MODULE excor USE GlobalMath IMPLICIT NONE ! CHARACTER(132) :: exctype INTEGER, PRIVATE, SAVE :: itype INTEGER, PRIVATE, PARAMETER :: LDA_PW=14 INTEGER, PRIVATE, PARAMETER :: GGA_PBE=16 ! Parameters for the Perdew-Wang (PRB 45,13244 (1992)) LDA correlation REAL(8), PRIVATE, PARAMETER :: AA=0.0310907d0 REAL(8), PRIVATE, PARAMETER :: a1=0.21370d0 REAL(8), PRIVATE, PARAMETER :: b1=7.59570d0 REAL(8), PRIVATE, PARAMETER :: b2=3.58760d0 REAL(8), PRIVATE, PARAMETER :: b3=1.63820d0 REAL(8), PRIVATE, PARAMETER :: b4=0.49294d0 ! Parameters from PBE paper (and in program from web) REAL(8), PRIVATE, PARAMETER :: kappa= 0.804d0 REAL(8), PRIVATE, PARAMETER :: mu = 0.2195149727645171d0 REAL(8), PRIVATE, PARAMETER :: beta = 0.06672455060314922d0 REAL(8), PRIVATE, PARAMETER :: gamm = 0.03109069086965489503494086371273d0 REAL(8), PRIVATE, PARAMETER :: betabygamm=beta/gamm CONTAINS !*********************************************************** !* Initialization !*********************************************************** SUBROUTINE initexch ! choose form of exchange-correlation potential SELECT CASE(TRIM(exctype)) CASE default itype = LDA_PW exctype='LDA-PW' WRITE(6,*) 'Perdew-Wang correlation' CASE('LDA-PW') itype = LDA_PW WRITE(6,*) 'Perdew-Wang correlation' CASE('GGA-PBE') itype = GGA_PBE WRITE(6,*) 'Perdew - Burke - Ernzerhof GGA' END SELECT END SUBROUTINE initexch !********************************************************** ! Subroutine to calculate the LDA exchange correlation functionals ! using the form of Perdew and Wang (PRB 45, 13244 (1992) ! assuming no spin polarization ! Inside this routine, energies are in Hartree units SUBROUTINE pwldafunc(den,exc,vxc) REAL(8), INTENT(IN) :: den !density REAL(8), INTENT(OUT) :: exc,vxc ! Variables depending on den REAL(8) :: n,g,kf,rs,ks REAL(8) :: ex,ec,pprs,decdrs REAL(8) :: term n=den IF (n < machine_zero) THEN exc=0.d0; vxc=0.d0 RETURN ENDIF kf=(3.d0*(pi**2)*n)**0.3333333333333333333333333333d0 rs=(3.d0/(4.d0*pi*n))**0.3333333333333333333333333333d0 ks=SQRT(4.d0*kf/pi) ex=-3.d0*kf/(4.d0*pi) pprs=SQRT(rs)*(b1+b3*rs)+rs*(b2+b4*rs) !ec=-2.d0*AA*(1.d0+a1*rs)*ddlog(1.d0+1.d0/(2.d0*AA*pprs)) term=Logofterm(1.d0/(2.d0*AA*pprs)) ec=-2.d0*AA*(1.d0+a1*rs)*term exc=ex+ec !decdrs=-(2.d0*AA*a1)*ddlog(1.d0+1.d0/(2*AA*pprs)) & decdrs=-(2.d0*AA*a1)*term & +((1.d0+a1*rs)*((b1+3*b3*rs)/(2.d0*SQRT(rs))+b2+2*b4*rs))/(pprs*(pprs+1.d0/(2.d0*AA))) vxc = (4.d0/3.d0)*ex+ec-(decdrs*rs)/3.d0 IF ((ABS(exc).GT.1.d65).OR.(ABS(vxc).GT.1.d65)) THEN WRITE(6,*) 'Problem in PW',n,rs,ec ENDIF exc=2*exc; vxc=2*vxc ! change to Rydberg units RETURN END SUBROUTINE pwldafunc !********************************************************** ! Subroutine to calculate the exchange correlation functionals ! using the form of Perdew, Burke, and Ernzerhof (PRL 77, 3865 (1996)) ! assuming no spin polarization ! Inside this routine, energies are in Hartree units SUBROUTINE pbefunc(den,grad,fxc,dfxcdn,dfxcdgbg) REAL(8), INTENT(IN) :: den,grad !density, magnitude of grad(density) REAL(8), INTENT(OUT) :: fxc,dfxcdn,dfxcdgbg ! Variables depending on den & grad REAL(8) :: n,g,kf,rs,ks,s,t REAL(8) :: ex,ec,Fx,H,A,pprs,ppt,At2,dFds,dHdt,decdrs,dHdrs,dHdA,dAdrs REAL(8) :: term,dHdtbg,dFdsbg n=den IF (n < machine_zero) THEN fxc=0.d0; dfxcdn=0.d0; dfxcdgbg=0.d0 RETURN ENDIF g=grad IF (g < machine_zero) g=machine_zero kf=(3.d0*(pi**2)*n)**0.3333333333333333333333333333d0 rs=(3.d0/(4.d0*pi*n))**0.3333333333333333333333333333d0 ks=SQRT(4.d0*kf/pi) s=g/(2.d0*kf*n) t=g/(2.d0*ks*n) ex=-3.d0*kf/(4.d0*pi) pprs=SQRT(rs)*(b1+b3*rs)+rs*(b2+b4*rs) !ec=-2.d0*AA*(1.d0+a1*rs)*ddlog(1.d0+1.d0/(2.d0*AA*pprs)) term=Logofterm(1.d0/(2.d0*AA*pprs)) ec=-2.d0*AA*(1.d0+a1*rs)*term Fx=1.d0+kappa -kappa/(1.d0+(mu/kappa)*s*s) A=Aofec(ec) At2=A*t*t ppt=(1.d0+At2*(1.d0+At2)) !H=gamm*ddlog(1.d0+(betabygamm)*(t*t)*((1.d0+At2)/ppt)) H=gamm*Logofterm((betabygamm)*(t*t)*((1.d0+At2)/ppt)) fxc=n*(ex*Fx+ec+H) dFds = (2.d0*mu*s)/(1.d0+(mu/kappa)*(s**2))**2 dFdsbg = ((2.d0*mu)/(1.d0+(mu/kappa)*(s**2))**2)/(2.d0*kf*n) dHdt = (2.d0*t*beta*gamm*(1.d0+2.d0*At2))/((gamm*ppt+beta*t*t*(1.d0+At2))*ppt) dHdtbg = ((2.d0*beta*gamm*(1.d0+ & 2.d0*At2))/((gamm*ppt+beta*t*t*(1.d0+At2))*ppt))/(2.d0*ks*n) !decdrs=-(2.d0*AA*a1)*ddlog(1.d0+1.d0/(2*AA*pprs)) & decdrs=-(2.d0*AA*a1)*term & +((1.d0+a1*rs)*((b1+3*b3*rs)/(2.d0*SQRT(rs))+b2+2*b4*rs))/(pprs*(pprs+1.d0/(2.d0*AA))) dHdA=((2.d0+At2)*(At2*t*t*t*t*beta*gamm))/((gamm*ppt+beta*t*t*(1.d0+At2))*ppt) dAdrs=-ddexp(-ec/gamm)*A*A*decdrs/beta dHdrs=dHdA*dAdrs dfxcdn = (4.d0/3.d0)*ex*(Fx-dFds*s)+ec-(decdrs*rs)/3.d0+H-(dHdrs*rs)/3.d0 & - (7.d0/6.d0)*dHdt*t dfxcdgbg = ex*dFdsbg/(2.d0*kf) + dHdtbg/(2.d0*ks) IF ((ABS(fxc).GT.1.d65).OR.(ABS(dfxcdn).GT.1.d65).OR.(ABS(dfxcdgbg).GT.1.d65)) THEN WRITE(6,*) 'Problem in PBE',n,g,rs,s,t,ec,A,H ENDIF RETURN END SUBROUTINE pbefunc !******************************************************************* ! ! Function Logofterm -- needed to take care of behavior for small term ! Evaluates log(1+term) FUNCTION Logofterm(term) REAL(8) :: term, Logofterm IF (ABS(term)>machine_precision) THEN Logofterm=ddlog(1.d0+term) ELSE Logofterm=term ENDIF RETURN END FUNCTION Logofterm !******************************************************************* ! ! Function Aofec -- needed to take care of behavior for small ec FUNCTION Aofec(ec) REAL(8) :: ec, Aofec IF (ABS(ec)>machine_precision) THEN Aofec=betabygamm/(ddexp(-ec/gamm)-1.d0) ELSEIF (ABS(ec)>machine_zero) THEN Aofec=beta/(-ec) ELSE Aofec=-beta*DSIGN(machine_infinity,ec) ENDIF RETURN END FUNCTION Aofec !******************************************************************** ! ! Subroutine radialexcpbe ! Density(:) input on a uniform radial mesh of Npts ! R(:) input mesh points ! Exc - output integrated exchange correlation energy -- in Rydberg units ! vxc(:) -- output exchange correlation potential -- iin Rydberg units ! !******************************************************************** SUBROUTINE radialexcpbe(Npts,r,density,Exc,vxc) IMPLICIT NONE INTEGER, INTENT(IN) ::Npts REAL(8), INTENT(IN) :: r(Npts),density(Npts) REAL(8), INTENT(OUT) :: Exc, vxc(Npts) INTEGER :: i,j , ierr REAL(8) :: gradient(Npts),gradmag(Npts),gxc(Npts),dgxcdr(Npts) REAL(8) :: fxc(Npts),dfxcdn,dfxcdgbg,sgn REAL(8) :: h h=r(2)-r(1) CALL nderiv(h,density,gradient,Npts,ierr) IF (ierr/=0) THEN WRITE(6,*) 'Error in nderiv ', ierr STOP ENDIF gradmag=ABS(gradient) DO i=1,Npts CALL pbefunc(density(i),gradmag(i),fxc(i),dfxcdn,dfxcdgbg) vxc(i)=dfxcdn gxc(i)=dfxcdgbg*gradient(i) ENDDO CALL nderiv(h,gxc,dgxcdr,Npts,ierr) IF (ierr/=0) THEN WRITE(6,*) 'Error in nderiv ', ierr STOP ENDIF DO i=1,Npts fxc(i)=2*fxc(i)*4*pi*(r(i)*r(i)) !2* changes from Hartree to Rydberg vxc(i)=2*vxc(i)-2*dgxcdr(i) ! units ENDDO Exc = overint(Npts-1,h,fxc(2)) RETURN END SUBROUTINE radialexcpbe !******************************************************************* SUBROUTINE exch(n,h,den,rvxc,etxc,eexc) ! calculate exchange correlation potentials and energys ! for density functional theory from electron density ! den(n) is electron density * (4*pi*r**2) ! rvxc(n) is returned as vxc * r ! to the total energy ! eexc is the total exchange energy (int(den*exc)) ! etxc is eexc - int(den*vxc) INTEGER, INTENT(IN) :: n REAL(8), INTENT(IN) :: h,den(n) REAL(8), INTENT(INOUT) :: rvxc(n),etxc,eexc REAL(8), ALLOCATABLE :: tmpr(:),tmpd(:),tmpv(:) REAL(8) :: fpi INTEGER :: i REAL(8) :: r,r2,rho,exc,vxc fpi=4*pi rvxc=0 SELECT CASE(itype) CASE default WRITE(6,*) 'Error in exch -- itype not correct', itype STOP CASE(GGA_PBE) !!!!!!!PBE form!!!!!! ALLOCATE(tmpr(n+1),tmpd(n+1),tmpv(n+1)) DO i=1,n tmpr(i+1)=i*h tmpd(i+1)=den(i)/(fpi*(i*h*i*h)) ENDDO tmpr(1)=0 tmpd(1)=2*tmpd(2)-tmpd(3) CALL radialexcpbe(n+1,tmpr,tmpd,eexc,tmpv) WRITE(6,*) 'eexc',eexc DO i=1,n rvxc(i)=tmpv(i+1)*tmpr(i+1) tmpv(i)=tmpv(i+1)*den(i) ENDDO etxc=eexc-overint(n,h,tmpv) DEALLOCATE(tmpr,tmpd,tmpv) CASE (LDA_PW) !!! ! Perdew-Wang LDA !!!! ALLOCATE(tmpd(n),tmpv(n)) DO i=1,n r=i*h r2=r*r rho=den(i)/(fpi*r2) CALL pwldafunc(rho,exc,vxc) tmpd(i)=den(i)*(exc-vxc) tmpv(i)=den(i)*exc rvxc(i)=r*vxc ENDDO ! calculate exchange-correlation contribution to the potential ! etxc=dsum(n,a,1)*h ! eexc=dsum(n,b,1)*h etxc=overint(n,h,tmpd) eexc=overint(n,h,tmpv) WRITE(6,*) 'etxc,eexc = ',etxc,eexc DEALLOCATE(tmpd,tmpv) END SELECT END SUBROUTINE exch END MODULE excor