src/species/xc_pbe.f90

   1
   2 ! This routine is based on code written by K. Burke.
   3
   4 !BOP
   5 ! !ROUTINE: xc_pbe
   6 ! !INTERFACE:
   7 subroutine xc_pbe(n,kappa,mu,beta,rhoup,rhodn,grho,gup,gdn,g2up,g2dn,g3rho, &
   8  g3up,g3dn,ex,ec,vxup,vxdn,vcup,vcdn)
   9 ! !INPUT/OUTPUT PARAMETERS:
  10 !   n     : number of density points (in,integer)
  11 !   kappa : parameter for large-gradient limit (in,real)
  12 !   mu    : gradient expansion coefficient (in,real)
  13 !   beta  : gradient expansion coefficient (in,real)
  14 !   rhoup : spin-up charge density (in,real(n))
  15 !   rhodn : spin-down charge density (in,real(n))
  16 !   grho  : |grad rho| (in,real(n))
  17 !   gup   : |grad rhoup| (in,real(n))
  18 !   gdn   : |grad rhodn| (in,real(n))
  19 !   g2up  : grad^2 rhoup (in,real(n))
  20 !   g2dn  : grad^2 rhodn (in,real(n))
  21 !   g3rho : (grad rho).(grad |grad rho|) (in,real(n))
  22 !   g3up  : (grad rhoup).(grad |grad rhoup|) (in,real(n))
  23 !   g3dn  : (grad rhodn).(grad |grad rhodn|) (in,real(n))
  24 !   ex    : exchange energy density (out,real(n))
  25 !   ec    : correlation energy density (out,real(n))
  26 !   vxup  : spin-up exchange potential (out,real(n))
  27 !   vxdn  : spin-down exchange potential (out,real(n))
  28 !   vcup  : spin-up correlation potential (out,real(n))
  29 !   vcdn  : spin-down correlation potential (out,real(n))
  30 ! !DESCRIPTION:
  31 !   Spin-polarised exchange-correlation potential and energy of the generalised
  32 !   gradient approximation functional of J. P. Perdew, K. Burke and M. Ernzerhof
  33 !   {\it Phys. Rev. Lett.} {\bf 77}, 3865 (1996) and {\bf 78}, 1396(E) (1997).
  34 !   The parameter $\kappa$, which controls the large-gradient limit, can be set
  35 !   to $0.804$ or $1.245$ corresponding to the value in the original article or
  36 !   the revised version of Y. Zhang and W. Yang, {\it Phys. Rev. Lett.}
  37 !   {\bf 80}, 890 (1998).
  38 !
  39 ! !REVISION HISTORY:
  40 !   Modified routines written by K. Burke, October 2004 (JKD)
  41 !EOP
  42 !BOC
  43 implicit none
  44 ! arguments
  45 integer, intent(in) :: n
  46 real(8), intent(in) :: kappa
  47 real(8), intent(in) :: mu
  48 real(8), intent(in) :: beta
  49 real(8), intent(in) :: rhoup(n)
  50 real(8), intent(in) :: rhodn(n)
  51 real(8), intent(in) :: grho(n)
  52 real(8), intent(in) :: gup(n)
  53 real(8), intent(in) :: gdn(n)
  54 real(8), intent(in) :: g2up(n)
  55 real(8), intent(in) :: g2dn(n)
  56 real(8), intent(in) :: g3rho(n)
  57 real(8), intent(in) :: g3up(n)
  58 real(8), intent(in) :: g3dn(n)
  59 real(8), intent(out) :: ex(n)
  60 real(8), intent(out) :: ec(n)
  61 real(8), intent(out) :: vxup(n)
  62 real(8), intent(out) :: vxdn(n)
  63 real(8), intent(out) :: vcup(n)
  64 real(8), intent(out) :: vcdn(n)
  65 ! local variables
  66 integer i
  67 real(8), parameter :: thrd=1.d0/3.d0
  68 real(8), parameter :: thrd2=2.d0/3.d0
  69 real(8), parameter :: pi=3.1415926535897932385d0
  70 ! maximum allowed |grad rho|
  71 real(8), parameter :: gmax=1.d6
  72 ! maximum allowed grad^2 rho
  73 real(8), parameter :: g2max=1.d12
  74 ! maximum allowed (grad rho).(grad |grad rho|)
  75 real(8), parameter :: g3max=1.d14
  76 real(8) r,r2,kf,s,u,v
  77 real(8) rs,z,g,ks,ksg
  78 real(8) t,uu,vv,ww
  79 real(8) grho_,gup_,gdn_,g2up_,g2dn_,g2rho_
  80 real(8) g3rho_,g3up_,g3dn_
  81 real(8) exup,exdn
  82 do i=1,n
  83   if ((rhoup(i).gt.1.d-12).and.(rhodn(i).gt.1.d-12)) then
  84     grho_=grho(i)
  85     gup_=gup(i)
  86     gdn_=gdn(i)
  87     g2up_=g2up(i)
  88     g2dn_=g2dn(i)
  89     g3rho_=g3rho(i)
  90     g3up_=g3up(i)
  91     g3dn_=g3dn(i)
  92 ! check gradients are within range
  93     if (grho_.gt.gmax) grho_=gmax
  94     if (gup_.gt.gmax) gup_=gmax
  95     if (gdn_.gt.gmax) gdn_=gmax
  96     if (abs(g2up_).gt.g2max) g2up_=sign(g2max,g2up_)
  97     if (abs(g2dn_).gt.g2max) g2dn_=sign(g2max,g2dn_)
  98     if (abs(g3rho_).gt.g3max) g3rho_=sign(g3max,g3rho_)
  99     if (abs(g3up_).gt.g3max) g3up_=sign(g3max,g3up_)
 100     if (abs(g3dn_).gt.g3max) g3dn_=sign(g3max,g3dn_)
 101 ! exchange energy density and potential
 102 ! spin-up
 103     r=rhoup(i)
 104     r2=2.d0*r
 105     kf=(r2*3.d0*pi**2)**thrd
 106     s=gup_/(2.d0*kf*r)
 107     u=g3up_/((r**2)*(2.d0*kf)**3)
 108     v=g2up_/(r*(2.d0*kf)**2)
 109     call x_pbe(kappa,mu,r2,s,u,v,exup,vxup(i))
 110 ! spin-down
 111     r=rhodn(i)
 112     r2=2.d0*r
 113     kf=(r2*3.d0*pi**2)**thrd
 114     s=gdn_/(2.d0*kf*r)
 115     u=g3dn_/((r**2)*(2.d0*kf)**3)
 116     v=g2dn_/(r*(2.d0*kf)**2)
 117     call x_pbe(kappa,mu,r2,s,u,v,exdn,vxdn(i))
 118 ! total density
 119     r=rhoup(i)+rhodn(i)
 120 ! average exchange energy density
 121     ex(i)=(exup*rhoup(i)+exdn*rhodn(i))/r
 122 ! correlation
 123     rs=(3.d0/(4.d0*pi*r))**thrd
 124     z=(rhoup(i)-rhodn(i))/r
 125     g=((1.d0+z)**thrd2+(1.d0-z)**thrd2)/2.d0
 126     kf=(r*3.d0*pi**2)**thrd
 127     ks=sqrt(4.d0*kf/pi)
 128     ksg=2.d0*ks*g
 129     t=grho_/(ksg*r)
 130     uu=g3rho_/((r**2)*ksg**3)
 131     g2rho_=g2up_+g2dn_
 132     vv=g2rho_/(r*ksg**2)
 133     ww=(gup_**2-gdn_**2-z*grho_**2)/(r*r*ksg**2)
 134     call c_pbe(beta,rs,z,t,uu,vv,ww,ec(i),vcup(i),vcdn(i))
 135   else
 136     ex(i)=0.d0
 137     ec(i)=0.d0
 138     vxup(i)=0.d0
 139     vxdn(i)=0.d0
 140     vcup(i)=0.d0
 141     vcdn(i)=0.d0
 142   end if
 143 end do
 144 return
 145 end subroutine
 146 !EOC
 147