LCPQ/qp2
9
1
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-12-22 03:23:29 +01:00
Code Releases Activity
c960dfd245
qp2/src/dft_utils_one_e/exc_sr_lda.irp.f

1008 lines
33 KiB
Fortran
Raw Normal View History

Initial commit
2019-01-25 11:39:31 +01:00
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
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)
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
ex_a = - (rho_a_2*(z24*rho_a_2/pi)**f13) * (z3/z8-a*(sqpi*derf(f12/a)+(z2*a-z4*a3)*dexp(-f14/a2)-z3*a+z4*a3))
vx_a = -(z3*rho_a_2/pi)**f13 + z2*a*mu/pi*(dexp(-f14/a2)-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)
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
ex_b = - (rho_b_2*(z24*rho_b_2/pi)**f13)*(z3/z8-a*(sqpi*derf(f12/a)+(z2*a-z4*a3)*dexp(-f14/a2)-z3*a+z4*a3))
vx_b = -(z3*rho_b_2/pi)**f13+ z2*a*mu/pi*(dexp(-f14/a2)-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
fak=2.540118935556d0*dexp(-eta*a*a)
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
fak=2.540118935556d0*dexp(-eta*a*a)
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
d2anti=(q1a*rs+q2a*rs**2)*exp(-abs(q3a)*rs)/rs**2
d3anti=(t1a*rs+t2a*rs**2)*exp(-abs(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)/(sqrt(2.d0*pi)*rs**3)
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)
coe5=-9.d0/40.d0/(sqrt(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/(sqrt(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*sqrt(rs)/phi.lt.0.d0)then
print*,'phi',phi
print*,'mu ',mu
print*,'rs ',rs
stop -1
endif
eclr=(phi**3*Qrpa(mu*sqrt(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
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(-rs*t3a)
!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/sqrt(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*sqrt(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
g0f=(1.d0-(0.7317d0-D0f)*x+C0f*x**2+E0f*x**3+ &
F0f*x**4)*exp(-abs(D0f)*x)/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
expsum=exp(-Fg0*rs)
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
expsum=exp(-Fg0*rs)
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*(log(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)
!if(((1.d0+a2*x+b2*x**2+c2*x**3)/(1.d0+a2*x+d2*x**2)).le.0.d0)then
! print*,(1.d0+a2*x+b2*x**2+c2*x**3)/(1.d0+a2*x+d2*x**2)
! print*,(1.d0+a2*x+b2*x**2+c2*x**3),(1.d0+a2*x+d2*x**2)
! print*,x
! pause
!endif
!Qrpa=Acoul*log(dabs((1.d0+a2*x+b2*x**2+c2*x**3)/(1.d0+a2*x+d2*x**2)))
Qrpa=Acoul*log((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*(log(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*(log(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-log(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)
!SCD
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=sqrt(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
Copy Permalink