Skip to content

GitLab

fleur
fleur
Commit dea084ba authored by Matthias Redies's avatar Matthias Redies
Browse files
Options

Merge branch 'develop' of iffgit.fz-juelich.de:fleur/fleur into develop

parents 97362865 cf55ed4e
Expand all Hide whitespace changes
Inline Side-by-side
Showing with 435 additions and 801 deletions
fleurinput/types_wannier.f90
View file @ dea084ba
......@@ -15,31 +15,33 @@ MODULE m_types_wannier
TYPE,EXTENDS(t_fleurinput_base):: t_wann
!New parameters not handled correctly yet...
LOGICAL :: l_socmatvec=.FALSE.
INTEGER :: socmatvecfmt=1
LOGICAL :: l_socmatvecrs=.FALSE.
INTEGER :: socmatvecrsfmt=1
LOGICAL :: l_mmn0_unf_to_spn_unf=.FALSE.
LOGICAL :: l_mmn0_to_spn_unf=.FALSE.
LOGICAL :: l_mmn0_to_spn=.FALSE.
LOGICAL :: l_mmn0_to_spn2=.FALSE.
LOGICAL :: l_mmn0_unf_to_spn=.FALSE.
LOGICAL :: l_perpmag_unf_to_tor_unf=.FALSE.
LOGICAL :: l_perpmag_to_tor_unf=.FALSE.
LOGICAL :: l_perpmag_to_tor=.FALSE.
LOGICAL :: l_perpmag_unf_to_tor=.FALSE.
LOGICAL :: l_hsomtxvec_unf_to_lmpzsoc_unf=.FALSE.
LOGICAL :: l_hsomtxvec_to_lmpzsoc_unf=.FALSE.
LOGICAL :: l_hsomtxvec_to_lmpzsoc=.FALSE.
LOGICAL :: l_hsomtxvec_unf_to_lmpzsoc=.FALSE.
LOGICAL :: l_hsomtx_unf_to_hsoc_unf=.FALSE.
LOGICAL :: l_hsomtx_to_hsoc_unf=.FALSE.
LOGICAL :: l_hsomtx_to_hsoc=.FALSE.
LOGICAL :: l_hsomtx_unf_to_hsoc=.FALSE.
INTEGER :: perpmagl
LOGICAL :: l_perpmagatlres=.FALSE.
LOGICAL :: l_mmn0_to_spn_unf=.FALSE.
LOGICAL :: l_mmn0_to_spn=.FALSE.
LOGICAL :: l_mmn0_to_spn2=.FALSE.
LOGICAL :: l_mmn0_unf_to_spn=.FALSE.
LOGICAL :: l_perpmag_unf_to_tor_unf=.FALSE.
LOGICAL :: l_perpmag_to_tor_unf=.FALSE.
LOGICAL :: l_perpmag_to_tor=.FALSE.
LOGICAL :: l_perpmag_unf_to_tor=.FALSE.
LOGICAL :: l_hsomtxvec_unf_to_lmpzsoc_unf=.FALSE.
LOGICAL :: l_hsomtxvec_to_lmpzsoc_unf=.FALSE.
LOGICAL :: l_hsomtxvec_to_lmpzsoc=.FALSE.
LOGICAL :: l_hsomtxvec_unf_to_lmpzsoc=.FALSE.
LOGICAL :: l_hsomtx_unf_to_hsoc_unf=.FALSE.
LOGICAL :: l_hsomtx_to_hsoc_unf=.FALSE.
LOGICAL :: l_hsomtx_to_hsoc=.FALSE.
LOGICAL :: l_hsomtx_unf_to_hsoc=.FALSE.
INTEGER :: perpmagl
LOGICAL :: l_perpmagatlres=.FALSE.
INTEGER :: wan90version =3
INTEGER :: oc_num_orbs =0
INTEGER, ALLOCATABLE :: oc_orbs(:)
LOGICAL :: l_unformatted =.FALSE.
INTEGER :: wan90version =3
INTEGER :: oc_num_orbs =0
INTEGER, ALLOCATABLE :: oc_orbs(:)
LOGICAL :: l_unformatted =.FALSE.
LOGICAL :: l_oc_f=.FALSE.
LOGICAL :: l_ndegen=.FALSE.
LOGICAL :: l_orbitalmom=.FALSE.
......@@ -50,7 +52,9 @@ MODULE m_types_wannier
LOGICAL :: l_perturb=.FALSE.
LOGICAL :: l_nedrho=.FALSE.
LOGICAL :: l_anglmomrs=.FALSE.
INTEGER :: anglmomrsfmt=1
LOGICAL :: l_anglmom=.FALSE.
INTEGER :: anglmomfmt=1
LOGICAL :: l_spindisp=.FALSE.
LOGICAL :: l_spindisprs=.FALSE.
LOGICAL :: l_socspicom=.FALSE.
......@@ -58,17 +62,25 @@ MODULE m_types_wannier
LOGICAL :: l_offdiposoprs=.FALSE.
LOGICAL :: l_offdiposop=.FALSE.
LOGICAL :: l_torque=.FALSE.
INTEGER :: torquefmt=1
LOGICAL :: l_torquers=.FALSE.
INTEGER :: torquersfmt=1
LOGICAL :: l_atomlist=.FALSE.
INTEGER :: atomlist_num=0
INTEGER, ALLOCATABLE :: atomlist(:)
LOGICAL :: l_berry=.FALSE.
LOGICAL :: l_perpmagrs=.FALSE.
INTEGER :: perpmagrsfmt=1
LOGICAL :: l_perpmag=.FALSE.
INTEGER :: perpmagfmt=1
LOGICAL :: l_perpmagat=.FALSE.
INTEGER :: perpmagatfmt=1
LOGICAL :: l_perpmagatrs=.FALSE.
INTEGER :: perpmagatrsfmt=1
LOGICAL :: l_socmatrs=.FALSE.
INTEGER :: socmatrsfmt=1
LOGICAL :: l_socmat=.FALSE.
INTEGER :: socmatfmt=1
LOGICAL :: l_soctomom=.FALSE.
LOGICAL :: l_kptsreduc2=.FALSE.
LOGICAL :: l_nablapaulirs=.FALSE.
......@@ -81,13 +93,16 @@ MODULE m_types_wannier
LOGICAL :: l_nabla=.FALSE.
LOGICAL :: l_socodi=.FALSE.
LOGICAL :: l_pauli=.FALSE.
INTEGER :: paulifmt=1
LOGICAL :: l_pauliat=.FALSE.
INTEGER :: pauliatfmt=1
LOGICAL :: l_potmat=.FALSE.
LOGICAL :: l_projgen=.FALSE.
LOGICAL :: l_plot_symm=.FALSE.
LOGICAL :: l_socmmn0=.FALSE.
LOGICAL :: l_bzsym=.FALSE.
LOGICAL :: l_hopping=.FALSE.
INTEGER :: hoppingfmt=1
LOGICAL :: l_kptsreduc=.FALSE.
LOGICAL :: l_prepwan90=.FALSE.
LOGICAL :: l_plot_umdat=.FALSE.
......@@ -95,7 +110,9 @@ MODULE m_types_wannier
LOGICAL :: l_bynumber=.FALSE.
LOGICAL :: l_stopopt=.FALSE.
LOGICAL :: l_matrixmmn=.FALSE.
INTEGER :: matrixmmnfmt=1
LOGICAL :: l_matrixamn=.FALSE.
INTEGER :: matrixamnfmt=1
LOGICAL :: l_projmethod=.FALSE.
LOGICAL :: l_wannierize=.FALSE.
LOGICAL :: l_plotw90=.FALSE.
......@@ -111,7 +128,9 @@ MODULE m_types_wannier
LOGICAL :: l_dipole2=.FALSE.
LOGICAL :: l_dipole3=.FALSE.
LOGICAL :: l_mmn0=.FALSE.
INTEGER :: mmn0fmt=1
LOGICAL :: l_mmn0at=.FALSE.
INTEGER :: mmn0atfmt=1
LOGICAL :: l_manyfiles=.FALSE.
LOGICAL :: l_collectmanyfiles=.FALSE.
LOGICAL :: l_ldauwan=.FALSE.
......@@ -123,7 +142,9 @@ MODULE m_types_wannier
LOGICAL :: l_finishgwf=.FALSE.
LOGICAL :: l_skipkov=.FALSE.
LOGICAL :: l_matrixuHu=.FALSE.
INTEGER :: matrixuHufmt=1
LOGICAL :: l_matrixuHu_dmi=.FALSE.
INTEGER :: matrixuHudmifmt=1
INTEGER :: ikptstart=1
INTEGER :: band_min(1:2)=-1
INTEGER :: band_max(1:2)=-1
......@@ -171,7 +192,27 @@ CONTAINS
ELSE
rank=0
END IF
CALL mpi_bc(this%socmatvecfmt,rank,mpi_comm)
CALL mpi_bc(this%socmatvecrsfmt,rank,mpi_comm)
CALL mpi_bc(this%anglmomrsfmt,rank,mpi_comm)
CALL mpi_bc(this%anglmomfmt,rank,mpi_comm)
CALL mpi_bc(this%torquefmt,rank,mpi_comm)
CALL mpi_bc(this%torquersfmt,rank,mpi_comm)
CALL mpi_bc(this%perpmagrsfmt,rank,mpi_comm)
CALL mpi_bc(this%perpmagfmt,rank,mpi_comm)
CALL mpi_bc(this%perpmagatfmt,rank,mpi_comm)
CALL mpi_bc(this%perpmagatrsfmt,rank,mpi_comm)
CALL mpi_bc(this%socmatrsfmt,rank,mpi_comm)
CALL mpi_bc(this%socmatfmt,rank,mpi_comm)
CALL mpi_bc(this%paulifmt,rank,mpi_comm)
CALL mpi_bc(this%pauliatfmt,rank,mpi_comm)
CALL mpi_bc(this%hoppingfmt,rank,mpi_comm)
CALL mpi_bc(this%matrixmmnfmt,rank,mpi_comm)
CALL mpi_bc(this%matrixamnfmt,rank,mpi_comm)
CALL mpi_bc(this%mmn0fmt,rank,mpi_comm)
CALL mpi_bc(this%mmn0atfmt,rank,mpi_comm)
CALL mpi_bc(this%matrixuHufmt,rank,mpi_comm)
CALL mpi_bc(this%matrixuHudmifmt,rank,mpi_comm)
CALL mpi_bc(this%wan90version ,rank,mpi_comm)
CALL mpi_bc(this%oc_num_orbs ,rank,mpi_comm)
CALL mpi_bc(this%oc_orbs,rank,mpi_comm)
......@@ -330,9 +371,14 @@ CONTAINS
CHARACTER(len=100):: xPathA
CHARACTER(len=255):: valueString
CHARACTER(len=30):: jobname
CHARACTER(len=30):: param
INTEGER :: parampos
INTEGER :: stat
INTEGER :: numberNodes,i,n,numtokens
LOGICAL,ALLOCATABLE:: wannAtomList(:)
LOGICAL :: l_param
REAL :: version_real
xPathA = '/fleurInput/output/wannier'
numberNodes = xml%getNumberOfNodes(xPathA)
......@@ -383,6 +429,15 @@ CONTAINS
DO i = 1, numTokens
this%jobList(i) = xml%popFirstStringToken(valueString)
IF(this%jobList(i)(1:1).EQ.'!')cycle
parampos=index(this%jobList(i),'=')
if(parampos.gt.1)then
jobname=this%jobList(i)(1:parampos-1)
param=this%jobList(i)(parampos+1:)
l_param=.true.
else
jobname=this%jobList(i)
l_param=.false.
endif
IF(this%jobList(i).EQ.'socmat')THEN
this%l_socmat=.TRUE.
ELSEIF(this%jobList(i).EQ.'socmatvec')THEN
......@@ -458,8 +513,30 @@ CONTAINS
this%l_wann_plot=.TRUE.
ELSEIF(this%jobList(i).EQ.'bynumber')THEN
this%l_bynumber=.TRUE.
ELSEIF(this%jobList(i).EQ.'matrixmmn')THEN
this%l_matrixmmn=.TRUE.
ELSEIF(jobname.EQ.'matrixmmn')THEN
this%l_matrixmmn=.TRUE.
if(l_param)then
read(param,*,iostat=stat) this%matrixmmnfmt
if(stat/=0)then
CALL juDFT_error("problem with jobparam=",calledby="wann_read_inp")
endif
endif
ELSEIF(jobname.EQ.'perpmag')THEN
this%l_perpmag=.TRUE.
if(l_param)then
read(param,*,iostat=stat) this%perpmagfmt
if(stat/=0)then
CALL juDFT_error("problem with jobparam=",calledby="wann_read_inp")
endif
endif
ELSEIF(jobname.EQ.'perpmagrs')THEN
this%l_perpmagrs=.TRUE.
if(l_param)then
read(param,*,iostat=stat) this%perpmagrsfmt
if(stat/=0)then
CALL juDFT_error("problem with jobparam=",calledby="wann_read_inp")
endif
endif
ELSEIF(this%jobList(i).EQ.'projmethod')THEN
this%l_projmethod=.TRUE.
ELSEIF(this%jobList(i).EQ.'matrixamn')THEN
......@@ -520,20 +597,28 @@ CONTAINS
this%l_matrixuHu=.TRUE.
ELSEIF(this%jobList(i).EQ.'matrixuhu-dmi') THEN
this%l_matrixuHu_dmi=.TRUE.
!Not done
! ELSEIF(this%jobList(i).EQ.'wan90version')THEN
! backspace(916)
! read(916,*,iostat=ios)task,version_real
! if (ios /= 0) CALL judft_error("error reading wan90version", calledby="wann_read_inp")
! if(abs(version_real-1.1).lt.1.e-9)THEN
! this%wan90version=1
! ELSEIF(abs(version_real-1.2).lt.1.e-9)THEN
! this%wan90version=2
! ELSEIF(abs(version_real-2.0).lt.1.e-9)THEN
! this%wan90version=3
! ELSE
! CALL judft_error ("chosen w90 version unknown", calledby="wann_read_inp")
! endif
ELSEIF(jobname.EQ.'wan90version')THEN
if(l_param)then
read(param,*,iostat=stat) version_real
if(stat/=0)then
write(*,*)"problem with jobparam=",param
CALL juDFT_error ("problem with jobparam", calledby="wann_read_inp")
endif
else
CALL juDFT_error ("parameter needed in wan90version", calledby="wann_read_inp")
endif
if(abs(version_real-1.1).lt.1.e-9)THEN
this%wan90version=1
ELSEIF(abs(version_real-1.2).lt.1.e-9)THEN
this%wan90version=2
ELSEIF(abs(version_real-2.0).lt.1.e-9)THEN
this%wan90version=3
ELSE
CALL judft_error ("chosen w90 version unknown", calledby="wann_read_inp")
endif
!Not done
! ELSEIF(this%jobList(i).EQ.'ikptstart')THEN
! this%l_ikptstart=.TRUE.
......
greensf/CMakeLists.txt
View file @ dea084ba
......@@ -7,6 +7,7 @@ greensf/greensfCalcRealPart.F90
greensf/greensfUtils.f90
greensf/greensfPostProcess.F90
greensf/kkintgr.f90
greensf/lorentzian_smooth.f90
greensf/kk_cutoff.F90
greensf/hybridization.f90
greensf/occmtx.f90
......
greensf/greensfCalcRealPart.F90
View file @ dea084ba
......@@ -23,7 +23,7 @@ MODULE m_greensfCalcRealPart
IMPLICIT NONE
INTEGER, PARAMETER :: int_method(3) = (/method_direct,method_direct,method_maclaurin/)
INTEGER, PARAMETER :: int_method(3) = [method_direct,method_direct,method_maclaurin]
CONTAINS
......@@ -44,11 +44,12 @@ MODULE m_greensfCalcRealPart
INTEGER :: spin_cut,nn,natom,contourShape,dummy
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_skip
REAL :: fac,del,eb,et,fixedCutoff
REAL, ALLOCATABLE :: eMesh(:)
!Get the information on the real axis energy mesh
CALL gfinp%eMesh(ef,del,eb,et)
CALL gfinp%eMesh(ef,del_out=del,eb_out=eb,et_out=et,eMesh=eMesh)
nspins = MERGE(3,input%jspins,gfinp%l_mperp)
......@@ -99,13 +100,13 @@ MODULE m_greensfCalcRealPart
IF(i_gf.LT.1 .OR. i_gf.GT.gfinp%n) CYCLE !Make sure to not produce segfaults with mpi
!Get the information of ith current element
l = gfinp%elem(i_gf)%l
lp = gfinp%elem(i_gf)%lp
nType = gfinp%elem(i_gf)%atomType
l = gfinp%elem(i_gf)%l
lp = gfinp%elem(i_gf)%lp
nType = gfinp%elem(i_gf)%atomType
nTypep = gfinp%elem(i_gf)%atomTypep
contourShape = gfinp%contour(gfinp%elem(i_gf)%iContour)%shape
contourShape = gfinp%contour(gfinp%elem(i_gf)%iContour)%shape
l_fixedCutoffset = gfinp%elem(i_gf)%l_fixedCutoffset
fixedCutoff = gfinp%elem(i_gf)%fixedCutoff
fixedCutoff = gfinp%elem(i_gf)%fixedCutoff
CALL uniqueElements_gfinp(gfinp,dummy,ind=i_gf,indUnique=i_elem)
......@@ -124,11 +125,9 @@ MODULE m_greensfCalcRealPart
greensfImagPart%kkintgr_cutoff(i_gf,:,2) = INT((fixedCutoff+ef-eb)/del)+1
ELSE
!For all other elements we just use ef+elup as a hard cutoff
!(maybe give option to specify outside of changing the realAxis grid)
greensfImagPart%kkintgr_cutoff(i_gf,:,1) = 1
greensfImagPart%kkintgr_cutoff(i_gf,:,2) = gfinp%ne
ENDIF
!
!Perform the Kramers-Kronig-Integration if not already calculated
!
......@@ -136,23 +135,43 @@ MODULE m_greensfCalcRealPart
DO jspin = spin_start, spin_end
spin_cut = MERGE(1,jspin,jspin.GT.2)
kkcut = greensfImagPart%kkintgr_cutoff(i_gf,spin_cut,2)
!------------------------------------------------------------
! Set everything above the cutoff in the imaginary part to 0
! We do this explicitely because when we just use the hard cutoff index
! Things might get lost when the imaginary part is smoothed explicitely
!------------------------------------------------------------
IF(kkcut.ne.SIZE(eMesh)) THEN
greensfImagPart%sphavg(kkcut+1:,-l:l,-l:l,i_elem,jspin) = 0.0
ENDIF
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
!Don't waste time on empty elements
l_skip = .FALSE.
DO ie = 1, SIZE(eMesh)
IF(ABS(greensfImagPart%sphavg(ie,m,mp,i_elem,jspin)).GT.1e-12) EXIT
IF(ie==SIZE(eMesh)) l_skip = .TRUE.
ENDDO
IF(l_skip) THEN
g(i_gf)%gmmpMat(:,m,mp,jspin,ipm) = cmplx_0
CYCLE
ENDIF
IF(gfinp%l_sphavg) THEN
CALL kkintgr(greensfImagPart%sphavg(:,m,mp,i_elem,jspin),eb,del,kkcut,&
g(i_gf)%gmmpMat(:,m,mp,jspin,ipm),g(i_gf)%contour%e,(ipm.EQ.2),g(i_gf)%contour%nz,int_method(contourShape))
CALL kkintgr(greensfImagPart%sphavg(:,m,mp,i_elem,jspin),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
CALL kkintgr(greensfImagPart%uu(:,m,mp,i_elem,jspin),eb,del,kkcut,&
g(i_gf)%uu(:,m,mp,jspin,ipm),g(i_gf)%contour%e,(ipm.EQ.2),g(i_gf)%contour%nz,int_method(contourShape))
CALL kkintgr(greensfImagPart%dd(:,m,mp,i_elem,jspin),eb,del,kkcut,&
g(i_gf)%dd(:,m,mp,jspin,ipm),g(i_gf)%contour%e,(ipm.EQ.2),g(i_gf)%contour%nz,int_method(contourShape))
CALL kkintgr(greensfImagPart%du(:,m,mp,i_elem,jspin),eb,del,kkcut,&
g(i_gf)%du(:,m,mp,jspin,ipm),g(i_gf)%contour%e,(ipm.EQ.2),g(i_gf)%contour%nz,int_method(contourShape))
CALL kkintgr(greensfImagPart%ud(:,m,mp,i_elem,jspin),eb,del,kkcut,&
g(i_gf)%ud(:,m,mp,jspin,ipm),g(i_gf)%contour%e,(ipm.EQ.2),g(i_gf)%contour%nz,int_method(contourShape))
CALL kkintgr(greensfImagPart%uu(:,m,mp,i_elem,jspin),eMesh,g(i_gf)%contour%e,(ipm.EQ.2),&
g(i_gf)%uu(:,m,mp,jspin,ipm),int_method(contourShape))
CALL kkintgr(greensfImagPart%ud(:,m,mp,i_elem,jspin),eMesh,g(i_gf)%contour%e,(ipm.EQ.2),&
g(i_gf)%ud(:,m,mp,jspin,ipm),int_method(contourShape))
CALL kkintgr(greensfImagPart%du(:,m,mp,i_elem,jspin),eMesh,g(i_gf)%contour%e,(ipm.EQ.2),&
g(i_gf)%du(:,m,mp,jspin,ipm),int_method(contourShape))
CALL kkintgr(greensfImagPart%dd(:,m,mp,i_elem,jspin),eMesh,g(i_gf)%contour%e,(ipm.EQ.2),&
g(i_gf)%dd(:,m,mp,jspin,ipm),int_method(contourShape))
ENDIF
ENDDO
ENDDO
......
greensf/kk_cutoff.F90
View file @ dea084ba
......@@ -26,36 +26,19 @@ MODULE m_kk_cutoff
REAL, INTENT(IN) :: e_top
INTEGER, INTENT(INOUT) :: cutoff(:,:)
CHARACTER(len=5) :: filename
INTEGER :: i,m,n_c,ispin,spins_cut
INTEGER :: m,ispin,spins_cut
REAL :: lowerBound,upperBound,integral,n_states,scale,e_cut
REAL :: projDOS(ne,jspins)
projDOS = 0.0
REAL, ALLOCATABLE :: projDOS(:,:)
!Calculate the trace over m,mp of the Imaginary Part matrix to obtain the projected DOS
!n_f(e) = -1/pi * TR[Im(G_f(e))]
ALLOCATE(projDOS(ne,jspins),source=0.0)
DO ispin = 1, jspins
DO m = -l , l
DO i = 1, ne
projDOS(i,ispin) = projDOS(i,ispin) + im(i,m,m,ispin)
ENDDO
projDOS(:,ispin) = projDOS(:,ispin) - 1/pi_const * im(:,m,m,ispin)
ENDDO
ENDDO
projDOS = -1/pi_const*projDOS
!#ifdef CPP_DEBUG
!DO ispin = 1, jspins
! WRITE(filename,9010) ispin
!9010 FORMAT("projDOS",I1)
! OPEN(unit=1337,file=filename,status="replace")
! DO i = 1, ne
! WRITE(1337,"(2f14.8)") (i-1)*del+e_bot,projDOS(i,ispin)
! ENDDO
! CLOSE(unit=1337)
!ENDDO
!#endif
spins_cut = MERGE(1,jspins,noco%l_noco.AND.l_mperp)
n_states = (2*l+1) * MERGE(2.0,2.0/jspins,noco%l_noco.AND.l_mperp)
......@@ -69,6 +52,8 @@ MODULE m_kk_cutoff
!Check the integral up to the hard cutoff
!----------------------------------------
IF(spins_cut.EQ.1 .AND.jspins.EQ.2) projDOS(:,1) = projDOS(:,1) + projDOS(:,2)
!Initial complete integral
integral = trapz(projDOS(:,ispin),del,ne)
#ifdef CPP_DEBUG
......@@ -91,19 +76,23 @@ MODULE m_kk_cutoff
! If the integral is to small we terminate here to avoid problems
CALL juDFT_warn("Integral over DOS too small for f -> increase elup(<1htr) or numbands", calledby="kk_cutoff")
ENDIF
ELSE IF((integral.GT.n_states).AND.((integral-n_states).GT.0.00001)) THEN
ELSE
!IF the integral is bigger than 2l+1, search for the cutoff using the bisection method
lowerBound = e_bot
upperBound = e_top
DO WHILE(upperBound-lowerBound.GT.del)
DO WHILE(ABS(upperBound-lowerBound).GT.del/2.0)
e_cut = (lowerBound+upperBound)/2.0
cutoff(ispin,2) = INT((e_cut-e_bot)/del)+1
!Integrate the DOS up to the cutoff
integral = trapz(projDOS(:,ispin),del,cutoff(ispin,2))
IF(integral.LT.n_states) THEN
IF(ABS(integral-n_states).LT.1e-12) THEN
EXIT
ELSE IF(integral.LT.n_states) THEN
!integral to small -> choose the right interval
lowerBound = e_cut
ELSE IF(integral.GT.n_states) THEN
......@@ -118,7 +107,8 @@ MODULE m_kk_cutoff
WRITE(*,*) "CORRESPONDING ENERGY", e_cut
WRITE(*,*) "INTEGRAL OVER projDOS with cutoff: ", integral
#endif
IF(spins_cut.EQ.1.AND.jspins.EQ.2) cutoff(2,2) = cutoff(1,2)
!Copy cutoff to second spin if only one was calculated
IF(spins_cut.EQ.1 .AND. jspins.EQ.2) cutoff(2,2) = cutoff(1,2)
ENDIF
ENDDO
......
greensf/kkintgr.f90
View file @ dea084ba
......@@ -14,7 +14,6 @@ MODULE m_kkintgr
! TODO: Look at FFT for Transformation
! How to do changing imaginary parts
!------------------------------------------------------------------------------
USE ieee_arithmetic
USE m_constants
USE m_juDFT
......@@ -25,15 +24,11 @@ MODULE m_kkintgr
INTEGER, PARAMETER :: method_direct = 3
INTEGER, PARAMETER :: method_fft = 4
CHARACTER(len=10), PARAMETER :: smooth_method = 'lorentzian' !(or gaussian)
!PARAMETER FOR LORENTZIAN SMOOTHING
REAL, PARAMETER :: cut = 1e-8
CHARACTER(len=10), PARAMETER :: smooth_method = 'lorentzian' !(lorentzian or gaussian)
CONTAINS
SUBROUTINE kkintgr(im,eb,del,ne,g,ez,l_conjg,nz,method)
SUBROUTINE kkintgr(im,eMesh,ez,l_conjg,g,method)
!calculates the Kramer Kronig Transformation on the same contour where the imaginary part was calculated
!Re(G(E+i * delta)) = -1/pi * int_bot^top dE' P(1/(E-E')) * Im(G(E'+i*delta))
......@@ -41,33 +36,30 @@ MODULE m_kkintgr
!The dominant source of error for this routine is a insufficiently dense energy mesh on the real axis
!TODO: Some way to estimate the error (maybe search for the sharpest peak and estimate from width)
USE m_smooth
USE m_lorentzian_smooth
!Information about the integrand
REAL, INTENT(IN) :: im(:) !Imaginary part of the green's function on the real axis
REAL, INTENT(IN) :: eb !Bottom energy cutoff
REAL, INTENT(IN) :: del !Energy step on the real axis
INTEGER, INTENT(IN) :: ne !Number of energy points on the real axis
!Information about the complex energy contour
COMPLEX, INTENT(INOUT) :: g(:) !Green's function on the complex plane
REAL, INTENT(IN) :: eMesh(:) !Energy grid on the real axis
COMPLEX, INTENT(IN) :: ez(:) !Complex energy contour
LOGICAL, INTENT(IN) :: l_conjg !Switch determines wether we calculate g on the complex conjugate of the contour ez
INTEGER, INTENT(IN) :: nz !Number of energy points on the complex contour
!Information about the method
COMPLEX, INTENT(INOUT) :: g(:) !Green's function on the complex plane
INTEGER, INTENT(IN) :: method !Integer associated with the method to be used (definitions above)
INTEGER :: iz,izp,n1,n2,i
INTEGER :: iz,izp,n1,n2,i,ne,nz
INTEGER :: ismooth,nsmooth
REAL :: e(ne)
REAL :: eb,del
REAL :: re_n1,re_n2,im_n1,im_n2
INTEGER :: smoothInd(nz)
REAL :: sigma(nz)
REAL, ALLOCATABLE :: smoothed(:,:)
INTEGER, ALLOCATABLE :: smoothInd(:)
REAL, ALLOCATABLE :: sigma(:)
REAL, ALLOCATABLE :: smoothed(:,:)
DO i = 1, ne
e(i) = (i-1) * del + eb
ENDDO
nz = SIZE(ez)
ne = SIZE(eMesh)
eb = eMesh(1)
del = eMesh(2) - eMesh(1)
ALLOCATE(smoothInd(nz),source=0)
ALLOCATE(sigma(nz),source=0.0)
IF(method.NE.method_direct) THEN
CALL timestart("kkintgr: smoothing")
......@@ -85,19 +77,19 @@ MODULE m_kkintgr
smoothInd(iz) = nsmooth
sigma(nsmooth) = AIMAG(ez(iz))
ENDDO outer
ALLOCATE(smoothed(nsmooth,ne), source=0.0)
ALLOCATE(smoothed(ne,nsmooth), source=0.0)
!$OMP PARALLEL DEFAULT(none) &
!$OMP SHARED(nsmooth,smoothed,sigma,ne,e,im) &
!$OMP SHARED(nsmooth,smoothed,sigma,ne,eMesh,im) &
!$OMP PRIVATE(ismooth)
!$OMP DO
DO ismooth = 1, nsmooth
smoothed(ismooth,:) = im(:ne)
smoothed(:,ismooth) = im(:ne)
IF(ABS(sigma(ismooth)).LT.1e-12) CYCLE
SELECT CASE (TRIM(ADJUSTL(smooth_method)))
CASE('lorentzian')
CALL lorentzian_smooth(e,smoothed(ismooth,:),sigma(ismooth),ne)
CALL lorentzian_smooth(eMesh,smoothed(:,ismooth),sigma(ismooth),ne)
CASE('gaussian')
CALL smooth(e,smoothed(ismooth,:),sigma(ismooth),ne)
CALL smooth(eMesh,smoothed(:,ismooth),sigma(ismooth),ne)
CASE DEFAULT
CALL juDFT_error("No valid smooth_method set",&
hint="This is a bug in FLEUR, please report",&
......@@ -113,14 +105,14 @@ MODULE m_kkintgr
CALL timestart("kkintgr: integration")
!$OMP PARALLEL DEFAULT(none) &
!$OMP SHARED(nz,ne,method,del,eb,l_conjg) &
!$OMP SHARED(g,ez,im,e,smoothed,smoothInd) &
!$OMP SHARED(g,ez,im,smoothed,smoothInd) &
!$OMP PRIVATE(iz,n1,n2,re_n1,re_n2,im_n1,im_n2)
!$OMP DO
DO iz = 1, nz
SELECT CASE(method)
CASE(method_direct)
g(iz) = g_circle(im,ne,MERGE(conjg(ez(iz)),ez(iz),l_conjg),del,eb)
g(iz) = kk_direct(im,ne,MERGE(conjg(ez(iz)),ez(iz),l_conjg),del,eb)
CASE(method_maclaurin, method_deriv)
!Use the previously smoothed version and interpolate after
!Next point to the left
......@@ -128,31 +120,28 @@ MODULE m_kkintgr
!next point to the right
n2 = n1 + 1
!Here we perform the Kramers-kronig-Integration
re_n2 = re_ire(smoothed(smoothInd(iz),:),ne,n2,method)
re_n1 = re_ire(smoothed(smoothInd(iz),:),ne,n1,method)
re_n2 = kk_num(smoothed(:,smoothInd(iz)),ne,n2,method)
re_n1 = kk_num(smoothed(:,smoothInd(iz)),ne,n1,method)