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-2022 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 .
SOURCE
16 #if defined HAVE_CONFIG_H 17 #include "config.h" 18 #endif 19 20 #include "abi_common.h" 21 22 module m_xclda 23 24 use defs_basis 25 use m_errors 26 use m_abicore 27 28 use m_numeric_tools, only : invcb 29 30 implicit none 31 32 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)
SOURCE
988 subroutine xchelu(exc,npt,order,rspts,vxc,dvxc) ! dvxc is optional 989 990 !Arguments ------------------------------------ 991 !scalars 992 integer,intent(in) :: npt,order 993 !arrays 994 real(dp),intent(in) :: rspts(npt) 995 real(dp),intent(out) :: exc(npt),vxc(npt) 996 real(dp),intent(out),optional :: dvxc(npt) 997 998 !Local variables------------------------------- 999 !aa and cc are H-L fitting parameters A and C (C in hartree) 1000 !rs = (3/(4 Pi))**(1/3) * rho(r)**(-1/3). 1001 !scalars 1002 integer :: ipt 1003 real(dp),parameter :: aa=21_dp,c1_21=one/21_dp,c4_9=4.0_dp/9.0_dp,cc=0.0225_dp 1004 real(dp) :: dfac,efac,rs,rsm1,vfac,xx 1005 character(len=500) :: message 1006 1007 ! ************************************************************************* 1008 1009 !Checks the values of order 1010 if(order<0 .or. order>2)then 1011 write(message, '(a,a,a,i0)' )& 1012 & 'With Hedin-Lundqvist xc functional, the only',ch10,& 1013 & 'allowed values for order are 0, 1 or 2, while it is found to be',order 1014 ABI_BUG(message) 1015 end if 1016 1017 !Compute vfac=(3/(2*Pi))^(2/3) 1018 vfac=(1.5_dp/pi)**(2.0_dp/3.0_dp) 1019 !Compute efac=(3/4)*vfac 1020 efac=0.75_dp*vfac 1021 !Compute dfac=(4*Pi/9)*vfac 1022 dfac=(4.0_dp*pi/9.0_dp)*vfac 1023 !separate cases with respect to order 1024 if (order==2) then 1025 ! Loop over grid points 1026 do ipt=1,npt 1027 rs=rspts(ipt) 1028 rsm1=one/rs 1029 ! compute energy density exc (hartree) 1030 xx=rs*c1_21 1031 exc(ipt)=-cc*((one+xx**3)*log(one+one/xx)+& 1032 & half*xx-xx*xx-third) - efac*rsm1 1033 ! compute xc potential d(rho*exc)/d(rho) (hartree) 1034 vxc(ipt)=-cc*log(one+aa*rsm1)-vfac*rsm1 1035 ! compute d(vxc)/d(rho) (hartree*bohr^3) 1036 dvxc(ipt)=-(rs**2)*((c4_9*pi)*cc*rs/(one+xx) + dfac) 1037 end do 1038 else 1039 ! Loop over grid points 1040 do ipt=1,npt 1041 rs=rspts(ipt) 1042 rsm1=one/rs 1043 ! compute energy density exc (hartree) 1044 xx=rs*c1_21 1045 exc(ipt)=-cc*((one+xx**3)*log(one+one/xx)+& 1046 & half*xx-xx*xx-third) - efac*rsm1 1047 ! compute xc potential d(rho*exc)/d(rho) (hartree) 1048 vxc(ipt)=-cc*log(one+aa*rsm1)-vfac*rsm1 1049 end do 1050 end if 1051 ! 1052 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.
SOURCE
1175 subroutine xclb(grho2_updn,npts,nspden,rho_updn,vxci) 1176 1177 !Arguments ------------------------------------ 1178 !scalars 1179 integer,intent(in) :: npts,nspden 1180 !arrays 1181 real(dp),intent(in) :: grho2_updn(npts,2*nspden-1),rho_updn(npts,nspden) 1182 real(dp),intent(inout) :: vxci(npts,nspden) 1183 1184 !Local variables------------------------------- 1185 !scalars 1186 integer :: ipts,ispden 1187 real(dp),parameter :: beta=0.05_dp 1188 real(dp) :: density,density_gradient,density_t13,s_g_sq,scaled_gradient 1189 real(dp) :: scaling_factor,vx_lb 1190 1191 ! ************************************************************************* 1192 1193 !DEBUG 1194 !write(std_out,*) ' %xclb: enter' 1195 !ENDDEBUG 1196 1197 !scale the spin densities for evaluating spin up or down exchange 1198 scaling_factor=one 1199 if(nspden == 2) scaling_factor=two 1200 1201 do ispden=1,nspden 1202 1203 do ipts=1,npts 1204 1205 density= scaling_factor * rho_updn(ipts,ispden) 1206 density_gradient= scaling_factor * sqrt(grho2_updn(ipts,ispden)) 1207 1208 density_t13= density**third 1209 scaled_gradient= density_gradient/max(density*density_t13,1.e-12_dp) 1210 1211 s_g_sq= scaled_gradient*scaled_gradient 1212 1213 vx_lb= -beta*density_t13 * s_g_sq/ & 1214 & (one+3.d0*beta* scaled_gradient*log(scaled_gradient+sqrt(one+s_g_sq*s_g_sq))) 1215 1216 vxci(ipts,ispden)=vxci(ipts,ispden)+vx_lb 1217 end do 1218 1219 end do 1220 1221 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)
SOURCE
73 subroutine xcpzca(exc,npt,order,rhor,rspts,vxc,& !Mandatory arguments 74 & dvxc) !Optional arguments 75 76 !Arguments ------------------------------------ 77 !scalars 78 integer,intent(in) :: npt,order 79 !arrays 80 real(dp),intent(in) :: rhor(npt),rspts(npt) 81 real(dp),intent(out) :: exc(npt),vxc(npt) 82 real(dp),intent(out),optional :: dvxc(npt) 83 84 !Local variables------------------------------- 85 !Perdew-Zunger parameters a, b, b1, b2, c, d, gamma 86 !scalars 87 integer :: ipt 88 real(dp),parameter :: aa=0.0311_dp,b1=1.0529_dp,b2=0.3334_dp,bb=-0.048_dp 89 real(dp),parameter :: c4_3=4.0_dp/3.0_dp,c7_6=7.0_dp/6.0_dp,cc=0.0020_dp 90 real(dp),parameter :: dd=-0.0116_dp,ga=-0.1423_dp 91 real(dp) :: den,den3,dfac,efac,logrs,rs,rsm1,t1,t2,vfac 92 character(len=500) :: message 93 94 ! ************************************************************************* 95 96 !Compute vfac=(3/(2*Pi))^(2/3) 97 vfac=(1.5_dp/pi)**(2.0_dp/3.0_dp) 98 !Compute efac=(3/4)*vfac 99 efac=0.75_dp*vfac 100 !Compute dfac=(4*Pi/9)*vfac 101 dfac=(4.0_dp*pi/9.0_dp)*vfac 102 103 !Checks the values of order 104 if(order<0 .or. order>2)then 105 write(message, '(a,a,a,i0)' )& 106 & 'With Perdew-Zunger Ceperley-Alder xc functional, the only',ch10,& 107 & 'allowed values for order are 0, 1 or 2, while it is found to be',order 108 ABI_BUG(message) 109 end if 110 111 !Checks the compatibility between the order and the presence of the optional arguments 112 if(order <= 1 .and. present(dvxc))then 113 write(message, '(a,a,a,i0)' )& 114 & 'The order chosen does not need the presence',ch10,& 115 & 'of the vector dvxc, that is needed only with order=2 , while we have',order 116 ABI_BUG(message) 117 end if 118 119 !separate cases with respect to order 120 if(order==2) then 121 ! Loop over grid points 122 do ipt=1,npt 123 rs=rspts(ipt) 124 rsm1=1.0_dp/rs 125 ! Consider two regimes: rs<1 or rs>=1 126 if (rs<1._dp) then 127 logrs=log(rs) 128 ! compute energy density exc (hartree) 129 exc(ipt)=(aa+cc*rs)*logrs+dd*rs+bb-efac*rsm1 130 ! compute potential vxc=d(rho*exc)/d(rho) (hartree) 131 vxc(ipt)=(aa+two_thirds*cc*rs)*logrs+(dd+dd-cc)*rs*third+& 132 & (bb-aa*third)-vfac*rsm1 133 ! compute d(vxc)/d(rho) (hartree*bohr^3) 134 dvxc(ipt)=-(3._dp*aa+(cc+dd+dd)*rs+2._dp*cc*rs*logrs)& 135 & /(9._dp*rhor(ipt))-dfac*rs**2 136 else if (rs<1000._dp) then 137 t1=b1*sqrt(rs) 138 t2=b2*rs 139 den=1._dp/(1._dp+t1+t2) 140 exc(ipt)=ga*den-efac*rsm1 141 vxc(ipt)=ga*(1._dp+c7_6*t1+c4_3*t2)*den**2-vfac*rsm1 142 den3=den**3 143 dvxc(ipt)=(ga*den3/(36._dp*rhor(ipt)))*(5._dp*t1+8._dp*t2+& 144 & 7._dp*t1**2+16._dp*t2**2+21._dp*t1*t2)-dfac*rs**2 145 else 146 t1=b1*sqrt(rs) 147 t2=b2*rs 148 den=1._dp/(1._dp+t1+t2) 149 exc(ipt)=ga*den-efac*rsm1 150 vxc(ipt)=ga*(1._dp+c7_6*t1+c4_3*t2)*den**2-vfac*rsm1 151 dvxc(ipt)=0._dp 152 end if 153 end do 154 else 155 ! Loop over grid points 156 do ipt=1,npt 157 rs=rspts(ipt) 158 rsm1=1.0_dp/rs 159 ! Consider two regimes: rs<1 or rs>=1 160 if (rs<1._dp) then 161 logrs=log(rs) 162 ! compute energy density exc (hartree) 163 exc(ipt)=(aa+cc*rs)*logrs+dd*rs+bb-efac*rsm1 164 ! compute potential vxc=d(rho*exc)/d(rho) (hartree) 165 vxc(ipt)=(aa+two_thirds*cc*rs)*logrs+(dd+dd-cc)*rs*third+& 166 & (bb-aa*third)-vfac*rsm1 167 ! compute d(vxc)/d(rho) (hartree*bohr^3) 168 else 169 t1=b1*sqrt(rs) 170 t2=b2*rs 171 den=1._dp/(1._dp+t1+t2) 172 exc(ipt)=ga*den-efac*rsm1 173 vxc(ipt)=ga*(1._dp+c7_6*t1+c4_3*t2)*den**2-vfac*rsm1 174 end if 175 end do 176 end if 177 ! 178 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}$.
SOURCE
236 subroutine xcspol(exc,npts,nspden,order,rspts,vxc,zeta,ndvxc,& !Mandatory arguments 237 & dvxc) !Optional arguments 238 239 !Arguments ------------------------------------ 240 !scalars 241 integer,intent(in) :: ndvxc,npts,nspden,order 242 !arrays 243 real(dp),intent(in) :: rspts(npts),zeta(npts) 244 real(dp),intent(out) :: exc(npts),vxc(npts,nspden) 245 real(dp),intent(out),optional :: dvxc(npts,ndvxc) 246 247 !Local variables------------------------------- 248 !The generation of density from rs needs rsfac and rsfac^(-3) : 249 !rsfac=(3/(4 Pi))^(1/3) ; rsfacm3=4pi/3 250 !Mike Teter s parameters of 8 April 1993. 251 !New parameters which accomodate spin polarization (fit to P-W) 252 !Paramagnetic limit:a0p,...b4p 253 !(a0=(3/4)(3/(2 Pi))^(2/3)) (note that b1=1 is fixed) 254 !Differences, ferromagnetic - paramagnetic (delta params):da1,da2,da3,db1,db2,db3,db4 255 !scalars 256 integer :: ipts 257 real(dp),parameter :: a0p=.4581652932831429_dp,a1p=2.217058676663745_dp 258 real(dp),parameter :: a2p=0.7405551735357053_dp,a3p=0.01968227878617998_dp 259 real(dp),parameter :: alpha_zeta=one-1.0d-6,b1p=one,b2p=4.504130959426697_dp 260 real(dp),parameter :: b3p=1.110667363742916_dp,b4p=0.02359291751427506_dp 261 real(dp),parameter :: da0=.119086804055547_dp,da1=0.6157402568883345_dp 262 real(dp),parameter :: da2=0.1574201515892867_dp,da3=0.003532336663397157_dp 263 real(dp),parameter :: db1=zero,db2=0.2673612973836267_dp 264 real(dp),parameter :: db3=0.2052004607777787_dp,db4=0.004200005045691381_dp 265 real(dp),parameter :: ft=4._dp/3._dp,rsfac=0.6203504908994000_dp 266 real(dp),parameter :: rsfacm3=rsfac**(-3) 267 real(dp) :: a0,a1,a2,a3,b1,b2,b3,b4,d1,d1m1,d2d1drs2,d2d1drsdf,d2excdf2 268 real(dp) :: d2excdrs2,d2excdrsdf,d2excdz2,d2fxcdz2,d2n1drs2,d2n1drsdf,dd1df 269 real(dp) :: dd1drs,dexcdf,dexcdrs,dexcdz,dfxcdz,dn1df,dn1drs,dvxcdrs 270 real(dp) :: dvxcpdrho,dvxcpdz,excipt,fact,fxc,n1 271 real(dp) :: rhom1,rs,vxcp,zet,zetm,zetm_third 272 real(dp) :: zetp,zetp_third 273 character(len=500) :: message 274 !no_abirules 275 !Set a minimum rho below which terms are 0 276 real(dp),parameter :: rhotol=1.d-28 277 !real(dp) :: delta,rho,rho_dn,rho_dnm,rho_dnp,rho_up,rho_upm,rho_upp,zeta_mean 278 279 ! ************************************************************************* 280 281 !Checks the compatibility between the presence of dvxc and ndvxc 282 if(ndvxc /=0 .neqv. present(dvxc))then 283 message = 'If ndvxc/=0 there must be the optional argument dvxc' 284 ABI_BUG(message) 285 end if 286 287 !Checks the compatibility between the inputs and the presence of the optional arguments 288 if(abs(order) <= 1 .and. ndvxc /= 0)then 289 write(message, '(4a,i0)' )ch10,& 290 & 'The order chosen does not need the presence',ch10,& 291 & 'of the vector dvxc, that is needed only with |order|>1 , while we have',order 292 ABI_BUG(message) 293 end if 294 295 if(nspden == 1 .and. ndvxc /=0 .and. ndvxc /= 2)then 296 write(message,'(a,i0)')' Once nspden=1 we must have ndvxc=2, while we have',ndvxc 297 ABI_BUG(message) 298 end if 299 300 if(nspden == 2 .and. ndvxc /=0 .and. ndvxc /= 3)then 301 write(message, '(a,i0)' )' Once nspden=2 we must have ndvxc=3, while we have',ndvxc 302 ABI_BUG(message) 303 end if 304 305 306 !Although fact is parameter value, some compilers are not able to evaluate 307 !it at compile time. 308 fact=one/(two**(four*third)-two) 309 310 !DEBUG 311 !Finite-difference debugging, do not take away 312 !debug=1 313 !zeta_mean=0.1_dp 314 !delta=0.0001 315 !if(debug==1)then 316 !do ipts=1,npts,5 317 !rho=ipts*0.01_dp 318 !rho_up=rho*(one+zeta_mean)*half 319 !rho_dn=rho*(one-zeta_mean)*half 320 !rho_upp=rho_up+delta 321 !rho_upm=rho_up-delta 322 !rho_dnp=rho_dn+delta 323 !rho_dnm=rho_dn-delta 324 !First possibility : vary rho up , and then rho down 325 !zeta(ipts )=(rho_up -rho_dn )/(rho_up +rho_dn ) 326 !zeta(ipts+1)=(rho_upp-rho_dn )/(rho_upp+rho_dn ) 327 !zeta(ipts+2)=(rho_upm-rho_dn )/(rho_upm+rho_dn ) 328 !zeta(ipts+3)=(rho_up -rho_dnp)/(rho_up +rho_dnp) 329 !zeta(ipts+4)=(rho_up -rho_dnm)/(rho_up +rho_dnm) 330 !rspts(ipts )=rsfac*(rho_up +rho_dn )**(-third) 331 !rspts(ipts+1)=rsfac*(rho_upp+rho_dn )**(-third) 332 !rspts(ipts+2)=rsfac*(rho_upm+rho_dn )**(-third) 333 !rspts(ipts+3)=rsfac*(rho_up +rho_dnp)**(-third) 334 !rspts(ipts+4)=rsfac*(rho_up +rho_dnm)**(-third) 335 !DEBUGBUG : another possibility : vary rho and zeta 336 !zeta(ipts+1)=zeta(ipts ) 337 !zeta(ipts+2)=zeta(ipts ) 338 !zeta(ipts+3)=zeta(ipts )+delta 339 !zeta(ipts+4)=zeta(ipts )-delta 340 !rspts(ipts+1)=rsfac*(rho+delta)**(-third) 341 !rspts(ipts+2)=rsfac*(rho-delta )**(-third) 342 !rspts(ipts+3)=rspts(ipts ) 343 !rspts(ipts+4)=rspts(ipts ) 344 !ENDDEBUGBUG 345 !end do 346 !end if 347 !nspden=2 348 !order=2 349 !ENDDEBUG 350 351 if (nspden==1) then 352 ! separate cases with respect to order 353 if(order==-2) then 354 ! No spin-polarization so skip steps related to zeta not 0 355 do ipts=1,npts 356 357 rs=rspts(ipts) 358 n1=a0p+rs*(a1p+rs*(a2p+rs*a3p)) 359 d1=rs*(b1p+rs*(b2p+rs*(b3p+rs*b4p))) 360 d1m1=one/d1 361 362 ! Exchange-correlation energy 363 excipt=-n1*d1m1 364 exc(ipts)=excipt 365 366 ! Exchange-correlation potential 367 dn1drs=a1p+rs*(2._dp*a2p+rs*(3._dp*a3p)) 368 dd1drs=b1p+rs*(2._dp*b2p+rs*(3._dp*b3p+rs*(4._dp*b4p))) 369 370 ! dexcdrs is d(exc)/d(rs) 371 dexcdrs=-(dn1drs+excipt*dd1drs)*d1m1 372 vxc(ipts,1)=excipt-third*rs*dexcdrs 373 374 ! If the exchange-correlation kernel is needed 375 376 d2n1drs2=2._dp*a2p+rs*(6._dp*a3p) 377 d2d1drs2=2._dp*b2p+rs*(6._dp*b3p+rs*(12._dp*b4p)) 378 ! d2excdrs2 is d2(exc)/d(rs)2 379 d2excdrs2=-(d2n1drs2+2._dp*dexcdrs*dd1drs+excipt*d2d1drs2)*d1m1 380 dvxcdrs=third*(2.0_dp*dexcdrs-rs*d2excdrs2) 381 ! And d(vxc)/d(rho)=(-rs/(3*rho))*d(vxc)/d(rs) 382 dvxc(ipts,1)= -rs**4*rsfacm3*third*dvxcdrs 383 384 ! dn1df=d(n1)/d(fxc) and dd1df=d(d1)/d(fxc) 385 dn1df=da0+rs*(da1+rs*(da2+rs*da3)) 386 dd1df=rs*(db1+rs*(db2+rs*(db3+rs*db4))) 387 dexcdf=-(dn1df+excipt*dd1df)*d1m1 388 ! d2(fxc)/d(zeta)2 389 d2fxcdz2=ft*third*(alpha_zeta**2)*2._dp*fact 390 ! d2(exc)/d(zeta)2 391 d2excdz2=d2fxcdz2*dexcdf 392 rhom1=rsfacm3*rs**3 393 dvxc(ipts,2)= dvxc(ipts,1) - d2excdz2*rhom1 394 end do 395 else if(order**2>1) then 396 ! No spin-polarization so skip steps related to zeta not 0 397 do ipts=1,npts 398 399 rs=rspts(ipts) 400 n1=a0p+rs*(a1p+rs*(a2p+rs*a3p)) 401 d1=rs*(b1p+rs*(b2p+rs*(b3p+rs*b4p))) 402 d1m1=one/d1 403 404 ! Exchange-correlation energy 405 excipt=-n1*d1m1 406 exc(ipts)=excipt 407 408 ! Exchange-correlation potential 409 dn1drs=a1p+rs*(2._dp*a2p+rs*(3._dp*a3p)) 410 dd1drs=b1p+rs*(2._dp*b2p+rs*(3._dp*b3p+rs*(4._dp*b4p))) 411 412 ! dexcdrs is d(exc)/d(rs) 413 dexcdrs=-(dn1drs+excipt*dd1drs)*d1m1 414 vxc(ipts,1)=excipt-third*rs*dexcdrs 415 416 ! If the exchange-correlation kernel is needed 417 d2n1drs2=2._dp*a2p+rs*(6._dp*a3p) 418 d2d1drs2=2._dp*b2p+rs*(6._dp*b3p+rs*(12._dp*b4p)) 419 ! d2excdrs2 is d2(exc)/d(rs)2 420 d2excdrs2=-(d2n1drs2+2._dp*dexcdrs*dd1drs+excipt*d2d1drs2)*d1m1 421 dvxcdrs=third*(2.0_dp*dexcdrs-rs*d2excdrs2) 422 ! And d(vxc)/d(rho)=(-rs/(3*rho))*d(vxc)/d(rs) 423 dvxc(ipts,1)= -rs**4*rsfacm3*third*dvxcdrs 424 425 end do 426 else 427 ! No spin-polarization so skip steps related to zeta not 0 428 do ipts=1,npts 429 430 rs=rspts(ipts) 431 n1=a0p+rs*(a1p+rs*(a2p+rs*a3p)) 432 d1=rs*(b1p+rs*(b2p+rs*(b3p+rs*b4p))) 433 d1m1=one/d1 434 435 ! Exchange-correlation energy 436 excipt=-n1*d1m1 437 exc(ipts)=excipt 438 439 ! Exchange-correlation potential 440 dn1drs=a1p+rs*(2._dp*a2p+rs*(3._dp*a3p)) 441 dd1drs=b1p+rs*(2._dp*b2p+rs*(3._dp*b3p+rs*(4._dp*b4p))) 442 443 ! dexcdrs is d(exc)/d(rs) 444 dexcdrs=-(dn1drs+excipt*dd1drs)*d1m1 445 vxc(ipts,1)=excipt-third*rs*dexcdrs 446 end do 447 448 end if 449 450 451 ! Allows for nspden==1, in the case of testing nspden=1 against nspden=2 452 else if (nspden<=2) then 453 454 455 ! DEBUG 456 ! do not take away : allows to compare nspden=1 and nspden=2 coding 457 ! if (nspden==1)then 458 ! zeta(:)=zero 459 ! end if 460 ! ENDDEBUG 461 ! separate cases with respect to order 462 if(abs(order)>1) then 463 ! Allow for spin polarization. This part could be optimized for speed. 464 do ipts=1,npts 465 466 rs=rspts(ipts) 467 zet=zeta(ipts) 468 zetp=one+zet*alpha_zeta 469 zetm=one-zet*alpha_zeta 470 zetp_third=zetp**third 471 zetm_third=zetm**third 472 ! Exchange energy spin interpolation function f(zeta) 473 fxc=( zetp*zetp_third + zetm*zetm_third - two ) *fact 474 475 a0=a0p+fxc*da0 476 a1=a1p+fxc*da1 477 a2=a2p+fxc*da2 478 a3=a3p+fxc*da3 479 b1=b1p+fxc*db1 480 b2=b2p+fxc*db2 481 b3=b3p+fxc*db3 482 b4=b4p+fxc*db4 483 484 n1= a0+rs*(a1+rs*(a2+rs*a3)) 485 d1=rs*(b1+rs*(b2+rs*(b3+rs*b4))) 486 d1m1=one/d1 487 488 ! Exchange-correlation energy 489 excipt=-n1*d1m1 490 exc(ipts)=excipt 491 492 ! Exchange-correlation potential 493 dn1drs=a1+rs*(2._dp*a2+rs*(3._dp*a3)) 494 dd1drs=b1+rs*(2._dp*b2+rs*(3._dp*b3+rs*(4._dp*b4))) 495 ! dexcdrs is d(exc)/d(rs) 496 dexcdrs=-(dn1drs+excipt*dd1drs)*d1m1 497 498 ! Only vxcp contributes when paramagnetic 499 vxcp=excipt-third*rs*dexcdrs 500 501 ! d(fxc)/d(zeta) (which is 0 at zeta=0) 502 dfxcdz=ft*alpha_zeta*(zetp_third-zetm_third)*fact 503 504 ! dn1df=d(n1)/d(fxc) and dd1df=d(d1)/d(fxc) 505 dn1df=da0+rs*(da1+rs*(da2+rs*da3)) 506 dd1df=rs*(db1+rs*(db2+rs*(db3+rs*db4))) 507 508 ! dexcdz is d(exc)/d(zeta) 509 dexcdf=-(dn1df+excipt*dd1df)*d1m1 510 dexcdz=dfxcdz*dexcdf 511 512 ! Compute Vxc for both spin channels 513 514 vxc(ipts,1)=vxcp - (zet-one)*dexcdz 515 vxc(ipts,2)=vxcp - (zet+one)*dexcdz 516 517 ! DEBUG Allow to check the variation of rho and zeta 518 ! vxc(ipts,1)=vxcp 519 ! vxc(ipts,2)=dexcdz 520 ! ENDDEBUG 521 ! Compute second derivative with respect to rho 522 d2n1drs2=2._dp*a2+rs*(6._dp*a3) 523 d2d1drs2=2._dp*b2+rs*(6._dp*b3+rs*(12._dp*b4)) 524 ! d2excdrs2 is d2(exc)/d(rs)2 525 d2excdrs2=-(d2n1drs2+two*dexcdrs*dd1drs+excipt*d2d1drs2)*d1m1 526 dvxcdrs=third*(two*dexcdrs-rs*d2excdrs2) 527 ! And d(vxc)/d(rho) paramagnetic =(-rs/(3*rho))*d(vxcp)/d(rs) 528 ! remember : 1/rho=(4pi/3)*rs**3=rsfacm3*rs**3 529 rhom1=rsfacm3*rs**3 530 dvxcpdrho= -rs*rhom1*third * dvxcdrs 531 532 ! Compute mixed second derivative with respect to rho and zeta 533 d2n1drsdf=da1+rs*(2._dp*da2+rs*(3._dp*da3)) 534 d2d1drsdf=db1+rs*(2._dp*db2+rs*(3._dp*db3+rs*(4._dp*db4))) 535 ! d2excdrsdf is d2(exc)/d(rs)df 536 d2excdrsdf=-(d2n1drsdf+dexcdrs*dd1df+dexcdf*dd1drs+excipt*d2d1drsdf)*d1m1 537 ! d(vxc)/d(zeta) paramagnetic 538 dvxcpdz=dexcdz-third*rs*dfxcdz*d2excdrsdf 539 540 ! Compute second derivative with respect to zeta 541 ! the second derivative of n1 and d1 wrt f vanishes 542 d2excdf2=-(two*dexcdf*dd1df)*d1m1 543 ! d2(fxc)/d(zeta)2 544 d2fxcdz2=ft*third*(alpha_zeta**2)*(zetp_third**(-2)+zetm_third**(-2))*fact 545 ! d2(exc)/d(zeta)2 546 d2excdz2=d2fxcdz2*dexcdf+dfxcdz**2*d2excdf2 547 548 ! Compute now the three second derivatives of the Exc energy with respect 549 ! to : wrt twice spin-up ; wrt spin-up and spin-dn ; wrt twice spin-down 550 dvxc(ipts,1)= dvxcpdrho & 551 & +two*rhom1*( one-zet)*(dvxcpdz-dexcdz) & 552 & +d2excdz2*rhom1*(one-zet)**2 553 dvxc(ipts,2)= dvxcpdrho & 554 & +two*rhom1*( -zet)*(dvxcpdz-dexcdz) & 555 & +d2excdz2*rhom1*(one-zet)*(-one-zet) 556 ! if(nspden==2)then 557 dvxc(ipts,3)= dvxcpdrho & 558 & +two*rhom1*(-one-zet)*(dvxcpdz-dexcdz) & 559 & +d2excdz2*rhom1*(-one-zet)**2 560 ! else 561 ! ! For testing purposes, need the spin-averaged quantity 562 ! dvxc(ipts,1)= ( dvxc(ipts,1) + dvxc(ipts,2) ) * half 563 ! end if 564 565 ! DEBUG Allow to check the variation of rho and zeta 566 ! dvxc(ipts,1)=dvxcpdrho 567 ! dvxc(ipts,2)=d2excdz2 568 ! dvxc(ipts,3)=dvxcpdz 569 ! ENDDEBUG 570 end do 571 else 572 ! Allow for spin polarization. This part could be optimized for speed. 573 do ipts=1,npts 574 575 rs=rspts(ipts) 576 zet=zeta(ipts) 577 zetp=one+zet*alpha_zeta 578 zetm=one-zet*alpha_zeta 579 zetp_third=zetp**third 580 zetm_third=zetm**third 581 ! Exchange energy spin interpolation function f(zeta) 582 fxc=( zetp*zetp_third + zetm*zetm_third - two ) *fact 583 584 a0=a0p+fxc*da0 585 a1=a1p+fxc*da1 586 a2=a2p+fxc*da2 587 a3=a3p+fxc*da3 588 b1=b1p+fxc*db1 589 b2=b2p+fxc*db2 590 b3=b3p+fxc*db3 591 b4=b4p+fxc*db4 592 593 n1= a0+rs*(a1+rs*(a2+rs*a3)) 594 d1=rs*(b1+rs*(b2+rs*(b3+rs*b4))) 595 d1m1=one/d1 596 597 ! Exchange-correlation energy 598 excipt=-n1*d1m1 599 exc(ipts)=excipt 600 601 ! Exchange-correlation potential 602 dn1drs=a1+rs*(2._dp*a2+rs*(3._dp*a3)) 603 dd1drs=b1+rs*(2._dp*b2+rs*(3._dp*b3+rs*(4._dp*b4))) 604 ! dexcdrs is d(exc)/d(rs) 605 dexcdrs=-(dn1drs+excipt*dd1drs)*d1m1 606 607 ! Only vxcp contributes when paramagnetic 608 vxcp=excipt-third*rs*dexcdrs 609 610 ! d(fxc)/d(zeta) (which is 0 at zeta=0) 611 dfxcdz=ft*alpha_zeta*(zetp_third-zetm_third)*fact 612 613 ! dn1df=d(n1)/d(fxc) and dd1df=d(d1)/d(fxc) 614 dn1df=da0+rs*(da1+rs*(da2+rs*da3)) 615 dd1df=rs*(db1+rs*(db2+rs*(db3+rs*db4))) 616 617 ! dexcdz is d(exc)/d(zeta) 618 dexcdf=-(dn1df+excipt*dd1df)*d1m1 619 dexcdz=dfxcdz*dexcdf 620 621 ! Compute Vxc for both spin channels 622 623 vxc(ipts,1)=vxcp - (zet-one)*dexcdz 624 vxc(ipts,2)=vxcp - (zet+one)*dexcdz 625 626 ! DEBUG Allow to check the variation of rho and zeta 627 ! vxc(ipts,1)=vxcp 628 ! vxc(ipts,2)=dexcdz 629 ! ENDDEBUG 630 end do 631 end if 632 else 633 634 ! Disallowed value for nspden 635 write(message, '(3a,i0)' )& 636 & ' Argument nspden must be 1 or 2; ',ch10,& 637 & ' Value provided as argument was ',nspden 638 ABI_BUG(message) 639 end if 640 641 !DEBUG 642 !Finite-difference debugging, do not take away 643 !if(debug==1)then 644 !write(std_out,*)' delta =',delta 645 !do ipts=1,npts,5 646 !rho=(rspts(ipts)/rsfac)**(-3) 647 !write(std_out,'(a,i5,a,2es16.8)' ) ' Point number',ipts,' with rho,zeta=',rho,zeta(ipts) 648 !write(std_out,'(3es16.8)' )exc(ipts)*rho,vxc(ipts,1),vxc(ipts,2) 649 !write(std_out,'(3es16.8)' )dvxc(ipts,1),dvxc(ipts,3),dvxc(ipts,2) 650 !write(std_out,'(3es16.8)' )exc(ipts)*rho,& 651 !& ( exc(ipts+1)*(rho+delta) - exc(ipts+2)*(rho-delta) )/2._dp/delta,& 652 !& ( exc(ipts+3)*(rho+delta) - exc(ipts+4)*(rho-delta) )/2._dp/delta 653 !write(std_out,'(4es16.8)' )& 654 !& ( vxc(ipts+1,1) - vxc(ipts+2,1) )/2._dp/delta,& 655 !& ( vxc(ipts+3,2) - vxc(ipts+4,2) )/2._dp/delta,& 656 !& ( vxc(ipts+3,1) - vxc(ipts+4,1) )/2._dp/delta,& 657 !& ( vxc(ipts+1,2) - vxc(ipts+2,2) )/2._dp/delta 658 !end do 659 !stop 660 !end if 661 !ENDDEBUG 662 663 !DEBUG 664 !if(order==-2)then 665 !write(std_out,*)' xcspol : ipts,npts ',ipts,npts 666 !write(std_out,*)dvxcdrs,d2excdz2,d2fxcdz2,dexcdf 667 !write(std_out,*)rhom1 668 !write(std_out,*)dvxc(1000,1),dvxc(1000,2) 669 !stop 670 !end if 671 !ENDDEBUG 672 673 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)$.
SOURCE
726 subroutine xctetr(exc,npt,order,rhor,rspts,vxc,& !Mandatory arguments 727 & d2vxc,dvxc) !Optional arguments 728 729 !Arguments ------------------------------------ 730 !scalars 731 integer,intent(in) :: npt,order 732 !arrays 733 real(dp),intent(in) :: rhor(npt),rspts(npt) 734 real(dp),intent(out) :: exc(npt),vxc(npt) 735 real(dp),intent(out),optional :: d2vxc(npt),dvxc(npt) 736 737 !Local variables------------------------------- 738 !rsfac=(3/(4 Pi))^(1/3) 739 !Mike Teter s parameters: (keep 8 digits after decimal) 740 !(a0=(3/4)(3/(2 Pi))^(2/3) 741 !scalars 742 integer :: ipt 743 real(dp),parameter :: a0=.4581652932831429_dp,a1=2.40875407_dp,a2=.88642404_dp 744 real(dp),parameter :: a3=.02600342_dp,b1=1.0_dp,b2=4.91962865_dp 745 real(dp),parameter :: b3=1.34799453_dp,b4=.03120453_dp,c1=4._dp*a0*b1/3.0_dp 746 real(dp),parameter :: c2=5.0_dp*a0*b2/3.0_dp+a1*b1 747 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 748 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 749 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 750 real(dp),parameter :: c6=5.0_dp*a2*b4/3.0_dp+a3*b3,c7=4.0_dp*a3*b4/3.0_dp 751 real(dp),parameter :: rsfac=0.6203504908994000_dp 752 real(dp) :: d1,d1m1,d2vxcr,d3,d4,dvxcdr,n1,n2,n3,n4,rhom1,rs 753 character(len=500) :: message 754 755 ! ************************************************************************* 756 ! 757 !Checks the values of order 758 if(order<0 .or. order>3)then 759 write(message, '(a,a,a,i6)' )& 760 & 'With Teter 91 Ceperley-Alder xc functional, the only',ch10,& 761 & 'allowed values for order are 0, 1, 2 or 3, while it is found to be',order 762 ABI_BUG(message) 763 end if 764 765 !Checks the compatibility between the order and the presence of the optional arguments 766 if(order /=3 .and. present(d2vxc))then 767 write(message, '(a,a,a,i6)' )& 768 & 'The order chosen does not need the presence',ch10,& 769 & 'of the vector d2vxc, that is needed only with order=3, while we have',order 770 ABI_BUG(message) 771 end if 772 773 if(order <= 1 .and. present(dvxc))then 774 write(message, '(a,a,a,i6)' )& 775 & 'The order chosen does not need the presence',ch10,& 776 & 'of the vector dvxc, that is needed with order > 1, while we have',order 777 ABI_BUG(message) 778 end if 779 780 !separated cases with respect to order 781 782 if (order<=1) then 783 ! Loop over grid points 784 do ipt=1,npt 785 rs=rspts(ipt) 786 n1=-(a0+rs*(a1+rs*(a2+rs*a3))) 787 d1=rs*(b1+rs*(b2+rs*(b3+rs*b4))) 788 d1m1=1.0_dp/d1 789 n2=-rs*(c1+rs*(c2+rs*(c3+rs*(c4+rs*(c5+rs*(c6+rs*c7)))))) 790 ! 791 ! Exchange-correlation energy 792 exc(ipt)=n1*d1m1 793 ! 794 ! Exchange-correlation potential 795 vxc(ipt)=n2*d1m1**2 796 end do 797 else if (order>2) then 798 ! Loop over grid points 799 do ipt=1,npt 800 rs=rspts(ipt) 801 n1=-(a0+rs*(a1+rs*(a2+rs*a3))) 802 d1=rs*(b1+rs*(b2+rs*(b3+rs*b4))) 803 d1m1=1.0_dp/d1 804 n2=-rs*(c1+rs*(c2+rs*(c3+rs*(c4+rs*(c5+rs*(c6+rs*c7)))))) 805 ! 806 ! Exchange-correlation energy 807 exc(ipt)=n1*d1m1 808 ! 809 ! Exchange-correlation potential 810 vxc(ipt)=n2*d1m1**2 811 ! Assemble derivative of vxc wrt rs 812 n3=-(c1+rs*(2._dp*c2+rs*(3._dp*c3+rs*(4._dp*c4+rs*(5._dp*c5+& 813 & rs*(6._dp*c6+rs*(7._dp*c7))))))) 814 d3=b1+rs*(2._dp*b2+rs*(3._dp*b3+rs*(4._dp*b4))) 815 dvxcdr=(n3-d3*(2._dp*n2*d1m1))*d1m1**2 816 rhom1=1.0_dp/rhor(ipt) 817 ! 818 ! derivative of vxc wrt rho 819 dvxc(ipt)=-dvxcdr*rs*third*rhom1 820 ! 821 822 ! Assemble derivative d^2(Vxc)/d(rs)^2 823 n4=-(2.0_dp*c2+rs*(6.0_dp*c3+rs*(12.0_dp*c4+rs*(20.0_dp*c5+& 824 & rs*(30.0_dp*c6+rs*(42.0_dp*c7)))))) 825 d4=2.0_dp*b2+rs*(6.0_dp*b3+rs*(12.0_dp*b4)) 826 d2vxcr=(n4-2.0_dp*(2.0_dp*n3*d3+n2*d4-3.0_dp*n2*d3**2*d1m1)*d1m1)*d1m1**2 827 828 ! Derivative d^2(Vxc)/d(rho)^2 829 d2vxc(ipt)=(rs*third*rhom1)*(4.0_dp*dvxcdr+rs*d2vxcr)*third*rhom1 830 831 end do 832 else if (order>1) then 833 ! Loop over grid points 834 do ipt=1,npt 835 rs=rspts(ipt) 836 n1=-(a0+rs*(a1+rs*(a2+rs*a3))) 837 d1=rs*(b1+rs*(b2+rs*(b3+rs*b4))) 838 d1m1=1.0_dp/d1 839 n2=-rs*(c1+rs*(c2+rs*(c3+rs*(c4+rs*(c5+rs*(c6+rs*c7)))))) 840 ! 841 ! Exchange-correlation energy 842 exc(ipt)=n1*d1m1 843 ! 844 ! Exchange-correlation potential 845 vxc(ipt)=n2*d1m1**2 846 ! Assemble derivative of vxc wrt rs 847 n3=-(c1+rs*(2._dp*c2+rs*(3._dp*c3+rs*(4._dp*c4+rs*(5._dp*c5+& 848 & rs*(6._dp*c6+rs*(7._dp*c7))))))) 849 d3=b1+rs*(2._dp*b2+rs*(3._dp*b3+rs*(4._dp*b4))) 850 dvxcdr=(n3-d3*(2._dp*n2*d1m1))*d1m1**2 851 rhom1=1.0_dp/rhor(ipt) 852 ! 853 ! derivative of vxc wrt rho 854 dvxc(ipt)=-dvxcdr*rs*third*rhom1 855 ! 856 end do 857 end if 858 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) npts= number of points to be computed nspden=number if spin density component (necessarily 1 here) 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), 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
SOURCE
1252 subroutine xctfw(temp,exci,fxci,usefxc,rho_updn,vxci,npts,nspden,dvxcdgr,ndvxcdgr,grho2_updn) 1253 1254 !Arguments ------------------------------------ 1255 !scalars 1256 integer,intent(in) :: ndvxcdgr,npts,nspden,usefxc 1257 real(dp),intent(in) :: temp 1258 !arrays 1259 real(dp),intent(in) :: grho2_updn(npts,2*nspden-1),rho_updn(npts,nspden) 1260 real(dp),intent(inout) :: dvxcdgr(npts,ndvxcdgr),fxci(npts*usefxc),exci(npts),vxci(npts,nspden) 1261 1262 !Local variables------------------------------- 1263 !scalars 1264 integer :: iperrot,ipts 1265 logical :: has_fxc,has_dvxcdgr 1266 real(dp) :: etfw,rho,rho_inv,rhomot,yperrot0,vtfw 1267 real(dp) :: yperrot,uperrot,dyperrotdn,duperrotdyperrot 1268 real(dp) :: hperrot,dhperrotdyperrot,dhperrotdn,dhperrotduperrot 1269 !arrays 1270 real(dp) :: wpy(0:7), wpu(0:7) 1271 real(dp),allocatable :: rho_updnm1_3(:,:) 1272 1273 ! ************************************************************************* 1274 1275 has_fxc=(usefxc/=0) 1276 has_dvxcdgr=(ndvxcdgr/=0) 1277 1278 yperrot0=1.666081101_dp 1279 1280 wpy(0)=0.5_dp; wpy(1)=-0.1999176316_dp 1281 wpy(2)=0.09765615709_dp; wpy(3)=-0.06237609924_dp 1282 wpy(4)=0.05801466322_dp; wpy(5)=-0.04449287774_dp 1283 wpy(6)=0.01903211697_dp; wpy(7)=-0.003284096926_dp 1284 1285 wpu(0)=one/6._dp; wpu(1)=0.311590799_dp 1286 wpu(2)=3.295662439_dp; wpu(3)=-29.22038326_dp 1287 wpu(4)=116.1084531_dp; wpu(5)=-250.4543147_dp 1288 wpu(6)=281.433688_dp; wpu(7)=-128.8784806_dp 1289 1290 ABI_MALLOC(rho_updnm1_3,(npts,2)) 1291 1292 call invcb(rho_updn(:,1),rho_updnm1_3(:,1),npts) 1293 1294 do ipts=1,npts 1295 1296 rho =rho_updn(ipts,1) 1297 rhomot=rho_updnm1_3(ipts,1) 1298 rho_inv=rhomot*rhomot*rhomot 1299 1300 yperrot=pi*pi/sqrt2/temp**1.5*two*rho 1301 uperrot=yperrot**(2./3.) 1302 1303 dyperrotdn=pi*pi/sqrt2/temp**1.5*2.0_dp 1304 1305 hperrot=zero 1306 dhperrotdyperrot=zero 1307 dhperrotduperrot=zero 1308 if(yperrot<=yperrot0)then 1309 do iperrot=0,7 1310 hperrot=hperrot+wpy(iperrot)*yperrot**iperrot 1311 dhperrotdyperrot=dhperrotdyperrot+iperrot*wpy(iperrot)*yperrot**(iperrot-1) 1312 end do 1313 hperrot=one/12.0_dp*hperrot 1314 dhperrotdyperrot=one/12.0_dp*dhperrotdyperrot 1315 dhperrotdn=dhperrotdyperrot*dyperrotdn 1316 else 1317 do iperrot=0,7 1318 hperrot=hperrot+wpu(iperrot)/uperrot**(2*iperrot) 1319 dhperrotduperrot=dhperrotduperrot-2.*iperrot*wpu(iperrot)/uperrot**(2*iperrot+1) 1320 end do 1321 hperrot=one/12.0_dp*hperrot 1322 dhperrotduperrot=one/12.0_dp*dhperrotduperrot 1323 duperrotdyperrot=two/3._dp/yperrot**(1./3.) 1324 dhperrotdn=dhperrotduperrot*duperrotdyperrot*dyperrotdn 1325 end if 1326 1327 etfw=hperrot*grho2_updn(ipts,1)*rho_inv*rho_inv 1328 vtfw=-etfw + rho/hperrot*dhperrotdn*etfw 1329 1330 if(yperrot<=yperrot0)then 1331 exci(ipts) = exci(ipts) + etfw + 1.5_dp*yperrot*dhperrotdyperrot*grho2_updn(ipts,1)*rho_inv*rho_inv 1332 else 1333 exci(ipts) = exci(ipts) + etfw + uperrot*dhperrotduperrot*grho2_updn(ipts,1)*rho_inv*rho_inv 1334 end if 1335 vxci(ipts,1) = vxci(ipts,1) + vtfw 1336 if (has_fxc) fxci(ipts) = fxci(ipts) + etfw 1337 if (has_dvxcdgr) dvxcdgr(ipts,1)= dvxcdgr(ipts,1)+two*hperrot*rho_inv 1338 1339 end do 1340 1341 ABI_FREE(rho_updnm1_3) 1342 1343 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)
SOURCE
884 subroutine xcwign(exc,npt,order,rspts,vxc,& !Mandatory arguments 885 & dvxc) !Optional arguments 886 887 !Arguments ------------------------------------ 888 !scalars 889 integer,intent(in) :: npt,order 890 !arrays 891 real(dp),intent(in) :: rspts(npt) 892 real(dp),intent(out) :: exc(npt),vxc(npt) 893 real(dp),intent(out),optional :: dvxc(npt) 894 895 !Local variables------------------------------- 896 !c1 and c2 are the Wigner parameters in hartree and bohr resp. 897 !scalars 898 integer :: ipt 899 real(dp),parameter :: c1=0.44_dp,c2=7.8_dp,c4_3=4.0_dp/3.0_dp 900 real(dp),parameter :: c8_27=8.0_dp/27.0_dp 901 real(dp) :: dfac,efac,rs,rsc2m1,rsm1,vfac,vxcnum 902 character(len=500) :: message 903 904 ! ************************************************************************* 905 906 !Checks the values of order 907 if(order<0 .or. order>2)then 908 write(message, '(a,a,a,i0)' )& 909 & 'With Wigner xc functional, the only',ch10,& 910 & 'allowed values for order are 0, 1 or 2, while it is found to be',order 911 ABI_BUG(message) 912 end if 913 914 !Checks the compatibility between the order and the presence of the optional arguments 915 if(order <= 1 .and. present(dvxc))then 916 write(message, '(a,a,a,i3)' )& 917 & 'The order chosen does not need the presence',ch10,& 918 & 'of the vector dvxc, that is needed only with order=2 , while we have',order 919 ABI_BUG(message) 920 end if 921 922 !Compute vfac=(3/(2*Pi))^(2/3) 923 vfac=(1.5_dp/pi)**(2.0_dp/3.0_dp) 924 !Compute efac=(3/4)*vfac 925 efac=0.75_dp*vfac 926 !Compute dfac=(4*Pi/9)*vfac 927 dfac=(4.0_dp*pi/9.0_dp)*vfac 928 929 !separate cases with respect to order 930 if (order==2) then 931 932 ! Loop over grid points 933 do ipt=1,npt 934 rs=rspts(ipt) 935 rsm1=1.0_dp/rs 936 rsc2m1=1.0_dp/(rs+c2) 937 ! compute energy density (hartree) 938 exc(ipt)=-c1*rsc2m1-efac*rsm1 939 vxcnum=-(c4_3*rs+c2)*c1 940 ! compute potential (hartree) 941 vxc(ipt)=vxcnum*rsc2m1**2-vfac*rsm1 942 ! compute d(vxc)/d(rho) (hartree*bohr^3) 943 dvxc(ipt)=-(c8_27*pi)*(c1*rs**4)*(rs+rs+c2)*rsc2m1**3-dfac*rs**2 944 end do 945 else 946 947 ! Loop over grid points 948 do ipt=1,npt 949 rs=rspts(ipt) 950 rsm1=1.0_dp/rs 951 rsc2m1=1.0_dp/(rs+c2) 952 ! compute energy density (hartree) 953 exc(ipt)=-c1*rsc2m1-efac*rsm1 954 vxcnum=-(c4_3*rs+c2)*c1 955 ! compute potential (hartree) 956 vxc(ipt)=vxcnum*rsc2m1**2-vfac*rsm1 957 end do 958 959 end if 960 ! 961 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)
SOURCE
1078 subroutine xcxalp(exc,npt,order,rspts,vxc, dvxc) ! dvxc is optional 1079 1080 !Arguments ------------------------------------ 1081 !scalars 1082 integer,intent(in) :: npt,order 1083 !arrays 1084 real(dp),intent(in) :: rspts(npt) 1085 real(dp),intent(out) :: exc(npt),vxc(npt) 1086 real(dp),intent(out),optional :: dvxc(npt) 1087 1088 !Local variables------------------------------- 1089 !Set value of alpha in "X-alpha" method 1090 !scalars 1091 integer :: ipt 1092 real(dp),parameter :: alpha=1.0_dp 1093 real(dp) :: dfac,efac,rs,rsm1,vfac 1094 character(len=500) :: message 1095 1096 ! ************************************************************************* 1097 1098 !Checks the values of order 1099 if(order<0 .or. order>2)then 1100 write(message, '(a,a,a,i3)' )& 1101 & 'With X-alpha xc functional, the only',ch10,& 1102 & 'allowed values for order are 0, 1 or 2, while it is found to be',order 1103 ABI_BUG(message) 1104 end if 1105 1106 !Compute vfac=(3/(2*Pi))^(2/3) 1107 vfac=(1.5_dp/pi)**(2.0_dp/3.0_dp) 1108 !Compute efac=(3/4)*vfac 1109 efac=0.75_dp*vfac 1110 !Compute dfac=(4*Pi/9)*vfac 1111 dfac=(4.0_dp*pi/9.0_dp)*vfac 1112 1113 !separate cases with respect to order 1114 if(order==2) then 1115 ! Loop over grid points 1116 do ipt=1,npt 1117 rs=rspts(ipt) 1118 rsm1=1.0_dp/rs 1119 ! compute energy density (hartree) 1120 exc(ipt)=-alpha*efac*rsm1 1121 ! compute potential (hartree) 1122 vxc(ipt)=-alpha*vfac*rsm1 1123 ! compute d(vxc)/d(rho) (hartree*bohr^3) 1124 dvxc(ipt)=-alpha*dfac*rs**2 1125 end do 1126 else 1127 ! Loop over grid points 1128 do ipt=1,npt 1129 rs=rspts(ipt) 1130 rsm1=1.0_dp/rs 1131 ! compute energy density (hartree) 1132 exc(ipt)=-alpha*efac*rsm1 1133 ! compute potential (hartree) 1134 vxc(ipt)=-alpha*vfac*rsm1 1135 end do 1136 end if 1137 ! 1138 end subroutine xcxalp