Commit ee8f47c8 authored by Matthias Redies's avatar Matthias Redies

Merge branch 'develop' into MetaGGA

parents ae409a52 088e558e
MODULE m_cored MODULE m_cored
CONTAINS CONTAINS
SUBROUTINE cored(input, jspin, atoms, rho, DIMENSION, sphhar, vr, qint, rhc, tec, seig, EnergyDen) SUBROUTINE cored(&
! ******************************************************* & input,jspin,atoms,&
! ***** set up the core densities for compounds. ***** & rho,DIMENSION,&
! ***** d.d.koelling ***** & sphhar,&
! ******************************************************* & vr,&
USE m_juDFT & qint,rhc,tec,seig)
USE m_intgr, ONLY : intgr3,intgr0,intgr1
USE m_constants, ONLY : c_light,sfp_const ! *******************************************************
USE m_setcor ! ***** set up the core densities for compounds. *****
USE m_differ ! ***** d.d.koelling *****
USE m_types ! *******************************************************
USE m_xmlOutput USE m_juDFT
IMPLICIT NONE USE m_intgr, ONLY : intgr3,intgr0,intgr1
TYPE(t_dimension),INTENT(IN) :: DIMENSION USE m_constants, ONLY : c_light,sfp_const
TYPE(t_input),INTENT(IN) :: input USE m_setcor
TYPE(t_sphhar),INTENT(IN) :: sphhar USE m_differ
TYPE(t_atoms),INTENT(IN) :: atoms USE m_types
! USE m_xmlOutput
! .. Scalar Arguments .. IMPLICIT NONE
INTEGER, INTENT (IN) :: jspin TYPE(t_dimension),INTENT(IN) :: DIMENSION
REAL, INTENT (OUT) :: seig TYPE(t_input),INTENT(IN) :: input
! .. TYPE(t_sphhar),INTENT(IN) :: sphhar
! .. Array Arguments .. TYPE(t_atoms),INTENT(IN) :: atoms
REAL, INTENT(IN) :: vr(atoms%jmtd,atoms%ntype) !
REAL, INTENT(INOUT) :: rho(atoms%jmtd,0:sphhar%nlhd,atoms%ntype,input%jspins) ! .. Scalar Arguments ..
REAL, INTENT(INOUT) :: rhc(DIMENSION%msh,atoms%ntype,input%jspins) INTEGER, INTENT (IN) :: jspin
REAL, INTENT(INOUT) :: qint(atoms%ntype,input%jspins) REAL, INTENT (OUT) :: seig
REAL, INTENT(INOUT) :: tec(atoms%ntype,input%jspins) ! ..
REAL, INTENT(INOUT), OPTIONAL :: EnergyDen(atoms%jmtd,0:sphhar%nlhd,atoms%ntype,input%jspins) ! .. Array Arguments ..
! .. REAL, INTENT(IN) :: vr(atoms%jmtd,atoms%ntype)
! .. Local Scalars .. REAL, INTENT(INOUT) :: rho(atoms%jmtd,0:sphhar%nlhd,atoms%ntype,input%jspins)
REAL eig,fj,fl,fn,q,rad,rhos,rhs,sea,sume,t2 REAL, INTENT(INOUT) :: rhc(DIMENSION%msh,atoms%ntype,input%jspins)
REAL d,dxx,rn,rnot,z,t1,rr,r,lambd,c,bmu,weight, aux_weight REAL, INTENT(INOUT) :: qint(atoms%ntype,input%jspins)
INTEGER i,j,jatom,korb,n,ncmsh,nm,nm1,nst ,l,ierr REAL, INTENT(INOUT) :: tec(atoms%ntype,input%jspins)
! .. ! ..
! .. Local Arrays .. ! .. Local Scalars ..
REAL e,fj,fl,fn,q,rad,rhos,rhs,sea,sume,t2
REAL rhcs(DIMENSION%msh),rhoc(DIMENSION%msh),rhoss(DIMENSION%msh),vrd(DIMENSION%msh),f(0:3) REAL d,dxx,rn,rnot,z,t1,rr,r,lambd,c,bmu,weight
REAL rhcs_aux(DIMENSION%msh), rhoss_aux(DIMENSION%msh) !> quantities for energy density calculations INTEGER i,j,jatom,korb,n,ncmsh,nm,nm1,nst ,l,ierr
REAL occ(DIMENSION%nstd),a(DIMENSION%msh),b(DIMENSION%msh),ain(DIMENSION%msh),ahelp(DIMENSION%msh) ! ..
REAL occ_h(DIMENSION%nstd,2) ! .. Local Arrays ..
INTEGER kappa(DIMENSION%nstd),nprnc(DIMENSION%nstd)
CHARACTER(LEN=20) :: attributes(6) REAL rhcs(DIMENSION%msh),rhoc(DIMENSION%msh),rhoss(DIMENSION%msh),vrd(DIMENSION%msh),f(0:3)
REAL stateEnergies(29) REAL occ(DIMENSION%nstd),a(DIMENSION%msh),b(DIMENSION%msh),ain(DIMENSION%msh),ahelp(DIMENSION%msh)
! .. REAL occ_h(DIMENSION%nstd,2)
INTEGER kappa(DIMENSION%nstd),nprnc(DIMENSION%nstd)
c = c_light(1.0) CHARACTER(LEN=20) :: attributes(6)
seig = 0. REAL stateEnergies(29)
! ! ..
IF (input%frcor) THEN c = c_light(1.0)
DO n = 1,atoms%ntype seig = 0.
rnot = atoms%rmsh(1,n) ; dxx = atoms%dx(n) !
ncmsh = NINT( LOG( (atoms%rmt(n)+10.0)/rnot ) / dxx + 1 ) IF (input%frcor) THEN
ncmsh = MIN( ncmsh, DIMENSION%msh ) DO n = 1,atoms%ntype
! ---> update spherical charge density rnot = atoms%rmsh(1,n) ; dxx = atoms%dx(n)
DO i = 1,atoms%jri(n) ncmsh = NINT( LOG( (atoms%rmt(n)+10.0)/rnot ) / dxx + 1 )
rhoc(i) = rhc(i,n,jspin) ncmsh = MIN( ncmsh, DIMENSION%msh )
rho(i,0,n,jspin) = rho(i,0,n,jspin) + rhoc(i)/sfp_const ! ---> update spherical charge density
ENDDO DO i = 1,atoms%jri(n)
! ---> for total energy calculations, determine the sum of the rhoc(i) = rhc(i,n,jspin)
! ---> eigenvalues by requiring that the core kinetic energy rho(i,0,n,jspin) = rho(i,0,n,jspin) + rhoc(i)/sfp_const
! ---> remains constant. ENDDO
DO i = 1,atoms%jri(n) ! ---> for total energy calculations, determine the sum of the
rhoc(i) = rhoc(i)*vr(i,n)/atoms%rmsh(i,n) ! ---> eigenvalues by requiring that the core kinetic energy
ENDDO ! ---> remains constant.
nm = atoms%jri(n) DO i = 1,atoms%jri(n)
CALL intgr3(rhoc,atoms%rmsh(1,n),atoms%dx(n),nm,rhos) rhoc(i) = rhoc(i)*vr(i,n)/atoms%rmsh(i,n)
sea = tec(n,jspin) + rhos ENDDO
WRITE (16,FMT=8030) n,jspin,tec(n,jspin),sea nm = atoms%jri(n)
WRITE (6,FMT=8030) n,jspin,tec(n,jspin),sea CALL intgr3(rhoc,atoms%rmsh(1,n),atoms%dx(n),nm,rhos)
seig = seig + atoms%neq(n)*sea sea = tec(n,jspin) + rhos
ENDDO WRITE (16,FMT=8030) n,jspin,tec(n,jspin),sea
RETURN WRITE (6,FMT=8030) n,jspin,tec(n,jspin),sea
END IF seig = seig + atoms%neq(n)*sea
ENDDO
! ---> set up densities RETURN
DO jatom = 1,atoms%ntype END IF
sume = 0.
z = atoms%zatom(jatom) ! ---> set up densities
! rn = rmt(jatom) DO jatom = 1,atoms%ntype
dxx = atoms%dx(jatom) sume = 0.
bmu = 0.0 z = atoms%zatom(jatom)
CALL setcor(jatom,input%jspins,atoms,input,bmu,nst,kappa,nprnc,occ_h) ! rn = rmt(jatom)
IF ((bmu > 99.)) THEN dxx = atoms%dx(jatom)
occ(1:nst) = input%jspins * occ_h(1:nst,jspin) bmu = 0.0
ELSE CALL setcor(jatom,input%jspins,atoms,input,bmu,nst,kappa,nprnc,occ_h)
occ(1:nst) = occ_h(1:nst,1) IF ((bmu > 99.)) THEN
ENDIF occ(1:nst) = input%jspins * occ_h(1:nst,jspin)
rnot = atoms%rmsh(1,jatom) ELSE
d = EXP(atoms%dx(jatom)) occ(1:nst) = occ_h(1:nst,1)
ncmsh = NINT( LOG( (atoms%rmt(jatom)+10.0)/rnot ) / dxx + 1 ) ENDIF
ncmsh = MIN( ncmsh, DIMENSION%msh ) rnot = atoms%rmsh(1,jatom)
rn = rnot* (d** (ncmsh-1)) d = EXP(atoms%dx(jatom))
WRITE (6,FMT=8000) z,rnot,dxx,atoms%jri(jatom) ncmsh = NINT( LOG( (atoms%rmt(jatom)+10.0)/rnot ) / dxx + 1 )
WRITE (16,FMT=8000) z,rnot,dxx,atoms%jri(jatom) ncmsh = MIN( ncmsh, DIMENSION%msh )
DO j = 1,atoms%jri(jatom) rn = rnot* (d** (ncmsh-1))
rhoss(j) = 0.0 WRITE (6,FMT=8000) z,rnot,dxx,atoms%jri(jatom)
vrd(j) = vr(j,jatom) WRITE (16,FMT=8000) z,rnot,dxx,atoms%jri(jatom)
ENDDO DO j = 1,atoms%jri(jatom)
! rhoss(j) = 0.
IF (input%l_core_confpot) THEN vrd(j) = vr(j,jatom)
!---> linear extension of the potential with slope t1 / a.u. ENDDO
t1=0.125 !
t1 = MAX( (vrd(atoms%jri(jatom)) - vrd(atoms%jri(jatom)-1)*d)*& IF (input%l_core_confpot) THEN
d / (atoms%rmt(jatom)**2 * (d-1) ) , t1) !---> linear extension of the potential with slope t1 / a.u.
t2=vrd(atoms%jri(jatom))/atoms%rmt(jatom)-atoms%rmt(jatom)*t1 t1=0.125
rr = atoms%rmt(jatom) t1 = MAX( (vrd(atoms%jri(jatom)) - vrd(atoms%jri(jatom)-1)*d)*&
ELSE d / (atoms%rmt(jatom)**2 * (d-1) ) , t1)
t2 = vrd(atoms%jri(jatom)) / ( atoms%jri(jatom) - ncmsh ) t2=vrd(atoms%jri(jatom))/atoms%rmt(jatom)-atoms%rmt(jatom)*t1
ENDIF rr = atoms%rmt(jatom)
IF ( atoms%jri(jatom) < ncmsh) THEN ELSE
DO i = atoms%jri(jatom) + 1,ncmsh t2 = vrd(atoms%jri(jatom)) / ( atoms%jri(jatom) - ncmsh )
rhoss(i) = 0. ENDIF
IF (input%l_core_confpot) THEN IF ( atoms%jri(jatom) < ncmsh) THEN
rr = d*rr DO i = atoms%jri(jatom) + 1,ncmsh
vrd(i) = rr*( t2 + rr*t1 ) rhoss(i) = 0.
! vrd(i) = 2*vrd(jri(jatom)) - rr*( t2 + rr*t1 ) IF (input%l_core_confpot) THEN
ELSE rr = d*rr
vrd(i) = vrd(atoms%jri(jatom)) + t2* (i-atoms%jri(jatom)) vrd(i) = rr*( t2 + rr*t1 )
ENDIF ! vrd(i) = 2*vrd(jri(jatom)) - rr*( t2 + rr*t1 )
! ELSE
ENDDO vrd(i) = vrd(atoms%jri(jatom)) + t2* (i-atoms%jri(jatom))
END IF ENDIF
!
nst = atoms%ncst(jatom) ! for lda+U ENDDO
END IF
IF (input%gw==1 .OR. input%gw==3)&
& WRITE(15) nst,atoms%rmsh(1:atoms%jri(jatom),jatom) nst = atoms%ncst(jatom) ! for lda+U
IF (input%gw==1 .OR. input%gw==3)&
& WRITE(15) nst,atoms%rmsh(1:atoms%jri(jatom),jatom)
stateEnergies = 0.0
DO korb = 1,nst
IF (occ(korb) /= 0.0) THEN
fn = nprnc(korb)
fj = iabs(kappa(korb)) - .5e0
weight = 2*fj + 1.e0
IF (bmu > 99.) weight = occ(korb)
fl = fj + (.5e0)*isign(1,kappa(korb))
e = -2* (z/ (fn+fl))**2
CALL differ(fn,fl,fj,c,z,dxx,rnot,rn,d,ncmsh,vrd, e, a,b,ierr)
stateEnergies(korb) = e
WRITE (6,FMT=8010) fn,fl,fj,e,weight
WRITE (16,FMT=8010) fn,fl,fj,e,weight
IF (ierr/=0) CALL juDFT_error("error in core-level routine" ,calledby ="cored")
IF (input%gw==1 .OR. input%gw==3) WRITE (15) NINT(fl),weight,e,&
a(1:atoms%jri(jatom)),b(1:atoms%jri(jatom))
sume = sume + weight*e/input%jspins
DO j = 1,ncmsh
rhcs(j) = weight* (a(j)**2+b(j)**2)
rhoss(j) = rhoss(j) + rhcs(j)
ENDDO
ENDIF
ENDDO
stateEnergies = 0.0 ! ---->update spherical charge density rho with the core density.
DO korb = 1,nst ! ---->for spin-polarized (jspins=2), take only half the density
IF (occ(korb) /= 0.0) THEN nm = atoms%jri(jatom)
fn = nprnc(korb) DO j = 1,nm
rhoc(j) = rhoss(j)/input%jspins
rho(j,0,jatom,jspin) = rho(j,0,jatom,jspin) + rhoc(j)/sfp_const
ENDDO
rhc(1:ncmsh,jatom,jspin) = rhoss(1:ncmsh) / input%jspins
rhc(ncmsh+1:DIMENSION%msh,jatom,jspin) = 0.0
seig = seig + atoms%neq(jatom)*sume
DO i = 1,nm
rhoc(i) = rhoc(i)*vr(i,jatom)/atoms%rmsh(i,jatom)
ENDDO
CALL intgr3(rhoc,atoms%rmsh(1,jatom),atoms%dx(jatom),nm,rhs)
tec(jatom,jspin) = sume - rhs
WRITE (6,FMT=8030) jatom,jspin,tec(jatom,jspin),sume
WRITE (16,FMT=8030) jatom,jspin,tec(jatom,jspin),sume
! ---> simpson integration
rad = atoms%rmt(jatom)
q = rad*rhoss(nm)/2.
DO nm1 = nm + 1,ncmsh - 1,2
rad = d*rad
q = q + 2*rad*rhoss(nm1)
rad = d*rad
q = q + rad*rhoss(nm1+1)
ENDDO
q = 2*q*dxx/3
!+sb
WRITE (6,FMT=8020) q/input%jspins
WRITE (16,FMT=8020) q/input%jspins
!-sb
qint(jatom,jspin) = q*atoms%neq(jatom)
attributes = ''
WRITE(attributes(1),'(i0)') jatom
WRITE(attributes(2),'(i0)') NINT(z)
WRITE(attributes(3),'(i0)') jspin
WRITE(attributes(4),'(f18.10)') tec(jatom,jspin)
WRITE(attributes(5),'(f18.10)') sume
WRITE(attributes(6),'(f9.6)') q/input%jspins
CALL openXMLElementForm('coreStates',(/'atomType ','atomicNumber ','spin ','kinEnergy ',&
'eigValSum ','lostElectrons'/),&
attributes,RESHAPE((/8,12,4,9,9,13,6,3,1,18,18,9/),(/6,2/)))
DO korb = 1, atoms%ncst(jatom)
fj = iabs(kappa(korb)) - .5e0 fj = iabs(kappa(korb)) - .5e0
weight = 2*fj + 1.e0 weight = 2*fj + 1.e0
IF (bmu > 99.) weight = occ(korb) IF (bmu > 99.) weight = occ(korb)
fl = fj + (.5e0)*isign(1,kappa(korb)) fl = fj + (.5e0)*isign(1,kappa(korb))
eig = -2* (z/ (fn+fl))**2 attributes = ''
WRITE(attributes(1),'(i0)') nprnc(korb)
CALL differ(fn,fl,fj,c,z,dxx,rnot,rn,d,ncmsh,vrd, eig, a,b,ierr) WRITE(attributes(2),'(i0)') NINT(fl)
stateEnergies(korb) = eig WRITE(attributes(3),'(f4.1)') fj
WRITE (6,FMT=8010) fn,fl,fj,eig,weight WRITE(attributes(4),'(f20.10)') stateEnergies(korb)
WRITE (16,FMT=8010) fn,fl,fj,eig,weight WRITE(attributes(5),'(f15.10)') weight
IF (ierr/=0) CALL juDFT_error("error in core-level routine" ,calledby ="cored") CALL writeXMLElementForm('state',(/'n ','l ','j ','energy','weight'/),&
IF (input%gw==1 .OR. input%gw==3) WRITE (15) NINT(fl),weight,eig,& attributes(1:5),RESHAPE((/1,1,1,6,6,2,2,4,20,15/),(/5,2/)))
a(1:atoms%jri(jatom)),b(1:atoms%jri(jatom)) END DO
CALL closeXMLElement('coreStates')
sume = sume + weight*eig/input%jspins ENDDO
DO j = 1,ncmsh
rhcs(j) = weight* (a(j)**2+b(j)**2) RETURN
rhoss(j) = rhoss(j) + rhcs(j)
ENDDO 8000 FORMAT (/,/,10x,'z=',f4.0,5x,'r(1)=',e14.6,5x,'dx=',f9.6,5x,&
& 'm.t.index=',i4,/,15x,'n',4x,'l',5x,'j',4x,'energy',7x,&
IF(present(EnergyDen)) THEN & 'weight')
rhoss_aux = rhoss 8010 FORMAT (12x,2f5.0,f6.1,f10.4,f12.4)
DO j = 1,ncmsh 8020 FORMAT (f20.8,' electrons lost from core.')
! for energy density we want to multiply the weights 8030 FORMAT (10x,'atom type',i3,' (spin',i2,') ',/,10x,&
! with the eigenenergies & 'kinetic energy=',e20.12,5x,'sum of the eigenvalues=',&
rhoss_aux(j) = rhoss_aux(j) + (rhcs(j) * eig) & e20.12)
ENDDO END SUBROUTINE cored
ENDIF
ENDIF
ENDDO
! ---->update spherical charge density rho with the core density.
! ---->for spin-polarized (jspins=2), take only half the density
nm = atoms%jri(jatom)
DO j = 1,nm
rhoc(j) = rhoss(j)/input%jspins
rho(j,0,jatom,jspin) = rho(j,0,jatom,jspin) + rhoc(j)/sfp_const
ENDDO
IF(present(EnergyDen)) then
DO j = 1,nm
rhoc(j) = rhoss(j)/input%jspins
EnergyDen(j,0,jatom,jspin) = EnergyDen(j,0,jatom,jspin) &
+ rhoss_aux(j) /(input%jspins * sfp_const)
ENDDO
ENDIF
rhc(1:ncmsh,jatom,jspin) = rhoss(1:ncmsh) / input%jspins
rhc(ncmsh+1:DIMENSION%msh,jatom,jspin) = 0.0
seig = seig + atoms%neq(jatom)*sume
DO i = 1,nm
rhoc(i) = rhoc(i)*vr(i,jatom)/atoms%rmsh(i,jatom)
ENDDO
CALL intgr3(rhoc,atoms%rmsh(1,jatom),atoms%dx(jatom),nm,rhs)
tec(jatom,jspin) = sume - rhs
WRITE (6,FMT=8030) jatom,jspin,tec(jatom,jspin),sume
WRITE (16,FMT=8030) jatom,jspin,tec(jatom,jspin),sume
! ---> simpson integration
rad = atoms%rmt(jatom)
q = rad*rhoss(nm)/2.
DO nm1 = nm + 1,ncmsh - 1,2
rad = d*rad
q = q + 2*rad*rhoss(nm1)
rad = d*rad
q = q + rad*rhoss(nm1+1)
ENDDO
q = 2*q*dxx/3
!+sb
WRITE (6,FMT=8020) q/input%jspins
WRITE (16,FMT=8020) q/input%jspins
!-sb
qint(jatom,jspin) = q*atoms%neq(jatom)
attributes = ''
WRITE(attributes(1),'(i0)') jatom
WRITE(attributes(2),'(i0)') NINT(z)
WRITE(attributes(3),'(i0)') jspin
WRITE(attributes(4),'(f18.10)') tec(jatom,jspin)
WRITE(attributes(5),'(f18.10)') sume
WRITE(attributes(6),'(f9.6)') q/input%jspins
CALL openXMLElementForm('coreStates',(/'atomType ','atomicNumber ','spin ','kinEnergy ',&
'eigValSum ','lostElectrons'/),&
attributes,RESHAPE((/8,12,4,9,9,13,6,3,1,18,18,9/),(/6,2/)))
DO korb = 1, atoms%ncst(jatom)
fj = iabs(kappa(korb)) - .5e0
weight = 2*fj + 1.e0
IF (bmu > 99.) weight = occ(korb)
fl = fj + (.5e0)*isign(1,kappa(korb))
attributes = ''
WRITE(attributes(1),'(i0)') nprnc(korb)
WRITE(attributes(2),'(i0)') NINT(fl)
WRITE(attributes(3),'(f4.1)') fj
WRITE(attributes(4),'(f20.10)') stateEnergies(korb)
WRITE(attributes(5),'(f15.10)') weight
CALL writeXMLElementForm('state',(/'n ','l ','j ','energy','weight'/),&
attributes(1:5),RESHAPE((/1,1,1,6,6,2,2,4,20,15/),(/5,2/)))
END DO
CALL closeXMLElement('coreStates')
ENDDO
RETURN
8000 FORMAT (/,/,10x,'z=',f4.0,5x,'r(1)=',e14.6,5x,'dx=',f9.6,5x,&
& 'm.t.index=',i4,/,15x,'n',4x,'l',5x,'j',4x,'energy',7x,&
& 'weight')
8010 FORMAT (12x,2f5.0,f6.1,f10.4,f12.4)
8020 FORMAT (f20.8,' electrons lost from core.')
8030 FORMAT (10x,'atom type',i3,' (spin',i2,') ',/,10x,&
& 'kinetic energy=',e20.12,5x,'sum of the eigenvalues=',&
& e20.12)
END SUBROUTINE cored
END MODULE m_cored END MODULE m_cored
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment