Commit 248cdd03 authored by Henning Janssen's avatar Henning Janssen

Started integrating new input types gfinp and hub1inp

   Now remove the redundant parameters
parent 29b51239
......@@ -12,7 +12,7 @@ CONTAINS
SUBROUTINE cdnval(eig_id, mpi,kpts,jspin,noco,input,banddos,cell,atoms,enpara,stars,&
vacuum,sphhar,sym,vTot,oneD,cdnvalJob,den,regCharges,dos,results,&
moments,hub1,coreSpecInput,mcd,slab,orbcomp,greensfCoeffs,angle)
moments,gfinp,hub1inp,hub1data,coreSpecInput,mcd,slab,orbcomp,greensfCoeffs,angle)
!************************************************************************************
! This is the FLEUR valence density generator
......@@ -73,13 +73,15 @@ SUBROUTINE cdnval(eig_id, mpi,kpts,jspin,noco,input,banddos,cell,atoms,enpara,st
TYPE(t_kpts), INTENT(IN) :: kpts
TYPE(t_sphhar), INTENT(IN) :: sphhar
TYPE(t_atoms), INTENT(IN) :: atoms
TYPE(t_gfinp), INTENT(IN) :: gfinp
TYPE(t_hub1inp), INTENT(IN) :: hub1inp
TYPE(t_potden), INTENT(IN) :: vTot
TYPE(t_cdnvalJob), INTENT(IN) :: cdnvalJob
TYPE(t_potden), INTENT(INOUT) :: den
TYPE(t_regionCharges), INTENT(INOUT) :: regCharges
TYPE(t_dos), INTENT(INOUT) :: dos
TYPE(t_moments), INTENT(INOUT) :: moments
TYPE(t_hub1ham), OPTIONAL, INTENT(INOUT) :: hub1
TYPE(t_hub1data), OPTIONAL, INTENT(INOUT) :: hub1data
TYPE(t_coreSpecInput), OPTIONAL, INTENT(IN) :: coreSpecInput
TYPE(t_mcd), OPTIONAL, INTENT(INOUT) :: mcd
TYPE(t_slab), OPTIONAL, INTENT(INOUT) :: slab
......@@ -95,7 +97,7 @@ SUBROUTINE cdnval(eig_id, mpi,kpts,jspin,noco,input,banddos,cell,atoms,enpara,st
#endif
! Local Scalars
INTEGER :: ikpt,ikpt_i,jsp_start,jsp_end,ispin,jsp,ne
INTEGER :: ikpt,ikpt_i,jsp_start,jsp_end,ispin,jsp
INTEGER :: iErr,nbands,noccbd,iType
INTEGER :: skip_t,skip_tt,nbasfcn
LOGICAL :: l_orbcomprot, l_real, l_dosNdir, l_corespec
......@@ -123,7 +125,7 @@ SUBROUTINE cdnval(eig_id, mpi,kpts,jspin,noco,input,banddos,cell,atoms,enpara,st
l_real = sym%invs.AND.(.NOT.noco%l_soc).AND.(.NOT.noco%l_noco)
l_dosNdir = banddos%dos.AND.(banddos%ndir.EQ.-3)
IF (noco%l_mperp.OR.input%l_gfmperp) THEN
IF (noco%l_mperp.OR.gfinp%l_mperp) THEN
! when the off-diag. part of the desinsity matrix, i.e. m_x and
! m_y, is calculated inside the muffin-tins (l_mperp = T), cdnval
! is called only once. therefore, several spin loops have been
......@@ -170,9 +172,9 @@ SUBROUTINE cdnval(eig_id, mpi,kpts,jspin,noco,input,banddos,cell,atoms,enpara,st
DO iType = 1, atoms%ntype
DO ispin = jsp_start, jsp_end
CALL genMTBasis(atoms,enpara,vTot,mpi,iType,ispin,usdus,f(:,:,0:,ispin),g(:,:,0:,ispin),flo(:,:,:,ispin),input%l_dftspinpol)
CALL genMTBasis(atoms,enpara,vTot,mpi,iType,ispin,usdus,f(:,:,0:,ispin),g(:,:,0:,ispin),flo(:,:,:,ispin),hub1inp%l_dftspinpol)
END DO
IF (noco%l_mperp.OR.input%l_gfmperp) CALL denCoeffsOffdiag%addRadFunScalarProducts(atoms,f,g,flo,iType)
IF (noco%l_mperp.OR.gfinp%l_mperp) CALL denCoeffsOffdiag%addRadFunScalarProducts(atoms,f,g,flo,iType)
IF (banddos%l_mcd) CALL mcd_init(atoms,input,vTot%mt(:,0,:,:),g,f,mcd,iType,jspin)
IF (l_coreSpec) CALL corespec_rme(atoms,input,iType,29,input%jspins,jspin,results%ef,&
atoms%msh,vTot%mt(:,0,:,:),f,g)
......@@ -185,14 +187,14 @@ SUBROUTINE cdnval(eig_id, mpi,kpts,jspin,noco,input,banddos,cell,atoms,enpara,st
jsp = MERGE(1,jspin,noco%l_noco)
!Get the energy mesh for the tetrahedron method
IF(atoms%n_gf>0) CALL greensfCoeffs%eMesh(ne,eMesh)
IF(gfinp%n>0) CALL gfinp%eMesh(results%ef,eMesh=eMesh)
DO ikpt_i = 1,size(cdnvalJob%k_list)
ikpt=cdnvalJob%k_list(ikpt_i)
CALL lapw%init(input,noco, kpts,atoms,sym,ikpt,cell,.false., mpi)
skip_t = skip_tt
ev_list=cdnvaljob%compact_ev_list(ikpt_i,banddos%dos.OR.atoms%n_gf>0)
ev_list=cdnvaljob%compact_ev_list(ikpt_i,banddos%dos.OR.gfinp%n>0)
noccbd = SIZE(ev_list)
we = cdnvalJob%weights(ev_list,ikpt)
eig = results%eig(ev_list,ikpt,jsp)
......@@ -231,13 +233,13 @@ SUBROUTINE cdnval(eig_id, mpi,kpts,jspin,noco,input,banddos,cell,atoms,enpara,st
! valence density in the atomic spheres
CALL eigVecCoeffs%init(input,atoms,noco,jspin,noccbd)
IF (atoms%n_gf.GT.0.AND.(input%tria.OR.input%gfTet)) THEN
IF (gfinp%n.GT.0.AND.(input%tria.OR.input%gfTet)) THEN
CALL timestart("TetrahedronWeights")
ALLOCATE(dosWeights(greensfCoeffs%ne,noccbd))
ALLOCATE(resWeights(greensfCoeffs%ne,noccbd))
resWeights = 0.0
ALLOCATE(dosWeights(gfinp%ne,noccbd),source=0.0)
ALLOCATE(resWeights(gfinp%ne,noccbd),source=0.0)
ALLOCATE(dosBound(noccbd,2))
CALL tetrahedronInit(kpts,ikpt,results%eig(ev_list,:,jsp),noccbd,eMesh,ne,input%film,dosWeights,bounds=dosBound,dos=.TRUE.)
CALL tetrahedronInit(kpts,ikpt,results%eig(ev_list,:,jsp),noccbd,eMesh,gfinp%ne,&
input%film,dosWeights,bounds=dosBound,dos=.TRUE.)
CALL timestop("TetrahedronWeights")
ENDIF
DO ispin = jsp_start, jsp_end
......@@ -249,8 +251,8 @@ SUBROUTINE cdnval(eig_id, mpi,kpts,jspin,noco,input,banddos,cell,atoms,enpara,st
IF (atoms%n_u.GT.0.AND.noco%l_mperp.AND.(ispin==jsp_end)) CALL n_mat21(atoms,sym,angle,noccbd,we,denCoeffsOffdiag,&
eigVecCoeffs,den%mmpMat(:,:,:,3))
IF (atoms%n_gf.GT.0) CALL bzIntegrationGF(atoms,sym,input,angle,ispin,noccbd,dosWeights,resWeights,dosBound,kpts%wtkpt(ikpt),&
eig,denCoeffsOffdiag,usdus,eigVecCoeffs,greensfCoeffs,ispin==jsp_end)
IF (gfinp%n.GT.0) CALL bzIntegrationGF(atoms,gfinp,sym,input,angle,ispin,noccbd,dosWeights,resWeights,dosBound,kpts%wtkpt(ikpt),&
results%ef,eig,denCoeffsOffdiag,usdus,eigVecCoeffs,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
......@@ -277,7 +279,7 @@ SUBROUTINE cdnval(eig_id, mpi,kpts,jspin,noco,input,banddos,cell,atoms,enpara,st
IF(l_coreSpec) CALL corespec_dos(atoms,usdus,ispin,atoms%lmaxd*(atoms%lmaxd+2),kpts%nkpt,ikpt,input%neig,&
noccbd,results%ef,banddos%sig_dos,eig,we,eigVecCoeffs)
END DO ! end loop over ispin
IF (atoms%n_gf.GT.0.AND.(input%tria.OR.input%gfTet)) DEALLOCATE(dosWeights,resWeights,dosBound)
IF (gfinp%n.GT.0.AND.(input%tria.OR.input%gfTet)) DEALLOCATE(dosWeights,resWeights,dosBound)
IF (noco%l_mperp) CALL denCoeffsOffdiag%calcCoefficients(atoms,sphhar,sym,eigVecCoeffs,we,noccbd)
IF ((banddos%dos.OR.banddos%vacdos.OR.input%cdinf).AND.(banddos%ndir.GT.0)) THEN
......@@ -288,13 +290,13 @@ 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,&
CALL mpi_col_den(mpi,sphhar,atoms,oneD,stars,vacuum,input,noco,gfinp,ispin,regCharges,dos,&
results,denCoeffs,orb,denCoeffsOffdiag,den,mcd,slab,orbcomp,greensfCoeffs)
END DO
#endif
CALL cdnmt(mpi,input%jspins,input,atoms,sym,sphhar,noco,jsp_start,jsp_end,&
enpara,vTot%mt(:,0,:,:),denCoeffs,usdus,orb,denCoeffsOffdiag,moments,den%mt,hub1,input%l_dftspinpol)
CALL cdnmt(mpi,input%jspins,input,atoms,sym,sphhar,noco,jsp_start,jsp_end,enpara,&
vTot%mt(:,0,:,:),denCoeffs,usdus,orb,denCoeffsOffdiag,moments,den%mt,hub1inp,hub1data)
IF (mpi%irank==0) THEN
IF (l_coreSpec) CALL corespec_ddscs(jspin,input%jspins)
......
......@@ -11,7 +11,7 @@ MODULE m_cdnmt
!***********************************************************************
CONTAINS
SUBROUTINE cdnmt(mpi,jspd,input,atoms,sym,sphhar,noco,jsp_start,jsp_end,enpara,&
vr,denCoeffs,usdus,orb,denCoeffsOffdiag,moments,rho,hub1,l_dftspinpol)
vr,denCoeffs,usdus,orb,denCoeffsOffdiag,moments,rho,hub1inp,hub1data)
use m_constants,only: sfp_const
USE m_rhosphnlo
USE m_radfun
......@@ -25,11 +25,11 @@ CONTAINS
TYPE(t_noco), INTENT(IN) :: noco
TYPE(t_sphhar), INTENT(IN) :: sphhar
TYPE(t_atoms), INTENT(IN) :: atoms
TYPE(t_sym), INTENT(IN) :: sym
TYPE(t_sym), INTENT(IN) :: sym
TYPE(t_enpara), INTENT(IN) :: enpara
TYPE(t_hub1inp), INTENT(IN) :: hub1inp
TYPE(t_moments), INTENT(INOUT) :: moments
TYPE(t_hub1ham), OPTIONAL, INTENT(INOUT) :: hub1
LOGICAL, INTENT(IN) :: l_dftspinpol
TYPE(t_hub1data), OPTIONAL, INTENT(INOUT) :: hub1data
! .. Scalar Arguments ..
INTEGER, INTENT (IN) :: jsp_start,jsp_end,jspd
......@@ -69,7 +69,7 @@ CONTAINS
ENDIF
! !$OMP PARALLEL DEFAULT(none) &
! !$OMP SHARED(usdus,rho,moments,qmtl,hub1) &
! !$OMP SHARED(usdus,rho,moments,qmtl,hub1inp,hub1data) &
! !$OMP SHARED(atoms,jsp_start,jsp_end,enpara,vr,denCoeffs,sphhar,l_dftspinpol)&
! !$OMP SHARED(orb,noco,denCoeffsOffdiag,jspd)&
! !$OMP PRIVATE(itype,na,ispin,l,rho21,f,g,nodeu,noded,wronk,i,j,s,qmtllo,qmtt,nd,lh,lp,llp,llpb,cs)
......@@ -102,7 +102,7 @@ CONTAINS
!In the case of a spin-polarized calculation with Hubbard 1 we want to treat
!the correlated orbitals with a non-spin-polarized basis
IF(l_hia.AND.jspd.EQ.2.AND..NOT.l_dftspinpol) THEN
IF(l_hia.AND.jspd.EQ.2.AND..NOT.hub1inp%l_dftspinpol) THEN
vrTmp = (vr(:,itype,1) + vr(:,itype,2))/2.0
ELSE
vrTmp = vr(:,itype,ispin)
......@@ -149,11 +149,12 @@ CONTAINS
moments%chmom(itype,ispin) = qmtt
!Get the magnetic moment for the shells where we defined additional exchange splittings for DFT+Hubbard 1
IF(PRESENT(hub1)) THEN
IF(PRESENT(hub1data)) THEN
DO i_hia = 1, atoms%n_hia
IF(atoms%lda_u(atoms%n_u+i_hia)%atomType.NE.itype) CYCLE
DO i_exc = 1, hub1%n_exc_given(i_hia)
hub1%mag_mom(i_hia,i_exc) = hub1%mag_mom(i_hia,i_exc) + (-1)**(ispin-1) * qmtl(hub1%exc_l(i_hia,i_exc),ispin,itype)
DO i_exc = 1, hub1inp%n_exc(i_hia)
hub1data%mag_mom(i_hia,i_exc) = hub1data%mag_mom(i_hia,i_exc) + &
(-1)**(ispin-1) * qmtl(hub1inp%exc_l(i_hia,i_exc),ispin,itype)
ENDDO
ENDDO
ENDIF
......
......@@ -20,7 +20,7 @@ CONTAINS
!>@author D. Wortmann
SUBROUTINE eigen(mpi,stars,sphhar,atoms,xcpot,sym,kpts,vacuum,input,&
cell,enpara,banddos,noco,oneD,mpdata,hybinp,hybdat,&
iter,eig_id,results,inden,v,vx,hub1)
iter,eig_id,results,inden,v,vx,hub1inp,hub1data)
#include"cpp_double.h"
USE m_types
......@@ -64,7 +64,8 @@ CONTAINS
TYPE(t_sphhar),INTENT(IN) :: sphhar
TYPE(t_atoms),INTENT(IN) :: atoms
TYPE(t_potden),INTENT(IN) :: inden !
TYPE(t_hub1ham),INTENT(INOUT):: hub1
TYPE(t_hub1inp),INTENT(IN) :: hub1inp
TYPE(t_hub1data),INTENT(INOUT):: hub1data
TYPE(t_potden), INTENT(IN) :: vx
TYPE(t_potden),INTENT(INOUT) :: v !u_setup will modify the potential matrix
......@@ -129,7 +130,7 @@ CONTAINS
! Set up and solve the eigenvalue problem
! loop over spins
! set up k-point independent t(l'm',lm) matrices
CALL mt_setup(atoms,sym,sphhar,input,noco,enpara,hub1,inden,v,mpi,results,td,ud)
CALL mt_setup(atoms,sym,sphhar,input,noco,enpara,hub1inp,hub1data,inden,v,mpi,results,td,ud)
neigBuffer = 0
results%neig = 0
......
......@@ -7,7 +7,7 @@
MODULE m_mt_setup
CONTAINS
SUBROUTINE mt_setup(atoms,sym,sphhar,input,noco,enpara,hub1,inden,vTot,mpi,results,td,ud)
SUBROUTINE mt_setup(atoms,sym,sphhar,input,noco,enpara,hub1inp,hub1data,inden,vTot,mpi,results,td,ud)
USE m_types
USE m_usetup
USE m_tlmplm_cholesky
......@@ -27,7 +27,8 @@ CONTAINS
TYPE(t_potden),INTENT(INOUT) :: vTot
TYPE(t_tlmplm),INTENT(INOUT) :: td
TYPE(t_usdus),INTENT(INOUT) :: ud
TYPE(t_hub1ham),INTENT(INOUT) :: hub1
TYPE(t_hub1inp),INTENT(IN) :: hub1inp
TYPE(t_hub1data),INTENT(INOUT)::hub1data
INTEGER:: jsp
......@@ -50,7 +51,7 @@ CONTAINS
!Setup of soc parameters for first-variation SOC
IF (noco%l_soc.AND.noco%l_noco.AND..NOT.noco%l_ss) THEN
CALL spnorb(atoms,noco,input,mpi,enpara,vTot%mt,ud,td%rsoc,.FALSE.,hub1)
CALL spnorb(atoms,noco,input,mpi,enpara,vTot%mt,ud,td%rsoc,.FALSE.,hub1inp,hub1data)
END IF
END SUBROUTINE mt_setup
......
......@@ -21,7 +21,7 @@ MODULE m_eigenso
!
CONTAINS
SUBROUTINE eigenso(eig_id,mpi,stars,vacuum,atoms,sphhar,&
sym,cell,noco,input,kpts,oneD,vTot,enpara,results,hub1)
sym,cell,noco,input,kpts,oneD,vTot,enpara,results,hub1inp,hub1data)
USE m_types
USE m_eig66_io, ONLY : read_eig,write_eig
......@@ -48,7 +48,8 @@ CONTAINS
TYPE(t_potden),INTENT(IN) :: vTot
TYPE(t_enpara),INTENT(IN) :: enpara
TYPE(t_results),INTENT(INOUT) :: results
TYPE(t_hub1ham),OPTIONAL,INTENT(INOUT) :: hub1
TYPE(t_hub1inp),OPTIONAL,INTENT(IN) :: hub1inp
TYPE(t_hub1data),OPTIONAL,INTENT(INOUT) :: hub1data
#ifdef CPP_MPI
INCLUDE 'mpif.h'
......@@ -114,7 +115,7 @@ CONTAINS
! ..
!Get spin-orbit coupling matrix elements
CALL spnorb( atoms,noco,input,mpi, enpara,vTot%mt,usdus,rsoc,.TRUE.,hub1)
CALL spnorb( atoms,noco,input,mpi, enpara,vTot%mt,usdus,rsoc,.TRUE.,hub1inp,hub1data)
!
......
......@@ -12,7 +12,7 @@ MODULE m_spnorb
! using the functions anglso and sgml.
!*********************************************************************
CONTAINS
SUBROUTINE spnorb(atoms,noco,input,mpi, enpara, vr, usdus, rsoc,l_angles,hub1)
SUBROUTINE spnorb(atoms,noco,input,mpi, enpara, vr, usdus, rsoc,l_angles,hub1inp,hub1data)
USE m_sorad
USE m_constants, only : hartree_to_ev_const
USE m_types
......@@ -26,7 +26,8 @@ CONTAINS
TYPE(t_usdus),INTENT(INOUT) :: usdus
TYPE(t_rsoc),INTENT(OUT) :: rsoc
LOGICAL,INTENT(IN) :: l_angles
TYPE(t_hub1ham),OPTIONAL,INTENT(INOUT) :: hub1
TYPE(t_hub1inp),OPTIONAL,INTENT(IN) :: hub1inp
TYPE(t_hub1data),OPTIONAL,INTENT(INOUT) :: hub1data
! ..
! ..
! .. Array Arguments ..
......@@ -74,12 +75,12 @@ CONTAINS
ENDDO
!Read in SOC-parameter for shell with hubbard 1
IF(PRESENT(hub1).AND.mpi%irank.EQ.0) THEN
IF(PRESENT(hub1inp).AND.mpi%irank.EQ.0) THEN
DO i_hia = 1, atoms%n_hia
IF(hub1%l_soc_given(i_hia)) CYCLE
IF(hub1inp%l_soc_given(i_hia)) CYCLE
n = atoms%lda_u(atoms%n_u+i_hia)%atomType
l = atoms%lda_u(atoms%n_u+i_hia)%l
hub1%xi(i_hia) = 2.0*rsoc%rsopp(n,l,1,1)*hartree_to_ev_const
IF(PRESENT(hub1data)) hub1data%xi(i_hia) = 2.0*rsoc%rsopp(n,l,1,1)*hartree_to_ev_const
ENDDO
ENDIF
......
......@@ -366,7 +366,7 @@ CONTAINS
ENDDO
END FUNCTION find_gfelem
SUBROUTINE eMesh_gfinp(this,ef,del_out,eb_out,eMesh)
SUBROUTINE eMesh_gfinp(this,ef,del_out,eb_out,et_out,eMesh)
!Gives back the information for the energy mesh on the real axis
!Energies are shifted according to the fermi level
......@@ -375,6 +375,7 @@ CONTAINS
REAL, INTENT(IN) :: ef
REAL, OPTIONAL, INTENT(INOUT) :: del_out
REAL, OPTIONAL, INTENT(INOUT) :: eb_out
REAL, OPTIONAL, INTENT(INOUT) :: et_out
REAL, ALLOCATABLE, OPTIONAL, INTENT(INOUT) :: eMesh(:)
INTEGER :: ie
......@@ -384,6 +385,7 @@ CONTAINS
del = (this%elup-this%ellow)/REAL(this%ne-1)
IF(PRESENT(eb_out)) eb_out = eb
IF(PRESENT(et_out)) et_out = ef+this%elup
IF(PRESENT(del_out)) del_out = del
IF(PRESENT(eMesh)) THEN
......
......@@ -11,7 +11,6 @@ greensf/kk_cutoff.f90
greensf/gfcalc.f90
greensf/crystalfield.f90
greensf/j0.f90
greensf/ind_greensf.f90
greensf/hybridization.f90
greensf/rot_gf.f90
greensf/occmtx.f90
......
......@@ -15,7 +15,6 @@ MODULE m_crystalfield
USE m_types
USE m_constants
USE m_kkintgr
USE m_ind_greensf
USE m_sgml
USE m_anglso
......@@ -25,37 +24,35 @@ MODULE m_crystalfield
CONTAINS
SUBROUTINE crystal_field(atoms,input,noco,greensfCoeffs,hub1,v)
SUBROUTINE crystal_field(atoms,gfinp,hub1inp,input,noco,greensfCoeffs,v,ef,hub1data)
!calculates the crystal-field matrix for the local hamiltonian
IMPLICIT NONE
!-Type Arguments
TYPE(t_greensfCoeffs), INTENT(IN) :: greensfCoeffs
TYPE(t_atoms), INTENT(IN) :: atoms
TYPE(t_gfinp), INTENT(IN) :: gfinp
TYPE(t_input), INTENT(IN) :: input
TYPE(t_noco), INTENT(IN) :: noco
TYPE(t_hub1ham), INTENT(INOUT) :: hub1
!-Array Arguments
TYPE(t_hub1inp), INTENT(IN) :: hub1inp
TYPE(t_potden), INTENT(IN) :: v !LDA+U potential (should be removed from h_loc)
REAL, INTENT(IN) :: ef
TYPE(t_hub1data), INTENT(INOUT) :: hub1data
!-Local Scalars
INTEGER i_gf,l,nType,jspin,m,mp,ie,i_hia,kkcut,i_u,isp
REAL tr,xiSOC
REAL tr,xiSOC,del,eb
COMPLEX vso
!-Local Arrays
REAL :: h_loc(-lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const,atoms%n_hia,input%jspins)
REAL :: ex(-lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const)
REAL :: integrand(greensfCoeffs%ne), norm(greensfCoeffs%ne)
REAL :: integrand(gfinp%ne), norm(gfinp%ne)
h_loc = 0.0
DO i_hia = 1, atoms%n_hia
l = atoms%lda_u(atoms%n_u+i_hia)%l
nType = atoms%lda_u(atoms%n_u+i_hia)%atomType
i_gf = ind_greensf(atoms,l,nType)
i_gf = gfinp%find(l,nType)
!---------------------------------------------------------
! Perform the integration
!---------------------------------------------------------
......@@ -65,6 +62,7 @@ MODULE m_crystalfield
! connected to the imaginary part of the greens function
! with a factor -1/pi
!---------------------------------------------------------
CALL gfinp%eMesh(ef,del,eb)
DO jspin = 1, input%jspins
!Use the same cutoffs as in the kramer kronigs integration
kkcut = greensfCoeffs%kkintgr_cutoff(i_gf,jspin,2)
......@@ -73,11 +71,11 @@ MODULE m_crystalfield
DO mp = -l, l
integrand = 0.0
DO ie = 1, kkcut
integrand(ie) = -1.0/pi_const * ((ie-1) * greensfCoeffs%del+greensfCoeffs%e_bot) &
integrand(ie) = -1.0/pi_const * ((ie-1) * del+eb) &
* REAL(greensfCoeffs%projdos(ie,m,mp,0,i_gf,jspin)/(3.0-input%jspins))
IF(m.EQ.mp) norm(ie) = norm(ie) -1.0/pi_const * REAL(greensfCoeffs%projdos(ie,m,mp,0,i_gf,jspin))/(3.0-input%jspins)
ENDDO
h_loc(m,mp,i_hia,jspin) = trapz(integrand(1:kkcut),greensfCoeffs%del,kkcut)
h_loc(m,mp,i_hia,jspin) = trapz(integrand(1:kkcut),del,kkcut)
ENDDO
ENDDO
ENDDO
......@@ -101,12 +99,12 @@ MODULE m_crystalfield
DO m = -l, l
DO mp = -l, l
isp = 3.0-2.0*jspin !1,-1
IF((ABS(noco%theta).LT.0.00001).AND.(ABS(noco%phi).LT.0.00001)) THEN
IF((ABS(noco%theta).LT.1e-5).AND.(ABS(noco%phi).LT.1e-5)) THEN
vso = CMPLX(sgml(l,m,isp,l,mp,isp),0.0)
ELSE
vso = anglso(noco%theta,noco%phi,l,m,isp,l,mp,isp)
ENDIF
h_loc(m,mp,i_hia,jspin) = h_loc(m,mp,i_hia,jspin) - REAL(vso)/2.0 * hub1%xi(i_hia)/hartree_to_ev_const
h_loc(m,mp,i_hia,jspin) = h_loc(m,mp,i_hia,jspin) - REAL(vso)/2.0 * hub1data%xi(i_hia)/hartree_to_ev_const
ENDDO
ENDDO
ENDDO
......@@ -149,32 +147,32 @@ MODULE m_crystalfield
ENDIF
!Average over spins
hub1%ccfmat(i_hia,:,:) = 0.0
hub1data%ccfmat(i_hia,:,:) = 0.0
DO m = -l, l
DO mp = -l, l
hub1%ccfmat(i_hia,m,mp) = SUM(h_loc(m,mp,i_hia,:))/2.0
hub1data%ccfmat(i_hia,m,mp) = SUM(h_loc(m,mp,i_hia,:))/2.0
!For jspins.EQ.1 we need to take care of the fact that the spin-orbit coupling is opposite in spin 1/2
IF(input%jspins.EQ.1) hub1%ccfmat(i_hia,m,mp) = (h_loc(m,mp,i_hia,1)+h_loc(-m,-mp,i_hia,1))/2.0
IF(input%jspins.EQ.1) hub1data%ccfmat(i_hia,m,mp) = (h_loc(m,mp,i_hia,1)+h_loc(-m,-mp,i_hia,1))/2.0
ENDDO
ENDDO
IF(l_debug) THEN
WRITE(*,*) "Average"
WRITE(*,"(7f7.3)") hub1%ccfmat(i_hia,-3:3,-3:3)
WRITE(*,"(7f7.3)") hub1data%ccfmat(i_hia,-3:3,-3:3)
ENDIF
DO m = -l, l
DO mp = -l, l
hub1%ccfmat(i_hia,m,mp) = (hub1%ccfmat(i_hia,m,mp)+hub1%ccfmat(i_hia,-m,-mp))/2.0
hub1%ccfmat(i_hia,-m,-mp) = hub1%ccfmat(i_hia,m,mp)
hub1data%ccfmat(i_hia,m,mp) = (hub1data%ccfmat(i_hia,m,mp)+hub1data%ccfmat(i_hia,-m,-mp))/2.0
hub1data%ccfmat(i_hia,-m,-mp) = hub1data%ccfmat(i_hia,m,mp)
ENDDO
ENDDO
IF(l_debug) THEN
WRITE(*,*) "SOC"
WRITE(*,"(7f7.3)") hub1%ccfmat(i_hia,-3:3,-3:3)
WRITE(*,"(7f7.3)") hub1data%ccfmat(i_hia,-3:3,-3:3)
ENDIF
tr = 0.0
!calculate the trace
DO m = -l, l
tr = tr + hub1%ccfmat(i_hia,m,m)
tr = tr + hub1data%ccfmat(i_hia,m,m)
ENDDO
IF(l_debug) THEN
WRITE(*,*) "TRACE"
......@@ -182,12 +180,12 @@ MODULE m_crystalfield
ENDIF
!Remove trace
DO m = -l, l
hub1%ccfmat(i_hia,m,m) = hub1%ccfmat(i_hia,m,m) - tr/(2*l+1)
hub1data%ccfmat(i_hia,m,m) = hub1data%ccfmat(i_hia,m,m) - tr/(2*l+1)
ENDDO
IF(l_debug) THEN
WRITE(*,*) "TRACELESS (eV)"
WRITE(*,"(7f7.3)") hub1%ccfmat(i_hia,-3:3,-3:3)*hartree_to_ev_const
WRITE(*,"(7f7.3)") hub1data%ccfmat(i_hia,-3:3,-3:3)*hartree_to_ev_const
ENDIF
ENDDO
END SUBROUTINE crystal_field
......
......@@ -8,14 +8,13 @@ MODULE m_gfDOS
CONTAINS
SUBROUTINE gfDOS(g,l,nType,jobID,atoms,input,ef)
SUBROUTINE gfDOS(g,l,nType,jobID,gfinp,input,ef)
!Provides the spin up/down DOS as well as the high/low J DOS
!calculated from the greens function in gfDOS.jobID
TYPE(t_greensf), INTENT(IN) :: g
TYPE(t_atoms), INTENT(IN) :: atoms
TYPE(t_gfinp), INTENT(IN) :: gfinp
TYPE(t_input), INTENT(IN) :: input
INTEGER, INTENT(IN) :: l
INTEGER, INTENT(IN) :: nType
......@@ -38,7 +37,7 @@ MODULE m_gfDOS
IF(io_error.NE.0) CALL juDFT_error("IO-error",calledby="gfDOS")
!Write out warnings
IF(.NOT.PRESENT(ef)) WRITE(3456,"(A)") "This gfDOS is not corrected to have ef=0"
IF(g%mode.NE.3) WRITE(3456,"(A)") "You are using an energy contour where the gfDOS might not be very informative"
IF(gfinp%mode.NE.3) WRITE(3456,"(A)") "You are using an energy contour where the gfDOS might not be very informative"
!Calculate the transformation between |L,S> and |J,mj> basis
ns = 2*l+1
CALL cmat%init(.TRUE.,2*ns,2*ns)
......@@ -50,7 +49,7 @@ MODULE m_gfDOS
re = 0.0
DO ipm = 1, 2 !Sum over G^+ and G^-
!Get the full gf matrix at the energy point
CALL g%get(gmat,atoms,input,iz,l,nType,ipm.EQ.2)
CALL g%get(gmat,gfinp,input,iz,l,nType,ipm.EQ.2)
!Convert to eV^-1
gmat%data_c = gmat%data_c/hartree_to_eV_const
!Calculate up/down dos
......
......@@ -3,6 +3,7 @@ MODULE m_gfcalc
!Contains the main routines called from cdnval and cdngen
USE m_juDFT
USE m_constants
USE m_types
!These dependencies are here so that all these subroutines can be called by simply adding USE m_gfcalc to the file
USE m_j0
USE m_gfDOS
......@@ -10,33 +11,31 @@ MODULE m_gfcalc
USE m_hybridization
USE m_crystalfield
IMPLICIT NONE
CONTAINS
SUBROUTINE bzIntegrationGF(atoms,sym,input,angle,ispin,nbands,dosWeights,resWeights,indBound,wtkpt,eig,denCoeffsOffdiag,&
usdus,eigVecCoeffs,greensfCoeffs,l21)
SUBROUTINE bzIntegrationGF(atoms,gfinp,sym,input,angle,ispin,nbands,dosWeights,resWeights,&
indBound,wtkpt,ef,eig,denCoeffsOffdiag,usdus,eigVecCoeffs,greensfCoeffs,l21)
USE m_greensfImag
USE m_greensfImag21
IMPLICIT NONE
!-Type Arguments
TYPE(t_atoms), INTENT(IN) :: atoms
TYPE(t_gfinp), INTENT(IN) :: gfinp
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
!-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
REAL, INTENT(IN) :: ef
LOGICAL, INTENT(IN) :: l21 !Calculate spin off-diagonal part ?
!-Array Arguments
REAL, INTENT(IN) :: resWeights(:,:)
REAL, INTENT(IN) :: dosWeights(:,:) !Precalculated tetrahedron weights for the current k-point
INTEGER, INTENT(IN) :: indBound(:,:) !Gives the range where the tetrahedron weights are non-zero
......@@ -44,9 +43,11 @@ MODULE m_gfcalc
REAL, INTENT(IN) :: angle(:) !Phases for spin-offdiagonal part
CALL timestart("Greens Function: Imaginary Part")
CALL greensfImag(atoms,sym,input,ispin,nbands,dosWeights,resWeights,indBound,wtkpt,eig,usdus,eigVecCoeffs,greensfCoeffs)
IF(input%l_gfmperp.AND.l21) THEN
CALL greensfImag21(atoms,sym,angle,input,nbands,dosWeights,resWeights,indBound,wtkpt,eig,denCoeffsOffdiag,eigVecCoeffs,greensfCoeffs)
CALL greensfImag(atoms,gfinp,sym,input,ispin,nbands,dosWeights,resWeights,indBound,&
wtkpt,ef,eig,usdus,eigVecCoeffs,greensfCoeffs)
IF(gfinp%l_mperp.AND.l21) THEN
CALL greensfImag21(atoms,gfinp,sym,angle,input,nbands,dosWeights,resWeights,indBound,&
wtkpt,ef,eig,denCoeffsOffdiag,eigVecCoeffs,greensfCoeffs)
ENDIF
CALL timestop("Greens Function: Imaginary Part")
......@@ -54,54 +55,56 @@ MODULE m_gfcalc
SUBROUTINE postProcessGF(greensf,greensfCoeffs,atoms,input,sym,noco,vTot,hub1,results,angle)
SUBROUTINE postProcessGF(greensf,greensfCoeffs,atoms,gfinp,input,sym,noco,vTot,hub1inp,hub1data,results,angle)
!contains all the modules for calculating properties from the greens function
USE m_onsite
USE m_rot_gf
TYPE(t_atoms), INTENT(IN) :: atoms
TYPE(t_gfinp), INTENT(IN) :: gfinp
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_hub1inp), INTENT(IN) :: hub1inp
TYPE(t_results), INTENT(IN) :: results
TYPE(t_potden), INTENT(IN) :: vTot
REAL, INTENT(IN) :: angle(:)
TYPE(t_hub1data), INTENT(INOUT) :: hub1data
TYPE(t_greensfCoeffs), INTENT(INOUT) :: greensfCoeffs
TYPE(t_greensf), INTENT(INOUT) :: greensf
INTEGER i_gf,l,nType
COMPLEX mmpmat(-lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const,atoms%n_gf,3)
COMPLEX mmpmat(-lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const,gfinp%n,3)
LOGICAL err
CALL timestart("Green's Function: Postprocess")
CALL rot_projDOS(sym,atoms,input,angle,greensfCoeffs)
CALL rot_projDOS(sym,atoms,gfinp,input,angle,greensfCoeffs)
!Perform the Kramer-Kronigs-Integration if we only have the imaginary part at this point
CALL calc_onsite(atoms,input,sym,noco,angle,greensfCoeffs,greensf)
CALL calc_onsite(atoms,gfinp,input,sym,noco,results%ef,angle,greensfCoeffs,greensf)
!-------------------------------------------------------------
! 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)) THEN
CALL crystal_field(atoms,input,noco,greensfCoeffs,hub1,vTot)
IF(atoms%n_hia.GT.0.AND.ANY(ABS(hub1inp%ccf(:)).GT.1e-12)) THEN
CALL crystal_field(atoms,gfinp,hub1inp,input,noco,greensfCoeffs,vTot,results%ef,hub1data)
ENDIF
IF(input%jspins.EQ.2) THEN
CALL eff_excinteraction(greensf,atoms,input,results%ef,greensfCoeffs)
CALL eff_excinteraction(greensf,gfinp,input,results%ef,greensfCoeffs)
ENDIF
CALL timestart("Green's Function: Occupation/DOS")
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
DO i_gf = 1, gfinp%n
l = gfinp%elem(i_gf)%l
nType = gfinp%elem(i_gf)%atomType
IF(l.NE.gfinp%elem(i_gf)%lp) CYCLE
IF(nType.NE.gfinp%elem(i_gf)%atomTypep) CYCLE
!Density of states from Greens function
CALL gfDOS(greensf,l,nType,i_gf,atoms,input,results%ef)
CALL gfDOS(greensf,l,nType,i_gf,gfinp,input,results%ef)
!Occupation matrix