Commit 7abfde37 authored by Henning Janssen's avatar Henning Janssen

Reset multiple commits

parent 9c7af07f
......@@ -12,7 +12,7 @@ CONTAINS
SUBROUTINE cdnval(eig_id, mpi,kpts,jspin,noco,input,banddos,cell,atoms,enpara,stars,&
vacuum,dimension,sphhar,sym,vTot,oneD,cdnvalJob,den,regCharges,dos,results,&
moments,hub1,coreSpecInput,mcd,slab,orbcomp,greensfCoeffs,angle,ntria,as,itria,atr)
moments,hub1,coreSpecInput,mcd,slab,orbcomp,greensfCoeffs,greensf,angle,ntria,as,itria,atr)
!************************************************************************************
! This is the FLEUR valence density generator
......@@ -38,10 +38,8 @@ SUBROUTINE cdnval(eig_id, mpi,kpts,jspin,noco,input,banddos,cell,atoms,enpara,st
USE m_pwden
USE m_forcea8
USE m_checkdopall
USE m_onsite ! calculates the on-site green's function
USE m_greensfImag
USE m_greensfImag21
USE m_tetra_weights
USE m_tetrahedronInit
USE m_gfcalc
USE m_cdnmt ! calculate the density and orbital moments etc.
USE m_orbmom ! coeffd for orbital moments
USE m_qmtsl ! These subroutines divide the input%film into vacuum%layers
......@@ -85,6 +83,7 @@ SUBROUTINE cdnval(eig_id, mpi,kpts,jspin,noco,input,banddos,cell,atoms,enpara,st
TYPE(t_slab), OPTIONAL, INTENT(INOUT) :: slab
TYPE(t_orbcomp), OPTIONAL, INTENT(INOUT) :: orbcomp
TYPE(t_greensfCoeffs), OPTIONAL, INTENT(INOUT) :: greensfCoeffs
TYPE(t_greensf), OPTIONAL, INTENT(INOUT) :: greensf
INTEGER, OPTIONAL, INTENT(IN) :: ntria
REAL, OPTIONAL, INTENT(IN) :: as
INTEGER, OPTIONAL, INTENT(IN) :: itria(3,2*kpts%nkpt)
......@@ -108,7 +107,7 @@ SUBROUTINE cdnval(eig_id, mpi,kpts,jspin,noco,input,banddos,cell,atoms,enpara,st
REAL,ALLOCATABLE :: we(:),eig(:)
INTEGER,ALLOCATABLE :: ev_list(:)
REAL, ALLOCATABLE :: f(:,:,:,:),g(:,:,:,:),flo(:,:,:,:) ! radial functions
REAL, ALLOCATABLE :: tetweights(:,:)
COMPLEX, ALLOCATABLE :: tetweights(:,:)
INTEGER, ALLOCATABLE :: tet_ind(:,:)
TYPE (t_lapw) :: lapw
......@@ -235,17 +234,10 @@ SUBROUTINE cdnval(eig_id, mpi,kpts,jspin,noco,input,banddos,cell,atoms,enpara,st
! valence density in the atomic spheres
CALL eigVecCoeffs%init(input,DIMENSION,atoms,noco,jspin,noccbd)
IF (atoms%n_gf.GT.0.AND.(input%tria.OR.input%gfTet)) THEN
IF(input%film) THEN
CALL timestart("TriangularWeights")
tetweights = 0.0
CALL tria_weights(ikpt,kpts,ntria,as,itria,atr,SIZE(ev_list),results%eig(ev_list,:,jsp),greensfCoeffs,tetweights,tet_ind,results%ef)
CALL timestop("TriangularWeights")
ELSE
CALL timestart("TetrahedronWeights")
tetweights = 0.0
CALL tetra_weights(ikpt,kpts,SIZE(ev_list),results%eig(ev_list,:,jsp),greensfCoeffs,tetweights,tet_ind,results%ef)
CALL timestop("TetrahedronWeights")
ENDIF
CALL timestart("TetrahedronWeights")
CALL tetrahedronInit(ikpt,kpts,input,SIZE(ev_list),results%eig(ev_list,:,jsp),&
greensfCoeffs,results%ef,greensf%e,greensf%nz,tetweights,tet_ind)
CALL timestop("TetrahedronWeights")
ENDIF
DO ispin = jsp_start, jsp_end
IF (input%l_f) CALL force%init2(noccbd,input,atoms)
......@@ -254,16 +246,8 @@ SUBROUTINE cdnval(eig_id, mpi,kpts,jspin,noco,input,banddos,cell,atoms,enpara,st
eigVecCoeffs%ccof(-atoms%llod:,:,:,:,ispin),zMat,eig,force)
IF (atoms%n_u.GT.0) CALL n_mat(atoms,sym,noccbd,usdus,ispin,we,eigVecCoeffs,den%mmpMat(:,:,:,jspin))
IF (atoms%n_gf.GT.0) THEN
CALL timestart("Greens Function: Imaginary Part")
CALL greensfImag(atoms,sym,input,ispin,noccbd,tetweights,tet_ind,kpts%wtkpt(ikpt),eig,usdus,eigVecCoeffs,greensfCoeffs)
IF(input%l_gfmperp.AND.ispin==jsp_end) THEN
CALL greensfImag21(atoms,sym,angle,input,noccbd,tetweights,tet_ind,kpts%wtkpt(ikpt),eig,denCoeffsOffdiag,eigVecCoeffs,greensfCoeffs)
ENDIF
CALL timestop("Greens Function: Imaginary Part")
ENDIF
IF (atoms%n_gf.GT.0) CALL bzIntegrationGF(atoms,sym,input,angle,ispin,noccbd,tetweights,tet_ind,kpts%wtkpt(ikpt),&
eig,denCoeffsOffdiag,usdus,eigVecCoeffs,greensf,greensfCoeffs,ispin==jsp_end)
! perform Brillouin zone integration and summation over the
! bands in order to determine the energy parameters for each atom and angular momentum
......@@ -301,7 +285,8 @@ SUBROUTINE cdnval(eig_id, mpi,kpts,jspin,noco,input,banddos,cell,atoms,enpara,st
#ifdef CPP_MPI
DO ispin = jsp_start,jsp_end
CALL mpi_col_den(mpi,sphhar,atoms,oneD,stars,vacuum,input,noco,ispin,regCharges,dos,&
results,denCoeffs,orb,denCoeffsOffdiag,den,den%mmpMat(:,:,:,jspin),mcd,slab,orbcomp,greensfCoeffs)
results,denCoeffs,orb,denCoeffsOffdiag,den,den%mmpMat(:,:,:,jspin),&
mcd,slab,orbcomp,greensfCoeffs,greensf)
END DO
#endif
......
......@@ -46,7 +46,11 @@ CONTAINS
IF (atoms%jri(ntyp)>atoms%jmtd) CALL juDFT_error("atoms%jri(ntyp).GT.atoms%jmtd",calledby ="sorad")
ALLOCATE ( p(atoms%jmtd,2),pd(atoms%jmtd,2),q(atoms%jmtd,2),plo(atoms%jmtd,2),fint(atoms%jmtd),&
& qlo(atoms%jmtd,2),plop(atoms%jmtd,2),qd(atoms%jmtd,2),v0(atoms%jmtd),vso(atoms%jmtd,2) )
!
p = 0.0 ; pd = 0.0 ; q = 0.0 ; plo = 0.0 ; fint = 0.0
qlo = 0.0 ; plop = 0.0 ; qd = 0.0 ; v0 = 0.0 ; vso = 0.0
DO l = 0,atoms%lmax(ntyp)
DO jspin = 1,input%jspins
......@@ -130,6 +134,7 @@ CONTAINS
!+apw+lo
IF (atoms%l_dulo(ilo,ntyp).OR.atoms%ulo_der(ilo,ntyp).GE.1) THEN ! calculate energy derivative (of order atoms%ulo_der) at e
ALLOCATE (glo(atoms%jmtd,2),pqlo(atoms%jmtd,2),filo(atoms%jmtd,2))
glo = 0.0 ; pqlo = 0.0 ; filo = 0.0
pqlo(1:atoms%jri(ntyp),1)=plo(1:atoms%jri(ntyp),jspin)
pqlo(1:atoms%jri(ntyp),2)=qlo(1:atoms%jri(ntyp),jspin)
i = atoms%ulo_der(ilo,ntyp)
......@@ -190,6 +195,7 @@ CONTAINS
!+apw+lo
IF (atoms%l_dulo(ilo,ntyp).OR.atoms%ulo_der(ilo,ntyp).GE.1) THEN ! calculate orthogonal energy derivative at e
ALLOCATE (glo(atoms%jmtd,2),pqlo(atoms%jmtd,2),filo(atoms%jmtd,2))
glo = 0.0 ; pqlo = 0.0 ; filo = 0.0
pqlo(1:atoms%jri(ntyp),1)=plop(1:atoms%jri(ntyp),jspin)
pqlo(1:atoms%jri(ntyp),2)=q(1:atoms%jri(ntyp),1)
i = atoms%ulo_der(ilo,ntyp)
......
......@@ -7,8 +7,7 @@ MODULE m_fermie
! or fermi-function p.kurz
!----------------------------------------------------------------------
CONTAINS
SUBROUTINE fermie(eig_id, mpi,kpts,obsolete,&
input, noco,e_min,cell,results)
SUBROUTINE fermie(eig_id,mpi,kpts,input,noco,e_min,cell,results)
!---------------------------------------------------f--------------------
!
......@@ -39,13 +38,12 @@ CONTAINS
USE m_types
USE m_xmlOutput
IMPLICIT NONE
TYPE(t_results),INTENT(INOUT) :: results
TYPE(t_mpi),INTENT(IN) :: mpi
TYPE(t_obsolete),INTENT(IN) :: obsolete
TYPE(t_input),INTENT(IN) :: input
TYPE(t_noco),INTENT(IN) :: noco
TYPE(t_cell),INTENT(IN) :: cell
TYPE(t_kpts),INTENT(IN) :: kpts
TYPE(t_results), INTENT(INOUT) :: results
TYPE(t_mpi), INTENT(IN) :: mpi
TYPE(t_input), INTENT(IN) :: input
TYPE(t_noco), INTENT(IN) :: noco
TYPE(t_cell), INTENT(IN) :: cell
TYPE(t_kpts), INTENT(IN) :: kpts
! ..
! .. Scalar Arguments ..
INTEGER, INTENT (IN) :: eig_id
......@@ -61,8 +59,9 @@ CONTAINS
! ..
! .. Local Arrays ..
!
INTEGER, ALLOCATABLE :: idxeig(:),idxjsp(:),idxkpt(:),INDEX(:)
REAL, ALLOCATABLE :: e(:),eig(:,:,:),we(:)
INTEGER :: idxeig(SIZE(results%w_iks)),idxjsp(SIZE(results%w_iks)),idxkpt(SIZE(results%w_iks)),INDEX(SIZE(results%w_iks))
REAL :: e(SIZE(results%w_iks)),we(SIZE(results%w_iks))
REAL, ALLOCATABLE :: eig(:,:,:)
INTEGER ne(kpts%nkpt,SIZE(results%w_iks,3))
CHARACTER(LEN=20) :: attributes(5)
......@@ -94,15 +93,12 @@ CONTAINS
!***********************************************************************
! .. Data statements ..
DATA del/1.0e-6/
! ..
n=SIZE(results%w_iks) !size of list of all eigenvalues
ALLOCATE (idxeig(n),idxjsp(n),idxkpt(n),INDEX(n),e(n),we(n) )
ALLOCATE (eig(SIZE(results%w_iks,1),SIZE(results%w_iks,2),SIZE(results%w_iks,3)))
! initiliaze e
e = 0
IF ( mpi%irank == 0 ) WRITE (6,FMT=8000)
8000 FORMAT (/,/,1x,'fermi energy and band-weighting factors:')
!
......@@ -126,17 +122,19 @@ CONTAINS
DO jsp = 1,nspins
DO k = 1,kpts%nkpt
CALL read_eig(eig_id,k,jsp,neig=ne(k,jsp),eig=eig(:,k,jsp))
IF ( mpi%irank == 0 ) THEN
IF (mpi%irank == 0) THEN
WRITE (6,'(a2,3f10.5,f12.6)') 'at',kpts%bk(:,k),kpts%wtkpt(k)
WRITE (6,'(i5,a14)') ne(k,jsp),' eigenvalues :'
WRITE (6,'(8f12.6)') (eig(i,k,jsp),i=1,ne(k,jsp))
attributes = ''
WRITE(attributes(1),'(i0)') jsp
WRITE(attributes(2),'(i0)') k
WRITE(attributes(3),'(f15.8)') kpts%bk(1,k)
WRITE(attributes(4),'(f15.8)') kpts%bk(2,k)
WRITE(attributes(5),'(f15.8)') kpts%bk(3,k)
CALL writeXMLElementPoly('eigenvaluesAt',(/'spin','ikpt','k_x ','k_y ','k_z '/),attributes,eig(1:ne(k,jsp),k,jsp))
IF(.NOT.judft_was_argument("-minimalOutput")) THEN
attributes = ''
WRITE(attributes(1),'(i0)') jsp
WRITE(attributes(2),'(i0)') k
WRITE(attributes(3),'(f15.8)') kpts%bk(1,k)
WRITE(attributes(4),'(f15.8)') kpts%bk(2,k)
WRITE(attributes(5),'(f15.8)') kpts%bk(3,k)
CALL writeXMLElementPoly('eigenvaluesAt',(/'spin','ikpt','k_x ','k_y ','k_z '/),attributes,eig(1:ne(k,jsp),k,jsp))
END IF
END IF
END DO
ENDDO
......@@ -253,7 +251,7 @@ CONTAINS
ENDIF
efermi = results%ef
enddo
DEALLOCATE ( idxeig,idxjsp,idxkpt,index,e,eig,we )
DEALLOCATE (eig)
attributes = ''
WRITE(attributes(1),'(f20.10)') results%ef
......
......@@ -4,10 +4,14 @@ set(fleur_F77 ${fleur_F77}
set(fleur_F90 ${fleur_F90}
greensf/greensfImag.f90
greensf/greensfImag21.f90
greensf/greensfRes.f90
greensf/greensfRes21.f90
greensf/angles.f90
greensf/onsite.f90
greensf/kkintgr.f90
greensf/tetra_weights.f90
greensf/dosWeights.f90
greensf/resWeights.f90
greensf/tetrahedronInit.f90
greensf/kk_cutoff.f90
greensf/gfcalc.f90
greensf/crystalfield.f90
......
......@@ -10,7 +10,6 @@ MODULE m_calc_tetra
USE m_types
USE m_juDFT
USE m_constants
USE m_gen_bz
IMPLICIT NONE
......@@ -25,32 +24,6 @@ MODULE m_calc_tetra
REAL sumvol,volbz
CALL timestart("Calculation of Tetrahedra")
!We calculate the kpts in the whole bz
IF(kpts%nkptf.EQ.0) THEN
WRITE(*,*) "Recalculating the kpoints over the full brillouin zone"
CALL gen_bz(kpts,sym)
ENDIF
kpts%nkpt = kpts%nkptf
DEALLOCATE(kpts%bk,kpts%wtkpt)
ALLOCATE(kpts%bk(3,kpts%nkptf),kpts%wtkpt(kpts%nkptf))
kpts%bk(:,:) = kpts%bkf(:,:)
IF (kpts%nkpt3(1)*kpts%nkpt3(2)*kpts%nkpt3(3).NE.kpts%nkptf) THEN
IF(kpts%l_gamma) THEN
kpts%wtkpt = 1.0 / (kpts%nkptf-1)
DO ikpt = 1, kpts%nkptf
IF(ALL(kpts%bk(:,ikpt).EQ.0.0)) THEN
kpts%wtkpt(ikpt) = 0.0
END IF
END DO
ELSE
CALL juDFT_error("nkptf does not match product of nkpt3(i).",calledby="calc_tetra")
END IF
ELSE
kpts%wtkpt = 1.0 / kpts%nkptf
END IF
volbz = sqrt(sum(cell%bmat(:,1)*cell%bmat(:,1))*sum(cell%bmat(:,2)*cell%bmat(:,2))*sum(cell%bmat(:,3)*cell%bmat(:,3)))
!Choose the tetrahedra decomposition along the shortest diagonal
......@@ -58,8 +31,8 @@ MODULE m_calc_tetra
!MISSING: Generate all kpts
!Set up pointer array for the kpts
p = 0
DO ikpt = 1, kpts%nkpt
iarr = nint(kpts%bk(:,ikpt)*kpts%nkpt3)
DO ikpt = 1, kpts%nkptf
iarr = nint(kpts%bkf(:,ikpt)*kpts%nkpt3)
p(iarr(1),iarr(2),iarr(3)) = ikpt
ENDDO
!Wrap around at the end
......@@ -73,7 +46,7 @@ MODULE m_calc_tetra
p(:,:,0) = p(:,:,kpts%nkpt3(3))
ENDIF
!Check for invalid indices
IF(ANY(p<=0).OR.ANY(p>kpts%nkpt)) CALL juDFT_error("Invalid kpoint index in pointer array",calledby="calc_tetra")
IF(ANY(p<=0).OR.ANY(p>kpts%nkptf)) CALL juDFT_error("Invalid kpoint index in pointer array",calledby="calc_tetra")
IF (ALLOCATED(kpts%ntetra)) THEN
DEALLOCATE(kpts%ntetra)
......@@ -81,8 +54,8 @@ MODULE m_calc_tetra
IF (ALLOCATED(kpts%voltet)) THEN
DEALLOCATE(kpts%voltet)
END IF
ALLOCATE(kpts%ntetra(4,kpts%nkpt*6))
ALLOCATE(kpts%voltet(kpts%nkpt*6))
ALLOCATE(kpts%ntetra(4,kpts%nkptf*6))
ALLOCATE(kpts%voltet(kpts%nkptf*6))
kpts%ntet = 0
sumvol = 0.0
......@@ -110,7 +83,7 @@ MODULE m_calc_tetra
ENDDO
ENDDO
ENDDO
IF(ABS(sumvol-volbz).GT.0.001) CALL juDFT_error("calc_tetra failed", calledby="calc_tetra")
IF(ABS(sumvol-volbz).GT.1E-12) CALL juDFT_error("calc_tetra failed", calledby="calc_tetra")
input%gfTet = .TRUE.
kpts%voltet = kpts%voltet*kpts%ntet/volbz
CALL timestop("Calculation of Tetrahedra")
......
MODULE m_gfcalc
!contains all the modules for calculating properties from the greens function
!Contains the main routines called from cdnval and cdngen
USE m_juDFT
USE m_constants
!These dependencies are here so that all these subroutines can be called by simply adding USE m_gfcalc tot hte file
USE m_j0
USE m_gfDOS
USE m_occmtx
USE m_hybridization
USE m_crystalfield
CONTAINS
SUBROUTINE bzIntegrationGF(atoms,sym,input,angle,ispin,nbands,tetWeights,tetInd,wtkpt,eig,denCoeffsOffdiag,&
usdus,eigVecCoeffs,greensf,greensfCoeffs,l21)
USE m_greensfImag
USE m_greensfImag21
USE m_greensfRes
USE m_greensfRes21
IMPLICIT NONE
!-Type Arguments
TYPE(t_atoms), INTENT(IN) :: atoms
TYPE(t_sym), INTENT(IN) :: sym
TYPE(t_input), INTENT(IN) :: input
TYPE(t_eigVecCoeffs), INTENT(IN) :: eigVecCoeffs
TYPE(t_usdus), INTENT(IN) :: usdus
TYPE(t_denCoeffsOffDiag), INTENT(IN) :: denCoeffsOffdiag
TYPE(t_greensfCoeffs), INTENT(INOUT) :: greensfCoeffs
TYPE(t_greensf), INTENT(INOUT) :: greensf
!-Scalar Arguments
INTEGER, INTENT(IN) :: ispin !Current spin index
INTEGER, INTENT(IN) :: nbands !Number of bands to be considered
REAL, INTENT(IN) :: wtkpt !Weight of the current k-point
LOGICAL, INTENT(IN) :: l21 !Calculate spin off-diagonal part ?
!-Array Arguments
COMPLEX, INTENT(IN) :: tetWeights(greensfCoeffs%ne,nbands) !Precalculated tetrahedron weights for the current k-point
INTEGER, INTENT(IN) :: tetInd(nbands,2) !Gives the range where the tetrahedron weights are non-zero
REAL, INTENT(IN) :: eig(nbands) !Eigenvalues for the current k-point
REAL, INTENT(IN) :: angle(sym%nop) !Phases for spin-offdiagonal part
IF(input%l_resolvent) THEN
!Calculate greens function directly (very time consuming)
CALL timestart("Greens Function: Resolvent")
CALL greensfRes(atoms,sym,input,ispin,nbands,tetWeights,tetInd,wtkpt,eig,usdus,eigVecCoeffs,greensf)
IF(input%l_gfmperp.AND.l21) THEN
CALL greensfRes21(atoms,sym,angle,input,nbands,tetWeights,tetInd,wtkpt,eig,denCoeffsOffdiag,eigVecCoeffs,greensf)
ENDIF
CALL timestop("Greens Function: Resolvent")
ELSE
CALL timestart("Greens Function: Imaginary Part")
CALL greensfImag(atoms,sym,input,ispin,nbands,tetWeights,tetInd,wtkpt,eig,usdus,eigVecCoeffs,greensfCoeffs)
IF(input%l_gfmperp.AND.l21) THEN
CALL greensfImag21(atoms,sym,angle,input,nbands,tetWeights,tetInd,wtkpt,eig,denCoeffsOffdiag,eigVecCoeffs,greensfCoeffs)
ENDIF
CALL timestop("Greens Function: Imaginary Part")
ENDIF
END SUBROUTINE bzIntegrationGF
SUBROUTINE postProcessGF(greensf,greensfCoeffs,atoms,input,sym,noco,vTot,hub1,results)
!contains all the modules for calculating properties from the greens function
USE m_onsite
TYPE(t_atoms), INTENT(IN) :: atoms
TYPE(t_input), INTENT(IN) :: input
TYPE(t_sym), INTENT(IN) :: sym
TYPE(t_noco), INTENT(IN) :: noco
TYPE(t_greensfCoeffs), INTENT(INOUT) :: greensfCoeffs
TYPE(t_greensf), INTENT(INOUT) :: greensf
TYPE(t_hub1ham), INTENT(INOUT) :: hub1
TYPE(t_results), INTENT(IN) :: results
TYPE(t_potden), INTENT(IN) :: vTot
INTEGER i_gf,l,nType
COMPLEX mmpmat(-lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const,atoms%n_gf,3)
CALL timestart("Green's Function: Postprocess")
!Perform the Kramer-Kronigs-Integration if we only have he imaginary part at this point
IF(.NOT.input%l_resolvent) THEN
CALL calc_onsite(atoms,input,sym,noco,greensfCoeffs,greensf)
ENDIF
!-------------------------------------------------------------
! Calculate various properties from the greens function
!-------------------------------------------------------------
!calculate the crystal field contribution to the local hamiltonian in LDA+Hubbard 1
IF(atoms%n_hia.GT.0.AND.ANY(hub1%ccf(:).NE.0.0).AND..NOT.input%l_resolvent) THEN
CALL crystal_field(atoms,input,greensfCoeffs,hub1,vTot)
ENDIF
IF(input%jspins.EQ.2.AND..NOT.input%l_resolvent) THEN
CALL eff_excinteraction(greensf,atoms,input,greensfCoeffs)
ENDIF
DO i_gf = 1, atoms%n_gf
l = atoms%gfelem(i_gf)%l
nType = atoms%gfelem(i_gf)%atomType
IF(l.NE.atoms%gfelem(i_gf)%lp) CYCLE
IF(nType.NE.atoms%gfelem(i_gf)%atomTypep) CYCLE
!Occupation matrix
CALL occmtx(greensf,l,nType,atoms,sym,input,mmpmat(:,:,i_gf,:),l_write=.TRUE.)
!Hybridization function
CALL hybridization(greensf,l,nType,atoms,input,results%ef)
!Density of states from Greens function
CALL gfDOS(greensf,l,nType,i_gf,atoms,input)
ENDDO
CALL timestop("Green's Function: Postprocess")
END SUBROUTINE postProcessGF
END MODULE m_gfcalc
\ No newline at end of file
......@@ -30,7 +30,7 @@ SUBROUTINE greensfImag(atoms,sym,input,ispin,nbands,tetweights,ind,wtkpt,eig,usd
REAL, INTENT(IN) :: wtkpt !Weight of the current k-point
!-Array Arguments
REAL, INTENT(IN) :: tetweights(greensfCoeffs%ne,nbands) !Precalculated tetrahedron weights for the current k-point
COMPLEX, INTENT(IN) :: tetweights(greensfCoeffs%ne,nbands) !Precalculated tetrahedron weights for the current k-point
INTEGER, INTENT(IN) :: ind(nbands,2) !Gives the range where the tetrahedron weights are non-zero
REAL, INTENT(IN) :: eig(nbands) !Eigenvalues for the current k-point
......@@ -63,7 +63,7 @@ SUBROUTINE greensfImag(atoms,sym,input,ispin,nbands,tetweights,ind,wtkpt,eig,usd
l_zero = .true.
IF(l_tria) THEN
!TETRAHEDRON METHOD: check if the weight for this eigenvalue is non zero
IF(ANY(tetweights(ind(ib,1):ind(ib,2),ib).NE.0.0)) l_zero = .false.
IF(ANY(ABS(tetweights(ind(ib,1):ind(ib,2),ib)).NE.0.0)) l_zero = .false.
ELSE
!HISTOGRAM METHOD: check if eigenvalue is inside the energy range
j = NINT((eig(ib)-greensfCoeffs%e_bot)/greensfCoeffs%del)+1
......@@ -79,7 +79,7 @@ SUBROUTINE greensfImag(atoms,sym,input,ispin,nbands,tetweights,ind,wtkpt,eig,usd
!Choose the relevant energy points depending on the bz-integration method
DO ie = MERGE(ind(ib,1),j,l_tria), MERGE(ind(ib,2),j,l_tria)
!weight for the bz-integration including spin-degeneracy
weight = 2.0/input%jspins * MERGE(tetweights(ie,ib),wtkpt/greensfCoeffs%del,l_tria)
weight = 2.0/input%jspins * MERGE(REAL(tetweights(ie,ib)),wtkpt/greensfCoeffs%del,l_tria)
!
!Contribution from states
!
......
......@@ -37,7 +37,7 @@ MODULE m_greensfImag21
INTEGER, INTENT(IN) :: nbands
REAL, INTENT(IN) :: wtkpt
REAL, INTENT(IN) :: tetweights(greensfCoeffs%ne,nbands)
COMPLEX, INTENT(IN) :: tetweights(greensfCoeffs%ne,nbands)
INTEGER, INTENT(IN) :: ind(nbands,2)
REAL, INTENT(IN) :: eig(nbands)
......@@ -61,7 +61,7 @@ MODULE m_greensfImag21
im = 0.0
!Loop through bands
!$OMP PARALLEL DEFAULT(none) &
!$OMP SHARED(natom,l,nType,wtkpt,i_gf,nbands) &
!$OMP SHARED(natom,l,nType,wtkpt,i_gf,nbands,l_tria) &
!$OMP SHARED(atoms,im,input,eigVecCoeffs,greensfCoeffs,denCoeffsOffDiag,eig,tetweights,ind) &
!$OMP PRIVATE(ie,m,mp,lm,lmp,weight,ib,j,l_zero,ilo,ilop)
!$OMP DO
......@@ -70,7 +70,7 @@ MODULE m_greensfImag21
l_zero = .true.
IF(l_tria) THEN
!TETRAHEDRON METHOD: check if the weight for this eigenvalue is non zero
IF(ANY(tetweights(:,ib).NE.0.0)) l_zero = .false.
IF(ANY(ABS(tetweights(:,ib)).NE.0.0)) l_zero = .false.
ELSE
!HISTOGRAM METHOD: check if eigenvalue is inside the energy range
j = NINT((eig(ib)-greensfCoeffs%e_bot)/greensfCoeffs%del)+1
......@@ -85,7 +85,7 @@ MODULE m_greensfImag21
lmp = l*(l+1) + mp
DO ie = MERGE(ind(ib,1),j,l_tria), MERGE(ind(ib,2),j,l_tria)
weight = MERGE(tetweights(ie,ib),wtkpt/greensfCoeffs%del,l_tria)
weight = MERGE(REAL(tetweights(ie,ib)),wtkpt/greensfCoeffs%del,l_tria)
!
!Contribution from states
!
......
MODULE m_greensfRes
USE m_juDFT
USE m_types
USE m_constants
CONTAINS
SUBROUTINE greensfRes(atoms,sym,input,ispin,nbands,tetweights,ind,wtkpt,eig,usdus,eigVecCoeffs,greensf)
!This Subroutine calculates the contribution to the imaginary part of the Matrix elements G^[n \sigma]_{Lm Lm'}(E+i*sigma)
!of the current k-Point (it is called in cdnval) inside the MT-sphere
!and sums over the Brillouin-Zone using the histogram method or linear tetrahedron method
!It is essentially the l-density of states in a (m,mp) matrix with an additional factor - pi
IMPLICIT NONE
!-Type Arguments
TYPE(t_atoms), INTENT(IN) :: atoms
TYPE(t_sym), INTENT(IN) :: sym
TYPE(t_eigVecCoeffs), INTENT(IN) :: eigVecCoeffs
TYPE(t_usdus), INTENT(IN) :: usdus
TYPE(t_greensf), INTENT(INOUT) :: greensf
TYPE(t_input), INTENT(IN) :: input
!-Scalar Arguments
INTEGER, INTENT(IN) :: ispin !Current spin index
INTEGER, INTENT(IN) :: nbands !Number of bands to be considered
REAL, INTENT(IN) :: wtkpt !Weight of the current k-point
!-Array Arguments
COMPLEX, INTENT(IN) :: tetweights(greensf%nz,nbands) !Precalculated tetrahedron weights for the current k-point
INTEGER, INTENT(IN) :: ind(nbands,2) !Gives the range where the tetrahedron weights are non-zero
REAL, INTENT(IN) :: eig(nbands) !Eigenvalues for the current k-point
!-Local Scalars
LOGICAL l_zero,l_tria
INTEGER i_gf,ib,iz,j,nType,natom,l,m,mp,lm,lmp,ilo,ilop,imat,it,is,isi
REAL fac
COMPLEX weight
COMPLEX g(greensf%nz,-lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const,MERGE(1,5,input%l_gfsphavg))
COMPLEX d_mat(-lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const),calc_mat(-lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const)
l_tria = input%tria.OR.input%gfTet
IF(.NOT.l_tria) CALL juDFT_error("Only available with tetra",calledby="greensfRes")
IF(.NOT.input%l_gfsphavg) CALL juDFT_error("No radial dependence",calledby="greensfRes")
!Loop through the gf elements to be calculated
DO i_gf = 1, atoms%n_gf
l = atoms%gfelem(i_gf)%l
nType = atoms%gfelem(i_gf)%atomType
!Loop through equivalent atoms
DO natom = SUM(atoms%neq(:nType-1)) + 1, SUM(atoms%neq(:nType))
g = 0.0
!Loop through bands
!$OMP PARALLEL DEFAULT(none) &
!$OMP SHARED(natom,l,nType,ispin,wtkpt,i_gf,nbands,l_tria) &
!$OMP SHARED(atoms,g,input,eigVecCoeffs,usdus,greensf,eig,tetweights,ind) &
!$OMP PRIVATE(iz,m,mp,lm,lmp,ilo,ilop,weight,ib,j,l_zero)
!$OMP DO
DO ib = 1, nbands
DO m = -l, l
lm = l*(l+1)+m
DO mp = -l,l
lmp = l*(l+1)+mp
!Choose the relevant energy points depending on the bz-integration method
DO iz = 1, greensf%nz
!weight for the bz-integration including spin-degeneracy
weight = 2.0/input%jspins * tetweights(iz,ib)
!
!Contribution from states
!
g(iz,m,mp,1) = g(iz,m,mp,1) + weight *&
(conjg(eigVecCoeffs%acof(ib,lmp,natom,ispin))*eigVecCoeffs%acof(ib,lm,natom,ispin) +&
conjg(eigVecCoeffs%bcof(ib,lmp,natom,ispin))*eigVecCoeffs%bcof(ib,lm,natom,ispin) *&
usdus%ddn(l,nType,ispin))
!IF(.NOT.input%l_gfsphavg) THEN
! im(ie,m,mp,2) = im(ie,m,mp,2) - pi_const * weight *&
! conjg(eigVecCoeffs%acof(ib,lm,natom,ispin))*eigVecCoeffs%acof(ib,lmp,natom,ispin)
! im(ie,m,mp,3) = im(ie,m,mp,3) - pi_const * weight *&
! conjg(eigVecCoeffs%bcof(ib,lm,natom,ispin))*eigVecCoeffs%bcof(ib,lmp,natom,ispin)
! im(ie,m,mp,4) = im(ie,m,mp,4) - pi_const * weight *&
! conjg(eigVecCoeffs%acof(ib,lm,natom,ispin))*eigVecCoeffs%bcof(ib,lmp,natom,ispin)
! im(ie,m,mp,5) = im(ie,m,mp,5) - pi_const * weight *&
! conjg(eigVecCoeffs%bcof(ib,lm,natom,ispin))*eigVecCoeffs%acof(ib,lmp,natom,ispin)
!END IF
!
! add local orbital contribution (not implemented for radial dependence yet and not tested for average)
!
DO ilo = 1, atoms%nlo(nType)
IF(atoms%llo(ilo,nType).EQ.l) THEN
g(iz,m,mp,1) = g(iz,m,mp,1) + weight * ( usdus%uulon(ilo,nType,ispin) * (&
conjg(eigVecCoeffs%acof(ib,lmp,natom,ispin))*eigVecCoeffs%ccof(m,ib,ilo,natom,ispin) +&
conjg(eigVecCoeffs%ccof(mp,ib,ilo,natom,ispin))*eigVecCoeffs%acof(ib,lm,natom,ispin) )&
+ usdus%dulon(ilo,nType,ispin) * (&
conjg(eigVecCoeffs%bcof(ib,lmp,natom,ispin))*eigVecCoeffs%ccof(m,ib,ilo,natom,ispin) +&
conjg(eigVecCoeffs%ccof(mp,ib,ilo,natom,ispin))*eigVecCoeffs%bcof(ib,lm,natom,ispin)))
DO ilop = 1, atoms%nlo(nType)
IF (atoms%llo(ilop,nType).EQ.l) THEN
g(iz,m,mp,1) = g(iz,m,mp,1) + weight * usdus%uloulopn(ilo,ilop,nType,ispin) *&
conjg(eigVecCoeffs%ccof(mp,ib,ilop,natom,ispin)) *eigVecCoeffs%ccof(m,ib,ilo,natom,ispin)
ENDIF
ENDDO
ENDIF
ENDDO
ENDDO! ie
ENDDO !mp
ENDDO !m
ENDDO !ib
!$OMP END DO
!$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 available",calledby="greensfImag")
DO imat = 1, MERGE(1,5,input%l_gfsphavg)
DO iz = 1, greensf%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) ) , &
g(iz,-lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const,imat))
calc_mat = matmul( calc_mat, d_mat )
DO m = -l,l
DO mp = -l,l
IF(imat.EQ.1) THEN
greensf%gmmpmat(iz,i_gf,m,mp,ispin,1) = greensf%gmmpmat(iz,i_gf,m,mp,ispin,1) + fac * calc_mat(m,mp)
greensf%gmmpmat(iz,i_gf,m,mp,ispin,2) = greensf%gmmpmat(iz,i_gf,m,mp,ispin,2) + fac * conjg(calc_mat(m,mp))
!ELSE IF(imat.EQ.2) THEN
! greensfCoeffs%uu(ie,i_gf,m,mp,ispin) = greensfCoeffs%uu(ie,i_gf,m,mp,ispin) + fac * conjg(calc_mat(m,mp))
!ELSE IF(imat.EQ.3) THEN
! greensfCoeffs%dd(ie,i_gf,m,mp,ispin) = greensfCoeffs%dd(ie,i_gf,m,mp,ispin) + fac * conjg(calc_mat(m,mp))
!ELSE IF(imat.EQ.4) THEN
! greensfCoeffs%ud(ie,i_gf,m,mp,ispin) = greensfCoeffs%ud(ie,i_gf,m,mp,ispin) + fac * conjg(calc_mat(m,mp))
!ELSE IF(imat.EQ.5) THEN
! greensfCoeffs%du(ie,i_gf,m,mp,ispin) = greensfCoeffs%du(ie,i_gf,m,mp,ispin) + fac * conjg(calc_mat(m,mp))
ENDIF
ENDDO
ENDDO
ENDDO!it
ENDDO!ie
ENDDO!imat
ENDDO!natom
ENDDO !i_gf
END SUBROUTINE greensfRes
END MODULE m_greensfRes
\ No newline at end of file
MODULE m_greensfRes21
USE m_juDFT
USE m_types
USE m_constants
CONTAINS
SUBROUTINE greensfRes21(atoms,sym,angle,input,nbands,tetweights,ind,wtkpt,eig,denCoeffsOffDiag,eigVecCoeffs,greensf)
!This Subroutine calculates the contribution to the imaginary part of the Matrix elements G^[n \sigma]_{Lm Lm'}(E+i*sigma)