...
 
Commits (148)
......@@ -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,coreSpecInput,mcd,slab,orbcomp)
moments,hub1,coreSpecInput,mcd,slab,orbcomp,greensfCoeffs,angle,ntria,as,itria,atr)
!************************************************************************************
! This is the FLEUR valence density generator
......@@ -38,6 +38,10 @@ 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_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
......@@ -75,10 +79,17 @@ SUBROUTINE cdnval(eig_id, mpi,kpts,jspin,noco,input,banddos,cell,atoms,enpara,st
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_coreSpecInput), OPTIONAL, INTENT(IN) :: coreSpecInput
TYPE(t_mcd), OPTIONAL, INTENT(INOUT) :: mcd
TYPE(t_slab), OPTIONAL, INTENT(INOUT) :: slab
TYPE(t_orbcomp), OPTIONAL, INTENT(INOUT) :: orbcomp
TYPE(t_greensfCoeffs), OPTIONAL, INTENT(INOUT) :: greensfCoeffs
INTEGER, OPTIONAL, INTENT(IN) :: ntria
REAL, OPTIONAL, INTENT(IN) :: as
INTEGER, OPTIONAL, INTENT(IN) :: itria(3,2*kpts%nkpt)
REAL, OPTIONAL, INTENT(IN) :: atr(2*kpts%nkpt)
REAL, OPTIONAL, INTENT(IN) :: angle(sym%nop)
! Scalar Arguments
INTEGER, INTENT(IN) :: eig_id, jspin
......@@ -97,6 +108,8 @@ SUBROUTINE cdnval(eig_id, mpi,kpts,jspin,noco,input,banddos,cell,atoms,enpara,st
REAL :: we(MAXVAL(cdnvalJob%noccbd(:kpts%nkpt)))
REAL :: eig(MAXVAL(cdnvalJob%noccbd(:kpts%nkpt)))
REAL, ALLOCATABLE :: f(:,:,:,:),g(:,:,:,:),flo(:,:,:,:) ! radial functions
REAL, ALLOCATABLE :: tetweights(:,:)
INTEGER, ALLOCATABLE :: tet_ind(:,:)
TYPE (t_lapw) :: lapw
TYPE (t_orb) :: orb
......@@ -110,10 +123,11 @@ SUBROUTINE cdnval(eig_id, mpi,kpts,jspin,noco,input,banddos,cell,atoms,enpara,st
CALL timestart("cdnval")
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) THEN
IF (noco%l_mperp.OR.input%l_gfmperp) 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
......@@ -124,10 +138,12 @@ SUBROUTINE cdnval(eig_id, mpi,kpts,jspin,noco,input,banddos,cell,atoms,enpara,st
jsp_start = jspin
jsp_end = jspin
END IF
ALLOCATE (f(atoms%jmtd,2,0:atoms%lmaxd,jsp_start:jsp_end)) ! Deallocation before mpi_col_den
ALLOCATE (g(atoms%jmtd,2,0:atoms%lmaxd,jsp_start:jsp_end))
ALLOCATE (flo(atoms%jmtd,2,atoms%nlod,input%jspins))
ALLOCATE(tetweights(greensfCoeffs%ne,dimension%neigd))
ALLOCATE(tet_ind(dimension%neigd,2))
! Initializations
CALL usdus%init(atoms,input%jspins)
......@@ -162,7 +178,7 @@ SUBROUTINE cdnval(eig_id, mpi,kpts,jspin,noco,input,banddos,cell,atoms,enpara,st
DO ispin = jsp_start, jsp_end
CALL genMTBasis(atoms,enpara,vTot,mpi,iType,ispin,usdus,f(:,:,0:,ispin),g(:,:,0:,ispin),flo(:,:,:,ispin))
END DO
IF (noco%l_mperp) CALL denCoeffsOffdiag%addRadFunScalarProducts(atoms,f,g,flo,iType)
IF (noco%l_mperp.OR.input%l_gfmperp) CALL denCoeffsOffdiag%addRadFunScalarProducts(atoms,f,g,flo,iType)
IF (banddos%l_mcd) CALL mcd_init(atoms,input,dimension,vTot%mt(:,0,:,:),g,f,mcd,iType,jspin)
IF (l_coreSpec) CALL corespec_rme(atoms,input,iType,dimension%nstd,input%jspins,jspin,results%ef,&
dimension%msh,vTot%mt(:,0,:,:),f,g)
......@@ -225,6 +241,19 @@ 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) THEN
IF(input%film) THEN
CALL timestart("TriangularWeights")
tetweights = 0.0
CALL tria_weights(ikpt,kpts,ntria,as,itria,atr,results%neig(:,jsp),results%eig(:,:,jsp),greensfCoeffs,tetweights,tet_ind,results%ef)
CALL timestop("TriangularWeights")
ELSE
CALL timestart("TetrahedronWeights")
tetweights = 0.0
CALL tetra_weights(ikpt,kpts,results%neig(:,jsp),results%eig(:,:,jsp),greensfCoeffs,tetweights,tet_ind,results%ef)
CALL timestop("TetrahedronWeights")
ENDIF
ENDIF
DO ispin = jsp_start, jsp_end
IF (input%l_f) CALL force%init2(noccbd,input,atoms)
CALL abcof(input,atoms,sym,cell,lapw,noccbd,usdus,noco,ispin,oneD,&
......@@ -232,6 +261,18 @@ 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 onsite_coeffs(atoms,input,ispin,noccbd,tetweights,tet_ind,kpts%wtkpt(ikpt),eig,usdus,eigVecCoeffs,greensfCoeffs)
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
! perform Brillouin zone integration and summation over the
! bands in order to determine the energy parameters for each atom and angular momentum
CALL eparas(ispin,atoms,noccbd,mpi,ikpt,noccbd,we,eig,&
......@@ -239,6 +280,7 @@ SUBROUTINE cdnval(eig_id, mpi,kpts,jspin,noco,input,banddos,cell,atoms,enpara,st
IF (noco%l_mperp.AND.(ispin==jsp_end)) CALL qal_21(dimension,atoms,input,noccbd,noco,eigVecCoeffs,denCoeffsOffdiag,ikpt,dos)
! layer charge of each valence state in this k-point of the SBZ from the mt-sphere region of the film
IF (l_dosNdir) THEN
IF (PRESENT(slab)) CALL q_mt_sl(ispin,atoms,noccbd,ikpt,noccbd,skip_t,noccbd,eigVecCoeffs,usdus,slab)
......@@ -267,12 +309,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,&
results,denCoeffs,orb,denCoeffsOffdiag,den,den%mmpMat(:,:,:,jspin),mcd,slab,orbcomp)
results,denCoeffs,orb,denCoeffsOffdiag,den,den%mmpMat(:,:,:,jspin),mcd,slab,orbcomp,greensfCoeffs)
END DO
#endif
CALL cdnmt(mpi,input%jspins,atoms,sphhar,noco,jsp_start,jsp_end,&
enpara,vTot%mt(:,0,:,:),denCoeffs,usdus,orb,denCoeffsOffdiag,moments,den%mt)
enpara,vTot%mt(:,0,:,:),denCoeffs,usdus,orb,denCoeffsOffdiag,moments,den%mt,hub1)
IF (mpi%irank==0) THEN
IF (l_coreSpec) CALL corespec_ddscs(jspin,input%jspins)
DO ispin = jsp_start,jsp_end
......
......@@ -11,7 +11,7 @@ MODULE m_cdnmt
!***********************************************************************
CONTAINS
SUBROUTINE cdnmt(mpi,jspd,atoms,sphhar,noco,jsp_start,jsp_end,enpara,&
vr,denCoeffs,usdus,orb,denCoeffsOffdiag,moments,rho)
vr,denCoeffs,usdus,orb,denCoeffsOffdiag,moments,rho,hub1)
use m_constants,only: sfp_const
USE m_rhosphnlo
USE m_radfun
......@@ -26,6 +26,7 @@ CONTAINS
TYPE(t_atoms), INTENT(IN) :: atoms
TYPE(t_enpara), INTENT(IN) :: enpara
TYPE(t_moments), INTENT(INOUT) :: moments
TYPE(t_hub1ham), OPTIONAL, INTENT(INOUT) :: hub1
! .. Scalar Arguments ..
INTEGER, INTENT (IN) :: jsp_start,jsp_end,jspd
......@@ -39,12 +40,13 @@ CONTAINS
! ..
! .. Local Scalars ..
INTEGER itype,na,nd,l,lp,llp ,lh,j,ispin,noded,nodeu
INTEGER ilo,ilop,i
INTEGER ilo,ilop,i,i_hia,i_exc
REAL s,wronk,sumlm,qmtt
COMPLEX cs
LOGICAL l_hia
! ..
! .. Local Arrays ..
REAL qmtl(0:atoms%lmaxd,jspd,atoms%ntype),qmtllo(0:atoms%lmaxd)
REAL qmtl(0:atoms%lmaxd,jspd,atoms%ntype),qmtllo(0:atoms%lmaxd),vrav(atoms%jmtd)
CHARACTER(LEN=20) :: attributes(6)
! ..
......@@ -64,10 +66,10 @@ CONTAINS
ENDIF
!$OMP PARALLEL DEFAULT(none) &
!$OMP SHARED(usdus,rho,moments,qmtl) &
!$OMP SHARED(usdus,rho,moments,qmtl,hub1) &
!$OMP SHARED(atoms,jsp_start,jsp_end,enpara,vr,denCoeffs,sphhar)&
!$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,cs)
!$OMP PRIVATE(itype,na,ispin,l,rho21,f,g,nodeu,noded,wronk,i,j,s,qmtllo,qmtt,nd,lh,lp,llp,cs,l_hia,vrav,i_hia)
IF (noco%l_mperp) THEN
ALLOCATE ( f(atoms%jmtd,2,0:atoms%lmaxd,jspd),g(atoms%jmtd,2,0:atoms%lmaxd,jspd) )
ELSE
......@@ -86,8 +88,24 @@ CONTAINS
!---> spherical component
DO ispin = jsp_start,jsp_end
DO l = 0,atoms%lmax(itype)
CALL radfun(l,itype,ispin,enpara%el0(l,itype,ispin),vr(1,itype,ispin),atoms,&
f(1,1,l,ispin),g(1,1,l,ispin),usdus, nodeu,noded,wronk)
!Check if the orbital is to be treated with Hubbard 1
l_hia=.FALSE.
DO i = 1, atoms%n_hia
IF(atoms%lda_hia(i)%atomType.EQ.itype.AND.atoms%lda_hia(i)%l.EQ.l) THEN
l_hia=.TRUE.
ENDIF
ENDDO
IF(l_hia.AND.jspd.EQ.2) THEN
!In the case of a spin-polarized calculation with Hubbard 1 we want to treat
!the correlated orbitals with a non-spin-polarized basis
vrav = (vr(:,itype,1) + vr(:,itype,2))/2.0
CALL radfun(l,itype,ispin,enpara%el0(l,itype,ispin),vrav,atoms,&
f(1,1,l,ispin),g(1,1,l,ispin),usdus, nodeu,noded,wronk)
ELSE
CALL radfun(l,itype,ispin,enpara%el0(l,itype,ispin),vr(1,itype,ispin),atoms,&
f(1,1,l,ispin),g(1,1,l,ispin),usdus, nodeu,noded,wronk)
ENDIF
DO j = 1,atoms%jri(itype)
s = denCoeffs%uu(l,itype,ispin)*( f(j,1,l,ispin)*f(j,1,l,ispin)+f(j,2,l,ispin)*f(j,2,l,ispin) )&
+ denCoeffs%dd(l,itype,ispin)*( g(j,1,l,ispin)*g(j,1,l,ispin)+g(j,2,l,ispin)*g(j,2,l,ispin) )&
......@@ -122,6 +140,16 @@ CONTAINS
END DO
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
DO i_hia = 1, atoms%n_hia
IF(atoms%lda_hia(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)
ENDDO
ENDDO
ENDIF
!+soc
!---> spherical angular component
IF (noco%l_soc) THEN
......@@ -152,6 +180,7 @@ CONTAINS
ENDDO
ENDDO ! end of spin loop (ispin = jsp_start,jsp_end)
IF (noco%l_mperp) THEN
!---> calculate off-diagonal integrated density
......
......@@ -16,6 +16,7 @@ include("cmake/tests/test_GPU.cmake")
include("cmake/tests/test_MAGMA.cmake")
include("cmake/tests/test_LibXC.cmake")
include("cmake/tests/test_ELPA_onenode.cmake")
include("cmake/tests/test_EDSOLVER.cmake")
if (FLEUR_USE_MPI)
include("cmake/tests/test_SCALAPACK.cmake")
include("cmake/tests/test_ELPA.cmake")
......
......@@ -40,6 +40,9 @@ include(wannier/uhu/CMakeLists.txt)
include(forcetheorem/CMakeLists.txt)
include(rdmft/CMakeLists.txt)
include(kpoints/CMakeLists.txt)
include(ldahia/CMakeLists.txt)
include(intersite/CMakeLists.txt)
include(greensf/CMakeLists.txt)
include(tests/tests_old.cmake)
......@@ -63,7 +66,7 @@ inpgen/closure.f90 math/intgr.F90
io/w_inpXML.f90 kpoints/julia.f90 global/utility.F90
init/compile_descr.F90 kpoints/kpoints.f90 io/xmlOutput.F90 kpoints/brzone2.f90 cdn/slab_dim.f90 cdn/slabgeom.f90 dos/nstm3.f90 cdn/int_21.f90
cdn/int_21lo.f90 cdn_mt/rhomt21.f90 cdn_mt/rhonmt21.f90 force/force_a21.F90 force/force_a21_lo.f90 force/force_a21_U.f90 force/force_a12.f90
eigen/tlmplm_store.F90 xc-pot/gaunt.f90 kpoints/unfoldBandKPTS.f90)
eigen/tlmplm_store.F90 xc-pot/gaunt.f90 kpoints/unfoldBandKPTS.f90 math/ExpSave.f90 greensf/ind_greensf.f90)
set(fleur_SRC ${fleur_F90} ${fleur_F77})
......
......@@ -25,6 +25,7 @@ message("${Green} HDF5 Library found : ${CReset} ${FLEUR_USE_HDF5
message("${Green} Wannier90 1.2 Library found : ${CReset} ${FLEUR_USE_WANN}")
message("${Green} Wannier90-4 Library found : ${CReset} ${FLEUR_USE_WANN4}")
message("${Green} Wannier90-5 Library found : ${CReset} ${FLEUR_USE_WANN5}")
message("${Green} EDSolver Library found : ${CReset} ${FLEUR_USE_EDSOLVER}")
message("${Green} MAGMA Library found : ${CReset} ${FLEUR_USE_MAGMA}")
message("${Green} MPI Library found : ${CReset} ${FLEUR_USE_MPI}")
if (FLEUR_USE_MPI)
......
......@@ -18,6 +18,7 @@ do
-external) shift;external_lib="$external_lib $1";;
-hdf5) shift; CLI_USE_HDF5=$1;;
-wannier) shift; CLI_USE_WANNIER=$1;;
-edsolver) shift; CLI_USE_EDSOLVER=$1;;
-mpi) shift; CLI_USE_MPI=$1;;
-magma) shift; CLI_USE_MAGMA=$1;;
-gpu) shift; CLI_USE_GPU=$1;;
......@@ -83,10 +84,11 @@ General options:
Command line options to switch on/off features. These options overwrite the results of
the test and might lead to the configuration to fail.
-hdf5 [TRUE|FALSE] : use HDF5 library
-wannier [TRUE|FALSE] : use Wannier90 library
-mpi [TRUE|FALSE] : compile the MPI parallel version
-libxc [TRUE|FALSE] : use libxc library
-hdf5 [TRUE|FALSE] : use HDF5 library
-wannier [TRUE|FALSE] : use Wannier90 library
-mpi [TRUE|FALSE] : compile the MPI parallel version
-libxc [TRUE|FALSE] : use libxc library
-edsolver [TRUE|FALSE] : use the Exact Diagonalization library by Jindrich Kolorenc
Command line option to compile external libraries:
-external # : download and compile external libraries before building FLEUR
......
......@@ -86,6 +86,11 @@ then
echo "set(CLI_FLEUR_USE_WANNIER $CLI_USE_WANNIER)" >>config.cmake
fi
if [ "$CLI_USE_EDSOLVER" ]
then
echo "set(CLI_FLEUR_USE_EDSOLVER $CLI_USE_EDSOLVER)" >>config.cmake
fi
if [ "$CLI_USE_CHASE" ]
then
echo "set(CLI_FLEUR_USE_CHASE $CLI_USE_CHASE)" >>config.cmake
......
#first try if EdSolver already works
try_compile(FLEUR_USE_EDSOLVER ${CMAKE_BINARY_DIR} ${CMAKE_SOURCE_DIR}/cmake/tests/test_EDSOLVER.f90
LINK_LIBRARIES ${FLEUR_LIBRARIES})
message("EDSolver Library found:${FLEUR_USE_EDSOLVER}")
if (DEFINED CLI_FLEUR_USE_EDSOLVER)
if (${CLI_FLEUR_USE_EDSOLVER})
if (NOT FLEUR_USE_EDSOLVER)
message(FATAL ERROR "You explicitly asked for EDSolver but cmake couldn't find it")
endif()
else()
if (FLEUR_USE_EDSOLVER)
message("EDSolver library found, but you explicitely asked not to use it")
set(FLEUR_USE_EDSOLVER FALSE)
endif()
endif()
endif()
if (FLEUR_USE_EDSOLVER)
set(FLEUR_DEFINITIONS ${FLEUR_DEFINITIONS} "CPP_EDSOLVER")
set(FLEUR_MPI_DEFINITIONS ${FLEUR_MPI_DEFINITIONS} "CPP_EDSOLVER")
endif()
! ==========================================================================
! Example calls of the exact-diagonalization DMFT solver. The files
! selfen_from_bathGF_diag.dat and selfen_from_cfg.dat produced by this
! program should be identical. See module EDsolver for details.
! --------------------------------------------------------------------------
! Copyright (c) J. Kolorenc <kolorenc@fzu.cz> 2013, 2014
! All rights reserved.
! ==========================================================================
program EDdemo
use EDsolver, only: EDsolver_from_bathGF_diag, &
EDsolver_from_bathGF_offdiag, EDsolver_from_cfg
implicit none
integer, parameter :: dp=kind(1.0d0)
integer, parameter :: dpc=kind((1.0d0,1.0d0))
integer, parameter :: Ndos=1001
integer, parameter :: maxNbpar=3
integer, parameter :: bathMethod=1 ! large-z assymptotics
integer, parameter :: rediagonalize=0 ! no diagonalization in the second call
real(dp), parameter :: Emin =-15.0_dp
real(dp), parameter :: Emax = 5.0_dp
real(dp), parameter :: eps = 0.1_dp
real(dp), parameter :: eDC = 58.0_dp
complex(dpc), dimension(:,:,:), allocatable :: G0, G0save, selfen
complex(dpc), dimension(:), allocatable :: z, wt
real(dp), dimension(:,:), allocatable :: re, im
real(dp), dimension(:,:), allocatable :: bpar
integer, dimension(:), allocatable :: Nbpar
real(dp) :: a, b, c, d
integer :: Ncirc, dim
integer :: i
real(dp) :: time0, time1
real(dp), dimension(3) :: time
! precomputed bath Green's function on a semicircle in the complex plane
open(unit=11,file="G0circ32.dat",action="read")
read(unit=11,fmt=*) Ncirc, dim
allocate( G0(dim,dim,Ncirc+Ndos), G0save(dim,dim,Ncirc), &
selfen(dim,dim,Ncirc+Ndos), z(Ncirc+Ndos), wt(Ncirc), &
Nbpar(dim), bpar(maxNbpar,dim), re(dim,dim), im(dim,dim) )
Nbpar=3
do i=1, Ncirc
read(unit=11,fmt=*) a, b, c, d
z(i)=cmplx(a,b,dpc)
wt(i)=cmplx(c,d,dpc)
read(unit=11,fmt=*) re
read(unit=11,fmt=*) im
G0(:,:,i)=cmplx(re,im,dpc)
end do
close(unit=11)
! the energies where we want the selfenergy to be computed
do i=1, Ndos
z(Ncirc+i)=cmplx( Emin + (Emax-Emin)/(Ndos-1)*(i-1), eps, dpc )
end do
! EDsolver_from_bathGF_diag() first constructs the discrete impurity model
! from the supplied bath Green's function and then calculates the selfenergy
G0save=G0(:,:,1:Ncirc) ! G0 is destroyed by the subroutine
call cpu_time(time0)
call EDsolver_from_bathGF_diag(dim,Ncirc+Ndos,z,G0,Ncirc,wt,1,1,bathMethod, &
0.0_dp,Nbpar,bpar,eDC,selfen)
call cpu_time(time1)
time(1)=time1-time0
call write_selfen(z(Ncirc+1:Ncirc+Ndos),selfen(:,:,Ncirc+1:Ncirc+Ndos), &
"selfen_from_bathGF_diag.dat")
! EDsolver_from_cfg() takes an existing impurity model from configuration
! files stored on the disk and then calculates the selfenergy
call cpu_time(time0)
call EDsolver_from_cfg( dim, Ndos, z(Ncirc+1:Ncirc+Ndos), &
selfen(:,:,Ncirc+1:Ncirc+Ndos), rediagonalize )
call cpu_time(time1)
time(2)=time1-time0
call write_selfen(z(Ncirc+1:Ncirc+Ndos),selfen(:,:,Ncirc+1:Ncirc+Ndos), &
"selfen_from_cfg.dat")
! EDsolver_from_bathGF_offdiag() first constructs the discrete impurity model
! from the supplied bath Green's function and then calculates the selfenergy
G0(:,:,1:Ncirc)=G0save ! restore the bath Green's function
call cpu_time(time0)
call EDsolver_from_bathGF_offdiag(dim,Ncirc+Ndos,z,G0,Ncirc,wt, &
Nbpar(1),bpar,eDC,0,selfen)
call cpu_time(time1)
time(3)=time1-time0
call write_selfen(z(Ncirc+1:Ncirc+Ndos),selfen(:,:,Ncirc+1:Ncirc+Ndos), &
"selfen_from_bathGF_offdiag.dat")
call write_selfen_full_mtrx(z(Ncirc+1:Ncirc+Ndos), &
selfen(:,:,Ncirc+1:Ncirc+Ndos), "selfen_full_mtrx.dat")
! from_cfg is faster because the diagonalization was not performed again,
! from_bathGF_offdiag is a lot slower because it does not explore any
! symmetry
print *
print '(1x,a,f9.4)', "time in EDsolver_from_bathGF_diag() :", time(1)
print '(1x,a,f9.4)', "time in EDsolver_from_cfg() :", time(2)
print '(1x,a,f9.4)', "time in EDsolver_from_bathGF_offdiag():", time(3)
contains
subroutine write_selfen(z,selfen,filename)
! {{{
complex(dpc), dimension(:), intent(in) :: z
complex(dpc), dimension(:,:,:), intent(in) :: selfen
character(len=*), intent(in) :: filename
integer :: i
open(unit=11,file=filename,action="write")
do i=1, size(z)
write(unit=11,fmt="(5f15.7)") real(z(i),dp), &
real(selfen(3,3,i),dp), aimag(selfen(3,3,i)), &
real(selfen(2,2,i),dp), aimag(selfen(2,2,i))
end do
close(unit=11)
! }}}
end subroutine write_selfen
subroutine write_selfen_full_mtrx(z,selfen,filename)
! {{{
complex(dpc), dimension(:), intent(in) :: z
complex(dpc), dimension(:,:,:), intent(in) :: selfen
character(len=*), intent(in) :: filename
integer :: i
open(unit=11,file=filename,action="write")
do i=1, size(z)
write(unit=11,fmt="(2f15.7)") real(z(i),dp), aimag(z(i))
write(unit=11,fmt="(10f15.7)") real(selfen(:,:,i),dp)
write(unit=11,fmt="(10f15.7)") aimag(selfen(:,:,i))
end do
close(unit=11)
! }}}
end subroutine write_selfen_full_mtrx
end program EDdemo
! Local variables:
! folded-file: t
! End:
......@@ -19,7 +19,8 @@ CONTAINS
!> The matrices generated and diagonalized here are of type m_mat as defined in m_types_mat.
!>@author D. Wortmann
SUBROUTINE eigen(mpi,stars,sphhar,atoms,xcpot,sym,kpts,DIMENSION,vacuum,input,&
cell,enpara,banddos,noco,oneD,hybrid,iter,eig_id,results,inden,v,vx)
cell,enpara,banddos,noco,oneD,hybrid,iter,eig_id,results,inden,v,vx,gOnsite,hub1)
#include"cpp_double.h"
USE m_types
......@@ -62,6 +63,8 @@ CONTAINS
TYPE(t_atoms),INTENT(IN) :: atoms
TYPE(t_potden),INTENT(IN) :: inden,vx
TYPE(t_potden),INTENT(INOUT) :: v !u_setup will modify the potential matrix
TYPE(t_greensf),INTENT(INOUT):: gOnsite
TYPE(t_hub1ham),INTENT(INOUT) :: hub1
#ifdef CPP_MPI
INCLUDE 'mpif.h'
......@@ -123,7 +126,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,inden,v,mpi,results,DIMENSION,td,ud)
CALL mt_setup(atoms,sym,sphhar,input,noco,enpara,gOnsite,hub1,inden,v,mpi,results,DIMENSION,td,ud)
neigBuffer = 0
results%neig = 0
......
......@@ -7,12 +7,13 @@
MODULE m_mt_setup
CONTAINS
SUBROUTINE mt_setup(atoms,sym,sphhar,input,noco,enpara,inden,vTot,mpi,results,DIMENSION,td,ud)
SUBROUTINE mt_setup(atoms,sym,sphhar,input,noco,enpara,gOnsite,hub1,inden,vTot,mpi,results,DIMENSION,td,ud)
USE m_types
USE m_usetup
USE m_tlmplm_cholesky
USE m_tlmplm_store
USE m_spnorb
USE m_hubbard1_setup
IMPLICIT NONE
TYPE(t_results),INTENT(INOUT):: results
TYPE(t_mpi),INTENT(IN) :: mpi
......@@ -23,16 +24,21 @@ CONTAINS
TYPE(t_sym),INTENT(IN) :: sym
TYPE(t_sphhar),INTENT(IN) :: sphhar
TYPE(t_atoms),INTENT(IN) :: atoms
TYPE(t_greensf),INTENT(IN) :: gOnsite
TYPE(t_potden),INTENT(IN) :: inDen
TYPE(t_potden),INTENT(INOUT) :: vTot
TYPE(t_tlmplm),INTENT(INOUT) :: td
TYPE(t_usdus),INTENT(INOUT) :: ud
TYPE(t_hub1ham),INTENT(INOUT) :: hub1
INTEGER:: jsp
IF ((atoms%n_u.GT.0)) THEN
CALL u_setup(sym,atoms,sphhar,input,enpara%el0(0:,:,:),inDen,vTot,mpi,results)
END IF
IF((atoms%n_hia.GT.0)) THEN
CALL hubbard1_setup(atoms,hub1,sym,mpi,noco,input,ud,vTot,gOnsite,results)
END IF
CALL timestart("tlmplm")
CALL td%init(DIMENSION%lmplmd,DIMENSION%lmd,atoms%ntype,atoms%lmaxd,atoms%llod,SUM(atoms%nlo),&
......@@ -42,13 +48,13 @@ CONTAINS
DO jsp=1,MERGE(4,input%jspins,noco%l_mtNocoPot)
!CALL tlmplm_cholesky(sphhar,atoms,DIMENSION,enpara, jsp,1,mpi,vTot%mt(:,0,1,jsp),input,vTot%mmpMat, td,ud)
CALL tlmplm_cholesky(sphhar,atoms,noco,enpara,jsp,jsp,mpi,vTot,input,td,ud)
IF (input%l_f) CALL write_tlmplm(td,vTot%mmpMat,atoms%n_u>0,jsp,jsp,input%jspins)
IF (input%l_f) CALL write_tlmplm(td,vTot%mmpMat,atoms%n_u+atoms%n_hia>0,jsp,jsp,input%jspins)
END DO
CALL timestop("tlmplm")
!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.)
CALL spnorb(atoms,noco,input,mpi,enpara,vTot%mt,ud,td%rsoc,.FALSE.,hub1)
END IF
END SUBROUTINE mt_setup
......
......@@ -36,6 +36,7 @@ CONTAINS
REAL :: temp
INTEGER i,l,l2,lamda,lh,lm,lmin,lmin0,lmp,lmpl,lmplm,lmx,lmxx,lp,info,in
INTEGER lp1,lpl ,mem,mems,mp,mu,nh,na,m,nsym,s,i_u,jspin1,jspin2
LOGICAL l_remove
ALLOCATE( dvd(0:atoms%lmaxd*(atoms%lmaxd+3)/2,0:sphhar%nlhd ))
ALLOCATE( dvu(0:atoms%lmaxd*(atoms%lmaxd+3)/2,0:sphhar%nlhd ))
......@@ -106,6 +107,18 @@ CONTAINS
CALL intgr3(x,atoms%rmsh(1,n),atoms%dx(n),atoms%jri(n),temp)
dvd(lpl,lh) = temp
END IF
l_remove=.FALSE.
DO i = 1, atoms%n_hia
IF(atoms%lda_hia(i)%atomType.EQ.n.AND.atoms%lda_hia(i)%l.EQ.l.AND.lp.EQ.l) THEN
l_remove=.TRUE.
ENDIF
ENDDO
IF(l_remove) THEN
uvu(lpl,lh) = 0.0
dvd(lpl,lh) = 0.0
uvd(lpl,lh) = 0.0
dvu(lpl,lh) = 0.0
ENDIF
END DO
END DO
END DO
......
......@@ -13,6 +13,7 @@ CONTAINS
USE m_tlmplm
USE m_types
USE m_gaunt, ONLY: gaunt2
USE m_constants, ONLY: lmaxU_const
IMPLICIT NONE
TYPE(t_mpi),INTENT(IN) :: mpi
TYPE(t_noco),INTENT(IN) :: noco
......@@ -52,7 +53,6 @@ CONTAINS
td%tuloulo(:,:,:,jsp) = CMPLX(0.0,0.0)
td%h_off=0.0
!$ call gaunt2(atoms%lmaxd)
!$OMP PARALLEL DO DEFAULT(NONE)&
......@@ -96,12 +96,12 @@ CONTAINS
END DO
END DO
ENDDO
!Include contribution from LDA+U
DO i_u=1,atoms%n_u
IF (n.NE.atoms%lda_u(i_u)%atomtype) CYCLE
!Include contribution from LDA+U and LDA+HIA (latter are behind LDA+U contributions)
DO i_u=1,atoms%n_u+atoms%n_hia
IF (n.NE.MERGE(atoms%lda_u(i_u)%atomType,atoms%lda_hia(i_u)%atomType,i_u.LE.atoms%n_u)) CYCLE
!Found a "U" for this atom type
l=atoms%lda_u(i_u)%l
lp=atoms%lda_u(i_u)%l
l=MERGE(atoms%lda_u(i_u)%l,atoms%lda_hia(i_u)%l,i_u.LE.atoms%n_u)
lp=MERGE(atoms%lda_u(i_u)%l,atoms%lda_hia(i_u)%l,i_u.LE.atoms%n_u)
DO m = -l,l
lm = l* (l+1) + m
DO mp = -lp,lp
......@@ -111,7 +111,8 @@ CONTAINS
ENDDO
ENDDO
END DO
!Now add diagonal contribution to matrices
IF (jsp<3) THEN
!Create Cholesky decomposition of local hamiltonian
......
......@@ -67,7 +67,7 @@ CONTAINS
! read from eigenvalue and -vector file
!
l_real=sym%invs.and..not.noco%l_noco.and..not.(noco%l_soc.and.atoms%n_u>0)
l_real=sym%invs.and..not.noco%l_noco.and..not.(noco%l_soc.and.atoms%n_u+atoms%n_hia>0)
zmat%l_real=l_real
zMat(1:input%jspins)%matsize1=lapw%nv(1:input%jspins)+atoms%nlotot
zmat%matsize2=dimension%neigd
......
......@@ -21,7 +21,7 @@ MODULE m_eigenso
!
CONTAINS
SUBROUTINE eigenso(eig_id,mpi,DIMENSION,stars,vacuum,atoms,sphhar,&
obsolete,sym,cell,noco,input,kpts,oneD,vTot,enpara,results)
obsolete,sym,cell,noco,input,kpts,oneD,vTot,enpara,results,hub1)
USE m_types
USE m_eig66_io, ONLY : read_eig,write_eig
......@@ -49,6 +49,7 @@ 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
#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.)
CALL spnorb( atoms,noco,input,mpi, enpara,vTot%mt,usdus,rsoc,.TRUE.,hub1)
!
......
......@@ -12,8 +12,9 @@ MODULE m_spnorb
! using the functions anglso and sgml.
!*********************************************************************
CONTAINS
SUBROUTINE spnorb(atoms,noco,input,mpi, enpara, vr, usdus, rsoc,l_angles)
SUBROUTINE spnorb(atoms,noco,input,mpi, enpara, vr, usdus, rsoc,l_angles,hub1)
USE m_sorad
USE m_constants, only : hartree_to_ev_const
USE m_types
IMPLICIT NONE
......@@ -25,13 +26,14 @@ 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
! ..
! ..
! .. Array Arguments ..
REAL, INTENT (IN) :: vr(:,0:,:,:) !(atoms%jmtd,0:sphhar%nlhd,atoms%ntype,input%jspins)
! ..
! .. Local Scalars ..
INTEGER is1,is2,jspin1,jspin2,l,l1,l2,m1,m2,n
INTEGER is1,is2,jspin1,jspin2,l,l1,l2,m1,m2,n,i_hia
LOGICAL, SAVE :: first_k = .TRUE.
! ..
......@@ -70,6 +72,16 @@ CONTAINS
ENDIF
ENDDO
!Read in SOC-parameter for shell with hubbard 1
IF(PRESENT(hub1).AND.mpi%irank.EQ.0) THEN
DO i_hia = 1, atoms%n_hia
IF(hub1%l_soc_given(i_hia)) CYCLE
n = atoms%lda_hia(i_hia)%atomType
l = atoms%lda_hia(i_hia)%l
hub1%xi(i_hia) = (rsoc%rsopp(n,l,1,1)+rsoc%rsopp(n,l,2,2))*hartree_to_ev_const
ENDDO
ENDIF
!DO some IO into out file
IF ((first_k).AND.(mpi%irank.EQ.0)) THEN
DO n = 1,atoms%ntype
......
if [ ! -r libEDsolver_FLEUR ]
then
if [ -r libEDsolver_FLEUR.tar.gz ]
then
tar xzf libEDsolver_FLEUR.tar.gz
else
echo "No valid version of EDsolver library packaged for FLEUR available"
fi
fi
if [ -r libEDsolver_FLEUR ]
then
cd libEDsolver_FLEUR
#First add ARPACK
if [ ! -r ARPACK ]
then
echo "ARPACK is needed for the EDsolver"
else
FLEUR_LIBDIR="$PWD/ARPACK $FLEUR_LIBDIR"
fi
#Now look at the EDsolver
edsolver_version="2019-Jun-10"
if [ ! -r libEDsolver.${edsolver_version} ]
then
if [ -r libEDsolver.${edsolver_version}.tar.gz ]
then
tar xzf libEDsolver.${edsolver_version}.tar.gz
cp default.mk libEDsolver.${edsolver_version}/make/
cd libEDsolver.${edsolver_version}.tar.gz
make lib
else
echo "libEDsolver not present"
fi
else
cp default.mk libEDsolver.${edsolver_version}/make/
cd libEDsolver.${edsolver_version}
if [ ! -r "libEDsolver.a" ]
then
make lib
fi
fi
#Store the installation location
FLEUR_LIBDIR="$PWD/ $FLEUR_LIBDIR"
FLEUR_INCLUDEDIR="$PWD/ $FLEUR_INCLUDEDIR"
fi
\ No newline at end of file
......@@ -33,7 +33,9 @@ CONTAINS
REAL :: wronk
LOGICAL :: l_write
LOGICAL :: l_write,l_hia
REAL vrav(atoms%jmtd)
INTEGER :: i
l_write=mpi%irank==0
!$ l_write = l_write .and. omp_get_num_threads()==1
......@@ -41,8 +43,24 @@ CONTAINS
IF (l_write) WRITE (6,FMT=8000) iType
DO l = 0,atoms%lmax(iType)
CALL radfun(l,iType,jspin,enpara%el0(l,iType,jspin),vTot%mt(:,0,iType,jspin),atoms,&
f(1,1,l),g(1,1,l),usdus,nodeu,noded,wronk)
!Check if the orbital is to be treated with Hubbard 1
l_hia=.FALSE.
DO i = 1, atoms%n_hia
IF(atoms%lda_hia(i)%atomType.EQ.itype.AND.atoms%lda_hia(i)%l.EQ.l) THEN
l_hia=.TRUE.
ENDIF
ENDDO
IF(l_hia.AND.SIZE(vTot%mt,4).GT.1) THEN
!In the case of a spin-polarized calculation with Hubbard 1 we want to treat
!the correlated orbitals with a non-spin-polarized basis
vrav = (vTot%mt(:,0,iType,1) + vTot%mt(:,0,iType,2))/2.0
CALL radfun(l,iType,jspin,enpara%el0(l,iType,jspin),vrav,atoms,&
f(1,1,l),g(1,1,l),usdus, nodeu,noded,wronk)
ELSE
CALL radfun(l,iType,jspin,enpara%el0(l,iType,jspin),vTot%mt(:,0,iType,jspin),atoms,&
f(1,1,l),g(1,1,l),usdus,nodeu,noded,wronk)
ENDIF
IF (l_write) THEN
WRITE (6,FMT=8010) l,enpara%el0(l,iType,jspin),usdus%us(l,iType,jspin),usdus%dus(l,iType,jspin),&
nodeu,usdus%uds(l,iType,jspin),usdus%duds(l,iType,jspin),noded,usdus%ddn(l,iType,jspin),wronk
......
set(fleur_F77 ${fleur_F77}
)
set(fleur_F90 ${fleur_F90}
greensf/greensfImag.f90
greensf/greensfImag21.f90
greensf/angles.f90
greensf/onsite.f90
greensf/kkintgr.f90
greensf/tetra_weights.f90
greensf/gfcalc.f90
greensf/j0.f90
greensf/kk_cutoff.f90
greensf/ind_greensf.f90
)
MODULE m_angles
CONTAINS
!Calculate the correct angles for the rotation in spin-space (needed for LDA+U and Greens functions with noco)
SUBROUTINE angles(sym,angle)
USE m_types
USE m_constants
IMPLICIT NONE
TYPE(t_sym), INTENT(IN) :: sym
REAL, INTENT(OUT) :: angle(sym%nop)
INTEGER iop,d,t
angle = 0.0
DO iop = 1, sym%nop
d = det(sym%mrot(:,:,iop))
t = (sym%mrot(1,1,iop)+sym%mrot(2,2,iop)+sym%mrot(3,3,iop)) * d
IF(t.EQ.-1) angle(iop) = 1.0
IF(t.EQ.0) angle(iop) = 2.0/3.0
IF(t.EQ.1) angle(iop) = 1.0/2.0
IF(t.EQ.2) angle(iop) = 1.0/3.0
IF(t.EQ.3) angle(iop) = 0.0
angle(iop) = d*angle(iop)*pi_const
ENDDO
END SUBROUTINE
INTEGER FUNCTION det(m)
INTEGER m(3,3)
det = m(1,1)*m(2,2)*m(3,3) + m(1,2)*m(2,3)*m(3,1) + &
m(2,1)*m(3,2)*m(1,3) - m(1,3)*m(2,2)*m(3,1) - &
m(2,3)*m(3,2)*m(1,1) - m(2,1)*m(1,2)*m(3,3)
END FUNCTION det
END MODULE m_angles
\ No newline at end of file
MODULE m_gfcalc
!------------------------------------------------------------------------------
!
! MODULE: m_gfcalc
!
!> @author
!> Henning Janßen
!
! DESCRIPTION:
!> This module contains various subroutines calculating properties from the
!> green's functions:
!> -calculates the effective exchange interaction from the onsite
!> green's function according to Condens. Matter 26 (2014) 476003 EQ.1
!> -calculates the occupation matrix
!> -calculates the crystal-field-contrbution for the local hamiltonian
!
! REVISION HISTORY:
! February 2019 - Initial Version
! March 2019 - Changed calculation of the onsite exchange matrix
!------------------------------------------------------------------------------
USE m_juDFT
USE m_types
USE m_constants
IMPLICIT NONE
CONTAINS
SUBROUTINE occmtx(g,l,nType,atoms,sym,input,mmpMat,lp,nTypep)
USE m_ind_greensf
!calculates the occupation of a orbital treated with DFT+HIA from the related greens function
!The Greens-function should already be prepared on a energy contour ending at e_fermi
!The occupation is calculated with:
!
! n^sigma_mm' = -1/2pi int^Ef dz (G^+(z)^sigma_mm'-G^-(z)^sigma_mm')
IMPLICIT NONE
TYPE(t_greensf), INTENT(IN) :: g
TYPE(t_atoms), INTENT(IN) :: atoms
TYPE(t_sym), INTENT(IN) :: sym
TYPE(t_input), INTENT(IN) :: input
COMPLEX, INTENT(OUT) :: mmpMat(-lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const,MERGE(3,input%jspins,input%l_gfmperp))
INTEGER, INTENT(IN) :: l
INTEGER, INTENT(IN) :: nType
INTEGER, OPTIONAL, INTENT(IN) :: lp
INTEGER, OPTIONAL, INTENT(IN) :: nTypep
INTEGER ind1,ind2,ipm,iz,ispin,m,mp,lp_loop
LOGICAL l_vertcorr
REAL re,imag
TYPE(t_mat) :: gmat
l_vertcorr = .false. !Enables/Disables a correction for the vertical parts of the rectangular contour
mmpMat(:,:,:) = CMPLX(0.0,0.0)
IF(.NOT.PRESENT(lp)) THEN
lp_loop = l
ELSE
lp_loop = lp
ENDIF
!REPLACE: input%jspins --> MERGE(3,input%jspins,input%l_gfmperp)
DO ispin = 1, MERGE(3,input%jspins,input%l_gfmperp)
DO ipm = 1, 2
!Integrate over the contour:
DO iz = 1, g%nz
!get the corresponding gf-matrix
CALL g%get_gf(gmat,atoms,input,iz,l,nType,ipm.EQ.2,spin=ispin,lp=lp,nTypep=nTypep)
ind1 = 0
DO m = -l, l
ind1 = ind1 + 1
ind2 = 0
DO mp = -lp_loop,lp_loop
ind2 = ind2 + 1
mmpMat(m,mp,ispin) = mmpMat(m,mp,ispin) - 1/(2.0*pi_const*ImagUnit) * (-1)**(ipm-1) * gmat%data_c(ind1,ind2) &
* MERGE(g%de(iz),conjg(g%de(iz)),ipm.EQ.1)
ENDDO
ENDDO
CALL gmat%free()
ENDDO
ENDDO
ENDDO
END SUBROUTINE occmtx
SUBROUTINE crystal_field(atoms,input,greensfCoeffs,hub1,vu)
!calculates the crystal-field matrix for the local hamiltonian
!In addition we calculate the onsite exchange splitting for the calculation of j0
USE m_kkintgr
USE m_ind_greensf
IMPLICIT NONE
!-Type Arguments
TYPE(t_greensfCoeffs), INTENT(IN) :: greensfCoeffs
TYPE(t_atoms), INTENT(IN) :: atoms
TYPE(t_input), INTENT(IN) :: input
TYPE(t_hub1ham), INTENT(INOUT) :: hub1
!-Array Arguments
COMPLEX, INTENT(IN) :: vu(-lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const,atoms%n_hia,input%jspins) !LDA+U potential (should be removed from h_loc)
!-Local Scalars
INTEGER i_gf,l,nType,jspin,m,mp,ie,i_hia,kkcut,spin_cut
!-Local Arrays
REAL :: h_loc(-lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const,atoms%n_hia,input%jspins)
REAL :: integrand(greensfCoeffs%ne), norm(input%jspins),tr_hloc(atoms%n_hia,input%jspins)
h_loc = 0.0
tr_hloc = 0.0
DO i_hia = 1, atoms%n_hia
l = atoms%lda_hia(i_hia)%l
nType = atoms%lda_hia(i_hia)%atomType
i_gf = ind_greensf(atoms,l,nType)
!Perform the integration
!
! \int_{E_b}^{E_c} dE E * N_LL'(E)
!
norm = 0.0
DO jspin = 1, input%jspins
spin_cut = MERGE(1,jspin,jspin.GT.2)
kkcut = greensfCoeffs%kkintgr_cutoff(i_gf,spin_cut,2)
DO m = -l, l
DO mp = -l, l
integrand = 0.0
DO ie = 1, kkcut
integrand(ie) = ((ie-1) * greensfCoeffs%del+greensfCoeffs%e_bot) * greensfCoeffs%projdos(ie,i_gf,m,mp,jspin)
ENDDO
h_loc(m,mp,i_hia,jspin) = trapz(integrand,greensfCoeffs%del,greensfCoeffs%ne)
ENDDO
!trace of the integrated E*projdos
tr_hloc(i_hia,jspin) = tr_hloc(i_hia,jspin) + h_loc(m,m,i_hia,jspin) - REAL(vu(m,m,i_hia,jspin))
ENDDO
ENDDO
!Average over spins
hub1%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,:))/REAL(input%jspins) &
- SUM(vu(m,mp,i_hia,:))/REAL(input%jspins)
ENDDO
hub1%ccfmat(i_hia,m,m) = hub1%ccfmat(i_hia,m,m) - SUM(tr_hloc(i_hia,:))/REAL(input%jspins*(2*l+1))
ENDDO
ENDDO
END SUBROUTINE crystal_field
SUBROUTINE ldosmtx(app,g,i_gf,atoms,sym,input)
!calculates the l-dos from the onsite green's function
!If mode = 2 this may not make sense
USE m_intgr
USE m_ExpSave
CHARACTER(len=*), INTENT(IN) :: app
TYPE(t_greensf), INTENT(IN) :: g
TYPE(t_atoms), INTENT(IN) :: atoms
TYPE(t_sym), INTENT(IN) :: sym
TYPE(t_input), INTENT(IN) :: input
INTEGER, INTENT(IN) :: i_gf
REAL :: dos(-lmaxU_const:lmaxU_const,g%nz,input%jspins)
REAL :: re(-lmaxU_const:lmaxU_const,g%nz,input%jspins)
INTEGER ispin,m,j,l,n,i,ipm
REAL imag
dos = 0.0
re = 0.0
IF(.NOT.input%l_gfsphavg) CALL juDFT_warn("Not implemented for radial dependence",calledby="ldosmtx")
!n(E) = 1/2pii * Tr(G^+-G^-)
l = atoms%gfelem(i_gf)%l
n = atoms%gfelem(i_gf)%atomType
DO ispin = 1, input%jspins
DO m = -l , l
DO ipm = 1, 2
DO j = 1, g%nz
imag = AIMAG(g%gmmpMat(j,i_gf,m,m,ispin,ipm))
dos(m,j,ispin) = dos(m,j,ispin) + (-1)**(ipm-1) * imag
re(m,j,ispin) = re(m,j,ispin) + REAL(g%gmmpMat(j,i_gf,m,m,ispin,ipm))
ENDDO
ENDDO
ENDDO
ENDDO
OPEN(1337,file="lDOS_up_" // TRIM(ADJUSTL(app)) // ".txt",action="write",status="replace")
DO i = 1, g%nz
WRITE(1337,"(9f15.8)") REAL(g%e(i)), SUM(dos(-l:l,i,1)), (dos(m,i,1), m = -lmaxU_const, lmaxU_const)
ENDDO
CLOSE(unit = 1337)
IF(input%jspins.EQ.2) THEN
OPEN(1337,file="lDOS_dwn_" // TRIM(ADJUSTL(app)) // ".txt",action="write",status="replace")
DO i = 1, g%nz
WRITE(1337,"(9f15.8)") REAL(g%e(i)),SUM(dos(-l:l,i,2)), (dos(m,i,2), m = -lmaxU_const, lmaxU_const)
ENDDO
CLOSE(unit = 1337)
ENDIF
OPEN(1337,file="Re_up_" // TRIM(ADJUSTL(app)) // ".txt",action="write",status="replace")
DO i = 1, g%nz
WRITE(1337,"(9f15.8)") REAL(g%e(i)), SUM(re(-l:l,i,1)), (re(m,i,1), m = -lmaxU_const, lmaxU_const)
ENDDO
CLOSE(unit = 1337)
IF(input%jspins.EQ.2) THEN
OPEN(1337,file="Re_dwn_" // TRIM(ADJUSTL(app)) // ".txt",action="write",status="replace")
DO i = 1, g%nz
WRITE(1337,"(9f15.8)") REAL(g%e(i)),SUM(re(-l:l,i,2)), (re(m,i,2), m = -lmaxU_const, lmaxU_const)
ENDDO
CLOSE(unit = 1337)
ENDIF
END SUBROUTINE
END MODULE m_gfcalc
\ No newline at end of file
MODULE m_greensfImag
USE m_juDFT
USE m_types
USE m_constants
CONTAINS
SUBROUTINE greensfImag(atoms,sym,input,ispin,nbands,tetweights,ind,wtkpt,eig,usdus,eigVecCoeffs,greensfCoeffs)
!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_greensfCoeffs), INTENT(INOUT) :: greensfCoeffs
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
REAL, 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
!-Local Scalars
LOGICAL l_zero
INTEGER i_gf,ib,ie,j,nType,natom,l,m,mp,lm,lmp,ilo,ilop,imat,it,is,isi
REAL weight,fac
COMPLEX im(greensfCoeffs%ne,-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)
!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))
im = 0.0
!Loop through bands
!$OMP PARALLEL DEFAULT(none) &
!$OMP SHARED(natom,l,nType,ispin,wtkpt,i_gf,nbands) &
!$OMP SHARED(atoms,im,input,eigVecCoeffs,usdus,greensfCoeffs,eig,tetweights,ind) &
!$OMP PRIVATE(ie,m,mp,lm,lmp,ilo,ilop,weight,ib,j,l_zero)
!$OMP DO
DO ib = 1, nbands
!Check wether there is a non-zero weight for the energy window
l_zero = .true.
IF(input%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.
ELSE
!HISTOGRAM METHOD: check if eigenvalue is inside the energy range
j = NINT((eig(ib)-greensfCoeffs%e_bot)/greensfCoeffs%del)+1
IF( (j.LE.greensfCoeffs%ne).AND.(j.GE.1) ) l_zero = .false.
END IF
IF(l_zero) CYCLE
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 ie = MERGE(ind(ib,1),j,input%tria), MERGE(ind(ib,2),j,input%tria)
!weight for the bz-integration including spin-degeneracy
weight = 2.0/input%jspins * MERGE(tetweights(ie,ib),wtkpt/greensfCoeffs%del,input%tria)
!
!Contribution from states
!
im(ie,m,mp,1) = im(ie,m,mp,1) - pi_const * 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
im(ie,m,mp,1) = im(ie,m,mp,1) - pi_const * 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,