1c VS98 exchange functional 2c META GGA 3C utilizes ingredients: 4c rho - density 5c delrho - gradient of density 6c tau - K.S kinetic energy density 7c tauU - uniform-gas KE density 8c ijzy - 1 VS98 9c ijzy - 2 M06-L 10c ijzy - 3 M06-HF 11c ijzy - 4 M06 12c ijzy - 6 revM06-L 13c ijzy - 7 revM06 14c 15c References: 16c 17c [a] T. V. Voorhis and G. E. Scuseria, J. Chem. Phys. 109, 400 (1998). 18c [b] Y. Zhao and D. G. Truhlar, J. Chem. Phys. 125, 194101 (2006). 19 20 21 22 Subroutine xc_xvs98(tol_rho, fac,lfac,nlfac, rho, delrho, 23 & Amat, Cmat, nq, ipol, Ex, 24 & qwght, ldew, func, tau, Mmat, ijzy) 25 26 27c 28c$Id$ 29c 30 implicit none 31c 32#include "errquit.fh" 33c 34 double precision fac, Ex 35 integer nq, ipol 36 logical lfac, nlfac,ldew 37 double precision func(*) ! value of the functional [output] 38c 39c Charge Density & Its Cube Root 40c 41 double precision rho(nq,ipol*(ipol+1)/2) 42c 43c Charge Density Gradient 44c 45 double precision delrho(nq,3,ipol) 46c 47c Quadrature Weights 48c 49 double precision qwght(nq) 50c 51c Sampling Matrices for the XC Potential & Energy 52c 53 double precision Amat(nq,ipol), Cmat(nq,*) 54 55 56 double precision tol_rho, pi 57 58c 59 integer n, ijzy 60 double precision rrho, rho43, rho13, rhoo, rho53, rho83 61 double precision Gamma 62c 63c kinetic energy density or tau 64c 65 double precision tau(nq,ipol), Mmat(nq,*) 66 67 double precision tauN,tauu,DTol 68 double precision Tiny, f13, f43, f53, f83, f113 69 double precision gx, gg, x, z, kx,xk,zk 70 double precision One, Two, Three, Four, Five, Six, Seven, Eight 71 double precision Nine, F10, F11 72 double precision cf, Axlsda, r1, r2, r3, r4, r5, r6 73 74c functional derivatives below FFFFFFFFFFFF 75 76 double precision dxdr, dxdg, dzdr, dzdt, dgdx, dgdz 77 78c functional derivatives above FFFFFFFFFFFF 79 80 81cedo parameter( pi = 3.1415926535897932384626433832795d0 ) 82 83 parameter (cf = 9.115599720d0, Axlsda = -0.9305257363491d0 ) 84 parameter (gg = 0.00186726d0) 85 parameter (f13=1.d0/3.d0,f43=4.0d0/3.0d0,f53=5.0d0/3.0d0) 86 parameter (f83=8.d0/3.0d0, F113=11.0d0/3.d0) 87 parameter (One=1.0d0, Two=2.0d0, Three=3.0d0, Four=4.0d0, 88 & Five=5.0d0,Six=6.0d0, Seven=7.0d0, 89 & Eight=8.0d0, Nine=9.0d0,F10=10.d0, F11=11.d0) 90 pi=acos(-1d0) 91 92 93 if (ijzy.eq.1) then 94c 95c Parameters for VS98 96c 97 r1= -9.800683d-01 98 r2= -3.556788d-03 99 r3= 6.250326d-03 100 r4= -2.354518d-05 101 r5= -1.282732d-04 102 r6= 3.574822d-04 103 elseif (ijzy.eq.2) then 104c 105c Parameters for M06-L 106c 107 r1 = 6.012244D-01*Axlsda 108 r2 = 4.748822D-03*Axlsda 109 r3 = -8.635108D-03*Axlsda 110 r4 = -9.308062D-06*Axlsda 111 r5 = 4.482811D-05*Axlsda 112 r6 = 0.000000D+00 113 elseif (ijzy.eq.3) then 114c 115c Parameters for M06-HF 116c 117 r1 = -1.179732D-01*Axlsda 118 r2 = -2.500000D-03*Axlsda 119 r3 = -1.180065D-02*Axlsda 120 r4 = 0.000000D+00 121 r5 = 0.000000D+00 122 r6 = 0.000000D+00 123 elseif (ijzy.eq.4) then 124c 125c Parameters for M06 126c 127 r1 = 1.422057D-01*Axlsda 128 r2 = 7.370319D-04*Axlsda 129 r3 = -1.601373D-02*Axlsda 130 r4 = 0.000000D+00 131 r5 = 0.000000D+00 132 r6 = 0.000000D+00 133 elseif (ijzy.eq.6) then 134c 135c Parameters for revM06-L 136c 137 r1 = -4.23227252D-01*Axlsda 138 r2 = 0.000000D+00*Axlsda 139 r3 = 3.724234D-03*Axlsda 140 r4 = 0.000000D+00*Axlsda 141 r5 = 0.000000D+00*Axlsda 142 r6 = 0.000000D+00 143 elseif (ijzy.eq.7) then 144c 145c Parameters for revM06 146c 147 r1 = -5.523940140D-02*Axlsda 148 r2 = 0.000000D+00*Axlsda 149 r3 = -3.782631233D-03*Axlsda 150 r4 = 0.000000D+00*Axlsda 151 r5 = 0.000000D+00*Axlsda 152 r6 = 0.000000D+00 153 else 154 call errquit("xc_xvs98: illegal value of ijzy",ijzy,UERR) 155 endif 156 157 DTol = tol_rho 158 159 160c 161 162 163c 164 if (ipol.eq.1 )then 165c 166c ======> SPIN-RESTRICTED <====== 167c or 168c SPIN-UNPOLARIZED 169c 170c 171 do 10 n = 1, nq 172 rhoo = rho(n,1)/Two ! rho_sigma 173 if (rhoo.lt.DTol) goto 10 174 rho43 = rhoo**f43 175 rrho = 1d0/rhoo ! reciprocal of rho 176 rho13 = rho43*rrho 177 rho53 = rhoo**f53 178 rho83 = rho53*rhoo 179c 180 181 tauN = tau(n,1) 182 if(taun.lt.dtol) goto 10 183 tauu=tauN 184 gamma = delrho(n,1,1)*delrho(n,1,1) + 185 & delrho(n,2,1)*delrho(n,2,1) + 186 & delrho(n,3,1)*delrho(n,3,1) 187 if(sqrt(gamma).lt.dtol) goto 10 188 gamma = gamma/Four 189 x = gamma/rho83 190 dxdr = -f83*x*rrho 191 dxdg = One/rho83 192 z = tauu/rho53 - cf 193 dzdr = -f53 * tauu/rho83 194 dzdt = One/rho53 195 kx = One + gg*x + gg*z 196 xk = x/kx 197 zk = z/kx 198 call gvt4(gx,dgdx,dgdz,xk,zk,kx,gg,gg,r1,r2,r3,r4,r5,r6) 199 200 Ex = Ex + Two*rho43*gx*qwght(n) 201 if(ldew) func(n)=func(n)+ Two*rho43*gx 202c 203c functional derivatives 204c 205 Amat(n,1) = Amat(n,1) + f43*rho13*gx + 206 & rho43*(dgdx*dxdr + dgdz*dzdr) 207 Cmat(n,1)= Cmat(n,1) + rho43*(dgdx*dxdg) 208 Mmat(n,1)= Mmat(n,1) + rho43*(dgdz*dzdt) 209 21010 continue 211 212 213c UUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUnrestricted 214 else ! ipol=2 215c 216c ======> SPIN-UNRESTRICTED <====== 217 218 219 do 20 n = 1, nq 220 if (rho(n,1).lt.DTol) goto 20 221c 222c Alpha ALPHA ALPHA 223c 224 if (rho(n,2).lt.DTol) goto 25 225 rhoo = rho(n,2) 226 rho43 = rhoo**f43 227 rrho = 1.0d0/rhoo ! reciprocal of rho 228 rho13 = rho43*rrho 229 rho53 = rhoo**f53 230 rho83 = rho53*rhoo 231 232c 233 234 tauN = tau(n,1) 235 tauu = tauN*Two 236 237 Gamma = delrho(n,1,1)*delrho(n,1,1) + 238 & delrho(n,2,1)*delrho(n,2,1) + 239 & delrho(n,3,1)*delrho(n,3,1) 240 241 x = gamma/rho83 242 dxdr = -f83*x*rrho 243 dxdg = One/rho83 244 z = tauu/rho53 - cf 245 dzdr = -f53 * tauu/rho83 246 dzdt = One/rho53 247 kx = One + gg*x + gg*z 248 xk = x/kx 249 zk = z/kx 250 call gvt4(gx,dgdx,dgdz,xk,zk,kx,gg,gg,r1,r2,r3,r4,r5,r6) 251 Ex = Ex + rho43*gx*qwght(n) 252 if(ldew) func(n)=func(n)+ rho43*gx 253c 254c functional derivatives 255c 256 Amat(n,1) = Amat(n,1) + f43*rho13*gx + 257 & rho43*(dgdx*dxdr + dgdz*dzdr) 258 Cmat(n,1)= Cmat(n,1) + rho43*(dgdx*dxdg) 259 Mmat(n,1)= Mmat(n,1) + rho43*(dgdz*dzdt) 260 261 262c 263c Beta BETA BETA 264c 265 26625 continue 267 268c 269c Beta 270c 271 if (rho(n,3).lt.DTol) goto 20 272 rhoo = rho(n,3) 273 rho43 = rhoo**f43 274 rrho = 1.0d0/rhoo ! reciprocal of rho 275 rho13 = rho43*rrho 276 rho53 = rhoo**f53 277 rho83 = rho53*rhoo 278 279c 280 281 tauN = tau(n,2) 282 tauu = Two*tauN 283 284 285 Gamma = delrho(n,1,2)*delrho(n,1,2) + 286 & delrho(n,2,2)*delrho(n,2,2) + 287 & delrho(n,3,2)*delrho(n,3,2) 288 289 x = gamma/rho83 290 dxdr = -f83*x*rrho 291 dxdg = One/rho83 292 z = tauu/rho53 - cf 293 dzdr = -f53 * tauu/rho83 294 dzdt = One/rho53 295 kx = One + gg*x + gg*z 296 xk = x/kx 297 zk = z/kx 298 call gvt4(gx,dgdx,dgdz,xk,zk,kx,gg,gg,r1,r2,r3,r4,r5,r6) 299 300 Ex = Ex + rho43*gx*qwght(n) 301 if(ldew) func(n)=func(n)+ rho43*gx 302c 303c functional derivatives 304c 305 Amat(n,2) = Amat(n,2) + f43*rho13*gx + 306 & rho43*(dgdx*dxdr + dgdz*dzdr) 307 Cmat(n,3)= Cmat(n,3) + rho43*(dgdx*dxdg) 308 Mmat(n,2)= Mmat(n,2) + rho43*(dgdz*dzdt) 309 310c 31120 continue 312 endif 313c 314 return 315 end 316 317 318 Subroutine xc_xvs98_d2() 319 implicit none 320 call errquit(' xvs98: d2 not coded ',0,0) 321 return 322 end 323 324 Subroutine gvt4(g,dgdx,dgdz,xk,zk,k,c,ct,r1,r2,r3,r4,r5,r6) 325 Implicit none 326c 327c Evaluate the GVT4 form in the VS98 functional 328c 329c 330c some working variables 331 double precision g,dgdx,dgdz,xk,zk,k,c,ct,r1,r2,r3,r4,r5,r6 332 double precision sxk,szk,sk,sc,sct,sr1,sr2,sr3,sr4,sr5,sr6,sk2 333 double precision One, Two, Three, Four, Six 334 Data One/1.0d0/, Two/2.0d0/, Three/3.0d0/, Four/4.0d0/, Six/6.0d0/ 335C 336 sxk = xk 337 szk = zk 338 sk = k 339 sc = c 340 sct = ct 341 sr1 = r1 342 sr2 = r2 343 sr3 = r3 344 sr4 = r4 345 sr5 = r5 346 sr6 = r6 347 sk2 = sk*sk 348 g = (sr1 + sr2*sxk + sr3*szk 349 $ +sr4*sxk*sxk + sr5*szk*sxk + sr6*szk*szk)/sk 350 dgdx = (-sr1*sc 351 $ +sr2*(One-Two*sc*sxk) 352 $ -Two*sr3*szk*sc 353 $ +sr4*(Two*sxk-Three*sxk*sxk*sc) 354 $ +sr5*(szk -Three*szk*sxk*sc) 355 $ -Three*sr6*szk*szk*sc )/sk2 356 dgdz = (-sr1*sct 357 $ -Two*sr2*sxk*sct 358 $ +sr3*(One-Two*szk*sct) 359 $ -Three*sr4*sxk*sxk*sct 360 $ +sr5*(sxk-Three*sxk*szk*sct) 361 $ +sr6*(Two*szk-Three*szk*szk*sct))/sk2 362 363 return 364 end 365 366 367