Commit 911deba0 authored by Matthias Redies's avatar Matthias Redies

Merge branch 'develop' into MetaGGA

parents e4721869 3fbad4c3
......@@ -55,7 +55,7 @@ io/calculator.f global/ss_sym.f global/soc_sym.f math/inv3.f io/rw_symfile.f
kpoints/kptgen_hybrid.f kpoints/od_kptsgen.f kpoints/bravais.f kpoints/divi.f kpoints/brzone.f
kpoints/kptmop.f kpoints/kpttet.f init/bandstr1.F kpoints/ordstar.f kpoints/fulstar.f kpoints/kprep.f
kpoints/tetcon.f kpoints/kvecon.f init/boxdim.f global/radsra.f math/intgr.F global/differ.f math/inwint.f
math/outint.f xc-pot/gaunt.f math/grule.f
math/outint.f math/grule.f
)
set(inpgen_F90 ${inpgen_F90} global/constants.f90 io/xsf_io.f90
......@@ -65,7 +65,7 @@ inpgen/closure.f90 inpgen/inpgen_arguments.F90
juDFT/info.F90 juDFT/stop.F90 juDFT/args.F90 juDFT/time.F90 juDFT/init.F90 juDFT/sysinfo.F90 juDFT/string.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 kpoints/unfoldBandKPTS.f90)
eigen/tlmplm_store.F90 xc-pot/gaunt.f90 kpoints/unfoldBandKPTS.f90)
set(fleur_SRC ${fleur_F90} ${fleur_F77})
......
......@@ -4,59 +4,50 @@ xc-pot/metagga.F90
xc-pot/potl0.f90
xc-pot/mkgl0.f90
xc-pot/grdchlh.f90
xc-pot/easypbe.f90
xc-pot/corg91.f90
xc-pot/corl91.f90
xc-pot/corpbe.f90
xc-pot/excepbe.f90
xc-pot/exchpbe.f90
xc-pot/vxcepbe.f90
xc-pot/excl91.f90
xc-pot/excpw91.f90
xc-pot/gaunt.f90
xc-pot/excwb91.f90
xc-pot/pbecor2.f90
xc-pot/relcor.f90
xc-pot/vxcl91.f90
xc-pot/vxcpw91.f90
xc-pot/vxcwb91.f90
xc-pot/xcbh.f90
xc-pot/xcpz.f90
xc-pot/xcvwn.f90
xc-pot/xcwgn.f90
xc-pot/xcxal.f90
xc-pot/xch91.f90
)
set(fleur_F77 ${fleur_F77}
xc-pot/corg91.F
xc-pot/corl91.f
xc-pot/corpbe.f
xc-pot/easypbe.f
xc-pot/excepbe.f
xc-pot/exchpbe.F
xc-pot/excl91.f
xc-pot/excpw91.f
xc-pot/excwb91.f
xc-pot/gaunt.f
xc-pot/pbecor2.f
xc-pot/relcor.f
xc-pot/vxcepbe.f
xc-pot/vxcl91.f
xc-pot/vxcpw91.f
xc-pot/vxcwb91.f
xc-pot/xcbh.f
xc-pot/xch91.F
xc-pot/xcpz.f
xc-pot/xcvwn.f
xc-pot/xcwgn.f
xc-pot/xcxal.f
set(inpgen_F90 ${inpgen_F90}
xc-pot/easypbe.f90
xc-pot/corg91.f90
xc-pot/corl91.f90
xc-pot/corpbe.f90
xc-pot/excepbe.f90
xc-pot/exchpbe.f90
xc-pot/vxcepbe.f90
xc-pot/excl91.f90
xc-pot/excpw91.f90
xc-pot/excwb91.f90
xc-pot/pbecor2.f90
xc-pot/relcor.f90
xc-pot/vxcl91.f90
xc-pot/vxcpw91.f90
xc-pot/vxcwb91.f90
xc-pot/xcbh.f90
xc-pot/xcpz.f90
xc-pot/xcvwn.f90
xc-pot/xcwgn.f90
xc-pot/xcxal.f90
xc-pot/xch91.f90
)
set(ipngen_F90 ${inpgen_F90}
xc-pot/grdchlh.f90
)
set(inpgen_F77 ${inpgen_F77}
xc-pot/corg91.F
xc-pot/corl91.f
xc-pot/corpbe.f
xc-pot/easypbe.f
xc-pot/excepbe.f
xc-pot/exchpbe.F
xc-pot/excl91.f
xc-pot/excpw91.f
xc-pot/excwb91.f
xc-pot/pbecor2.f
xc-pot/relcor.f
xc-pot/vxcepbe.f
xc-pot/vxcl91.f
xc-pot/vxcpw91.f
xc-pot/vxcwb91.f
xc-pot/xcbh.f
xc-pot/xch91.F
xc-pot/xcpz.f
xc-pot/xcvwn.f
xc-pot/xcwgn.f
xc-pot/xcxal.f
)
MODULE m_corg91
c.....-----------------------------------------------------------------
c gga91 correlation
c.....-----------------------------------------------------------------
CONTAINS
SUBROUTINE corg91(fk,sk,gz,ec,ecrs,eczta,rs,zta,t,uu,vv,ww,h,
+ dvcup,dvcdn)
c.....-----------------------------------------------------------------
c input
c rs: seitz radius
c zta: relative spin polarization
c t: abs(grad d)/(d*2.*ks*gz)
c uu: (grad d)*grad(abs(grad d))/(d**2 * (2*ks*gz)**3)
c vv: (laplacian d)/(d * (2*ks*gz)**2)
c ww: (grad d)*(gradzta)/(d * (2*ks*gz)**2
c output
c h: nonlocal part of correlation energy per electron
c dvcup,-dn: nonlocal parts of correlation potentials.
MODULE m_corg91
!.....-----------------------------------------------------------------
! gga91 correlation
!.....-----------------------------------------------------------------
CONTAINS
SUBROUTINE corg91(fk,sk,gz,ec,ecrs,eczta,rs,zta,t,uu,vv,ww,h, &
dvcup,dvcdn)
!.....-----------------------------------------------------------------
! input
! rs: seitz radius
! zta: relative spin polarization
! t: abs(grad d)/(d*2.*ks*gz)
! uu: (grad d)*grad(abs(grad d))/(d**2 * (2*ks*gz)**3)
! vv: (laplacian d)/(d * (2*ks*gz)**2)
! ww: (grad d)*(gradzta)/(d * (2*ks*gz)**2
! output
! h: nonlocal part of correlation energy per electron
! dvcup,-dn: nonlocal parts of correlation potentials.
c with ks=sqrt(4*kf/pai), gz=[(1+zta)**(2/3)+(1-zta)**(2/3)]/2, &
c kf=cbrt(3*pai**2*d).
c.....-----------------------------------------------------------------
c.....-----------------------------------------------------------------
c .. previously untyped names ..
! with ks=sqrt(4*kf/pai), gz=[(1+zta)**(2/3)+(1-zta)**(2/3)]/2, &
! kf=cbrt(3*pai**2*d).
!.....-----------------------------------------------------------------
!.....-----------------------------------------------------------------
! .. previously untyped names ..
IMPLICIT NONE
REAL dvcdn,dvcup,ec,ecrs,eczta,fk,gz,h,rs,sk,t,uu,vv,ww,zta,a4,
+ alf,argmx,b,b2,b2fac,bec,bet,bg,c1,c2,c3,c4,c5,c6,cc,cc0,
+ ccrs,coeff,comm,cx,delt,fact0,fact1,fact2,fact3,fact4,fact5,
+ gm,gz3,gz4,h0,h0b,h0bt,h0rs,h0rst,h0t,h0tt,h0z,h0zt,h1,h1rs,
+ h1rst,h1t,h1tt,h1z,h1zt,hrs,hrst,ht,htt,hz,hzt,pon,pref,q4,
+ q5,q6,q7,q8,q9,r0,r1,r2,r3,r4,rs2,rs3,rsthrd,t2,t4,t6,thrd2,
+ thrdm,xnu,r1t2,expsm
c ..
REAL :: dvcdn,dvcup,ec,ecrs,eczta,fk,gz,h,rs,sk,t,uu,vv,ww,zta,a4, &
alf,argmx,b,b2,b2fac,bec,bet,bg,c1,c2,c3,c4,c5,c6,cc,cc0, &
ccrs,coeff,comm,cx,delt,fact0,fact1,fact2,fact3,fact4,fact5, &
gm,gz3,gz4,h0,h0b,h0bt,h0rs,h0rst,h0t,h0tt,h0z,h0zt,h1,h1rs, &
h1rst,h1t,h1tt,h1z,h1zt,hrs,hrst,ht,htt,hz,hzt,pon,pref,q4, &
q5,q6,q7,q8,q9,r0,r1,r2,r3,r4,rs2,rs3,rsthrd,t2,t4,t6,thrd2, &
thrdm,xnu,r1t2,expsm
! ..
DATA xnu,cc0,cx,alf/15.75592e0,0.00423500,-0.001667212e0,0.0900/
DATA c1,c2,c3,c4/0.00256800,0.02326600,7.389e-6,8.723e0/
DATA c5,c6,a4/0.472e0,7.389e-2,100.00/
DATA thrdm,thrd2/-0.333333333333e0,0.666666666667e0/
c.....-----------------------------------------------------------------
c expsm: argument of exponential-smallest.
!.....-----------------------------------------------------------------
! expsm: argument of exponential-smallest.
expsm=minexponent(expsm)/1.5
argmx = 174.0
bet = xnu*cc0
......@@ -45,10 +45,10 @@ c expsm: argument of exponential-smallest.
gz3 = gz**3
gz4 = gz3*gz
pon = -delt*ec/ (gz3*bet)
IF (pon.gt.argmx) THEN
b = 0.
IF (pon > argmx) THEN
b = 0.
ELSE
b = delt/ (exp(pon)-1.00)
b = delt/ (exp(pon)-1.00)
ENDIF
b2 = b*b
t2 = t*t
......@@ -65,26 +65,26 @@ c expsm: argument of exponential-smallest.
r1 = a4*r0*gz4
coeff = cc - cc0 - 3.e0*cx/7.00
r2 = xnu*coeff*gz3
c+tagu
! tagu
r1t2=max(-r1*t2,expsm)
r3 = exp(r1t2)
c-tagu
! tagu
h0 = gz3* (bet/delt)*log(1.00+delt*q4*t2/q5)
h1 = r3*r2*t2
h = h0 + h1
c local correlation option:
c h=0.0e0
! local correlation option:
! h=0.0e0
c energy done. now the potential:
! energy done. now the potential:
ccrs = (c2+2.*c3*rs)/q7 - q6* (c4+2.*c5*rs+3.*c6*rs2)/q7**2
rsthrd = rs/3.e0
r4 = rsthrd*ccrs/coeff
gm = ((1.00+zta)**thrdm- (1.00-zta)**thrdm)/3.00
IF (pon.gt.argmx) THEN
b2fac = 0.
IF (pon > argmx) THEN
b2fac = 0.
ELSE
b2fac = b2* (delt/b+1.00)
b2fac = b2* (delt/b+1.00)
ENDIF
bg = -3.e0*ec*b2fac/ (bet*gz4)
bec = b2fac/ (bet*gz3)
......@@ -122,9 +122,9 @@ c energy done. now the potential:
dvcup = comm + pref
dvcdn = comm - pref
c local correlation option:
c dvcup=0.0e0
c dvcdn=0.0e0
! local correlation option:
! dvcup=0.0e0
! dvcdn=0.0e0
END SUBROUTINE corg91
END MODULE m_corg91
END SUBROUTINE corg91
END MODULE m_corg91
MODULE m_corl91
c.....-----------------------------------------------------------------
c uniform-gas correlation of perdew and wang 1991
c.....-----------------------------------------------------------------
CONTAINS
SUBROUTINE corl91(
> rs,zta,
< ec,vcup,vcdn,ecrs,eczta,alfc)
c.....-----------------------------------------------------------------
c input: seitz radius (rs), relative spin polarization (zta)
c output: correlation energy per electron (ec),
c up- and down-spin potentials (vcup,vcdn),
c derivatives of ec wrt rs (ecrs) &zta (eczta).
c output: correlation contribution (alfc) to the spin stiffness
c.....-----------------------------------------------------------------
IMPLICIT NONE
REAL, INTENT (IN) :: rs,zta
REAL, INTENT (OUT) :: alfc,ec,ecrs,eczta,vcdn,vcup
REAL :: alfm,alfrsm,comm,ep,eprs,eu,eurs,f,fz,z4
REAL :: fzz,gam,thrd,thrd4
c.....-----------------------------------------------------------------
c ..
DATA gam,fzz/0.5198421,1.709921/
DATA thrd,thrd4/0.333333333333e0,1.333333333333e0/
c.....-----------------------------------------------------------------
f = ((1.0+zta)**thrd4+ (1.0-zta)**thrd4-2.e0)/gam
CALL gcor91(0.0310907,0.21370,7.5957,3.5876,1.6382,
+ 0.49294,1.00,rs,eu,eurs)
CALL gcor91(0.01554535,0.20548,14.1189,6.1977,3.3662,
+ 0.62517,1.00,rs,ep,eprs)
CALL gcor91(0.0168869,0.11125,10.357,3.6231,0.88026,
+ 0.49671,1.00,rs,alfm,alfrsm)
c alfm is minus the spin stiffness alfc
alfc = -alfm
z4 = zta**4
ec = eu* (1.0-f*z4) + ep*f*z4 - alfm*f* (1.0-z4)/fzz
c energy done. now the potential:
ecrs = eurs* (1.0-f*z4) + eprs*f*z4 - alfrsm*f* (1.0-z4)/fzz
fz = thrd4* ((1.0+zta)**thrd- (1.0-zta)**thrd)/gam
eczta = 4.e0* (zta**3)*f* (ep-eu+alfm/fzz) +
+ fz* (z4*ep-z4*eu- (1.0-z4)*alfm/fzz)
comm = ec - rs*ecrs/3.e0 - zta*eczta
vcup = comm + eczta
vcdn = comm - eczta
END SUBROUTINE corl91
c.....-----------------------------------------------------------------
c called by corl91
c.....-----------------------------------------------------------------
SUBROUTINE gcor91(
> a,a1,b1,b2,b3,b4,p,rs,
< gg,ggrs)
IMPLICIT NONE
REAL, INTENT (IN) :: a,a1,b1,b2,b3,b4,p,rs
REAL, INTENT (OUT) :: gg,ggrs
REAL :: p1,q0,q1,q2,q3,rs12,rs32,rsp
c.....-----------------------------------------------------------------
c ..
p1 = p + 1.0
q0 = -2.e0*a* (1.0+a1*rs)
rs12 = sqrt(rs)
rs32 = rs12**3
rsp = rs**p
q1 = 2.e0*a* (b1*rs12+b2*rs+b3*rs32+b4*rs*rsp)
q2 = log(1.0+1.0/q1)
gg = q0*q2
q3 = a* (b1/rs12+2.e0*b2+3.e0*b3*rs12+2.e0*b4*p1*rsp)
ggrs = -2.e0*a*a1*q2 - q0*q3/ (q1**2+q1)
c
END SUBROUTINE gcor91
END MODULE m_corl91
MODULE m_corl91
!.....-----------------------------------------------------------------
! uniform-gas correlation of perdew and wang 1991
!.....-----------------------------------------------------------------
CONTAINS
SUBROUTINE corl91( &
rs,zta, &
ec,vcup,vcdn,ecrs,eczta,alfc)
!.....-----------------------------------------------------------------
! input: seitz radius (rs), relative spin polarization (zta)
! output: correlation energy per electron (ec),
! up- and down-spin potentials (vcup,vcdn),
! derivatives of ec wrt rs (ecrs) &zta (eczta).
! output: correlation contribution (alfc) to the spin stiffness
!.....-----------------------------------------------------------------
IMPLICIT NONE
REAL, INTENT (IN) :: rs,zta
REAL, INTENT (OUT) :: alfc,ec,ecrs,eczta,vcdn,vcup
REAL :: alfm,alfrsm,comm,ep,eprs,eu,eurs,f,fz,z4
REAL :: fzz,gam,thrd,thrd4
!.....-----------------------------------------------------------------
! ..
DATA gam,fzz/0.5198421,1.709921/
DATA thrd,thrd4/0.333333333333e0,1.333333333333e0/
!.....-----------------------------------------------------------------
f = ((1.0+zta)**thrd4+ (1.0-zta)**thrd4-2.e0)/gam
CALL gcor91(0.0310907,0.21370,7.5957,3.5876,1.6382, &
& 0.49294,1.00,rs,eu,eurs)
CALL gcor91(0.01554535,0.20548,14.1189,6.1977,3.3662, &
& 0.62517,1.00,rs,ep,eprs)
CALL gcor91(0.0168869,0.11125,10.357,3.6231,0.88026, &
& 0.49671,1.00,rs,alfm,alfrsm)
! alfm is minus the spin stiffness alfc
alfc = -alfm
z4 = zta**4
ec = eu* (1.0-f*z4) + ep*f*z4 - alfm*f* (1.0-z4)/fzz
! energy done. now the potential:
ecrs = eurs* (1.0-f*z4) + eprs*f*z4 - alfrsm*f* (1.0-z4)/fzz
fz = thrd4* ((1.0+zta)**thrd- (1.0-zta)**thrd)/gam
eczta = 4.e0* (zta**3)*f* (ep-eu+alfm/fzz) + &
fz* (z4*ep-z4*eu- (1.0-z4)*alfm/fzz)
comm = ec - rs*ecrs/3.e0 - zta*eczta
vcup = comm + eczta
vcdn = comm - eczta
END SUBROUTINE corl91
!.....-----------------------------------------------------------------
! called by corl91
!.....-----------------------------------------------------------------
SUBROUTINE gcor91( &
a,a1,b1,b2,b3,b4,p,rs, &
gg,ggrs)
IMPLICIT NONE
REAL, INTENT (IN) :: a,a1,b1,b2,b3,b4,p,rs
REAL, INTENT (OUT) :: gg,ggrs
REAL :: p1,q0,q1,q2,q3,rs12,rs32,rsp
!.....-----------------------------------------------------------------
! ..
p1 = p + 1.0
q0 = -2.e0*a* (1.0+a1*rs)
rs12 = sqrt(rs)
rs32 = rs12**3
rsp = rs**p
q1 = 2.e0*a* (b1*rs12+b2*rs+b3*rs32+b4*rs*rsp)
q2 = log(1.0+1.0/q1)
gg = q0*q2
q3 = a* (b1/rs12+2.e0*b2+3.e0*b3*rs12+2.e0*b4*p1*rsp)
ggrs = -2.e0*a*a1*q2 - q0*q3/ (q1**2+q1)
END SUBROUTINE gcor91
END MODULE m_corl91
MODULE m_corpbe
c----------------------------------------------------------------------
c official pbe correlation code. k. burke, may 14, 1996.
c----------------------------------------------------------------------
c references:
c [a] j.p.~perdew, k.~burke, and m.~ernzerhof,
c {\sl generalized gradient approximation made simple}, sub.
c to phys. rev.lett. may 1996.
c [b] j. p. perdew, k. burke, and y. wang, {\sl real-space cutoff
c construction of a generalized gradient approximation: the pw91
c density functional}, submitted to phys. rev. b, feb. 1996.
c [c] j. p. perdew and y. wang, phys. rev. b {\bf 45}, 13244 (1992).
c----------------------------------------------------------------------
CONTAINS
SUBROUTINE corpbe(
> l_pbes,rs,zet,t,uu,vv,ww,lgga,lpot,
< ec,vcup,vcdn,h,dvcup,dvcdn)
c----------------------------------------------------------------------
c input: rs=seitz radius=(3/4pi rho)^(1/3)
c : zet=relative spin polarization = (rhoup-rhodn)/rho
c : t=abs(grad rho)/(rho*2.*ks*g) -- only needed for pbe
c : uu=(grad rho)*grad(abs(grad rho))/(rho**2 * (2*ks*g)**3)
c : vv=(laplacian rho)/(rho * (2*ks*g)**2)
c : ww=(grad rho)*(grad zet)/(rho * (2*ks*g)**2
c : uu,vv,ww, only needed for pbe potential
c : lgga=flag to do gga (0=>lsd only)
c : lpot=flag to do potential (0=>energy only)
c output: ec=lsd correlation energy from [a]
c : vcup=lsd up correlation potential
c : vcdn=lsd dn correlation potential
c : h=nonlocal part of correlation energy per electron
c : dvcup=nonlocal correction to vcup
c : dvcdn=nonlocal correction to vcdn
c----------------------------------------------------------------------
MODULE m_corpbe
!----------------------------------------------------------------------
! official pbe correlation code. k. burke, may 14, 1996.
!----------------------------------------------------------------------
! references:
! [a] j.p.~perdew, k.~burke, and m.~ernzerhof,
! {\sl generalized gradient approximation made simple}, sub.
! to phys. rev.lett. may 1996.
! [b] j. p. perdew, k. burke, and y. wang, {\sl real-space cutoff
! construction of a generalized gradient approximation: the pw91
! density functional}, submitted to phys. rev. b, feb. 1996.
! [c] j. p. perdew and y. wang, phys. rev. b {\bf 45}, 13244 (1992).
!----------------------------------------------------------------------
CONTAINS
SUBROUTINE corpbe( &
l_pbes,rs,zet,t,uu,vv,ww,lgga,lpot, &
ec,vcup,vcdn,h,dvcup,dvcdn)
!----------------------------------------------------------------------
! input: rs=seitz radius=(3/4pi rho)^(1/3)
! : zet=relative spin polarization = (rhoup-rhodn)/rho
! : t=abs(grad rho)/(rho*2.*ks*g) -- only needed for pbe
! : uu=(grad rho)*grad(abs(grad rho))/(rho**2 * (2*ks*g)**3)
! : vv=(laplacian rho)/(rho * (2*ks*g)**2)
! : ww=(grad rho)*(grad zet)/(rho * (2*ks*g)**2
! : uu,vv,ww, only needed for pbe potential
! : lgga=flag to do gga (0=>lsd only)
! : lpot=flag to do potential (0=>energy only)
! output: ec=lsd correlation energy from [a]
! : vcup=lsd up correlation potential
! : vcdn=lsd dn correlation potential
! : h=nonlocal part of correlation energy per electron
! : dvcup=nonlocal correction to vcup
! : dvcdn=nonlocal correction to vcdn
!----------------------------------------------------------------------
USE m_pbecor2
IMPLICIT NONE
......@@ -41,15 +41,15 @@ c----------------------------------------------------------------------
REAL, INTENT (IN) :: rs,zet,t,uu,vv,ww
REAL, INTENT (OUT) :: dvcdn,dvcup,ec,h,vcdn,vcup
c thrd*=various multiples of 1/3
c numbers for use in lsd energy spin-interpolation formula, [c](9).
c gam= 2^(4/3)-2
c fzz=f''(0)= 8/(9*gam)
c numbers for construction of pbe
c gamma=(1-log(2))/pi^2
c bet=coefficient in gradient expansion for correlation, [a](4).
c eta=small number to stop d phi/ dzeta from blowing up at
c |zeta|=1.
! thrd*=various multiples of 1/3
! numbers for use in lsd energy spin-interpolation formula, [c](9).
! gam= 2^(4/3)-2
! fzz=f''(0)= 8/(9*gam)
! numbers for construction of pbe
! gamma=(1-log(2))/pi^2
! bet=coefficient in gradient expansion for correlation, [a](4).
! eta=small number to stop d phi/ dzeta from blowing up at
! |zeta|=1.
REAL, PARAMETER :: thrd=1.e0/3.e0
REAL, PARAMETER :: thrdm=-thrd
......@@ -60,61 +60,61 @@ c |zeta|=1.
REAL, PARAMETER :: fzz=8.e0/ (9.e0*gam)
REAL, PARAMETER :: gamma=0.03109069086965489503494086371273e0
REAL, PARAMETER :: eta=1.e-12
c----------------------------------------------------------------------
c----------------------------------------------------------------------
c find lsd energy contributions, using [c](10) and table i[c].
c eu=unpolarized lsd correlation energy
c eurs=deu/drs
c ep=fully polarized lsd correlation energy
c eprs=dep/drs
c alfm=-spin stiffness, [c](3).
c alfrsm=-dalpha/drs
c f=spin-scaling factor from [c](9).
c construct ec, using [c](8)
!----------------------------------------------------------------------
!----------------------------------------------------------------------
! find lsd energy contributions, using [c](10) and table i[c].
! eu=unpolarized lsd correlation energy
! eurs=deu/drs
! ep=fully polarized lsd correlation energy
! eprs=dep/drs
! alfm=-spin stiffness, [c](3).
! alfrsm=-dalpha/drs
! f=spin-scaling factor from [c](9).
! construct ec, using [c](8)
REAL alfm,alfrsm,b,b2,bec,bg,comm,ecrs,eczet,ep,eprs,eu,eurs,
+ f,fac,fact0,fact1,fact2,fact3,fact5,fz,g,g3,g4,gz,hb,hbt,hrs,
+ hrst,ht,htt,hz,hzt,pon,pref,q4,q5,q8,q9,rsthrd,rtrs,
+ t2,t4,t6,z4,delt,bet
c ..
REAL :: alfm,alfrsm,b,b2,bec,bg,comm,ecrs,eczet,ep,eprs,eu,eurs, &
f,fac,fact0,fact1,fact2,fact3,fact5,fz,g,g3,g4,gz,hb,hbt,hrs, &
hrst,ht,htt,hz,hzt,pon,pref,q4,q5,q8,q9,rsthrd,rtrs, &
t2,t4,t6,z4,delt,bet
! ..
IF (l_pbes) THEN ! PBE_sol
bet=0.046e0
bet=0.046e0
ELSE
bet=0.06672455060314922e0
bet=0.06672455060314922e0
ENDIF
delt=bet/gamma
rtrs = sqrt(rs)
CALL pbecor2(0.0310907,0.21370,7.5957,3.5876,1.6382,
+ 0.49294,rtrs,eu,eurs)
CALL pbecor2(0.01554535,0.20548,14.1189,6.1977,3.3662,
+ 0.62517,rtrs,ep,eprs)
CALL pbecor2(0.0168869,0.11125,10.357,3.6231,0.88026,
+ 0.49671,rtrs,alfm,alfrsm)
CALL pbecor2(0.0310907,0.21370,7.5957,3.5876,1.6382, &
& 0.49294,rtrs,eu,eurs)
CALL pbecor2(0.01554535,0.20548,14.1189,6.1977,3.3662, &
& 0.62517,rtrs,ep,eprs)
CALL pbecor2(0.0168869,0.11125,10.357,3.6231,0.88026, &
& 0.49671,rtrs,alfm,alfrsm)
z4 = zet**4
f = ((1.e0+zet)**thrd4+ (1.e0-zet)**thrd4-2.e0)/gam
ec = eu* (1.e0-f*z4) + ep*f*z4 - alfm*f* (1.e0-z4)/fzz
c----------------------------------------------------------------------
c----------------------------------------------------------------------
c lsd potential from [c](a1)
c ecrs = dec/drs [c](a2)
c eczet=dec/dzeta [c](a3)
c fz = df/dzeta [c](a4)
!----------------------------------------------------------------------
!----------------------------------------------------------------------
! lsd potential from [c](a1)
! ecrs = dec/drs [c](a2)
! eczet=dec/dzeta [c](a3)
! fz = df/dzeta [c](a4)
ecrs = eurs* (1.e0-f*z4) + eprs*f*z4 - alfrsm*f* (1.e0-z4)/fzz
fz = thrd4* ((1.e0+zet)**thrd- (1.e0-zet)**thrd)/gam
eczet = 4.e0* (zet**3)*f* (ep-eu+alfm/fzz) +
+ fz* (z4*ep-z4*eu- (1.e0-z4)*alfm/fzz)
eczet = 4.e0* (zet**3)*f* (ep-eu+alfm/fzz) + &
fz* (z4*ep-z4*eu- (1.e0-z4)*alfm/fzz)
comm = ec - rs*ecrs/3.e0 - zet*eczet
vcup = comm + eczet
vcdn = comm - eczet
IF (lgga.EQ.0) RETURN
c----------------------------------------------------------------------
c----------------------------------------------------------------------
c pbe correlation energy
c g=phi(zeta), given after [a](3)
c delt=bet/gamma
c b=a of [a](8)
IF (lgga == 0) RETURN
!----------------------------------------------------------------------
!----------------------------------------------------------------------
! pbe correlation energy
! g=phi(zeta), given after [a](3)
! delt=bet/gamma
! b=a of [a](8)
g = ((1.e0+zet)**thrd2+ (1.e0-zet)**thrd2)/2.e0
g3 = g**3
pon = -ec/ (g3*gamma)
......@@ -125,15 +125,15 @@ c b=a of [a](8)
q4 = 1.e0 + b*t2
q5 = 1.e0 + b*t2 + b2*t4
h = g3* (bet/delt)*log(1.e0+delt*q4*t2/q5)
IF (lpot.EQ.0) RETURN
c----------------------------------------------------------------------
c----------------------------------------------------------------------
c energy done. now the potential, using appendix e of [b].
IF (lpot == 0) RETURN
!----------------------------------------------------------------------
!----------------------------------------------------------------------
! energy done. now the potential, using appendix e of [b].
g4 = g3*g
t6 = t4*t2
rsthrd = rs/3.e0
gz = (((1.e0+zet)**2+eta)**sixthm- ((1.e0-zet)**2+eta)**sixthm)/
+ 3.e0
gz = (((1.e0+zet)**2+eta)**sixthm- ((1.e0-zet)**2+eta)**sixthm)/ &
& 3.e0
fac = delt/b + 1.e0
bg = -3.e0*b2*ec*fac/ (bet*g4)
bec = b2*fac/ (bet*g3)
......@@ -157,6 +157,6 @@ c energy done. now the potential, using appendix e of [b].
comm = comm - pref*zet - uu*htt - vv*ht - ww* (hzt-fact5)
dvcup = comm + pref
dvcdn = comm - pref
END SUBROUTINE corpbe
END MODULE m_corpbe
END SUBROUTINE corpbe
END MODULE m_corpbe
MODULE m_easypbe
c----------------------------------------------------------------------
c easypbe---exchpbe
c --corpbe --- pbecor2
c----------------------------------------------------------------------
c easypbe is a driver for the pbe subroutines, using simple inputs
c k. burke, may 14, 1996.
c----------------------------------------------------------------------
CONTAINS
SUBROUTINE easypbe(xcpot,
> up,