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