TABLE OF CONTENTS
- ABINIT/m_xclda
- ABINIT/xchelu
- ABINIT/xclb
- ABINIT/xcpzca
- ABINIT/xcspol
- ABINIT/xctetr
- ABINIT/xctfw
- ABINIT/xcwign
- ABINIT/xcxalp
ABINIT/m_xclda [ Modules ]
NAME
m_xclda
FUNCTION
LDA or LDA-like XC functionals.
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (DCA, XG, GMR, LG, MF, JFD, LK) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
PARENTS
CHILDREN
SOURCE
21 #if defined HAVE_CONFIG_H 22 #include "config.h" 23 #endif 24 25 #include "abi_common.h" 26 27 module m_xclda 28 29 use defs_basis 30 use m_errors 31 use m_abicore 32 33 use m_numeric_tools, only : invcb 34 35 implicit none 36 37 private
ABINIT/xchelu [ Functions ]
NAME
xchelu
FUNCTION
Returns exc, vxc, and eventually d(vxc)/d($\rho$) from input rho.
NOTES
Hedin-Lundqvist exchange and correlation (xc)-- L. Hedin and B.I. Lundqvist, J. Phys. C. 4, 2064 (1971) [[cite:Hedin1971]]
INPUTS
npt=number of real space points on which density is provided order=gives the maximal derivative of Exc computed. rspts(npt)=Wigner-Seitz radii at each point
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
1054 subroutine xchelu(exc,npt,order,rspts,vxc,dvxc) ! dvxc is optional 1055 1056 1057 !This section has been created automatically by the script Abilint (TD). 1058 !Do not modify the following lines by hand. 1059 #undef ABI_FUNC 1060 #define ABI_FUNC 'xchelu' 1061 !End of the abilint section 1062 1063 implicit none 1064 1065 !Arguments ------------------------------------ 1066 !scalars 1067 integer,intent(in) :: npt,order 1068 !arrays 1069 real(dp),intent(in) :: rspts(npt) 1070 real(dp),intent(out) :: exc(npt),vxc(npt) 1071 real(dp),intent(out),optional :: dvxc(npt) 1072 1073 !Local variables------------------------------- 1074 !aa and cc are H-L fitting parameters A and C (C in hartree) 1075 !rs = (3/(4 Pi))**(1/3) * rho(r)**(-1/3). 1076 !scalars 1077 integer :: ipt 1078 real(dp),parameter :: aa=21_dp,c1_21=one/21_dp,c4_9=4.0_dp/9.0_dp,cc=0.0225_dp 1079 real(dp) :: dfac,efac,rs,rsm1,vfac,xx 1080 character(len=500) :: message 1081 1082 ! ************************************************************************* 1083 1084 !Checks the values of order 1085 if(order<0 .or. order>2)then 1086 write(message, '(a,a,a,i0)' )& 1087 & 'With Hedin-Lundqvist xc functional, the only',ch10,& 1088 & 'allowed values for order are 0, 1 or 2, while it is found to be',order 1089 MSG_BUG(message) 1090 end if 1091 1092 !Compute vfac=(3/(2*Pi))^(2/3) 1093 vfac=(1.5_dp/pi)**(2.0_dp/3.0_dp) 1094 !Compute efac=(3/4)*vfac 1095 efac=0.75_dp*vfac 1096 !Compute dfac=(4*Pi/9)*vfac 1097 dfac=(4.0_dp*pi/9.0_dp)*vfac 1098 !separate cases with respect to order 1099 if (order==2) then 1100 ! Loop over grid points 1101 do ipt=1,npt 1102 rs=rspts(ipt) 1103 rsm1=one/rs 1104 ! compute energy density exc (hartree) 1105 xx=rs*c1_21 1106 exc(ipt)=-cc*((one+xx**3)*log(one+one/xx)+& 1107 & half*xx-xx*xx-third) - efac*rsm1 1108 ! compute xc potential d(rho*exc)/d(rho) (hartree) 1109 vxc(ipt)=-cc*log(one+aa*rsm1)-vfac*rsm1 1110 ! compute d(vxc)/d(rho) (hartree*bohr^3) 1111 dvxc(ipt)=-(rs**2)*((c4_9*pi)*cc*rs/(one+xx) + dfac) 1112 end do 1113 else 1114 ! Loop over grid points 1115 do ipt=1,npt 1116 rs=rspts(ipt) 1117 rsm1=one/rs 1118 ! compute energy density exc (hartree) 1119 xx=rs*c1_21 1120 exc(ipt)=-cc*((one+xx**3)*log(one+one/xx)+& 1121 & half*xx-xx*xx-third) - efac*rsm1 1122 ! compute xc potential d(rho*exc)/d(rho) (hartree) 1123 vxc(ipt)=-cc*log(one+aa*rsm1)-vfac*rsm1 1124 end do 1125 end if 1126 ! 1127 end subroutine xchelu
ABINIT/xclb [ Functions ]
NAME
xclb
FUNCTION
Computes the GGA like part (vx_lb) of the Leeuwen-Baerends exchange-correlation potential (vxc_lb) and adds it to the lda exchange-correlation potential (vxc_lda) which must be provided as input, vxci <-- vxc_lb =: vxc_lda + vx_lb R van Leeuwen and EJ Baerends, Phys Rev A 49, 2421 (1994) [[cite:VanLeeuwen1994]] With respect to spin, the van Leeuwen-Baerends potential is "exchange-like" : separate contributions from spin up and spin down.
INPUTS
npts= number of points to be computed
nspden=1 for unpolarized, 2 for spin-polarized
grho2_updn(npts,2*nspden-1)=square of the gradient of the spin-up,
and, if nspden==2, spin-down, and total density (Hartree/Bohr**2)
rho_updn(npts,nspden)=spin-up and spin-down density (Hartree/bohr**3)
OUTPUT
(see side effects)
SIDE EFFECTS
Input/Output: vxci(npts,nspden)=input xc potential to which Leeuwen-Baerends correction is added at output.
PARENTS
drivexc
CHILDREN
SOURCE
1269 subroutine xclb(grho2_updn,npts,nspden,rho_updn,vxci) 1270 1271 1272 !This section has been created automatically by the script Abilint (TD). 1273 !Do not modify the following lines by hand. 1274 #undef ABI_FUNC 1275 #define ABI_FUNC 'xclb' 1276 !End of the abilint section 1277 1278 implicit none 1279 1280 !Arguments ------------------------------------ 1281 !scalars 1282 integer,intent(in) :: npts,nspden 1283 !arrays 1284 real(dp),intent(in) :: grho2_updn(npts,2*nspden-1),rho_updn(npts,nspden) 1285 real(dp),intent(inout) :: vxci(npts,nspden) 1286 1287 !Local variables------------------------------- 1288 !scalars 1289 integer :: ipts,ispden 1290 real(dp),parameter :: beta=0.05_dp 1291 real(dp) :: density,density_gradient,density_t13,s_g_sq,scaled_gradient 1292 real(dp) :: scaling_factor,vx_lb 1293 1294 ! ************************************************************************* 1295 1296 !DEBUG 1297 !write(std_out,*) ' %xclb: enter' 1298 !ENDDEBUG 1299 1300 !scale the spin densities for evaluating spin up or down exchange 1301 scaling_factor=one 1302 if(nspden == 2) scaling_factor=two 1303 1304 do ispden=1,nspden 1305 1306 do ipts=1,npts 1307 1308 density= scaling_factor * rho_updn(ipts,ispden) 1309 density_gradient= scaling_factor * sqrt(grho2_updn(ipts,ispden)) 1310 1311 density_t13= density**third 1312 scaled_gradient= density_gradient/max(density*density_t13,1.e-12_dp) 1313 1314 s_g_sq= scaled_gradient*scaled_gradient 1315 1316 vx_lb= -beta*density_t13 * s_g_sq/ & 1317 & (one+3.d0*beta* scaled_gradient*log(scaled_gradient+sqrt(one+s_g_sq*s_g_sq))) 1318 1319 vxci(ipts,ispden)=vxci(ipts,ispden)+vx_lb 1320 end do 1321 1322 end do 1323 1324 end subroutine xclb
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) [[cite:Perdew1981]] D.M. Ceperley and B.J. Alder, Phys. Rev. Lett. 45, 566 (1980) [[cite:Ceperley1980]]
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
83 subroutine xcpzca(exc,npt,order,rhor,rspts,vxc,& !Mandatory arguments 84 & dvxc) !Optional arguments 85 86 87 !This section has been created automatically by the script Abilint (TD). 88 !Do not modify the following lines by hand. 89 #undef ABI_FUNC 90 #define ABI_FUNC 'xcpzca' 91 !End of the abilint section 92 93 implicit none 94 95 !Arguments ------------------------------------ 96 !scalars 97 integer,intent(in) :: npt,order 98 !arrays 99 real(dp),intent(in) :: rhor(npt),rspts(npt) 100 real(dp),intent(out) :: exc(npt),vxc(npt) 101 real(dp),intent(out),optional :: dvxc(npt) 102 103 !Local variables------------------------------- 104 !Perdew-Zunger parameters a, b, b1, b2, c, d, gamma 105 !scalars 106 integer :: ipt 107 real(dp),parameter :: aa=0.0311_dp,b1=1.0529_dp,b2=0.3334_dp,bb=-0.048_dp 108 real(dp),parameter :: c4_3=4.0_dp/3.0_dp,c7_6=7.0_dp/6.0_dp,cc=0.0020_dp 109 real(dp),parameter :: dd=-0.0116_dp,ga=-0.1423_dp 110 real(dp) :: den,den3,dfac,efac,logrs,rs,rsm1,t1,t2,vfac 111 character(len=500) :: message 112 113 ! ************************************************************************* 114 115 !Compute vfac=(3/(2*Pi))^(2/3) 116 vfac=(1.5_dp/pi)**(2.0_dp/3.0_dp) 117 !Compute efac=(3/4)*vfac 118 efac=0.75_dp*vfac 119 !Compute dfac=(4*Pi/9)*vfac 120 dfac=(4.0_dp*pi/9.0_dp)*vfac 121 122 !Checks the values of order 123 if(order<0 .or. order>2)then 124 write(message, '(a,a,a,i0)' )& 125 & 'With Perdew-Zunger Ceperley-Alder xc functional, the only',ch10,& 126 & 'allowed values for order are 0, 1 or 2, while it is found to be',order 127 MSG_BUG(message) 128 end if 129 130 !Checks the compatibility between the order and the presence of the optional arguments 131 if(order <= 1 .and. present(dvxc))then 132 write(message, '(a,a,a,i0)' )& 133 & 'The order chosen does not need the presence',ch10,& 134 & 'of the vector dvxc, that is needed only with order=2 , while we have',order 135 MSG_BUG(message) 136 end if 137 138 !separate cases with respect to order 139 if(order==2) then 140 ! Loop over grid points 141 do ipt=1,npt 142 rs=rspts(ipt) 143 rsm1=1.0_dp/rs 144 ! Consider two regimes: rs<1 or rs>=1 145 if (rs<1._dp) then 146 logrs=log(rs) 147 ! compute energy density exc (hartree) 148 exc(ipt)=(aa+cc*rs)*logrs+dd*rs+bb-efac*rsm1 149 ! compute potential vxc=d(rho*exc)/d(rho) (hartree) 150 vxc(ipt)=(aa+two_thirds*cc*rs)*logrs+(dd+dd-cc)*rs*third+& 151 & (bb-aa*third)-vfac*rsm1 152 ! compute d(vxc)/d(rho) (hartree*bohr^3) 153 dvxc(ipt)=-(3._dp*aa+(cc+dd+dd)*rs+2._dp*cc*rs*logrs)& 154 & /(9._dp*rhor(ipt))-dfac*rs**2 155 else if (rs<1000._dp) then 156 t1=b1*sqrt(rs) 157 t2=b2*rs 158 den=1._dp/(1._dp+t1+t2) 159 exc(ipt)=ga*den-efac*rsm1 160 vxc(ipt)=ga*(1._dp+c7_6*t1+c4_3*t2)*den**2-vfac*rsm1 161 den3=den**3 162 dvxc(ipt)=(ga*den3/(36._dp*rhor(ipt)))*(5._dp*t1+8._dp*t2+& 163 & 7._dp*t1**2+16._dp*t2**2+21._dp*t1*t2)-dfac*rs**2 164 else 165 t1=b1*sqrt(rs) 166 t2=b2*rs 167 den=1._dp/(1._dp+t1+t2) 168 exc(ipt)=ga*den-efac*rsm1 169 vxc(ipt)=ga*(1._dp+c7_6*t1+c4_3*t2)*den**2-vfac*rsm1 170 dvxc(ipt)=0._dp 171 end if 172 end do 173 else 174 ! Loop over grid points 175 do ipt=1,npt 176 rs=rspts(ipt) 177 rsm1=1.0_dp/rs 178 ! Consider two regimes: rs<1 or rs>=1 179 if (rs<1._dp) then 180 logrs=log(rs) 181 ! compute energy density exc (hartree) 182 exc(ipt)=(aa+cc*rs)*logrs+dd*rs+bb-efac*rsm1 183 ! compute potential vxc=d(rho*exc)/d(rho) (hartree) 184 vxc(ipt)=(aa+two_thirds*cc*rs)*logrs+(dd+dd-cc)*rs*third+& 185 & (bb-aa*third)-vfac*rsm1 186 ! compute d(vxc)/d(rho) (hartree*bohr^3) 187 else 188 t1=b1*sqrt(rs) 189 t2=b2*rs 190 den=1._dp/(1._dp+t1+t2) 191 exc(ipt)=ga*den-efac*rsm1 192 vxc(ipt)=ga*(1._dp+c7_6*t1+c4_3*t2)*den**2-vfac*rsm1 193 end if 194 end do 195 end if 196 ! 197 end subroutine xcpzca
ABINIT/xcspol [ Functions ]
NAME
xcspol
FUNCTION
Spin-polarized exchange and correlation, parameterized by Mike Teter of Corning Incorporated.
INPUTS
nspden=number of spin-density components npts= number of points to be computed order=its absolute value gives the maximal derivative of Exc to be computed. rspts(npts)=Seitz electron radius (bohr) zeta(npts)=$(\rho\uparrow-\rho\downarrow)/(\rho\uparrow+\rho\downarrow)$=degree of polarization (ignored if nspden=1, in which case zeta should be 0)
OUTPUT
if(abs(order)>1) dvxc(npts,1+nspden)= (Hartree*bohr^3)
if(nspden=1 .and. order==2): dvxc(:,1)=dvxc/d$\rho$ , dvxc(:,2) empty
if(nspden=1 .and. order==-2): also compute dvxc(:,2)=dvxc($\uparrow$)/d$\rho(\downarrow)$
if(nspden=2): dvxc(:,1)=dvxc($\uparrow$)/d$\rho(\uparrow)$,
dvxc(:,2)=dvxc($\uparrow$)/d$\rho(\downarrow)$, dvxc(:,3)=dvxc($\downarrow$)/d$\rho(\downarrow)$
exc(npts)=exchange-correlation energy density (hartree)
vxc(npts,nspden)=xc potent. (d($\rho$*exc)/d($\rho\uparrow$)) and d/d($\rho\downarrow$) (ha)
(only overall potential d($\rho$*exc)/d($\rho$) returned in vxc(1) for nspden=1)
ndvxc= size of dvxc(npts,ndvxc)
Normalization: Exc=$\int(exc(r)*\rho(r) d^3 r)$ for $\rho$(r)=electron density.
TODO
To be added later d2vxc=derivative $d^2 (Vxc)/d(rho)^2$ (hartree*bohr^6)
NOTES
This form is based on Mike Teter s rational polynomial exc=-(a0+a1*rs+a2*rs**2+a3*rs**3)/(b1*rs+b2*rs**2+b3*rs**3+b4*rs**4) where the parameters are fit to reproduce (in this case) the Perdew-Wang parameterization of the correlation energy given in Phys. Rev. B 45, 13244-13249 (1992) [[cite:Perdew1992]]. Each parameter is interpolated between zeta=0 and 1 by a_i(zeta)=a_i(0)+(a_i(1)-a_i(0))*f_x(zeta) and f_x(zeta)=[(1+zeta)$^{4/3}$+(1-zeta)$^{4/3}$-2]/(2*(2$^{1/3}$-1)). Beware : in this expression, zeta is actually replaced by zeta*alpha_zeta, where alpha_zeta is very close to 1, but slightly lower. This is to remove the singularity in the derivatives when abs(zeta) is 1 Below, a_i(1)-a_i(0) is called "da" for delta a, same for b s. rs = $(3/(4\pi))^{1/3} * \rho(r)^{-1/3}$ zeta = $(\rho\uparrow-\rho\downarrow)/(\rho\uparrow+\rho\downarrow)$ b1 must be 1 and a0 must be $(3/4)(3/(2\pi))^{2/3}$.
PARENTS
drivexc
CHILDREN
SOURCE
260 subroutine xcspol(exc,npts,nspden,order,rspts,vxc,zeta,ndvxc,& !Mandatory arguments 261 & dvxc) !Optional arguments 262 263 264 !This section has been created automatically by the script Abilint (TD). 265 !Do not modify the following lines by hand. 266 #undef ABI_FUNC 267 #define ABI_FUNC 'xcspol' 268 !End of the abilint section 269 270 implicit none 271 272 !Arguments ------------------------------------ 273 !scalars 274 integer,intent(in) :: ndvxc,npts,nspden,order 275 !arrays 276 real(dp),intent(in) :: rspts(npts),zeta(npts) 277 real(dp),intent(out) :: exc(npts),vxc(npts,nspden) 278 real(dp),intent(out),optional :: dvxc(npts,ndvxc) 279 280 !Local variables------------------------------- 281 !The generation of density from rs needs rsfac and rsfac^(-3) : 282 !rsfac=(3/(4 Pi))^(1/3) ; rsfacm3=4pi/3 283 !Mike Teter s parameters of 8 April 1993. 284 !New parameters which accomodate spin polarization (fit to P-W) 285 !Paramagnetic limit:a0p,...b4p 286 !(a0=(3/4)(3/(2 Pi))^(2/3)) (note that b1=1 is fixed) 287 !Differences, ferromagnetic - paramagnetic (delta params):da1,da2,da3,db1,db2,db3,db4 288 !scalars 289 integer :: ipts 290 real(dp),parameter :: a0p=.4581652932831429_dp,a1p=2.217058676663745_dp 291 real(dp),parameter :: a2p=0.7405551735357053_dp,a3p=0.01968227878617998_dp 292 real(dp),parameter :: alpha_zeta=one-1.0d-6,b1p=one,b2p=4.504130959426697_dp 293 real(dp),parameter :: b3p=1.110667363742916_dp,b4p=0.02359291751427506_dp 294 real(dp),parameter :: da0=.119086804055547_dp,da1=0.6157402568883345_dp 295 real(dp),parameter :: da2=0.1574201515892867_dp,da3=0.003532336663397157_dp 296 real(dp),parameter :: db1=zero,db2=0.2673612973836267_dp 297 real(dp),parameter :: db3=0.2052004607777787_dp,db4=0.004200005045691381_dp 298 real(dp),parameter :: ft=4._dp/3._dp,rsfac=0.6203504908994000_dp 299 real(dp),parameter :: rsfacm3=rsfac**(-3) 300 real(dp) :: a0,a1,a2,a3,b1,b2,b3,b4,d1,d1m1,d2d1drs2,d2d1drsdf,d2excdf2 301 real(dp) :: d2excdrs2,d2excdrsdf,d2excdz2,d2fxcdz2,d2n1drs2,d2n1drsdf,dd1df 302 real(dp) :: dd1drs,dexcdf,dexcdrs,dexcdz,dfxcdz,dn1df,dn1drs,dvxcdrs 303 real(dp) :: dvxcpdrho,dvxcpdz,excipt,fact,fxc,n1 304 real(dp) :: rhom1,rs,vxcp,zet,zetm,zetm_third 305 real(dp) :: zetp,zetp_third 306 character(len=500) :: message 307 !no_abirules 308 !Set a minimum rho below which terms are 0 309 real(dp),parameter :: rhotol=1.d-28 310 !real(dp) :: delta,rho,rho_dn,rho_dnm,rho_dnp,rho_up,rho_upm,rho_upp,zeta_mean 311 312 ! ************************************************************************* 313 314 !Checks the compatibility between the presence of dvxc and ndvxc 315 if(ndvxc /=0 .neqv. present(dvxc))then 316 message = 'If ndvxc/=0 there must be the optional argument dvxc' 317 MSG_BUG(message) 318 end if 319 320 !Checks the compatibility between the inputs and the presence of the optional arguments 321 if(abs(order) <= 1 .and. ndvxc /= 0)then 322 write(message, '(4a,i0)' )ch10,& 323 & 'The order chosen does not need the presence',ch10,& 324 & 'of the vector dvxc, that is needed only with |order|>1 , while we have',order 325 MSG_BUG(message) 326 end if 327 328 if(nspden == 1 .and. ndvxc /=0 .and. ndvxc /= 2)then 329 write(message,'(a,i0)')' Once nspden=1 we must have ndvxc=2, while we have',ndvxc 330 MSG_BUG(message) 331 end if 332 333 if(nspden == 2 .and. ndvxc /=0 .and. ndvxc /= 3)then 334 write(message, '(a,i0)' )' Once nspden=2 we must have ndvxc=3, while we have',ndvxc 335 MSG_BUG(message) 336 end if 337 338 339 !Although fact is parameter value, some compilers are not able to evaluate 340 !it at compile time. 341 fact=one/(two**(four*third)-two) 342 343 !DEBUG 344 !Finite-difference debugging, do not take away 345 !debug=1 346 !zeta_mean=0.1_dp 347 !delta=0.0001 348 !if(debug==1)then 349 !do ipts=1,npts,5 350 !rho=ipts*0.01_dp 351 !rho_up=rho*(one+zeta_mean)*half 352 !rho_dn=rho*(one-zeta_mean)*half 353 !rho_upp=rho_up+delta 354 !rho_upm=rho_up-delta 355 !rho_dnp=rho_dn+delta 356 !rho_dnm=rho_dn-delta 357 !First possibility : vary rho up , and then rho down 358 !zeta(ipts )=(rho_up -rho_dn )/(rho_up +rho_dn ) 359 !zeta(ipts+1)=(rho_upp-rho_dn )/(rho_upp+rho_dn ) 360 !zeta(ipts+2)=(rho_upm-rho_dn )/(rho_upm+rho_dn ) 361 !zeta(ipts+3)=(rho_up -rho_dnp)/(rho_up +rho_dnp) 362 !zeta(ipts+4)=(rho_up -rho_dnm)/(rho_up +rho_dnm) 363 !rspts(ipts )=rsfac*(rho_up +rho_dn )**(-third) 364 !rspts(ipts+1)=rsfac*(rho_upp+rho_dn )**(-third) 365 !rspts(ipts+2)=rsfac*(rho_upm+rho_dn )**(-third) 366 !rspts(ipts+3)=rsfac*(rho_up +rho_dnp)**(-third) 367 !rspts(ipts+4)=rsfac*(rho_up +rho_dnm)**(-third) 368 !DEBUGBUG : another possibility : vary rho and zeta 369 !zeta(ipts+1)=zeta(ipts ) 370 !zeta(ipts+2)=zeta(ipts ) 371 !zeta(ipts+3)=zeta(ipts )+delta 372 !zeta(ipts+4)=zeta(ipts )-delta 373 !rspts(ipts+1)=rsfac*(rho+delta)**(-third) 374 !rspts(ipts+2)=rsfac*(rho-delta )**(-third) 375 !rspts(ipts+3)=rspts(ipts ) 376 !rspts(ipts+4)=rspts(ipts ) 377 !ENDDEBUGBUG 378 !end do 379 !end if 380 !nspden=2 381 !order=2 382 !ENDDEBUG 383 384 if (nspden==1) then 385 ! separate cases with respect to order 386 if(order==-2) then 387 ! No spin-polarization so skip steps related to zeta not 0 388 do ipts=1,npts 389 390 rs=rspts(ipts) 391 n1=a0p+rs*(a1p+rs*(a2p+rs*a3p)) 392 d1=rs*(b1p+rs*(b2p+rs*(b3p+rs*b4p))) 393 d1m1=one/d1 394 395 ! Exchange-correlation energy 396 excipt=-n1*d1m1 397 exc(ipts)=excipt 398 399 ! Exchange-correlation potential 400 dn1drs=a1p+rs*(2._dp*a2p+rs*(3._dp*a3p)) 401 dd1drs=b1p+rs*(2._dp*b2p+rs*(3._dp*b3p+rs*(4._dp*b4p))) 402 403 ! dexcdrs is d(exc)/d(rs) 404 dexcdrs=-(dn1drs+excipt*dd1drs)*d1m1 405 vxc(ipts,1)=excipt-third*rs*dexcdrs 406 407 ! If the exchange-correlation kernel is needed 408 409 d2n1drs2=2._dp*a2p+rs*(6._dp*a3p) 410 d2d1drs2=2._dp*b2p+rs*(6._dp*b3p+rs*(12._dp*b4p)) 411 ! d2excdrs2 is d2(exc)/d(rs)2 412 d2excdrs2=-(d2n1drs2+2._dp*dexcdrs*dd1drs+excipt*d2d1drs2)*d1m1 413 dvxcdrs=third*(2.0_dp*dexcdrs-rs*d2excdrs2) 414 ! And d(vxc)/d(rho)=(-rs/(3*rho))*d(vxc)/d(rs) 415 dvxc(ipts,1)= -rs**4*rsfacm3*third*dvxcdrs 416 417 ! dn1df=d(n1)/d(fxc) and dd1df=d(d1)/d(fxc) 418 dn1df=da0+rs*(da1+rs*(da2+rs*da3)) 419 dd1df=rs*(db1+rs*(db2+rs*(db3+rs*db4))) 420 dexcdf=-(dn1df+excipt*dd1df)*d1m1 421 ! d2(fxc)/d(zeta)2 422 d2fxcdz2=ft*third*(alpha_zeta**2)*2._dp*fact 423 ! d2(exc)/d(zeta)2 424 d2excdz2=d2fxcdz2*dexcdf 425 rhom1=rsfacm3*rs**3 426 dvxc(ipts,2)= dvxc(ipts,1) - d2excdz2*rhom1 427 end do 428 else if(order**2>1) then 429 ! No spin-polarization so skip steps related to zeta not 0 430 do ipts=1,npts 431 432 rs=rspts(ipts) 433 n1=a0p+rs*(a1p+rs*(a2p+rs*a3p)) 434 d1=rs*(b1p+rs*(b2p+rs*(b3p+rs*b4p))) 435 d1m1=one/d1 436 437 ! Exchange-correlation energy 438 excipt=-n1*d1m1 439 exc(ipts)=excipt 440 441 ! Exchange-correlation potential 442 dn1drs=a1p+rs*(2._dp*a2p+rs*(3._dp*a3p)) 443 dd1drs=b1p+rs*(2._dp*b2p+rs*(3._dp*b3p+rs*(4._dp*b4p))) 444 445 ! dexcdrs is d(exc)/d(rs) 446 dexcdrs=-(dn1drs+excipt*dd1drs)*d1m1 447 vxc(ipts,1)=excipt-third*rs*dexcdrs 448 449 ! If the exchange-correlation kernel is needed 450 d2n1drs2=2._dp*a2p+rs*(6._dp*a3p) 451 d2d1drs2=2._dp*b2p+rs*(6._dp*b3p+rs*(12._dp*b4p)) 452 ! d2excdrs2 is d2(exc)/d(rs)2 453 d2excdrs2=-(d2n1drs2+2._dp*dexcdrs*dd1drs+excipt*d2d1drs2)*d1m1 454 dvxcdrs=third*(2.0_dp*dexcdrs-rs*d2excdrs2) 455 ! And d(vxc)/d(rho)=(-rs/(3*rho))*d(vxc)/d(rs) 456 dvxc(ipts,1)= -rs**4*rsfacm3*third*dvxcdrs 457 458 end do 459 else 460 ! No spin-polarization so skip steps related to zeta not 0 461 do ipts=1,npts 462 463 rs=rspts(ipts) 464 n1=a0p+rs*(a1p+rs*(a2p+rs*a3p)) 465 d1=rs*(b1p+rs*(b2p+rs*(b3p+rs*b4p))) 466 d1m1=one/d1 467 468 ! Exchange-correlation energy 469 excipt=-n1*d1m1 470 exc(ipts)=excipt 471 472 ! Exchange-correlation potential 473 dn1drs=a1p+rs*(2._dp*a2p+rs*(3._dp*a3p)) 474 dd1drs=b1p+rs*(2._dp*b2p+rs*(3._dp*b3p+rs*(4._dp*b4p))) 475 476 ! dexcdrs is d(exc)/d(rs) 477 dexcdrs=-(dn1drs+excipt*dd1drs)*d1m1 478 vxc(ipts,1)=excipt-third*rs*dexcdrs 479 end do 480 481 end if 482 483 484 ! Allows for nspden==1, in the case of testing nspden=1 against nspden=2 485 else if (nspden<=2) then 486 487 488 ! DEBUG 489 ! do not take away : allows to compare nspden=1 and nspden=2 coding 490 ! if (nspden==1)then 491 ! zeta(:)=zero 492 ! end if 493 ! ENDDEBUG 494 ! separate cases with respect to order 495 if(abs(order)>1) then 496 ! Allow for spin polarization. This part could be optimized for speed. 497 do ipts=1,npts 498 499 rs=rspts(ipts) 500 zet=zeta(ipts) 501 zetp=one+zet*alpha_zeta 502 zetm=one-zet*alpha_zeta 503 zetp_third=zetp**third 504 zetm_third=zetm**third 505 ! Exchange energy spin interpolation function f(zeta) 506 fxc=( zetp*zetp_third + zetm*zetm_third - two ) *fact 507 508 a0=a0p+fxc*da0 509 a1=a1p+fxc*da1 510 a2=a2p+fxc*da2 511 a3=a3p+fxc*da3 512 b1=b1p+fxc*db1 513 b2=b2p+fxc*db2 514 b3=b3p+fxc*db3 515 b4=b4p+fxc*db4 516 517 n1= a0+rs*(a1+rs*(a2+rs*a3)) 518 d1=rs*(b1+rs*(b2+rs*(b3+rs*b4))) 519 d1m1=one/d1 520 521 ! Exchange-correlation energy 522 excipt=-n1*d1m1 523 exc(ipts)=excipt 524 525 ! Exchange-correlation potential 526 dn1drs=a1+rs*(2._dp*a2+rs*(3._dp*a3)) 527 dd1drs=b1+rs*(2._dp*b2+rs*(3._dp*b3+rs*(4._dp*b4))) 528 ! dexcdrs is d(exc)/d(rs) 529 dexcdrs=-(dn1drs+excipt*dd1drs)*d1m1 530 531 ! Only vxcp contributes when paramagnetic 532 vxcp=excipt-third*rs*dexcdrs 533 534 ! d(fxc)/d(zeta) (which is 0 at zeta=0) 535 dfxcdz=ft*alpha_zeta*(zetp_third-zetm_third)*fact 536 537 ! dn1df=d(n1)/d(fxc) and dd1df=d(d1)/d(fxc) 538 dn1df=da0+rs*(da1+rs*(da2+rs*da3)) 539 dd1df=rs*(db1+rs*(db2+rs*(db3+rs*db4))) 540 541 ! dexcdz is d(exc)/d(zeta) 542 dexcdf=-(dn1df+excipt*dd1df)*d1m1 543 dexcdz=dfxcdz*dexcdf 544 545 ! Compute Vxc for both spin channels 546 547 vxc(ipts,1)=vxcp - (zet-one)*dexcdz 548 vxc(ipts,2)=vxcp - (zet+one)*dexcdz 549 550 ! DEBUG Allow to check the variation of rho and zeta 551 ! vxc(ipts,1)=vxcp 552 ! vxc(ipts,2)=dexcdz 553 ! ENDDEBUG 554 ! Compute second derivative with respect to rho 555 d2n1drs2=2._dp*a2+rs*(6._dp*a3) 556 d2d1drs2=2._dp*b2+rs*(6._dp*b3+rs*(12._dp*b4)) 557 ! d2excdrs2 is d2(exc)/d(rs)2 558 d2excdrs2=-(d2n1drs2+two*dexcdrs*dd1drs+excipt*d2d1drs2)*d1m1 559 dvxcdrs=third*(two*dexcdrs-rs*d2excdrs2) 560 ! And d(vxc)/d(rho) paramagnetic =(-rs/(3*rho))*d(vxcp)/d(rs) 561 ! remember : 1/rho=(4pi/3)*rs**3=rsfacm3*rs**3 562 rhom1=rsfacm3*rs**3 563 dvxcpdrho= -rs*rhom1*third * dvxcdrs 564 565 ! Compute mixed second derivative with respect to rho and zeta 566 d2n1drsdf=da1+rs*(2._dp*da2+rs*(3._dp*da3)) 567 d2d1drsdf=db1+rs*(2._dp*db2+rs*(3._dp*db3+rs*(4._dp*db4))) 568 ! d2excdrsdf is d2(exc)/d(rs)df 569 d2excdrsdf=-(d2n1drsdf+dexcdrs*dd1df+dexcdf*dd1drs+excipt*d2d1drsdf)*d1m1 570 ! d(vxc)/d(zeta) paramagnetic 571 dvxcpdz=dexcdz-third*rs*dfxcdz*d2excdrsdf 572 573 ! Compute second derivative with respect to zeta 574 ! the second derivative of n1 and d1 wrt f vanishes 575 d2excdf2=-(two*dexcdf*dd1df)*d1m1 576 ! d2(fxc)/d(zeta)2 577 d2fxcdz2=ft*third*(alpha_zeta**2)*(zetp_third**(-2)+zetm_third**(-2))*fact 578 ! d2(exc)/d(zeta)2 579 d2excdz2=d2fxcdz2*dexcdf+dfxcdz**2*d2excdf2 580 581 ! Compute now the three second derivatives of the Exc energy with respect 582 ! to : wrt twice spin-up ; wrt spin-up and spin-dn ; wrt twice spin-down 583 dvxc(ipts,1)= dvxcpdrho & 584 & +two*rhom1*( one-zet)*(dvxcpdz-dexcdz) & 585 & +d2excdz2*rhom1*(one-zet)**2 586 dvxc(ipts,2)= dvxcpdrho & 587 & +two*rhom1*( -zet)*(dvxcpdz-dexcdz) & 588 & +d2excdz2*rhom1*(one-zet)*(-one-zet) 589 ! if(nspden==2)then 590 dvxc(ipts,3)= dvxcpdrho & 591 & +two*rhom1*(-one-zet)*(dvxcpdz-dexcdz) & 592 & +d2excdz2*rhom1*(-one-zet)**2 593 ! else 594 ! ! For testing purposes, need the spin-averaged quantity 595 ! dvxc(ipts,1)= ( dvxc(ipts,1) + dvxc(ipts,2) ) * half 596 ! end if 597 598 ! DEBUG Allow to check the variation of rho and zeta 599 ! dvxc(ipts,1)=dvxcpdrho 600 ! dvxc(ipts,2)=d2excdz2 601 ! dvxc(ipts,3)=dvxcpdz 602 ! ENDDEBUG 603 end do 604 else 605 ! Allow for spin polarization. This part could be optimized for speed. 606 do ipts=1,npts 607 608 rs=rspts(ipts) 609 zet=zeta(ipts) 610 zetp=one+zet*alpha_zeta 611 zetm=one-zet*alpha_zeta 612 zetp_third=zetp**third 613 zetm_third=zetm**third 614 ! Exchange energy spin interpolation function f(zeta) 615 fxc=( zetp*zetp_third + zetm*zetm_third - two ) *fact 616 617 a0=a0p+fxc*da0 618 a1=a1p+fxc*da1 619 a2=a2p+fxc*da2 620 a3=a3p+fxc*da3 621 b1=b1p+fxc*db1 622 b2=b2p+fxc*db2 623 b3=b3p+fxc*db3 624 b4=b4p+fxc*db4 625 626 n1= a0+rs*(a1+rs*(a2+rs*a3)) 627 d1=rs*(b1+rs*(b2+rs*(b3+rs*b4))) 628 d1m1=one/d1 629 630 ! Exchange-correlation energy 631 excipt=-n1*d1m1 632 exc(ipts)=excipt 633 634 ! Exchange-correlation potential 635 dn1drs=a1+rs*(2._dp*a2+rs*(3._dp*a3)) 636 dd1drs=b1+rs*(2._dp*b2+rs*(3._dp*b3+rs*(4._dp*b4))) 637 ! dexcdrs is d(exc)/d(rs) 638 dexcdrs=-(dn1drs+excipt*dd1drs)*d1m1 639 640 ! Only vxcp contributes when paramagnetic 641 vxcp=excipt-third*rs*dexcdrs 642 643 ! d(fxc)/d(zeta) (which is 0 at zeta=0) 644 dfxcdz=ft*alpha_zeta*(zetp_third-zetm_third)*fact 645 646 ! dn1df=d(n1)/d(fxc) and dd1df=d(d1)/d(fxc) 647 dn1df=da0+rs*(da1+rs*(da2+rs*da3)) 648 dd1df=rs*(db1+rs*(db2+rs*(db3+rs*db4))) 649 650 ! dexcdz is d(exc)/d(zeta) 651 dexcdf=-(dn1df+excipt*dd1df)*d1m1 652 dexcdz=dfxcdz*dexcdf 653 654 ! Compute Vxc for both spin channels 655 656 vxc(ipts,1)=vxcp - (zet-one)*dexcdz 657 vxc(ipts,2)=vxcp - (zet+one)*dexcdz 658 659 ! DEBUG Allow to check the variation of rho and zeta 660 ! vxc(ipts,1)=vxcp 661 ! vxc(ipts,2)=dexcdz 662 ! ENDDEBUG 663 end do 664 end if 665 else 666 667 ! Disallowed value for nspden 668 write(message, '(3a,i0)' )& 669 & ' Argument nspden must be 1 or 2; ',ch10,& 670 & ' Value provided as argument was ',nspden 671 MSG_BUG(message) 672 end if 673 674 !DEBUG 675 !Finite-difference debugging, do not take away 676 !if(debug==1)then 677 !write(std_out,*)' delta =',delta 678 !do ipts=1,npts,5 679 !rho=(rspts(ipts)/rsfac)**(-3) 680 !write(std_out,'(a,i5,a,2es16.8)' ) ' Point number',ipts,' with rho,zeta=',rho,zeta(ipts) 681 !write(std_out,'(3es16.8)' )exc(ipts)*rho,vxc(ipts,1),vxc(ipts,2) 682 !write(std_out,'(3es16.8)' )dvxc(ipts,1),dvxc(ipts,3),dvxc(ipts,2) 683 !write(std_out,'(3es16.8)' )exc(ipts)*rho,& 684 !& ( exc(ipts+1)*(rho+delta) - exc(ipts+2)*(rho-delta) )/2._dp/delta,& 685 !& ( exc(ipts+3)*(rho+delta) - exc(ipts+4)*(rho-delta) )/2._dp/delta 686 !write(std_out,'(4es16.8)' )& 687 !& ( vxc(ipts+1,1) - vxc(ipts+2,1) )/2._dp/delta,& 688 !& ( vxc(ipts+3,2) - vxc(ipts+4,2) )/2._dp/delta,& 689 !& ( vxc(ipts+3,1) - vxc(ipts+4,1) )/2._dp/delta,& 690 !& ( vxc(ipts+1,2) - vxc(ipts+2,2) )/2._dp/delta 691 !end do 692 !stop 693 !end if 694 !ENDDEBUG 695 696 !DEBUG 697 !if(order==-2)then 698 !write(std_out,*)' xcspol : ipts,npts ',ipts,npts 699 !write(std_out,*)dvxcdrs,d2excdz2,d2fxcdz2,dexcdf 700 !write(std_out,*)rhom1 701 !write(std_out,*)dvxc(1000,1),dvxc(1000,2) 702 !stop 703 !end if 704 !ENDDEBUG 705 706 end subroutine xcspol
ABINIT/xctetr [ Functions ]
NAME
xctetr
FUNCTION
Returns exc, vxc, and d(vxc)/d($\rho$) from input $\rho$. Also returns $d^2(Vxc)/d(\rho)^2$ as needed for third-order DFPT
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) if(order>2) d2vxc(npt)=derivative d$^2$(Vxc)/d$(\rho)^2$ (hartree*bohr^6)
NOTES
Teter exchange and correlation (xc)--Mike Teter s fit to Ceperly-Alder electron gas energy data. Data from D.M. Ceperley and B.J. Alder, Phys. Rev. Lett. 45, 566 (1980) [[cite:Ceperley1980]] and private communication from authors. This form is based on Mike Teter s rational polynomial exc=-(a0+a1*rs+a2*rs**2+a3*rs**3)/(b1*rs+b2*rs**2+b3*rs**3+b4*rs**4) where the parameters of the fit are fit to reproduce Ceperley-Alder data and the high density limit (rs->0) of the electron gas (pure exchange). rs = $(3/(4\pi))^{1/3} * \rho(r)^{-1/3}$. b1 must be 1 and a0 must be $(3/4)(3/(2\pi))^{2/3}$. Fit is by Mike Teter, Corning Incorporated. Note that d(vxc)/d($\rho$) gets a little wild at small rho. d$^2$(Vxc)/d$(\rho)^2$ is probably wilder. Some notation: (XG 990224, sign convention should be changed, see xcspol.f) $Exc = N1/D1$ with $N1=-(a0+a1*rs+...)$ given above and $D1= (b1*rs+b2*rs^2+...)$ also given above. $Vxc = N2/D1^2$ with $N2=d(N1)/d(rs)$. $d(Vxc)/d(rs)=(N3-D3*(2*N2/D1))/D1^2 with N3=d(N2)/d(rs)$ and $D3=d(D1)/d(rs)$. $d(Vxc)/d(\rho) = (-rs/(3*\rho))* d(Vxc)/d(rs)$. $d^2(Vxc)/d(rs)^2 = (N4-2*(2*N3*D3+N2*D4-3*N2*D3^2/D1)/D1)/D1^2$ with $N4=d(N3)/d(rs), D4=d(D3)/d(rs)$. $d^2(Vxc)/d(\rho)^2= rs/(3*\rho)^2)*(4*d(Vxc)/d(rs)+rs*d^2(Vxc)/d(rs)^2)$.
PARENTS
drivexc
CHILDREN
SOURCE
764 subroutine xctetr(exc,npt,order,rhor,rspts,vxc,& !Mandatory arguments 765 & d2vxc,dvxc) !Optional arguments 766 767 768 !This section has been created automatically by the script Abilint (TD). 769 !Do not modify the following lines by hand. 770 #undef ABI_FUNC 771 #define ABI_FUNC 'xctetr' 772 !End of the abilint section 773 774 implicit none 775 776 !Arguments ------------------------------------ 777 !scalars 778 integer,intent(in) :: npt,order 779 !arrays 780 real(dp),intent(in) :: rhor(npt),rspts(npt) 781 real(dp),intent(out) :: exc(npt),vxc(npt) 782 real(dp),intent(out),optional :: d2vxc(npt),dvxc(npt) 783 784 !Local variables------------------------------- 785 !rsfac=(3/(4 Pi))^(1/3) 786 !Mike Teter s parameters: (keep 8 digits after decimal) 787 !(a0=(3/4)(3/(2 Pi))^(2/3) 788 !scalars 789 integer :: ipt 790 real(dp),parameter :: a0=.4581652932831429_dp,a1=2.40875407_dp,a2=.88642404_dp 791 real(dp),parameter :: a3=.02600342_dp,b1=1.0_dp,b2=4.91962865_dp 792 real(dp),parameter :: b3=1.34799453_dp,b4=.03120453_dp,c1=4._dp*a0*b1/3.0_dp 793 real(dp),parameter :: c2=5.0_dp*a0*b2/3.0_dp+a1*b1 794 real(dp),parameter :: c3=2.0_dp*a0*b3+4.0_dp*a1*b2/3.0_dp+2.0_dp*a2*b1/3.0_dp 795 real(dp),parameter :: c4=7.0_dp*a0*b4/3.0_dp+5.0_dp*a1*b3/3.0_dp+a2*b2+a3*b1/3.0_dp 796 real(dp),parameter :: c5=2.0_dp*a1*b4+4.0_dp*a2*b3/3.0_dp+2.0_dp*a3*b2/3.0_dp 797 real(dp),parameter :: c6=5.0_dp*a2*b4/3.0_dp+a3*b3,c7=4.0_dp*a3*b4/3.0_dp 798 real(dp),parameter :: rsfac=0.6203504908994000_dp 799 real(dp) :: d1,d1m1,d2vxcr,d3,d4,dvxcdr,n1,n2,n3,n4,rhom1,rs 800 character(len=500) :: message 801 802 ! ************************************************************************* 803 ! 804 !Checks the values of order 805 if(order<0 .or. order>3)then 806 write(message, '(a,a,a,i6)' )& 807 & 'With Teter 91 Ceperley-Alder xc functional, the only',ch10,& 808 & 'allowed values for order are 0, 1, 2 or 3, while it is found to be',order 809 MSG_BUG(message) 810 end if 811 812 !Checks the compatibility between the order and the presence of the optional arguments 813 if(order /=3 .and. present(d2vxc))then 814 write(message, '(a,a,a,i6)' )& 815 & 'The order chosen does not need the presence',ch10,& 816 & 'of the vector d2vxc, that is needed only with order=3, while we have',order 817 MSG_BUG(message) 818 end if 819 820 if(order <= 1 .and. present(dvxc))then 821 write(message, '(a,a,a,i6)' )& 822 & 'The order chosen does not need the presence',ch10,& 823 & 'of the vector dvxc, that is needed with order > 1, while we have',order 824 MSG_BUG(message) 825 end if 826 827 !separated cases with respect to order 828 829 if (order<=1) then 830 ! Loop over grid points 831 do ipt=1,npt 832 rs=rspts(ipt) 833 n1=-(a0+rs*(a1+rs*(a2+rs*a3))) 834 d1=rs*(b1+rs*(b2+rs*(b3+rs*b4))) 835 d1m1=1.0_dp/d1 836 n2=-rs*(c1+rs*(c2+rs*(c3+rs*(c4+rs*(c5+rs*(c6+rs*c7)))))) 837 ! 838 ! Exchange-correlation energy 839 exc(ipt)=n1*d1m1 840 ! 841 ! Exchange-correlation potential 842 vxc(ipt)=n2*d1m1**2 843 end do 844 else if (order>2) then 845 ! Loop over grid points 846 do ipt=1,npt 847 rs=rspts(ipt) 848 n1=-(a0+rs*(a1+rs*(a2+rs*a3))) 849 d1=rs*(b1+rs*(b2+rs*(b3+rs*b4))) 850 d1m1=1.0_dp/d1 851 n2=-rs*(c1+rs*(c2+rs*(c3+rs*(c4+rs*(c5+rs*(c6+rs*c7)))))) 852 ! 853 ! Exchange-correlation energy 854 exc(ipt)=n1*d1m1 855 ! 856 ! Exchange-correlation potential 857 vxc(ipt)=n2*d1m1**2 858 ! Assemble derivative of vxc wrt rs 859 n3=-(c1+rs*(2._dp*c2+rs*(3._dp*c3+rs*(4._dp*c4+rs*(5._dp*c5+& 860 & rs*(6._dp*c6+rs*(7._dp*c7))))))) 861 d3=b1+rs*(2._dp*b2+rs*(3._dp*b3+rs*(4._dp*b4))) 862 dvxcdr=(n3-d3*(2._dp*n2*d1m1))*d1m1**2 863 rhom1=1.0_dp/rhor(ipt) 864 ! 865 ! derivative of vxc wrt rho 866 dvxc(ipt)=-dvxcdr*rs*third*rhom1 867 ! 868 869 ! Assemble derivative d^2(Vxc)/d(rs)^2 870 n4=-(2.0_dp*c2+rs*(6.0_dp*c3+rs*(12.0_dp*c4+rs*(20.0_dp*c5+& 871 & rs*(30.0_dp*c6+rs*(42.0_dp*c7)))))) 872 d4=2.0_dp*b2+rs*(6.0_dp*b3+rs*(12.0_dp*b4)) 873 d2vxcr=(n4-2.0_dp*(2.0_dp*n3*d3+n2*d4-3.0_dp*n2*d3**2*d1m1)*d1m1)*d1m1**2 874 875 ! Derivative d^2(Vxc)/d(rho)^2 876 d2vxc(ipt)=(rs*third*rhom1)*(4.0_dp*dvxcdr+rs*d2vxcr)*third*rhom1 877 878 end do 879 else if (order>1) then 880 ! Loop over grid points 881 do ipt=1,npt 882 rs=rspts(ipt) 883 n1=-(a0+rs*(a1+rs*(a2+rs*a3))) 884 d1=rs*(b1+rs*(b2+rs*(b3+rs*b4))) 885 d1m1=1.0_dp/d1 886 n2=-rs*(c1+rs*(c2+rs*(c3+rs*(c4+rs*(c5+rs*(c6+rs*c7)))))) 887 ! 888 ! Exchange-correlation energy 889 exc(ipt)=n1*d1m1 890 ! 891 ! Exchange-correlation potential 892 vxc(ipt)=n2*d1m1**2 893 ! Assemble derivative of vxc wrt rs 894 n3=-(c1+rs*(2._dp*c2+rs*(3._dp*c3+rs*(4._dp*c4+rs*(5._dp*c5+& 895 & rs*(6._dp*c6+rs*(7._dp*c7))))))) 896 d3=b1+rs*(2._dp*b2+rs*(3._dp*b3+rs*(4._dp*b4))) 897 dvxcdr=(n3-d3*(2._dp*n2*d1m1))*d1m1**2 898 rhom1=1.0_dp/rhor(ipt) 899 ! 900 ! derivative of vxc wrt rho 901 dvxc(ipt)=-dvxcdr*rs*third*rhom1 902 ! 903 end do 904 end if 905 end subroutine xctetr
ABINIT/xctfw [ Functions ]
NAME
xctfw
FUNCTION
Add gradient part of the Thomas-Fermi-Weizsacker functional Perrot F., Phys. Rev. A20, 586-594 (1979) [[cite:Perrot1979]]
INPUTS
ndvxcdgr= size of dvxcdgr(npts,ndvxcdgr)
ngr2= size of grho2_updn(npts,ngr2)
npts= number of points to be computed
nspden=number if spin density component (necessarily 1 here)
grho2_updn(npts,ngr2)=square of the gradient of the spin-up,
and, if nspden==2, spin-down, and total density (Hartree/Bohr**2),
only used if gradient corrected functional (option=2,-2,-4 and 4 or beyond)
rho_updn(npts,nspden)=spin-up and spin-down density (Hartree/bohr**3)
temp= electronic temperature
usefxc=1 if free energy fxc is used
SIDE EFFECTS
The following arrays are modified (gradient correction added): dvxcdgr(npts,3)=partial derivative of the XC energy divided by the norm of the gradient exci(npts)=exchange-correlation energy density fxci(npts)=free energy energy density vxci(npts,nspden)=exchange-correlation potential
PARENTS
rhotoxc
CHILDREN
invcb
SOURCE
1362 subroutine xctfw(temp,exci,fxci,usefxc,rho_updn,vxci,npts,nspden,dvxcdgr,ndvxcdgr,grho2_updn,ngr2) 1363 1364 1365 !This section has been created automatically by the script Abilint (TD). 1366 !Do not modify the following lines by hand. 1367 #undef ABI_FUNC 1368 #define ABI_FUNC 'xctfw' 1369 !End of the abilint section 1370 1371 implicit none 1372 1373 !Arguments ------------------------------------ 1374 !scalars 1375 integer,intent(in) :: ndvxcdgr,ngr2,npts,nspden,usefxc 1376 real(dp),intent(in) :: temp 1377 !arrays 1378 real(dp),intent(in) :: grho2_updn(npts,ngr2),rho_updn(npts,nspden) 1379 real(dp),intent(inout) :: dvxcdgr(npts,ndvxcdgr),fxci(npts*usefxc),exci(npts),vxci(npts,nspden) 1380 1381 !Local variables------------------------------- 1382 !scalars 1383 integer :: iperrot,ipts 1384 logical :: has_fxc,has_dvxcdgr 1385 real(dp) :: etfw,rho,rho_inv,rhomot,yperrot0,vtfw 1386 real(dp) :: yperrot,uperrot,dyperrotdn,duperrotdyperrot 1387 real(dp) :: hperrot,dhperrotdyperrot,dhperrotdn,dhperrotduperrot 1388 !arrays 1389 real(dp) :: wpy(0:7), wpu(0:7) 1390 real(dp),allocatable :: rho_updnm1_3(:,:) 1391 1392 ! ************************************************************************* 1393 1394 !DBG_ENTER('COLL') 1395 1396 has_fxc=(usefxc/=0) 1397 has_dvxcdgr=(ndvxcdgr/=0) 1398 1399 yperrot0=1.666081101_dp 1400 1401 wpy(0)=0.5_dp; wpy(1)=-0.1999176316_dp 1402 wpy(2)=0.09765615709_dp; wpy(3)=-0.06237609924_dp 1403 wpy(4)=0.05801466322_dp; wpy(5)=-0.04449287774_dp 1404 wpy(6)=0.01903211697_dp; wpy(7)=-0.003284096926_dp 1405 1406 wpu(0)=one/6._dp; wpu(1)=0.311590799_dp 1407 wpu(2)=3.295662439_dp; wpu(3)=-29.22038326_dp 1408 wpu(4)=116.1084531_dp; wpu(5)=-250.4543147_dp 1409 wpu(6)=281.433688_dp; wpu(7)=-128.8784806_dp 1410 1411 ABI_ALLOCATE(rho_updnm1_3,(npts,2)) 1412 1413 call invcb(rho_updn(:,1),rho_updnm1_3(:,1),npts) 1414 1415 do ipts=1,npts 1416 1417 rho =rho_updn(ipts,1) 1418 rhomot=rho_updnm1_3(ipts,1) 1419 rho_inv=rhomot*rhomot*rhomot 1420 1421 yperrot=pi*pi/sqrt2/temp**1.5*two*rho 1422 uperrot=yperrot**(2./3.) 1423 1424 dyperrotdn=pi*pi/sqrt2/temp**1.5*2.0_dp 1425 1426 hperrot=zero 1427 dhperrotdyperrot=zero 1428 dhperrotduperrot=zero 1429 if(yperrot<=yperrot0)then 1430 do iperrot=0,7 1431 hperrot=hperrot+wpy(iperrot)*yperrot**iperrot 1432 dhperrotdyperrot=dhperrotdyperrot+iperrot*wpy(iperrot)*yperrot**(iperrot-1) 1433 end do 1434 hperrot=one/12.0_dp*hperrot 1435 dhperrotdyperrot=one/12.0_dp*dhperrotdyperrot 1436 dhperrotdn=dhperrotdyperrot*dyperrotdn 1437 else 1438 do iperrot=0,7 1439 hperrot=hperrot+wpu(iperrot)/uperrot**(2*iperrot) 1440 dhperrotduperrot=dhperrotduperrot-2.*iperrot*wpu(iperrot)/uperrot**(2*iperrot+1) 1441 end do 1442 hperrot=one/12.0_dp*hperrot 1443 dhperrotduperrot=one/12.0_dp*dhperrotduperrot 1444 duperrotdyperrot=two/3._dp/yperrot**(1./3.) 1445 dhperrotdn=dhperrotduperrot*duperrotdyperrot*dyperrotdn 1446 end if 1447 1448 etfw=hperrot*grho2_updn(ipts,1)*rho_inv*rho_inv 1449 vtfw=-etfw + rho/hperrot*dhperrotdn*etfw 1450 1451 if(yperrot<=yperrot0)then 1452 exci(ipts) = exci(ipts) + etfw + 1.5_dp*yperrot*dhperrotdyperrot*grho2_updn(ipts,1)*rho_inv*rho_inv 1453 else 1454 exci(ipts) = exci(ipts) + etfw + uperrot*dhperrotduperrot*grho2_updn(ipts,1)*rho_inv*rho_inv 1455 end if 1456 vxci(ipts,1) = vxci(ipts,1) + vtfw 1457 if (has_fxc) fxci(ipts) = fxci(ipts) + etfw 1458 if (has_dvxcdgr) dvxcdgr(ipts,1)= dvxcdgr(ipts,1)+two*hperrot*rho_inv 1459 1460 end do 1461 1462 ABI_DEALLOCATE(rho_updnm1_3) 1463 1464 !DBG_EXIT('COLL') 1465 1466 end subroutine xctfw
ABINIT/xcwign [ Functions ]
NAME
xcwign
FUNCTION
Returns exc, vxc, and eventually d(vxc)/d($\rho$) from input $\rho$. Wigner exchange and correlation (xc)--see e.g. David Pines, Elementary Excitations in Solids, p. 94, NY 1964. Expression is exc=-(0.44)/(rs+7.8)-efac/rs (hartree), efac below. rs = $(3/(4\pi))^{1/3}* \rho (r)^{-1/3}$.
INPUTS
npt=number of real space points on which density is provided order=gives the maximal derivative of Exc computed. 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
936 subroutine xcwign(exc,npt,order,rspts,vxc,& !Mandatory arguments 937 & dvxc) !Optional arguments 938 939 940 !This section has been created automatically by the script Abilint (TD). 941 !Do not modify the following lines by hand. 942 #undef ABI_FUNC 943 #define ABI_FUNC 'xcwign' 944 !End of the abilint section 945 946 implicit none 947 948 !Arguments ------------------------------------ 949 !scalars 950 integer,intent(in) :: npt,order 951 !arrays 952 real(dp),intent(in) :: rspts(npt) 953 real(dp),intent(out) :: exc(npt),vxc(npt) 954 real(dp),intent(out),optional :: dvxc(npt) 955 956 !Local variables------------------------------- 957 !c1 and c2 are the Wigner parameters in hartree and bohr resp. 958 !scalars 959 integer :: ipt 960 real(dp),parameter :: c1=0.44_dp,c2=7.8_dp,c4_3=4.0_dp/3.0_dp 961 real(dp),parameter :: c8_27=8.0_dp/27.0_dp 962 real(dp) :: dfac,efac,rs,rsc2m1,rsm1,vfac,vxcnum 963 character(len=500) :: message 964 965 ! ************************************************************************* 966 967 !Checks the values of order 968 if(order<0 .or. order>2)then 969 write(message, '(a,a,a,i0)' )& 970 & 'With Wigner xc functional, the only',ch10,& 971 & 'allowed values for order are 0, 1 or 2, while it is found to be',order 972 MSG_BUG(message) 973 end if 974 975 !Checks the compatibility between the order and the presence of the optional arguments 976 if(order <= 1 .and. present(dvxc))then 977 write(message, '(a,a,a,i3)' )& 978 & 'The order chosen does not need the presence',ch10,& 979 & 'of the vector dvxc, that is needed only with order=2 , while we have',order 980 MSG_BUG(message) 981 end if 982 983 !Compute vfac=(3/(2*Pi))^(2/3) 984 vfac=(1.5_dp/pi)**(2.0_dp/3.0_dp) 985 !Compute efac=(3/4)*vfac 986 efac=0.75_dp*vfac 987 !Compute dfac=(4*Pi/9)*vfac 988 dfac=(4.0_dp*pi/9.0_dp)*vfac 989 990 !separate cases with respect to order 991 if (order==2) then 992 993 ! Loop over grid points 994 do ipt=1,npt 995 rs=rspts(ipt) 996 rsm1=1.0_dp/rs 997 rsc2m1=1.0_dp/(rs+c2) 998 ! compute energy density (hartree) 999 exc(ipt)=-c1*rsc2m1-efac*rsm1 1000 vxcnum=-(c4_3*rs+c2)*c1 1001 ! compute potential (hartree) 1002 vxc(ipt)=vxcnum*rsc2m1**2-vfac*rsm1 1003 ! compute d(vxc)/d(rho) (hartree*bohr^3) 1004 dvxc(ipt)=-(c8_27*pi)*(c1*rs**4)*(rs+rs+c2)*rsc2m1**3-dfac*rs**2 1005 end do 1006 else 1007 1008 ! Loop over grid points 1009 do ipt=1,npt 1010 rs=rspts(ipt) 1011 rsm1=1.0_dp/rs 1012 rsc2m1=1.0_dp/(rs+c2) 1013 ! compute energy density (hartree) 1014 exc(ipt)=-c1*rsc2m1-efac*rsm1 1015 vxcnum=-(c4_3*rs+c2)*c1 1016 ! compute potential (hartree) 1017 vxc(ipt)=vxcnum*rsc2m1**2-vfac*rsm1 1018 end do 1019 1020 end if 1021 ! 1022 end subroutine xcwign
ABINIT/xcxalp [ Functions ]
NAME
xcxalp
FUNCTION
Returns exc, vxc, and eventually d(vxc)/d($\rho$) from input $\rho$.
"X$\alpha$" method is used in this subroutine:
a single fixed value is chosen for "alpha", set below.
Expression is exc=-alpha*efac/rs (hartree), efac below.
rs = $(3/(4\pi))^{1/3}* \rho (r)^{-1/3}$.
INPUTS
npt=number of real space points on which density is provided order=gives the maximal derivative of Exc computed. rspts(npt)=Wigner-Seitz radii, at each point
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
1158 subroutine xcxalp(exc,npt,order,rspts,vxc, dvxc) ! dvxc is optional 1159 1160 1161 !This section has been created automatically by the script Abilint (TD). 1162 !Do not modify the following lines by hand. 1163 #undef ABI_FUNC 1164 #define ABI_FUNC 'xcxalp' 1165 !End of the abilint section 1166 1167 implicit none 1168 1169 !Arguments ------------------------------------ 1170 !scalars 1171 integer,intent(in) :: npt,order 1172 !arrays 1173 real(dp),intent(in) :: rspts(npt) 1174 real(dp),intent(out) :: exc(npt),vxc(npt) 1175 real(dp),intent(out),optional :: dvxc(npt) 1176 1177 !Local variables------------------------------- 1178 !Set value of alpha in "X-alpha" method 1179 !scalars 1180 integer :: ipt 1181 real(dp),parameter :: alpha=1.0_dp 1182 real(dp) :: dfac,efac,rs,rsm1,vfac 1183 character(len=500) :: message 1184 1185 ! ************************************************************************* 1186 1187 !Checks the values of order 1188 if(order<0 .or. order>2)then 1189 write(message, '(a,a,a,i3)' )& 1190 & 'With X-alpha xc functional, the only',ch10,& 1191 & 'allowed values for order are 0, 1 or 2, while it is found to be',order 1192 MSG_BUG(message) 1193 end if 1194 1195 !Compute vfac=(3/(2*Pi))^(2/3) 1196 vfac=(1.5_dp/pi)**(2.0_dp/3.0_dp) 1197 !Compute efac=(3/4)*vfac 1198 efac=0.75_dp*vfac 1199 !Compute dfac=(4*Pi/9)*vfac 1200 dfac=(4.0_dp*pi/9.0_dp)*vfac 1201 1202 !separate cases with respect to order 1203 if(order==2) then 1204 ! Loop over grid points 1205 do ipt=1,npt 1206 rs=rspts(ipt) 1207 rsm1=1.0_dp/rs 1208 ! compute energy density (hartree) 1209 exc(ipt)=-alpha*efac*rsm1 1210 ! compute potential (hartree) 1211 vxc(ipt)=-alpha*vfac*rsm1 1212 ! compute d(vxc)/d(rho) (hartree*bohr^3) 1213 dvxc(ipt)=-alpha*dfac*rs**2 1214 end do 1215 else 1216 ! Loop over grid points 1217 do ipt=1,npt 1218 rs=rspts(ipt) 1219 rsm1=1.0_dp/rs 1220 ! compute energy density (hartree) 1221 exc(ipt)=-alpha*efac*rsm1 1222 ! compute potential (hartree) 1223 vxc(ipt)=-alpha*vfac*rsm1 1224 end do 1225 end if 1226 ! 1227 end subroutine xcxalp