TABLE OF CONTENTS
ABINIT/m_xchcth [ Modules ]
NAME
m_xchcth
FUNCTION
Treat XC GGA functional of HCTH type. See Hamprecht, Cohen, Tozer and Handy, J. Chem. Phys. 109, 6264 (1998) [[cite:Hamprecht1998]] for HCTH-93 Boese, Doltsinis, Handy and Sprik, J. Chem. Phys. 112, 1670 (2000) [[cite:Boese2000]] for HCTH-120 and HCTH-147. Boese and Handy , J. Chem. Phys. 114, 5497 (2001) [[cite:Boese2001]] for HCTH-407.
COPYRIGHT
Copyright (C) 2002-2024 ABINIT group (XG,LG) 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
19 #if defined HAVE_CONFIG_H 20 #include "config.h" 21 #endif 22 23 #include "abi_common.h" 24 25 module m_xchcth 26 27 use defs_basis 28 use m_errors 29 use m_abicore 30 31 use m_numeric_tools, only : invcb 32 33 implicit none 34 35 private
ABINIT/xchcth [ Functions ]
NAME
xchcth
FUNCTION
Treat XC GGA functional of HCTH type. See Hamprecht, Cohen, Tozer and Handy, J. Chem. Phys. 109, 6264 (1998) [[cite:Hamprecht1998]] for HCTH-93. Boese, Doltsinis, Handy and Sprik, J. Chem. Phys. 112, 1670 (2000) [[cite:Boese2000]] for HCTH-120 and HCTH-147. Boese and Handy , J. Chem. Phys. 114, 5497 (2001) [[cite:Boese2001]] for HCTH-407. For a series of values of the density and the square of the gradient of the density, return the associated Exc energy, potential, and, in case of response-function, functions needed to build the XC kernel.
INPUTS
ixc=number of the XC functional : 16 for HCTH-93, 17 for HCTH-120, 26 for HCTH-147 and 27 for HCTH-407.
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)
order=maximal derivative of Exc to be computed
(1 => energy and potential, or 2 => also XC kernel )
Warning : order=2 not yet available
rho_updn(npts,nspden)=spin-up and spin-down density (Hartree/bohr**3)
OUTPUT
dvxcdgr(npts,3)=partial derivative of the exchange-correlation energy (exci*$\rho$) with respect to the spin-up (dvxcdgr(:,1)), spin-down (dvxcdgr(:,2)) square of gradients of the density exci(npts)=exchange-correlation energy density (hartree) vxci(npts,nspden)=partial derivative of the exchange-correlation energy (exci*$\rho$) with respect to the spin-down (vxci(:,1)) and spin-up (vxci(:,2) densities Normalization: Exc=$\int (exc(r)*\rho (r) d^3 r)$ for $\rho$(r)=electron density.
TODO
Response function not coded yet, but part of it are already present
SOURCE
86 subroutine xchcth(dvxcdgr,exci,grho2_updn,ixc,npts,nspden,order,rho_updn,vxci) 87 88 !Arguments ------------------------------------ 89 !scalars 90 integer,intent(in) :: ixc,npts,nspden,order 91 !arrays 92 real(dp),intent(in) :: grho2_updn(npts,2*nspden-1),rho_updn(npts,nspden) 93 real(dp),intent(out) :: dvxcdgr(npts,2),exci(npts),vxci(npts,nspden) 94 95 !Local variables------------------------------- 96 !scalars 97 integer,save :: initialized=0 98 integer :: ipts,ispden 99 real(dp),parameter :: alpha_zeta2=one-1e-6_dp,alpha_zeta=one-1e-6_dp 100 real(dp),parameter :: fsec_inv=one/1.709921_dp,gammacab=0.006_dp 101 real(dp),parameter :: gammacsig=0.2_dp,gammax=0.004_dp 102 real(dp),parameter :: rsfac=0.6203504908994000_dp 103 real(dp),save :: factf_zeta,factfp_zeta,sixpi2_1_3,sixpi2m1_3,sq_rsfac 104 real(dp),save :: sq_rsfac_inv,threefourth_divpi,twom1_3 105 real(dp) :: ccab0,ccab1,ccab2,ccab3,ccab4,ccsig0,ccsig1,ccsig2,ccsig3,ccsig4 106 real(dp) :: coeffss,cxsig0,cxsig1,cxsig2,cxsig3,cxsig4,d2ecrs0_drs2 107 real(dp) :: d2ecrs1_drs2,d2ecrs_drs2,d2ecrs_drsdzeta,d2ecrs_dzeta2 108 real(dp) :: d2fzeta4_dzeta2,d2gcrs_drs2,d2macrs_drs2,decrs0_drs,decrs1_drho 109 real(dp) :: decrs1_drs,decrs_drs,decrs_dzeta,dfzeta4_dzeta,dgcabdss 110 real(dp) :: dgcrs_drs,dgcsigdss,dgxsigdss,divcab,divcsig,divx,dmacrs_drs 111 real(dp) :: drhoecab_drhodn,drhoecab_drhoup,drhoecrs1_drhodn,drhoecrs1_drhoup 112 real(dp) :: drsdrho,dssdg,dssdndg,dssdndrho,dssdrho,dssupdg,dssupdrho,ducabdss 113 real(dp) :: ducsigdss,duxsigdss,ec0_a1,ec0_aa,ec0_b1,ec0_b2,ec0_b3,ec0_b4 114 real(dp) :: ec0_den,ec0_log,ec0_q0,ec0_q1,ec0_q1p,ec0_q1pp,ec1_a1,ec1_aa 115 real(dp) :: ec1_b1,ec1_b2,ec1_b3,ec1_b4,ec1_den,ec1_log,ec1_q0,ec1_q1,ec1_q1p 116 real(dp) :: ec1_q1pp,ecrs,ecrs0,ecrs1,ex_lsd,exc,f_zeta,factfpp_zeta 117 real(dp) :: fp_zeta,fpp_zeta,gcab,gcrs,gcsig,gxsig,mac_a1,mac_aa,mac_b1 118 real(dp) :: mac_b2,mac_b3,mac_b4,mac_den,mac_log,mac_q0,mac_q1,mac_q1p 119 real(dp) :: mac_q1pp,macrs,rho,rho_inv 120 real(dp) :: rhoecab,rhoecrs1_dn,rhoecrs1_up,rhomo6,rhomot,rhoo6 121 real(dp) :: rhotmo6,rhotmot,rhoto6,rhotot,rhotot_inv,rs,rsm1_2,sqr_rs,ss,ss_dn 122 real(dp) :: ss_up,ssavg,ucab,ucsig,uxsig,vxcadd,zeta,zeta4,zetm_1_3 123 real(dp) :: zetp_1_3 124 character(len=500) :: message 125 !arrays 126 real(dp),allocatable :: rho_updnm1_3(:,:),rhoarr(:),rhom1_3(:),zetm(:) 127 real(dp),allocatable :: zetmm1_3(:),zetp(:),zetpm1_3(:) 128 !no_abirules 129 !real(dp) :: delta,factor,grr,rho_dn,rho_dnm,rho_dnp,rho_up,rho_upm,rho_upp,zeta_mean 130 131 ! ************************************************************************* 132 133 !DEBUG 134 !write(std_out,*)' xchcth : enter' 135 !write(std_out,*)' nspden=',nspden 136 !ENDDEBUG 137 138 if (order/=1) then 139 write(message, '(a,i0)' )' Order must be 1 ; argument was ',order 140 ABI_BUG(message) 141 end if 142 143 if(initialized==0)then 144 twom1_3=two**(-third) 145 sixpi2_1_3=(six*pi**2)**third 146 sixpi2m1_3=one/sixpi2_1_3 147 threefourth_divpi=three_quarters*piinv 148 factf_zeta= one / ( two**(four/three)-two ) 149 factfp_zeta= four_thirds * factf_zeta * alpha_zeta2 150 sq_rsfac=sqrt(rsfac) 151 sq_rsfac_inv=one/sq_rsfac 152 initialized=1 153 end if 154 155 if(ixc==16)then 156 ! HCTH-93 parameters from table II of JCP 109, 6264 (1998). Note that there is 157 ! an error in tables III of JCP 112, 1670 (2000) and JCP 114, 5497 (2001) for the coefficient ccab1 158 cxsig0=1.09320_dp ; cxsig1=-0.744056_dp ; cxsig2=5.59920_dp ; cxsig3=-6.78549_dp ; cxsig4=4.49357_dp 159 ccsig0=0.222601_dp; ccsig1=-0.0338622_dp; ccsig2=-0.0125170_dp; ccsig3=-0.802496_dp; ccsig4=1.55396_dp 160 ccab0=0.729974_dp ; ccab1=3.35287_dp ; ccab2=-11.543_dp ; ccab3=8.08564_dp ; ccab4=-4.47857_dp 161 else if(ixc==17)then 162 ! HCTH-120 parameters from table III of JCP 112, 1670 (2000) 163 ! Note the correction of the sign of cxsig1 and ccsig1, as there is a misprint in the HCTH paper ! 164 ! see the exchange of mail with MMarques 23 Dec 2008, and 6 Jan 2009. 165 ! Now, the functional agrees with the libXC results 166 cxsig0=1.09163_dp; cxsig1=-0.7472_dp; cxsig2=5.0783_dp; cxsig3=-4.1075_dp; cxsig4=1.1717_dp 167 ccsig0=0.48951_dp; ccsig1=-0.2607_dp; ccsig2=0.4329_dp; ccsig3=-1.9925_dp; ccsig4=2.4853_dp 168 ccab0=0.51473_dp ; ccab1=6.9298_dp ; ccab2=-24.707_dp; ccab3=23.110_dp ; ccab4=-11.323_dp 169 ! Exactly same values than in lib_xc 170 cxsig0=1.09163_dp; cxsig1=-0.747215_dp; cxsig2=5.07833_dp; cxsig3=-4.10746_dp; cxsig4=1.17173_dp 171 ccsig0=0.489508_dp; ccsig1=-0.260699_dp; ccsig2=0.432917_dp; ccsig3=-1.99247_dp; ccsig4=2.48531_dp 172 ccab0=0.51473_dp ; ccab1=6.92982_dp ; ccab2=-24.7073_dp; ccab3=23.1098_dp ; ccab4=-11.3234_dp 173 else if(ixc==26)then 174 ! HCTH-147 parameters from table III of JCP 112, 1670 (2000) 175 cxsig0=1.09025_dp; cxsig1=-0.7992_dp; cxsig2=5.5721_dp ; cxsig3=-5.8676_dp; cxsig4=3.0454_dp 176 ccsig0=0.56258_dp; ccsig1=-0.0171_dp; ccsig2=-1.3064_dp; ccsig3= 1.0575_dp; ccsig4=0.8854_dp 177 ccab0=0.54235_dp ; ccab1=7.0146_dp ; ccab2=-28.382_dp ; ccab3=35.033_dp ; ccab4=-20.428_dp 178 ! Exactly same values than in lib_xc 179 cxsig0=1.09025_dp; cxsig1=-0.799194_dp; cxsig2=5.57212_dp ; cxsig3=-5.86760_dp; cxsig4=3.04544_dp 180 ccsig0=0.562576_dp; ccsig1= 0.0171436_dp; ccsig2=-1.30636_dp; ccsig3= 1.05747_dp; ccsig4=0.885429_dp 181 ccab0=0.542352_dp ; ccab1=7.01464_dp ; ccab2=-28.3822_dp ; ccab3=35.0329_dp ; ccab4=-20.4284_dp 182 else if(ixc==27)then 183 ! HCTH-407 parameters from table IV of JCP 114, 5497 (2001) 184 cxsig0=1.08184_dp; cxsig1=-0.5183_dp; cxsig2=3.4256_dp; cxsig3=-2.6290_dp; cxsig4=2.2886_dp 185 ccsig0=1.18777_dp; ccsig1=-2.4029_dp; ccsig2=5.6174_dp; ccsig3=-9.1792_dp; ccsig4=6.2480_dp 186 ccab0=0.58908_dp ; ccab1=4.4237_dp ; ccab2=-19.222_dp; ccab3=42.572_dp ; ccab4=-42.005_dp 187 ! Exactly same values than in lib_xc 188 cxsig0=1.08184_dp; cxsig1=-0.518339_dp; cxsig2=3.42562_dp; cxsig3=-2.62901_dp; cxsig4=2.28855_dp 189 ccsig0=1.18777_dp; ccsig1=-2.40292_dp; ccsig2=5.61741_dp; ccsig3=-9.17923_dp; ccsig4=6.24798_dp 190 ccab0=0.589076_dp ; ccab1=4.42374_dp ; ccab2=-19.2218_dp; ccab3=42.5721_dp ; ccab4=-42.0052_dp 191 else 192 write(message, '(a,i0)' )' xchcth : ixc must be 16, 17, 26, or 27 ; argument was ',ixc 193 ABI_BUG(message) 194 end if 195 196 !Parameters for the Perdew-Wang 92 LSD as well as LSD-RPA, 197 !see Table I of Phys.Rev.B 45,13244 (1992) [[cite:Perdew1992a]] 198 ec0_aa=0.031091_dp ; ec1_aa=0.015545_dp ; mac_aa=0.016887_dp 199 ec0_a1=0.21370_dp ; ec1_a1=0.20548_dp ; mac_a1=0.11125_dp 200 ec0_b1=7.5957_dp ; ec1_b1=14.1189_dp ; mac_b1=10.357_dp 201 ec0_b2=3.5876_dp ; ec1_b2=6.1977_dp ; mac_b2=3.6231_dp 202 ec0_b3=1.6382_dp ; ec1_b3=3.3662_dp ; mac_b3=0.88026_dp 203 ec0_b4=0.49294_dp ; ec1_b4=0.62517_dp ; mac_b4=0.49671_dp 204 205 !DEBUG 206 !Finite-difference debugging, do not take away 207 !Note : here work with collinear gradients. Might be generalized ... 208 !debug=2 ! Choose 1 (rho grads) or 2 (grho grads) 209 !factor=one 210 !zeta_mean=0.98_dp 211 !zeta_mean=zero 212 !delta=0.000025*factor 213 !delta=0.0000125*factor 214 !if(debug/=0)then 215 !do ipts=1,npts,5 216 !rho=ipts*0.zero*factor 217 !rho_up=rho*(one+zeta_mean)*half 218 !rho_dn=rho*(one-zeta_mean)*half 219 !rho_upp=rho_up+delta 220 !rho_upm=rho_up-delta 221 !rho_dnp=rho_dn+delta 222 !rho_dnm=rho_dn-delta 223 !! Here, vary rho 224 !if(debug==1)then 225 !rho_updn(ipts ,1)=rho_up ; rho_updn(ipts ,2)=rho_dn 226 !rho_updn(ipts+1,1)=rho_upp; rho_updn(ipts+1,2)=rho_dn 227 !rho_updn(ipts+2,1)=rho_upm; rho_updn(ipts+2,2)=rho_dn 228 !rho_updn(ipts+3,1)=rho_up ; rho_updn(ipts+3,2)=rho_dnp 229 !rho_updn(ipts+4,1)=rho_up ; rho_updn(ipts+4,2)=rho_dnm 230 !grho2_updn(ipts:ipts+4,1)=(0.2_dp*factor)**2 ! grad2 of spin up density 231 !grho2_updn(ipts:ipts+4,2)=(0.2_dp*factor)**2 ! grad2 of spin down density 232 !grho2_updn(ipts:ipts+4,3)=(0.3_dp*factor)**2 ! grad2 of total density 233 !else 234 !! Here, vary grho (interchange rho and grho) 235 !grho2_updn(ipts ,1)=rho_up**2 ; grho2_updn(ipts ,2)=rho_dn**2 236 !grho2_updn(ipts+1,1)=rho_upp**2; grho2_updn(ipts+1,2)=rho_dn**2 237 !grho2_updn(ipts+2,1)=rho_upm**2; grho2_updn(ipts+2,2)=rho_dn**2 238 !grho2_updn(ipts+3,1)=rho_up**2 ; grho2_updn(ipts+3,2)=rho_dnp**2 239 !grho2_updn(ipts+4,1)=rho_up**2 ; grho2_updn(ipts+4,2)=rho_dnm**2 240 !grho2_updn(ipts ,3)=(ipts*0.zero*factor)**2 241 !grho2_updn(ipts+1,3)=(ipts*0.zero*factor+delta)**2 242 !grho2_updn(ipts+2,3)=(ipts*0.zero*factor-delta)**2 243 !grho2_updn(ipts+3,3)=(ipts*0.zero*factor+delta)**2 ! identical to ipts+1 244 !grho2_updn(ipts+4,3)=(ipts*0.zero*factor-delta)**2 ! identical to ipts+2 245 !rho_updn(ipts:ipts+4,1)=0.2_dp*factor*(one+zeta_mean)*half ! spin up density 246 !rho_updn(ipts:ipts+4,2)=0.2_dp*factor*(one-zeta_mean)*half ! spin down density 247 !end if 248 !end do 249 !end if 250 !Usual option : 251 !nspden=2 ; order=2 252 !GGA 253 !nspden=2 ; order=1 254 !Might take also, although finite difference later is meaningless 255 !nspden=1 ; order=-2 256 !ENDDEBUG 257 258 if(order**2 >1)then 259 factfpp_zeta= third * factfp_zeta * alpha_zeta2 260 end if 261 262 ABI_MALLOC(rhoarr,(npts)) 263 ABI_MALLOC(rhom1_3,(npts)) 264 ABI_MALLOC(rho_updnm1_3,(npts,2)) 265 do ispden=1,nspden 266 call invcb(rho_updn(:,ispden),rho_updnm1_3(:,ispden),npts) 267 end do 268 if(nspden==1)then 269 rhoarr(:)=two*rho_updn(:,1) 270 rhom1_3(:)=twom1_3*rho_updnm1_3(:,1) 271 rho_updnm1_3(:,2)=rho_updnm1_3(:,1) 272 else 273 rhoarr(:)=rho_updn(:,1)+rho_updn(:,2) 274 call invcb(rhoarr,rhom1_3,npts) 275 ABI_MALLOC(zetm,(npts)) 276 ABI_MALLOC(zetmm1_3,(npts)) 277 ABI_MALLOC(zetp,(npts)) 278 ABI_MALLOC(zetpm1_3,(npts)) 279 do ipts=1,npts 280 rhotmot=rhom1_3(ipts) 281 rhotot_inv=rhotmot*rhotmot*rhotmot 282 zeta=(rho_updn(ipts,1)-rho_updn(ipts,2))*rhotot_inv 283 zetp(ipts)=one+zeta*alpha_zeta 284 zetm(ipts)=one-zeta*alpha_zeta 285 end do 286 call invcb(zetp,zetpm1_3,npts) 287 call invcb(zetm,zetmm1_3,npts) 288 end if 289 290 if (nspden==1) then 291 292 if(order==-2) then 293 294 do ipts=1,npts 295 296 ! ----------------------------------------------------------------------- 297 ! First take care of the spin-split part of the functional 298 exc=zero 299 ispden=1 300 rho =rho_updn(ipts,ispden) 301 rhomot=rho_updnm1_3(ipts,ispden) 302 rho_inv=rhomot*rhomot*rhomot 303 304 ! Exchange part 305 ex_lsd= - threefourth_divpi * sixpi2_1_3*rhomot*rhomot*rho 306 ! Note that this definition differs from the PBE one 307 coeffss=rho_inv*rho_inv*rhomot*rhomot 308 ss=grho2_updn(ipts,ispden)*coeffss 309 dssdrho=-eight*third*ss*rho_inv 310 dssdg=two*coeffss 311 312 divx=one/(one+gammax*ss) 313 uxsig=gammax*ss*divx 314 duxsigdss=gammax*divx*(one-ss*gammax*divx) 315 316 gxsig=cxsig0+uxsig*(cxsig1+uxsig*(cxsig2+uxsig*(cxsig3+uxsig*cxsig4))) 317 dgxsigdss=(cxsig1+uxsig*(two*cxsig2+uxsig*(three*cxsig3+uxsig*four*cxsig4)))& 318 & *duxsigdss 319 320 exc=exc+ex_lsd*rho*gxsig 321 vxci(ipts,ispden)=ex_lsd*(four_thirds*gxsig+rho*dgxsigdss*dssdrho) 322 dvxcdgr(ipts,ispden)=ex_lsd*rho*dgxsigdss*dssdg 323 324 ! Spin parallel correlation part 325 ! Note that this definition is for fully spin-polarized quantities 326 rs=rsfac*rhomot 327 rhomo6=sqrt(rhomot) 328 sqr_rs=sq_rsfac*rhomo6 329 rhoo6=rho*rhomot*rhomot*rhomo6 330 rsm1_2=sq_rsfac_inv*rhoo6 331 drsdrho=-third*rs*rho_inv 332 333 ec1_q0=-two*ec1_aa*(one+ec1_a1*rs) 334 ec1_q1=two*ec1_aa*(ec1_b1*sqr_rs+ec1_b2*rs+ec1_b3*rs*sqr_rs+ec1_b4*rs*rs) 335 ec1_q1p=ec1_aa*(ec1_b1*rsm1_2+two*ec1_b2+three*ec1_b3*sqr_rs+four*ec1_b4*rs) 336 ec1_den=one/(ec1_q1*ec1_q1+ec1_q1) 337 ! ec1_log=log( one + one / ec1_q1 ) 338 ec1_log=-log( ec1_q1*ec1_q1*ec1_den ) 339 ecrs1=ec1_q0*ec1_log 340 decrs1_drs= -two*ec1_aa*ec1_a1*ec1_log - ec1_q0*ec1_q1p *ec1_den 341 decrs1_drho=ecrs1+decrs1_drs*drsdrho*rho 342 343 ! Store the LSDA corr energy and density derivative 344 rhoecrs1_up=rho*ecrs1 345 drhoecrs1_drhoup=decrs1_drho 346 ss_up=ss 347 dssupdrho=dssdrho 348 dssupdg=dssdg 349 350 divcsig=one/(one+gammacsig*ss) 351 ucsig=gammacsig*ss*divcsig 352 ducsigdss=gammacsig*divcsig*(one-ss*gammacsig*divcsig) 353 354 gcsig=ccsig0+ucsig*(ccsig1+ucsig*(ccsig2+ucsig*(ccsig3+ucsig*ccsig4))) 355 dgcsigdss=(ccsig1+ucsig*(two*ccsig2+ucsig*(three*ccsig3+ucsig*four*ccsig4)))& 356 & *ducsigdss 357 358 ! DEBUG 359 ! gcsig=zero 360 ! dgcsigdss=zero 361 ! ENDDEBUG 362 exc=exc+ecrs1*rho*gcsig 363 vxci(ipts,ispden)=vxci(ipts,ispden)+decrs1_drho*gcsig+& 364 & ecrs1*rho*dgcsigdss*dssdrho 365 dvxcdgr(ipts,ispden)=dvxcdgr(ipts,ispden)+ecrs1*rho*dgcsigdss*dssdg 366 367 rhoecrs1_dn=rhoecrs1_up 368 drhoecrs1_drhodn=drhoecrs1_drhoup 369 ss_dn=ss_up 370 dssdndrho=dssupdrho 371 dssdndg=dssupdg 372 exc=exc*2 373 374 ! ----------------------------------------------------------------------------- 375 ! Then takes care of the LSD correlation part of the functional 376 377 rhotot=rhoarr(ipts) 378 rhotmot=rhom1_3(ipts) 379 rhotot_inv=rhotmot*rhotmot*rhotmot 380 rhotmo6=sqrt(rhotmot) 381 rhoto6=rhotot*rhotmot*rhotmot*rhotmo6 382 383 ! From now, the coding of the PW92 functional start. It is identical in xcpbe.f 384 rs=rsfac*rhotmot 385 sqr_rs=sq_rsfac*rhotmo6 386 rsm1_2=sq_rsfac_inv*rhoto6 387 388 ! Formulas A6-A8 of PW92LSD 389 ec0_q0=-2.0d0*ec0_aa*(one+ec0_a1*rs) 390 ec0_q1=2.0d0*ec0_aa*(ec0_b1*sqr_rs+ec0_b2*rs+ec0_b3*rs*sqr_rs+ec0_b4*rs*rs) 391 ec0_q1p=ec0_aa*(ec0_b1*rsm1_2+2.d0*ec0_b2+3.d0*ec0_b3*sqr_rs+4.d0*ec0_b4*rs) 392 ec0_den=one/(ec0_q1*ec0_q1+ec0_q1) 393 ! ec0_log=log( one + one / ec0_q1 ) 394 ec0_log=-log( ec0_q1*ec0_q1*ec0_den ) 395 ecrs0=ec0_q0*ec0_log 396 decrs0_drs= -2.0d0*ec0_aa*ec0_a1*ec0_log - ec0_q0*ec0_q1p *ec0_den 397 398 ec0_q1pp=half*ec0_aa*(-ec0_b1*rsm1_2**3+3.d0*ec0_b3*rsm1_2+8.d0*ec0_b4) 399 d2ecrs0_drs2= 4.0d0*ec0_aa*ec0_a1*ec0_q1p*ec0_den & 400 & -ec0_q0*ec0_q1pp*ec0_den & 401 & +ec0_q0*ec0_q1p**2*ec0_den**2*(2.d0*ec0_q1+one) 402 403 404 405 mac_q0=-2.0d0*mac_aa*(one+mac_a1*rs) 406 mac_q1=2.0d0*mac_aa*(mac_b1*sqr_rs+mac_b2*rs+mac_b3*rs*sqr_rs+mac_b4*rs*rs) 407 mac_q1p=mac_aa*(mac_b1*rsm1_2+2.d0*mac_b2+3.d0*mac_b3*sqr_rs+4.d0*mac_b4*rs) 408 mac_den=one/(mac_q1*mac_q1+mac_q1) 409 mac_log=-log( mac_q1*mac_q1*mac_den ) 410 macrs=mac_q0*mac_log 411 dmacrs_drs= -2.0d0*mac_aa*mac_a1*mac_log - mac_q0*mac_q1p*mac_den 412 413 414 ecrs=ecrs0 415 decrs_drs=decrs0_drs 416 decrs_dzeta=zero 417 418 d2ecrs_drs2=d2ecrs0_drs2 419 420 d2ecrs_dzeta2=alpha_zeta**2*(-macrs) 421 422 zeta=zero 423 424 ! At this point, the coding of the PW92 functional finishes. 425 426 ! The correlation between different spin from HCTH is now computed 427 ! First, the part without gradient correction factor 428 rhoecab=ecrs*rhotot-rhoecrs1_dn-rhoecrs1_up 429 vxcadd=ecrs-rs*third*decrs_drs-zeta*decrs_dzeta 430 drhoecab_drhoup=vxcadd+decrs_dzeta-drhoecrs1_drhoup 431 drhoecab_drhodn=vxcadd-decrs_dzeta-drhoecrs1_drhodn 432 433 ! Now, the gradient correction factor 434 ssavg=half*(ss_up+ss_dn) 435 divcab=one/(one+gammacab*ssavg) 436 ucab=gammacab*ssavg*divcab 437 ducabdss=gammacab*divcab*(one-ssavg*gammacab*divcab) 438 439 gcab=ccab0+ucab*(ccab1+ucab*(ccab2+ucab*(ccab3+ucab*ccab4))) 440 dgcabdss=(ccab1+ucab*(two*ccab2+ucab*(three*ccab3+ucab*four*ccab4)))& 441 & *ducabdss 442 443 exc=exc+rhoecab*gcab 444 445 vxci(ipts,1)=vxci(ipts,1)+drhoecab_drhoup*gcab+rhoecab*dgcabdss*half*dssupdrho 446 dvxcdgr(ipts,1)=dvxcdgr(ipts,1)+rhoecab*dgcabdss*half*dssupdg 447 ! If non spin-polarized, treat spin down contribution now, similar to spin up 448 ! vxci(ipts,2)=vxci(ipts,1) 449 dvxcdgr(ipts,2)=dvxcdgr(ipts,1) 450 451 ! Final division by the total density, to give the energy density 452 exci(ipts)=exc*rhotot_inv 453 454 end do ! ipts=1,npts 455 456 else if(order**2>1) then 457 458 do ipts=1,npts 459 460 ! ----------------------------------------------------------------------- 461 ! First take care of the spin-split part of the functional 462 exc=zero 463 ispden=1 464 rho =rho_updn(ipts,ispden) 465 rhomot=rho_updnm1_3(ipts,ispden) 466 rho_inv=rhomot*rhomot*rhomot 467 468 ! Exchange part 469 ex_lsd= - threefourth_divpi * sixpi2_1_3*rhomot*rhomot*rho 470 ! Note that this definition differs from the PBE one 471 coeffss=rho_inv*rho_inv*rhomot*rhomot 472 ss=grho2_updn(ipts,ispden)*coeffss 473 dssdrho=-eight*third*ss*rho_inv 474 dssdg=two*coeffss 475 476 divx=one/(one+gammax*ss) 477 uxsig=gammax*ss*divx 478 duxsigdss=gammax*divx*(one-ss*gammax*divx) 479 480 gxsig=cxsig0+uxsig*(cxsig1+uxsig*(cxsig2+uxsig*(cxsig3+uxsig*cxsig4))) 481 dgxsigdss=(cxsig1+uxsig*(two*cxsig2+uxsig*(three*cxsig3+uxsig*four*cxsig4)))& 482 & *duxsigdss 483 484 exc=exc+ex_lsd*rho*gxsig 485 vxci(ipts,ispden)=ex_lsd*(four_thirds*gxsig+rho*dgxsigdss*dssdrho) 486 dvxcdgr(ipts,ispden)=ex_lsd*rho*dgxsigdss*dssdg 487 488 ! Spin parallel correlation part 489 ! Note that this definition is for fully spin-polarized quantities 490 rs=rsfac*rhomot 491 rhomo6=sqrt(rhomot) 492 sqr_rs=sq_rsfac*rhomo6 493 rhoo6=rho*rhomot*rhomot*rhomo6 494 rsm1_2=sq_rsfac_inv*rhoo6 495 drsdrho=-third*rs*rho_inv 496 497 ec1_q0=-two*ec1_aa*(one+ec1_a1*rs) 498 ec1_q1=two*ec1_aa*(ec1_b1*sqr_rs+ec1_b2*rs+ec1_b3*rs*sqr_rs+ec1_b4*rs*rs) 499 ec1_q1p=ec1_aa*(ec1_b1*rsm1_2+two*ec1_b2+three*ec1_b3*sqr_rs+four*ec1_b4*rs) 500 ec1_den=one/(ec1_q1*ec1_q1+ec1_q1) 501 ! ec1_log=log( one + one / ec1_q1 ) 502 ec1_log=-log( ec1_q1*ec1_q1*ec1_den ) 503 ecrs1=ec1_q0*ec1_log 504 decrs1_drs= -two*ec1_aa*ec1_a1*ec1_log - ec1_q0*ec1_q1p *ec1_den 505 decrs1_drho=ecrs1+decrs1_drs*drsdrho*rho 506 507 ! Store the LSDA corr energy and density derivative 508 rhoecrs1_up=rho*ecrs1 509 drhoecrs1_drhoup=decrs1_drho 510 ss_up=ss 511 dssupdrho=dssdrho 512 dssupdg=dssdg 513 514 divcsig=one/(one+gammacsig*ss) 515 ucsig=gammacsig*ss*divcsig 516 ducsigdss=gammacsig*divcsig*(one-ss*gammacsig*divcsig) 517 518 gcsig=ccsig0+ucsig*(ccsig1+ucsig*(ccsig2+ucsig*(ccsig3+ucsig*ccsig4))) 519 dgcsigdss=(ccsig1+ucsig*(two*ccsig2+ucsig*(three*ccsig3+ucsig*four*ccsig4)))& 520 & *ducsigdss 521 522 ! DEBUG 523 ! gcsig=zero 524 ! dgcsigdss=zero 525 ! ENDDEBUG 526 exc=exc+ecrs1*rho*gcsig 527 vxci(ipts,ispden)=vxci(ipts,ispden)+decrs1_drho*gcsig+& 528 & ecrs1*rho*dgcsigdss*dssdrho 529 dvxcdgr(ipts,ispden)=dvxcdgr(ipts,ispden)+ecrs1*rho*dgcsigdss*dssdg 530 531 rhoecrs1_dn=rhoecrs1_up 532 drhoecrs1_drhodn=drhoecrs1_drhoup 533 ss_dn=ss_up 534 dssdndrho=dssupdrho 535 dssdndg=dssupdg 536 exc=exc*2 537 538 ! ----------------------------------------------------------------------------- 539 ! Then takes care of the LSD correlation part of the functional 540 541 rhotot=rhoarr(ipts) 542 rhotmot=rhom1_3(ipts) 543 rhotot_inv=rhotmot*rhotmot*rhotmot 544 rhotmo6=sqrt(rhotmot) 545 rhoto6=rhotot*rhotmot*rhotmot*rhotmo6 546 547 ! From now, the coding of the PW92 functional start. It is identical in xcpbe.f 548 rs=rsfac*rhotmot 549 sqr_rs=sq_rsfac*rhotmo6 550 rsm1_2=sq_rsfac_inv*rhoto6 551 552 ! Formulas A6-A8 of PW92LSD 553 ec0_q0=-2.0d0*ec0_aa*(one+ec0_a1*rs) 554 ec0_q1=2.0d0*ec0_aa*(ec0_b1*sqr_rs+ec0_b2*rs+ec0_b3*rs*sqr_rs+ec0_b4*rs*rs) 555 ec0_q1p=ec0_aa*(ec0_b1*rsm1_2+2.d0*ec0_b2+3.d0*ec0_b3*sqr_rs+4.d0*ec0_b4*rs) 556 ec0_den=one/(ec0_q1*ec0_q1+ec0_q1) 557 ! ec0_log=log( one + one / ec0_q1 ) 558 ec0_log=-log( ec0_q1*ec0_q1*ec0_den ) 559 ecrs0=ec0_q0*ec0_log 560 decrs0_drs= -2.0d0*ec0_aa*ec0_a1*ec0_log - ec0_q0*ec0_q1p *ec0_den 561 ec0_q1pp=half*ec0_aa*(-ec0_b1*rsm1_2**3+3.d0*ec0_b3*rsm1_2+8.d0*ec0_b4) 562 d2ecrs0_drs2= 4.0d0*ec0_aa*ec0_a1*ec0_q1p*ec0_den & 563 & -ec0_q0*ec0_q1pp*ec0_den & 564 & +ec0_q0*ec0_q1p**2*ec0_den**2*(2.d0*ec0_q1+one) 565 566 567 ecrs=ecrs0 568 decrs_drs=decrs0_drs 569 decrs_dzeta=zero 570 d2ecrs_drs2=d2ecrs0_drs2 571 zeta=zero 572 573 ! At this point, the coding of the PW92 functional finishes. 574 575 ! The correlation between different spin from HCTH is now computed 576 ! First, the part without gradient correction factor 577 rhoecab=ecrs*rhotot-rhoecrs1_dn-rhoecrs1_up 578 vxcadd=ecrs-rs*third*decrs_drs-zeta*decrs_dzeta 579 drhoecab_drhoup=vxcadd+decrs_dzeta-drhoecrs1_drhoup 580 drhoecab_drhodn=vxcadd-decrs_dzeta-drhoecrs1_drhodn 581 582 ! Now, the gradient correction factor 583 ssavg=half*(ss_up+ss_dn) 584 divcab=one/(one+gammacab*ssavg) 585 ucab=gammacab*ssavg*divcab 586 ducabdss=gammacab*divcab*(one-ssavg*gammacab*divcab) 587 588 gcab=ccab0+ucab*(ccab1+ucab*(ccab2+ucab*(ccab3+ucab*ccab4))) 589 dgcabdss=(ccab1+ucab*(two*ccab2+ucab*(three*ccab3+ucab*four*ccab4)))& 590 & *ducabdss 591 592 exc=exc+rhoecab*gcab 593 594 vxci(ipts,1)=vxci(ipts,1)+drhoecab_drhoup*gcab+rhoecab*dgcabdss*half*dssupdrho 595 dvxcdgr(ipts,1)=dvxcdgr(ipts,1)+rhoecab*dgcabdss*half*dssupdg 596 ! If non spin-polarized, treat spin down contribution now, similar to spin up 597 ! vxci(ipts,2)=vxci(ipts,1) 598 dvxcdgr(ipts,2)=dvxcdgr(ipts,1) 599 600 ! Final division by the total density, to give the energy density 601 exci(ipts)=exc*rhotot_inv 602 603 end do ! ipts=1,npts 604 else 605 606 do ipts=1,npts 607 608 ! ----------------------------------------------------------------------- 609 ! First take care of the spin-split part of the functional 610 exc=zero 611 ispden=1 612 rho =rho_updn(ipts,ispden) 613 rhomot=rho_updnm1_3(ipts,ispden) 614 rho_inv=rhomot*rhomot*rhomot 615 616 ! Exchange part 617 ex_lsd= - threefourth_divpi * sixpi2_1_3*rhomot*rhomot*rho 618 ! Note that this definition differs from the PBE one 619 coeffss=rho_inv*rho_inv*rhomot*rhomot 620 ss=grho2_updn(ipts,ispden)*coeffss 621 dssdrho=-eight*third*ss*rho_inv 622 dssdg=two*coeffss 623 624 divx=one/(one+gammax*ss) 625 uxsig=gammax*ss*divx 626 duxsigdss=gammax*divx*(one-ss*gammax*divx) 627 628 gxsig=cxsig0+uxsig*(cxsig1+uxsig*(cxsig2+uxsig*(cxsig3+uxsig*cxsig4))) 629 dgxsigdss=(cxsig1+uxsig*(two*cxsig2+uxsig*(three*cxsig3+uxsig*four*cxsig4)))& 630 & *duxsigdss 631 632 exc=exc+ex_lsd*rho*gxsig 633 vxci(ipts,ispden)=ex_lsd*(four_thirds*gxsig+rho*dgxsigdss*dssdrho) 634 dvxcdgr(ipts,ispden)=ex_lsd*rho*dgxsigdss*dssdg 635 636 ! Spin parallel correlation part 637 ! Note that this definition is for fully spin-polarized quantities 638 rs=rsfac*rhomot 639 rhomo6=sqrt(rhomot) 640 sqr_rs=sq_rsfac*rhomo6 641 rhoo6=rho*rhomot*rhomot*rhomo6 642 rsm1_2=sq_rsfac_inv*rhoo6 643 drsdrho=-third*rs*rho_inv 644 645 ec1_q0=-two*ec1_aa*(one+ec1_a1*rs) 646 ec1_q1=two*ec1_aa*(ec1_b1*sqr_rs+ec1_b2*rs+ec1_b3*rs*sqr_rs+ec1_b4*rs*rs) 647 ec1_q1p=ec1_aa*(ec1_b1*rsm1_2+two*ec1_b2+three*ec1_b3*sqr_rs+four*ec1_b4*rs) 648 ec1_den=one/(ec1_q1*ec1_q1+ec1_q1) 649 ! ec1_log=log( one + one / ec1_q1 ) 650 ec1_log=-log( ec1_q1*ec1_q1*ec1_den ) 651 ecrs1=ec1_q0*ec1_log 652 decrs1_drs= -two*ec1_aa*ec1_a1*ec1_log - ec1_q0*ec1_q1p *ec1_den 653 decrs1_drho=ecrs1+decrs1_drs*drsdrho*rho 654 655 ! Store the LSDA corr energy and density derivative 656 rhoecrs1_up=rho*ecrs1 657 drhoecrs1_drhoup=decrs1_drho 658 ss_up=ss 659 dssupdrho=dssdrho 660 dssupdg=dssdg 661 662 divcsig=one/(one+gammacsig*ss) 663 ucsig=gammacsig*ss*divcsig 664 ducsigdss=gammacsig*divcsig*(one-ss*gammacsig*divcsig) 665 666 gcsig=ccsig0+ucsig*(ccsig1+ucsig*(ccsig2+ucsig*(ccsig3+ucsig*ccsig4))) 667 dgcsigdss=(ccsig1+ucsig*(two*ccsig2+ucsig*(three*ccsig3+ucsig*four*ccsig4)))& 668 & *ducsigdss 669 670 ! DEBUG 671 ! gcsig=zero 672 ! dgcsigdss=zero 673 ! ENDDEBUG 674 exc=exc+ecrs1*rho*gcsig 675 vxci(ipts,ispden)=vxci(ipts,ispden)+decrs1_drho*gcsig+& 676 & ecrs1*rho*dgcsigdss*dssdrho 677 dvxcdgr(ipts,ispden)=dvxcdgr(ipts,ispden)+ecrs1*rho*dgcsigdss*dssdg 678 679 rhoecrs1_dn=rhoecrs1_up 680 drhoecrs1_drhodn=drhoecrs1_drhoup 681 ss_dn=ss_up 682 dssdndrho=dssupdrho 683 dssdndg=dssupdg 684 exc=exc*2 685 686 ! ----------------------------------------------------------------------------- 687 ! Then takes care of the LSD correlation part of the functional 688 689 rhotot=rhoarr(ipts) 690 rhotmot=rhom1_3(ipts) 691 rhotot_inv=rhotmot*rhotmot*rhotmot 692 rhotmo6=sqrt(rhotmot) 693 rhoto6=rhotot*rhotmot*rhotmot*rhotmo6 694 695 ! From now, the coding of the PW92 functional start. It is identical in xcpbe.f 696 rs=rsfac*rhotmot 697 sqr_rs=sq_rsfac*rhotmo6 698 rsm1_2=sq_rsfac_inv*rhoto6 699 700 ! Formulas A6-A8 of PW92LSD 701 ec0_q0=-2.0d0*ec0_aa*(one+ec0_a1*rs) 702 ec0_q1=2.0d0*ec0_aa*(ec0_b1*sqr_rs+ec0_b2*rs+ec0_b3*rs*sqr_rs+ec0_b4*rs*rs) 703 ec0_q1p=ec0_aa*(ec0_b1*rsm1_2+2.d0*ec0_b2+3.d0*ec0_b3*sqr_rs+4.d0*ec0_b4*rs) 704 ec0_den=one/(ec0_q1*ec0_q1+ec0_q1) 705 ! ec0_log=log( one + one / ec0_q1 ) 706 ec0_log=-log( ec0_q1*ec0_q1*ec0_den ) 707 ecrs0=ec0_q0*ec0_log 708 decrs0_drs= -2.0d0*ec0_aa*ec0_a1*ec0_log - ec0_q0*ec0_q1p *ec0_den 709 710 ecrs=ecrs0 711 decrs_drs=decrs0_drs 712 decrs_dzeta=zero 713 zeta=zero 714 715 ! At this point, the coding of the PW92 functional finishes. 716 717 ! The correlation between different spin from HCTH is now computed 718 ! First, the part without gradient correction factor 719 rhoecab=ecrs*rhotot-rhoecrs1_dn-rhoecrs1_up 720 vxcadd=ecrs-rs*third*decrs_drs-zeta*decrs_dzeta 721 drhoecab_drhoup=vxcadd+decrs_dzeta-drhoecrs1_drhoup 722 drhoecab_drhodn=vxcadd-decrs_dzeta-drhoecrs1_drhodn 723 724 ! Now, the gradient correction factor 725 ssavg=half*(ss_up+ss_dn) 726 divcab=one/(one+gammacab*ssavg) 727 ucab=gammacab*ssavg*divcab 728 ducabdss=gammacab*divcab*(one-ssavg*gammacab*divcab) 729 730 gcab=ccab0+ucab*(ccab1+ucab*(ccab2+ucab*(ccab3+ucab*ccab4))) 731 dgcabdss=(ccab1+ucab*(two*ccab2+ucab*(three*ccab3+ucab*four*ccab4)))& 732 & *ducabdss 733 734 exc=exc+rhoecab*gcab 735 736 vxci(ipts,1)=vxci(ipts,1)+drhoecab_drhoup*gcab+rhoecab*dgcabdss*half*dssupdrho 737 dvxcdgr(ipts,1)=dvxcdgr(ipts,1)+rhoecab*dgcabdss*half*dssupdg 738 ! If non spin-polarized, treat spin down contribution now, similar to spin up 739 ! vxci(ipts,2)=vxci(ipts,1) 740 dvxcdgr(ipts,2)=dvxcdgr(ipts,1) 741 742 ! Final division by the total density, to give the energy density 743 exci(ipts)=exc*rhotot_inv 744 745 end do ! ipts=1,npts 746 end if 747 748 else if(nspden==2) then 749 750 if(order**2>1) then 751 752 do ipts=1,npts 753 754 ! ----------------------------------------------------------------------- 755 ! First take care of the spin-split part of the functional 756 exc=zero 757 do ispden=1,nspden 758 rho =rho_updn(ipts,ispden) 759 rhomot=rho_updnm1_3(ipts,ispden) 760 rho_inv=rhomot*rhomot*rhomot 761 762 ! Exchange part 763 ex_lsd= - threefourth_divpi * sixpi2_1_3*rhomot*rhomot*rho 764 ! Note that this definition differs from the PBE one 765 coeffss=rho_inv*rho_inv*rhomot*rhomot 766 ss=grho2_updn(ipts,ispden)*coeffss 767 dssdrho=-eight*third*ss*rho_inv 768 dssdg=two*coeffss 769 770 divx=one/(one+gammax*ss) 771 uxsig=gammax*ss*divx 772 duxsigdss=gammax*divx*(one-ss*gammax*divx) 773 774 gxsig=cxsig0+uxsig*(cxsig1+uxsig*(cxsig2+uxsig*(cxsig3+uxsig*cxsig4))) 775 dgxsigdss=(cxsig1+uxsig*(two*cxsig2+uxsig*(three*cxsig3+uxsig*four*cxsig4)))& 776 & *duxsigdss 777 778 exc=exc+ex_lsd*rho*gxsig 779 vxci(ipts,ispden)=ex_lsd*(four_thirds*gxsig+rho*dgxsigdss*dssdrho) 780 dvxcdgr(ipts,ispden)=ex_lsd*rho*dgxsigdss*dssdg 781 782 ! Spin parallel correlation part 783 ! Note that this definition is for fully spin-polarized quantities 784 rs=rsfac*rhomot 785 rhomo6=sqrt(rhomot) 786 sqr_rs=sq_rsfac*rhomo6 787 rhoo6=rho*rhomot*rhomot*rhomo6 788 rsm1_2=sq_rsfac_inv*rhoo6 789 drsdrho=-third*rs*rho_inv 790 791 ec1_q0=-two*ec1_aa*(one+ec1_a1*rs) 792 ec1_q1=two*ec1_aa*(ec1_b1*sqr_rs+ec1_b2*rs+ec1_b3*rs*sqr_rs+ec1_b4*rs*rs) 793 ec1_q1p=ec1_aa*(ec1_b1*rsm1_2+two*ec1_b2+three*ec1_b3*sqr_rs+four*ec1_b4*rs) 794 ec1_den=one/(ec1_q1*ec1_q1+ec1_q1) 795 ! ec1_log=log( one + one / ec1_q1 ) 796 ec1_log=-log( ec1_q1*ec1_q1*ec1_den ) 797 ecrs1=ec1_q0*ec1_log 798 decrs1_drs= -two*ec1_aa*ec1_a1*ec1_log - ec1_q0*ec1_q1p *ec1_den 799 decrs1_drho=ecrs1+decrs1_drs*drsdrho*rho 800 801 ! Store the LSDA corr energy and density derivative 802 if(ispden==1)then 803 rhoecrs1_up=rho*ecrs1 804 drhoecrs1_drhoup=decrs1_drho 805 ss_up=ss 806 dssupdrho=dssdrho 807 dssupdg=dssdg 808 else 809 rhoecrs1_dn=rho*ecrs1 810 drhoecrs1_drhodn=decrs1_drho 811 ss_dn=ss 812 dssdndrho=dssdrho 813 dssdndg=dssdg 814 end if 815 816 divcsig=one/(one+gammacsig*ss) 817 ucsig=gammacsig*ss*divcsig 818 ducsigdss=gammacsig*divcsig*(one-ss*gammacsig*divcsig) 819 820 gcsig=ccsig0+ucsig*(ccsig1+ucsig*(ccsig2+ucsig*(ccsig3+ucsig*ccsig4))) 821 dgcsigdss=(ccsig1+ucsig*(two*ccsig2+ucsig*(three*ccsig3+ucsig*four*ccsig4)))& 822 & *ducsigdss 823 824 ! DEBUG 825 ! gcsig=zero 826 ! dgcsigdss=zero 827 ! ENDDEBUG 828 exc=exc+ecrs1*rho*gcsig 829 vxci(ipts,ispden)=vxci(ipts,ispden)+decrs1_drho*gcsig+& 830 & ecrs1*rho*dgcsigdss*dssdrho 831 dvxcdgr(ipts,ispden)=dvxcdgr(ipts,ispden)+ecrs1*rho*dgcsigdss*dssdg 832 833 end do 834 835 ! ----------------------------------------------------------------------------- 836 ! Then takes care of the LSD correlation part of the functional 837 838 rhotot=rhoarr(ipts) 839 rhotmot=rhom1_3(ipts) 840 rhotot_inv=rhotmot*rhotmot*rhotmot 841 rhotmo6=sqrt(rhotmot) 842 rhoto6=rhotot*rhotmot*rhotmot*rhotmo6 843 844 ! From now, the coding of the PW92 functional start. It is identical in xcpbe.f 845 rs=rsfac*rhotmot 846 sqr_rs=sq_rsfac*rhotmo6 847 rsm1_2=sq_rsfac_inv*rhoto6 848 849 ! Formulas A6-A8 of PW92LSD 850 ec0_q0=-2.0d0*ec0_aa*(one+ec0_a1*rs) 851 ec0_q1=2.0d0*ec0_aa*(ec0_b1*sqr_rs+ec0_b2*rs+ec0_b3*rs*sqr_rs+ec0_b4*rs*rs) 852 ec0_q1p=ec0_aa*(ec0_b1*rsm1_2+2.d0*ec0_b2+3.d0*ec0_b3*sqr_rs+4.d0*ec0_b4*rs) 853 ec0_den=one/(ec0_q1*ec0_q1+ec0_q1) 854 ! ec0_log=log( one + one / ec0_q1 ) 855 ec0_log=-log( ec0_q1*ec0_q1*ec0_den ) 856 ecrs0=ec0_q0*ec0_log 857 decrs0_drs= -2.0d0*ec0_aa*ec0_a1*ec0_log - ec0_q0*ec0_q1p *ec0_den 858 ec0_q1pp=half*ec0_aa*(-ec0_b1*rsm1_2**3+3.d0*ec0_b3*rsm1_2+8.d0*ec0_b4) 859 d2ecrs0_drs2= 4.0d0*ec0_aa*ec0_a1*ec0_q1p*ec0_den & 860 & -ec0_q0*ec0_q1pp*ec0_den & 861 & +ec0_q0*ec0_q1p**2*ec0_den**2*(2.d0*ec0_q1+one) 862 863 mac_q0=-2.0d0*mac_aa*(one+mac_a1*rs) 864 mac_q1=2.0d0*mac_aa*(mac_b1*sqr_rs+mac_b2*rs+mac_b3*rs*sqr_rs+mac_b4*rs*rs) 865 mac_q1p=mac_aa*(mac_b1*rsm1_2+2.d0*mac_b2+3.d0*mac_b3*sqr_rs+4.d0*mac_b4*rs) 866 mac_den=one/(mac_q1*mac_q1+mac_q1) 867 mac_log=-log( mac_q1*mac_q1*mac_den ) 868 macrs=mac_q0*mac_log 869 dmacrs_drs= -2.0d0*mac_aa*mac_a1*mac_log - mac_q0*mac_q1p*mac_den 870 871 zeta=(rho_updn(ipts,1)-rho_updn(ipts,2))*rhotot_inv 872 ec1_q0=-two*ec1_aa*(one+ec1_a1*rs) 873 ec1_q1=two*ec1_aa*(ec1_b1*sqr_rs+ec1_b2*rs+ec1_b3*rs*sqr_rs+ec1_b4*rs*rs) 874 ec1_q1p=ec1_aa*(ec1_b1*rsm1_2+two*ec1_b2+three*ec1_b3*sqr_rs+four*ec1_b4*rs) 875 ec1_den=one/(ec1_q1*ec1_q1+ec1_q1) 876 ! ec1_log=log( one + one / ec1_q1 ) 877 ec1_log=-log( ec1_q1*ec1_q1*ec1_den ) 878 ecrs1=ec1_q0*ec1_log 879 decrs1_drs= -two*ec1_aa*ec1_a1*ec1_log - ec1_q0*ec1_q1p *ec1_den 880 881 ! alpha_zeta is introduced in order to remove singularities for fully 882 ! polarized systems. 883 zetp_1_3=(one+zeta*alpha_zeta)*zetpm1_3(ipts)**2 884 zetm_1_3=(one-zeta*alpha_zeta)*zetmm1_3(ipts)**2 885 886 f_zeta=( (one+zeta*alpha_zeta2)*zetp_1_3 + & 887 & (one-zeta*alpha_zeta2)*zetm_1_3 - 2.0d0 ) * factf_zeta 888 fp_zeta=( zetp_1_3 - zetm_1_3 ) * factfp_zeta 889 zeta4=zeta**4 890 gcrs=ecrs1-ecrs0+macrs*fsec_inv 891 ! ecrs=ecrs0+f_zeta*(-macrs*(one-zeta4)*fsec_inv+(ecrs1-ecrs0)*zeta4) 892 ecrs=ecrs0+f_zeta*(zeta4*gcrs-macrs*fsec_inv) 893 894 dgcrs_drs=decrs1_drs-decrs0_drs+dmacrs_drs*fsec_inv 895 ! decrs_drs=decrs0_drs+f_zeta*& 896 ! & (-dmacrs_drs*(one-zeta4)*fsec_inv+(decrs1_drs-decrs0_drs)*zeta4) 897 decrs_drs=decrs0_drs+f_zeta*(zeta4*dgcrs_drs-dmacrs_drs*fsec_inv) 898 dfzeta4_dzeta=4.0d0*zeta**3*f_zeta+fp_zeta*zeta4 899 decrs_dzeta=dfzeta4_dzeta*gcrs-fp_zeta*macrs*fsec_inv 900 901 ec1_q1pp=half*ec1_aa*(-ec1_b1*rsm1_2**3+3.d0*ec1_b3*rsm1_2+8.d0*ec1_b4) 902 d2ecrs1_drs2= 4.0d0*ec1_aa*ec1_a1*ec1_q1p*ec1_den & 903 & -ec1_q0*ec1_q1pp*ec1_den & 904 & +ec1_q0*ec1_q1p**2*ec1_den**2*(2.d0*ec1_q1+one) 905 906 mac_q1pp=half*mac_aa*(-mac_b1*rsm1_2**3+3.d0*mac_b3*rsm1_2+8.d0*mac_b4) 907 d2macrs_drs2= 4.0d0*mac_aa*mac_a1*mac_q1p*mac_den & 908 & -mac_q0*mac_q1pp*mac_den & 909 & +mac_q0*mac_q1p**2*mac_den**2*(2.d0*mac_q1+one) 910 911 d2gcrs_drs2=d2ecrs1_drs2-d2ecrs0_drs2+d2macrs_drs2*fsec_inv 912 fpp_zeta=(zetpm1_3(ipts)**2+zetmm1_3(ipts)**2) * factfpp_zeta 913 d2fzeta4_dzeta2=12.0d0*zeta**2*f_zeta & 914 & + 8.0d0*zeta**3*fp_zeta & 915 & + zeta4 *fpp_zeta 916 917 d2ecrs_drs2=d2ecrs0_drs2+& 918 & f_zeta*(zeta4*d2gcrs_drs2-d2macrs_drs2*fsec_inv) 919 d2ecrs_drsdzeta=dfzeta4_dzeta*dgcrs_drs-fp_zeta*dmacrs_drs*fsec_inv 920 d2ecrs_dzeta2=d2fzeta4_dzeta2*gcrs-fpp_zeta*macrs*fsec_inv 921 922 923 ! At this point, the coding of the PW92 functional finishes. 924 925 ! The correlation between different spin from HCTH is now computed 926 ! First, the part without gradient correction factor 927 rhoecab=ecrs*rhotot-rhoecrs1_dn-rhoecrs1_up 928 vxcadd=ecrs-rs*third*decrs_drs-zeta*decrs_dzeta 929 drhoecab_drhoup=vxcadd+decrs_dzeta-drhoecrs1_drhoup 930 drhoecab_drhodn=vxcadd-decrs_dzeta-drhoecrs1_drhodn 931 932 ! Now, the gradient correction factor 933 ssavg=half*(ss_up+ss_dn) 934 divcab=one/(one+gammacab*ssavg) 935 ucab=gammacab*ssavg*divcab 936 ducabdss=gammacab*divcab*(one-ssavg*gammacab*divcab) 937 938 gcab=ccab0+ucab*(ccab1+ucab*(ccab2+ucab*(ccab3+ucab*ccab4))) 939 dgcabdss=(ccab1+ucab*(two*ccab2+ucab*(three*ccab3+ucab*four*ccab4)))& 940 & *ducabdss 941 942 exc=exc+rhoecab*gcab 943 944 vxci(ipts,1)=vxci(ipts,1)+drhoecab_drhoup*gcab+rhoecab*dgcabdss*half*dssupdrho 945 dvxcdgr(ipts,1)=dvxcdgr(ipts,1)+rhoecab*dgcabdss*half*dssupdg 946 vxci(ipts,2)=vxci(ipts,2)+drhoecab_drhodn*gcab+rhoecab*dgcabdss*half*dssdndrho 947 dvxcdgr(ipts,2)=dvxcdgr(ipts,2)+rhoecab*dgcabdss*half*dssdndg 948 949 ! Final division by the total density, to give the energy density 950 exci(ipts)=exc*rhotot_inv 951 952 end do ! ipts=1,npts 953 954 else 955 956 do ipts=1,npts 957 958 ! ----------------------------------------------------------------------- 959 ! First take care of the spin-split part of the functional 960 exc=zero 961 do ispden=1,nspden 962 rho =rho_updn(ipts,ispden) 963 rhomot=rho_updnm1_3(ipts,ispden) 964 rho_inv=rhomot*rhomot*rhomot 965 966 ! Exchange part 967 ex_lsd= - threefourth_divpi * sixpi2_1_3*rhomot*rhomot*rho 968 ! Note that this definition differs from the PBE one 969 coeffss=rho_inv*rho_inv*rhomot*rhomot 970 ss=grho2_updn(ipts,ispden)*coeffss 971 dssdrho=-eight*third*ss*rho_inv 972 dssdg=two*coeffss 973 974 divx=one/(one+gammax*ss) 975 uxsig=gammax*ss*divx 976 duxsigdss=gammax*divx*(one-ss*gammax*divx) 977 978 gxsig=cxsig0+uxsig*(cxsig1+uxsig*(cxsig2+uxsig*(cxsig3+uxsig*cxsig4))) 979 dgxsigdss=(cxsig1+uxsig*(two*cxsig2+uxsig*(three*cxsig3+uxsig*four*cxsig4)))& 980 & *duxsigdss 981 982 exc=exc+ex_lsd*rho*gxsig 983 vxci(ipts,ispden)=ex_lsd*(four_thirds*gxsig+rho*dgxsigdss*dssdrho) 984 dvxcdgr(ipts,ispden)=ex_lsd*rho*dgxsigdss*dssdg 985 986 ! Spin parallel correlation part 987 ! Note that this definition is for fully spin-polarized quantities 988 rs=rsfac*rhomot 989 rhomo6=sqrt(rhomot) 990 sqr_rs=sq_rsfac*rhomo6 991 rhoo6=rho*rhomot*rhomot*rhomo6 992 rsm1_2=sq_rsfac_inv*rhoo6 993 drsdrho=-third*rs*rho_inv 994 995 ec1_q0=-two*ec1_aa*(one+ec1_a1*rs) 996 ec1_q1=two*ec1_aa*(ec1_b1*sqr_rs+ec1_b2*rs+ec1_b3*rs*sqr_rs+ec1_b4*rs*rs) 997 ec1_q1p=ec1_aa*(ec1_b1*rsm1_2+two*ec1_b2+three*ec1_b3*sqr_rs+four*ec1_b4*rs) 998 ec1_den=one/(ec1_q1*ec1_q1+ec1_q1) 999 ! ec1_log=log( one + one / ec1_q1 ) 1000 ec1_log=-log( ec1_q1*ec1_q1*ec1_den ) 1001 ecrs1=ec1_q0*ec1_log 1002 decrs1_drs= -two*ec1_aa*ec1_a1*ec1_log - ec1_q0*ec1_q1p *ec1_den 1003 decrs1_drho=ecrs1+decrs1_drs*drsdrho*rho 1004 1005 ! Store the LSDA corr energy and density derivative 1006 if(ispden==1)then 1007 rhoecrs1_up=rho*ecrs1 1008 drhoecrs1_drhoup=decrs1_drho 1009 ss_up=ss 1010 dssupdrho=dssdrho 1011 dssupdg=dssdg 1012 else 1013 rhoecrs1_dn=rho*ecrs1 1014 drhoecrs1_drhodn=decrs1_drho 1015 ss_dn=ss 1016 dssdndrho=dssdrho 1017 dssdndg=dssdg 1018 end if 1019 1020 divcsig=one/(one+gammacsig*ss) 1021 ucsig=gammacsig*ss*divcsig 1022 ducsigdss=gammacsig*divcsig*(one-ss*gammacsig*divcsig) 1023 1024 gcsig=ccsig0+ucsig*(ccsig1+ucsig*(ccsig2+ucsig*(ccsig3+ucsig*ccsig4))) 1025 dgcsigdss=(ccsig1+ucsig*(two*ccsig2+ucsig*(three*ccsig3+ucsig*four*ccsig4)))& 1026 & *ducsigdss 1027 1028 ! DEBUG 1029 ! gcsig=zero 1030 ! dgcsigdss=zero 1031 ! ENDDEBUG 1032 exc=exc+ecrs1*rho*gcsig 1033 vxci(ipts,ispden)=vxci(ipts,ispden)+decrs1_drho*gcsig+& 1034 & ecrs1*rho*dgcsigdss*dssdrho 1035 dvxcdgr(ipts,ispden)=dvxcdgr(ipts,ispden)+ecrs1*rho*dgcsigdss*dssdg 1036 1037 end do 1038 1039 ! ----------------------------------------------------------------------------- 1040 ! Then takes care of the LSD correlation part of the functional 1041 1042 rhotot=rhoarr(ipts) 1043 rhotmot=rhom1_3(ipts) 1044 rhotot_inv=rhotmot*rhotmot*rhotmot 1045 rhotmo6=sqrt(rhotmot) 1046 rhoto6=rhotot*rhotmot*rhotmot*rhotmo6 1047 1048 ! From now, the coding of the PW92 functional start. It is identical in xcpbe.f 1049 rs=rsfac*rhotmot 1050 sqr_rs=sq_rsfac*rhotmo6 1051 rsm1_2=sq_rsfac_inv*rhoto6 1052 1053 ! Formulas A6-A8 of PW92LSD 1054 ec0_q0=-2.0d0*ec0_aa*(one+ec0_a1*rs) 1055 ec0_q1=2.0d0*ec0_aa*(ec0_b1*sqr_rs+ec0_b2*rs+ec0_b3*rs*sqr_rs+ec0_b4*rs*rs) 1056 ec0_q1p=ec0_aa*(ec0_b1*rsm1_2+2.d0*ec0_b2+3.d0*ec0_b3*sqr_rs+4.d0*ec0_b4*rs) 1057 ec0_den=one/(ec0_q1*ec0_q1+ec0_q1) 1058 ! ec0_log=log( one + one / ec0_q1 ) 1059 ec0_log=-log( ec0_q1*ec0_q1*ec0_den ) 1060 ecrs0=ec0_q0*ec0_log 1061 decrs0_drs= -2.0d0*ec0_aa*ec0_a1*ec0_log - ec0_q0*ec0_q1p *ec0_den 1062 1063 mac_q0=-2.0d0*mac_aa*(one+mac_a1*rs) 1064 mac_q1=2.0d0*mac_aa*(mac_b1*sqr_rs+mac_b2*rs+mac_b3*rs*sqr_rs+mac_b4*rs*rs) 1065 mac_q1p=mac_aa*(mac_b1*rsm1_2+2.d0*mac_b2+3.d0*mac_b3*sqr_rs+4.d0*mac_b4*rs) 1066 mac_den=one/(mac_q1*mac_q1+mac_q1) 1067 mac_log=-log( mac_q1*mac_q1*mac_den ) 1068 macrs=mac_q0*mac_log 1069 dmacrs_drs= -2.0d0*mac_aa*mac_a1*mac_log - mac_q0*mac_q1p*mac_den 1070 1071 zeta=(rho_updn(ipts,1)-rho_updn(ipts,2))*rhotot_inv 1072 ec1_q0=-two*ec1_aa*(one+ec1_a1*rs) 1073 ec1_q1=two*ec1_aa*(ec1_b1*sqr_rs+ec1_b2*rs+ec1_b3*rs*sqr_rs+ec1_b4*rs*rs) 1074 ec1_q1p=ec1_aa*(ec1_b1*rsm1_2+two*ec1_b2+three*ec1_b3*sqr_rs+four*ec1_b4*rs) 1075 ec1_den=one/(ec1_q1*ec1_q1+ec1_q1) 1076 ! ec1_log=log( one + one / ec1_q1 ) 1077 ec1_log=-log( ec1_q1*ec1_q1*ec1_den ) 1078 ecrs1=ec1_q0*ec1_log 1079 decrs1_drs= -two*ec1_aa*ec1_a1*ec1_log - ec1_q0*ec1_q1p *ec1_den 1080 1081 ! alpha_zeta is introduced in order to remove singularities for fully 1082 ! polarized systems. 1083 zetp_1_3=(one+zeta*alpha_zeta)*zetpm1_3(ipts)**2 1084 zetm_1_3=(one-zeta*alpha_zeta)*zetmm1_3(ipts)**2 1085 1086 f_zeta=( (one+zeta*alpha_zeta2)*zetp_1_3 + & 1087 & (one-zeta*alpha_zeta2)*zetm_1_3 - 2.0d0 ) * factf_zeta 1088 fp_zeta=( zetp_1_3 - zetm_1_3 ) * factfp_zeta 1089 zeta4=zeta**4 1090 gcrs=ecrs1-ecrs0+macrs*fsec_inv 1091 ! ecrs=ecrs0+f_zeta*(-macrs*(one-zeta4)*fsec_inv+(ecrs1-ecrs0)*zeta4) 1092 ecrs=ecrs0+f_zeta*(zeta4*gcrs-macrs*fsec_inv) 1093 1094 dgcrs_drs=decrs1_drs-decrs0_drs+dmacrs_drs*fsec_inv 1095 ! decrs_drs=decrs0_drs+f_zeta*& 1096 ! & (-dmacrs_drs*(one-zeta4)*fsec_inv+(decrs1_drs-decrs0_drs)*zeta4) 1097 decrs_drs=decrs0_drs+f_zeta*(zeta4*dgcrs_drs-dmacrs_drs*fsec_inv) 1098 dfzeta4_dzeta=4.0d0*zeta**3*f_zeta+fp_zeta*zeta4 1099 decrs_dzeta=dfzeta4_dzeta*gcrs-fp_zeta*macrs*fsec_inv 1100 1101 ! At this point, the coding of the PW92 functional finishes. 1102 1103 ! The correlation between different spin from HCTH is now computed 1104 ! First, the part without gradient correction factor 1105 rhoecab=ecrs*rhotot-rhoecrs1_dn-rhoecrs1_up 1106 vxcadd=ecrs-rs*third*decrs_drs-zeta*decrs_dzeta 1107 drhoecab_drhoup=vxcadd+decrs_dzeta-drhoecrs1_drhoup 1108 drhoecab_drhodn=vxcadd-decrs_dzeta-drhoecrs1_drhodn 1109 1110 ! Now, the gradient correction factor 1111 ssavg=half*(ss_up+ss_dn) 1112 divcab=one/(one+gammacab*ssavg) 1113 ucab=gammacab*ssavg*divcab 1114 ducabdss=gammacab*divcab*(one-ssavg*gammacab*divcab) 1115 1116 gcab=ccab0+ucab*(ccab1+ucab*(ccab2+ucab*(ccab3+ucab*ccab4))) 1117 dgcabdss=(ccab1+ucab*(two*ccab2+ucab*(three*ccab3+ucab*four*ccab4)))& 1118 & *ducabdss 1119 1120 exc=exc+rhoecab*gcab 1121 1122 vxci(ipts,1)=vxci(ipts,1)+drhoecab_drhoup*gcab+rhoecab*dgcabdss*half*dssupdrho 1123 dvxcdgr(ipts,1)=dvxcdgr(ipts,1)+rhoecab*dgcabdss*half*dssupdg 1124 vxci(ipts,2)=vxci(ipts,2)+drhoecab_drhodn*gcab+rhoecab*dgcabdss*half*dssdndrho 1125 dvxcdgr(ipts,2)=dvxcdgr(ipts,2)+rhoecab*dgcabdss*half*dssdndg 1126 1127 ! Final division by the total density, to give the energy density 1128 exci(ipts)=exc*rhotot_inv 1129 1130 end do ! ipts=1,npts 1131 1132 end if 1133 1134 else 1135 ! Disallowed value for nspden 1136 write(message, '(3a,i0)' )& 1137 & ' Argument nspden must be 1 or 2; ',ch10,& 1138 & ' Value provided as argument was ',nspden 1139 ABI_BUG(message) 1140 end if 1141 1142 !DEBUG 1143 !Finite-difference debugging, do not take away 1144 !Beware that dvxcdgr(:,3) no longer exists 1145 !if(debug/=0)then 1146 !do ipts=1,npts,5 1147 1148 !rho=rho_updn(ipts,1)+rho_updn(ipts,2) 1149 !write(std_out,'(a,i5,a,es16.8)' ) ' Point number',ipts,' with rho=',rho 1150 1151 !! For rho 1152 !if(debug==1)then 1153 !write(std_out,'(3es16.8)' )exci(ipts)*rho,vxci(ipts,1),vxci(ipts,2) 1154 !else 1155 !! For grho2 1156 !write(std_out,'(4es16.8)' )exci(ipts)*rho,dvxcdgr(ipts,1),& 1157 !& dvxcdgr(ipts,2),dvxcdgr(ipts,3) 1158 !end if 1159 1160 !if(debug==1)then 1161 !! For rho 1162 !write(std_out,'(3es16.8)' )exci(ipts)*rho,& 1163 !& ( exci(ipts+1)*(rho+delta) - exci(ipts+2)*(rho-delta) )/2.d0/delta,& 1164 !& ( exci(ipts+3)*(rho+delta) - exci(ipts+4)*(rho-delta) )/2.d0/delta 1165 !write(std_out,'(3es16.8)' )& 1166 !& ( vxci(ipts+1,1) - vxci(ipts+2,1) )/2.d0/delta,& 1167 !& ( vxci(ipts+3,1) - vxci(ipts+4,1) )/2.d0/delta,& 1168 !& ( vxci(ipts+3,2) - vxci(ipts+4,2) )/2.d0/delta 1169 !write(std_out,'(4es16.8)' )& 1170 !& ( dvxcdgr(ipts+1,1) - dvxcdgr(ipts+2,1) )/2.d0/delta,& 1171 !& ( dvxcdgr(ipts+3,2) - dvxcdgr(ipts+4,2) )/2.d0/delta,& 1172 !& ( dvxcdgr(ipts+1,3) - dvxcdgr(ipts+2,3) )/2.d0/delta,& 1173 !& ( dvxcdgr(ipts+3,3) - dvxcdgr(ipts+4,3) )/2.d0/delta 1174 !else 1175 !! For grho2 (should distinguish exchange and correlation ...) 1176 !grr=sqrt(grho2_updn(ipts,1)) ! Analysis of exchange 1177 !grr=sqrt(grho2_updn(ipts,3)) ! Analysis of correlation 1178 !write(std_out,'(3es16.8)' )exci(ipts)*rho,& 1179 !& ( exci(ipts+1)*rho - exci(ipts+2)*rho )/2.d0/delta/grr,& 1180 !& ( exci(ipts+3)*rho - exci(ipts+4)*rho )/2.d0/delta/grr 1181 !write(std_out,'(3es16.8)' )& 1182 !& ( vxci(ipts+1,1) - vxci(ipts+2,1) )/2.d0/delta/grr,& 1183 !& ( vxci(ipts+3,1) - vxci(ipts+4,1) )/2.d0/delta/grr,& 1184 !& ( vxci(ipts+3,2) - vxci(ipts+4,2) )/2.d0/delta/grr 1185 !write(std_out,'(4es16.8)' )& 1186 !& ( dvxcdgr(ipts+1,1) - dvxcdgr(ipts+2,1) )/2.d0/delta/grr,& 1187 !& ( dvxcdgr(ipts+3,2) - dvxcdgr(ipts+4,2) )/2.d0/delta/grr,& 1188 !& ( dvxcdgr(ipts+1,3) - dvxcdgr(ipts+2,3) )/2.d0/delta/grr,& 1189 !& ( dvxcdgr(ipts+3,3) - dvxcdgr(ipts+4,3) )/2.d0/delta/grr 1190 !end if 1191 !end do 1192 !stop 1193 !end if 1194 !ENDDEBUG 1195 1196 ABI_FREE(rhoarr) 1197 ABI_FREE(rhom1_3) 1198 ABI_FREE(rho_updnm1_3) 1199 if(nspden==2) then 1200 ABI_FREE(zetm) 1201 ABI_FREE(zetmm1_3) 1202 ABI_FREE(zetp) 1203 ABI_FREE(zetpm1_3) 1204 end if 1205 1206 !DEBUG 1207 !write(std_out,*)' xchcth : exit' 1208 !write(std_out,*)' nspden=',nspden 1209 !if(order==2)stop 1210 !ENDDEBUG 1211 1212 end subroutine xchcth