1 Subroutine xc_xm11(tol_rho, fac,lfac,nlfac, rho, delrho, 2 & Amat, Cmat, nq, ipol, Ex, 3 & qwght, ldew, func, tau, Mmat, ijzy) 4 5 6c 7c$Id$ 8c 9c 10c 11c**********************************************************************c 12c c 13c xc_xm11 evaluates the exchange part of the M08 and M11 suite of c 14c functionals on the grid. c 15c !!! Second derivatives are not available yet. c 16c c 17c Ref: (a) Zhao, Y. and Truhlar, D. G. JCTC, 2008, 4 , 1849 c 18c (b) Peverati, R. and Truhlar, D. G. J.P.C.Lett. 2011, 2, 2810 c 19c (c) Peverati, R. and Truhlar, D. G. J.P.C.Lett. 2011, 3, 117 c 20c c 21c ijzy - 1 M08-HX (a) c 22c ijzy - 2 M08-SO (a) c 23c ijzy - 3 M11 (b) c 24c ijzy - 4 M11-L (c) c 25c c 26c Coded by Roberto Peverati (12/11) c 27c c 28c**********************************************************************c 29c 30 implicit none 31c 32#include "errquit.fh" 33c 34 double precision fac, Ex 35 integer nq, ipol 36 logical lfac, nlfac,ldew, uselc 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,*), Mmat(nq,*) 54 double precision tol_rho, pi 55c 56c kinetic energy density or tau 57c 58 double precision tau(nq,ipol) 59 double precision tauN,tauu,DTol 60c 61c functional derivatives 62c 63 double precision dWdT, dTdR, dTdTau 64c 65c Intermediate derivative results, etc. 66c 67 integer n, ijzy 68c 69 double precision Ax, s, s2 70 double precision kapa,kapas,mu,mus 71c 72 double precision f0,f1,f2,f3,f4,f5,f6,f7,f8,f9,f10,f11 73 double precision F1o3,F2o3,F3o5,F4o3,F5o3,F48,F81 74 double precision Fsig1, Fsig2, Fsig3, Fsig4, Fx1, Fx2 75 double precision ElSR, ElLR 76 double precision PDUM 77 double precision GGA1, GGA2, GGA3, GGA4 78 double precision Emu, X, deno 79 double precision ds2drho, ds2dg, dfx1ds2 80 double precision dfx2ds2, df1dw 81 double precision dfx1drho,dfx1dg,dfx2drho,dfx2dg,df2dw 82 double precision df3dw, df4dw, delsrdr, dellrdr 83 double precision dgga1dr, dgga2dr, dgga3dr, dgga4dr 84 double precision df1dr, df1dtau, df2dr, df2dtau 85 double precision df3dr, df3dtau, df4dr, df4dtau 86 double precision dgga1dg, dgga2dg, dgga3dg, dgga4dg 87c 88 double precision at00,at01,at02,at03,at04,at05,at06 89 double precision at07,at08,at09,at10,at11 90 double precision bt00,bt01,bt02,bt03,bt04,bt05,bt06 91 double precision bt07,bt08,bt09,bt10,bt11 92 double precision ct00,ct01,ct02,ct03,ct04,ct05,ct06 93 double precision ct07,ct08,ct09,ct10,ct11 94 double precision dt00,dt01,dt02,dt03,dt04,dt05,dt06 95 double precision dt07,dt08,dt09,dt10,dt11 96 double precision rho43, rho13, rhoo, rho53 97 double precision Gamma2, Gamma 98 double precision TauUEG, Tsig, Wsig, W1, W2, W3, W4, W5, W6 99 double precision W7, W8, W9, W10, W11 100c 101 parameter( F0=0.0D+00, F1=1.0D+00, F2=2.0D+00, 102 $ F3=3.0D+00, F4=4.0D+00, F5=5.0D+00, 103 $ F6=6.0D+00, F7=7.0D+00, F8=8.0D+00, 104 $ F9=9.0D+00, F10=10.0D+00,F11=11.0D+00) 105c 106 pi=acos(-1d0) 107c 108 ct00= 0D+00 109 ct01= 0D+00 110 ct02= 0D+00 111 ct03= 0D+00 112 ct04= 0D+00 113 ct05= 0D+00 114 ct06= 0D+00 115 ct07= 0D+00 116 ct08= 0D+00 117 ct09= 0D+00 118 ct10= 0D+00 119 ct11= 0D+00 120C 121 dt00= 0D+00 122 dt01= 0D+00 123 dt02= 0D+00 124 dt03= 0D+00 125 dt04= 0D+00 126 dt05= 0D+00 127 dt06= 0D+00 128 dt07= 0D+00 129 dt08= 0D+00 130 dt09= 0D+00 131 dt10= 0D+00 132 dt11= 0D+00 133 at00=0.000000D+00 134 at01=0.000000D+00 135 at02=0.000000D+00 136 at03=0.000000D+00 137 at04=0.000000D+00 138 at05=0.000000D+00 139 at06=0.000000D+00 140 at07=0.000000D+00 141 at08=0.000000D+00 142 at09=0.000000D+00 143 at10=0.000000D+00 144 at11=0.000000D+00 145 bt00=0.000000D+00 146 bt01=0.000000D+00 147 bt02=0.000000D+00 148 bt03=0.000000D+00 149 bt04=0.000000D+00 150 bt05=0.000000D+00 151 bt06=0.000000D+00 152 bt07=0.000000D+00 153 bt08=0.000000D+00 154 bt09=0.000000D+00 155 bt10=0.000000D+00 156 bt11=0.000000D+00 157 UseLC=.False. 158 Emu =0.00D+00 159C 160 if (ijzy.eq.1) then 161C Parameters for M08-HX 162 at00= 1.3340172D+00 163 at01= -9.4751087D+00 164 at02= -1.2541893D+01 165 at03= 9.1369974D+00 166 at04= 3.4717204D+01 167 at05= 5.8831807D+01 168 at06= 7.1369574D+01 169 at07= 2.3312961D+01 170 at08= 4.8314679D+00 171 at09= -6.5044167D+00 172 at10= -1.4058265D+01 173 at11= 1.2880570D+01 174 175 bt00= -8.5631823D-01 176 bt01= 9.2810354D+00 177 bt02= 1.2260749D+01 178 bt03= -5.5189665D+00 179 bt04= -3.5534989D+01 180 bt05= -8.2049996D+01 181 bt06= -6.8586558D+01 182 bt07= 3.6085694D+01 183 bt08= -9.3740983D+00 184 bt09= -5.9731688D+01 185 bt10= 1.6587868D+01 186 bt11= 1.3993203D+01 187C 188 UseLC=.False. 189C 190 elseif (ijzy.eq.2) then 191C Parameters for M08-SO 192 at00= -3.4888428D-01 193 at01= -5.8157416D+00 194 at02= 3.7550810D+01 195 at03= 6.3727406D+01 196 at04= -5.3742313D+01 197 at05= -9.8595529D+01 198 at06= 1.6282216D+01 199 at07= 1.7513468D+01 200 at08= -6.7627553D+00 201 at09= 1.1106658D+01 202 at10= 1.5663545D+00 203 at11= 8.7603470D+00 204 205 bt00= 7.8098428D-01 206 bt01= 5.4538178D+00 207 bt02= -3.7853348D+01 208 bt03= -6.2295080D+01 209 bt04= 4.6713254D+01 210 bt05= 8.7321376D+01 211 bt06= 1.6053446D+01 212 bt07= 2.0126920D+01 213 bt08= -4.0343695D+01 214 bt09= -5.8577565D+01 215 bt10= 2.0890272D+01 216 bt11= 1.0946903D+01 217C 218 UseLC=.False. 219C 220 elseif (ijzy.eq.3) then 221C Parameters for M11 222 at00= -0.18399900D+00 223 at01= -1.39046703D+01 224 at02= 1.18206837D+01 225 at03= 3.10098465D+01 226 at04= -5.19625696D+01 227 at05= 1.55750312D+01 228 at06= -6.94775730D+00 229 at07= -1.58465014D+02 230 at08= -1.48447565D+00 231 at09= 5.51042124D+01 232 at10= -1.34714184D+01 233 at11= 0.00000000D+00 234 235 bt00= 0.75599900D+00 236 bt01= 1.37137944D+01 237 bt02= -1.27998304D+01 238 bt03= -2.93428814D+01 239 bt04= 5.91075674D+01 240 bt05= -2.27604866D+01 241 bt06= -1.02769340D+01 242 bt07= 1.64752731D+02 243 bt08= 1.85349258D+01 244 bt09= -5.56825639D+01 245 bt10= 7.47980859D+00 246 bt11= 0.00000000D+00 247C 248 UseLC=.True. 249 Emu =0.25D+00 250C 251 elseif (ijzy.eq.4) then 252C Parameters for M11-L 253 at00= 8.121131D-01 254 at01= 1.738124D+01 255 at02= 1.154007D+00 256 at03= 6.869556D+01 257 at04= 1.016864D+02 258 at05= -5.887467D+00 259 at06= 4.517409D+01 260 at07= -2.773149D+00 261 at08= -2.617211D+01 262 at09= 0.000000D+00 263 at10= 0.000000D+00 264 at11= 0.000000D+00 265C 266 bt00= 1.878869D-01 267 bt01= -1.653877D+01 268 bt02= 6.755753D-01 269 bt03= -7.567572D+01 270 bt04= -1.040272D+02 271 bt05= 1.831853D+01 272 bt06= -5.573352D+01 273 bt07= -3.520210D+00 274 bt08= 3.724276D+01 275 bt09= 0.000000D+00 276 bt10= 0.000000D+00 277 bt11= 0.000000D+00 278C 279 ct00= -4.386615D-01 280 ct01= -1.214016D+02 281 ct02= -1.393573D+02 282 ct03= -2.046649D+00 283 ct04= 2.804098D+01 284 ct05= -1.312258D+01 285 ct06= -6.361819D+00 286 ct07= -8.055758D-01 287 ct08= 3.736551D+00 288 ct09= 0.000000D+00 289 ct10= 0.000000D+00 290 ct11= 0.000000D+00 291C 292 dt00= 1.438662D+00 293 dt01= 1.209465D+02 294 dt02= 1.328252D+02 295 dt03= 1.296355D+01 296 dt04= 5.854866D+00 297 dt05= -3.378162D+00 298 dt06= -4.423393D+01 299 dt07= 6.844475D+00 300 dt08= 1.949541D+01 301 dt09= 0.000000D+00 302 dt10= 0.000000D+00 303 dt11= 0.000000D+00 304C 305 UseLC=.True. 306 Emu =0.25D+00 307C 308 else 309 call errquit("xc_xm11: illegal value of ijzy",ijzy,UERR) 310 endif 311 DTol=tol_rho 312 F1o3 = F1/F3 313 F2o3 = F2/F3 314 F3o5 = F3/F5 315 F4o3 = F4/F3 316 F5o3 = F5/F3 317 F48 = 48.0d0 318 F81 = 81.0d0 319 Ax = -(F3/F2) * (F4o3*Pi)**(-F1o3) 320C RPBE parameters 321 Mus = F10/F81 322 kapas = 0.552d0 323C PBE parameters 324 Mu = 0.21951d0 325 kapa = 0.804d0 326c 327 if (ipol.eq.1 )then 328c 329c ======> SPIN-RESTRICTED <====== 330c or 331c SPIN-UNPOLARIZED 332c 333c 334 do 10 n = 1, nq 335 if (rho(n,1).lt.DTol) goto 10 336 rhoo = rho(n,1)/F2 337 rho43 = rhoo**F4o3 338 rho13 = rho43/rhoo 339 rho53 = rhoo**F5o3 340c 341 tauN = tau(n,1) 342 if(taun.lt.dtol) goto 10 343 tauu=tauN 344 TAUUEG=F3O5*((F6*PI*PI)**F2O3)*RHO53 345 Tsig =TauUEG/tauu 346 Wsig =(Tsig - F1)/(Tsig + F1) 347 W1=Wsig 348 W2=Wsig*W1 349 W3=Wsig*W2 350 W4=Wsig*W3 351 W5=Wsig*W4 352 W6=Wsig*W5 353 W7=Wsig*W6 354 W8=Wsig*W7 355 W9=Wsig*W8 356 W10=Wsig*W9 357 W11=Wsig*W10 358 Fsig1 =(at00 + at01*W1 + at02*W2 + at03*W3 359 $ + at04*W4 + at05*W5 + at06*W6 + at07*W7 360 $ + at08*W8 + at09*W9 + at10*W10+ at11*W11) 361 Fsig2 =(bt00 + bt01*W1 + bt02*W2 + bt03*W3 362 $ + bt04*W4 + bt05*W5 + bt06*W6 + bt07*W7 363 $ + bt08*W8 + bt09*W9 + bt10*W10+ bt11*W11) 364 Fsig3 =(ct00 + ct01*W1 + ct02*W2 + ct03*W3 365 $ + ct04*W4 + ct05*W5 + ct06*W6 + ct07*W7 366 $ + ct08*W8 + ct09*W9 + ct10*W10+ ct11*W11) 367 Fsig4 =(dt00 + dt01*W1 + dt02*W2 + dt03*W3 368 $ + dt04*W4 + dt05*W5 + dt06*W6 + dt07*W7 369 $ + dt08*W8 + dt09*W9 + dt10*W10+ dt11*W11) 370 371 Gamma2 =(delrho(n,1,1)*delrho(n,1,1) + 372 & delrho(n,2,1)*delrho(n,2,1) + 373 & delrho(n,3,1)*delrho(n,3,1))/F4 374 Gamma = dsqrt(Gamma2) 375 if(gamma.lt.dtol) goto 10 376 377 X = GAMMA/RHO43 378 S = X/(F48*PI*PI)**F1o3 379 s2 = s*s 380 Deno = (F1 + Mu*s2/kapa) 381 fx1=F1+kapa*(F1-F1/Deno) 382 fx2=F1+kapas*(F1-DExp(-Mus*s2/kapas)) 383 If(UseLC) then 384 CALL LRCLSDA(EMU,RHOO,ElSR,PDUM) 385 ElLR = Ax*Rho43-ElSR 386 else 387 ElSR = Ax*Rho43 388 ElLR = F0 389 endIf 390 GGA1 = ElSR*fx1 391 GGA2 = ElSR*fx2 392 GGA3 = ElLR*fx1 393 GGA4 = ElLR*fx2 394C 395 Ex = Ex +F2*(GGA1*Fsig1 + GGA2*Fsig2 396 $ + GGA3*Fsig3 + GGA4*Fsig4)*qwght(n) 397 if(ldew) func(n)=func(n)+F2*(GGA1*Fsig1+GGA2*Fsig2 398 $ + GGA3*Fsig3+GGA4*Fsig4) 399 400c 401c functional derivatives 402c 403 ds2dRho = -(F8/F3) * s2/rhoo 404 ds2dG = s2/Gamma2 405C 406 dfx1ds2 = Mu*(F1/(Deno*Deno)) 407 dfx1dRho = dfx1ds2*ds2dRho 408 dfx1dG = dfx1ds2*ds2dG 409C 410 dfx2ds2 = Mus*DExp(-Mus*s2/kapas) 411 dfx2dRho = dfx2ds2*ds2dRho 412 dfx2dG = dfx2ds2*ds2dG 413c 414 dF1dW = (at01 + F2*at02*W1 + F3*at03*W2 415 $ + F4*at04*W3 + F5*at05*W4 416 $ + F6*at06*W5 + F7*at07*W6 417 $ + F8*at08*W7 + F9*at09*W8 418 $ + F10*at10*W9+F11*at11*W10) 419 dF2dW = (bt01 + F2*bt02*W1 + F3*bt03*W2 420 $ + F4*bt04*W3 + F5*bt05*W4 421 $ + F6*bt06*W5 + F7*bt07*W6 422 $ + F8*bt08*W7 + F9*bt09*W8 423 $ + F10*Bt10*W9+F11*Bt11*W10) 424 dF3dW = (ct01 + F2*ct02*W1 + F3*ct03*W2 425 $ + F4*ct04*W3 + F5*ct05*W4 426 $ + F6*ct06*W5 + F7*ct07*W6 427 $ + F8*ct08*W7 + F9*ct09*W8 428 $ + F10*ct10*W9+F11*ct11*W10) 429 dF4dW = (dt01 + F2*dt02*W1 + F3*dt03*W2 430 $ + F4*dt04*W3 + F5*dt05*W4 431 $ + F6*dt06*W5 + F7*dt07*W6 432 $ + F8*dt08*W7 + F9*dt09*W8 433 $ + F10*dt10*W9+F11*dt11*W10) 434c 435 dWdT = F2/((F1 + Tsig)**F2) 436 dTdR = ((F6*PI*PI)**F2o3)*(rhoo**F2o3)/tauN 437 dTdTau = -TauUEG/tauN**F2 438C 439 If(UseLC) then 440 dElSRdR = PDUM 441 dElLRdR = Ax*F4o3*Rho13-PDUM 442 else 443 dElSRdR=Ax*F4o3*Rho13 444 dElLRdR=F0 445 endIf 446 dGGA1dR = dElSRdR*fx1 + ElSR*dfx1dRho 447 dGGA2dR = dElSRdR*fx2 + ElSR*dfx2dRho 448 dGGA3dR = dElLRdR*fx1 + ElLR*dfx1dRho 449 dGGA4dR = dElLRdR*fx2 + ElLR*dfx2dRho 450c 451 dF1dR = dF1dW*dWdT*dTdR 452 dF1dTau=dF1dW*dWdT*dTdTau 453 dF2dR = dF2dW*dWdT*dTdR 454 dF2dTau=dF2dW*dWdT*dTdTau 455 dF3dR = dF3dW*dWdT*dTdR 456 dF3dTau=dF3dW*dWdT*dTdTau 457 dF4dR = dF4dW*dWdT*dTdR 458 dF4dTau=dF4dW*dWdT*dTdTau 459c 460 dGGA1dG = ElSR*dfx1dG 461 dGGA2dG = ElSR*dfx2dG 462 dGGA3dG = ElLR*dfx1dG 463 dGGA4dG = ElLR*dfx2dG 464c 465 Amat(n,1) = Amat(n,1) +dGGA1dR*Fsig1 + GGA1*dF1dR 466 $ +dGGA2dR*Fsig2 + GGA2*dF2dR 467 $ +dGGA3dR*Fsig3 + GGA3*dF3dR 468 $ +dGGA4dR*Fsig4 + GGA4*dF4dR 469 Cmat(n,1)= Cmat(n,1) +dGGA1dG*Fsig1 + dGGA2dG*Fsig2 470 $ +dGGA3dG*Fsig3 + dGGA4dG*Fsig4 471 Mmat(n,1)= Mmat(n,1) +GGA1*dF1dTau + GGA2*dF2dTau 472 $ +GGA3*dF3dTau + GGA4*dF4dTau 473c 47410 continue 475c 476c UUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUnrestricted 477 else ! ipol=2 478c 479c ======> SPIN-UNRESTRICTED <====== 480 481c 482c use spin density functional theory ie n-->2n 483c Ex=(1/2)Ex[2*alpha] + (1/2)Ex[2*beta] 484c 485c Alpha ALPHA ALPHA 486c 487 do 20 n = 1, nq 488 if (rho(n,1).lt.DTol) goto 20 489 if (rho(n,2).lt.DTol) goto 25 490 rhoo = rho(n,2) 491 rho43 = rhoo**F4o3 492 rho13 = rho43/rhoo 493 rho53 = rhoo**F5o3 494c 495 tauN = tau(n,1)*F2 496 497 if(taun.lt.dtol) goto 25 498 tauu=tauN 499 TAUUEG=F3O5*((F6*PI*PI)**F2O3)*RHO53 500 Tsig =TauUEG/tauu 501 Wsig =(Tsig - F1)/(Tsig + F1) 502 W1=Wsig 503 W2=Wsig*W1 504 W3=Wsig*W2 505 W4=Wsig*W3 506 W5=Wsig*W4 507 W6=Wsig*W5 508 W7=Wsig*W6 509 W8=Wsig*W7 510 W9=Wsig*W8 511 W10=Wsig*W9 512 W11=Wsig*W10 513 Fsig1 =(at00 + at01*W1 + at02*W2 + at03*W3 514 $ + at04*W4 + at05*W5 + at06*W6 + at07*W7 515 $ + at08*W8 + at09*W9 + at10*W10+ at11*W11) 516 Fsig2 =(bt00 + bt01*W1 + bt02*W2 + bt03*W3 517 $ + bt04*W4 + bt05*W5 + bt06*W6 + bt07*W7 518 $ + bt08*W8 + bt09*W9 + bt10*W10+ bt11*W11) 519 Fsig3 =(ct00 + ct01*W1 + ct02*W2 + ct03*W3 520 $ + ct04*W4 + ct05*W5 + ct06*W6 + ct07*W7 521 $ + ct08*W8 + ct09*W9 + ct10*W10+ ct11*W11) 522 Fsig4 =(dt00 + dt01*W1 + dt02*W2 + dt03*W3 523 $ + dt04*W4 + dt05*W5 + dt06*W6 + dt07*W7 524 $ + dt08*W8 + dt09*W9 + dt10*W10+ dt11*W11) 525 526 Gamma2 =(delrho(n,1,1)*delrho(n,1,1) + 527 & delrho(n,2,1)*delrho(n,2,1) + 528 & delrho(n,3,1)*delrho(n,3,1)) 529 Gamma = dsqrt(Gamma2) 530 if(gamma.lt.dtol) goto 25 531 532 X = GAMMA/RHO43 533 S = X/(F48*PI*PI)**F1o3 534 s2 = s*s 535 Deno = (F1 + Mu*s2/kapa) 536 fx1=F1+kapa*(F1-F1/Deno) 537 fx2=F1+kapas*(F1-DExp(-Mus*s2/kapas)) 538 If(UseLC) then 539 CALL LRCLSDA(EMU,RHOO,ElSR,PDUM) 540 ElLR = Ax*Rho43-ElSR 541 else 542 ElSR = Ax*Rho43 543 ElLR = F0 544 endIf 545 GGA1 = ElSR*fx1 546 GGA2 = ElSR*fx2 547 GGA3 = ElLR*fx1 548 GGA4 = ElLR*fx2 549C 550 Ex = Ex + (GGA1*Fsig1 + GGA2*Fsig2 551 $ + GGA3*Fsig3 + GGA4*Fsig4)*qwght(n) 552 if(ldew) func(n)=func(n)+ (GGA1*Fsig1+GGA2*Fsig2 553 $ + GGA3*Fsig3+GGA4*Fsig4) 554c 555c functional derivatives 556c 557 ds2dRho = -(F8/F3) * s2/rhoo 558 ds2dG = s2/Gamma2 559C 560 dfx1ds2 = Mu*(F1/(Deno*Deno)) 561 dfx1dRho = dfx1ds2*ds2dRho 562 dfx1dG = dfx1ds2*ds2dG 563C 564 dfx2ds2 = Mus*DExp(-Mus*s2/kapas) 565 dfx2dRho = dfx2ds2*ds2dRho 566 dfx2dG = dfx2ds2*ds2dG 567c 568 dF1dW = (at01 + F2*at02*W1 + F3*at03*W2 569 $ + F4*at04*W3 + F5*at05*W4 570 $ + F6*at06*W5 + F7*at07*W6 571 $ + F8*at08*W7 + F9*at09*W8 572 $ + F10*at10*W9+F11*at11*W10) 573 dF2dW = (bt01 + F2*bt02*W1 + F3*bt03*W2 574 $ + F4*bt04*W3 + F5*bt05*W4 575 $ + F6*bt06*W5 + F7*bt07*W6 576 $ + F8*bt08*W7 + F9*bt09*W8 577 $ + F10*Bt10*W9+F11*Bt11*W10) 578 dF3dW = (ct01 + F2*ct02*W1 + F3*ct03*W2 579 $ + F4*ct04*W3 + F5*ct05*W4 580 $ + F6*ct06*W5 + F7*ct07*W6 581 $ + F8*ct08*W7 + F9*ct09*W8 582 $ + F10*ct10*W9+F11*ct11*W10) 583 dF4dW = (dt01 + F2*dt02*W1 + F3*dt03*W2 584 $ + F4*dt04*W3 + F5*dt05*W4 585 $ + F6*dt06*W5 + F7*dt07*W6 586 $ + F8*dt08*W7 + F9*dt09*W8 587 $ + F10*dt10*W9+F11*dt11*W10) 588 589 dWdT = F2/((F1 + Tsig)**F2) 590 dTdR = ((F6*PI*PI)**F2o3)*(rhoo**F2o3)/tauN 591 dTdTau = -TauUEG/tauN**F2 592C 593 If(UseLC) then 594 dElSRdR = PDUM 595 dElLRdR = Ax*F4o3*Rho13-PDUM 596 else 597 dElSRdR=Ax*F4o3*Rho13 598 dElLRdR=F0 599 endIf 600 dGGA1dR = dElSRdR*fx1 + ElSR*dfx1dRho 601 dGGA2dR = dElSRdR*fx2 + ElSR*dfx2dRho 602 dGGA3dR = dElLRdR*fx1 + ElLR*dfx1dRho 603 dGGA4dR = dElLRdR*fx2 + ElLR*dfx2dRho 604c 605 dF1dR = dF1dW*dWdT*dTdR 606 dF1dTau=dF1dW*dWdT*dTdTau 607 dF2dR = dF2dW*dWdT*dTdR 608 dF2dTau=dF2dW*dWdT*dTdTau 609 dF3dR = dF3dW*dWdT*dTdR 610 dF3dTau=dF3dW*dWdT*dTdTau 611 dF4dR = dF4dW*dWdT*dTdR 612 dF4dTau=dF4dW*dWdT*dTdTau 613c 614 dGGA1dG = ElSR*dfx1dG 615 dGGA2dG = ElSR*dfx2dG 616 dGGA3dG = ElLR*dfx1dG 617 dGGA4dG = ElLR*dfx2dG 618c 619 Amat(n,1) = Amat(n,1) +dGGA1dR*Fsig1 + GGA1*dF1dR 620 $ +dGGA2dR*Fsig2 + GGA2*dF2dR 621 $ +dGGA3dR*Fsig3 + GGA3*dF3dR 622 $ +dGGA4dR*Fsig4 + GGA4*dF4dR 623 Cmat(n,1)= Cmat(n,1) +dGGA1dG*Fsig1 + dGGA2dG*Fsig2 624 $ +dGGA3dG*Fsig3 + dGGA4dG*Fsig4 625 Mmat(n,1)= Mmat(n,1) +GGA1*dF1dTau + GGA2*dF2dTau 626 $ +GGA3*dF3dTau + GGA4*dF4dTau 627c 62825 continue 629c 630c Beta BETA BETA 631c 632 if (rho(n,3).lt.DTol) goto 20 633 rhoo = rho(n,3) 634 rho43 = rhoo**F4o3 635 rho13 = rho43/rhoo 636 rho53 = rhoo**F5o3 637c 638 639 tauN = tau(n,2)*F2 640 641 if(taun.lt.dtol) goto 20 642 tauu=tauN 643 TAUUEG=F3O5*((F6*PI*PI)**F2O3)*RHO53 644 Tsig =TauUEG/tauu 645 Wsig =(Tsig - F1)/(Tsig + F1) 646 W1=Wsig 647 W2=Wsig*W1 648 W3=Wsig*W2 649 W4=Wsig*W3 650 W5=Wsig*W4 651 W6=Wsig*W5 652 W7=Wsig*W6 653 W8=Wsig*W7 654 W9=Wsig*W8 655 W10=Wsig*W9 656 W11=Wsig*W10 657 Fsig1 =(at00 + at01*W1 + at02*W2 + at03*W3 658 $ + at04*W4 + at05*W5 + at06*W6 + at07*W7 659 $ + at08*W8 + at09*W9 + at10*W10+ at11*W11) 660 Fsig2 =(bt00 + bt01*W1 + bt02*W2 + bt03*W3 661 $ + bt04*W4 + bt05*W5 + bt06*W6 + bt07*W7 662 $ + bt08*W8 + bt09*W9 + bt10*W10+ bt11*W11) 663 Fsig3 =(ct00 + ct01*W1 + ct02*W2 + ct03*W3 664 $ + ct04*W4 + ct05*W5 + ct06*W6 + ct07*W7 665 $ + ct08*W8 + ct09*W9 + ct10*W10+ ct11*W11) 666 Fsig4 =(dt00 + dt01*W1 + dt02*W2 + dt03*W3 667 $ + dt04*W4 + dt05*W5 + dt06*W6 + dt07*W7 668 $ + dt08*W8 + dt09*W9 + dt10*W10+ dt11*W11) 669 670 Gamma2 =(delrho(n,1,2)*delrho(n,1,2) + 671 & delrho(n,2,2)*delrho(n,2,2) + 672 & delrho(n,3,2)*delrho(n,3,2)) 673 Gamma = dsqrt(Gamma2) 674 if(gamma.lt.dtol) goto 20 675 676 X = GAMMA/RHO43 677 S = X/(F48*PI*PI)**F1o3 678 s2 = s*s 679 Deno = (F1 + Mu*s2/kapa) 680 fx1=F1+kapa*(F1-F1/Deno) 681 fx2=F1+kapas*(F1-DExp(-Mus*s2/kapas)) 682 If(UseLC) then 683 CALL LRCLSDA(EMU,RHOO,ElSR,PDUM) 684 ElLR = Ax*Rho43-ElSR 685 else 686 ElSR = Ax*Rho43 687 ElLR = F0 688 endIf 689 GGA1 = ElSR*fx1 690 GGA2 = ElSR*fx2 691 GGA3 = ElLR*fx1 692 GGA4 = ElLR*fx2 693C 694 Ex = Ex + (GGA1*Fsig1 + GGA2*Fsig2 695 $ + GGA3*Fsig3 + GGA4*Fsig4)*qwght(n) 696 if(ldew) func(n)=func(n)+ (GGA1*Fsig1+GGA2*Fsig2 697 $ + GGA3*Fsig3+GGA4*Fsig4) 698 699c 700c functional derivatives 701c 702 ds2dRho = -(F8/F3) * s2/rhoo 703 ds2dG = s2/Gamma2 704C 705 dfx1ds2 = Mu*(F1/(Deno*Deno)) 706 dfx1dRho = dfx1ds2*ds2dRho 707 dfx1dG = dfx1ds2*ds2dG 708C 709 dfx2ds2 = Mus*DExp(-Mus*s2/kapas) 710 dfx2dRho = dfx2ds2*ds2dRho 711 dfx2dG = dfx2ds2*ds2dG 712c 713 dF1dW = (at01 + F2*at02*W1 + F3*at03*W2 714 $ + F4*at04*W3 + F5*at05*W4 715 $ + F6*at06*W5 + F7*at07*W6 716 $ + F8*at08*W7 + F9*at09*W8 717 $ + F10*at10*W9+F11*at11*W10) 718 dF2dW = (bt01 + F2*bt02*W1 + F3*bt03*W2 719 $ + F4*bt04*W3 + F5*bt05*W4 720 $ + F6*bt06*W5 + F7*bt07*W6 721 $ + F8*bt08*W7 + F9*bt09*W8 722 $ + F10*Bt10*W9+F11*Bt11*W10) 723 dF3dW = (ct01 + F2*ct02*W1 + F3*ct03*W2 724 $ + F4*ct04*W3 + F5*ct05*W4 725 $ + F6*ct06*W5 + F7*ct07*W6 726 $ + F8*ct08*W7 + F9*ct09*W8 727 $ + F10*ct10*W9+F11*ct11*W10) 728 dF4dW = (dt01 + F2*dt02*W1 + F3*dt03*W2 729 $ + F4*dt04*W3 + F5*dt05*W4 730 $ + F6*dt06*W5 + F7*dt07*W6 731 $ + F8*dt08*W7 + F9*dt09*W8 732 $ + F10*dt10*W9+F11*dt11*W10) 733 734 dWdT = F2/((F1 + Tsig)**F2) 735 dTdR = ((F6*PI*PI)**F2o3)*(rhoo**F2o3)/tauN 736 dTdTau = -TauUEG/tauN**F2 737C 738 If(UseLC) then 739 dElSRdR = PDUM 740 dElLRdR = Ax*F4o3*Rho13-PDUM 741 else 742 dElSRdR=Ax*F4o3*Rho13 743 dElLRdR=F0 744 endIf 745 dGGA1dR = dElSRdR*fx1 + ElSR*dfx1dRho 746 dGGA2dR = dElSRdR*fx2 + ElSR*dfx2dRho 747 dGGA3dR = dElLRdR*fx1 + ElLR*dfx1dRho 748 dGGA4dR = dElLRdR*fx2 + ElLR*dfx2dRho 749c 750 dF1dR = dF1dW*dWdT*dTdR 751 dF1dTau=dF1dW*dWdT*dTdTau 752 dF2dR = dF2dW*dWdT*dTdR 753 dF2dTau=dF2dW*dWdT*dTdTau 754 dF3dR = dF3dW*dWdT*dTdR 755 dF3dTau=dF3dW*dWdT*dTdTau 756 dF4dR = dF4dW*dWdT*dTdR 757 dF4dTau=dF4dW*dWdT*dTdTau 758c 759 dGGA1dG = ElSR*dfx1dG 760 dGGA2dG = ElSR*dfx2dG 761 dGGA3dG = ElLR*dfx1dG 762 dGGA4dG = ElLR*dfx2dG 763c 764 Amat(n,2) = Amat(n,2) +dGGA1dR*Fsig1 + GGA1*dF1dR 765 $ +dGGA2dR*Fsig2 + GGA2*dF2dR 766 $ +dGGA3dR*Fsig3 + GGA3*dF3dR 767 $ +dGGA4dR*Fsig4 + GGA4*dF4dR 768 Cmat(n,3)= Cmat(n,3) +dGGA1dG*Fsig1 + dGGA2dG*Fsig2 769 $ +dGGA3dG*Fsig3 + dGGA4dG*Fsig4 770 Mmat(n,2)= Mmat(n,2) +GGA1*dF1dTau + GGA2*dF2dTau 771 $ +GGA3*dF3dTau + GGA4*dF4dTau 772c 77320 continue 774 endif 775 return 776 end 777c 778 Subroutine xc_xm11_d2() 779 call errquit(' not coded ',0,0) 780 return 781 end 782c 783 SUBROUTINE LRCLSDA(Emu,Rho,F,D1F) 784c 785c*********************************************** 786c 787c INPUT: 788c Emu - Value of mu (or omega) 789c Rho - Spin density 790c 791c OUTPUT: 792c F - Functional value 793c D1F - First derivative 794c 795c*********************************************** 796c 797 IMPLICIT REAL*8 (a-h,o-z) 798 Save F1, F2, F3, F4, F5, F6, F7, F8, F9 799 DATA F1/1.0D+00/,F2/2.0D+00/,F3/3.0D+00/,F4/4.0D+00/,F5/5.0D+00/, 800 $ F6/6.0D+00/,F7/7.0D+00/,F8/8.0D+00/,F9/9.0D+00/ 801C 802 PARAMETER( PI = 3.1415926535897932384626433832795D+00 ) 803C 804 F1o2 = F1 / F2 805 F1o3 = F1 / F3 806 F1o4 = F1 / F4 807 F4o3 = F4 / F3 808 F8o3 = F8 / F3 809 PI12 = SQRT(Pi) 810C 811 AX = -(F3/F2) * (F4o3*PI)**(-F1o3) 812 Cmu = (F6*Pi**F2)**F1o3 813C 814 Rho13 = Rho**F1o3 815 Rho43 = Rho**F4o3 816c 817 tmu = Emu/(F2*Cmu*Rho13) 818 tmu2 = tmu*tmu 819 tmu3 = tmu*tmu2 820c 821 W = DExp(-F1o4/tmu2) 822 ERFV = DErf( F1o2/tmu) 823 dtmudR = -F1o3*tmu / Rho 824c 825 Fsr = F1-F4o3*tmu*(-F6*tmu+F8*tmu3+W* 826 $ (F4*tmu-F8*tmu3)+F2*PI12*ERFV) 827 dFsrdtmu = F8o3*(F2*tmu*(F3-F8*tmu2+W* 828 $ (-F1+F8*tmu2))-PI12*ERFV) 829c 830 F = Ax*Rho43*Fsr 831 D1F = Ax*F4o3*Rho13*Fsr + Ax*Rho43*(dFsrdtmu*dtmudR) 832c 833 RETURN 834 END 835 836