Commit 04fc83d6 authored by Henning Janssen's avatar Henning Janssen

Changed calculation of offdiagonal green's function elements

Now the offdiagonal contributions are calculated on each atom
and only after the Kramers-Kronig transformation they are rotatet with the d_wgn matrices
parent f5f52a63
......@@ -67,7 +67,7 @@ MODULE m_gfcalc
SUBROUTINE postProcessGF(greensf,greensfCoeffs,atoms,input,sym,noco,vTot,hub1,results)
SUBROUTINE postProcessGF(greensf,greensfCoeffs,atoms,input,sym,noco,vTot,hub1,results,angle)
!contains all the modules for calculating properties from the greens function
USE m_onsite
......@@ -81,6 +81,7 @@ MODULE m_gfcalc
TYPE(t_hub1ham), INTENT(INOUT) :: hub1
TYPE(t_results), INTENT(IN) :: results
TYPE(t_potden), INTENT(IN) :: vTot
REAL, INTENT(IN) :: angle(sym%nop)
INTEGER i_gf,l,nType
COMPLEX mmpmat(-lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const,atoms%n_gf,3)
......@@ -88,7 +89,7 @@ MODULE m_gfcalc
CALL timestart("Green's Function: Postprocess")
!Perform the Kramer-Kronigs-Integration if we only have he imaginary part at this point
CALL calc_onsite(atoms,input,sym,noco,greensfCoeffs,greensf)
CALL calc_onsite(atoms,input,sym,noco,angle,greensfCoeffs,greensf)
!-------------------------------------------------------------
! Calculate various properties from the greens function
!-------------------------------------------------------------
......
......@@ -42,7 +42,7 @@ MODULE m_greensfImag21
INTEGER, INTENT(IN) :: ind(nbands,2)
REAL, INTENT(IN) :: eig(nbands)
INTEGER i_gf,nType,l,natom,ib,j,ie,m,lm,mp,lmp,imat,it,is,isi,ilo,ilop
INTEGER i_gf,nType,l,natom,ib,j,ie,m,lm,mp,lmp,imat,it,is,isi,ilo,ilop,nn
REAL fac
COMPLEX phase,weight
LOGICAL l_zero,l_tria
......@@ -58,12 +58,13 @@ MODULE m_greensfImag21
l = atoms%gfelem(i_gf)%l
DO natom = SUM(atoms%neq(:nType-1)) + 1, SUM(atoms%neq(:nType))
DO nn = 1, atoms%neq(nType)
natom = SUM(atoms%neq(:nType-1)) + nn
im = 0.0
!Loop through bands
!$OMP PARALLEL DEFAULT(none) &
!$OMP SHARED(natom,l,nType,wtkpt,i_gf,nbands,l_tria) &
!$OMP SHARED(natom,nn,l,nType,wtkpt,i_gf,nbands,l_tria) &
!$OMP SHARED(atoms,im,input,eigVecCoeffs,greensfCoeffs,denCoeffsOffDiag,eig,dosWeights,resWeights,ind) &
!$OMP PRIVATE(ie,m,mp,lm,lmp,weight,ib,j,l_zero,ilo,ilop)
!$OMP DO
......@@ -138,41 +139,47 @@ MODULE m_greensfImag21
!$OMP END PARALLEL
!Rotate the eqivalent atom into the irreducible brillouin zone
fac = 1.0/(sym%invarind(natom)*atoms%neq(nType))
IF(sym%invarind(natom).EQ.0) CALL juDFT_error("No symmetry operations",calledby="greensfImag")
DO imat = 1, MERGE(1,5,input%l_gfsphavg)
DO ie = 1, greensfCoeffs%ne
DO it = 1, sym%invarind(natom)
is = sym%invarop(natom,it)
isi = sym%invtab(is)
d_mat(:,:) = cmplx(0.0,0.0)
DO m = -l,l
DO mp = -l,l
d_mat(m,mp) = sym%d_wgn(m,mp,l,isi)
ENDDO
ENDDO
calc_mat = matmul( transpose( conjg(d_mat) ) , &
im(ie,-lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const,imat))
calc_mat = matmul( calc_mat, d_mat )
phase = exp(ImagUnit * angle(isi))
DO m = -l,l
DO mp = -l,l
IF(imat.EQ.1) THEN
greensfCoeffs%projdos(ie,i_gf,m,mp,3) = greensfCoeffs%projdos(ie,i_gf,m,mp,3) - AIMAG(fac * phase * conjg(calc_mat(m,mp)))
ELSE IF(imat.EQ.2) THEN
greensfCoeffs%uu(ie,i_gf,m,mp,3) = greensfCoeffs%uu(ie,i_gf,m,mp,3) - AIMAG(fac * phase * conjg(calc_mat(m,mp)))
ELSE IF(imat.EQ.3) THEN
greensfCoeffs%dd(ie,i_gf,m,mp,3) = greensfCoeffs%dd(ie,i_gf,m,mp,3) - AIMAG(fac * phase * conjg(calc_mat(m,mp)))
ELSE IF(imat.EQ.4) THEN
greensfCoeffs%ud(ie,i_gf,m,mp,3) = greensfCoeffs%ud(ie,i_gf,m,mp,3) - AIMAG(fac * phase * conjg(calc_mat(m,mp)))
ELSE IF(imat.EQ.5) THEN
greensfCoeffs%du(ie,i_gf,m,mp,3) = greensfCoeffs%du(ie,i_gf,m,mp,3) - AIMAG(fac * phase * conjg(calc_mat(m,mp)))
ENDIF
ENDDO
ENDDO
ENDDO!it
ENDDO!ie
ENDDO!imat
!fac = 1.0/(sym%invarind(natom)*atoms%neq(nType))
!IF(sym%invarind(natom).EQ.0) CALL juDFT_error("No symmetry operations",calledby="greensfImag")
!DO imat = 1, MERGE(1,5,input%l_gfsphavg)
! DO ie = 1, greensfCoeffs%ne
! DO it = 1, sym%invarind(natom)
! is = sym%invarop(natom,it)
! isi = sym%invtab(is)
! d_mat(:,:) = cmplx(0.0,0.0)
! DO m = -l,l
! DO mp = -l,l
! d_mat(m,mp) = sym%d_wgn(m,mp,l,isi)
! ENDDO
! ENDDO
! calc_mat = matmul( transpose( conjg(d_mat) ) , &
! im(ie,-lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const,imat))
! calc_mat = matmul( calc_mat, d_mat )
! !phase = exp(ImagUnit * angle(isi))
! DO m = -l,l
! DO mp = -l,l
! IF(imat.EQ.1) THEN
DO ie = 1, greensfCoeffs%ne
DO m = -l,l
DO mp = -l,l
greensfCoeffs%projdos21(ie,i_gf,nn,m,mp) = greensfCoeffs%projdos21(ie,i_gf,nn,m,mp) - AIMAG(im(ie,m,mp,1))
ENDDO
ENDDO
ENDDO
!ELSE IF(imat.EQ.2) THEN
! greensfCoeffs%uu(ie,i_gf,m,mp,3) = greensfCoeffs%uu(ie,i_gf,m,mp,3) - AIMAG(fac * conjg(calc_mat(m,mp)))
!ELSE IF(imat.EQ.3) THEN
! greensfCoeffs%dd(ie,i_gf,m,mp,3) = greensfCoeffs%dd(ie,i_gf,m,mp,3) - AIMAG(fac * conjg(calc_mat(m,mp)))
!ELSE IF(imat.EQ.4) THEN
! greensfCoeffs%ud(ie,i_gf,m,mp,3) = greensfCoeffs%ud(ie,i_gf,m,mp,3) - AIMAG(fac * conjg(calc_mat(m,mp)))
!ELSE IF(imat.EQ.5) THEN
! greensfCoeffs%du(ie,i_gf,m,mp,3) = greensfCoeffs%du(ie,i_gf,m,mp,3) - AIMAG(fac * conjg(calc_mat(m,mp)))
!ENDIF
!ENDDO
!ENDDO
!ENDDO!it
!ENDDO!ie
!ENDDO!imat
ENDDO !natom
ENDDO !i_gf
......
......@@ -24,7 +24,7 @@ INTEGER, PARAMETER :: int_method(3) = (/3,3,1/)
CONTAINS
SUBROUTINE calc_onsite(atoms,input,sym,noco,greensfCoeffs,g)
SUBROUTINE calc_onsite(atoms,input,sym,noco,angle,greensfCoeffs,g)
USE m_kkintgr
USE m_kk_cutoff
......@@ -38,10 +38,15 @@ SUBROUTINE calc_onsite(atoms,input,sym,noco,greensfCoeffs,g)
TYPE(t_sym), INTENT(IN) :: sym
TYPE(t_noco), INTENT(IN) :: noco
TYPE(t_input), INTENT(IN) :: input
REAL, INTENT(IN) :: angle(sym%nop)
!-Local Scalars
INTEGER i_gf,ie,l,m,mp,nType,jspin,ipm,kkcut,lp,nTypep,spin_cut
INTEGER i_gf,ie,l,m,mp,nType,jspin,ipm,kkcut,lp,nTypep,spin_cut,nn,natom
REAL fac
INTEGER it,is,isi
COMPLEX phase
COMPLEX d_mat(-lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const),calc_mat(-lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const)
COMPLEX g21(g%nz,-lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const)
DO i_gf = 1, atoms%n_gf
l = atoms%gfelem(i_gf)%l
......@@ -67,7 +72,7 @@ SUBROUTINE calc_onsite(atoms,input,sym,noco,greensfCoeffs,g)
!Perform the Kramers-Kronig-Integration if not already calculated
!
CALL timestart("On-Site: Kramer-Kronigs-Integration")
DO jspin = 1, MERGE(3,input%jspins,input%l_gfmperp)
DO jspin = 1, input%jspins
spin_cut = MERGE(1,jspin,jspin.GT.2)
kkcut = greensfCoeffs%kkintgr_cutoff(i_gf,spin_cut,2)
DO m= -l,l
......@@ -91,6 +96,45 @@ SUBROUTINE calc_onsite(atoms,input,sym,noco,greensfCoeffs,g)
ENDDO
ENDDO
ENDDO
IF(input%l_gfmperp) THEN
kkcut = greensfCoeffs%kkintgr_cutoff(i_gf,1,2)
DO nn = 1, atoms%neq(nType)
natom = SUM(atoms%neq(:nType-1)) + nn
DO ipm = 1, 2 !upper or lower half of the complex plane (G(E \pm i delta))
DO m= -l,l
DO mp= -lp,lp
CALL kkintgr(greensfCoeffs%projdos21(1:kkcut,i_gf,nn,m,mp),greensfCoeffs%e_bot,greensfCoeffs%del,kkcut,&
g21(:,m,mp),g%e,(ipm.EQ.2),g%mode,g%nz,int_method(g%mode))
ENDDO
ENDDO
fac = 1.0/(sym%invarind(natom)*atoms%neq(nType))
IF(sym%invarind(natom).EQ.0) CALL juDFT_error("No symmetry operations available",calledby="greensfImag")
DO ie = 1, g%nz
DO it = 1, sym%invarind(natom)
is = sym%invarop(natom,it)
isi = sym%invtab(is)
d_mat(:,:) = cmplx(0.0,0.0)
DO m = -l,l
DO mp = -l,l
d_mat(m,mp) = sym%d_wgn(m,mp,l,isi)
ENDDO
ENDDO
calc_mat = matmul( transpose( conjg(d_mat) ) , &
g21(ie,-lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const))
calc_mat = matmul( calc_mat, d_mat )
phase = exp(ImagUnit*angle(isi))
DO m = -l,l
DO mp = -l,l
g%gmmpMat(ie,i_gf,m,mp,3,ipm) =&
g%gmmpMat(ie,i_gf,m,mp,3,ipm) + phase * fac * conjg(calc_mat(m,mp))
ENDDO
ENDDO
ENDDO!it
ENDDO!ie
ENDDO
ENDDO!natom
ENDIF
CALL timestop("On-Site: Kramer-Kronigs-Integration")
!IF(input%l_gfmperp) THEN
! CALL rot_gf_mat(g,noco)
......
......@@ -351,7 +351,7 @@ MODULE m_cdn_io
IF (iofl < 0) EXIT
numLines = numLines + 1
END DO
IF (MOD(numLines,14*input%jspins).NE.0) THEN
IF (MOD(numLines,14*SIZE(den%mmpMat,4)).NE.0) THEN
WRITE(*,*) 'The n_mmp_mat file could not be read.'
WRITE(*,*) 'Was it an old style file with linear mixing parameter specification'
WRITE(*,*) 'in the last line? Linear mixing for the density matrix can now be'
......
......@@ -146,7 +146,7 @@ SUBROUTINE cdngen(eig_id,mpi,input,banddos,sliceplot,vacuum,&
IF(PRESENT(gOnsite).AND.mpi%irank.EQ.0) THEN
IF(atoms%n_gf.GT.0) THEN
CALL postProcessGF(gOnsite,greensfCoeffs,atoms,input,sym,noco,vTot,hub1,results)
CALL postProcessGF(gOnsite,greensfCoeffs,atoms,input,sym,noco,vTot,hub1,results,angle)
ENDIF
ENDIF
......
......@@ -451,10 +451,10 @@ CONTAINS
ENDIF
DEALLOCATE(r_b)
IF(input%l_gfmperp.AND.jspin.EQ.1) THEN
n = greensfCoeffs%ne*atoms%n_gf*(2*lmaxU_const+1)**2
n = greensfCoeffs%ne*atoms%n_gf*(2*lmaxU_const+1)**2*MAXVAL(atoms%neq)
ALLOCATE(r_b(n))
CALL MPI_REDUCE(greensfCoeffs%projdos(:,:,-lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const,3),r_b,n,CPP_MPI_REAL,MPI_SUM,0,MPI_COMM_WORLD,ierr)
IF(mpi%irank.EQ.0) CALL CPP_BLAS_scopy(n,r_b,1,greensfCoeffs%projdos(:,:,-lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const,3),1)
CALL MPI_REDUCE(greensfCoeffs%projdos21(:,:,:,-lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const),r_b,n,CPP_MPI_REAL,MPI_SUM,0,MPI_COMM_WORLD,ierr)
IF(mpi%irank.EQ.0) CALL CPP_BLAS_scopy(n,r_b,1,greensfCoeffs%projdos21(:,:,:,-lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const),1)
IF(.NOT.input%l_gfsphavg) THEN
CALL MPI_REDUCE(greensfCoeffs%uu(:,:,-lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const,3),r_b,n,CPP_MPI_REAL,MPI_SUM,0,MPI_COMM_WORLD,ierr)
IF(mpi%irank.EQ.0) CALL CPP_BLAS_scopy(n,r_b,1,greensfCoeffs%uu(:,:,-lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const,3),1)
......
......@@ -41,6 +41,7 @@ MODULE m_types_greensfCoeffs
!If we look at the Green's function that only depends on Energy and not on spatial arguments
!the imaginary part is equal to the proected density of states
REAL, ALLOCATABLE :: projdos(:,:,:,:,:)
REAL, ALLOCATABLE :: projdos21(:,:,:,:,:)
! These arrays are only used in the case we want the green's function with radial dependence
REAL, ALLOCATABLE :: uu(:,:,:,:,:)
......@@ -92,8 +93,10 @@ MODULE m_types_greensfCoeffs
IF(atoms%n_gf.GT.0) THEN
ALLOCATE(thisGREENSFCOEFFS%kkintgr_cutoff(atoms%n_gf,input%jspins,2))
ALLOCATE (thisGREENSFCOEFFS%projdos(thisGREENSFCOEFFS%ne,MAX(1,atoms%n_gf),-lmax:lmax,-lmax:lmax,spin_dim))
ALLOCATE (thisGREENSFCOEFFS%projdos(thisGREENSFCOEFFS%ne,MAX(1,atoms%n_gf),-lmax:lmax,-lmax:lmax,input%jspins))
ALLOCATE (thisGREENSFCOEFFS%projdos21(thisGREENSFCOEFFS%ne,MAX(1,atoms%n_gf),MAXVAL(atoms%neq),-lmax:lmax,-lmax:lmax))
thisGREENSFCOEFFS%projdos = 0.0
thisGREENSFCOEFFS%projdos21 = 0.0
IF(.NOT.input%l_gfsphavg) THEN
ALLOCATE (thisGREENSFCOEFFS%uu(thisGREENSFCOEFFS%ne,MAX(1,atoms%n_gf),-lmax:lmax,-lmax:lmax,spin_dim))
ALLOCATE (thisGREENSFCOEFFS%dd(thisGREENSFCOEFFS%ne,MAX(1,atoms%n_gf),-lmax:lmax,-lmax:lmax,spin_dim))
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment