1058 lines
34 KiB
Fortran
1058 lines
34 KiB
Fortran
subroutine ex_lda(rho_a,rho_b,ex,vx_a,vx_b)
|
|
include 'constants.include.F'
|
|
implicit none
|
|
double precision, intent(in) :: rho_a,rho_b
|
|
double precision, intent(out) :: ex,vx_a,vx_b
|
|
double precision :: tmp_a,tmp_b
|
|
tmp_a = rho_a**(c_1_3)
|
|
tmp_b = rho_b**(c_1_3)
|
|
ex = cst_lda * (tmp_a*tmp_a*tmp_a*tmp_a + tmp_b*tmp_b*tmp_b*tmp_b)
|
|
vx_a = cst_lda * c_4_3 * tmp_a
|
|
vx_b = cst_lda * c_4_3 * tmp_b
|
|
|
|
end
|
|
|
|
|
|
subroutine ec_lda(rho_a,rho_b,ec,vc_a,vc_b)
|
|
implicit none
|
|
include 'constants.include.F'
|
|
double precision, intent(out) :: ec
|
|
double precision, intent(out) :: vc_a,vc_b
|
|
double precision, intent(in) :: rho_a,rho_b
|
|
|
|
! Double precision numbers
|
|
|
|
double precision :: rsfac,rho,rs,rhoa,rhob,z
|
|
double precision :: eccoul, ecd, ecz, ecdd, eczd
|
|
double precision :: vcup,vcdown
|
|
rsfac = (3.0d0/(4.0d0*pi))**c_1_3
|
|
|
|
! Test on density
|
|
rho = rho_a + rho_b
|
|
if (dabs(rho).ge.1.d-10) then
|
|
|
|
rs=rsfac/(rho**c_1_3)
|
|
rhoa=max(rho_a,1.0d-15)
|
|
rhob=max(rho_b,1.0d-15)
|
|
z=(rhoa-rhob)/(rhoa+rhob)
|
|
|
|
call ecPW(rs,z,eccoul,ecd,ecz,ecdd,eczd)
|
|
ec=(eccoul)*rho
|
|
|
|
|
|
vcup=eccoul-rs/3.d0*ecd-(z-1.d0)*ecz
|
|
vcdown=eccoul-rs/3.d0*ecd-(z+1.d0)*ecz
|
|
vc_a = vcup
|
|
vc_b = vcdown
|
|
else
|
|
ec = 1.d-15
|
|
vc_a = 1.d-15
|
|
vc_b = 1.d-15
|
|
|
|
endif
|
|
|
|
end
|
|
|
|
subroutine ec_lda_sr(mu,rho_a,rho_b,ec,vc_a,vc_b)
|
|
implicit none
|
|
include 'constants.include.F'
|
|
double precision, intent(out) :: ec
|
|
double precision, intent(out) :: vc_a,vc_b
|
|
double precision, intent(in) :: mu,rho_a,rho_b
|
|
|
|
! Double precision numbers
|
|
|
|
double precision :: rsfac,rho,rs,rhoa,rhob,z
|
|
double precision :: eccoul, ecd, ecz, ecdd, eczd
|
|
double precision :: eclr,vcup,vcdown,vclrup,vclrdown,vclrupd,vclrdownd
|
|
rsfac = (3.0d0/(4.0d0*pi))**c_1_3
|
|
|
|
ec = 0.d0
|
|
vc_a = 0.d0
|
|
vc_b = 0.d0
|
|
! Test on density
|
|
rho = rho_a + rho_b
|
|
if (dabs(rho).ge.1.d-12) then
|
|
|
|
rs=rsfac/(rho**c_1_3)
|
|
rhoa=max(rho_a,1.0d-15)
|
|
rhob=max(rho_b,1.0d-15)
|
|
z=(rhoa-rhob)/(rhoa+rhob)
|
|
|
|
call ecPW(rs,z,eccoul,ecd,ecz,ecdd,eczd)
|
|
call ecorrlr(rs,z,mu,eclr)
|
|
ec=(eccoul-eclr)*rho
|
|
|
|
|
|
vcup=eccoul-rs/3.d0*ecd-(z-1.d0)*ecz
|
|
vcdown=eccoul-rs/3.d0*ecd-(z+1.d0)*ecz
|
|
call vcorrlr(rs,z,mu,vclrup,vclrdown,vclrupd,vclrdownd)
|
|
vc_a = vcup-vclrup
|
|
vc_b = vcdown-vclrdown
|
|
|
|
else
|
|
ec = 1.d-15
|
|
vc_a = 1.d-15
|
|
vc_b = 1.d-15
|
|
|
|
endif
|
|
|
|
end
|
|
|
|
subroutine ex_lda_sr(mu,rho_a,rho_b,ex,vx_a,vx_b)
|
|
include 'constants.include.F'
|
|
implicit none
|
|
double precision, intent(out) :: ex
|
|
double precision, intent(out) :: vx_a,vx_b
|
|
double precision, intent(in) :: rho_a,rho_b,mu
|
|
|
|
|
|
double precision :: rho_a_2,rho_b_2
|
|
double precision :: z0,z1,z2,z3,z4,z6,z8,z16,z24,z96,z12
|
|
double precision :: ex_a,ex_b
|
|
|
|
double precision :: f12,f13,f14,f32,f23,f43,f16
|
|
double precision :: ckf
|
|
double precision :: a, akf,a2, a3
|
|
double precision :: exp_f14a2
|
|
|
|
z0 = 0.D0
|
|
z1 = 1.D0
|
|
z2 = 2.D0
|
|
z3 = 3.D0
|
|
z4 = 4.D0
|
|
z6 = 6.D0
|
|
z8 = 8.D0
|
|
z12 = 12.D0
|
|
z16 = 16.D0
|
|
z24 = 24.D0
|
|
z96 = 96.D0
|
|
f12 = 0.5d0
|
|
f13 = 0.3333333333333333d0
|
|
f14 = 0.25d0
|
|
f32 = 1.5d0
|
|
f23 = 0.6666666666666666d0
|
|
f43 = 1.3333333333333333d0
|
|
f16 = 0.16666666666666666d0
|
|
ckf = 3.0936677262801355d0
|
|
|
|
!Density and kF
|
|
rho_a_2=rho_a*2.D0
|
|
akf = ckf*(rho_a_2**f13)
|
|
! Avoid division by zero
|
|
if (akf == 0.d0) akf = 1.d-20
|
|
a = mu/(z2*akf)
|
|
a2 = a*a
|
|
a3 = a2*a
|
|
|
|
!Test on the value of a
|
|
|
|
!Limit for small a (expansion not so important as for large a)
|
|
if (a.lt.1.d-9) then
|
|
ex_a = -z3/z8*rho_a_2*(z24*rho_a_2/pi)**f13
|
|
vx_a = - ((z3/pi)*rho_a_2)**f13
|
|
|
|
!Intermediate values of a
|
|
elseif (a.le.100d0) then
|
|
if(dabs(f14/a2).lt.50.d0)then
|
|
exp_f14a2 = dexp(-f14/a2)
|
|
else
|
|
exp_f14a2 = 0.d0
|
|
endif
|
|
ex_a = - (rho_a_2*(z24*rho_a_2/pi)**f13) * (z3/z8-a*(sqpi*derf(f12/a)+(z2*a-z4*a3)* exp_f14a2 -z3*a+z4*a3))
|
|
vx_a = -(z3*rho_a_2/pi)**f13 + z2*a*mu/pi*(exp_f14a2 - z1)+mu/sqpi * derf(f12/a)
|
|
|
|
|
|
!Expansion for large a
|
|
elseif (a.lt.1.d+9) then
|
|
ex_a = -(rho_a_2*(z24*rho_a_2/pi)**f13) * z1/(z96*a2)
|
|
vx_a = -pi*rho_a_2/(z2*mu*mu)
|
|
|
|
!Limit for large a
|
|
else
|
|
ex_a = 0.d0
|
|
vx_a = 0.d0
|
|
end if
|
|
|
|
!Density and kF
|
|
rho_b_2= rho_b * 2.d0
|
|
akf = ckf*(rho_b_2**f13)
|
|
if (akf == 0.d0) akf = 1.d-20
|
|
a = mu/(z2*akf)
|
|
a2 = a*a
|
|
a3 = a2*a
|
|
|
|
!Test on the value of a
|
|
|
|
!Limit for small a (expansion not so important as for large a)
|
|
if (a.lt.1.d-9) then
|
|
ex_b = -z3/z8*rho_b_2*(z24*rho_b_2/pi)**f13
|
|
vx_b = - ((z3/pi)*rho_b_2)**f13
|
|
|
|
!Intermediate values of a
|
|
elseif (a.le.100d0) then
|
|
if(dabs(f14/a2).lt.50.d0)then
|
|
exp_f14a2 = dexp(-f14/a2)
|
|
else
|
|
exp_f14a2 = 0.d0
|
|
endif
|
|
ex_b = - (rho_b_2*(z24*rho_b_2/pi)**f13)*(z3/z8-a*(sqpi*derf(f12/a)+(z2*a-z4*a3)*exp_f14a2-z3*a+z4*a3))
|
|
vx_b = -(z3*rho_b_2/pi)**f13+ z2*a*mu/pi*(exp_f14a2-z1)+mu/sqpi* derf(f12/a)
|
|
|
|
!Expansion for large a
|
|
elseif (a.lt.1.d+9) then
|
|
ex_b = - (rho_b_2*(z24*rho_b_2/pi)**f13) *z1/(z96*a2)
|
|
vx_b = - pi*rho_b_2/(z2*mu*mu)
|
|
|
|
!Limit for large a
|
|
else
|
|
ex_b = z0
|
|
vx_b = 0.d0
|
|
end if
|
|
|
|
ex = (ex_a+ex_b) * 0.5d0
|
|
end
|
|
|
|
|
|
subroutine ec_only_lda_sr(mu,rho_a,rho_b,ec)
|
|
implicit none
|
|
include 'constants.include.F'
|
|
double precision, intent(out) :: ec
|
|
double precision, intent(in) :: mu,rho_a,rho_b
|
|
|
|
! Double precision numbers
|
|
|
|
double precision :: rsfac,rho,rs,rhoa,rhob,z
|
|
double precision :: eccoul, ecd, ecz, ecdd, eczd
|
|
double precision :: eclr
|
|
rsfac = (3.0d0/(4.0d0*pi))**c_1_3
|
|
|
|
ec = 0.d0
|
|
! Test on density
|
|
rho = rho_a + rho_b
|
|
if (dabs(rho).ge.1.d-12) then
|
|
|
|
rs=rsfac/(rho**c_1_3)
|
|
rhoa=max(rho_a,1.0d-15)
|
|
rhob=max(rho_b,1.0d-15)
|
|
z=(rhoa-rhob)/(rhoa+rhob)
|
|
|
|
call ecPW(rs,z,eccoul,ecd,ecz,ecdd,eczd)
|
|
call ecorrlr(rs,z,mu,eclr)
|
|
ec=(eccoul-eclr)*rho
|
|
|
|
endif
|
|
|
|
end
|
|
|
|
!-------------------------------------------
|
|
function berf(a)
|
|
!-------------------------------------------
|
|
! Second-order exchange gradient expansion coefficient for erf
|
|
! interaction
|
|
! a = mu/(2*kF)
|
|
!
|
|
! Author : J. Toulouse
|
|
! Date : 10-03-04
|
|
!-------------------------------------------
|
|
implicit none
|
|
include 'constants.include.F'
|
|
|
|
double precision a
|
|
double precision eta,fak,berf,berf_dexp
|
|
|
|
! function
|
|
double precision derf
|
|
|
|
eta=19.0d0
|
|
if(dabs(eta*a*a).lt.50.d0)then
|
|
fak=2.540118935556d0*dexp(-eta*a*a)
|
|
else
|
|
fak=0.d0
|
|
endif
|
|
|
|
if(a .lt. 0.075d0) then
|
|
! expansion for small mu to avoid numerical problems
|
|
! denominator becomes zero for a approximately 0.4845801308
|
|
! (and for one negative and two complex values of a)
|
|
berf = (-7d0+72.d0*a*a)/(27.d0*(-3d0-24.d0*a*a+32.d0*a**4+8d0*dsqrt(pi)*a))
|
|
|
|
else if(a .gt. 50.d0) then
|
|
berf = 1.d0/(72.d0*a*a)-1.d0/(17280.d0*a**4)- 23.d0/(358400.d0*a**6)
|
|
|
|
else
|
|
|
|
|
|
! Code generated by Mathematica
|
|
berf_dexp=dexp(2.5d-1/a**2)
|
|
berf = (1.851851851851851851851852d-2*(-1.d0 + 1.44d2*a**4*(-1.d0 &
|
|
+ berf_dexp) - 2.d0*a**2*(1.1d1 + 7.d0*berf_dexp &
|
|
)))/(a**2*(3.2d1*a**4*(-1.d0 + berf_dexp) - 3.d0*berf_dexp &
|
|
+ 1.417963080724412821838534d1*a*derf(5.d-1/a)*berf_dexp &
|
|
- 8.d0*a**2*(-2.d0 + 3.d0*berf_dexp)))
|
|
|
|
end if
|
|
|
|
berf=berf*fak
|
|
|
|
return
|
|
end
|
|
|
|
!-------------------------------------------
|
|
function dberfda(a)
|
|
!-------------------------------------------
|
|
! Derivative of second-order exchange gradient
|
|
! expansion coefficient for erf interaction
|
|
! a = mu/(2*kF)
|
|
!
|
|
! Author : J. Toulouse
|
|
! Date : 10-03-04
|
|
!-------------------------------------------
|
|
implicit none
|
|
include 'constants.include.F'
|
|
|
|
double precision a
|
|
double precision eta,fak,dfakda,berf,dberfda,berf_dexp
|
|
double precision t1,t2,tdexp,t3,t4,t5
|
|
|
|
eta=19.0d0
|
|
if(dabs(eta*a*a).lt.50.d0)then
|
|
fak=2.540118935556d0*dexp(-eta*a*a)
|
|
else
|
|
fak=0.d0
|
|
endif
|
|
dfakda=-2.0d0*eta*a*fak
|
|
|
|
if(a .lt. 0.075d0) then
|
|
! expansion for small mu to avoid numerical problems
|
|
! denominator becomes zero for a approximately 0.4845801308
|
|
! (and for one negative and two complex values of a)
|
|
berf = (-7d0+72.d0*a*a)/(27.d0*(-3d0-24.d0*a*a+32.d0*a**4+8d0*dsqrt(pi)*a))
|
|
dberfda = (8d0*(-96.d0*a + 112.d0*a**3 - 576.d0*a**5 &
|
|
+ 7d0*dsqrt(pi) + 72.d0*a**2*dsqrt(pi)))/ &
|
|
(27.d0*(3d0 + 24.d0*a**2 - 32.d0*a**4 - 8d0*a*dsqrt(pi))**2)
|
|
|
|
else if(a .gt. 50.d0) then
|
|
berf = 1.d0/(72.d0*a*a)-1.d0/(17280.d0*a**4)- 23.d0/(358400.d0*a**6)
|
|
dberfda = - 1.d0/(36.d0*a**3) + 1.d0/(4320.d0*a**5)+ 69.d0/(179200.d0*a**7)
|
|
|
|
|
|
else
|
|
|
|
! Code generated by Mathematica
|
|
berf_dexp=dexp(2.5d-1/a**2)
|
|
|
|
berf = (1.851851851851851851851852d-2*(-1.d0 + 1.44d2*a**4*(-1.d0 + berf_dexp) - 2.d0*a**2*(1.1d1 + 7.d0*berf_dexp )))/(a**2*(3.2d1*a**4*(-1.d0 + berf_dexp) - 3.d0*berf_dexp + 1.417963080724412821838534d1*a*derf(5.d-1/a)*berf_dexp - 8.d0*a**2*(-2.d0 + 3.d0*berf_dexp)))
|
|
|
|
tdexp=dexp(2.5d-1/a**2)
|
|
t1 = (1.851851851851851851851852d-2*(5.76d2*a**3*(-1.d0 + tdexp ) + (7.d0*tdexp)/a - 7.2d1*a*tdexp - 4.d0*a*(1.1d1 + 7.d0*tdexp)))/(a**2*(3.2d1*a**4*(-1.d0 + tdexp) - 3.d0*tdexp + 1.417963080724412821838534d1*a*derf(5.d-1/a)*tdexp - 8.d0*a**2*(-2.d0 + 3.d0*tdexp)))
|
|
t2 = -1.851851851851851851851852d-2/a**2
|
|
t3 = -8.d0/a + 1.28d2*a**3*(-1.d0 + tdexp) + (1.5d0*tdexp)/a**3 + (1.2d1*tdexp)/a - 1.6d1*a* tdexp + 1.417963080724412821838534d1*derf(5.d-1/a)*tdexp - (7.08981540362206410919267d0*derf(5.d-1/a)*tdexp)/a**2 - 1.6d1*a*(-2.d0 + 3.d0*tdexp)
|
|
t4 = (-1.d0 + 1.44d2*a**4*(-1.d0 + tdexp) - 2.d0*a**2*(1.1d1 + 7.d0*tdexp))/(3.2d1*a**4*(-1.d0 + tdexp) - 3.d0*tdexp + 1.417963080724412821838534d1*a*derf(5.d-1/a)*tdexp - 8.d0*a**2*(-2.d0 + 3.d0*tdexp))**2
|
|
t5 = (-3.703703703703703703703704d-2*(-1.d0 + 1.44d2*a**4*(-1.d0 + tdexp) - 2.d0*a**2*(1.1d1 + 7.d0*tdexp )))/(a**3*(3.2d1*a**4*(-1.d0 + tdexp) - 3.d0*tdexp+ 1.417963080724412821838534d1*a*derf(5.d-1/a)*tdexp- 8.d0*a**2*(-2.d0 + 3.d0*tdexp)))
|
|
dberfda = t1 + t2*t3*t4 + t5
|
|
|
|
end if
|
|
|
|
dberfda=dberfda*fak+berf*dfakda
|
|
|
|
return
|
|
end
|
|
|
|
|
|
subroutine ecorrlr(rs,z,mu,eclr)
|
|
!cc Hartree atomic units used
|
|
!cc for given density parameter rs, spin polarization z
|
|
!cc and cutoff parameter mu
|
|
!cc gives the correlation energy of the LR gas
|
|
!cc => eclr
|
|
implicit none
|
|
double precision rs,z,mu,eclr,ec,ecd,ecz
|
|
double precision pi,alpha,cf,phi
|
|
double precision g0f,dpol,d2anti,d3anti,Qrpa
|
|
double precision coe2,coe3,coe4,coe5
|
|
double precision a1,a2,a3,a4,b0
|
|
double precision q1a,q2a,q3a,t1a,t2a,t3a,adib
|
|
!SCD
|
|
double precision ecdd,eczd
|
|
!SCF
|
|
pi=dacos(-1.d0)
|
|
alpha=(4.d0/9.d0/pi)**(1.d0/3.d0)
|
|
cf=1.d0/alpha
|
|
|
|
phi=((1.d0+z)**(2.d0/3.d0)+(1.d0-z)**(2.d0/3.d0))/2.d0
|
|
!c parameters from the fit
|
|
adib = 0.784949d0
|
|
q1a = -0.388d0
|
|
q2a = 0.676d0
|
|
q3a = 0.547d0
|
|
t1a = -4.95d0
|
|
t2a = 1.d0
|
|
t3a = 0.31d0
|
|
|
|
b0=adib*rs
|
|
|
|
double precision :: exp_q3a_rs
|
|
if(dabs(q3a*rs).lt.50.d0)then
|
|
exp_q3a_rs = dexp(-dabs(q3a)*rs)
|
|
else
|
|
exp_q3a_rs = 0.d0
|
|
endif
|
|
d2anti=(q1a*rs+q2a*rs**2)*exp_q3a_rs/rs**2
|
|
double precision :: exp_t3a_rs
|
|
if(dabs(t3a*rs).lt.50.d0)then
|
|
exp_t3a_rs = dexp(-dabs(t3a)*rs)
|
|
else
|
|
exp_t3a_rs = 0.d0
|
|
endif
|
|
d3anti=(t1a*rs+t2a*rs**2)*exp_t3a_rs/rs**3
|
|
|
|
coe2=-3.d0/8.d0/rs**3*(1.d0-z**2)*(g0f(rs)-0.5d0)
|
|
|
|
coe3=-(1.d0-z**2)*g0f(rs)/(dsqrt(2.d0*pi)*rs**3)
|
|
|
|
if(dabs(z).eq.1.d0) then
|
|
|
|
coe4=-9.d0/64.d0/rs**3*(dpol(rs) -cf**2*2d0**(5.d0/3.d0)/5.d0/rs**2)
|
|
coe5=-9.d0/40.d0/(dsqrt(2.d0*pi)*rs**3)*dpol(rs)
|
|
|
|
else
|
|
|
|
coe4=-9.d0/64.d0/rs**3*(((1.d0+z)/2.d0)**2* &
|
|
dpol(rs*(2d0/(1.d0+z))**(1.d0/3.d0))+((1.d0-z)/2.d0)**2 &
|
|
*dpol(rs*(2.d0/(1.d0-z))**(1.d0/3.d0))+ &
|
|
(1.-z**2)*d2anti-cf**2/10.d0*((1.d0+z)**(8.d0/3.d0) &
|
|
+(1.-z)**(8.d0/3.d0))/rs**2)
|
|
|
|
coe5=-9.d0/40.d0/(dsqrt(2.d0*pi)*rs**3)*(((1.d0+z)/2.d0)**2 &
|
|
*dpol(rs*(2.d0/(1.d0+z))**(1.d0/3.d0))+((1.d0-z)/2.d0)**2 &
|
|
*dpol(rs*(2.d0/(1.d0-z))**(1.d0/3.d0))+(1.d0-z**2)* &
|
|
d3anti)
|
|
end if
|
|
|
|
! call ecPW(rs,z,ec,ecd,ecz)
|
|
!SCD
|
|
call ecPW(rs,z,ec,ecd,ecz,ecdd,eczd)
|
|
!SCF
|
|
|
|
a1=4.d0*b0**6*coe3+b0**8*coe5
|
|
a2=4.d0*b0**6*coe2+b0**8*coe4+6.d0*b0**4*ec
|
|
a3=b0**8*coe3
|
|
a4=b0**6*(b0**2*coe2+4.d0*ec)
|
|
|
|
if(mu*dsqrt(rs)/phi.lt.0.d0)then
|
|
print*,'phi',phi
|
|
print*,'mu ',mu
|
|
print*,'rs ',rs
|
|
stop -1
|
|
endif
|
|
eclr=(phi**3*Qrpa(mu*dsqrt(rs)/phi)+a1*mu**3+a2*mu**4+a3*mu**5+ &
|
|
a4*mu**6+b0**8*mu**8*ec)/((1.d0+b0**2*mu**2)**4)
|
|
|
|
return
|
|
end
|
|
|
|
subroutine vcorrlr(rs,z,mu,vclrup,vclrdown,vclrupd,vclrdownd)
|
|
!SCF
|
|
!cc Hartree atomic units used
|
|
!cc for given density parameter rs, spin polarization z
|
|
!cc and cutoff mu it gives the correlation LSD potential for LR interaction
|
|
!cc => vclrup (spin-up electrons), vclrdown (spin-down electrons)
|
|
implicit none
|
|
double precision rs,z,mu,eclr,eclrrs,eclrz,vclrup,vclrdown
|
|
double precision ec,ecd,ecz
|
|
double precision pi,alpha,cf,phi
|
|
double precision g0f,dpol,d2anti,d3anti,Qrpa
|
|
double precision g0d,dpold,d2antid,d3antid,Qrpad,x
|
|
double precision coe2,coe3,coe4,coe5
|
|
double precision coe2rs,coe3rs,coe4rs,coe5rs
|
|
double precision coe2z,coe3z,coe4z,coe5z
|
|
double precision a1,a2,a3,a4,a5,b0,a1rs,a2rs,a3rs,a4rs,a5rs,b0rs,a1z,a2z,a3z,a4z,a5z,b0z
|
|
double precision q1a,q2a,q3a,t1a,t2a,t3a,adib
|
|
!SCD
|
|
double precision coe2rsd,coe3rsd,coe4rsd,coe5rsd,f23
|
|
double precision coe2zd,coe3zd,coe4zd,coe5zd
|
|
double precision g0dd,dpoldd,d2antidd,d3antidd
|
|
double precision a1rsd,a2rsd,a3rsd,a4rsd,a5rsd,a1zd,a2zd,a3zd,a4zd,a5zd
|
|
double precision ecdd,eczd,eclrrsd,vclrupd,vclrdownd
|
|
double precision u,du,ddu,v,dv,ddv,Qrpadd,eclrzd
|
|
!SCF
|
|
double precision sqrt2pi
|
|
pi=dacos(-1.d0)
|
|
alpha=(4.d0/9.d0/pi)**(1.d0/3.d0)
|
|
cf=1.d0/alpha
|
|
! sqrt2pi=sqrt(2.d0*pi)
|
|
sqrt2pi=2.5066282746310002d0
|
|
|
|
phi=((1.d0+z)**(2.d0/3.d0)+(1.d0-z)**(2.d0/3.d0))/2.d0
|
|
!c parameters from the fit
|
|
adib = 0.784949d0
|
|
q1a = -0.388d0
|
|
q2a = 0.676d0
|
|
q3a = 0.547d0
|
|
t1a = -4.95d0
|
|
t2a = 1.d0
|
|
t3a = 0.31d0
|
|
!SCD
|
|
f23 = 2.d0/3.d0
|
|
!SCF
|
|
|
|
b0=adib*rs
|
|
double precision :: exp_q3a_rs,exp_t3a_rs
|
|
if(dabs(q3a*rs).lt.50.d0)then
|
|
exp_q3a_rs = dexp(-q3a*rs)
|
|
else
|
|
exp_q3a_rs = 0.d0
|
|
endif
|
|
if(dabs(t3a*rs).lt.50.d0)then
|
|
exp_t3a_rs = dexp(-t3a*rs)
|
|
else
|
|
exp_t3a_rs = 0.d0
|
|
endif
|
|
|
|
d2anti=(q1a+q2a*rs)*exp_q3a_rs/rs
|
|
d3anti=(t1a+t2a*rs)*exp_t3a_rs/rs**2
|
|
|
|
d2antid=-((q1a + q1a*q3a*rs + q2a*q3a*rs**2)/rs**2)*exp_q3a_rs
|
|
d3antid=-((rs*t2a*(1d0 + rs*t3a) + t1a*(2d0 + rs*t3a))/rs**3)*exp_t3a_rs
|
|
|
|
!SCD
|
|
d2antidd = exp_q3a_rs/rs**3*( &
|
|
q3a**2*q1a*rs**2+q2a*q3a**2*rs**3 &
|
|
+2.d0*q3a*q1a*rs+2.d0*q1a)
|
|
d3antidd = exp_t3a_rs/rs**4* &
|
|
(2.d0*t3a*t2a*rs**2 + 2.d0*t2a*rs &
|
|
+ t1a*t3a**2*rs**2 + t2a*t3a**2*rs**3 &
|
|
+ 4.d0*t1a*t3a*rs + 6.d0*t1a)
|
|
!SCF
|
|
coe2=-3.d0/8.d0/rs**3*(1.d0-z**2)*(g0f(rs)-0.5d0)
|
|
coe2rs=-3.d0/8.d0/rs**3*(1.d0-z**2)*g0d(rs)+ &
|
|
9.d0/8.d0/rs**4*(1.d0-z**2)*(g0f(rs)-0.5d0)
|
|
coe2z=-3.d0/8.d0/rs**3*(-2.d0*z)*(g0f(rs)-0.5d0)
|
|
!SCD
|
|
coe2rsd=(1.d0-z**2)*(9.d0/4.d0/rs**4*g0d(rs) &
|
|
-3.d0/8.d0/rs**3*g0dd(rs) &
|
|
-9.d0/2.d0/rs**5*(g0f(rs)-0.5d0))
|
|
! coe2zd=3.d0/4.d0/rs**3*(g0f(rs)-0.5d0)
|
|
coe2zd=0.d0
|
|
!SCF
|
|
|
|
coe3=-(1.d0-z**2)*g0f(rs)/(sqrt2pi*rs**3)
|
|
coe3rs=-(1.d0-z**2)*g0d(rs)/(sqrt2pi*rs**3)+ &
|
|
3.d0*(1.d0-z**2)*g0f(rs)/(sqrt2pi*rs**4)
|
|
coe3z=2.d0*z*g0f(rs)/(sqrt2pi*rs**3)
|
|
!SCD
|
|
coe3rsd=(1.d0-z**2)/(sqrt2pi*rs**5) &
|
|
*(6.d0*rs*g0d(rs)-12.d0*g0f(rs) &
|
|
- g0dd(rs)*rs**2)
|
|
! coe3zd=2.d0*g0f(rs)/(sqrt2pi*rs**3)
|
|
coe3zd=0.d0
|
|
!SCF
|
|
|
|
if(abs(z).eq.1.d0) then
|
|
|
|
coe4=-9.d0/64.d0/rs**3*(dpol(rs) &
|
|
-cf**2*2d0**(5.d0/3.d0)/5.d0/rs**2)
|
|
coe4rs=-3.d0/rs*coe4-9.d0/64.d0/rs**3*(dpold(rs) &
|
|
+2.d0*cf**2*2d0**(5.d0/3.d0)/5.d0/rs**3)
|
|
coe4z=-9.d0/64.d0/rs**3*(dpol(rs)-rs/6.d0*dpold(rs)-2.d0*d2anti &
|
|
-4.d0/15.d0/rs**2*cf**2*2.d0**(5.d0/3.d0))*z
|
|
coe5=-9.d0/40.d0/(sqrt2pi*rs**3)*dpol(rs)
|
|
coe5rs=-3.d0/rs*coe5-9.d0/40.d0/(sqrt2pi*rs**3)*dpold(rs)
|
|
coe5z=-9.d0/40.d0/(sqrt2pi*rs**3)*(dpol(rs)-rs/6.d0* &
|
|
dpold(rs)-2.d0*d3anti)*z
|
|
!SCD
|
|
coe4rsd = -9.d0/64.d0/rs**7*(12.d0*dpol(rs)*rs**2 &
|
|
-12.d0*cf**2*2d0**(f23) &
|
|
-6.d0*dpold(rs)*rs**3 &
|
|
+dpoldd(rs)*rs**4)
|
|
coe4zd = 0.d0
|
|
|
|
coe5rsd = -9.d0/40.d0/dsqrt(2.d0/pi)/rs**5* &
|
|
(12.d0*dpol(rs)-6.d0*rs*dpold(rs) &
|
|
+rs**2*dpoldd(rs))
|
|
coe5zd = 0.d0
|
|
!SCF
|
|
|
|
else
|
|
|
|
coe4=-9.d0/64.d0/rs**3*(((1.d0+z)/2.d0)**2* &
|
|
dpol(rs*(2d0/(1.d0+z))**(1.d0/3.d0))+((1.d0-z)/2.d0)**2 &
|
|
*dpol(rs*(2.d0/(1.d0-z))**(1.d0/3.d0))+ &
|
|
(1.-z**2)*d2anti-cf**2/10.d0*((1.d0+z)**(8.d0/3.d0) &
|
|
+(1.-z)**(8.d0/3.d0))/rs**2)
|
|
coe4rs=-3.d0/rs*coe4-9.d0/64.d0/rs**3*( &
|
|
((1.d0+z)/2.d0)**(5.d0/3.d0)*dpold(rs*(2d0/(1.d0+z))** &
|
|
(1.d0/3.d0))+((1.d0-z)/2.d0)**(5.d0/3.d0)* &
|
|
dpold(rs*(2d0/(1.d0-z))**(1.d0/3.d0))+(1.d0-z**2)* &
|
|
d2antid+cf**2/5.d0*((1.d0+z)**(8.d0/3.d0) &
|
|
+(1.d0-z)**(8.d0/3.d0))/rs**3)
|
|
coe4z=-9.d0/64.d0/rs**3*(1.d0/2.d0*(1.d0+z)* &
|
|
dpol(rs*(2d0/(1.d0+z))**(1.d0/3.d0))-1.d0/2.d0*(1.d0-z)* &
|
|
dpol(rs*(2d0/(1.d0-z))**(1.d0/3.d0))-rs/6.d0* &
|
|
((1.d0+z)/2.d0)**(2.d0/3.d0)*dpold(rs*(2d0/(1.d0+z)) &
|
|
**(1.d0/3.d0))+rs/6.d0*((1.d0-z)/2.d0)**(2.d0/3.d0) &
|
|
*dpold(rs*(2d0/(1.d0-z))**(1.d0/3.d0))-2.d0*z*d2anti- &
|
|
4.d0/15.d0/rs**2*cf**2*((1.d0+z)**(5.d0/3.d0)- &
|
|
(1.d0-z)**(5.d0/3.d0)))
|
|
|
|
coe5=-9.d0/40.d0/(sqrt2pi*rs**3)*(((1.d0+z)/2.d0)**2 &
|
|
*dpol(rs*(2.d0/(1.d0+z))**(1.d0/3.d0))+((1.d0-z)/2.d0)**2 &
|
|
*dpol(rs*(2.d0/(1.d0-z))**(1.d0/3.d0))+(1.d0-z**2)* &
|
|
d3anti)
|
|
coe5rs=-3.d0/rs*coe5-9.d0/(40.d0*sqrt2pi*rs**3)*( &
|
|
((1.d0+z)/2.d0)**(5.d0/3.d0)*dpold(rs*(2d0/(1.d0+z))** &
|
|
(1.d0/3.d0))+((1.d0-z)/2.d0)**(5.d0/3.d0)* &
|
|
dpold(rs*(2d0/(1.d0-z))**(1.d0/3.d0))+(1.d0-z**2)* &
|
|
d3antid)
|
|
coe5z=-9.d0/40.d0/(sqrt2pi*rs**3)*(1.d0/2.d0*(1.d0+z)* &
|
|
dpol(rs*(2d0/(1.d0+z))**(1.d0/3.d0))-1.d0/2.d0*(1.d0-z)* &
|
|
dpol(rs*(2d0/(1.d0-z))**(1.d0/3.d0))-rs/6.d0* &
|
|
((1.d0+z)/2.d0)**(2.d0/3.d0)*dpold(rs*(2d0/(1.d0+z)) &
|
|
**(1.d0/3.d0))+rs/6.d0*((1.d0-z)/2.d0)**(2.d0/3.d0) &
|
|
*dpold(rs*(2d0/(1.d0-z))**(1.d0/3.d0))-2.d0*z*d3anti)
|
|
!SCD
|
|
! coe4rsd=+3.d0/rs**2*coe4-3.d0/rs*coe4rs+27.d0/64.d0/rs**4*(
|
|
! S ((1.d0+z)/2.d0)**(5.d0/3.d0)*dpold(rs*(2/(1.d0+z))**
|
|
! S (1.d0/3.d0))+((1.d0-z)/2.d0)**(5.d0/3.d0)*
|
|
! S dpold(rs*(2/(1.d0-z))**(1.d0/3.d0))+(1.d0-z**2)*
|
|
! S d2antid+cf**2/5.d0*((1.d0+z)**(8.d0/3.d0)
|
|
! S +(1.d0-z)**(8.d0/3.d0))/rs**3)-9.d0/64.d0/rs**3*(
|
|
! S ((1.d0+z)/2.d0)**(4.d0/3.d0)*dpoldd(rs*(2/(1.d0+z))**
|
|
! S (1.d0/3.d0))+((1.d0-z)/2.d0)**(4.d0/3.d0)*
|
|
! S dpoldd(rs*(2/(1.d0-z))**(1.d0/3.d0))+(1.d0-z**2)*
|
|
! S d2antidd-3.d0*cf**2/5.d0*((1.d0+z)**(8.d0/3.d0)
|
|
! S +(1.d0-z)**(8.d0/3.d0))/rs**4)
|
|
! Case where z=0
|
|
coe4rsd = -3.d0*coe4rs/rs + 3.d0*coe4/rs**2 &
|
|
+ 27.d0/64.d0/rs**4*(2d0**(-2.d0/3.d0)* &
|
|
dpold(2d0**(1.d0/3.d0)*rs)+d2antid &
|
|
+ 2.d0/5.d0/rs**3*cf**2) &
|
|
-9.d0/64.d0/rs**3*(2d0**(-1.d0/3.d0) &
|
|
* dpoldd(2d0**(1.d0/3.d0)*rs) &
|
|
+d2antidd - 6.d0/5.d0*cf**2/rs**4)
|
|
coe4zd = 0.d0
|
|
|
|
! coe5rsd = 3.d0/rs**2*coe5-3.d0/rs*coe5rs
|
|
! > +27.d0/40.d0/(sqrt2pi*rs**4)*(
|
|
! $ ((1.d0+z)/2.d0)**(5.d0/3.d0)*dpold(rs*(2/(1.d0+z))**
|
|
! $ (1.d0/3.d0))+((1.d0-z)/2.d0)**(5.d0/3.d0)*
|
|
! $ dpold(rs*(2/(1.d0-z))**(1.d0/3.d0))+(1.d0-z**2)*
|
|
! $ d3antid)-9.d0/40.d0/(sqrt2pi*rs**3)*(
|
|
! $ ((1.d0+z)/2.d0)**(4.d0/3.d0)*dpoldd(rs*(2/(1.d0+z))**
|
|
! $ (1.d0/3.d0))+((1.d0-z)/2.d0)**(4.d0/3.d0)*
|
|
! $ dpoldd(rs*(2/(1.d0-z))**(1.d0/3.d0))+(1.d0-z**2)*
|
|
! $ d3antidd)
|
|
! Case were z=0
|
|
coe5rsd = -3.d0*coe5rs/rs + 3.d0*coe5/rs**2 &
|
|
+27.d0/(40.d0*sqrt2pi*rs**4)* &
|
|
(2d0**(-2.d0/3.d0)*dpold(2d0**(1.d0/3.d0)*rs)+d3antid) &
|
|
-9.d0/(40.d0*sqrt2pi*rs**3)*(2d0**(-1.d0/3.d0)* &
|
|
dpoldd(2d0**(1.d0/3.d0)*rs)+d3antidd)
|
|
coe5zd = 0.d0
|
|
!SCF
|
|
|
|
end if
|
|
|
|
! call ecPW(rs,z,ec,ecd,ecz)
|
|
!SCD
|
|
call ecPW(rs,z,ec,ecd,ecz,ecdd,eczd)
|
|
!SCF
|
|
|
|
a1=4.d0*b0**6*coe3+b0**8*coe5
|
|
a1rs=24.d0*adib*b0**5*coe3+4.d0*b0**6*coe3rs+8.d0*adib*b0**7* coe5+b0**8*coe5rs
|
|
a1z=4.d0*b0**6*coe3z+b0**8*coe5z
|
|
|
|
a2=4.d0*b0**6*coe2+b0**8*coe4+6.d0*b0**4*ec
|
|
a2rs=24.d0*adib*b0**5*coe2+4.d0*b0**6*coe2rs+8.d0*adib*b0**7* &
|
|
coe4+b0**8*coe4rs+24.d0*adib*b0**3*ec+6.d0*b0**4*ecd
|
|
a2z=4.d0*b0**6*coe2z+b0**8*coe4z+6.d0*b0**4*ecz
|
|
|
|
a3=b0**8*coe3
|
|
a3rs=8.d0*adib*b0**7*coe3+b0**8*coe3rs
|
|
a3z=b0**8*coe3z
|
|
|
|
a4=b0**6*(b0**2*coe2+4.d0*ec)
|
|
a4rs=8.d0*adib*b0**7*coe2+b0**8*coe2rs+24.d0*adib*b0**5*ec+ &
|
|
4.d0*b0**6*ecd
|
|
a4z=b0**6*(b0**2*coe2z+4.d0*ecz)
|
|
|
|
a5=b0**8*ec
|
|
a5rs=8.d0*adib*b0**7*ec+b0**8*ecd
|
|
a5z=b0**8*ecz
|
|
!SCD
|
|
a1rsd = 120.d0*adib**2*b0**4*coe3 + 48.d0*adib*b0**5*coe3rs &
|
|
+ 4.d0*b0**6*coe3rsd + 56.d0*adib**2*b0**6*coe5 &
|
|
+ 16.d0*adib*b0**7*coe5rs + b0**8*coe5rsd
|
|
! a1zd = 4.d0*b0**6*coe3zd+b0**8*coe5zd
|
|
a1zd = 0.d0
|
|
!
|
|
a2rsd = 120.d0*adib**2*b0**4*coe2 + 48.d0*adib*b0**5*coe2rs &
|
|
+ 4.d0*b0**6*coe2rsd + 56.d0*b0**6*adib**2*coe4 &
|
|
+ 16.d0*b0**7*adib*coe4rs + b0**8*coe4rsd &
|
|
+ 72.d0*b0**2*adib**2*ec + 48.d0*b0**3*adib*ecd &
|
|
+ 6.d0*b0**4*ecdd
|
|
! a2zd = 4.d0*b0**6*coe2zd+b0**8*coe4zd+6.d0*b0**4*eczd
|
|
a2zd = 0.d0
|
|
!
|
|
a3rsd = 56.d0*adib**2*b0**6*coe3 + 16.d0*adib*b0**7*coe3rs &
|
|
+ b0**8*coe3rsd
|
|
! a3zd = b0**8*coe3zd
|
|
a3zd = 0.d0
|
|
!
|
|
a4rsd = 56.d0*adib**2*b0**6*coe2 + 16.d0*adib*b0**7*coe2rs &
|
|
+ b0**8*coe2rsd + 120.d0*adib**2*b0**4*ec &
|
|
+ 48.d0*adib*b0**5*ecd + 4.d0*b0**6*ecdd
|
|
! a4zd = b0**6*(b0**2*coe2zd+4.d0*eczd)
|
|
a4zd = 0.d0
|
|
!
|
|
a5rsd = 56.d0*adib**2*b0**6*ec + 16.d0*adib*b0**7*ecd &
|
|
+ b0**8*ecdd
|
|
! a5zd=b0**8*eczd
|
|
a5zd= 0.d0
|
|
!SCF
|
|
|
|
x=mu*dsqrt(rs)/phi
|
|
|
|
eclr=(phi**3*Qrpa(x)+a1*mu**3+a2*mu**4+a3*mu**5+ &
|
|
a4*mu**6+a5*mu**8)/((1.d0+b0**2*mu**2)**4)
|
|
|
|
eclrrs=-4.d0/(1.d0+b0**2*mu**2)*2.d0*adib*b0*mu**2*eclr+ &
|
|
1.d0/((1.d0+b0**2*mu**2)**4)*(phi**2*mu/(2.d0*sqrt(rs)) &
|
|
*Qrpad(x)+ &
|
|
a1rs*mu**3+a2rs*mu**4+a3rs*mu**5+a4rs*mu**6+a5rs*mu**8)
|
|
!SCD
|
|
! u=
|
|
! > (phi**2*mu/(2.d0*sqrt(rs))
|
|
! > *Qrpad(x)+
|
|
! > a1rs*mu**3+a2rs*mu**4+a3rs*mu**5+a4rs*mu**6+a5rs*mu**8)
|
|
! du=
|
|
! > (-phi**2*mu/(4.d0*rs**(3.d0/2.d0))*Qrpad(x)
|
|
! > +mu**2*phi/(4.d0*rs)*Qrpadd(x)*+
|
|
! > a1rsd*mu**3+a2rsd*mu**4+a3rsd*mu**5+a4rsd*mu**6+a5rsd*mu**8)
|
|
! v = (1.d0+b0**2*mu**2)**4
|
|
! dv= 8.d0*(1.d0+(b0*mu)**2)**3*b0*adib*mu**2
|
|
! eclrrsd= -8.d0*adib*b0*mu**2*eclrrs/(1.d0+b0**2*mu**2)
|
|
! > -8.d0*(adib*mu)**2/(1.d0+b0**2*mu**2)*eclr
|
|
! > +16.d0*(adib*mu)**4*rs**2/((1.d0+(b0*mu)**2))**2*eclr
|
|
! > +du/v-u*dv/v**2
|
|
u = (phi**3*Qrpa(x)+a1*mu**3+a2*mu**4+a3*mu**5+a4*mu**6+a5*mu**8)
|
|
du = (phi**2*mu/(2.d0*sqrt(rs))*Qrpad(x)+a1rs*mu**3+a2rs*mu**4 &
|
|
+a3rs*mu**5+a4rs*mu**6+a5rs*mu**8)
|
|
ddu = - phi**2*mu/(4.d0*rs**(3.d0/2.d0))*Qrpad(x) &
|
|
+ phi*mu**2/(4.d0*rs)*Qrpadd(x)+a1rsd*mu**3+a2rsd*mu**4 &
|
|
+ a3rsd*mu**5+a4rsd*mu**6+a5rsd*mu**8
|
|
v = ((1.d0+b0**2*mu**2)**4)
|
|
dv = 8.d0*(1.d0+b0**2*mu**2)**3*(adib**2*mu**2*rs)
|
|
ddv = 48.d0*(1.d0+b0**2*mu**2)**2*(adib**2*mu**2*rs)**2 &
|
|
+ 8.d0*(1.d0+b0**2*mu**2)**3*(adib**2*mu**2)
|
|
! eclrrsd = ddu/v - du*dv/v**2 - dv/v*eclrrs
|
|
! > - eclr*(ddv/v - (dv/v)**2)
|
|
eclrrsd = ddu/v - 2.d0*du*dv/v**2 - u*ddv/v**2 &
|
|
+ 2.d0*u*dv**2/v**3
|
|
|
|
!SCF
|
|
|
|
|
|
if(z.eq.1.d0) then
|
|
vclrup=eclr-rs/3.d0*eclrrs
|
|
vclrdown=0.d0
|
|
!SCD
|
|
vclrupd = eclrrs-1.d0/3.d0*eclrrs -rs/3.d0*eclrrsd
|
|
vclrdownd = 0.d0
|
|
!SCF
|
|
elseif(z.eq.-1.d0) then
|
|
vclrup=0.d0
|
|
vclrdown=eclr-rs/3.d0*eclrrs
|
|
!SCD
|
|
vclrupd = 0.d0
|
|
vclrdownd = eclrrs-1.d0/3.d0*eclrrs &
|
|
-rs/3.d0*eclrrsd
|
|
!SCF
|
|
else
|
|
|
|
eclrz=(phi**2*((1.d0+z)**(-1.d0/3.d0)-(1.d0-z)**(-1.d0/3.d0)) &
|
|
*Qrpa(x)-phi*Qrpad(x)*mu*sqrt(rs)*((1.d0+z)**(-1.d0/3.d0) &
|
|
-(1.d0-z)**(-1.d0/3.d0))/3.d0+ &
|
|
a1z*mu**3+a2z*mu**4+a3z*mu**5+ &
|
|
a4z*mu**6+a5z*mu**8)/((1.d0+b0**2*mu**2)**4)
|
|
!SCD
|
|
eclrzd=0.d0
|
|
!CSF
|
|
|
|
vclrup=eclr-rs/3.d0*eclrrs-(z-1.d0)*eclrz
|
|
vclrdown=eclr-rs/3.d0*eclrrs-(z+1.d0)*eclrz
|
|
!SCD
|
|
vclrupd = 2.d0/3.d0*eclrrs - rs/3.d0*eclrrsd
|
|
vclrdownd = 2.d0/3.d0*eclrrs - rs/3.d0*eclrrsd
|
|
!SCF
|
|
end if
|
|
return
|
|
end
|
|
|
|
|
|
double precision function g0f(x)
|
|
!cc on-top pair-distribution function
|
|
!cc Gori-Giorgi and Perdew, PRB 64, 155102 (2001)
|
|
!cc x -> rs
|
|
implicit none
|
|
double precision C0f,D0f,E0f,F0f,x
|
|
C0f = 0.0819306d0
|
|
D0f = 0.752411d0
|
|
E0f = -0.0127713d0
|
|
F0f = 0.00185898d0
|
|
double precision :: exp_d0fx
|
|
if(dabs(D0f*x).lt.50.d0)then
|
|
exp_d0fx = dexp(-dabs(D0f)*x)
|
|
else
|
|
exp_d0fx = 0.d0
|
|
endif
|
|
g0f=(1.d0-(0.7317d0-D0f)*x+C0f*x**2+E0f*x**3+ &
|
|
F0f*x**4)*exp_d0fx/2.d0
|
|
return
|
|
end
|
|
|
|
double precision function g0d(rs)
|
|
!cc derivative of on-top pair-distribution function
|
|
!cc Gori-Giorgi and Perdew, PRB 64, 155102 (2001)
|
|
implicit none
|
|
double precision Bg0,Cg0,Dg0,Eg0,Fg0,rs,expsum
|
|
Cg0 = 0.0819306d0
|
|
Fg0 = 0.752411d0
|
|
Dg0 = -0.0127713d0
|
|
Eg0 = 0.00185898d0
|
|
Bg0 =0.7317d0-Fg0
|
|
if(dabs(Fg0*rs).lt.50.d0)then
|
|
expsum=dexp(-Fg0*rs)
|
|
else
|
|
expsum = 0.d0
|
|
endif
|
|
g0d=(-Bg0+2d0*Cg0*rs+3d0*Dg0*rs**2+4d0*Eg0*rs**3)/2.d0 &
|
|
*expsum &
|
|
- (Fg0*(1d0 - Bg0*rs + Cg0*rs**2 + Dg0*rs**3 + Eg0*rs**4))/ &
|
|
2.d0*expsum
|
|
return
|
|
end
|
|
!SCD
|
|
double precision function g0dd(rs)
|
|
!cc derivative of g0d
|
|
implicit none
|
|
double precision Bg0,Cg0,Dg0,Eg0,Fg0,rs,expsum
|
|
Cg0 = 0.0819306d0
|
|
Fg0 = 0.752411d0
|
|
Dg0 = -0.0127713d0
|
|
Eg0 = 0.00185898d0
|
|
Bg0 = 0.7317d0-Fg0
|
|
if(dabs(Fg0*rs).lt.50.d0)then
|
|
expsum=dexp(-Fg0*rs)
|
|
else
|
|
expsum=0.d0
|
|
endif
|
|
g0dd = (2.d0*Cg0+6.d0*Dg0*rs+12.d0*Eg0*rs**2)/2.d0* &
|
|
expsum &
|
|
- (-Bg0+2.d0*Cg0*rs+3.d0*Dg0*rs**2+4.d0*Eg0*rs**3)*Fg0* &
|
|
expsum &
|
|
+ (1.d0-Bg0*rs+Cg0*rs**2+Dg0*rs**3+Eg0*rs**4)*Fg0**2* &
|
|
expsum/(2.d0)
|
|
return
|
|
end
|
|
!SCF
|
|
|
|
double precision function dpol(rs)
|
|
implicit none
|
|
double precision cf,pi,rs,p2p,p3p
|
|
pi=dacos(-1.d0)
|
|
cf=(9.d0*pi/4.d0)**(1.d0/3.d0)
|
|
p2p = 0.04d0
|
|
p3p = 0.4319d0
|
|
dpol=2.d0**(5.d0/3.d0)/5.d0*cf**2/rs**2*(1.d0+(p3p-0.454555d0)*rs) &
|
|
/(1.d0+p3p*rs+p2p*rs**2)
|
|
return
|
|
end
|
|
|
|
double precision function dpold(rs)
|
|
implicit none
|
|
double precision cf,pi,rs,p2p,p3p
|
|
pi=dacos(-1.d0)
|
|
cf=(9.d0*pi/4.d0)**(1.d0/3.d0)
|
|
p2p = 0.04d0
|
|
p3p = 0.4319d0
|
|
dpold=2.d0**(5.d0/3.d0)/5.d0*cf**2* &
|
|
(-2.d0 + (0.454555d0 - 4.d0*p3p)*rs + &
|
|
(-4.d0*p2p + &
|
|
(0.90911d0 - 2.d0*p3p)*p3p)*rs**2 &
|
|
+ p2p*(1.363665d0 - 3.d0*p3p)* &
|
|
rs**3)/ &
|
|
(rs**3*(1.d0 + p3p*rs + p2p*rs**2)**2)
|
|
return
|
|
end
|
|
|
|
!SCD
|
|
double precision function dpoldd(rs)
|
|
implicit none
|
|
double precision cf,pi,rs,p2p,p3p,p4p
|
|
pi=dacos(-1.d0)
|
|
cf=(9.d0*pi/4.d0)**(1.d0/3.d0)
|
|
p2p = 0.04d0
|
|
p3p = 0.4319d0
|
|
p4p = 0.454555d0
|
|
dpoldd = 4.d0/5.d0*2d0**(2.d0/3.d0)*cf**2*( &
|
|
9.d0*p2p*rs**2 + 8.d0*p3p**2*rs**4*p2p &
|
|
+ 6.d0*p3p*rs**5*p2p**2 - 3.d0*rs**3*p4p*p3p**2 &
|
|
- 6.d0*rs**5*p4p*p2p**2 - 3.d0*rs**2*p3p*p4p &
|
|
- 3.d0*rs**3*p2p*p4p + 10.d0*p2p**2*rs**4 &
|
|
+ 9.d0*p3p*rs + 9.d0*p3p**2*rs**2 + 3.d0 &
|
|
+ 3.d0*p3p**3*rs**3 - 8.d0*rs**4*p2p*p3p*p4p &
|
|
+ 18.d0*p3p*p2p*rs**3 - rs*p4p)/ &
|
|
(rs**4*(1.d0+p3p*rs+p2p*rs**2)**3)
|
|
return
|
|
end
|
|
!SCF
|
|
double precision function Qrpa(x)
|
|
implicit none
|
|
double precision pi,a2,b2,c2,d2,x,Acoul
|
|
pi=dacos(-1.d0)
|
|
Acoul=2.d0*(dlog(2.d0)-1.d0)/pi**2
|
|
a2 = 5.84605d0
|
|
c2 = 3.91744d0
|
|
d2 = 3.44851d0
|
|
b2=d2-3.d0/(2.d0*pi*Acoul)*(4.d0/(9.d0*pi))**(1.d0/3.d0)
|
|
Qrpa=Acoul*dlog((1.d0+a2*x+b2*x**2+c2*x**3)/(1.d0+a2*x+d2*x**2))
|
|
return
|
|
end
|
|
|
|
double precision function Qrpad(x)
|
|
implicit none
|
|
double precision pi,a2,b2,c2,d2,x,Acoul
|
|
pi=dacos(-1.d0)
|
|
Acoul=2.d0*(dlog(2.d0)-1.d0)/pi**2
|
|
a2 = 5.84605d0
|
|
c2 = 3.91744d0
|
|
d2 = 3.44851d0
|
|
b2=d2-3.d0/(2.d0*pi*Acoul)*(4.d0/(9.d0*pi))**(1.d0/3.d0)
|
|
Qrpad=Acoul*((x*(b2*(2.d0 + a2*x) + &
|
|
c2*x*(3.d0 + 2.d0*a2*x) + &
|
|
d2*(-2.d0 - a2*x + c2*x**3)))/ &
|
|
((1.d0 + a2*x + d2*x**2)* &
|
|
(1.d0 + a2*x + b2*x**2 + c2*x**3)))
|
|
return
|
|
end
|
|
!SCD
|
|
double precision function Qrpadd(x)
|
|
implicit none
|
|
double precision pi,a2,b2,c2,d2,x,Acoul
|
|
double precision uQ,duQ,dduQ,vQ,dvQ,ddvQ
|
|
pi=dacos(-1.d0)
|
|
Acoul=2.d0*(dlog(2.d0)-1.d0)/pi**2
|
|
a2 = 5.84605d0
|
|
c2 = 3.91744d0
|
|
d2 = 3.44851d0
|
|
b2=d2-3.d0/(2.d0*pi*Acoul)*(4.d0/(9.d0*pi))**(1.d0/3.d0)
|
|
uQ = 1.d0 + a2*x + b2*x**2 + c2*x**3
|
|
duQ = a2 + 2.d0*b2*x + 3.d0*c2*x**2
|
|
dduQ= 2.d0*b2 + 6.d0*c2*x
|
|
vQ = 1.d0 + a2*x + d2*x**2
|
|
dvQ = a2 + 2.d0*d2*x
|
|
ddvQ= 2.d0*d2
|
|
Qrpadd = Acoul*(dduQ/uQ - (duQ/uQ)**2 -ddvQ/vQ +(dvQ/vQ)**2)
|
|
return
|
|
end
|
|
!SCF
|
|
|
|
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
|
|
! correlation energy and its derivative w.r.t. rs and z at mu=infinity
|
|
! Perdew & Wang PRB 45, 13244 (1992)
|
|
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
|
|
! subroutine ecPW(x,y,ec,ecd,ecz)
|
|
!SCD
|
|
subroutine ecPW(x,y,ec,ecd,ecz,ecdd,eczd)
|
|
!SCF
|
|
! in Hartree; ec=ec(rs,zeta)
|
|
! x -> rs; y -> zeta
|
|
!cc ecd is d/drs ec
|
|
!cc ecz is d/dz ec
|
|
implicit none
|
|
double precision pi,f02,ff,x,y,ec,ecd,ec0,ec0d,ec1,ec1d,aaa,G,Gd,alfac,alfacd,ecz
|
|
!SCD
|
|
double precision alfacdd,ec0dd,ecdd,ec1dd,Gdd,eczd
|
|
!SCF
|
|
pi=dacos(-1.d0)
|
|
|
|
f02=4.d0/(9.d0*(2.d0**(1.d0/3.d0)-1.d0))
|
|
|
|
ff=((1.d0+y)**(4.d0/3.d0)+(1.d0-y)**(4.d0/3.d0)- &
|
|
2.d0)/(2.d0**(4.d0/3.d0)-2.d0)
|
|
|
|
aaa=(1.d0-dlog(2.d0))/pi**2
|
|
call GPW(x,aaa,0.21370d0,7.5957d0,3.5876d0, &
|
|
1.6382d0,0.49294d0,G,Gd,Gdd)
|
|
ec0=G
|
|
ec0d=Gd
|
|
ec0dd=Gdd
|
|
|
|
aaa=aaa/2.d0
|
|
call GPW(x,aaa,0.20548d0,14.1189d0,6.1977d0, &
|
|
3.3662d0,0.62517d0,G,Gd,Gdd)
|
|
ec1=G
|
|
ec1d=Gd
|
|
ec1dd=Gdd
|
|
call GPW(x,0.016887d0,0.11125d0,10.357d0,3.6231d0,0.88026d0,0.49671d0,G,Gd,Gdd)
|
|
alfac=-G
|
|
alfacd=-Gd
|
|
alfacdd=-Gdd
|
|
|
|
ec=ec0+alfac*ff/f02*(1.d0-y**4)+(ec1-ec0)*ff*y**4
|
|
ecd=ec0d+alfacd*ff/f02*(1.d0-y**4)+(ec1d-ec0d)* &
|
|
ff*y**4
|
|
ecz=alfac*(-4.d0*y**3)*ff/f02+alfac*(1.d0-y**4)/f02* &
|
|
4.d0/3.d0*((1.d0+y)**(1.d0/3.d0)-(1.d0-y)**(1.d0/3.d0))/ &
|
|
(2.d0**(4.d0/3.d0)-2.d0)+(ec1-ec0)*(4.d0*y**3*ff+ &
|
|
4.d0/3.d0*((1.d0+y)**(1.d0/3.d0)-(1.d0-y)**(1.d0/3.d0))/ &
|
|
(2.d0**(4.d0/3.d0)-2.d0)*y**4)
|
|
!SCD
|
|
ecdd = ec0dd + alfacdd*ff/f02*(1.D0-y**4) + (ec1dd - ec0dd)*ff*y**4
|
|
|
|
eczd = 0.d0
|
|
!SCF
|
|
|
|
return
|
|
end
|
|
|
|
subroutine GPW(x,Ac,alfa1,beta1,beta2,beta3,beta4,G,Gd,Gdd)
|
|
!SCF
|
|
!cc Gd is d/drs G
|
|
!cc Gdd is d/drs Gd
|
|
implicit none
|
|
double precision G,Gd,Ac,alfa1,beta1,beta2,beta3,beta4,x
|
|
!SCD
|
|
double precision f32,f34,f12,f14,Gdd
|
|
double precision A,dA,ddA,B
|
|
!SCF
|
|
double precision sqrtx
|
|
sqrtx=dsqrt(x)
|
|
G=-2.d0*Ac*(1.d0+alfa1*x)*dlog(1.d0+1.d0/(2.d0* &
|
|
Ac*(beta1*x**0.5d0+ &
|
|
beta2*x+beta3*x**1.5d0+beta4*x**2)))
|
|
Gd=(1.d0+alfa1*x)*(beta2+beta1/(2.d0*sqrtx)+3.d0*beta3* &
|
|
sqrtx/2.d0+2.d0*beta4*x)/((beta1*sqrtx+beta2*x+ &
|
|
beta3*x**(3.d0/2.d0)+beta4*x**2)**2*(1.d0+1.d0/ &
|
|
(2.d0*Ac*(beta1*sqrtx+beta2*x+beta3*x**(3.d0/2.d0)+&
|
|
beta4*x**2))))-2.d0*Ac*alfa1*dlog(1.d0+1.d0/(2.d0*Ac*&
|
|
(beta1*sqrtx+beta2*x+beta3*x**(3.d0/2.d0)+&
|
|
beta4*x**2)))
|
|
!SCD
|
|
f12=(1.d0)/(2.d0)
|
|
f14=(1.d0)/(4.d0)
|
|
f32=(3.d0)/(2.d0)
|
|
f34=(3.d0)/(4.d0)
|
|
A = beta1*sqrtx + beta2*x + beta3*x**(3.d0/2.d0) + beta4*x**2
|
|
dA = f12*beta1/sqrtx + beta2 + f32*beta3*sqrtx + 2.d0*beta4*x
|
|
ddA = -f14*beta1*x**(-f32) + f34*beta3/sqrtx + 2.d0*beta4
|
|
B = 1.d0 + 1.d0/(2.d0*Ac*A)
|
|
Gdd = 2.d0*alfa1*dA/(A**2*B) &
|
|
- 2.d0*(1.d0+alfa1*x)*dA**2/(A**3*B) &
|
|
+ (1.d0+alfa1*x)*ddA/(A**2*B) &
|
|
+ (1.d0+alfa1*x)*dA**2/(A**4*B**2*Ac*2.d0)
|
|
return
|
|
end
|