Commit 2533034f authored by Henning Janssen's avatar Henning Janssen

radial dependence calculation of Green's Function can now be switched between individual elements

parent f53fc1b0
......@@ -78,7 +78,6 @@ MODULE m_fertetra
!Now that the fermi energy is guaranteed to be in the interval [lowBound,upperBound]
!We use the bisection method to find it
!-----------------------------------------------------------------------------------
DO WHILE(upperBound-lowBound.GT.1e-10)
ef = (lowBound+upperBound)/2.0
......@@ -105,6 +104,7 @@ MODULE m_fertetra
!---------------------------------------------------------------------
!Calculate final occupation and weights for individual kpts and bands
!---------------------------------------------------------------------
ef = (lowBound+upperBound)/2.0
dfermi = 0.0
DO jspin = 1, jspins
!-------------------------------------------------------
......
......@@ -59,7 +59,6 @@ MODULE m_types_gfinp
TYPE, EXTENDS(t_fleurinput_base):: t_gfinp
!General logical switches
LOGICAL :: l_sphavg = .TRUE.
LOGICAL :: l_mperp = .FALSE.
REAL :: minCalcDistance=-1.0 !This distance has to be reached before green's functions are calculated
!Negative means it is evaluated at every iteration
......@@ -87,6 +86,8 @@ MODULE m_types_gfinp
PROCEDURE :: add => add_gfelem
PROCEDURE :: uniqueElements => uniqueElements_gfinp
PROCEDURE :: eMesh => eMesh_gfinp
PROCEDURE :: checkRadial => checkRadial_gfinp
PROCEDURE :: checkSphavg => checkSphavg_gfinp
PROCEDURE :: addNearestNeighbours => addNearestNeighbours_gfelem
END TYPE t_gfinp
......@@ -106,7 +107,6 @@ CONTAINS
ELSE
rank = 0
END IF
CALL mpi_bc(this%l_sphavg,rank,mpi_comm)
CALL mpi_bc(this%l_mperp,rank,mpi_comm)
CALL mpi_bc(this%minCalcDistance,rank,mpi_comm)
CALL mpi_bc(this%n,rank,mpi_comm)
......@@ -169,7 +169,7 @@ CONTAINS
REAL :: fixedCutoff
CHARACTER(len=200) :: xPathA,xPathS,label,cutoffArg,str
CHARACTER(len=1),PARAMETER :: spdf(0:3) = ['s','p','d','f']
LOGICAL :: l_gfinfo_given,l_fixedCutoffset
LOGICAL :: l_gfinfo_given,l_fixedCutoffset,l_sphavg
LOGICAL :: lp_calc(0:3,0:3)
xPathA = '/fleurInput/calculationSetup/greensFunction'
......@@ -177,7 +177,6 @@ CONTAINS
l_gfinfo_given = numberNodes.EQ.1
IF (l_gfinfo_given) THEN
this%l_sphavg=evaluateFirstBoolOnly(xml%GetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@l_sphavg'))
this%l_mperp=evaluateFirstBoolOnly(xml%GetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@l_mperp'))
this%minCalcDistance=evaluateFirstOnly(xml%GetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@minCalcDistance'))
......@@ -280,6 +279,7 @@ CONTAINS
WRITE(xPathA,'(a,i0,a)') TRIM(ADJUSTL(xPathS))//'/greensfCalculation[',i,']'
label = TRIM(ADJUSTL(xml%GetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@label')))
l_sphavg = evaluateFirstBoolOnly(xml%GetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@l_sphavg'))
iContour = this%find_contour(TRIM(ADJUSTL(label)))
cutoffArg = TRIM(ADJUSTL(xml%GetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@kkintgrCutoff')))
nshells = evaluateFirstIntOnly(xml%GetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@nshells'))
......@@ -310,7 +310,7 @@ CONTAINS
READ(str,'(4l2)') (lp_calc(lp,l),lp=0,3)
DO lp = 0,lmaxU_const
IF(.NOT.lp_calc(lp,l)) CYCLE
i_gf = this%add(l,itype,iContour,.TRUE.,lp=lp,l_fixedCutoffset=l_fixedCutoffset,&
i_gf = this%add(l,itype,iContour,l_sphavg,lp=lp,l_fixedCutoffset=l_fixedCutoffset,&
fixedCutoff=fixedCutoff,l_inter=(nshells/=0))
ENDDO
ENDDO
......@@ -323,7 +323,7 @@ CONTAINS
DO l = 0,lmaxU_const
lp_calc(l,l) = evaluateFirstBoolOnly(xml%GetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@'//spdf(l)))
IF(.NOT.lp_calc(l,l)) CYCLE
i_gf = this%add(l,itype,iContour,.TRUE.,l_fixedCutoffset=l_fixedCutoffset,&
i_gf = this%add(l,itype,iContour,l_sphavg,l_fixedCutoffset=l_fixedCutoffset,&
fixedCutoff=fixedCutoff,l_inter=(nshells/=0))
ENDDO
ENDIF
......@@ -331,11 +331,11 @@ CONTAINS
!Find the reference element
IF(refL /= -1) THEN
!Find the element
refGF = this%find(refL,itype,iContour,.TRUE.)
refGF = this%find(refL,itype,iContour,l_sphavg)
DO l = 0,lmaxU_const
DO lp = 0,lmaxU_const
IF(.NOT.lp_calc(lp,l)) CYCLE
i_gf = this%find(l,itype,iContour,.TRUE.,lp=lp)
i_gf = this%find(l,itype,iContour,l_sphavg,lp=lp)
IF(i_gf==refGF) CYCLE
this%elem(i_gf)%refCutoff = refGF
ENDDO
......@@ -357,6 +357,7 @@ CONTAINS
fixedCutoff = evaluateFirstOnly(TRIM(ADJUSTL(cutoffArg)))
l_fixedCutoffset = .TRUE.
ENDIF
!Hubbard 1 GF has to be spherically averaged
i_gf = this%add(l,itype,iContour,.TRUE.,l_fixedCutoffset=l_fixedCutoffset,&
fixedCutoff=fixedCutoff)
n_hia = n_hia + 1
......@@ -394,7 +395,8 @@ CONTAINS
READ(str,'(4l2)') (lp_calc(lp,l),lp=0,3)
DO lp = 0,lmaxU_const
IF(.NOT.lp_calc(lp,l)) CYCLE
i_gf = this%add(l,itype,iContour,.TRUE.,lp=lp,l_fixedCutoffset=l_fixedCutoffset,&
!Torgue GF has to have radial dependence
i_gf = this%add(l,itype,iContour,.FALSE.,lp=lp,l_fixedCutoffset=l_fixedCutoffset,&
fixedCutoff=fixedCutoff)
this%numTorgueElems(itype) = this%numTorgueElems(itype) + 1
this%torgueElem(itype,this%numTorgueElems(itype)) = i_gf
......@@ -404,11 +406,11 @@ CONTAINS
!Find the reference element
IF(refL /= -1) THEN
!Find the element
refGF = this%find(refL,itype,iContour,.TRUE.)
refGF = this%find(refL,itype,iContour,.FALSE.)
DO l = 0,lmaxU_const
DO lp = 0,lmaxU_const
IF(.NOT.lp_calc(lp,l)) CYCLE
i_gf = this%find(l,itype,iContour,.TRUE.,lp=lp)
i_gf = this%find(l,itype,iContour,.FALSE.,lp=lp)
IF(i_gf==refGF) CYCLE
this%elem(i_gf)%refCutoff = refGF
ENDDO
......@@ -452,7 +454,7 @@ CONTAINS
LOGICAL, INTENT(IN) :: l_write
INTEGER :: i_gf,l,lp,atomType,atomTypep,iContour
LOGICAL :: l_inter,l_offd
LOGICAL :: l_inter,l_offd,l_sphavg,l_interAvg,l_offdAvg
INTEGER :: hiaElem(atoms%n_hia)
IF(this%n==0) RETURN !Nothing to do here
......@@ -488,13 +490,22 @@ CONTAINS
l_inter = .FALSE.
l_offd = .FALSE.
l_interAvg = .FALSE.
l_offdAvg = .FALSE.
DO i_gf = 1, this%n
l = this%elem(i_gf)%l
lp = this%elem(i_gf)%lp
atomType = this%elem(i_gf)%atomType
atomTypep = this%elem(i_gf)%atomTypep
IF(atomType.NE.atomTypep) l_inter = .TRUE.
IF(l.NE.lp) l_offd = .TRUE.
l_sphavg = this%elem(i_gf)%l_sphavg
IF(atomType.NE.atomTypep) THEN
l_inter = .TRUE.
IF(l_sphavg) l_interAvg = .TRUE.
ENDIF
IF(l.NE.lp) THEN
l_offd = .TRUE.
IF(l_sphavg) l_offdAvg = .TRUE.
ENDIF
ENDDO
......@@ -502,7 +513,7 @@ CONTAINS
IF(sym%nop>1) THEN
CALL juDFT_error("Symmetries and intersite Green's Function not implemented",&
calledby="init_gfinp")
ELSE IF(this%l_sphavg) THEN
ELSE IF(l_interAvg) THEN
CALL juDFT_error("Spherical average and intersite Green's Function not implemented",&
calledby="init_gfinp")
ENDIF
......@@ -512,7 +523,7 @@ CONTAINS
IF(sym%nop>1) THEN
CALL juDFT_warn("Symmetries and l-offdiagonal Green's Function not correctly implemented",&
calledby="init_gfinp")
ELSE IF(this%l_sphavg) THEN
ELSE IF(l_offdAvg) THEN
CALL juDFT_error("Spherical average and l-offdiagonal Green's Function not implemented",&
calledby="init_gfinp")
ENDIF
......@@ -566,11 +577,6 @@ CONTAINS
IF(iUnique == i_gf) uniqueElements = uniqueElements +1
ENDDO
IF(uniqueElements==0 .AND. maxGF/=0) THEN
CALL juDFT_error("No unique GF elements",hint="This is a bug in FLEUR please report",&
calledby="uniqueElements_gfinp")
ENDIF
IF(PRESENT(indUnique)) THEN
IF(.NOT.PRESENT(ind)) CALL juDFT_error("ind and indUnique have to be provided at the same time",&
calledby="uniqueElements_gfinp")
......@@ -690,7 +696,8 @@ CONTAINS
IF(ABS(dist(natomp)-minDist).LT.1e-12) THEN
!Add the element to the gfinp%elem array
this%elem(this%n)%atomDiff = atoms%taual(:,natomp) - atoms%taual(:,refAtom)
i_gf = this%add(l,refAtom,iContour,.TRUE.,lp=lp,nTypep=natomp,l_fixedCutoffset=l_fixedCutoffset,&
!cannot be spherically averaged
i_gf = this%add(l,refAtom,iContour,.FALSE.,lp=lp,nTypep=natomp,l_fixedCutoffset=l_fixedCutoffset,&
fixedCutoff=fixedCutoff)
this%elem(i_gf)%refCutoff = refCutoff
dist(natomp) = 9e99 !Eliminate from the list
......@@ -840,4 +847,32 @@ CONTAINS
END SUBROUTINE eMesh_gfinp
PURE LOGICAL FUNCTION checkRadial_gfinp(this)
!Check if there are any elements with radial dependence
CLASS(t_gfinp), INTENT(IN) :: this
INTEGER :: i_gf
checkRadial_gfinp = .FALSE.
DO i_gf = 1, this%n
IF(.NOT.this%elem(i_gf)%l_sphavg) checkRadial_gfinp = .TRUE.
ENDDO
END FUNCTION checkRadial_gfinp
PURE LOGICAL FUNCTION checkSphavg_gfinp(this)
!Check if there are any elements with spherical averaging
CLASS(t_gfinp), INTENT(IN) :: this
INTEGER :: i_gf
checkSphavg_gfinp = .FALSE.
DO i_gf = 1, this%n
IF(this%elem(i_gf)%l_sphavg) checkSphavg_gfinp = .TRUE.
ENDDO
END FUNCTION checkSphavg_gfinp
END MODULE m_types_gfinp
......@@ -17,13 +17,13 @@ MODULE m_excSplitting
REAL, INTENT(IN) :: ef
INTEGER :: i_gf,i_elem,indUnique,ispin,kkcut,m,l,lp,atomType,atomTypep,ie
LOGICAL :: l_sphavg
REAL :: excSplit,del
REAL, ALLOCATABLE :: eMesh(:), imag(:)
REAL, ALLOCATABLE :: intCOM(:,:), intNorm(:,:)
CHARACTER(LEN=20) :: attributes(4)
IF(.NOT.gfinp%l_sphavg) RETURN
IF(.NOT.gfinp%checkSphavg()) RETURN
CALL gfinp%eMesh(ef,del=del,eMesh=eMesh)
......@@ -43,13 +43,14 @@ MODULE m_excSplitting
lp = gfinp%elem(i_gf)%lp
atomType = gfinp%elem(i_gf)%atomType
atomTypep = gfinp%elem(i_gf)%atomTypep
l_sphavg = gfinp%elem(i_gf)%l_sphavg
!Only onsite exchange splitting
IF(l /= lp) CYCLE
IF(atomType /= atomTypep) CYCLE
IF(.NOT.gfinp%l_sphavg) CYCLE
IF(.NOT.l_sphavg) CYCLE
i_elem = gfinp%uniqueElements(ind=i_gf,indUnique=indUnique)
i_elem = gfinp%uniqueElements(ind=i_gf,l_sphavg=l_sphavg,indUnique=indUnique)
IF(i_gf /= indUnique) CYCLE
!-------------------------------------------------
......@@ -64,7 +65,7 @@ MODULE m_excSplitting
intCOM = 0.0
intNorm = 0.0
DO m = -l, l
imag = greensfImagPart%applyCutoff(i_elem,i_gf,m,m,ispin)
imag = greensfImagPart%applyCutoff(i_elem,i_gf,m,m,ispin,l_sphavg)
intCOM(:,ispin) = intCOM(:,ispin) + eMesh*imag
intNorm(:,ispin) = intNorm(:,ispin) + imag
ENDDO
......
......@@ -34,6 +34,7 @@ MODULE m_greensfBZint
INTEGER :: spin1,spin2,ispin,spin_start,spin_end
COMPLEX :: phase
REAL :: atomFactor
LOGICAL :: l_sphavg
COMPLEX, ALLOCATABLE :: im(:,:,:,:,:)
spin_start = MERGE(1,jspin,gfinp%l_mperp)
......@@ -44,9 +45,6 @@ MODULE m_greensfBZint
CALL timestart("Green's Function: Brillouin-Zone-Integration")
ALLOCATE(im(-lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const,nBands,&
MERGE(1,4,gfinp%l_sphavg),spin_start:spin_end),source=cmplx_0)
DO i_gf = 1, gfinp%n
!Get the information about the current element
......@@ -54,13 +52,17 @@ MODULE m_greensfBZint
lp = gfinp%elem(i_gf)%lp
atomType = gfinp%elem(i_gf)%atomType
atomTypep = gfinp%elem(i_gf)%atomTypep
l_sphavg = gfinp%elem(i_gf)%l_sphavg
atomFactor = MERGE(1.0,1.0/atoms%neq(atomType),l.NE.lp)
atomFactor = MERGE(1.0,atomFactor,atomType.NE.atomTypep)
i_elem = gfinp%uniqueElements(ind=i_gf,indUnique=indUnique)
i_elem = gfinp%uniqueElements(ind=i_gf,l_sphavg=l_sphavg,indUnique=indUnique)
IF(i_gf/=indUnique) CYCLE
ALLOCATE(im(-lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const,nBands,&
MERGE(1,4,l_sphavg),spin_start:spin_end),source=cmplx_0)
natom_start = SUM(atoms%neq(:atomType-1)) + 1
natom_end = MERGE(SUM(atoms%neq(:atomType-1)) + 1,SUM(atoms%neq(:atomType)),l.NE.lp)
!Loop over equivalent atoms
......@@ -99,20 +101,20 @@ MODULE m_greensfBZint
!Spin diagonal elements
IF(spin1==spin2) THEN
CALL greensfSpinDiag(nBands,l,lp,natom,natomp,atomType,atomTypep,spin1,&
gfinp%l_sphavg,atoms,usdus,eigVecCoeffs,im(:,:,:,:,ispin))
l_sphavg,atoms,usdus,eigVecCoeffs,im(:,:,:,:,ispin))
ELSE
!Spin offdiagonal elements
CALL greensfSpinOffDiag(nBands,l,lp,natom,natomp,atomType,atomTypep,spin1,spin2,&
gfinp%l_sphavg,atoms,denCoeffsOffdiag,eigVecCoeffs,im(:,:,:,:,ispin))
l_sphavg,atoms,denCoeffsOffdiag,eigVecCoeffs,im(:,:,:,:,ispin))
ENDIF
CALL greensfSym(ikpt_i,i_elem,natom,l,natom.EQ.natomp.AND.l.EQ.lp,gfinp%l_sphavg,&
CALL greensfSym(ikpt_i,i_elem,natom,l,natom.EQ.natomp.AND.l.EQ.lp,l_sphavg,&
ispin,sym,atomFactor,phase,im(:,:,:,:,ispin),greensfBZintCoeffs)
ENDDO
ENDDO !natomp
ENDDO !natom
DEALLOCATE(im)
ENDDO !i_gf
CALL timestop("Green's Function: Brillouin-Zone-Integration")
......
......@@ -29,7 +29,7 @@ MODULE m_greensfCalcImagPart
INTEGER :: ikpt_i,ikpt,nBands,jsp,i_gf
INTEGER :: l,lp,m,mp,iBand,ie,j,eGrid_start,eGrid_end
INTEGER :: indUnique,i_elem
LOGICAL :: l_zero
LOGICAL :: l_zero,l_sphavg
REAL :: del,eb,wtkpt
COMPLEX :: fac,weight
INTEGER, ALLOCATABLE :: ev_list(:)
......@@ -86,18 +86,19 @@ MODULE m_greensfCalcImagPart
!Get the information about the current element
l = gfinp%elem(i_gf)%l
lp = gfinp%elem(i_gf)%lp
l_sphavg = gfinp%elem(i_gf)%l_sphavg
i_elem = gfinp%uniqueElements(ind=i_gf,indUnique=indUnique)
i_elem = gfinp%uniqueElements(ind=i_gf,l_sphavg=l_sphavg,indUnique=indUnique)
IF(i_gf/=indUnique) CYCLE
!$OMP parallel default(none) &
!$OMP shared(gfinp,input,greensfBZintCoeffs,greensfImagPart) &
!$OMP shared(i_elem,l,lp,ikpt_i,nBands,eMesh)&
!$OMP shared(input,greensfBZintCoeffs,greensfImagPart) &
!$OMP shared(i_elem,l,lp,ikpt_i,nBands,eMesh,l_sphavg)&
!$OMP shared(del,eb,eig,weights,indBound,fac,wtkpt,spin_ind) &
!$OMP private(ie,m,mp,iBand,j,eGrid_start,eGrid_end,weight,imag,imagReal,l_zero)
ALLOCATE(imag(SIZE(eMesh),MERGE(1,4,gfinp%l_sphavg)),source=cmplx_0)
ALLOCATE(imagReal(SIZE(eMesh),MERGE(1,4,gfinp%l_sphavg)),source=0.0)
ALLOCATE(imag(SIZE(eMesh),MERGE(1,4,l_sphavg)),source=cmplx_0)
ALLOCATE(imagReal(SIZE(eMesh),MERGE(1,4,l_sphavg)),source=0.0)
!$OMP do collapse(2)
DO mp = -lp, lp
DO m = -l, l
......@@ -132,7 +133,7 @@ MODULE m_greensfCalcImagPart
CASE DEFAULT
END SELECT
IF(gfinp%l_sphavg) THEN
IF(l_sphavg) THEN
imag(ie,1) = imag(ie,1) + weight * greensfBZintCoeffs%sphavg(iBand,m,mp,i_elem,ikpt_i,spin_ind)
ELSE
imag(ie,1) = imag(ie,1) + weight * greensfBZintCoeffs%uu(iBand,m,mp,i_elem,ikpt_i,spin_ind)
......@@ -143,7 +144,7 @@ MODULE m_greensfCalcImagPart
ELSE IF(eGrid_start==1 .AND. eGrid_end==SIZE(eMesh)) THEN!Here we always use the tetrahedron method
!We can only use the BLAS routine on the full array
IF(gfinp%l_sphavg) THEN
IF(l_sphavg) THEN
CALL CPP_BLAS_caxpy(SIZE(eMesh),greensfBZintCoeffs%sphavg(iBand,m,mp,i_elem,ikpt_i,spin_ind),&
weights(:,iBand),1,imag(:,1),1)
ELSE
......@@ -157,7 +158,7 @@ MODULE m_greensfCalcImagPart
weights(:,iBand),1,imag(:,4),1)
ENDIF
ELSE
IF(gfinp%l_sphavg) THEN
IF(l_sphavg) THEN
imag(eGrid_start:eGrid_end,1) = imag(eGrid_start:eGrid_end,1) + weights(eGrid_start:eGrid_end,iBand)&
* greensfBZintCoeffs%sphavg(iBand,m,mp,i_elem,ikpt_i,spin_ind)
ELSE
......@@ -174,7 +175,7 @@ MODULE m_greensfCalcImagPart
ENDDO!ib
imagReal = AIMAG(imag)
IF(gfinp%l_sphavg) THEN
IF(l_sphavg) THEN
CALL CPP_BLAS_saxpy(SIZE(eMesh),1.0,imagReal(:,1),1,greensfImagPart%sphavg(:,m,mp,i_elem,spin_ind),1)
ELSE
CALL CPP_BLAS_saxpy(SIZE(eMesh),1.0,imagReal(:,1),1,greensfImagPart%uu(:,m,mp,i_elem,spin_ind),1)
......
......@@ -47,7 +47,7 @@ MODULE m_greensfCalcRealPart
INTEGER :: spin_cut,contourShape
INTEGER :: i_gf_start,i_gf_end,spin_start,spin_end
INTEGER :: n_gf_task,extra
LOGICAL :: l_onsite,l_fixedCutoffset
LOGICAL :: l_onsite,l_fixedCutoffset,l_sphavg
REAL :: fac,del,eb,et,fixedCutoff
REAL, ALLOCATABLE :: eMesh(:),imag(:)
......@@ -65,11 +65,12 @@ MODULE m_greensfCalcRealPart
lp = g(i_gf)%elem%lp
nType = g(i_gf)%elem%atomType
nTypep = g(i_gf)%elem%atomTypep
l_sphavg = g(i_gf)%elem%l_sphavg
l_fixedCutoffset = g(i_gf)%elem%l_fixedCutoffset
fixedCutoff = g(i_gf)%elem%fixedCutoff
refCutoff = g(i_gf)%elem%refCutoff
i_elem = gfinp%uniqueElements(ind=i_gf,indUnique=indUnique)
i_elem = gfinp%uniqueElements(ind=i_gf,l_sphavg=l_sphavg,indUnique=indUnique)
IF(i_gf /= indUnique.AND..NOT.l_fixedCutoffset.AND.refCutoff==-1&
.AND..NOT.g(indUnique)%elem%l_fixedCutoffset.AND.g(indUnique)%elem%refCutoff==-1) THEN
......@@ -82,15 +83,15 @@ MODULE m_greensfCalcRealPart
!
!Check the integral over the fDOS to define a cutoff for the Kramer-Kronigs-Integration
!
IF(gfinp%l_sphavg) THEN
IF(l_sphavg) THEN
CALL kk_cutoff(greensfImagPart%sphavg(:,:,:,i_elem,:),noco,gfinp%l_mperp,l,input%jspins,&
eMesh,greensfImagPart%kkintgr_cutoff(i_gf,:,:),greensfImagPart%scalingFactor(i_elem,:))
eMesh,greensfImagPart%kkintgr_cutoff(i_gf,:,:),greensfImagPart%scalingFactorSphavg(i_elem,:))
ELSE
!Onsite element with radial dependence
CALL kk_cutoffRadial(greensfImagPart%uu(:,:,:,i_elem,:),greensfImagPart%ud(:,:,:,i_elem,:),&
greensfImagPart%du(:,:,:,i_elem,:),greensfImagPart%dd(:,:,:,i_elem,:),&
noco,atoms,vTot,enpara,fmpi,hub1inp,gfinp%l_mperp,l,nType,input,eMesh,&
greensfImagPart%kkintgr_cutoff(i_gf,:,:),greensfImagPart%scalingFactor(i_elem,:))
greensfImagPart%kkintgr_cutoff(i_gf,:,:),greensfImagPart%scalingFactorRadial(i_elem,:))
ENDIF
ELSE IF (l_fixedCutoffset) THEN
greensfImagPart%kkintgr_cutoff(i_gf,:,1) = 1
......@@ -113,13 +114,14 @@ MODULE m_greensfCalcRealPart
l_fixedCutoffset = g(i_gf)%elem%l_fixedCutoffset
fixedCutoff = g(i_gf)%elem%fixedCutoff
refCutoff = g(i_gf)%elem%refCutoff
i_elem = gfinp%uniqueElements(ind=i_gf,indUnique=indUnique)
l_sphavg = g(i_gf)%elem%l_sphavg
i_elem = gfinp%uniqueElements(ind=i_gf,l_sphavg=l_sphavg,indUnique=indUnique)
IF(refCutoff/=-1) THEN
!Overwrite cutoff with reference from other elements
greensfImagPart%kkintgr_cutoff(i_gf,:,:) = greensfImagPart%kkintgr_cutoff(refCutoff,:,:)
ENDIF
CALL greensfImagPart%scale(i_elem)
CALL greensfImagPart%scale(i_elem,l_sphavg)
ENDDO
CALL timestop("Green's Function: Integration Cutoff")
ENDIF
......@@ -179,9 +181,10 @@ MODULE m_greensfCalcRealPart
lp = g(i_gf)%elem%lp
nType = g(i_gf)%elem%atomType
nTypep = g(i_gf)%elem%atomTypep
l_sphavg = g(i_gf)%elem%l_sphavg
contourShape = gfinp%contour(g(i_gf)%elem%iContour)%shape
i_elem = gfinp%uniqueElements(ind=i_gf)
i_elem = gfinp%uniqueElements(ind=i_gf,l_sphavg=l_sphavg)
CALL timestart("Green's Function: Kramer-Kronigs-Integration")
DO jspin = spin_start, spin_end
......@@ -189,28 +192,28 @@ MODULE m_greensfCalcRealPart
DO m= -l,l
DO mp= -lp,lp
IF(greensfImagPart%checkEmpty(i_elem,m,mp,jspin)) THEN
IF(greensfImagPart%checkEmpty(i_elem,m,mp,jspin,l_sphavg)) THEN
CALL g(i_gf)%resetSingleElem(m,mp,jspin,ipm)
CYCLE
ENDIF
IF(gfinp%l_sphavg) THEN
imag = greensfImagPart%applyCutoff(i_elem,i_gf,m,mp,jspin)
IF(l_sphavg) THEN
imag = greensfImagPart%applyCutoff(i_elem,i_gf,m,mp,jspin,l_sphavg)
CALL kkintgr(imag,eMesh,g(i_gf)%contour%e,(ipm.EQ.2),&
g(i_gf)%gmmpMat(:,m,mp,jspin,ipm),int_method(contourShape))
ELSE
! In the case of radial dependence we perform the kramers-kronig-integration seperately for uu,dd,etc.
! We can do this because the radial functions are independent of E
imag = greensfImagPart%applyCutoff(i_elem,i_gf,m,mp,jspin,imat=1)
imag = greensfImagPart%applyCutoff(i_elem,i_gf,m,mp,jspin,l_sphavg,imat=1)
CALL kkintgr(imag,eMesh,g(i_gf)%contour%e,(ipm.EQ.2),&
g(i_gf)%uu(:,m,mp,jspin,ipm),int_method(contourShape))
imag = greensfImagPart%applyCutoff(i_elem,i_gf,m,mp,jspin,imat=2)
imag = greensfImagPart%applyCutoff(i_elem,i_gf,m,mp,jspin,l_sphavg,imat=2)
CALL kkintgr(imag,eMesh,g(i_gf)%contour%e,(ipm.EQ.2),&
g(i_gf)%dd(:,m,mp,jspin,ipm),int_method(contourShape))
imag = greensfImagPart%applyCutoff(i_elem,i_gf,m,mp,jspin,imat=3)
imag = greensfImagPart%applyCutoff(i_elem,i_gf,m,mp,jspin,l_sphavg,imat=3)
CALL kkintgr(imag,eMesh,g(i_gf)%contour%e,(ipm.EQ.2),&
g(i_gf)%ud(:,m,mp,jspin,ipm),int_method(contourShape))
imag = greensfImagPart%applyCutoff(i_elem,i_gf,m,mp,jspin,imat=4)
imag = greensfImagPart%applyCutoff(i_elem,i_gf,m,mp,jspin,l_sphavg,imat=4)
CALL kkintgr(imag,eMesh,g(i_gf)%contour%e,(ipm.EQ.2),&
g(i_gf)%du(:,m,mp,jspin,ipm),int_method(contourShape))
ENDIF
......
......@@ -39,6 +39,7 @@ MODULE m_greensfPostProcess
INTEGER i_gf,nType,l,lp,atomType,atomTypep,i_elem,indUnique,jspin,ierr
COMPLEX mmpmat(-lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const,gfinp%n,3)
LOGICAL l_sphavg
REAL :: torgue(3)
REAL, ALLOCATABLE :: u(:,:,:,:,:,:),udot(:,:,:,:,:,:)
......@@ -71,7 +72,7 @@ MODULE m_greensfPostProcess
CALL excSplitting(gfinp,input,greensfImagPart,results%ef)
IF(.NOT.gfinp%l_sphavg) THEN
IF(gfinp%checkRadial()) THEN
CALL timestart("Green's Function: Radial Functions")
ALLOCATE (f(atoms%jmtd,2,0:atoms%lmaxd),source=0.0)
ALLOCATE (g(atoms%jmtd,2,0:atoms%lmaxd),source=0.0)
......@@ -88,8 +89,10 @@ MODULE m_greensfPostProcess
lp = gfinp%elem(i_gf)%lp
atomType = gfinp%elem(i_gf)%atomType
atomTypep = gfinp%elem(i_gf)%atomTypep
l_sphavg = gfinp%elem(i_gf)%l_sphavg
IF(l_sphavg) CYCLE
i_elem = gfinp%uniqueElements(ind=i_gf,indUnique=indUnique)
i_elem = gfinp%uniqueElements(ind=i_gf,l_sphavg=l_sphavg,indUnique=indUnique)
IF(i_gf/=indUnique) THEN
u(:,:,:,:,:,i_gf) = u(:,:,:,:,:,indUnique)
......@@ -136,9 +139,10 @@ MODULE m_greensfPostProcess
lp = greensFunction(i_gf)%elem%lp
atomType = greensFunction(i_gf)%elem%atomType
atomTypep = greensFunction(i_gf)%elem%atomTypep
l_sphavg = gfinp%elem(i_gf)%l_sphavg
IF(l.NE.lp) CYCLE
IF(atomType.NE.atomTypep) CYCLE
IF(gfinp%l_sphavg) THEN
IF(l_sphavg) THEN
CALL occmtx(greensFunction(i_gf),gfinp,input,mmpmat(:,:,i_gf,:),l_write=.TRUE.,check=.TRUE.)
ELSE IF(.NOT.gfinp%l_mperp) THEN
CALL occmtx(greensFunction(i_gf),gfinp,input,mmpmat(:,:,i_gf,:),&
......
......@@ -70,7 +70,6 @@ MODULE m_greensf_io
CALL h5gcreate_f(fileID, '/general', generalGroupID, hdfError)
CALL io_write_attint0(generalGroupID,'spins',input%jspins)
CALL io_write_attreal0(generalGroupID,'FermiEnergy',eFermiPrev)
CALL io_write_attlog0(generalGroupID,'sphavg',gfinp%l_sphavg)
CALL io_write_attlog0(generalGroupID,'mperp',gfinp%l_mperp)
CALL h5gclose_f(generalGroupID, hdfError)
......@@ -126,6 +125,7 @@ MODULE m_greensf_io
INTEGER :: i_elem,n_elem
INTEGER(HSIZE_T) :: dims(7)
REAL :: trc(MERGE(3,input%jspins,gfinp%l_mperp))
LOGICAL :: l_sphavg
SELECT CASE(archiveType)
......@@ -168,6 +168,7 @@ MODULE m_greensf_io
lp = greensf(i_elem)%elem%lp
atomType = greensf(i_elem)%elem%atomType
atomTypep = greensf(i_elem)%elem%atomTypep
l_sphavg = greensf(i_elem)%elem%l_sphavg
iContour = greensf(i_elem)%elem%iContour
CALL h5gcreate_f(elementsGroupID, elementName, currentelementGroupID, hdfError)
......@@ -175,6 +176,7 @@ MODULE m_greensf_io
CALL io_write_attint0(currentelementGroupID,"lp",lp)
CALL io_write_attint0(currentelementGroupID,"atomType",atomType)
CALL io_write_attint0(currentelementGroupID,"atomTypep",atomTypep)
CALL io_write_attlog0(currentelementGroupID,'l_sphavg',l_sphavg)
CALL io_write_attint0(currentelementGroupID,'nz',greensf(i_elem)%contour%nz)
......@@ -233,7 +235,7 @@ MODULE m_greensf_io
CALL h5dclose_f(mmpmatSetID, hdfError)
!Spherically averaged greensfData
IF(gfinp%l_sphavg) THEN
IF(l_sphavg) THEN
dims(:6)=[2,greensf(i_elem)%contour%nz,2*lmaxU_Const+1,2*lmaxU_Const+1,jspinsOut,2]
dimsInt=dims
......@@ -243,7 +245,7 @@ MODULE m_greensf_io
CALL io_write_complex5(sphavgDataSetID,[-1,1,1,1,1,1],dimsInt(:6),greensf(i_elem)%gmmpmat)
CALL h5dclose_f(sphavgDataSetID, hdfError)
ELSE IF(.NOT.gfinp%l_sphavg.AND.archiveType.NE.GREENSF_HUBBARD_CONST) THEN
ELSE IF(.NOT.l_sphavg.AND.archiveType.NE.GREENSF_HUBBARD_CONST) THEN
!uu
dims(:6)=[2,greensf(i_elem)%contour%nz,2*lmaxU_Const+1,2*lmaxU_Const+1,jspinsOut,2]
......
......@@ -234,8 +234,8 @@ CONTAINS
WRITE (fileNum, 195) input%ldauLinMix, input%ldauMixParam, input%ldauSpinf
IF(l_gfOpt) THEN
205 FORMAT(' <greensFunction l_sphavg="', l1, '" l_mperp="', l1'">')