1#include "dft2drv.fh" 2c VSXC correlation functional 3c META GGA 4C utilizes ingredients: 5c rho - density 6c delrho - gradient of density 7c tau (tauN)- K.S kinetic energy density 8c ijzy - 1 VS98 9c ijzy - 2 M06-L 10c ijzy - 3 M06-HF 11c ijzy - 4 M06 12c ijzy - 5 M06-2X 13c ijzy - 6 revM06 14c ijzy - 7 revM06-L 15c 16 Subroutine xc_cvs98(tol_rho, cfac, lcfac, nlcfac, rho, delrho, 17 & nq, ipol, Ec, qwght, ldew, func, 18 & tau, Amat, Cmat, Mmat, ijzy) 19 20 21c 22c$Id$ 23c 24c Reference 25c [a] T. V. Voorhis and G. E. Scuseria, J. Chem. Phys. 109, 400 (1998). 26c [b] Y. Zhao and D. G. Truhlar, J. Chem. Phys. 125, 194101 (2006). 27 28c 29 implicit none 30c 31#include "errquit.fh" 32c 33c 34c Input and other parameters 35c 36 integer ipol, nq 37 38 double precision cfac 39 logical lcfac, nlcfac 40 41 logical lfac, nlfac 42 double precision fac 43 double precision tol_rho 44c 45 logical ldew 46 double precision func(*) 47c 48c Threshold parameters 49c 50 double precision DTol,F1, F2, F3, F4, gab, cf 51 Data F1/1.0d0/,F2/2.0d0/, 52 & F3/3.0d0/,F4/4.0d0/,gab/0.00304966d0/, 53 & cf/9.115599720d0/ 54c 55c Correlation energy 56c 57 double precision Ec 58c 59c Charge Density 60c 61 double precision rho(nq,ipol*(ipol+1)/2) 62c 63c Charge Density Gradient 64c 65 double precision delrho(nq,3,ipol) 66 67c 68c Kinetic Energy Density 69c 70 double precision tau(nq,ipol), tauN 71 72c 73c Quadrature Weights 74c 75 double precision qwght(nq) 76c 77c Sampling Matrices for the XC Potential 78c 79 double precision Amat(nq,ipol), Cmat(nq,*) 80 double precision Mmat(nq,*) 81 82 integer n, ijzy 83 84c call to the vs98css subroutine 85 double precision PA,GAA,TA,FA,FPA,FGA,FTA,EUA,EUEGA,ChiA,EUPA 86 &,ChiAP,ChiAG,ZA,ZAP,ZAT 87 double precision PB,GBB,TB,FB,FPB,FGB,FTB,EUB,EUEGB,ChiB,EUPB 88 &,ChiBP,ChiBG,ZB,ZBP,ZBT 89c 90 double precision Pi, F43, F13, Pi34, F6, PotLC, 91 &RS,RSP,Zeta,dZdA,dZdB,dLdS,dLdZ,d2LdSS,d2LdSZ,d2LdZZ 92 double precision P, EUEG, ZAB, XAB, kab, xk, zk 93 double precision dgdx,dgdz,dgdPA,dgdGA,dgdTA,dgdPB,dgdGB,dgdTB 94 double precision EUEGPA,EUEGPB,gcab 95 double precision r7, r8, r9, r10, r11, r12 96 97 98c 99c ======> BOTH SPIN-RESTRICETED AND UNRESTRICTED <====== 100c 101c DTol=1.0d-7 102 dtol=tol_rho 103C Parameters for VS98 104 if (ijzy.eq.1) then 105 r7= 7.035010d-01 106 r8= 7.694574d-03 107 r9= 5.152765d-02 108 r10= 3.394308d-05 109 r11= -1.269420d-03 110 r12= 1.296118d-03 111C Parameters for M06-L 112 elseif (ijzy.eq.2) then 113 r7= 3.957626D-01 114 r8= -5.614546D-01 115 r9= 1.403963D-02 116 r10= 9.831442D-04 117 r11= -3.577176D-03 118 r12= 0.000000D+00 119C Parameters for M06-HF 120 elseif (ijzy.eq.3) then 121 r7= -6.746338D-01 122 r8= -1.534002D-01 123 r9= -9.021521D-02 124 r10= -1.292037D-03 125 r11= -2.352983D-04 126 r12= 0.000000D+00 127 128C Parameters for M06 129 elseif (ijzy.eq.4) then 130 r7= -2.741539D+00 131 r8= -6.720113D-01 132 r9= -7.932688D-02 133 r10=1.918681D-03 134 r11=-2.032902D-03 135 r12=0.000000D+00 136 137C Parameters for M06-2X 138 elseif (ijzy.eq.5) then 139 r7= 1.166404D-01 140 r8= -9.120847D-02 141 r9= -6.726189D-02 142 r10= 6.720580D-05 143 r11= 8.448011D-04 144 r12= 0.000000D+00 145C Parameters for revM06-L 146 elseif (ijzy.eq.6) then 147 r7= 4.007146D-01 148 r8= 1.5796569D-02 149 r9= -3.2680984D-02 150 r10= 0.000000D+00 151 r11= 0.000000D+00 152 r12= 1.260132D-03 153C Parameters for revM06 154 elseif (ijzy.eq.7) then 155 r7= -3.390666720D-01 156 r8= 3.790156384D-03 157 r9= -2.762485975D-02 158 r10= 0.000000D+00 159 r11= 0.000000D+00 160 r12= 4.076285162D-04 161 else 162 call errquit("xc_cvs98: illegal value of ijzy",ijzy,UERR) 163 endif 164 Pi = F4*ATan(F1) 165 F6=6.0d0 166 F43 = F4 / F3 167 Pi34 = F3 / (F4*Pi) 168 F13 = F1 / F3 169 170 171 do 20 n = 1, nq 172 if (rho(n,1).lt.DTol) goto 20 173 if (ipol.eq.1) then 174c 175c get the density, gradient, and tau for the alpha spin from the total 176c 177 PA = rho(n,1)/F2 178 GAA = ( delrho(n,1,1)*delrho(n,1,1) + 179 & delrho(n,2,1)*delrho(n,2,1) + 180 & delrho(n,3,1)*delrho(n,3,1))/4.0d0 181c In the bc95css subroutine, we use 2*TA as the tau, so we do not divide 182c the tau by 2 here 183 184 TA = tau(n,1) 185 186 Call vs98ss(tol_rho,PA,GAA,TA,FA,FPA,FGA,FTA,EUA,ZA, 187 & ChiA,EUPA,ChiAP,ChiAG,ZAP,ZAT,ijzy) 188 PB = PA 189 GBB = GAA 190 TB = TA 191 FB = FA 192 FPB = FPA 193 FGB = FGA 194 FTB = FTA 195 EUB = EUA 196 ZB = ZA 197 ChiB = ChiA 198 EUPB = EUPA 199 ChiBP = ChiAP 200 ChiBG = ChiAG 201 ZBP = ZAP 202 ZBT = ZAT 203 204 Ec = Ec + 2.0d0*FA*qwght(n) !factor of 2 account for both spin 205 if(ldew) func(n)=func(n)+ 2.0d0*FA 206 Amat(n,1)=Amat(n,1)+ FPA 207 Cmat(n,D1_GAA)= Cmat(n,D1_GAA) + FGA 208 Mmat(n,1)= Mmat(n,1) + FTA 209c if 0 210c write (0,'(A,3F20.6)') " PA,EUA",PA,EUA 211c endif 212 213 214c UUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUnrestricted 215 else ! ipol=2 216c 217c ======> SPIN-UNRESTRICTED <====== 218c 219c 220c alpha 221c 222 223 PA = rho(n,2) 224 if (PA.le.DTol) go to 25 225 GAA = delrho(n,1,1)*delrho(n,1,1) + 226 & delrho(n,2,1)*delrho(n,2,1) + 227 & delrho(n,3,1)*delrho(n,3,1) 228c 229c In the bc95css subroutine, we use 2*TA as the tau 230c 231 TA = tau(n,1)*2.0d0 232 233 Call vs98ss(tol_rho,PA,GAA,TA,FA,FPA,FGA,FTA,EUA,ZA, 234 & ChiA,EUPA,ChiAP,ChiAG,ZAP,ZAT,ijzy) 235 Ec = Ec + FA*qwght(n) 236 if(ldew) func(n)=func(n)+ FA 237 Amat(n,1)=Amat(n,1)+ FPA 238 Cmat(n,D1_GAA)= Cmat(n,D1_GAA) + FGA 239c 2*0.5=1.0 for Mmat 240 Mmat(n,1)= Mmat(n,1) + FTA 241#if 0 242 write (0,'(A,3F20.6)') "AAmat Cmat Mmat",FPA,FGA,FTA 243#endif 244 245c 246c In the vs98ss subroutine, we use 2*TA as the tau, 247c 248c 249c Beta 250c 251 25 continue 252 PB = rho(n,3) 253 GBB = delrho(n,1,2)*delrho(n,1,2) + 254 & delrho(n,2,2)*delrho(n,2,2) + 255 & delrho(n,3,2)*delrho(n,3,2) 256 257 TB = tau(n,2)*2.0d0 258 259 Call vs98ss(tol_rho,PB,GBB,TB,FB,FPB,FGB,FTB,EUB,ZB, 260 & ChiB,EUPB,ChiBP,ChiBG,ZBP,ZBT,ijzy) 261 Ec = Ec + FB*qwght(n) 262 if(ldew) func(n)=func(n)+ FB 263 Amat(n,2)= Amat(n,2)+ FPB 264 Cmat(n,D1_GBB)= Cmat(n,D1_GBB) + FGB 265 Mmat(n,2)= Mmat(n,2) + FTB 266 267#if 0 268 write (0,'(A,3F20.6)') "BAmat Cmat Mmat",FPB,FGB,FTB 269#endif 270 endif 271 30 continue 272 P = rho(n,1) 273 If(PA.gt.DTol.and.PB.gt.DTol) then 274 RS = (Pi34/P) ** F13 275 RSP = -RS/(F3*P) 276 Zeta = (PA-PB)/P 277 dZdA = (F1-Zeta)/P 278 dZdB = (-F1-Zeta)/P 279 Call lsdac(tol_rho, 280 A RS,Zeta,PotLC,dLdS,dLdZ,d2LdSS,d2LdSZ, 281 $ d2LdZZ) 282 EUEG = P*PotLC - EUA - EUB 283 ZAB = ZA + ZB 284 XAB = ChiA+ChiB 285 kab = F1 + gab*(XAB+ZAB) 286 xk = XAB/kab 287 zk = ZAB/kab 288 call gvt4(gcab,dgdx,dgdz,xk,zk,kab,gab,gab,r7,r8,r9,r10,r11,r12) 289 Ec = Ec + gcab*EUEG*qwght(n) 290 if(ldew) func(n)=func(n)+ gcab*EUEG 291 dgdPA = dgdx*ChiAP + dgdz*ZAP 292 dgdGA = dgdx*ChiAG 293 dgdTA = dgdz*ZAT 294 dgdPB = dgdx*ChiBP + dgdz*ZBP 295 dgdGB = dgdx*ChiBG 296 dgdTB = dgdz*ZBT 297 EUEGPA = PotLC + P*dLdS*RSP + P*dLdZ*dZdA - EUPA 298 EUEGPB = PotLC + P*dLdS*RSP + P*dLdZ*dZdB - EUPB 299 if (ipol.eq.1) then 300 Amat(n,1) = Amat(n,1) + (EUEGPA*gcab + EUEG*dgdPA) 301 Cmat(n,D1_GAA)= Cmat(n,D1_GAA) + EUEG*dgdGA 302 Mmat(n,1)= Mmat(n,1) + EUEG*dgdTA 303 else 304 Amat(n,1) = Amat(n,1) + (EUEGPA*gcab + EUEG*dgdPA) 305 Amat(n,2) = Amat(n,2) + (EUEGPB*gcab + EUEG*dgdPB) 306 Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + EUEG*dgdGA 307 Cmat(n,D1_GBB) = Cmat(n,D1_GBB) + EUEG*dgdGB 308 Mmat(n,1)= Mmat(n,1) + EUEG*dgdTA 309 Mmat(n,2)= Mmat(n,2) + EUEG*dgdTB 310 endif 311 endIf 312c write (*,*) "Amat(n,1),Cmat(n,1),Mmat(n,1)",Amat(n,1),Cmat(n,1) 313c & ,Mmat(n,1),ipol 314c stop 31520 continue 316 end 317 318 Subroutine xc_cvs98_d2() 319 call errquit(' cvs98: d2 not coded ',0,0) 320 return 321 end 322 323 Subroutine vs98ss(tol_rho,PX,GX,TX,F,FP,FG,FT,EUEG,Z,Chi,EUEGP, 324 & ChiP,ChiG,ZP,ZT,ijzy) 325 Implicit none 326c 327#include "errquit.fh" 328C 329C Compute the same-spin part of the vs98 correlation functional for one grid 330C point and one spin-case. 331C 332 333 integer ijzy 334 double precision tol_rho 335 double precision r13, r14, r15, r16, r17, r18 336 double precision PX, GX, TX, F, FP, FG, FT, DTol, Z, ZP, ZT 337 double precision EUEG, Chi, EUEGP, ChiP, ChiG, cf, gcc 338 double precision Zero, Pt25, F1, F2, F3, F4, F5, F6, F8, F11 339 double precision Pi, Pi34, F13, F23, F43, F53, F83, F113 340 double precision RS, D, RSP, PotLC, DX, DZ, dgdP, dgdG, dgdT 341 double precision E,DP, DG, DT, rhoo, rho43, rho53, rho83 342 double precision rrho, F4o3, rho13, kc, xk, zk, gc, dgdx, dgdz 343 double precision d2LdSS, d2LdSZ, d2LdZZ, dLdS, dLdZ 344 345 Data Zero/0.0d0/, Pt25/0.25d0/, F1/1.0d0/, F2/2.0d0/, F3/3.0d0/, 346 $ F4/4.0d0/, F5/5.0d0/, F6/6.0d0/, F8/8.0d0/, F11/11.0d0/, 347 $ gcc/0.00515088d0/,cf/9.115599720d0/ 348 349 350 F4o3 = 4.0d0/3.0d0 351C Parameters for VS98 352 if (ijzy.eq.1) then 353 r13= 3.270912d-01 354 r14= -3.228915d-02 355 r15= -2.942406d-02 356 r16= 2.134222d-03 357 r17= -5.451559d-03 358 r18= 1.577575d-02 359C Parameters for M06-L 360 elseif (ijzy.eq.2) then 361 r13= 4.650534D-01 362 r14= 1.617589D-01 363 r15= 1.833657D-01 364 r16= 4.692100D-04 365 r17= -4.990573D-03 366 r18= 0.000000D+00 367C Parameters for M06-HF 368 elseif (ijzy.eq.3) then 369 r13= 8.976746D-01 370 r14= -2.345830D-01 371 r15= 2.368173D-01 372 r16= -9.913890D-04 373 r17= -1.146165D-02 374 r18= 0.000000D+00 375C Parameters for M06 376 elseif (ijzy.eq.4) then 377 r13= 4.905945D-01 378 r14= -1.437348D-01 379 r15= 2.357824D-01 380 r16= 1.871015D-03 381 r17= -3.788963D-03 382 r18= 0.000000D+00 383C Parameters for M06-2X 384 elseif (ijzy.eq.5) then 385 r13= 6.902145D-01 386 r14= 9.847204D-02 387 r15= 2.214797D-01 388 r16= -1.968264D-03 389 r17= -6.775479D-03 390 r18= 0.000000D+00 391C Parameters for revM06-L 392 elseif (ijzy.eq.6) then 393 r13= -5.38821292D-01 394 r14= -2.829603D-02 395 r15= 2.3889696D-02 396 r16= 0.000000D+00 397 r17= 0.000000D+00 398 r18= -2.437902D-03 399C Parameters for revM06 400 elseif (ijzy.eq.7) then 401 r13= -1.467095900D-01 402 r14= -1.832187007D-04 403 r15= 8.484372430D-02 404 r16= 0.000000D+00 405 r17= 0.000000D+00 406 r18= 2.280677172D-04 407 else 408 call errquit("vs98ss: illegal value of ijzy",ijzy,UERR) 409 endif 410 411 dtol=tol_rho 412 If(PX.le.DTol.or.gx.le.dtol.or.tx.le.dtol) then 413 EUEG = Zero 414 Chi = Zero 415 EUEGP = Zero 416 ChiP = Zero 417 ChiG = Zero 418 PX = Zero 419 GX = Zero 420 TX = Zero 421 F = Zero 422 FP = Zero 423 FG = Zero 424 FT = Zero 425 Z = Zero 426 ZP = Zero 427 ZT = Zero 428 else 429 Pi = F4*ATan(F1) 430 Pi34 = F3 / (F4*Pi) 431 F13 = F1 / F3 432 F23 = F2 / F3 433 F43 = F2 * F23 434 F53 = F5 / F3 435 F83 = F8 / F3 436 F113 = F11 / F3 437 rhoo = PX 438 rrho = 1.0d0/rhoo 439 rho43 = rhoo**F4o3 440 rho13 = rho43*rrho 441 rho53 = rhoo**F53 442 rho83 = rho53*rhoo 443 444 RS = (Pi34/PX) ** F13 445 Call lsdac(tol_rho, 446 A RS,F1,PotLC,dLdS,dLdZ,d2LdSS,d2LdSZ,d2LdZZ) 447 EUEG = PX*PotLC 448 Chi = GX/rho83 449 Z = (TX/rho53) - cf 450 kc = F1 + gcc*(Chi + Z) 451 xk = Chi/kc 452 zk = Z/kc 453 D = F1 - Chi/(F4*(Z + cf)) 454 call gvt4(gc,dgdx,dgdz,xk,zk,kc,gcc,gcc,r13,r14,r15,r16,r17,r18) 455 E = D*EUEG*gc 456c write (*,*) "Chi, Z, gc", CHi, Z, gc 457 F = E 458c 459 RSP = -RS/(F3*Px) 460 ChiG = F1/PX**F83 461 ChiP = -F83*Chi/PX 462 ZP = -F53 * TX/rho83 463 ZT = F1/rho53 464 DZ = Chi/(F4*(Z + cf)*(Z + cf)) 465 DX = -F1/(F4*(Z + cf)) 466 DP = DZ*ZP + DX*ChiP 467 DG = DX*ChiG 468 DT = DZ*ZT 469 dgdP = dgdx*ChiP + dgdz*ZP 470 dgdG = dgdx*ChiG 471 dgdT = dgdz*ZT 472 EUEGP = PotLC + PX*dLdS*RSP 473 FP = DP*EUEG*gc + D*EUEGP*gc + D*EUEG*dgdP 474 FG = DG*EUEG*gc + D*EUEG*dgdG 475 FT = DT*EUEG*gc + D*EUEG*dgdT 476 Endif 477 Return 478 End 479 480 481