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