...
 
Commits (377)
......@@ -13,6 +13,7 @@ cdn/int_21.f90
cdn/int_21lo.f90
cdn/m_perp.f90
cdn/n_mat.f90
cdn/n_mat21.f90
cdn/od_abvac.f90
cdn/prp_qfft_map.f90
cdn/pwden.F90
......
......@@ -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)
!************************************************************************************
! This is the FLEUR valence density generator
......@@ -34,10 +34,13 @@ SUBROUTINE cdnval(eig_id, mpi,kpts,jspin,noco,input,banddos,cell,atoms,enpara,st
USE m_qal21 ! off-diagonal part of partial charges
USE m_abcof
USE m_nmat ! calculate density matrix for LDA + U
USE m_nmat21
USE m_vacden
USE m_pwden
USE m_forcea8
USE m_checkdopall
USE m_tetrahedronInit
USE m_gfcalc
USE m_cdnmt ! calculate the density and orbital moments etc.
USE m_orbmom ! coeffd for orbital moments
USE m_qmtsl ! These subroutines divide the input%film into vacuum%layers
......@@ -75,10 +78,13 @@ 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
REAL, OPTIONAL, INTENT(IN) :: angle(:)
! Scalar Arguments
INTEGER, INTENT(IN) :: eig_id, jspin
......@@ -97,6 +103,9 @@ SUBROUTINE cdnval(eig_id, mpi,kpts,jspin,noco,input,banddos,cell,atoms,enpara,st
REAL,ALLOCATABLE :: we(:),eig(:)
INTEGER,ALLOCATABLE :: ev_list(:)
REAL, ALLOCATABLE :: f(:,:,:,:),g(:,:,:,:),flo(:,:,:,:) ! radial functions
REAL, ALLOCATABLE :: resWeights(:,:)
REAL, ALLOCATABLE :: dosWeights(:,:)
INTEGER, ALLOCATABLE :: dosBound(:,:)
TYPE (t_lapw) :: lapw
TYPE (t_orb) :: orb
......@@ -110,10 +119,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
......@@ -160,9 +170,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))
CALL genMTBasis(atoms,enpara,vTot,mpi,iType,ispin,usdus,f(:,:,0:,ispin),g(:,:,0:,ispin),flo(:,:,:,ispin),input%l_dftspinpol)
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)
......@@ -179,7 +189,7 @@ SUBROUTINE cdnval(eig_id, mpi,kpts,jspin,noco,input,banddos,cell,atoms,enpara,st
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)
ev_list=cdnvaljob%compact_ev_list(ikpt_i,banddos%dos.OR.input%l_gf)
noccbd = SIZE(ev_list)
we = cdnvalJob%weights(ev_list,ikpt)
eig = results%eig(ev_list,ikpt,jsp)
......@@ -218,12 +228,26 @@ SUBROUTINE cdnval(eig_id, mpi,kpts,jspin,noco,input,banddos,cell,atoms,enpara,st
! valence density in the atomic spheres
CALL eigVecCoeffs%init(input,DIMENSION,atoms,noco,jspin,noccbd)
IF (atoms%n_gf.GT.0.AND.(input%tria.OR.input%gfTet)) THEN
CALL timestart("TetrahedronWeights")
ALLOCATE(dosWeights(greensfCoeffs%ne,noccbd))
ALLOCATE(dosBound(noccbd,2))
ALLOCATE(resWeights(greensfCoeffs%ne,noccbd))
CALL tetrahedronInit(ikpt,kpts,input,noccbd,results%eig(ev_list,:,jsp),&
greensfCoeffs,results%ef,resWeights,dosWeights,dosBound)
CALL timestop("TetrahedronWeights")
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,&
eigVecCoeffs%acof(:,0:,:,ispin),eigVecCoeffs%bcof(:,0:,:,ispin),&
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_u.GT.0) CALL n_mat(atoms,sym,noccbd,usdus,ispin,we,eigVecCoeffs,den%mmpMat(:,:,:,ispin))
IF (atoms%n_u.GT.0.AND.noco%l_mperp.AND.(ispin==jsp_end)) CALL n_mat21(atoms,sym,noccbd,denCoeffsOffdiag,&
we,eigVecCoeffs,angle,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)
! perform Brillouin zone integration and summation over the
! bands in order to determine the energy parameters for each atom and angular momentum
......@@ -232,6 +256,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)
......@@ -249,6 +274,7 @@ SUBROUTINE cdnval(eig_id, mpi,kpts,jspin,noco,input,banddos,cell,atoms,enpara,st
IF(l_coreSpec) CALL corespec_dos(atoms,usdus,ispin,dimension%lmd,kpts%nkpt,ikpt,dimension%neigd,&
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 (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
......@@ -260,12 +286,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,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,input%l_dftspinpol)
IF (mpi%irank==0) THEN
IF (l_coreSpec) CALL corespec_ddscs(jspin,input%jspins)
DO ispin = jsp_start,jsp_end
......
MODULE m_nmat21
! ************************************************************
! This subroutine calculates the density matrix n^{s}_{m,m'}
! for a given atom 'n' and l-quantum number 'l'. The l's for
! all atoms are stored in lda_u(), if lda_u()<0, no +U is used.
! For details see Eq.(12) of Shick et al. PRB 60, 10765 (1999)
! Part of the LDA+U package G.B., Oct. 2000
! ************************************************************
CONTAINS
SUBROUTINE n_mat21(atoms,sym, ne,denCoeffsOffdiag,we,eigVecCoeffs,angle,n_mmp)
!
USE m_types
USE m_constants
IMPLICIT NONE
TYPE(t_sym),INTENT(IN) :: sym
TYPE(t_atoms),INTENT(IN) :: atoms
TYPE(t_eigVecCoeffs),INTENT(IN) :: eigVecCoeffs
TYPE(t_denCoeffsOffDiag), INTENT(IN) :: denCoeffsOffDiag
REAL, INTENT(IN) :: angle(sym%nop)
! ..
! .. Scalar Arguments ..
INTEGER, INTENT (IN) :: ne
! ..
! .. Array Arguments ..
REAL, INTENT (IN) :: we(:)!(dimension%neigd)
COMPLEX, INTENT (INOUT) :: n_mmp(-lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const,atoms%n_u)
! ..
! .. Local Scalars ..
COMPLEX c_0,phase
INTEGER i,j,k,l,m,mp,n,it,is,isi,natom,n_ldau,lp
INTEGER ilo,ilop,ll1,nn,lmp,lm,i_u,natomTemp
REAL fac
! ..
! .. Local Arrays ..
COMPLEX n_tmp(-3:3,-3:3),nr_tmp(-3:3,-3:3),d_tmp(-3:3,-3:3)
COMPLEX n1_tmp(-3:3,-3:3)
! ..
! calculate n_mat:
!
natom = 0
i_u = 1
DO n = 1,atoms%ntype
DO WHILE (i_u.LE.atoms%n_u)
IF (atoms%lda_u(i_u)%atomType.GT.n) EXIT
natomTemp = natom
n_tmp(:,:) = cmplx(0.0,0.0)
l = atoms%lda_u(i_u)%l
ll1 = (l+1)*l
DO nn = 1, atoms%neq(n)
natomTemp = natomTemp + 1
!
! prepare n_mat in local frame (in noco-calculations this depends
! also on alpha(n) and beta(n) )
!
DO m = -l,l
lm = ll1+m
DO mp = -l,l
lmp = ll1+mp
c_0 = cmplx(0.0,0.0)
DO i = 1,ne
c_0 = c_0 + we(i) * ( &
+ conjg(eigVecCoeffs%acof(i,lmp,natomTemp,2))*eigVecCoeffs%acof(i,lm,natomTemp,1) * denCoeffsOffdiag%uu21n(l,n) &
+ conjg(eigVecCoeffs%acof(i,lmp,natomTemp,2))*eigVecCoeffs%bcof(i,lm,natomTemp,1) * denCoeffsOffdiag%ud21n(l,n) &
+ conjg(eigVecCoeffs%bcof(i,lmp,natomTemp,2))*eigVecCoeffs%acof(i,lm,natomTemp,1) * denCoeffsOffdiag%du21n(l,n) &
+ conjg(eigVecCoeffs%bcof(i,lmp,natomTemp,2))*eigVecCoeffs%bcof(i,lm,natomTemp,1) * denCoeffsOffdiag%dd21n(l,n))
ENDDO
n_tmp(m,mp) = c_0
ENDDO
ENDDO
!
! add local orbital contribution (if there is one) (untested so far)
!
DO ilo = 1, atoms%nlo(n)
IF (atoms%llo(ilo,n).EQ.l) THEN
DO m = -l,l
lm = ll1+m
DO mp = -l,l
lmp = ll1+mp
c_0 = cmplx(0.0,0.0)
DO i = 1,ne
c_0 = c_0 + we(i) * ( &
conjg(eigVecCoeffs%acof(i,lmp,natomTemp,2))*eigVecCoeffs%ccof(m,i,ilo,natomTemp,1) * denCoeffsOffdiag%uulo21n(l,n) + &
conjg(eigVecCoeffs%ccof(mp,i,ilo,natomTemp,2))*eigVecCoeffs%acof(i,lm,natomTemp,1) * denCoeffsOffdiag%ulou21n(l,n) + &
conjg(eigVecCoeffs%bcof(i,lmp,natomTemp,2))*eigVecCoeffs%ccof(m,i,ilo,natomTemp,1) * denCoeffsOffdiag%dulo21n(l,n) + &
conjg(eigVecCoeffs%ccof(mp,i,ilo,natomTemp,2))*eigVecCoeffs%bcof(i,lm,natomTemp,1) * denCoeffsOffdiag%ulod21n(l,n))
ENDDO
DO ilop = 1, atoms%nlo(n)
IF (atoms%llo(ilop,n).EQ.l) THEN
DO i = 1,ne
c_0 = c_0 + we(i) * denCoeffsOffdiag%uloulop21n(ilo,ilop,n) *conjg(eigVecCoeffs%ccof(mp,i,ilop,natomTemp,2)) *eigVecCoeffs%ccof(m ,i,ilo ,natomTemp,1)
ENDDO
ENDIF
ENDDO
n_tmp(m,mp) = n_tmp(m,mp) + c_0
ENDDO
ENDDO
ENDIF
ENDDO
!
! n_mmp should be rotated by D_mm' ; compare force_a21
!
DO it = 1, sym%invarind(natomTemp)
fac = 1.0 / ( sym%invarind(natomTemp) * atoms%neq(n) )
is = sym%invarop(natomTemp,it)
isi = sym%invtab(is)
d_tmp(:,:) = cmplx(0.0,0.0)
DO m = -l,l
DO mp = -l,l
d_tmp(m,mp) = sym%d_wgn(m,mp,l,isi)
ENDDO
ENDDO
nr_tmp = matmul( transpose( conjg(d_tmp) ) , n_tmp)
n1_tmp = matmul( nr_tmp, d_tmp )
phase = exp((0.0,1.0)*angle(isi))
DO m = -l,l
DO mp = -l,l
n_mmp(m,mp,i_u) = n_mmp(m,mp,i_u) + conjg(n1_tmp(m,mp)) *fac*phase
ENDDO
ENDDO
ENDDO
ENDDO ! sum over equivalent
i_u = i_u +1
ENDDO
natom = natom + atoms%neq(n)
ENDDO ! loop over atom types
! do m=-l,l
! write(*,'(14f12.6)') (n_mmp21(m,mp,1),mp=-l,l)
! enddo
RETURN
END SUBROUTINE n_mat21
END MODULE m_nmat21
......@@ -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,l_dftspinpol)
use m_constants,only: sfp_const
USE m_rhosphnlo
USE m_radfun
......@@ -26,6 +26,8 @@ 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
LOGICAL, INTENT(IN) :: l_dftspinpol
! .. Scalar Arguments ..
INTEGER, INTENT (IN) :: jsp_start,jsp_end,jspd
......@@ -39,12 +41,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),vrTmp(atoms%jmtd)
CHARACTER(LEN=20) :: attributes(6)
! ..
......@@ -64,10 +67,10 @@ CONTAINS
ENDIF
!$OMP PARALLEL DEFAULT(none) &
!$OMP SHARED(usdus,rho,moments,qmtl) &
!$OMP SHARED(atoms,jsp_start,jsp_end,enpara,vr,denCoeffs,sphhar)&
!$OMP SHARED(usdus,rho,moments,qmtl,hub1) &
!$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,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,vrTmp,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 +89,26 @@ 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 treated with Hubbard 1
l_hia=.FALSE.
DO i = atoms%n_u+1, atoms%n_u+atoms%n_hia
IF(atoms%lda_u(i)%atomType.EQ.itype.AND.atoms%lda_u(i)%l.EQ.l) THEN
l_hia=.TRUE.
ENDIF
ENDDO
!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
vrTmp = (vr(:,itype,1) + vr(:,itype,2))/2.0
ELSE
vrTmp = vr(:,itype,ispin)
ENDIF
CALL radfun(l,itype,ispin,enpara%el0(l,itype,ispin),vrTmp,atoms,&
f(1,1,l,ispin),g(1,1,l,ispin),usdus, nodeu,noded,wronk)
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 +143,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_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)
ENDDO
ENDDO
ENDIF
!+soc
!---> spherical angular component
IF (noco%l_soc) THEN
......@@ -152,6 +183,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:
......@@ -52,7 +52,7 @@
! locals
INTEGER, PARAMETER :: lmax= 4, ned = 1301
INTEGER i,s,v,index,jspin,k,l,l1,l2,ln,n,nl,ntb,ntria,ntetra
INTEGER icore,qdim,n_orb,ncored
INTEGER icore,qdim,n_orb,ncored,jsp
REAL as,de,efermi,emax,emin,qmt,sigma,totdos,efermiPrev
REAL e_up,e_lo,e_test1,e_test2,fac,sumwei,dk,eFermiCorrection
LOGICAL l_tria,l_orbcomp,l_error
......@@ -146,7 +146,8 @@
qal(:,:,k) = 0.0
qval(:,:,k) = 0.0
ntb = max(ntb,results%neig(k,jspin))
jsp = MERGE(1,jspin,noco%l_noco)
ntb = max(ntb,results%neig(k,jsp))
IF (l_mcd) mcd_local(:,:,k) = RESHAPE(mcd%mcd(:,1:ncored,:,k,jspin),(/3*atoms%ntype*ncored,dimension%neigd/))
IF (.NOT.l_orbcomp) THEN
qal(1:lmax*atoms%ntype,:,k)=reshape(dos%qal(0:,:,:,k,jspin),(/lmax*atoms%ntype,size(dos%qal,3)/))
......@@ -159,10 +160,10 @@
qal(slab%nsld+1:2*slab%nsld,:,k) = slab%qmtsl(:,:,k,jspin)
ELSE
DO i = 1, 23
DO l = 1, results%neig(k,jspin)
DO l = 1, results%neig(k,jsp)
qal(i,l,k) = orbcomp%comp(l,i,n_orb,k,jspin)*orbcomp%qmtp(l,n_orb,k,jspin)/10000.
END DO
DO l = results%neig(k,jspin)+1, dimension%neigd
DO l = results%neig(k,jsp)+1, dimension%neigd
qal(i,l,k) = 0.0
END DO
END DO
......@@ -178,7 +179,7 @@
qal(lmax*atoms%ntype+3,n,k) = 0.0
ENDDO
ELSEIF ( banddos%vacdos .and. input%film ) THEN
DO i = 1,results%neig(k,jspin)
DO i = 1,results%neig(k,jsp)
DO v = 1,vacuum%nvac
DO l = 1,vacuum%layers
index = (l-1)*vacuum%nstars + (v-1)*(vacuum%nstars*vacuum%layers) + 1
......@@ -213,10 +214,10 @@
!
!---- > convert eigenvalues to ev and shift them by efermi
!
DO i = 1 , results%neig(k,jspin)
ev(i,k) = results%eig(i,k,jspin)*hartree_to_ev_const - efermi
DO i = 1 , results%neig(k,jsp)
ev(i,k) = results%eig(i,k,jsp)*hartree_to_ev_const - efermi
ENDDO
DO i = results%neig(k,jspin) + 1, dimension%neigd
DO i = results%neig(k,jsp) + 1, dimension%neigd
ev(i,k) = 9.9e+99
ENDDO
!
......@@ -282,6 +283,7 @@
67 CONTINUE ! tetrahedron-information read or created
ENDIF
ENDIF
!
IF ( .not.l_mcd ) THEN
ALLOCATE (g(ned,qdim))
......@@ -304,7 +306,7 @@
ELSE
write(*,*) efermi
CALL tetra_dos(lmax,atoms%ntype,dimension%neigd,ned,ntetra,kpts%nkpt,&
itetra,efermi,voltet,e,results%neig(:,jspin), ev,qal, g)
itetra,efermi,voltet,e,results%neig(:,jsp), ev,qal, g)
IF (input%jspins.EQ.1) g(:,:) = 2 * g(:,:)
ENDIF
ELSE
......@@ -313,10 +315,10 @@
!
IF ( .not.l_mcd ) THEN
CALL dos_bin(input%jspins,qdim,ned,emin,emax,dimension%neigd,kpts%nkpt,&
results%neig(:,jspin),kpts%wtkpt(1:kpts%nkpt),ev,qal, g)
results%neig(:,jsp),kpts%wtkpt(1:kpts%nkpt),ev,qal, g)
ELSE
CALL dos_bin(input%jspins,3*atoms%ntype*ncored,ned,emin,emax,ntb,kpts%nkpt,&
results%neig(:,jspin),kpts%wtkpt(1:kpts%nkpt),ev(1:ntb,1:kpts%nkpt), mcd_local(1:3*atoms%ntype*ncored,1:ntb,1:kpts%nkpt), g)
results%neig(:,jsp),kpts%wtkpt(1:kpts%nkpt),ev(1:ntb,1:kpts%nkpt), mcd_local(1:3*atoms%ntype*ncored,1:ntb,1:kpts%nkpt), g)
ENDIF
ENDIF
!
......@@ -353,6 +355,7 @@
!**** write out DOS
OPEN (18,FILE='DOS'//spin12(jspin))
IF (atoms%ntype >= 20) OPEN (1337,FILE="PROJDOS"//spin12(jspin))
DO i = 1 , ned
totdos = 0.0
......@@ -369,6 +372,7 @@
ELSE
WRITE (18,99001) e(i),totdos,g(i,lmax*atoms%ntype+1), &
g(i,lmax*atoms%ntype+2),g(i,lmax*atoms%ntype+3), (gpart(i,l),l=1,atoms%ntype)
WRITE (1337,99001) e(i),(g(i,l),l=1,lmax)
ENDIF
ELSEIF (n_orb == 0) THEN
DO nl = 1, slab%nsld
......@@ -383,6 +387,7 @@
ENDIF
ENDDO
CLOSE (18)
IF (atoms%ntype >= 20) CLOSE (1337)
ELSE
write(*,'(4f15.8)') ((mcd%e_mcd(n,jspin,i),n=1,atoms%ntype),i=1,ncored)
......
......@@ -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,hub1)
#include"cpp_double.h"
USE m_types
......@@ -60,7 +61,9 @@ CONTAINS
TYPE(t_kpts),INTENT(INOUT) :: kpts
TYPE(t_sphhar),INTENT(IN) :: sphhar
TYPE(t_atoms),INTENT(IN) :: atoms
TYPE(t_potden),INTENT(IN) :: inden,vx
TYPE(t_potden),INTENT(IN) :: inden !
TYPE(t_hub1ham),INTENT(INOUT):: hub1
TYPE(t_potden), INTENT(IN) :: vx
TYPE(t_potden),INTENT(INOUT) :: v !u_setup will modify the potential matrix
#ifdef CPP_MPI
......@@ -106,7 +109,7 @@ CONTAINS
ALLOCATE(bkpt(3))
ALLOCATE(eigBuffer(DIMENSION%neigd,kpts%nkpt,input%jspins))
l_real=sym%invs.AND..NOT.noco%l_noco
l_real=sym%invs.AND..NOT.noco%l_noco.AND..NOT.(noco%l_soc.AND.atoms%n_u+atoms%n_hia>0)
! check if z-reflection trick can be used
l_zref=(sym%zrfs.AND.(SUM(ABS(kpts%bk(3,:kpts%nkpt))).LT.1e-9).AND..NOT.noco%l_noco)
......@@ -122,7 +125,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,hub1,inden,v,mpi,results,DIMENSION,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,inden,vTot,mpi,results,DIMENSION,td,ud)
SUBROUTINE mt_setup(atoms,sym,sphhar,input,noco,enpara,hub1,inden,vTot,mpi,results,DIMENSION,td,ud)
USE m_types
USE m_usetup
USE m_tlmplm_cholesky
......@@ -27,13 +27,15 @@ 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
INTEGER:: jsp
IF ((atoms%n_u.GT.0)) THEN
CALL u_setup(sym,atoms,sphhar,input,enpara%el0(0:,:,:),inDen,vTot,mpi,results)
IF (atoms%n_u+atoms%n_hia>0) THEN
CALL u_setup(sym,atoms,sphhar,input,noco,enpara%el0(0:,:,:),inDen,vTot,mpi,results)
END IF
CALL timestart("tlmplm")
CALL td%init(DIMENSION%lmplmd,DIMENSION%lmd,atoms%ntype,atoms%lmaxd,atoms%llod,SUM(atoms%nlo),&
DOT_PRODUCT(atoms%nlo,atoms%nlo+1)/2,MERGE(4,input%jspins,noco%l_mtNocoPot),&
......@@ -42,13 +44,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
......
......@@ -31,18 +31,36 @@ MODULE m_radovlp
REAL, INTENT (OUT):: dun21(0:atoms%lmaxd,atoms%ntype),ddn21(0:atoms%lmaxd,atoms%ntype)
! ..
! .. Local Scalars ..
INTEGER itype,l,ispin,noded,nodeu
INTEGER itype,l,ispin,noded,nodeu,i
REAL wronk
LOGICAL l_hia
! ..
! .. Local Arrays ..
REAL :: f(atoms%jmtd,2,0:atoms%lmaxd,input%jspins)
REAL :: g(atoms%jmtd,2,0:atoms%lmaxd,input%jspins)
REAL :: vrTmp(atoms%jmtd)
! ..
DO itype = 1,atoms%ntype
DO l = 0,atoms%lmax(itype)
l_hia=.FALSE.
DO i = atoms%n_u+1, atoms%n_u+atoms%n_hia
IF(atoms%lda_u(i)%atomType.EQ.itype.AND.atoms%lda_u(i)%l.EQ.l) THEN
l_hia=.TRUE.
ENDIF
ENDDO
DO ispin = 1,input%jspins
CALL radfun(l,itype,ispin,epar(l,itype,ispin),vr(:,0,itype,ispin), atoms,&
f(1,1,l,ispin),g(1,1,l,ispin),usdus, nodeu,noded,wronk)
!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.input%jspins.EQ.2.AND..NOT.input%l_dftspinpol) THEN
vrTmp = (vr(:,0,iType,1) + vr(:,0,iType,2))/2.0
ELSE
vrTmp = vr(:,0,iType,ispin)
ENDIF
CALL radfun(l,iType,ispin,epar(l,iType,ispin),vrTmp,atoms,&
f(1,1,l,ispin),g(1,1,l,ispin),usdus,nodeu,noded,wronk)
ENDDO
CALL int_21(f,g,atoms,itype,l,uun21,udn21,dun21,ddn21)
ENDDO
......
......@@ -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 ))
......@@ -53,7 +54,7 @@ CONTAINS
IF (jsp<3) vr0(:,0)=0.0
DO i=MERGE(1,jspin,jspin>2),MERGE(2,jspin,jspin>2)
CALL genMTBasis(atoms,enpara,v,mpi,n,i,ud,f(:,:,:,i),g(:,:,:,i),flo)
CALL genMTBasis(atoms,enpara,v,mpi,n,i,ud,f(:,:,:,i),g(:,:,:,i),flo,input%l_dftspinpol)
ENDDO
IF (jspin>2) THEN
jspin1=1
......@@ -72,6 +73,15 @@ CONTAINS
lp1 = (lp* (lp+1))/2
DO l = 0,lp
lpl = lp1 + l
!----------------------------------------------------------------------------
! Remove non-spherical components for the orbitals treated with DFT+Hubbard-1
!----------------------------------------------------------------------------
l_remove=.FALSE.
IF(l.EQ.lp) THEN
DO i = atoms%n_u+1, atoms%n_u+atoms%n_hia
IF(atoms%lda_u(i)%atomType.EQ.n.AND.atoms%lda_u(i)%l.EQ.l) l_remove=.TRUE.
ENDDO
ENDIF
!---> loop over non-spherical components of the potential: must
!---> satisfy the triangular conditions and that l'+l+lamda even
!---> (conditions from the gaunt coefficient)
......@@ -79,7 +89,7 @@ CONTAINS
lamda = sphhar%llh(lh,nsym)
lmin = lp - l
lmx = lp + l
IF ((mod(lamda+lmx,2).EQ.1) .OR. (lamda.LT.lmin) .OR. (lamda.GT.lmx)) THEN
IF ((mod(lamda+lmx,2).EQ.1) .OR. (lamda.LT.lmin) .OR. (lamda.GT.lmx) .OR. l_remove) THEN
uvu(lpl,lh) = 0.0
dvd(lpl,lh) = 0.0
uvd(lpl,lh) = 0.0
......@@ -183,7 +193,7 @@ CONTAINS
!
!---> set up the t-matrices for the local orbitals,
!---> if there are any
IF (atoms%nlo(n).GE.1) THEN
IF (atoms%nlo(n).GE.1.AND.jspin<3) THEN
CALL tlo(atoms,sphhar,jspin,jsp,n,enpara,1,input,v%mt(1,0,n,jsp),&
na,flo,f(:,:,:,jspin),g(:,:,:,jspin),ud, ud%uuilon(:,:,jspin),ud%duilon(:,:,jspin),ud%ulouilopn(:,:,:,jspin), td)
......
......@@ -12,7 +12,9 @@ CONTAINS
jspin,jsp,mpi,v,input,td,ud)
USE m_tlmplm
USE m_types
USE m_radovlp
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
......@@ -36,7 +38,7 @@ CONTAINS
LOGICAL OK
! ..
! .. Local Arrays ..
REAL, ALLOCATABLE :: uun21(:,:),udn21(:,:),dun21(:,:),ddn21(:,:)
REAL,PARAMETER:: e_shift_min=0.5
REAL,PARAMETER:: e_shift_max=65.0
......@@ -46,20 +48,26 @@ CONTAINS
td%e_shift(:,jsp)=0.0
IF (jsp<3) td%e_shift(:,jsp)=e_shift_min
!Calculate overlap integrals
!For the off-diagonal LDA+U contributions
IF(jsp>=3) THEN
ALLOCATE(uun21(0:atoms%lmaxd,atoms%ntype),udn21(0:atoms%lmaxd,atoms%ntype),&
dun21(0:atoms%lmaxd,atoms%ntype),ddn21(0:atoms%lmaxd,atoms%ntype) )
CALL rad_ovlp(atoms,ud,input,v%mt,enpara%el0, uun21,udn21,dun21,ddn21)
ENDIF
td%tdulo(:,:,:,jsp) = CMPLX(0.0,0.0)
td%tuulo(:,:,:,jsp) = CMPLX(0.0,0.0)
td%tuloulo(:,:,:,jsp) = CMPLX(0.0,0.0)
td%h_off=0.0
!$ call gaunt2(atoms%lmaxd)
!$OMP PARALLEL DO DEFAULT(NONE)&
!$OMP PRIVATE(temp,i,l,lm,lmin,lmin0,lmp)&
!$OMP PRIVATE(lmplm,lp,m,mp,n)&
!$OMP PRIVATE(OK,s,in,info)&
!$OMP SHARED(atoms,jspin,jsp,sphhar,enpara,td,ud,v,mpi,input)
!$OMP SHARED(atoms,jspin,jsp,sphhar,enpara,td,ud,v,mpi,input,uun21,udn21,dun21,ddn21)
DO n = 1,atoms%ntype
CALL tlmplm(n,sphhar,atoms,enpara,jspin,jsp,mpi,v,input,td,ud)
OK=.FALSE.
......@@ -96,9 +104,9 @@ 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.atoms%lda_u(i_u)%atomType) CYCLE
!Found a "U" for this atom type
l=atoms%lda_u(i_u)%l
lp=atoms%lda_u(i_u)%l
......@@ -106,12 +114,30 @@ CONTAINS
lm = l* (l+1) + m
DO mp = -lp,lp
lmp = lp* (lp+1) + mp
td%h_loc(lm,lmp,n,jsp) =td%h_loc(lm,lmp,n,jsp) + v%mmpMat(m,mp,i_u,jsp)
td%h_loc(lm+s,lmp+s,n,jsp) =td%h_loc(lm+s,lmp+s,n,jsp)+ v%mmpMat(m,mp,i_u,jsp)*ud%ddn(lp,n,jsp)
!------------------------------------------------------------------------
! Currently for jsp >= 3 the convention is:
! -jsp=3 => real part of the off-diagonal hamiltonian
! -jsp=4 => imaginary part of the off-diagonal hamiltonian
!------------------------------------------------------------------------
IF (jsp < 3) THEN
td%h_loc(lm,lmp,n,jsp) =td%h_loc(lm,lmp,n,jsp) + v%mmpMat(m,mp,i_u,jsp)
td%h_loc(lm+s,lmp+s,n,jsp) =td%h_loc(lm+s,lmp+s,n,jsp) + v%mmpMat(m,mp,i_u,jsp) * ud%ddn(lp,n,jsp)
ELSE IF(jsp.EQ.3) THEN
td%h_loc(lm,lmp,n,jsp) =td%h_loc(lm,lmp,n,jsp) + REAL(v%mmpMat(m,mp,i_u,3)) * uun21(l,n)
td%h_loc(lm+s,lmp,n,jsp) =td%h_loc(lm+s,lmp,n,jsp) + REAL(v%mmpMat(m,mp,i_u,3)) * dun21(l,n)
td%h_loc(lm,lmp+s,n,jsp) =td%h_loc(lm,lmp+s,n,jsp) + REAL(v%mmpMat(m,mp,i_u,3)) * udn21(l,n)
td%h_loc(lm+s,lmp+s,n,jsp) =td%h_loc(lm+s,lmp+s,n,jsp) + REAL(v%mmpMat(m,mp,i_u,3)) * ddn21(l,n)
ELSE
td%h_loc(lm,lmp,n,jsp) =td%h_loc(lm,lmp,n,jsp) + AIMAG(v%mmpMat(m,mp,i_u,3)) * uun21(l,n)
td%h_loc(lm+s,lmp,n,jsp) =td%h_loc(lm+s,lmp,n,jsp) + AIMAG(v%mmpMat(m,mp,i_u,3)) * dun21(l,n)
td%h_loc(lm,lmp+s,n,jsp) =td%h_loc(lm,lmp+s,n,jsp) + AIMAG(v%mmpMat(m,mp,i_u,3)) * udn21(l,n)
td%h_loc(lm+s,lmp+s,n,jsp) =td%h_loc(lm+s,lmp+s,n,jsp) + AIMAG(v%mmpMat(m,mp,i_u,3)) * ddn21(l,n)
ENDIF
ENDDO
ENDDO
END DO
!Now add diagonal contribution to matrices
IF (jsp<3) THEN
!Create Cholesky decomposition of local hamiltonian
......@@ -199,9 +225,9 @@ CONTAINS
DO l=0,atoms%lnonsph(n)
c=(-0.5)*CMPLX(noco%b_con(1,n),noco%b_con(2,n))
td%h_off(l ,l ,n,1) =td%h_off(l ,l ,n,1) + uun21(l,n)*c
td%h_loc(l ,l+s,n,1) =td%h_off(l ,l+s,n,1) + udn21(l,n)*c
td%h_loc(l+s,l ,n,1) =td%h_off(l+s,l ,n,1) + dun21(l,n)*c
td%h_loc(l+s,l+s,n,1) =td%h_off(l+s,l+s,n,1) + ddn21(l,n)*c
td%h_off(l ,l+s,n,1) =td%h_off(l ,l+s,n,1) + udn21(l,n)*c
td%h_off(l+s,l ,n,1) =td%h_off(l+s,l ,n,1) + dun21(l,n)*c
td%h_off(l+s,l+s,n,1) =td%h_off(l+s,l+s,n,1) + ddn21(l,n)*c
ENDDO
......@@ -209,9 +235,9 @@ CONTAINS
DO l=0,atoms%lnonsph(n)
c=(-0.5)*CMPLX(noco%b_con(1,n),-noco%b_con(2,n))
td%h_off(l ,l ,n,2) =td%h_off(l ,l ,n,2) + uun21(l,n)*c
td%h_loc(l ,l+s,n,2) =td%h_off(l ,l+s,n,2) + udn21(l,n)*c
td%h_loc(l+s,l ,n,2) =td%h_off(l+s,l ,n,2) + dun21(l,n)*c
td%h_loc(l+s,l+s,n,2) =td%h_off(l+s,l+s,n,2) + ddn21(l,n)*c
td%h_off(l ,l+s,n,2) =td%h_off(l ,l+s,n,2) + udn21(l,n)*c
td%h_off(l+s,l ,n,2) =td%h_off(l+s,l ,n,2) + dun21(l,n)*c
td%h_off(l+s,l+s,n,2) =td%h_off(l+s,l+s,n,2) + ddn21(l,n)*c
ENDDO
END DO
END SUBROUTINE tlmplm_constrained
......
......@@ -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)
!
......
......@@ -35,38 +35,49 @@ CONTAINS
! .. Local Scalars ..
REAL ddn1,e ,ulops,dulops,duds1
INTEGER i,j,ir,jspin,l,noded,nodeu,ilo,ilop
LOGICAL l_hia
! ..
! .. Local Arrays ..
REAL, ALLOCATABLE :: p(:,:),pd(:,:),q(:,:),qd(:,:),plo(:,:)
REAL, ALLOCATABLE :: plop(:,:),glo(:,:),fint(:),pqlo(:,:)
REAL, ALLOCATABLE :: filo(:,:)
REAL, ALLOCATABLE :: v0(:),vso(:,:),qlo(:,:)
REAL, ALLOCATABLE :: v0(:),vso(:,:),qlo(:,:),vrTmp(:)
! ..
IF (atoms%jri(ntyp)>atoms%jmtd) CALL juDFT_error("atoms%jri(ntyp).GT.atoms%jmtd",calledby ="sorad")
ALLOCATE ( p(atoms%jmtd,2),pd(atoms%jmtd,2),q(atoms%jmtd,2),plo(atoms%jmtd,2),fint(atoms%jmtd),&
& qlo(atoms%jmtd,2),plop(atoms%jmtd,2),qd(atoms%jmtd,2),v0(atoms%jmtd),vso(atoms%jmtd,2) )
& qlo(atoms%jmtd,2),plop(atoms%jmtd,2),qd(atoms%jmtd,2),v0(atoms%jmtd),vso(atoms%jmtd,2),vrTmp(atoms%jmtd) )
p = 0.0 ; pd = 0.0 ; q = 0.0 ; plo = 0.0 ; fint = 0.0
qlo = 0.0 ; plop = 0.0 ; qd = 0.0 ; v0 = 0.0 ; vso = 0.0
qlo = 0.0 ; plop = 0.0 ; qd = 0.0 ; v0 = 0.0 ; vso = 0.0; vrTmp = 0.0
DO l = 0,atoms%lmax(ntyp)
l_hia=.FALSE.
DO i = atoms%n_u+1, atoms%n_u+atoms%n_hia
IF(atoms%lda_u(i)%atomType.EQ.ntyp.AND.atoms%lda_u(i)%l.EQ.l) THEN
l_hia=.TRUE.
ENDIF
ENDDO
DO jspin = 1,input%jspins
IF(l_hia.AND.input%jspins.EQ.2.AND..NOT.input%l_dftspinpol) THEN
vrTmp = (vr(:,1)+vr(:,2))/2.0
ELSE
vrTmp = vr(:,jspin)
ENDIF
!
!---> calculate normalized function at e: p and q
!
e = enpara%el0(l,ntyp,jspin)
CALL radsra(&
e,l,vr(:,jspin),atoms%rmsh(1,ntyp),atoms%dx(ntyp),atoms%jri(ntyp),c_light(1.0),&
e,l,vrTmp,atoms%rmsh(1,ntyp),atoms%dx(ntyp),atoms%jri(ntyp),c_light(1.0),&
usdus%us(l,ntyp,jspin),usdus%dus(l,ntyp,jspin),&
nodeu,p(:,jspin),q(:,jspin))