TABLE OF CONTENTS
ABINIT/xcpzca [ 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