1#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 2#if !defined(NWAD_PRINT) 3C> \ingroup nwxc 4C> @{ 5C> 6C> \file nwxc_c_m05.F 7C> Implementation of the M05 correlation functional 8C> 9C> @} 10#endif 11#endif 12C> 13C> \ingroup nwxc_priv 14C> @{ 15C> 16C> \brief The M05 correlation functional 17C> 18C> The M05 functional [1,2] is a meta-GGA of which this evaluates 19C> the correlation component. 20C> 21C> Due to the form of the meta-GGAs we need to screen on the kinetic 22C> energy density to ensure that LDA will be obtained when the kinetic 23C> energy density goes to zero [3]. 24C> 25C> ### References ### 26C> 27C> [1] Y Zhao, NE Schultz, DG Truhlar, 28C> "Exchange-correlation functional with broad accuracy for 29C> metallic and nonmetallic compounds, kinetics, and 30C> noncovalent interactions", 31C> J.Chem.Phys. <b>123</b>, 161103-161106 (2005), DOI: 32C> <a href="https://doi.org/10.1063/1.2126975"> 33C> 10.1063/1.2126975</a>. 34C> 35C> [2] Y Zhao, NE Schultz, DG Truhlar, 36C> "Design of density functionals by combining the method of 37C> constraint satisfaction parametrization for thermochemistry, 38C> thermochemical kinetics, and noncovalent interactions", 39C> J.Chem.Theory.Comput. <b>2</b>, 364-382 (2006), DOI: 40C> <a href="https://doi.org/10.1021/ct0502763"> 41C> 10.1021/ct0502763</a>. 42C> 43C> [3] J. Gräfenstein, D. Izotov, D. Cremer, 44C> "Avoiding singularity problems associated with meta-GGA exchange 45C> and correlation functionals containing the kinetic energy 46C> density", J. Chem. Phys. <b>127</b>, 214103 (2007), DOI: 47C> <a href="https://doi.org/10.1063/1.2800011"> 48C> 10.1063/1.2800011</a>. 49C> 50c M05 or M05-2X exchange functional 51c META GGA 52C utilizes ingredients: 53c rho - density 54c delrho - gradient of density 55c tau - K.S kinetic energy density 56c tauU - uniform-gas KE density 57c ijzy - 1 M05 58c ijzy - 2 M05-2X 59c References: 60c [a] Zhao, Y.; Schultz, N. E.; Truhlar, D. G. J. Chem. Phys. 2005, 123, 161103; 61c Note that in this communication we interchanged cCab,i and cCss,i in Table 1. 62c [b] Zhao, Y.; Schultz, N. E.; Truhlar, D. G. J. Chem. Theory Comput. 2006, in press. 63 64#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 65#if defined(NWAD_PRINT) 66 Subroutine nwxc_x_m05_p(param, tol_rho, ipol, nq, wght, rho, 67 & rgamma, tau, func) 68#else 69 Subroutine nwxc_x_m05(param, tol_rho, ipol, nq, wght, rho, rgamma, 70 & tau, func) 71#endif 72#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV) 73 Subroutine nwxc_x_m05_d2(param, tol_rho, ipol, nq, wght, rho, 74 & rgamma, tau, func) 75#else 76 Subroutine nwxc_x_m05_d3(param, tol_rho, ipol, nq, wght, rho, 77 & rgamma, tau, func) 78#endif 79c 80c$Id$ 81c 82#include "nwad.fh" 83c 84 implicit none 85c 86#include "nwxc_param.fh" 87c 88c 89c Input and other parameters 90c 91#if defined(NWAD_PRINT) 92#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 93 type(nwad_dble)::param(*) !< [Input] Parameters of functional 94 type(nwad_dble)::at1, at2, at3, at4, at5, at6, at7, at8, at9 95 type(nwad_dble)::at10, at11 96#else 97 double precision param(*) !< [Input] Parameters of functional 98 double precision at1, at2, at3, at4, at5, at6, at7, at8, at9 99 double precision at10, at11 100#endif 101#else 102 double precision param(*) !< [Input] Parameters of functional 103 !< - param(1): \f$ a_{1} \f$ 104 !< - param(2): \f$ a_{2} \f$ 105 !< - param(3): \f$ a_{3} \f$ 106 !< - param(4): \f$ a_{4} \f$ 107 !< - param(5): \f$ a_{5} \f$ 108 !< - param(6): \f$ a_{6} \f$ 109 !< - param(7): \f$ a_{7} \f$ 110 !< - param(8): \f$ a_{8} \f$ 111 !< - param(9): \f$ a_{9} \f$ 112 !< - param(10): \f$ a_{10} \f$ 113 !< - param(11): \f$ a_{11} \f$ 114 double precision at1, at2, at3, at4, at5, at6, at7, at8, at9 115 double precision at10, at11 116#endif 117 double precision tol_rho !< [Input] The lower limit on the density 118 integer nq !< [Input] The number of points 119 integer ipol !< [Input] The number of spin channels 120 double precision wght !< [Input] The weight of the functional 121c 122c Charge Density 123c 124 type(nwad_dble)::rho(nq,*) !< [Input] The density 125c 126c Charge Density Gradient Norm 127c 128 type(nwad_dble)::rgamma(nq,*) !< [Input] The density gradient norm 129c 130c Kinetic Energy Density 131c 132 type(nwad_dble)::tau(nq,*) !< [Input] The kinetic energy density 133c 134c Functional values 135c 136 type(nwad_dble)::func(*) !< [Output] The functional value 137c 138c Sampling Matrices for the XC Potential 139c 140c double precision Amat(nq,*) !< [Output] Derivative wrt density 141c double precision Cmat(nq,*) !< [Output] Derivative wrt rgamma 142c double precision Mmat(nq,*) !< [Output] Derivative wrt tau 143c 144 double precision pi 145c 146 integer n 147 double precision C1, C2, fL, fNL, at 148 type(nwad_dble)::rrho, rho43, rho13, rhoo, rho53 149 type(nwad_dble)::Gamma2, Gamma 150 type(nwad_dble)::TauUEG, Tsig, Wsig, W1, W2, W3, W4, W5, W6 151 type(nwad_dble)::W7, W8, W9, W10, W11, Fsig 152c 153c kinetic energy density or tau 154c 155 type(nwad_dble)::tauN,tauu 156 double precision DTol 157 double precision F83, F23, F53, F43, F13, F1o3 158 double precision F1o4, F2o3, F3o2, F4o3, F4o9, F3o5 159 double precision One, Two, Three, Four, Five, Six, Seven, Eight 160 double precision Nine, F10, F11 161 type(nwad_dble)::E,En,Ed,s,s2 162 double precision Cs, Ax, P32, dE, dEn, dEd 163 164c functional derivatives below FFFFFFFFFFFF 165 166 double precision dFdW, dWdT, dTdR, dTdTau, dGGAdR, dFdR 167 double precision dFdTau, dGGAdG,tauW 168 169c functional derivatives above FFFFFFFFFFFF 170 171 172 parameter( pi = 3.1415926535897932384626433832795d0 ) 173 174 175 parameter (F1o3=1.d0/3.d0, F1o4=1.d0/4.d0, F2o3=2.d0/3.d0, 176 & F3o2=3.d0/2.d0,F13=1.d0/3.d0) 177 parameter (F4o3=4.d0/3.d0, F4o9=4.d0/9.d0, F3o5=3.d0/5.d0) 178 parameter (F83=8.d0/3.0d0, F23=2.0d0/3.d0, F53=5.d0/3.d0) 179 parameter (One=1.0d0, Two=2.0d0, Three=3.0d0, Four=4.0d0, 180 & Five=5.0d0,Six=6.0d0, Seven=7.0d0, 181 & Eight=8.0d0, Nine=9.0d0,F10=10.d0, F11=11.d0) 182 183 184 at1= param(1) 185 at2= param(2) 186 at3= param(3) 187 at4= param(4) 188 at5= param(5) 189 at6= param(6) 190 at7= param(7) 191 at8= param(8) 192 at9= param(9) 193 at10= param(10) 194 at11= param(11) 195c if (ijzy.eq.1) then 196c Parameters for M05 197c at1= 0.08151d0 198c at2= -0.43956d0 199c at3= -3.22422d0 200c at4= 2.01819d0 201c at5= 8.79431d0 202c at6= -0.00295d0 203c at7= 9.82029d0 204c at8= -4.82351d0 205c at9= -48.17574d0 206c at10= 3.64802d0 207c at11= 34.02248d0 208c elseif (ijzy.eq.2) then 209c Parameters for M05-2X 210c at1= -0.56833d0 211c at2= -1.30057d0 212c at3= 5.50070d0 213c at4= 9.06402d0 214c at5= -32.21075d0 215c at6= -23.73298d0 216c at7= 70.22996d0 217c at8= 29.88614d0 218c at9= -60.25778d0 219c at10= -13.22205d0 220c at11= 15.23694d0 221c endif 222 223 at=1.0d0 224 C1 = 0.2195149727645171d0 225 C2 = C1/0.804d0 226 DTol=tol_rho 227C 228C Scale factors for local and non-local contributions. 229C 230 fL = wght 231 fNL = wght 232 Cs = 0.5d0/(3.0d0*pi*pi)**F13 233 P32 = (3.d0*pi**2)**F23 234 235c 236 Ax = (-0.75d0)*(3.0d0/pi)**F13 237c 238 if (ipol.eq.1 )then 239c 240c ======> SPIN-RESTRICTED <====== 241c or 242c SPIN-UNPOLARIZED 243c 244c 245 do 10 n = 1, nq 246 247 if (rho(n,R_T).lt.DTol) goto 10 248 rhoo = rho(n,R_T) 249 rho43 = rhoo**F4o3 250 rrho = 1d0/rhoo ! reciprocal of rho 251 rho13 = rho43*rrho 252 rho53 = rhoo**F53 253c 254 Gamma2 = rgamma(n,G_TT) 255c Gamma = sqrt(Gamma2) 256 tauN = tau(n,T_T) 257 tauu=tauN 258 TauUEG=0.3d0*P32*rho53 259c Tsig =TauUEG/tauu 260c Wsig =(Tsig-One)/(Tsig+One) 261 Wsig =(TauUEG-tauu)/(TauUEG+tauu) 262 W1=Wsig 263 W2=Wsig*W1 264 W3=Wsig*W2 265 W4=Wsig*W3 266 W5=Wsig*W4 267 W6=Wsig*W5 268 W7=Wsig*W6 269 W8=Wsig*W7 270 W9=Wsig*W8 271 W10=Wsig*W9 272 W11=Wsig*W10 273 Fsig =at*(at1*W1+ at2*W2 + at3*W3 274 & + at4*W4 + at5*W5 + at6*W6 + at7*W7 275 & + at8*W8 + at9*W9 + at10*W10 + at11*W11) 276 277c s = Cs*Gamma/rho43 278c s2 = s*s 279 s2 = Cs*Cs*Gamma2/(rho43*rho43) 280 En = C1*s2 281 Ed = One + C2*s2 282 E = En/Ed 283c Ex = Ex + rho43*Ax*(fL+fNL*E)*(One+Fsig)*qwght(n) 284 func(n)=func(n)+rho43*Ax*(fL+fNL*E)*(One+Fsig) 285c 286c functional derivatives 287c 288c dEn = Two*C1*s 289c dEd = Two*C2*s 290c dE = (dEn*Ed-En*dEd)/(Ed*Ed) 291c dFdW = at*(at1 + Two*at2*W1 + Three*at3*W2 292c & + Four*at4*W3 + Five*at5*W4 293c & + Six*at6*W5 + Seven*at7*W6 294c & + Eight*at8*W7 + Nine*at9*W8 295c & + F10*at10*W9 + F11*at11*W10) 296c dWdT = Two/((One + Tsig)**2) 297c dTdR = (0.5d0*P32*rho13*rho13)/tauu 298c dTdTau = -TauUEG/tauu**2 299c dGGAdR = F4o3*rho13*Ax*(fL+fNL*(E-s*dE)) 300c dFdR = dFdW*dWdT*dTdR 301c dFdTau=dFdW*dWdT*dTdTau 302c dGGAdG =(fNL*dE*s/(Two*Gamma2)) 303c Amat(n,D1_RA) = Amat(n,D1_RA) + dGGAdR*(One+Fsig) 304c & + (fL+fNL*E)*Ax*rho43*dFdR 305c Cmat(n,D1_GAA)= Cmat(n,D1_GAA) + 306c & Two*dGGAdG*Ax*rho43*(One+Fsig) 307c Mmat(n,D1_TA)= Mmat(n,D1_TA) 308c & + 0.5d0*rho43*Ax*(fL+fNL*E)*dFdTau 309 310 311 31210 continue 313 314 315c UUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUnrestricted 316 else ! ipol=2 317c 318c ======> SPIN-UNRESTRICTED <====== 319 320c 321c use spin density functional theory ie n-->2n 322c Ex=(1/2)Ex[2*alpha] + (1/2)Ex[2*beta] 323 324 do 20 n = 1, nq 325 if (rho(n,R_A)+rho(n,R_B).lt.DTol) goto 20 326c 327c Alpha ALPHA ALPHA 328c 329 if (rho(n,R_A).lt.0.5d0*DTol) goto 25 330 rhoo = Two*rho(n,R_A) 331 rho43 = rhoo**F4o3 332 rrho = 1.0d0/rhoo ! reciprocal of rho 333 rho13 = rho43*rrho 334 rho53 = rhoo**F53 335 336c 337 338 tauN = tau(n,T_A) 339 tauu = Two*tauN 340 TauUEG=0.3d0*P32*rho53 341c Tsig =TauUEG/tauu 342c Wsig =(Tsig-One)/(Tsig+One) 343 Wsig =(TauUEG-tauu)/(TauUEG+tauu) 344 W1=Wsig 345 W2=Wsig*W1 346 W3=Wsig*W2 347 W4=Wsig*W3 348 W5=Wsig*W4 349 W6=Wsig*W5 350 W7=Wsig*W6 351 W8=Wsig*W7 352 W9=Wsig*W8 353 W10=Wsig*W9 354 W11=Wsig*W10 355 Fsig =at*(at1*W1+ at2*W2 + at3*W3 356 & + at4*W4 + at5*W5 + at6*W6 + at7*W7 357 & + at8*W8 + at9*W9 + at10*W10 + at11*W11) 358 359 360c Gamma2 = delrho(n,1,1)*delrho(n,1,1) + 361c & delrho(n,2,1)*delrho(n,2,1) + 362c & delrho(n,3,1)*delrho(n,3,1) 363 Gamma2 = rgamma(n,G_AA) 364 Gamma2 = Four*Gamma2 365c Gamma = sqrt(Gamma2) 366 367c s = Cs*Gamma/rho43 368c s2 = s*s 369 s2 = Cs*Cs*Gamma2/(rho43*rho43) 370 En = C1*s2 371 Ed = One + C2*s2 372 E = En/Ed 373c Ex = Ex + rho43*Ax*(fL+fNL*E)*(One+Fsig)*0.5d0*qwght(n) 374 func(n)= func(n)+rho43*Ax*(fL+fNL*E)*(One+Fsig)*.5d0 375c 376c functional derivatives 377c 378c dEn = Two*C1*s 379c dEd = Two*C2*s 380c dE = (dEn*Ed-En*dEd)/(Ed*Ed) 381c dFdW = at*(at1 + Two*at2*W1 + Three*at3*W2 382c & + Four*at4*W3 + Five*at5*W4 383c & + Six*at6*W5 + Seven*at7*W6 384c & + Eight*at8*W7 + Nine*at9*W8 385c & + F10*at10*W9 + F11*at11*W10) 386c dWdT = Two/((One + Tsig)**2) 387c dTdR = Two*(0.5d0*P32*rho13*rho13)/tauu 388c dTdTau = -Two*TauUEG/tauu**2 389c dGGAdR = Two*F4o3*rho13*Ax*(fL+fNL*(E-s*dE)) 390c dFdR = dFdW*dWdT*dTdR 391c dFdTau=dFdW*dWdT*dTdTau 392c dGGAdG =Four*(fNL*dE*s/(Two*Gamma2)) 393 394c Amat(n,D1_RA) = Amat(n,D1_RA) + (dGGAdR*(One+Fsig) 395c & + (fL+fNL*E)*Ax*rho43*dFdR)*0.5d0 396c Cmat(n,D1_GAA)= Cmat(n,D1_GAA) + 397c & dGGAdG*Ax*rho43*(One+Fsig)*0.5d0 398c Mmat(n,D1_TA)= Mmat(n,D1_TA) + 399c & 0.5d0*0.5d0*rho43*Ax*(fL+fNL*E)*dFdTau 400 401c write (*,*) "Ex,Amat(n,1),Cmat(n,1)", 402c & Ex,Amat(n,1),Cmat(n,1) 403 404c 405c Beta BETA BETA 406c 407 40825 continue 409 410c 411c Beta 412c 413 if (rho(n,R_B).lt.0.5d0*DTol) goto 20 414 rhoo = Two*rho(n,R_B) 415 rho43 = rhoo**F4o3 416 rrho = 1.0d0/rhoo ! reciprocal of rho 417 rho13 = rho43*rrho 418 rho53 = rhoo**F53 419 420c 421 422 tauN = tau(n,T_B) 423 tauu = Two*tauN 424 TauUEG=0.3d0*P32*rho53 425c Tsig =TauUEG/tauu 426c Wsig =(Tsig-One)/(Tsig+One) 427 Wsig =(TauUEG-tauu)/(TauUEG+tauu) 428 W1=Wsig 429 W2=Wsig*W1 430 W3=Wsig*W2 431 W4=Wsig*W3 432 W5=Wsig*W4 433 W6=Wsig*W5 434 W7=Wsig*W6 435 W8=Wsig*W7 436 W9=Wsig*W8 437 W10=Wsig*W9 438 W11=Wsig*W10 439 Fsig =at*(at1*W1+ at2*W2 + at3*W3 440 & + at4*W4 + at5*W5 + at6*W6 + at7*W7 441 & + at8*W8 + at9*W9 + at10*W10 + at11*W11) 442 443 444c Gamma2 = delrho(n,1,2)*delrho(n,1,2) + 445c & delrho(n,2,2)*delrho(n,2,2) + 446c & delrho(n,3,2)*delrho(n,3,2) 447 Gamma2 = rgamma(n,G_BB) 448 Gamma2 = Four*Gamma2 449c Gamma = sqrt(Gamma2) 450c s = Cs*Gamma/rho43 451c s2 = s*s 452 s2 = Cs*Cs*Gamma2/(rho43*rho43) 453 En = C1*s2 454 Ed = One + C2*s2 455 E = En/Ed 456c Ex = Ex + rho43*Ax*(fL+fNL*E)*(One+Fsig)*0.5d0*qwght(n) 457 func(n)= func(n)+rho43*Ax*(fL+fNL*E)*(One+Fsig)*.5d0 458c 459c functional derivatives 460c 461c dEn = Two*C1*s 462c dEd = Two*C2*s 463c dE = (dEn*Ed-En*dEd)/(Ed*Ed) 464c dFdW = at*(at1 + Two*at2*W1 + Three*at3*W2 465c & + Four*at4*W3 + Five*at5*W4 466c & + Six*at6*W5 + Seven*at7*W6 467c & + Eight*at8*W7 + Nine*at9*W8 468c & + F10*at10*W9 + F11*at11*W10) 469c dWdT = Two/((One + Tsig)**2) 470c dTdR = Two*(0.5d0*P32*rho13*rho13)/tauu 471c dTdTau = -Two*TauUEG/tauu**2 472c dGGAdR = Two*F4o3*rho13*Ax*(fL+fNL*(E-s*dE)) 473c dFdR = dFdW*dWdT*dTdR 474c dFdTau=dFdW*dWdT*dTdTau 475c dGGAdG =Four*(fNL*dE*s/(Two*Gamma2)) 476 477c Amat(n,D1_RB) = Amat(n,D1_RB) + (dGGAdR*(One+Fsig) 478c & + (fL+fNL*E)*Ax*rho43*dFdR)*0.5d0 479c Cmat(n,D1_GBB)= Cmat(n,D1_GBB) + 480c & dGGAdG*Ax*rho43*(One+Fsig)*0.5d0 481c Mmat(n,D1_TB)= Mmat(n,D1_TB) + 482c & 0.5d0*0.5d0*rho43*Ax*(fL+fNL*E)*dFdTau 483 484 485c 48620 continue 487 endif 488c 489 return 490 end 491 492 493c---------------------------------------------------------------------- 494C> \brief Calculate the dlDF correlation functional 495C> 496C> Calculate the dlDF exchange functional [1]. 497C> 498C> ### References ### 499C> 500C> [1] K Pernal, R Podeszwa, K Patkowski, K Szalewicz, 501C> "Dispersionless density functional theory", 502C> Phys.Rev.Lett. <b>103</b>, 263201-263204 (2009), DOI: 503C> <a href="https://doi.org/10.1103/PhysRevLett.103.263201"> 504C> 10.1103/PhysRevLett.103.263201</a>. 505C> 506c dlDF exchange functional 507c META GGA 508C utilizes ingredients: 509c rho - density 510c delrho - gradient of density 511c tau - K.S kinetic energy density 512c tauU - uniform-gas KE density 513c 514c References: 515c [a] Pernal,Podeszwa,Patkowski,Szalewicz, PRL 103 263201 (2009) 516#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 517#if defined(NWAD_PRINT) 518 Subroutine nwxc_x_dldf_p(tol_rho, ipol, nq, wght, rho, rgamma, 519 & tau, func) 520#else 521 Subroutine nwxc_x_dldf(tol_rho, ipol, nq, wght, rho, rgamma, tau, 522 & func) 523#endif 524#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV) 525 Subroutine nwxc_x_dldf_d2(tol_rho, ipol, nq, wght, rho, rgamma, 526 & tau, func) 527#else 528 Subroutine nwxc_x_dldf_d3(tol_rho, ipol, nq, wght, rho, rgamma, 529 & tau, func) 530#endif 531c 532#include "nwad.fh" 533c 534 implicit none 535c 536#include "nwxc_param.fh" 537c 538c 539c Input and other parameters 540c 541 double precision tol_rho !< [Input] The lower limit on the density 542 integer nq !< [Input] The number of points 543 integer ipol !< [Input] The number of spin channels 544 double precision wght !< [Input] The weight of the functional 545c 546c Charge Density 547c 548 type(nwad_dble)::rho(nq,*) !< [Input] The density 549c 550c Charge Density Gradient Norm 551c 552 type(nwad_dble)::rgamma(nq,*) !< [Input] The density gradient norm 553c 554c Kinetic Energy Density 555c 556 type(nwad_dble)::tau(nq,*) !< [Input] The kinetic energy density 557c 558c Functional values 559c 560 type(nwad_dble)::func(*) !< [Output] The functional value 561c 562c Sampling Matrices for the XC Potential 563c 564c double precision Amat(nq,*) !< [Output] Derivative wrt density 565c double precision Cmat(nq,*) !< [Output] Derivative wrt rgamma 566c double precision Mmat(nq,*) !< [Output] Derivative wrt tau 567c 568 double precision pi 569c 570 integer n 571 double precision at1, at2, at3, at4, at5, at6, at7, at8, at9 572 double precision at, at10, at11, C1, C2, fL, fNL 573 type(nwad_dble)::rrho, rho43, rho13, rhoo, rho53 574c type(nwad_dble)::Gamma2, Gamma 575c type(nwad_dble)::TauUEG, Tsig, Wsig, W1, W2, W3, W4, W5, W6 576 type(nwad_dble)::Gamma2 577 type(nwad_dble)::TauUEG, Wsig, W1, W2, W3, W4, W5, W6 578 type(nwad_dble)::W7, W8, W9, W10, W11, Fsig 579c 580c kinetic energy density or tau 581c 582 type(nwad_dble)::tauN,tauu 583 double precision DTol 584 double precision F83, F23, F53, F43, F13, F1o3 585 double precision F1o4, F2o3, F3o2, F4o3, F4o9, F3o5 586 double precision One, Two, Three, Four, Five, Six, Seven, Eight 587 double precision Nine, F10, F11 588 type(nwad_dble)::E,En,Ed,s,s2 589 double precision Cs, Ax, P32, dE, dEn, dEd 590 591c functional derivatives below FFFFFFFFFFFF 592 593 double precision dFdW, dWdT, dTdR, dTdTau, dGGAdR, dFdR 594 double precision dFdTau, dGGAdG,tauW 595 596c functional derivatives above FFFFFFFFFFFF 597 598 599 parameter( pi = 3.1415926535897932384626433832795d0 ) 600 601 602 parameter (F1o3=1.d0/3.d0, F1o4=1.d0/4.d0, F2o3=2.d0/3.d0, 603 & F3o2=3.d0/2.d0,F13=1.d0/3.d0) 604 parameter (F4o3=4.d0/3.d0, F4o9=4.d0/9.d0, F3o5=3.d0/5.d0) 605 parameter (F83=8.d0/3.0d0, F23=2.0d0/3.d0, F53=5.d0/3.d0) 606 parameter (One=1.0d0, Two=2.0d0, Three=3.0d0, Four=4.0d0, 607 & Five=5.0d0,Six=6.0d0, Seven=7.0d0, 608 & Eight=8.0d0, Nine=9.0d0,F10=10.d0, F11=11.d0) 609 610c Parameters for dlDF 611 at1= -1.637571d-01 612 at2= -1.880028d-01 613 at3= -4.490609d-01 614 at4= -8.2359d-03 615 at5= 0.0d0 616 at6= 0.0d0 617 at7= 0.0d0 618 at8= 0.0d0 619 at9= 0.0d0 620 at10= 0.0d0 621 at11= 0.0d0 622 623 624 at=1.0d0 625 C1 = 0.3511128d0 626 C2 = C1/4.8827323d0 627 DTol=tol_rho 628C 629C Scale factors for local and non-local contributions. 630C 631 fL = wght 632 fNL = wght 633 Cs = 0.5d0/(3.0d0*pi*pi)**F13 634 P32 = (3.d0*pi**2)**F23 635 636c 637 Ax = (-0.75d0)*(3.0d0/pi)**F13 638 639 640c 641 if (ipol.eq.1 )then 642c 643c ======> SPIN-RESTRICTED <====== 644c or 645c SPIN-UNPOLARIZED 646c 647c 648 do 10 n = 1, nq 649 650 if (rho(n,R_T).lt.DTol) goto 10 651 rhoo = rho(n,R_T) 652 rho43 = rhoo**F4o3 653 rrho = 1d0/rhoo ! reciprocal of rho 654 rho13 = rho43*rrho 655 rho53 = rhoo**F53 656 657c 658 659 tauN = tau(n,T_T) 660 tauu=tauN 661cedo if (taun.lt.sqrt(DTol)) goto 10 662c Gamma2 = delrho(n,1,1)*delrho(n,1,1) + 663c & delrho(n,2,1)*delrho(n,2,1) + 664c & delrho(n,3,1)*delrho(n,3,1) 665 Gamma2 = rgamma(n,G_TT) 666c Gamma = sqrt(Gamma2) 667 TauUEG=0.3d0*P32*rho53 668c Tsig =TauUEG/tauu 669c Wsig =(Tsig-One)/(Tsig+One) 670 Wsig =(TauUEG-tauu)/(TauUEG+tauu) 671 W1=Wsig 672 W2=Wsig*W1 673 W3=Wsig*W2 674 W4=Wsig*W3 675 W5=Wsig*W4 676 W6=Wsig*W5 677 W7=Wsig*W6 678 W8=Wsig*W7 679 W9=Wsig*W8 680 W10=Wsig*W9 681 W11=Wsig*W10 682 Fsig =at*(at1*W1+ at2*W2 + at3*W3 683 & + at4*W4 + at5*W5 + at6*W6 + at7*W7 684 & + at8*W8 + at9*W9 + at10*W10 + at11*W11) 685 686c s = Cs*Gamma/rho43 687c s2 = s*s 688 s2 = Cs*Cs*Gamma2/(rho43*rho43) 689 En = C1*s2 690 Ed = One + C2*s2 691 E = En/Ed 692c Ex = Ex + rho43*Ax*(fL+fNL*E)*(One+Fsig)*qwght(n) 693 func(n)=func(n)+rho43*Ax*(fL+fNL*E)*(One+Fsig) 694c 695c functional derivatives 696c 697c dEn = Two*C1*s 698c dEd = Two*C2*s 699c dE = (dEn*Ed-En*dEd)/(Ed*Ed) 700c dFdW = at*(at1 + Two*at2*W1 + Three*at3*W2 701c & + Four*at4*W3 + Five*at5*W4 702c & + Six*at6*W5 + Seven*at7*W6 703c & + Eight*at8*W7 + Nine*at9*W8 704c & + F10*at10*W9 + F11*at11*W10) 705c dWdT = Two/((One + Tsig)**2) 706c dTdR = (0.5d0*P32*rho13*rho13)/tauu 707c dTdTau = -TauUEG/tauu**2 708c dGGAdR = F4o3*rho13*Ax*(fL+fNL*(E-s*dE)) 709c dFdR = dFdW*dWdT*dTdR 710c dFdTau=dFdW*dWdT*dTdTau 711c dGGAdG =(fNL*dE*s/(Two*Gamma2)) 712c Amat(n,D1_RA) = Amat(n,D1_RA) + dGGAdR*(One+Fsig) 713c & + (fL+fNL*E)*Ax*rho43*dFdR 714c Cmat(n,D1_GAA)= Cmat(n,D1_GAA) + 715c & Two*dGGAdG*Ax*rho43*(One+Fsig) 716c Mmat(n,D1_TA)= Mmat(n,D1_TA) 717c & + 0.5d0*rho43*Ax*(fL+fNL*E)*dFdTau 718 719 720 72110 continue 722 723 724c UUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUnrestricted 725 else ! ipol=2 726c 727c ======> SPIN-UNRESTRICTED <====== 728 729c 730c use spin density functional theory ie n-->2n 731c Ex=(1/2)Ex[2*alpha] + (1/2)Ex[2*beta] 732 733 do 20 n = 1, nq 734 if (rho(n,R_A)+rho(n,R_B).lt.DTol) goto 20 735c 736c Alpha ALPHA ALPHA 737c 738 if (rho(n,R_A).lt.0.5d0*DTol) goto 25 739 rhoo = Two*rho(n,R_A) 740 rho43 = rhoo**F4o3 741 rrho = 1.0d0/rhoo ! reciprocal of rho 742 rho13 = rho43*rrho 743 rho53 = rhoo**F53 744 745c 746 747 tauN = tau(n,T_A) 748 tauu = Two*tauN 749 TauUEG=0.3d0*P32*rho53 750c Tsig =TauUEG/tauu 751c Wsig =(Tsig-One)/(Tsig+One) 752 Wsig = (TauUEG-tauu)/(TauUEG+tauu) 753 W1=Wsig 754 W2=Wsig*W1 755 W3=Wsig*W2 756 W4=Wsig*W3 757 W5=Wsig*W4 758 W6=Wsig*W5 759 W7=Wsig*W6 760 W8=Wsig*W7 761 W9=Wsig*W8 762 W10=Wsig*W9 763 W11=Wsig*W10 764 Fsig =at*(at1*W1+ at2*W2 + at3*W3 765 & + at4*W4 + at5*W5 + at6*W6 + at7*W7 766 & + at8*W8 + at9*W9 + at10*W10 + at11*W11) 767 768 769c Gamma2 = delrho(n,1,1)*delrho(n,1,1) + 770c & delrho(n,2,1)*delrho(n,2,1) + 771c & delrho(n,3,1)*delrho(n,3,1) 772 Gamma2 = rgamma(n,G_AA) 773 Gamma2 = Four*Gamma2 774c Gamma = sqrt(Gamma2) 775c s = Cs*Gamma/rho43 776c s2 = s*s 777 s2 = Cs*Cs*Gamma2/(rho43*rho43) 778 En = C1*s2 779 Ed = One + C2*s2 780 E = En/Ed 781c Ex = Ex + rho43*Ax*(fL+fNL*E)*(One+Fsig)*0.5d0*qwght(n) 782 func(n)= func(n)+rho43*Ax*(fL+fNL*E)*(One+Fsig)*.5d0 783c 784c functional derivatives 785c 786c dEn = Two*C1*s 787c dEd = Two*C2*s 788c dE = (dEn*Ed-En*dEd)/(Ed*Ed) 789c dFdW = at*(at1 + Two*at2*W1 + Three*at3*W2 790c & + Four*at4*W3 + Five*at5*W4 791c & + Six*at6*W5 + Seven*at7*W6 792c & + Eight*at8*W7 + Nine*at9*W8 793c & + F10*at10*W9 + F11*at11*W10) 794c dWdT = Two/((One + Tsig)**2) 795c dTdR = Two*(0.5d0*P32*rho13*rho13)/tauu 796c dTdTau = -Two*TauUEG/tauu**2 797c dGGAdR = Two*F4o3*rho13*Ax*(fL+fNL*(E-s*dE)) 798c dFdR = dFdW*dWdT*dTdR 799c dFdTau=dFdW*dWdT*dTdTau 800c dGGAdG =Four*(fNL*dE*s/(Two*Gamma2)) 801 802c Amat(n,D1_RA) = Amat(n,D1_RA) + (dGGAdR*(One+Fsig) 803c & + (fL+fNL*E)*Ax*rho43*dFdR)*0.5d0 804c Cmat(n,D1_GAA)= Cmat(n,D1_GAA) + 805c & dGGAdG*Ax*rho43*(One+Fsig)*0.5d0 806c Mmat(n,D1_TA)= Mmat(n,D1_TA) + 807c & 0.5d0*0.5d0*rho43*Ax*(fL+fNL*E)*dFdTau 808 809c write (*,*) "Ex,Amat(n,1),Cmat(n,1)", 810c & Ex,Amat(n,1),Cmat(n,1) 811 812c 813c Beta BETA BETA 814c 815 81625 continue 817 818c 819c Beta 820c 821 if (rho(n,R_B).lt.0.5d0*DTol) goto 20 822 rhoo = Two*rho(n,R_B) 823 rho43 = rhoo**F4o3 824 rrho = 1.0d0/rhoo ! reciprocal of rho 825 rho13 = rho43*rrho 826 rho53 = rhoo**F53 827c 828 tauN = tau(n,T_B) 829 tauu = Two*tauN 830 TauUEG=0.3d0*P32*rho53 831c Tsig =TauUEG/tauu 832c Wsig =(Tsig-One)/(Tsig+One) 833 Wsig =(TauUEG-tauu)/(TauUEG+tauu) 834 W1=Wsig 835 W2=Wsig*W1 836 W3=Wsig*W2 837 W4=Wsig*W3 838 W5=Wsig*W4 839 W6=Wsig*W5 840 W7=Wsig*W6 841 W8=Wsig*W7 842 W9=Wsig*W8 843 W10=Wsig*W9 844 W11=Wsig*W10 845 Fsig =at*(at1*W1+ at2*W2 + at3*W3 846 & + at4*W4 + at5*W5 + at6*W6 + at7*W7 847 & + at8*W8 + at9*W9 + at10*W10 + at11*W11) 848 849 850c Gamma2 = delrho(n,1,2)*delrho(n,1,2) + 851c & delrho(n,2,2)*delrho(n,2,2) + 852c & delrho(n,3,2)*delrho(n,3,2) 853 Gamma2 = rgamma(n,G_BB) 854 Gamma2 = Four*Gamma2 855c Gamma = sqrt(Gamma2) 856c s = Cs*Gamma/rho43 857c s2 = s*s 858 s2 = Cs*Cs*Gamma2/(rho43*rho43) 859 En = C1*s2 860 Ed = One + C2*s2 861 E = En/Ed 862c Ex = Ex + rho43*Ax*(fL+fNL*E)*(One+Fsig)*0.5d0*qwght(n) 863 func(n)=func(n)+rho43*Ax*(fL+fNL*E)*(One+Fsig)*.5d0 864c 865c functional derivatives 866c 867c dEn = Two*C1*s 868c dEd = Two*C2*s 869c dE = (dEn*Ed-En*dEd)/(Ed*Ed) 870c dFdW = at*(at1 + Two*at2*W1 + Three*at3*W2 871c & + Four*at4*W3 + Five*at5*W4 872c & + Six*at6*W5 + Seven*at7*W6 873c & + Eight*at8*W7 + Nine*at9*W8 874c & + F10*at10*W9 + F11*at11*W10) 875c dWdT = Two/((One + Tsig)**2) 876c dTdR = Two*(0.5d0*P32*rho13*rho13)/tauu 877c dTdTau = -Two*TauUEG/tauu**2 878c dGGAdR = Two*F4o3*rho13*Ax*(fL+fNL*(E-s*dE)) 879c dFdR = dFdW*dWdT*dTdR 880c dFdTau=dFdW*dWdT*dTdTau 881c dGGAdG =Four*(fNL*dE*s/(Two*Gamma2)) 882 883c Amat(n,D1_RB) = Amat(n,D1_RB) + (dGGAdR*(One+Fsig) 884c & + (fL+fNL*E)*Ax*rho43*dFdR)*0.5d0 885c Cmat(n,D1_GBB)= Cmat(n,D1_GBB) + 886c & dGGAdG*Ax*rho43*(One+Fsig)*0.5d0 887c Mmat(n,D1_TB)= Mmat(n,D1_TB) + 888c & 0.5d0*0.5d0*rho43*Ax*(fL+fNL*E)*dFdTau 889 890 891c 89220 continue 893 endif 894c 895 return 896 end 897 898#if !defined(NWAD_PRINT) 899#define NWAD_PRINT 900c 901c Compile source again for Maxima 902c 903#include "nwxc_x_m05.F" 904#endif 905#if !defined(SECOND_DERIV) 906#define SECOND_DERIV 907c 908c Compile source again for the 2nd derivative case 909c 910#include "nwxc_x_m05.F" 911#endif 912#if !defined(THIRD_DERIV) 913#define THIRD_DERIV 914c 915c Compile source again for the 3rd derivative case 916c 917#include "nwxc_x_m05.F" 918#endif 919#undef NWAD_PRINT 920C> @} 921 922 923