1c M05 or M05-2X 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 M05 9c ijzy - 2 M05-2X 10c References: 11c [a] Zhao, Y.; Schultz, N. E.; Truhlar, D. G. J. Chem. Phys. 2005, 123, 161103; 12c Note that in this communication we interchanged cCab,i and cCss,i in Table 1. 13c [b] Zhao, Y.; Schultz, N. E.; Truhlar, D. G. J. Chem. Theory Comput. 2006, in press. 14 15 16 Subroutine xc_xm05(tol_rho, fac,lfac,nlfac, rho, delrho, 17 & Amat, Cmat, nq, ipol, Ex, 18 & qwght, ldew, func, tau, Mmat, ijzy) 19 20c 21c$Id$ 22c 23 implicit none 24c 25#include "errquit.fh" 26c 27 double precision fac, Ex 28 integer nq, ipol 29 logical lfac, nlfac,ldew 30 double precision func(*) ! value of the functional [output] 31c 32c Charge Density & Its Cube Root 33c 34 double precision rho(nq,ipol*(ipol+1)/2) 35c 36c Charge Density Gradient 37c 38 double precision delrho(nq,3,ipol) 39c 40c Quadrature Weights 41c 42 double precision qwght(nq) 43c 44c Sampling Matrices for the XC Potential & Energy 45c 46 double precision Amat(nq,ipol), Cmat(nq,*) 47 48 49 double precision tol_rho, pi 50 51c 52 integer n, ijzy 53 double precision at1, at2, at3, at4, at5, at6, at7, at8, at9 54 double precision at, at10, at11, C1, C2, fL, fNL 55 double precision rrho, rho43, rho13, rhoo, rho53 56 double precision Gamma2, Gamma 57 double precision TauUEG, Tsig, Wsig, W1, W2, W3, W4, W5, W6 58 double precision W7, W8, W9, W10, W11, Fsig 59c 60c kinetic energy density or tau 61c 62 double precision tau(nq,ipol), Mmat(nq,*) 63 64 double precision tauN,tauu,DTol 65 double precision F83, F23, F53, F43, F13, F1o3 66 double precision F1o4, F2o3, F3o2, F4o3, F4o9, F3o5 67 double precision One, Two, Three, Four, Five, Six, Seven, Eight 68 double precision Nine, F10, F11 69 double precision Cs, Ax, P32, s, s2, En, Ed, E, dE, dEn, dEd 70 71c functional derivatives below FFFFFFFFFFFF 72 73 double precision dFdW, dWdT, dTdR, dTdTau, dGGAdR, dFdR 74 double precision dFdTau, dGGAdG,tauW 75 76c functional derivatives above FFFFFFFFFFFF 77 78 79 parameter( pi = 3.1415926535897932384626433832795d0 ) 80 81 82 parameter (F1o3=1.d0/3.d0, F1o4=1.d0/4.d0, F2o3=2.d0/3.d0, 83 & F3o2=3.d0/2.d0,F13=1.d0/3.d0) 84 parameter (F4o3=4.d0/3.d0, F4o9=4.d0/9.d0, F3o5=3.d0/5.d0) 85 parameter (F83=8.d0/3.0d0, F23=2.0d0/3.d0, F53=5.d0/3.d0) 86 parameter (One=1.0d0, Two=2.0d0, Three=3.0d0, Four=4.0d0, 87 & Five=5.0d0,Six=6.0d0, Seven=7.0d0, 88 & Eight=8.0d0, Nine=9.0d0,F10=10.d0, F11=11.d0) 89 90 91 if (ijzy.eq.1) then 92c Parameters for M05 93 at1= 0.08151d0 94 at2= -0.43956d0 95 at3= -3.22422d0 96 at4= 2.01819d0 97 at5= 8.79431d0 98 at6= -0.00295d0 99 at7= 9.82029d0 100 at8= -4.82351d0 101 at9= -48.17574d0 102 at10= 3.64802d0 103 at11= 34.02248d0 104 elseif (ijzy.eq.2) then 105c Parameters for M05-2X 106 at1= -0.56833d0 107 at2= -1.30057d0 108 at3= 5.50070d0 109 at4= 9.06402d0 110 at5= -32.21075d0 111 at6= -23.73298d0 112 at7= 70.22996d0 113 at8= 29.88614d0 114 at9= -60.25778d0 115 at10= -13.22205d0 116 at11= 15.23694d0 117 else 118 call errquit("xc_xm05: illegal ijzy",ijzy,UERR) 119 endif 120 121 at=1.0d0 122 C1 = 0.2195149727645171d0 123 C2 = C1/0.804d0 124 DTol=tol_rho 125C 126C Scale factors for local and non-local contributions. 127C 128 fL = fac 129 fNL = fac 130 Cs = 0.5d0/(3.0d0*pi*pi)**F13 131 P32 = (3.d0*pi**2)**F23 132 133c 134 Ax = (-0.75d0)*(3.0d0/pi)**F13 135 136 137c 138 if (ipol.eq.1 )then 139c 140c ======> SPIN-RESTRICTED <====== 141c or 142c SPIN-UNPOLARIZED 143c 144c 145 do 10 n = 1, nq 146 147 if (rho(n,1).lt.DTol) goto 10 148 rhoo = rho(n,1) 149 rho43 = rhoo**F4o3 150 rrho = 1d0/rhoo ! reciprocal of rho 151 rho13 = rho43*rrho 152 rho53 = rhoo**F53 153 154c 155 156 tauN = tau(n,1) 157 tauu=tauN 158 Gamma2 = delrho(n,1,1)*delrho(n,1,1) + 159 & delrho(n,2,1)*delrho(n,2,1) + 160 & delrho(n,3,1)*delrho(n,3,1) 161 Gamma = dsqrt(Gamma2) 162 if (gamma.lt.DTol) goto 10 163 TauUEG=0.3d0*P32*rho53 164 Tsig =TauUEG/tauu 165 Wsig =(Tsig-One)/(Tsig+One) 166 W1=Wsig 167 W2=Wsig*W1 168 W3=Wsig*W2 169 W4=Wsig*W3 170 W5=Wsig*W4 171 W6=Wsig*W5 172 W7=Wsig*W6 173 W8=Wsig*W7 174 W9=Wsig*W8 175 W10=Wsig*W9 176 W11=Wsig*W10 177 Fsig =at*(at1*W1+ at2*W2 + at3*W3 178 & + at4*W4 + at5*W5 + at6*W6 + at7*W7 179 & + at8*W8 + at9*W9 + at10*W10 + at11*W11) 180 181 s = Cs*Gamma/rho43 182 s2 = s*s 183 En = C1*s2 184 Ed = One + C2*s2 185 E = En/Ed 186 Ex = Ex + rho43*Ax*(fL+fNL*E)*(One+Fsig)*qwght(n) 187 if(ldew) func(n)=func(n)+rho43*Ax*(fL+fNL*E)*(One+Fsig) 188c 189c functional derivatives 190c 191 dEn = Two*C1*s 192 dEd = Two*C2*s 193 dE = (dEn*Ed-En*dEd)/(Ed*Ed) 194 dFdW = at*(at1 + Two*at2*W1 + Three*at3*W2 195 & + Four*at4*W3 + Five*at5*W4 196 & + Six*at6*W5 + Seven*at7*W6 197 & + Eight*at8*W7 + Nine*at9*W8 198 & + F10*at10*W9 + F11*at11*W10) 199 dWdT = Two/((One + Tsig)**2) 200 dTdR = (0.5d0*P32*rho13*rho13)/tauu 201 dTdTau = -TauUEG/tauu**2 202 dGGAdR = F4o3*rho13*Ax*(fL+fNL*(E-s*dE)) 203 dFdR = dFdW*dWdT*dTdR 204 dFdTau=dFdW*dWdT*dTdTau 205 dGGAdG =(fNL*dE*s/(Two*Gamma2)) 206 Amat(n,1) = Amat(n,1) + dGGAdR*(One+Fsig) 207 & + (fL+fNL*E)*Ax*rho43*dFdR 208 Cmat(n,1)= Cmat(n,1) + 209 & Two*dGGAdG*Ax*rho43*(One+Fsig) 210 Mmat(n,1)= Mmat(n,1) + 0.5d0*rho43*Ax*(fL+fNL*E)*dFdTau 211 212 213 21410 continue 215 216 217c UUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUnrestricted 218 else ! ipol=2 219c 220c ======> SPIN-UNRESTRICTED <====== 221 222c 223c use spin density functional theory ie n-->2n 224c Ex=(1/2)Ex[2*alpha] + (1/2)Ex[2*beta] 225 226 do 20 n = 1, nq 227 if (rho(n,1).lt.DTol) goto 20 228c 229c Alpha ALPHA ALPHA 230c 231 if (rho(n,2).lt.DTol) goto 25 232 rhoo = Two*rho(n,2) 233 rho43 = rhoo**F4o3 234 rrho = 1.0d0/rhoo ! reciprocal of rho 235 rho13 = rho43*rrho 236 rho53 = rhoo**F53 237 238c 239 240 tauN = tau(n,1) 241 tauu = Two*tauN 242 TauUEG=0.3d0*P32*rho53 243 Tsig =TauUEG/tauu 244 Wsig =(Tsig-One)/(Tsig+One) 245 W1=Wsig 246 W2=Wsig*W1 247 W3=Wsig*W2 248 W4=Wsig*W3 249 W5=Wsig*W4 250 W6=Wsig*W5 251 W7=Wsig*W6 252 W8=Wsig*W7 253 W9=Wsig*W8 254 W10=Wsig*W9 255 W11=Wsig*W10 256 Fsig =at*(at1*W1+ at2*W2 + at3*W3 257 & + at4*W4 + at5*W5 + at6*W6 + at7*W7 258 & + at8*W8 + at9*W9 + at10*W10 + at11*W11) 259 260 261 Gamma2 = delrho(n,1,1)*delrho(n,1,1) + 262 & delrho(n,2,1)*delrho(n,2,1) + 263 & delrho(n,3,1)*delrho(n,3,1) 264 Gamma2 = Four*Gamma2 265 Gamma = dsqrt(Gamma2) 266 if (gamma.lt.DTol) goto 25 267 268 s = Cs*Gamma/rho43 269 s2 = s*s 270 En = C1*s2 271 Ed = One + C2*s2 272 E = En/Ed 273 Ex = Ex + rho43*Ax*(fL+fNL*E)*(One+Fsig)*0.5d0*qwght(n) 274 if(ldew) func(n)= 275 = func(n)+rho43*Ax*(fL+fNL*E)*(One+Fsig)*.5d0 276c 277c functional derivatives 278c 279 dEn = Two*C1*s 280 dEd = Two*C2*s 281 dE = (dEn*Ed-En*dEd)/(Ed*Ed) 282 dFdW = at*(at1 + Two*at2*W1 + Three*at3*W2 283 & + Four*at4*W3 + Five*at5*W4 284 & + Six*at6*W5 + Seven*at7*W6 285 & + Eight*at8*W7 + Nine*at9*W8 286 & + F10*at10*W9 + F11*at11*W10) 287 dWdT = Two/((One + Tsig)**2) 288 dTdR = Two*(0.5d0*P32*rho13*rho13)/tauu 289 dTdTau = -Two*TauUEG/tauu**2 290 dGGAdR = Two*F4o3*rho13*Ax*(fL+fNL*(E-s*dE)) 291 dFdR = dFdW*dWdT*dTdR 292 dFdTau=dFdW*dWdT*dTdTau 293 dGGAdG =Four*(fNL*dE*s/(Two*Gamma2)) 294 295 Amat(n,1) = Amat(n,1) + (dGGAdR*(One+Fsig) 296 & + (fL+fNL*E)*Ax*rho43*dFdR)*0.5d0 297 Cmat(n,1)= Cmat(n,1) + 298 & dGGAdG*Ax*rho43*(One+Fsig)*0.5d0 299 Mmat(n,1)= Mmat(n,1) + 300 & 0.5d0*0.5d0*rho43*Ax*(fL+fNL*E)*dFdTau 301 302c write (*,*) "Ex,Amat(n,1),Cmat(n,1)", 303c & Ex,Amat(n,1),Cmat(n,1) 304 305c 306c Beta BETA BETA 307c 308 30925 continue 310 311c 312c Beta 313c 314 if (rho(n,3).lt.DTol) goto 20 315 rhoo = Two*rho(n,3) 316 rho43 = rhoo**F4o3 317 rrho = 1.0d0/rhoo ! reciprocal of rho 318 rho13 = rho43*rrho 319 rho53 = rhoo**F53 320 321c 322 323 tauN = tau(n,2) 324 tauu = Two*tauN 325 TauUEG=0.3d0*P32*rho53 326 Tsig =TauUEG/tauu 327 Wsig =(Tsig-One)/(Tsig+One) 328 W1=Wsig 329 W2=Wsig*W1 330 W3=Wsig*W2 331 W4=Wsig*W3 332 W5=Wsig*W4 333 W6=Wsig*W5 334 W7=Wsig*W6 335 W8=Wsig*W7 336 W9=Wsig*W8 337 W10=Wsig*W9 338 W11=Wsig*W10 339 Fsig =at*(at1*W1+ at2*W2 + at3*W3 340 & + at4*W4 + at5*W5 + at6*W6 + at7*W7 341 & + at8*W8 + at9*W9 + at10*W10 + at11*W11) 342 343 344 Gamma2 = delrho(n,1,2)*delrho(n,1,2) + 345 & delrho(n,2,2)*delrho(n,2,2) + 346 & delrho(n,3,2)*delrho(n,3,2) 347 Gamma2 = Four*Gamma2 348 Gamma = dsqrt(Gamma2) 349 if (gamma.lt.DTol) goto 20 350 s = Cs*Gamma/rho43 351 s2 = s*s 352 En = C1*s2 353 Ed = One + C2*s2 354 E = En/Ed 355 Ex = Ex + rho43*Ax*(fL+fNL*E)*(One+Fsig)*0.5d0*qwght(n) 356 if(ldew) func(n)= 357 = func(n)+rho43*Ax*(fL+fNL*E)*(One+Fsig)*.5d0 358c 359c functional derivatives 360c 361 dEn = Two*C1*s 362 dEd = Two*C2*s 363 dE = (dEn*Ed-En*dEd)/(Ed*Ed) 364 dFdW = at*(at1 + Two*at2*W1 + Three*at3*W2 365 & + Four*at4*W3 + Five*at5*W4 366 & + Six*at6*W5 + Seven*at7*W6 367 & + Eight*at8*W7 + Nine*at9*W8 368 & + F10*at10*W9 + F11*at11*W10) 369 dWdT = Two/((One + Tsig)**2) 370 dTdR = Two*(0.5d0*P32*rho13*rho13)/tauu 371 dTdTau = -Two*TauUEG/tauu**2 372 dGGAdR = Two*F4o3*rho13*Ax*(fL+fNL*(E-s*dE)) 373 dFdR = dFdW*dWdT*dTdR 374 dFdTau=dFdW*dWdT*dTdTau 375 dGGAdG =Four*(fNL*dE*s/(Two*Gamma2)) 376 377 Amat(n,2) = Amat(n,2) + (dGGAdR*(One+Fsig) 378 & + (fL+fNL*E)*Ax*rho43*dFdR)*0.5d0 379 Cmat(n,3)= Cmat(n,3) + 380 & dGGAdG*Ax*rho43*(One+Fsig)*0.5d0 381 Mmat(n,2)= Mmat(n,2) + 382 & 0.5d0*0.5d0*rho43*Ax*(fL+fNL*E)*dFdTau 383 384 385c 38620 continue 387 endif 388c 389 return 390 end 391 392 393 394 395 Subroutine xc_xm05_d2() 396 call errquit(' xm05: d2 not coded ',0,0) 397 return 398 end 399 400 401c---------------------------------------------------------------------- 402c dlDF exchange functional 403c META GGA 404C utilizes ingredients: 405c rho - density 406c delrho - gradient of density 407c tau - K.S kinetic energy density 408c tauU - uniform-gas KE density 409c 410c References: 411c [a] Pernal,Podeszwa,Patkowski,Szalewicz, PRL 103 263201 (2009) 412 413 Subroutine xc_xdldf(tol_rho, fac,lfac,nlfac, rho, delrho, 414 & Amat, Cmat, nq, ipol, Ex, 415 & qwght, ldew, func, tau, Mmat) 416c 417c 418 implicit none 419c 420c 421 double precision fac, Ex 422 integer nq, ipol 423 logical lfac, nlfac,ldew 424 double precision func(*) ! value of the functional [output] 425c 426c Charge Density & Its Cube Root 427c 428 double precision rho(nq,ipol*(ipol+1)/2) 429c 430c Charge Density Gradient 431c 432 double precision delrho(nq,3,ipol) 433c 434c Quadrature Weights 435c 436 double precision qwght(nq) 437c 438c Sampling Matrices for the XC Potential & Energy 439c 440 double precision Amat(nq,ipol), Cmat(nq,*) 441 442 443 double precision tol_rho, pi 444 445c 446 integer n 447 double precision at1, at2, at3, at4, at5, at6, at7, at8, at9 448 double precision at, at10, at11, C1, C2, fL, fNL 449 double precision rrho, rho43, rho13, rhoo, rho53 450 double precision Gamma2, Gamma 451 double precision TauUEG, Tsig, Wsig, W1, W2, W3, W4, W5, W6 452 double precision W7, W8, W9, W10, W11, Fsig 453c 454c kinetic energy density or tau 455c 456 double precision tau(nq,ipol), Mmat(nq,*) 457 458 double precision tauN,tauu,DTol 459 double precision F83, F23, F53, F43, F13, F1o3 460 double precision F1o4, F2o3, F3o2, F4o3, F4o9, F3o5 461 double precision One, Two, Three, Four, Five, Six, Seven, Eight 462 double precision Nine, F10, F11 463 double precision Cs, Ax, P32, s, s2, En, Ed, E, dE, dEn, dEd 464 465c functional derivatives below FFFFFFFFFFFF 466 467 double precision dFdW, dWdT, dTdR, dTdTau, dGGAdR, dFdR 468 double precision dFdTau, dGGAdG,tauW 469 470c functional derivatives above FFFFFFFFFFFF 471 472 473 parameter( pi = 3.1415926535897932384626433832795d0 ) 474 475 476 parameter (F1o3=1.d0/3.d0, F1o4=1.d0/4.d0, F2o3=2.d0/3.d0, 477 & F3o2=3.d0/2.d0,F13=1.d0/3.d0) 478 parameter (F4o3=4.d0/3.d0, F4o9=4.d0/9.d0, F3o5=3.d0/5.d0) 479 parameter (F83=8.d0/3.0d0, F23=2.0d0/3.d0, F53=5.d0/3.d0) 480 parameter (One=1.0d0, Two=2.0d0, Three=3.0d0, Four=4.0d0, 481 & Five=5.0d0,Six=6.0d0, Seven=7.0d0, 482 & Eight=8.0d0, Nine=9.0d0,F10=10.d0, F11=11.d0) 483 484c Parameters for dlDF 485 at1= -1.637571d-01 486 at2= -1.880028d-01 487 at3= -4.490609d-01 488 at4= -8.2359d-03 489 at5= 0.0d0 490 at6= 0.0d0 491 at7= 0.0d0 492 at8= 0.0d0 493 at9= 0.0d0 494 at10= 0.0d0 495 at11= 0.0d0 496 497 498 at=1.0d0 499 C1 = 0.3511128d0 500 C2 = C1/4.8827323d0 501 DTol=tol_rho 502C 503C Scale factors for local and non-local contributions. 504C 505 fL = fac 506 fNL = fac 507 Cs = 0.5d0/(3.0d0*pi*pi)**F13 508 P32 = (3.d0*pi**2)**F23 509 510c 511 Ax = (-0.75d0)*(3.0d0/pi)**F13 512 513 514c 515 if (ipol.eq.1 )then 516c 517c ======> SPIN-RESTRICTED <====== 518c or 519c SPIN-UNPOLARIZED 520c 521c 522 do 10 n = 1, nq 523 524 if (rho(n,1).lt.DTol) goto 10 525 rhoo = rho(n,1) 526 rho43 = rhoo**F4o3 527 rrho = 1d0/rhoo ! reciprocal of rho 528 rho13 = rho43*rrho 529 rho53 = rhoo**F53 530 531c 532 533 tauN = tau(n,1) 534 tauu=tauN 535cedo if (taun.lt.sqrt(DTol)) goto 10 536 Gamma2 = delrho(n,1,1)*delrho(n,1,1) + 537 & delrho(n,2,1)*delrho(n,2,1) + 538 & delrho(n,3,1)*delrho(n,3,1) 539 Gamma = dsqrt(Gamma2) 540 if (gamma.lt.DTol) goto 10 541 TauUEG=0.3d0*P32*rho53 542 Tsig =TauUEG/tauu 543 Wsig =(Tsig-One)/(Tsig+One) 544 W1=Wsig 545 W2=Wsig*W1 546 W3=Wsig*W2 547 W4=Wsig*W3 548 W5=Wsig*W4 549 W6=Wsig*W5 550 W7=Wsig*W6 551 W8=Wsig*W7 552 W9=Wsig*W8 553 W10=Wsig*W9 554 W11=Wsig*W10 555 Fsig =at*(at1*W1+ at2*W2 + at3*W3 556 & + at4*W4 + at5*W5 + at6*W6 + at7*W7 557 & + at8*W8 + at9*W9 + at10*W10 + at11*W11) 558 559 s = Cs*Gamma/rho43 560 s2 = s*s 561 En = C1*s2 562 Ed = One + C2*s2 563 E = En/Ed 564 Ex = Ex + rho43*Ax*(fL+fNL*E)*(One+Fsig)*qwght(n) 565 if(ldew) func(n)=func(n)+rho43*Ax*(fL+fNL*E)*(One+Fsig) 566c 567c functional derivatives 568c 569 dEn = Two*C1*s 570 dEd = Two*C2*s 571 dE = (dEn*Ed-En*dEd)/(Ed*Ed) 572 dFdW = at*(at1 + Two*at2*W1 + Three*at3*W2 573 & + Four*at4*W3 + Five*at5*W4 574 & + Six*at6*W5 + Seven*at7*W6 575 & + Eight*at8*W7 + Nine*at9*W8 576 & + F10*at10*W9 + F11*at11*W10) 577 dWdT = Two/((One + Tsig)**2) 578 dTdR = (0.5d0*P32*rho13*rho13)/tauu 579 dTdTau = -TauUEG/tauu**2 580 dGGAdR = F4o3*rho13*Ax*(fL+fNL*(E-s*dE)) 581 dFdR = dFdW*dWdT*dTdR 582 dFdTau=dFdW*dWdT*dTdTau 583 dGGAdG =(fNL*dE*s/(Two*Gamma2)) 584 Amat(n,1) = Amat(n,1) + dGGAdR*(One+Fsig) 585 & + (fL+fNL*E)*Ax*rho43*dFdR 586 Cmat(n,1)= Cmat(n,1) + 587 & Two*dGGAdG*Ax*rho43*(One+Fsig) 588 Mmat(n,1)= Mmat(n,1) + 0.5d0*rho43*Ax*(fL+fNL*E)*dFdTau 589 590 591 59210 continue 593 594 595c UUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUnrestricted 596 else ! ipol=2 597c 598c ======> SPIN-UNRESTRICTED <====== 599 600c 601c use spin density functional theory ie n-->2n 602c Ex=(1/2)Ex[2*alpha] + (1/2)Ex[2*beta] 603 604 do 20 n = 1, nq 605 if (rho(n,1).lt.DTol) goto 20 606c 607c Alpha ALPHA ALPHA 608c 609 if (rho(n,2).lt.DTol) goto 25 610 rhoo = Two*rho(n,2) 611 rho43 = rhoo**F4o3 612 rrho = 1.0d0/rhoo ! reciprocal of rho 613 rho13 = rho43*rrho 614 rho53 = rhoo**F53 615 616c 617 618 tauN = tau(n,1) 619 tauu = Two*tauN 620 TauUEG=0.3d0*P32*rho53 621 Tsig =TauUEG/tauu 622 Wsig =(Tsig-One)/(Tsig+One) 623 W1=Wsig 624 W2=Wsig*W1 625 W3=Wsig*W2 626 W4=Wsig*W3 627 W5=Wsig*W4 628 W6=Wsig*W5 629 W7=Wsig*W6 630 W8=Wsig*W7 631 W9=Wsig*W8 632 W10=Wsig*W9 633 W11=Wsig*W10 634 Fsig =at*(at1*W1+ at2*W2 + at3*W3 635 & + at4*W4 + at5*W5 + at6*W6 + at7*W7 636 & + at8*W8 + at9*W9 + at10*W10 + at11*W11) 637 638 639 Gamma2 = delrho(n,1,1)*delrho(n,1,1) + 640 & delrho(n,2,1)*delrho(n,2,1) + 641 & delrho(n,3,1)*delrho(n,3,1) 642 Gamma2 = Four*Gamma2 643 Gamma = dsqrt(Gamma2) 644 if (gamma.lt.DTol) goto 25 645 646 s = Cs*Gamma/rho43 647 s2 = s*s 648 En = C1*s2 649 Ed = One + C2*s2 650 E = En/Ed 651 Ex = Ex + rho43*Ax*(fL+fNL*E)*(One+Fsig)*0.5d0*qwght(n) 652 if(ldew) func(n)= 653 = func(n)+rho43*Ax*(fL+fNL*E)*(One+Fsig)*.5d0 654c 655c functional derivatives 656c 657 dEn = Two*C1*s 658 dEd = Two*C2*s 659 dE = (dEn*Ed-En*dEd)/(Ed*Ed) 660 dFdW = at*(at1 + Two*at2*W1 + Three*at3*W2 661 & + Four*at4*W3 + Five*at5*W4 662 & + Six*at6*W5 + Seven*at7*W6 663 & + Eight*at8*W7 + Nine*at9*W8 664 & + F10*at10*W9 + F11*at11*W10) 665 dWdT = Two/((One + Tsig)**2) 666 dTdR = Two*(0.5d0*P32*rho13*rho13)/tauu 667 dTdTau = -Two*TauUEG/tauu**2 668 dGGAdR = Two*F4o3*rho13*Ax*(fL+fNL*(E-s*dE)) 669 dFdR = dFdW*dWdT*dTdR 670 dFdTau=dFdW*dWdT*dTdTau 671 dGGAdG =Four*(fNL*dE*s/(Two*Gamma2)) 672 673 Amat(n,1) = Amat(n,1) + (dGGAdR*(One+Fsig) 674 & + (fL+fNL*E)*Ax*rho43*dFdR)*0.5d0 675 Cmat(n,1)= Cmat(n,1) + 676 & dGGAdG*Ax*rho43*(One+Fsig)*0.5d0 677 Mmat(n,1)= Mmat(n,1) + 678 & 0.5d0*0.5d0*rho43*Ax*(fL+fNL*E)*dFdTau 679 680c write (*,*) "Ex,Amat(n,1),Cmat(n,1)", 681c & Ex,Amat(n,1),Cmat(n,1) 682 683c 684c Beta BETA BETA 685c 686 68725 continue 688 689c 690c Beta 691c 692 if (rho(n,3).lt.DTol) goto 20 693 rhoo = Two*rho(n,3) 694 rho43 = rhoo**F4o3 695 rrho = 1.0d0/rhoo ! reciprocal of rho 696 rho13 = rho43*rrho 697 rho53 = rhoo**F53 698 699c 700 701 tauN = tau(n,2) 702 tauu = Two*tauN 703 TauUEG=0.3d0*P32*rho53 704 Tsig =TauUEG/tauu 705 Wsig =(Tsig-One)/(Tsig+One) 706 W1=Wsig 707 W2=Wsig*W1 708 W3=Wsig*W2 709 W4=Wsig*W3 710 W5=Wsig*W4 711 W6=Wsig*W5 712 W7=Wsig*W6 713 W8=Wsig*W7 714 W9=Wsig*W8 715 W10=Wsig*W9 716 W11=Wsig*W10 717 Fsig =at*(at1*W1+ at2*W2 + at3*W3 718 & + at4*W4 + at5*W5 + at6*W6 + at7*W7 719 & + at8*W8 + at9*W9 + at10*W10 + at11*W11) 720 721 722 Gamma2 = delrho(n,1,2)*delrho(n,1,2) + 723 & delrho(n,2,2)*delrho(n,2,2) + 724 & delrho(n,3,2)*delrho(n,3,2) 725 Gamma2 = Four*Gamma2 726 Gamma = dsqrt(Gamma2) 727 if (gamma.lt.DTol) goto 20 728 s = Cs*Gamma/rho43 729 s2 = s*s 730 En = C1*s2 731 Ed = One + C2*s2 732 E = En/Ed 733 Ex = Ex + rho43*Ax*(fL+fNL*E)*(One+Fsig)*0.5d0*qwght(n) 734 if(ldew) func(n)= 735 = func(n)+rho43*Ax*(fL+fNL*E)*(One+Fsig)*.5d0 736c 737c functional derivatives 738c 739 dEn = Two*C1*s 740 dEd = Two*C2*s 741 dE = (dEn*Ed-En*dEd)/(Ed*Ed) 742 dFdW = at*(at1 + Two*at2*W1 + Three*at3*W2 743 & + Four*at4*W3 + Five*at5*W4 744 & + Six*at6*W5 + Seven*at7*W6 745 & + Eight*at8*W7 + Nine*at9*W8 746 & + F10*at10*W9 + F11*at11*W10) 747 dWdT = Two/((One + Tsig)**2) 748 dTdR = Two*(0.5d0*P32*rho13*rho13)/tauu 749 dTdTau = -Two*TauUEG/tauu**2 750 dGGAdR = Two*F4o3*rho13*Ax*(fL+fNL*(E-s*dE)) 751 dFdR = dFdW*dWdT*dTdR 752 dFdTau=dFdW*dWdT*dTdTau 753 dGGAdG =Four*(fNL*dE*s/(Two*Gamma2)) 754 755 Amat(n,2) = Amat(n,2) + (dGGAdR*(One+Fsig) 756 & + (fL+fNL*E)*Ax*rho43*dFdR)*0.5d0 757 Cmat(n,3)= Cmat(n,3) + 758 & dGGAdG*Ax*rho43*(One+Fsig)*0.5d0 759 Mmat(n,2)= Mmat(n,2) + 760 & 0.5d0*0.5d0*rho43*Ax*(fL+fNL*E)*dFdTau 761 762 763c 76420 continue 765 endif 766c 767 return 768 end 769 770 771 772 Subroutine xc_xdldf_d2() 773 call errquit('xdldf: d2 not coded',0,0) 774 return 775 end 776 777 778