TABLE OF CONTENTS


ABINIT/xcpzca [ Functions ]

[ Top ] [ Functions ]

NAME

 xcpzca

FUNCTION

 Returns exc, vxc, and d(vxc)/d($\rho$) from input rho.

 NOTE
 Perdew-Zunger parameterization of Ceperly-Alder electron gas
 energy data--
 J. Perdew and A. Zunger, Phys. Rev. B 23, 5048 (1981).
 D.M. Ceperly and B.J. Alder, Phys. Rev. Lett. 45, 566 (1980).

COPYRIGHT

 Copyright (C) 1998-2018 ABINIT group (DCA, XG, GMR)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .
 For the initials of contributors, see ~abinit/doc/developers/contributors.txt .

INPUTS

  npt=number of real space points on which density is provided
  order=gives the maximal derivative of Exc computed.
  rhor(npt)=electron number density (bohr^-3)
  rspts(npt)=corresponding Wigner-Seitz radii, precomputed

OUTPUT

  exc(npt)=exchange-correlation energy density (hartree)
  vxc(npt)=xc potential (d($\rho$*exc)/d($\rho$)) (hartree)
  if(order>1) dvxc(npt)=derivative d(vxc)/d($\rho$) (hartree*bohr^3)

PARENTS

      drivexc

CHILDREN

SOURCE

 40 #if defined HAVE_CONFIG_H
 41 #include "config.h"
 42 #endif
 43 
 44 #include "abi_common.h"
 45 
 46 
 47 subroutine xcpzca(exc,npt,order,rhor,rspts,vxc,&  !Mandatory arguments
 48 &                dvxc)                            !Optional arguments
 49 
 50  use defs_basis
 51  use m_errors
 52  use m_profiling_abi
 53 
 54 !This section has been created automatically by the script Abilint (TD).
 55 !Do not modify the following lines by hand.
 56 #undef ABI_FUNC
 57 #define ABI_FUNC 'xcpzca'
 58 !End of the abilint section
 59 
 60  implicit none
 61 
 62 !Arguments ------------------------------------
 63 !scalars
 64  integer,intent(in) :: npt,order
 65 !arrays
 66  real(dp),intent(in) :: rhor(npt),rspts(npt)
 67  real(dp),intent(out) :: exc(npt),vxc(npt)
 68  real(dp),intent(out),optional :: dvxc(npt)
 69 
 70 !Local variables-------------------------------
 71 !Perdew-Zunger parameters a, b, b1, b2, c, d, gamma
 72 !scalars
 73  integer :: ipt
 74  real(dp),parameter :: aa=0.0311_dp,b1=1.0529_dp,b2=0.3334_dp,bb=-0.048_dp
 75  real(dp),parameter :: c4_3=4.0_dp/3.0_dp,c7_6=7.0_dp/6.0_dp,cc=0.0020_dp
 76  real(dp),parameter :: dd=-0.0116_dp,ga=-0.1423_dp
 77  real(dp) :: den,den3,dfac,efac,logrs,rs,rsm1,t1,t2,vfac
 78  character(len=500) :: message
 79 
 80 ! *************************************************************************
 81 
 82 !Compute vfac=(3/(2*Pi))^(2/3)
 83  vfac=(1.5_dp/pi)**(2.0_dp/3.0_dp)
 84 !Compute efac=(3/4)*vfac
 85  efac=0.75_dp*vfac
 86 !Compute dfac=(4*Pi/9)*vfac
 87  dfac=(4.0_dp*pi/9.0_dp)*vfac
 88 
 89 !Checks the values of order
 90  if(order<0 .or. order>2)then
 91    write(message, '(a,a,a,i0)' )&
 92 &   'With Perdew-Zunger Ceperley-Alder xc functional, the only',ch10,&
 93 &   'allowed values for order are 0, 1 or 2, while it is found to be',order
 94    MSG_BUG(message)
 95  end if
 96 
 97 !Checks the compatibility between the order and the presence of the optional arguments
 98  if(order <= 1 .and. present(dvxc))then
 99    write(message, '(a,a,a,i0)' )&
100 &   'The order chosen does not need the presence',ch10,&
101 &   'of the vector dvxc, that is needed only with order=2 , while we have',order
102    MSG_BUG(message)
103  end if
104 
105 !separate cases with respect to order
106  if(order==2) then
107 !  Loop over grid points
108    do ipt=1,npt
109      rs=rspts(ipt)
110      rsm1=1.0_dp/rs
111 !    Consider two regimes: rs<1 or rs>=1
112      if (rs<1._dp) then
113        logrs=log(rs)
114 !      compute energy density exc (hartree)
115        exc(ipt)=(aa+cc*rs)*logrs+dd*rs+bb-efac*rsm1
116 !      compute potential vxc=d(rho*exc)/d(rho) (hartree)
117        vxc(ipt)=(aa+two_thirds*cc*rs)*logrs+(dd+dd-cc)*rs*third+&
118 &       (bb-aa*third)-vfac*rsm1
119 !      compute d(vxc)/d(rho) (hartree*bohr^3)
120        dvxc(ipt)=-(3._dp*aa+(cc+dd+dd)*rs+2._dp*cc*rs*logrs)&
121 &       /(9._dp*rhor(ipt))-dfac*rs**2
122      else if (rs<1000._dp) then
123        t1=b1*sqrt(rs)
124        t2=b2*rs
125        den=1._dp/(1._dp+t1+t2)
126        exc(ipt)=ga*den-efac*rsm1
127        vxc(ipt)=ga*(1._dp+c7_6*t1+c4_3*t2)*den**2-vfac*rsm1
128        den3=den**3
129        dvxc(ipt)=(ga*den3/(36._dp*rhor(ipt)))*(5._dp*t1+8._dp*t2+&
130 &       7._dp*t1**2+16._dp*t2**2+21._dp*t1*t2)-dfac*rs**2
131      else
132        t1=b1*sqrt(rs)
133        t2=b2*rs
134        den=1._dp/(1._dp+t1+t2)
135        exc(ipt)=ga*den-efac*rsm1
136        vxc(ipt)=ga*(1._dp+c7_6*t1+c4_3*t2)*den**2-vfac*rsm1
137        dvxc(ipt)=0._dp
138      end if
139    end do
140  else
141 !  Loop over grid points
142    do ipt=1,npt
143      rs=rspts(ipt)
144      rsm1=1.0_dp/rs
145 !    Consider two regimes: rs<1 or rs>=1
146      if (rs<1._dp) then
147        logrs=log(rs)
148 !      compute energy density exc (hartree)
149        exc(ipt)=(aa+cc*rs)*logrs+dd*rs+bb-efac*rsm1
150 !      compute potential vxc=d(rho*exc)/d(rho) (hartree)
151        vxc(ipt)=(aa+two_thirds*cc*rs)*logrs+(dd+dd-cc)*rs*third+&
152 &       (bb-aa*third)-vfac*rsm1
153 !      compute d(vxc)/d(rho) (hartree*bohr^3)
154      else
155        t1=b1*sqrt(rs)
156        t2=b2*rs
157        den=1._dp/(1._dp+t1+t2)
158        exc(ipt)=ga*den-efac*rsm1
159        vxc(ipt)=ga*(1._dp+c7_6*t1+c4_3*t2)*den**2-vfac*rsm1
160      end if
161    end do
162  end if
163 !
164 end subroutine xcpzca