9
1
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-10-09 17:33:05 +02:00
qp2/src/dft_utils_one_e/exc_sr_pbe.irp.f
2019-01-25 11:39:31 +01:00

514 lines
15 KiB
Fortran

subroutine ec_pbe_sr(mu,rhoc,rhoo,sigmacc,sigmaco,sigmaoo,ec,vrhoc,vrhoo,vsigmacc,vsigmaco,vsigmaoo)
BEGIN_DOC
! Short-range PBE correlation energy functional for erf interaction
!
! input : ==========
!
! mu = range separated parameter
!
! rhoc, rhoo = total density and spin density
!
! sigmacc = square of the gradient of the total density
!
! sigmaco = square of the gradient of the spin density
!
! sigmaoo = scalar product between the gradient of the total density and the one of the spin density
!
! output: ==========
!
! ec = correlation energy
!
! all variables v** are energy derivatives with respect to components of the density
!
! vrhoc = derivative with respect to the total density
!
! vrhoo = derivative with respect to spin density
!
! vsigmacc = derivative with respect to the square of the gradient of the total density
!
! vsigmaco = derivative with respect to scalar product between the gradients of total and spin densities
!
! vsigmaoo = derivative with respect to the square of the gradient of the psin density
END_DOC
include 'constants.include.F'
implicit none
double precision, intent(in) :: rhoc,rhoo,mu
double precision, intent(in) :: sigmacc,sigmaco,sigmaoo
double precision, intent(out) :: ec
double precision, intent(out) :: vrhoc,vrhoo
double precision, intent(out) :: vsigmacc,vsigmaco,vsigmaoo
double precision tol
parameter(tol=1d-12)
character(len=30) namedummy
double precision eccerflda
double precision vrhoccerflda
double precision vrhoocerflda
double precision ecclda
double precision vrhocclda
double precision vrhooclda
integer i,igrad
double precision rho,drho2,rhoa,rhob
double precision ecerflda,decerfldadrho
double precision eclda,decldadrho
double precision ecerfpbe,decerfpbedrho,decerfpbedrhoo
double precision decerfpbeddrho2
double precision arglog,arglogs,arglogss,alpha,beta,betas,gamma
double precision Aa,Ab,Ac,Aas,tq,tqs,tqss,decerfpur,decpur
double precision t1,t2,t3,t4,t5,t6,t7,t8,t9,t10
double precision t11,t12,t13,t14,t15,t16,t17,t18,t19
double precision zeta,phi,phi2,phi3,phi4,phis,arglogsc
double precision dlogarglog
double precision, parameter :: f13=0.333333333333333d0
! Parameter of the modified interaction
ec = 0.d0
vrhoc = 0.d0
vrhoo = 0.d0
vsigmacc = 0.d0
vsigmaco = 0.d0
vsigmaoo = 0.d0
! First-type gradient functional
igrad=1
alpha=2.78d0
gamma=3.1091d-2
! test on density
if (dabs(rhoc).lt.tol) return
double precision :: vc_a,vc_b
! Spin polarisation
rhoa=max((rhoc+rhoo)*.5d0,1.0d-15)
rhob=max((rhoc-rhoo)*.5d0,1.0d-15)
call ec_lda_sr(mu,rhoa,rhob,eccerflda,vc_a,vc_b)
ecerflda = eccerflda
vrhoccerflda = 0.5d0 * (vc_a + vc_b)
vrhoocerflda = 0.5d0 * (vc_a - vc_b)
! Density
rho = rhoc
! Square of density gradient
drho2 = sigmacc
zeta = (rhoa-rhob)/(rhoa+rhob)
! LDA energy density
double precision :: vc_a_lda,vc_b_lda
call ec_lda(rhoa,rhob,ecclda,vc_a_lda,vc_b_lda)
eclda = ecclda
if ((ecerflda/eclda).le.0d0) then
beta=0d0
else
beta=6.6725d-2*(ecerflda/eclda)**alpha
endif
phi=((1d0+zeta)**(2d0/3d0)+(1d0-zeta)**(2d0/3d0))/2d0
phi2=phi*phi
phi3=phi2*phi
phi4=phi3*phi
tq=drho2*6.346820607d-2*rho**(-7d0/3d0)/phi2
Ab=dexp(-ecerflda/(rho*gamma*phi3))-1d0
if (dabs(Ab).le.dabs(beta*tol)) then
ecerfpbe=ecerflda
else
Aa=beta/(gamma*Ab)
Ac=1d0+Aa*tq+Aa**2*tq**2
if (Aa.lt.tol) Aa=tol
arglog=1d0+beta*(1d0-1d0/Ac)/(gamma*Aa)
ecerfpbe=ecerflda+rho*phi3*gamma*dlog(arglog)
end if
ec = ecerfpbe
! Derive
! LDA energy density derivative
decerfldadrho = vrhoccerflda
decldadrho = 0.5d0 * (vc_a_lda+vc_b_lda)
decerfpur=(decerfldadrho-ecerflda/rho)/rho
decpur=(decldadrho-eclda/rho)/rho
betas=alpha*beta*(decerfpur*rho/ecerflda-decpur*rho/eclda)
phis=((rhoa - rhob)*((rhoa/(rhoa + rhob))**f13 - (rhob/(rhoa + rhob))**f13))/(3d0*2d0**f13*(rhoa/(rhoa + rhob))**f13*(rhob/(rhoa + rhob))**f13*(rhoa + rhob)**2)
if (dabs(Ab).le.dabs(beta*tol)) then
decerfpbedrho=decerfldadrho
else
Aas=betas/(gamma*Ab)+Aa*(1d0+1d0/Ab)*(decerfpur/phi3-3d0*phis*ecerflda/(rho*phi4))/gamma
tqs=-7d0*tq/(3d0*rho)-2d0*tq*phis/phi
arglogs=betas*tq*(1d0+Aa*tq)/(Ac*gamma)+beta*tqs*(1d0+Aa*tq)/(Ac*gamma)-beta*tq*Aa*tq*(Aas*tq+Aa*tqs)*(2d0+Aa*tq)/(Ac**2*gamma)
dlogarglog=dlog(arglog)
decerfpbedrho=decerfldadrho+gamma*(phi3*dlogarglog+3d0*rho*phis*phi2*dlogarglog+rho*phi3*arglogs/arglog)
end if
if (dabs(Ab).le.dabs(beta*tol)) then
decerfpbeddrho2=0.0d0
else
arglogsc=Ab*(Aa+2d0*Aa*Aa*tq)/(Ac*Ac)
tqss=6.346820607d-2*rho**(-7d0/3d0)/phi2
arglogss=tqss*arglogsc
decerfpbeddrho2=rho*gamma*phi3*arglogss/arglog
end if
! LDA energy density derivative
decerfldadrho = vrhoocerflda
decldadrho = 0.5d0 * (vc_a_lda-vc_b_lda)
decerfpur=decerfldadrho/rho
decpur=decldadrho/rho
betas=alpha*beta*(decerfpur*rho/ecerflda-decpur*rho/eclda)
phis=(rhob*(rhoa/(rhoa + rhob))**(2d0*f13)-rhoa*(rhob/(rhoa + rhob))**(2d0*f13))/(3d0*2d0**f13*rhoa*rhob)
if (dabs(Ab).le.dabs(beta*tol)) then
decerfpbedrhoo=decerfldadrho
else
Aas=betas/(gamma*Ab)+Aa*(1d0+1d0/Ab)*(decerfpur/phi3-3d0*phis*ecerflda/(rho*phi4))/gamma
tqs=-2d0*tq*phis/phi
arglogs=betas*tq*(1d0+Aa*tq)/(Ac*gamma)+beta*tqs*(1d0+Aa*tq)/(Ac*gamma)-beta*tq*Aa*tq*(Aas*tq+Aa*tqs)*(2d0+Aa*tq)/(Ac**2*gamma)
decerfpbedrhoo=decerfldadrho+gamma*(3d0*rho*phis*phi2*dlog(arglog)+rho*phi3*arglogs/arglog)
end if
! derivatives
vrhoc = vrhoc + decerfpbedrho
vrhoo = vrhoo + decerfpbedrhoo
vsigmacc = vsigmacc + decerfpbeddrho2
end
subroutine ex_pbe_sr(mu,rho_a,rho_b,grd_rho_a_2,grd_rho_b_2,grd_rho_a_b,ex,vx_rho_a,vx_rho_b,vx_grd_rho_a_2,vx_grd_rho_b_2,vx_grd_rho_a_b)
BEGIN_DOC
!mu = range separation parameter
!rho_a = density alpha
!rho_b = density beta
!grd_rho_a_2 = (gradient rho_a)^2
!grd_rho_b_2 = (gradient rho_b)^2
!grd_rho_a_b = (gradient rho_a).(gradient rho_b)
!ex = exchange energy density at the density and corresponding gradients of the density
!vx_rho_a = d ex / d rho_a
!vx_rho_b = d ex / d rho_b
!vx_grd_rho_a_2 = d ex / d grd_rho_a_2
!vx_grd_rho_b_2 = d ex / d grd_rho_b_2
!vx_grd_rho_a_b = d ex / d grd_rho_a_b
END_DOC
implicit none
! input
double precision, intent(in) :: mu,rho_a, rho_b
double precision, intent(in) :: grd_rho_a_2, grd_rho_b_2, grd_rho_a_b
! output
double precision, intent(out) :: ex
double precision, intent(out) :: vx_rho_a, vx_rho_b
double precision, intent(out) :: vx_grd_rho_a_2, vx_grd_rho_b_2, vx_grd_rho_a_b
! function
double precision berf
double precision dberfda
! local
double precision, parameter :: tol=1d-12
double precision, parameter :: f13=0.333333333333333d0
double precision exerflda,vxerflda_a,vxerflda_b
double precision dexerfldadrho
double precision exerfpbe_a, exerfpbe_b
double precision dexerfpbedrho_a, dexerfpbedrho_b
double precision dexerfpbeddrho2_a, dexerfpbeddrho2_b
double precision rho,drho2
double precision rho_a_2, rho_b_2
double precision t1,t2,t3,t4
double precision kappa,sq,sqs,sqss,fx,fxs,ksig
! Parameter of the modified interaction
! initialization
ex=0.d0
vx_rho_a=0.d0
vx_rho_b=0.d0
vx_grd_rho_a_2=0.d0
vx_grd_rho_b_2=0.d0
vx_grd_rho_a_b=0.d0
! spin scaling relation Ex[rho_a,rho_b] = (1/2) (Ex[2rho_a,2rho_a] + Ex[2rho_b,2rho_b])
! two times spin alpha density
rho = max(rho_a,tol)*2.d0
! test on density
if (rho >= tol) then
! call srLDA Ex[2*rho_a,2*rho_a]
call ex_lda_sr(mu,rho_a,rho_a,exerflda,vxerflda_a,vxerflda_b)
dexerfldadrho = (vxerflda_a + vxerflda_b)*0.5d0
! square of two times spin alpha density gradient
drho2=max(grd_rho_a_2,0d0)*4.0d0
kappa=0.804d0
sq=drho2*2.6121172985233599567768d-2*rho**(-8d0/3d0)
fx=1d0+kappa-kappa/(1d0+berf(1.616204596739954813d-1*mu*rho**(-f13))*sq/kappa)
exerfpbe_a=exerflda*fx
! Derivatives
sqs=-8d0*sq/(3d0*rho)
fxs=kappa**2*(-1.616204596739954813d-1*mu*rho**(-4d0*f13)/3d0*dberfda(1.616204596739954813d-1*mu*rho**(-f13))*sq+berf(1.616204596739954813d-1*mu*rho**(-f13))*sqs)/(kappa+berf(1.616204596739954813d-1*mu*rho**(-f13))*sq)**2
dexerfpbedrho_a=dexerfldadrho*fx+exerflda*fxs
sqss=2.6121172985233599567768d-2*rho**(-8d0/3d0)
dexerfpbeddrho2_a=exerflda*berf(1.616204596739954813d-1*mu*rho**(-1.d0/3.d0))*sqss*kappa**2/(kappa+berf(1.616204596739954813d-1*mu*rho**(-1.d0/3.d0))*sq)**2
endif
! two times spin beta density
rho = max(rho_b,tol)*2.d0
! test on density
if (rho >= tol) then
! call srLDA Ex[2*rho_b,2*rho_b]
call ex_lda_sr(mu,rho_b,rho_b,exerflda,vxerflda_a,vxerflda_b)
dexerfldadrho = (vxerflda_a + vxerflda_b)*0.5d0
! square of two times spin beta density gradient
drho2=max(grd_rho_b_2,0d0)*4.0d0
kappa=0.804d0
sq=drho2*2.6121172985233599567768d-2*rho**(-8d0/3d0)
fx=1d0+kappa-kappa/(1d0+berf(1.616204596739954813d-1*mu*rho**(-f13))*sq/kappa)
exerfpbe_b=exerflda*fx
! Derivatives
sqs=-8d0*sq/(3d0*rho)
fxs=kappa**2*(-1.616204596739954813d-1*mu*rho**(-4d0*f13)/3d0*dberfda(1.616204596739954813d-1*mu*rho**(-f13))*sq+berf(1.616204596739954813d-1*mu*rho**(-f13))*sqs)/(kappa+berf(1.616204596739954813d-1*mu*rho**(-f13))*sq)**2
dexerfpbedrho_b=dexerfldadrho*fx+exerflda*fxs
sqss=2.6121172985233599567768d-2*rho**(-8d0/3d0)
dexerfpbeddrho2_b=exerflda*berf(1.616204596739954813d-1*mu*rho**(-1.d0/3.d0))*sqss*kappa**2/(kappa+berf(1.616204596739954813d-1*mu*rho**(-1.d0/3.d0))*sq)**2
endif
ex = (exerfpbe_a+exerfpbe_b)*0.5d0
vx_rho_a = dexerfpbedrho_a
vx_rho_b = dexerfpbedrho_a
vx_grd_rho_a_2 = 2.d0*dexerfpbeddrho2_a
vx_grd_rho_b_2 = 2.d0*dexerfpbeddrho2_b
vx_grd_rho_a_b = 0.d0
end
subroutine ex_pbe_sr_only(mu,rho_a,rho_b,grd_rho_a_2,grd_rho_b_2,grd_rho_a_b,ex)
BEGIN_DOC
!rho_a = density alpha
!rho_b = density beta
!grd_rho_a_2 = (gradient rho_a)^2
!grd_rho_b_2 = (gradient rho_b)^2
!grd_rho_a_b = (gradient rho_a).(gradient rho_b)
!ex = exchange energy density at point r
END_DOC
implicit none
! input
double precision, intent(in) :: mu,rho_a, rho_b
double precision, intent(in) :: grd_rho_a_2, grd_rho_b_2, grd_rho_a_b
! output
double precision, intent(out) :: ex
! function
double precision berf
! local
double precision, parameter :: tol=1d-12
double precision, parameter :: f13=0.333333333333333d0
double precision exerflda,vxerflda_a,vxerflda_b
double precision exerfpbe_a, exerfpbe_b
double precision rho,drho2
double precision kappa,sq,fx
! initialization
ex=0.d0
! spin scaling relation Ex[rho_a,rho_b] = (1/2) (Ex[2rho_a,2rho_a] + Ex[2rho_b,2rho_b])
! two times spin alpha density
rho = max(rho_a,tol)*2.d0
! test on density
if (rho >= tol) then
! call srLDA Ex[2*rho_a,2*rho_a]
call ex_lda_sr(mu,rho_a,rho_a,exerflda,vxerflda_a,vxerflda_b)
! square of two times spin alpha density gradient
drho2=max(grd_rho_a_2,0d0)*4.0d0
kappa=0.804d0
sq=drho2*2.6121172985233599567768d-2*rho**(-8d0/3d0)
fx=1d0+kappa-kappa/(1d0+berf(1.616204596739954813d-1*mu*rho**(-f13))*sq/kappa)
exerfpbe_a=exerflda*fx
endif
! two times spin beta density
rho = max(rho_b,tol)*2.d0
! test on density
if (rho >= tol) then
! call srLDA Ex[2*rho_b,2*rho_b]
call ex_lda_sr(mu,rho_b,rho_b,exerflda,vxerflda_a,vxerflda_b)
! square of two times spin beta density gradient
drho2=max(grd_rho_b_2,0d0)*4.0d0
kappa=0.804d0
sq=drho2*2.6121172985233599567768d-2*rho**(-8d0/3d0)
fx=1d0+kappa-kappa/(1d0+berf(1.616204596739954813d-1*mu*rho**(-f13))*sq/kappa)
exerfpbe_b=exerflda*fx
endif
ex = (exerfpbe_a+exerfpbe_b)*0.5d0
end
subroutine ec_pbe_only(mu,rhoc,rhoo,sigmacc,sigmaco,sigmaoo,ec)
BEGIN_DOC
! Short-range PBE correlation energy functional for erf interaction
!
! input : ==========
!
! mu = range separated parameter
!
! rhoc, rhoo = total density and spin density
!
! sigmacc = square of the gradient of the total density
!
! sigmaco = square of the gradient of the spin density
!
! sigmaoo = scalar product between the gradient of the total density and the one of the spin density
!
! output: ==========
!
! ec = correlation energy
!
END_DOC
include 'constants.include.F'
implicit none
! input
double precision, intent(in) :: rhoc,rhoo,mu
double precision, intent(in) :: sigmacc,sigmaco,sigmaoo
! output
double precision, intent(out) :: ec
! local
double precision tol
parameter(tol=1d-12)
character(len=30) namedummy
double precision eccerflda
double precision vrhoccerflda
double precision vrhoocerflda
double precision ecclda
double precision vrhocclda
double precision vrhooclda
integer i,igrad
double precision rho,drho2,rhoa,rhob
double precision ecerflda
double precision eclda,decldadrho
double precision ecerfpbe
double precision arglog,alpha,beta,gamma
double precision Aa,Ab,Ac,tq
double precision zeta,phi,phi2,phi3,phi4
double precision, parameter :: f13=0.333333333333333d0
! Parameter of the modified interaction
ec = 0.d0
! First-type gradient functional
igrad=1
alpha=2.78d0
gamma=3.1091d-2
! test on density
if (dabs(rhoc).lt.tol) return
double precision :: vc_a,vc_b
! Spin polarisation
rhoa=max((rhoc+rhoo)*.5d0,1.0d-15)
rhob=max((rhoc-rhoo)*.5d0,1.0d-15)
call ec_lda_sr(mu,rhoa,rhob,eccerflda,vc_a,vc_b)
ecerflda = eccerflda
vrhoccerflda = 0.5d0 * (vc_a + vc_b)
vrhoocerflda = 0.5d0 * (vc_a - vc_b)
! Density
rho = rhoc
rho = max(rho,1.d-10)
! Square of density gradient
drho2 = sigmacc
zeta = (rhoa-rhob)/(rhoa+rhob)
zeta = max(zeta,1.d-10)
! LDA energy density
double precision :: vc_a_lda,vc_b_lda
call ec_lda(rhoa,rhob,ecclda,vc_a_lda,vc_b_lda)
eclda = ecclda
decldadrho = 0.5d0 * (vc_a_lda+vc_b_lda)
decldadrho = 0.5d0 * (vc_a_lda-vc_b_lda)
if ((ecerflda/eclda).le.0d0) then
beta=0d0
else
beta=6.6725d-2*(ecerflda/eclda)**alpha
endif
phi=((1d0+zeta)**(2d0/3d0)+(1d0-zeta)**(2d0/3d0))/2d0
phi2=phi*phi
phi3=phi2*phi
phi4=phi3*phi
tq=drho2*6.346820607d-2*rho**(-7d0/3d0)/phi2
Ab=dexp(-ecerflda/(rho*gamma*phi3))-1d0
if (dabs(Ab).le.dabs(beta*tol)) then
ecerfpbe=ecerflda
else
Aa=beta/(gamma*Ab)
Ac=1d0+Aa*tq+Aa**2*tq**2
if (Aa.lt.tol) Aa=tol
arglog=1d0+beta*(1d0-1d0/Ac)/(gamma*Aa)
arglog=max(arglog,1.d-10)
ecerfpbe=ecerflda+rho*phi3*gamma*dlog(arglog)
end if
ec = ecerfpbe
end