1#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 2#if !defined(NWAD_PRINT) 3C> \ingroup nwxc 4C> @{ 5C> 6C> \file nwxc_x_m06.F 7C> Implementation of the M06 exchange functional 8C> 9C> @} 10#endif 11#endif 12C> 13C> \ingroup nwxc_priv 14C> @{ 15C> 16C> \brief The M06 exchange functional 17C> 18C> The M06 functional [1,2] is a meta-GGA of which this evaluates 19C> the exchange component. 20C> 21C> ### References ### 22C> 23C> [1] Y Zhao, DG Truhlar, 24C> "A new local density functional for main-group thermochemistry, 25C> transition metal bonding, thermochemical kinetics, and noncovalent 26C> interactions", 27C> J. Chem. Phys. <b>125</b>, 194101 (2006), DOI: 28C> <a href="https://doi.org/10.1063/1.2370993"> 29C> 10.1063/1.2370993</a>. 30C> 31C> [2] Y Zhao, DG Truhlar, 32C> "Density functional for spectroscopy: No long-range self-interaction 33C> error, good performance for Rydberg and charge-transfer states, 34C> and better performance on average than B3LYP for ground states", 35C> J. Phys. Chem. A <b>110</b>, 13126-13130 (2006), DOI: 36C> <a href="https://doi.org/10.1021/jp066479k"> 37C> 10.1021/jp066479k</a>. 38C> 39c M06 suite exchange functional 40c META GGA 41C utilizes ingredients: 42c rho - density 43c delrho - gradient of density 44c tau - K.S kinetic energy density 45c tauU - uniform-gas KE density 46c ijzy - 1 M06-L 47c ijzy - 2 M06-HF 48c ijzy - 3 M06 49c ijzy - 4 M06-2X 50c References: 51c [a] Zhao, Y. and Truhlar, D. G. J. Chem. Phys. 2006, 125, 194101; 52c [b] Zhao, Y. and Truhlar, D. G. J. Phys. Chem. A (2006),110(49),13126-13130. 53 54#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 55#if defined(NWAD_PRINT) 56 Subroutine nwxc_x_m06_p(param, tol_rho, ipol, nq, wght, rho, 57 & rgamma, tau, func) 58#else 59 Subroutine nwxc_x_m06(param, tol_rho, ipol, nq, wght, rho, rgamma, 60 & tau, func) 61#endif 62#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV) 63 Subroutine nwxc_x_m06_d2(param, tol_rho, ipol, nq, wght, rho, 64 & rgamma, tau, func) 65#else 66 Subroutine nwxc_x_m06_d3(param, tol_rho, ipol, nq, wght, rho, 67 & rgamma, tau, func) 68#endif 69c 70c$Id$ 71c 72#include "nwad.fh" 73c 74 implicit none 75c 76#include "nwxc_param.fh" 77c 78#if defined(NWAD_PRINT) 79#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 80 type(nwad_dble)::param(*) !< [Input] Parameters of the functional 81 type(nwad_dble)::at1, at2, at3, at4, at5, at6, at7, at8, at9 82 type(nwad_dble)::at10, at11, at0 83#else 84 double precision param(*) !< [Input] Parameters of the functional 85 double precision at1, at2, at3, at4, at5, at6, at7, at8, at9 86 double precision at10, at11, at0 87#endif 88#else 89 double precision param(*) !< [Input] Parameters of the functional 90 !< - param( 1): \f$ d_0 \f$ 91 !< - param( 2): \f$ d_1 \f$ 92 !< - param( 3): \f$ d_2 \f$ 93 !< - param( 4): \f$ d_3 \f$ 94 !< - param( 5): \f$ d_4 \f$ 95 !< - param( 6): \f$ d_5 \f$ 96 !< - param( 7): \f$ a_0 \f$ 97 !< - param( 8): \f$ a_1 \f$ 98 !< - param( 9): \f$ a_2 \f$ 99 !< - param(10): \f$ a_3 \f$ 100 !< - param(11): \f$ a_4 \f$ 101 !< - param(12): \f$ a_5 \f$ 102 !< - param(13): \f$ a_6 \f$ 103 !< - param(14): \f$ a_7 \f$ 104 !< - param(15): \f$ a_8 \f$ 105 !< - param(16): \f$ a_9 \f$ 106 !< - param(17): \f$ a_10 \f$ 107 !< - param(18): \f$ a_11 \f$ 108 double precision at1, at2, at3, at4, at5, at6, at7, at8, at9 109 double precision at10, at11, at0 110#endif 111 double precision tol_rho !< [Input] The lower limit on the density 112 integer nq !< [Input] The number of points 113 integer ipol !< [Input] The number of spin channels 114 double precision wght !< [Input] The weight of the functional 115c 116c Charge Density 117c 118 type(nwad_dble)::rho(nq,*) !< [Input] The density 119c 120c Charge Density Gradient Norm 121c 122 type(nwad_dble)::rgamma(nq,*) !< [Input] The density gradient norm 123c 124c Kinetic Energy Density 125c 126 type(nwad_dble)::tau(nq,*) !< [Input] The kinetic energy density 127c 128c Functional values 129c 130 type(nwad_dble)::func(*) !< [Output] The functional value 131c 132c Sampling Matrices for the XC Potential 133c 134c double precision Amat(nq,*) !< [Output] Derivative wrt density 135c double precision Cmat(nq,*) !< [Output] Derivative wrt rgamma 136c double precision Mmat(nq,*) !< [Output] Derivative wrt tau 137c 138 double precision pi 139c 140 integer n 141 double precision at, C1, C2, fL, fNL 142 type(nwad_dble)::rrho, rho43, rho13, rhoo, rho53 143c type(nwad_dble)::Gamma2, Gamma 144c type(nwad_dble)::TauUEG, Tsig, Wsig, W1, W2, W3, W4, W5, W6 145 type(nwad_dble)::Gamma2 146 type(nwad_dble)::TauUEG, Wsig, W1, W2, W3, W4, W5, W6 147 type(nwad_dble)::W7, W8, W9, W10, W11, Fsig 148c 149 type(nwad_dble)::tauN,tauu 150 double precision DTol 151 double precision F83, F23, F53, F43, F13, F1o3 152 double precision F1o4, F2o3, F3o2, F4o3, F4o9, F3o5 153 double precision One, Two, Three, Four, Five, Six, Seven, Eight 154 double precision Nine, F10, F11 155 type(nwad_dble)::s,s2,En,Ed,E 156 double precision Cs, Ax, P32, dE, dEn, dEd 157 158c functional derivatives below FFFFFFFFFFFF 159 160 double precision dFdW, dWdT, dTdR, dTdTau, dGGAdR, dFdR 161 double precision dFdTau, dGGAdG,tauW 162 163c functional derivatives above FFFFFFFFFFFF 164 165 166cedo parameter( pi = 3.1415926535897932384626433832795d0 ) 167 168 169 parameter (F1o3=1.d0/3.d0, F1o4=1.d0/4.d0, F2o3=2.d0/3.d0, 170 & F3o2=3.d0/2.d0,F13=1.d0/3.d0) 171 parameter (F4o3=4.d0/3.d0, F4o9=4.d0/9.d0, F3o5=3.d0/5.d0) 172 parameter (F83=8.d0/3.0d0, F23=2.0d0/3.d0, F53=5.d0/3.d0) 173 parameter (One=1.0d0, Two=2.0d0, Three=3.0d0, Four=4.0d0, 174 & Five=5.0d0,Six=6.0d0, Seven=7.0d0, 175 & Eight=8.0d0, Nine=9.0d0,F10=10.d0, F11=11.d0) 176 pi=acos(-1d0) 177 178 at0 = param(7) 179 at1 = param(8) 180 at2 = param(9) 181 at3 = param(10) 182 at4 = param(11) 183 at5 = param(12) 184 at6 = param(13) 185 at7 = param(14) 186 at8 = param(15) 187 at9 = param(16) 188 at10 = param(17) 189 at11 = param(18) 190c if (ijzy.eq.1) then 191c at0= 3.987756D-01 192c at1= 2.548219D-01 193c at2= 3.923994D-01 194c at3= -2.103655D+00 195c at4= -6.302147D+00 196c at5= 1.097615D+01 197c at6= 3.097273D+01 198c at7= -2.318489D+01 199c at8= -5.673480D+01 200c at9= 2.160364D+01 201c at10= 3.421814D+01 202c at11= -9.049762D+00 203c elseif (ijzy.eq.2) then 204c Parameters for M06-HF 205c at0= 1.179732D-01 206c at1= -1.066708D+00 207c at2= -1.462405D-01 208c at3= 7.481848D+00 209c at4= 3.776679D+00 210c at5= -4.436118D+01 211c at6= -1.830962D+01 212c at7= 1.003903D+02 213c at8= 3.864360D+01 214c at9= -9.806018D+01 215c at10= -2.557716D+01 216c at11= 3.590404D+01 217c elseif (ijzy.eq.3) then 218c Parameters for M06 219c at0= 5.877943D-01 220c at1= -1.371776D-01 221c at2= 2.682367D-01 222c at3= -2.515898D+00 223c at4= -2.978892D+00 224c at5= 8.710679D+00 225c at6= 1.688195D+01 226c at7= -4.489724D+00 227c at8= -3.299983D+01 228c at9= -1.449050D+01 229c at10= 2.043747D+01 230c at11= 1.256504D+01 231c elseif (ijzy.eq.4) then 232c Parameters for M06-2X 233c at0= 4.600000D-01 234c at1= -2.206052D-01 235c at2= -9.431788D-02 236c at3= 2.164494D+00 237c at4= -2.556466D+00 238c at5= -1.422133D+01 239c at6= 1.555044D+01 240c at7= 3.598078D+01 241c at8= -2.722754D+01 242c at9= -3.924093D+01 243c at10= 1.522808D+01 244c at11= 1.522227D+01 245c endif 246 247 at=1.0d0 248#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 249#if defined(NWAD_PRINT) 250 call nwxc_x_vs98_p(param,tol_rho, ipol, nq, wght, rho, rgamma, 251 & tau, func) 252#else 253 call nwxc_x_vs98(param,tol_rho, ipol, nq, wght, rho, rgamma, tau, 254 & func) 255#endif 256#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV) 257 call nwxc_x_vs98_d2(param,tol_rho, ipol, nq, wght, rho, rgamma, 258 & tau, func) 259#else 260 call nwxc_x_vs98_d3(param,tol_rho, ipol, nq, wght, rho, rgamma, 261 & tau, func) 262#endif 263 264 265 C1 = 0.2195149727645171d0 266 C2 = C1/0.804d0 267cedo DTol=1.0D-8 268 DTol=tol_rho 269C 270C Scale factors for local and non-local contributions. 271C 272 fL = wght 273 fNL = wght 274 Cs = 0.5d0/(3.0d0*pi*pi)**F13 275 P32 = (3.d0*pi**2)**F23 276 277c 278 Ax = (-0.75d0)*(3.0d0/pi)**F13 279 280 281c 282 if (ipol.eq.1 )then 283c 284c ======> SPIN-RESTRICTED <====== 285c or 286c SPIN-UNPOLARIZED 287c 288c 289 do 10 n = 1, nq 290 if (rho(n,R_T).lt.DTol) goto 10 291 292 rhoo = rho(n,R_T) 293 rho43 = rhoo**F4o3 294 rrho = 1d0/rhoo ! reciprocal of rho 295 rho13 = rho43*rrho 296 rho53 = rhoo**F53 297 298c 299 300 tauN = tau(n,T_T) 301 tauu=tauN 302 TauUEG=0.3d0*P32*rho53 303c Tsig =TauUEG/tauu 304c Wsig =(Tsig-One)/(Tsig+One) 305 Wsig =(TauUEG-tauu)/(TauUEG+tauu) 306 W1=Wsig 307 W2=Wsig*W1 308 W3=Wsig*W2 309 W4=Wsig*W3 310 W5=Wsig*W4 311 W6=Wsig*W5 312 W7=Wsig*W6 313 W8=Wsig*W7 314 W9=Wsig*W8 315 W10=Wsig*W9 316 W11=Wsig*W10 317 Fsig =at*(at0 + at1*W1+ at2*W2 + at3*W3 318 & + at4*W4 + at5*W5 + at6*W6 + at7*W7 319 & + at8*W8 + at9*W9 + at10*W10 + at11*W11) 320 321c Gamma2 = delrho(n,1,1)*delrho(n,1,1) + 322c & delrho(n,2,1)*delrho(n,2,1) + 323c & delrho(n,3,1)*delrho(n,3,1) 324 Gamma2 = rgamma(n,G_TT) 325c Gamma = sqrt(Gamma2) 326c s = Cs*Gamma/rho43 327c s2 = s*s 328 s2 = Cs*Cs*Gamma2/(rho43*rho43) 329 En = C1*s2 330 Ed = One + C2*s2 331 E = En/Ed 332c Ex = Ex + rho43*Ax*(fL+fNL*E)*Fsig*qwght(n) 333 func(n)=func(n)+rho43*Ax*(fL+fNL*E)*Fsig 334c 335c functional derivatives 336c 337c dEn = Two*C1*s 338c dEd = Two*C2*s 339c dE = (dEn*Ed-En*dEd)/(Ed*Ed) 340c dFdW = at*(at1 + Two*at2*W1 + Three*at3*W2 341c & + Four*at4*W3 + Five*at5*W4 342c & + Six*at6*W5 + Seven*at7*W6 343c & + Eight*at8*W7 + Nine*at9*W8 344c & + F10*at10*W9 + F11*at11*W10) 345c dWdT = Two/((One + Tsig)**2) 346c dTdR = (0.5d0*P32*rho13*rho13)/tauu 347c dTdTau = -TauUEG/tauu**2 348c dGGAdR = F4o3*rho13*Ax*(fL+fNL*(E-s*dE)) 349c dFdR = dFdW*dWdT*dTdR 350c dFdTau=dFdW*dWdT*dTdTau 351c dGGAdG =(fNL*dE*s/(Two*Gamma2)) 352c Amat(n,D1_RA) = Amat(n,D1_RA) + dGGAdR*Fsig 353c & + (fL+fNL*E)*Ax*rho43*dFdR 354c Cmat(n,D1_GAA) = Cmat(n,D1_GAA) 355c & + Two*dGGAdG*Ax*rho43*Fsig 356c Mmat(n,D1_TA) = Mmat(n,D1_TA) 357c & + 0.5d0*rho43*Ax*(fL+fNL*E)*dFdTau 358 35910 continue 360 361 362c UUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUnrestricted 363 else ! ipol=2 364c 365c ======> SPIN-UNRESTRICTED <====== 366 367c 368c use spin density functional theory ie n-->2n 369c Ex=(1/2)Ex[2*alpha] + (1/2)Ex[2*beta] 370 371 do 20 n = 1, nq 372 if (rho(n,R_A)+rho(n,R_B).lt.DTol) goto 20 373c 374c Alpha ALPHA ALPHA 375c 376 if (rho(n,R_A).lt.0.5d0*DTol) goto 25 377 rhoo = Two*rho(n,R_A) 378 rho43 = rhoo**F4o3 379 rrho = 1.0d0/rhoo ! reciprocal of rho 380 rho13 = rho43*rrho 381 rho53 = rhoo**F53 382c 383 tauN = tau(n,T_A) 384 tauu = Two*tauN 385 TauUEG=0.3d0*P32*rho53 386c Tsig =TauUEG/tauu 387c Wsig =(Tsig-One)/(Tsig+One) 388 Wsig =(TauUEG-tauu)/(TauUEG+tauu) 389 W1=Wsig 390 W2=Wsig*W1 391 W3=Wsig*W2 392 W4=Wsig*W3 393 W5=Wsig*W4 394 W6=Wsig*W5 395 W7=Wsig*W6 396 W8=Wsig*W7 397 W9=Wsig*W8 398 W10=Wsig*W9 399 W11=Wsig*W10 400 Fsig =at*(at0 + at1*W1+ at2*W2 + at3*W3 401 & + at4*W4 + at5*W5 + at6*W6 + at7*W7 402 & + at8*W8 + at9*W9 + at10*W10 + at11*W11) 403 404 405c Gamma2 = delrho(n,1,1)*delrho(n,1,1) + 406c & delrho(n,2,1)*delrho(n,2,1) + 407c & delrho(n,3,1)*delrho(n,3,1) 408 Gamma2 = rgamma(n,G_AA) 409 Gamma2 = Four*Gamma2 410c Gamma = sqrt(Gamma2) 411c s = Cs*Gamma/rho43 412c s2 = s*s 413 s2 = Cs*Cs*Gamma2/(rho43*rho43) 414 En = C1*s2 415 Ed = One + C2*s2 416 E = En/Ed 417c Ex = Ex + rho43*Ax*(fL+fNL*E)*(Fsig)*0.5d0*qwght(n) 418 func(n)=func(n)+rho43*Ax*(fL+fNL*E)*(Fsig)*.5d0 419c 420c functional derivatives 421c 422c dEn = Two*C1*s 423c dEd = Two*C2*s 424c dE = (dEn*Ed-En*dEd)/(Ed*Ed) 425c dFdW = at*(at1 + Two*at2*W1 + Three*at3*W2 426c & + Four*at4*W3 + Five*at5*W4 427c & + Six*at6*W5 + Seven*at7*W6 428c & + Eight*at8*W7 + Nine*at9*W8 429c & + F10*at10*W9 + F11*at11*W10) 430c dWdT = Two/((One + Tsig)**2) 431c dTdR = Two*(0.5d0*P32*rho13*rho13)/tauu 432c dTdTau = -Two*TauUEG/tauu**2 433c dGGAdR = Two*F4o3*rho13*Ax*(fL+fNL*(E-s*dE)) 434c dFdR = dFdW*dWdT*dTdR 435c dFdTau=dFdW*dWdT*dTdTau 436c dGGAdG =Four*(fNL*dE*s/(Two*Gamma2)) 437 438c Amat(n,D1_RA) = Amat(n,D1_RA) + (dGGAdR*(Fsig) 439c & + (fL+fNL*E)*Ax*rho43*dFdR)*0.5d0 440c Cmat(n,D1_GAA) = Cmat(n,D1_GAA) 441c & + dGGAdG*Ax*rho43*(Fsig)*0.5d0 442c Mmat(n,D1_TA) = Mmat(n,D1_TA) 443c & + 0.5d0*0.5d0*rho43*Ax*(fL+fNL*E)*dFdTau 444 445c write (*,*) "Ex,Amat(n,1),Cmat(n,1)", 446c & Ex,Amat(n,1),Cmat(n,1) 447 448c 449c Beta BETA BETA 450c 451 45225 continue 453 454c 455c Beta 456c 457 if (rho(n,R_B).lt.0.5d0*DTol) goto 20 458 rhoo = Two*rho(n,R_B) 459 rho43 = rhoo**F4o3 460 rrho = 1.0d0/rhoo ! reciprocal of rho 461 rho13 = rho43*rrho 462 rho53 = rhoo**F53 463c 464 tauN = tau(n,T_B) 465 tauu = Two*tauN 466 TauUEG=0.3d0*P32*rho53 467c Tsig =TauUEG/tauu 468c Wsig =(Tsig-One)/(Tsig+One) 469 Wsig = (TauUEG-tauu)/(TauUEG+tauu) 470 W1=Wsig 471 W2=Wsig*W1 472 W3=Wsig*W2 473 W4=Wsig*W3 474 W5=Wsig*W4 475 W6=Wsig*W5 476 W7=Wsig*W6 477 W8=Wsig*W7 478 W9=Wsig*W8 479 W10=Wsig*W9 480 W11=Wsig*W10 481 Fsig =at*(at0+at1*W1+ at2*W2 + at3*W3 482 & + at4*W4 + at5*W5 + at6*W6 + at7*W7 483 & + at8*W8 + at9*W9 + at10*W10 + at11*W11) 484 485 486c Gamma2 = delrho(n,1,2)*delrho(n,1,2) + 487c & delrho(n,2,2)*delrho(n,2,2) + 488c & delrho(n,3,2)*delrho(n,3,2) 489 Gamma2 = rgamma(n,G_BB) 490 Gamma2 = Four*Gamma2 491c Gamma = sqrt(Gamma2) 492c s = Cs*Gamma/rho43 493c s2 = s*s 494 s2 = Cs*Cs*Gamma2/(rho43*rho43) 495 En = C1*s2 496 Ed = One + C2*s2 497 E = En/Ed 498c Ex = Ex + rho43*Ax*(fL+fNL*E)*(Fsig)*0.5d0*qwght(n) 499 func(n)= func(n)+rho43*Ax*(fL+fNL*E)*(Fsig)*.5d0 500c 501c functional derivatives 502c 503c dEn = Two*C1*s 504c dEd = Two*C2*s 505c dE = (dEn*Ed-En*dEd)/(Ed*Ed) 506c dFdW = at*(at1 + Two*at2*W1 + Three*at3*W2 507c & + Four*at4*W3 + Five*at5*W4 508c & + Six*at6*W5 + Seven*at7*W6 509c & + Eight*at8*W7 + Nine*at9*W8 510c & + F10*at10*W9 + F11*at11*W10) 511c dWdT = Two/((One + Tsig)**2) 512c dTdR = Two*(0.5d0*P32*rho13*rho13)/tauu 513c dTdTau = -Two*TauUEG/tauu**2 514c dGGAdR = Two*F4o3*rho13*Ax*(fL+fNL*(E-s*dE)) 515c dFdR = dFdW*dWdT*dTdR 516c dFdTau=dFdW*dWdT*dTdTau 517c dGGAdG =Four*(fNL*dE*s/(Two*Gamma2)) 518 519c Amat(n,D1_RB) = Amat(n,D1_RB) + (dGGAdR*(Fsig) 520c & + (fL+fNL*E)*Ax*rho43*dFdR)*0.5d0 521c Cmat(n,D1_GBB) = Cmat(n,D1_GBB) 522c & + dGGAdG*Ax*rho43*(Fsig)*0.5d0 523c Mmat(n,D1_TB) = Mmat(n,D1_TB) 524c & + 0.5d0*0.5d0*rho43*Ax*(fL+fNL*E)*dFdTau 525c 52620 continue 527 endif 528c 529 return 530 end 531C> 532C> \brief Evaluate the M06-2X exchange functional 533C> 534C> This routine evaluates the M06-2X exchange functional [1,2]. This functional 535C> is closely related to the M06, M06-HF and M06-L exchange functionals 536C> except that it does not use the VS98 exchange functional, hence the 537C> parameter list is defined differently. 538C> 539C> ### References ### 540C> 541C> [1] Y Zhao, DG Truhlar, 542C> "A new local density functional for main-group thermochemistry, 543C> transition metal bonding, thermochemical kinetics, and noncovalent 544C> interactions", 545C> J. Chem. Phys. <b>125</b>, 194101 (2006), DOI: 546C> <a href="https://doi.org/10.1063/1.2370993"> 547C> 10.1063/1.2370993</a>. 548C> 549C> [2] Y Zhao, DG Truhlar, 550C> "Density functional for spectroscopy: No long-range self-interaction 551C> error, good performance for Rydberg and charge-transfer states, 552C> and better performance on average than B3LYP for ground states", 553C> J. Phys. Chem. A <b>110</b>, 13126-13130 (2006), DOI: 554C> <a href="https://doi.org/10.1021/jp066479k"> 555C> 10.1021/jp066479k</a>. 556C> 557#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 558#if defined(NWAD_PRINT) 559 Subroutine nwxc_x_m06_2x_p(param, tol_rho, ipol, nq, wght, 560 & rho, rgamma, tau, func) 561#else 562 Subroutine nwxc_x_m06_2x(param, tol_rho, ipol, nq, wght, 563 & rho, rgamma, tau, func) 564#endif 565#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV) 566 Subroutine nwxc_x_m06_2x_d2(param, tol_rho, ipol, nq, wght, 567 & rho, rgamma, tau, func) 568#else 569 Subroutine nwxc_x_m06_2x_d3(param, tol_rho, ipol, nq, wght, 570 & rho, rgamma, tau, func) 571#endif 572c 573c$Id$ 574c 575#include "nwad.fh" 576c 577 implicit none 578c 579#include "nwxc_param.fh" 580c 581#if defined(NWAD_PRINT) 582#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 583 type(nwad_dble)::param(*) !< [Input] Parameters of the functional 584 type(nwad_dble)::at1, at2, at3, at4, at5, at6, at7, at8, at9 585 type(nwad_dble)::at10, at11, at0 586#else 587 double precision param(*) !< [Input] Parameters of the functional 588 double precision at1, at2, at3, at4, at5, at6, at7, at8, at9 589 double precision at10, at11, at0 590#endif 591#else 592 double precision param(*) !< [Input] Parameters of the functional 593 !< - param( 1): \f$ a_0 \f$ 594 !< - param( 2): \f$ a_1 \f$ 595 !< - param( 3): \f$ a_2 \f$ 596 !< - param( 4): \f$ a_3 \f$ 597 !< - param( 5): \f$ a_4 \f$ 598 !< - param( 6): \f$ a_5 \f$ 599 !< - param( 7): \f$ a_6 \f$ 600 !< - param( 8): \f$ a_7 \f$ 601 !< - param( 9): \f$ a_8 \f$ 602 !< - param(10): \f$ a_9 \f$ 603 !< - param(11): \f$ a_10 \f$ 604 !< - param(12): \f$ a_11 \f$ 605 double precision at1, at2, at3, at4, at5, at6, at7, at8, at9 606 double precision at10, at11, at0 607#endif 608 double precision tol_rho !< [Input] The lower limit on the density 609 integer nq !< [Input] The number of points 610 integer ipol !< [Input] The number of spin channels 611 double precision wght !< [Input] The weight of the functional 612c 613c Charge Density 614c 615 type(nwad_dble)::rho(nq,*) !< [Input] The density 616c 617c Charge Density Gradient Norm 618c 619 type(nwad_dble)::rgamma(nq,*) !< [Input] The density gradient norm 620c 621c Kinetic Energy Density 622c 623 type(nwad_dble)::tau(nq,*) !< [Input] The kinetic energy density 624c 625c Functional values 626c 627 type(nwad_dble)::func(*) !< [Output] The functional value 628c 629c Sampling Matrices for the XC Potential 630c 631c double precision Amat(nq,*) !< [Output] Derivative wrt density 632c double precision Cmat(nq,*) !< [Output] Derivative wrt rgamma 633c double precision Mmat(nq,*) !< [Output] Derivative wrt tau 634c 635 double precision pi 636c 637 integer n 638 double precision at, C1, C2, fL, fNL 639 type(nwad_dble)::rrho, rho43, rho13, rhoo, rho53 640c type(nwad_dble)::Gamma2, Gamma 641c type(nwad_dble)::TauUEG, Tsig, Wsig, W1, W2, W3, W4, W5, W6 642 type(nwad_dble)::Gamma2 643 type(nwad_dble)::TauUEG, Wsig, W1, W2, W3, W4, W5, W6 644 type(nwad_dble)::W7, W8, W9, W10, W11, Fsig 645c 646 type(nwad_dble)::tauN,tauu 647 double precision DTol 648 double precision F83, F23, F53, F43, F13, F1o3 649 double precision F1o4, F2o3, F3o2, F4o3, F4o9, F3o5 650 double precision One, Two, Three, Four, Five, Six, Seven, Eight 651 double precision Nine, F10, F11 652 type(nwad_dble)::s,s2,En,Ed,E 653 double precision Cs, Ax, P32, dE, dEn, dEd 654 655c functional derivatives below FFFFFFFFFFFF 656 657 double precision dFdW, dWdT, dTdR, dTdTau, dGGAdR, dFdR 658 double precision dFdTau, dGGAdG,tauW 659 660c functional derivatives above FFFFFFFFFFFF 661 662 663cedo parameter( pi = 3.1415926535897932384626433832795d0 ) 664 665 666 parameter (F1o3=1.d0/3.d0, F1o4=1.d0/4.d0, F2o3=2.d0/3.d0, 667 & F3o2=3.d0/2.d0,F13=1.d0/3.d0) 668 parameter (F4o3=4.d0/3.d0, F4o9=4.d0/9.d0, F3o5=3.d0/5.d0) 669 parameter (F83=8.d0/3.0d0, F23=2.0d0/3.d0, F53=5.d0/3.d0) 670 parameter (One=1.0d0, Two=2.0d0, Three=3.0d0, Four=4.0d0, 671 & Five=5.0d0,Six=6.0d0, Seven=7.0d0, 672 & Eight=8.0d0, Nine=9.0d0,F10=10.d0, F11=11.d0) 673 pi=acos(-1d0) 674 675 at0 = param(1) 676 at1 = param(2) 677 at2 = param(3) 678 at3 = param(4) 679 at4 = param(5) 680 at5 = param(6) 681 at6 = param(7) 682 at7 = param(8) 683 at8 = param(9) 684 at9 = param(10) 685 at10 = param(11) 686 at11 = param(12) 687c if (ijzy.eq.1) then 688c at0= 3.987756D-01 689c at1= 2.548219D-01 690c at2= 3.923994D-01 691c at3= -2.103655D+00 692c at4= -6.302147D+00 693c at5= 1.097615D+01 694c at6= 3.097273D+01 695c at7= -2.318489D+01 696c at8= -5.673480D+01 697c at9= 2.160364D+01 698c at10= 3.421814D+01 699c at11= -9.049762D+00 700c elseif (ijzy.eq.2) then 701c Parameters for M06-HF 702c at0= 1.179732D-01 703c at1= -1.066708D+00 704c at2= -1.462405D-01 705c at3= 7.481848D+00 706c at4= 3.776679D+00 707c at5= -4.436118D+01 708c at6= -1.830962D+01 709c at7= 1.003903D+02 710c at8= 3.864360D+01 711c at9= -9.806018D+01 712c at10= -2.557716D+01 713c at11= 3.590404D+01 714c elseif (ijzy.eq.3) then 715c Parameters for M06 716c at0= 5.877943D-01 717c at1= -1.371776D-01 718c at2= 2.682367D-01 719c at3= -2.515898D+00 720c at4= -2.978892D+00 721c at5= 8.710679D+00 722c at6= 1.688195D+01 723c at7= -4.489724D+00 724c at8= -3.299983D+01 725c at9= -1.449050D+01 726c at10= 2.043747D+01 727c at11= 1.256504D+01 728c elseif (ijzy.eq.4) then 729c Parameters for M06-2X 730c at0= 4.600000D-01 731c at1= -2.206052D-01 732c at2= -9.431788D-02 733c at3= 2.164494D+00 734c at4= -2.556466D+00 735c at5= -1.422133D+01 736c at6= 1.555044D+01 737c at7= 3.598078D+01 738c at8= -2.722754D+01 739c at9= -3.924093D+01 740c at10= 1.522808D+01 741c at11= 1.522227D+01 742c endif 743 744 at=1.0d0 745 746 C1 = 0.2195149727645171d0 747 C2 = C1/0.804d0 748cedo DTol=1.0D-8 749 DTol=tol_rho 750C 751C Scale factors for local and non-local contributions. 752C 753 fL = wght 754 fNL = wght 755 Cs = 0.5d0/(3.0d0*pi*pi)**F13 756 P32 = (3.d0*pi**2)**F23 757 758c 759 Ax = (-0.75d0)*(3.0d0/pi)**F13 760 761 762c 763 if (ipol.eq.1 )then 764c 765c ======> SPIN-RESTRICTED <====== 766c or 767c SPIN-UNPOLARIZED 768c 769c 770 do 10 n = 1, nq 771 if (rho(n,R_T).lt.DTol) goto 10 772 773 rhoo = rho(n,R_T) 774 rho43 = rhoo**F4o3 775 rrho = 1d0/rhoo ! reciprocal of rho 776 rho13 = rho43*rrho 777 rho53 = rhoo**F53 778c 779 tauN = tau(n,T_T) 780 tauu=tauN 781 TauUEG=0.3d0*P32*rho53 782c Tsig =TauUEG/tauu 783c Wsig =(Tsig-One)/(Tsig+One) 784 Wsig =(TauUEG-tauu)/(TauUEG+tauu) 785 W1=Wsig 786 W2=Wsig*W1 787 W3=Wsig*W2 788 W4=Wsig*W3 789 W5=Wsig*W4 790 W6=Wsig*W5 791 W7=Wsig*W6 792 W8=Wsig*W7 793 W9=Wsig*W8 794 W10=Wsig*W9 795 W11=Wsig*W10 796 Fsig =at*(at0 + at1*W1+ at2*W2 + at3*W3 797 & + at4*W4 + at5*W5 + at6*W6 + at7*W7 798 & + at8*W8 + at9*W9 + at10*W10 + at11*W11) 799 800c Gamma2 = delrho(n,1,1)*delrho(n,1,1) + 801c & delrho(n,2,1)*delrho(n,2,1) + 802c & delrho(n,3,1)*delrho(n,3,1) 803 Gamma2 = rgamma(n,G_TT) 804c Gamma = sqrt(Gamma2) 805c s = Cs*Gamma/rho43 806c s2 = s*s 807 s2 = Cs*Cs*Gamma2/(rho43*rho43) 808 En = C1*s2 809 Ed = One + C2*s2 810 E = En/Ed 811c Ex = Ex + rho43*Ax*(fL+fNL*E)*Fsig*qwght(n) 812 func(n)=func(n)+rho43*Ax*(fL+fNL*E)*Fsig 813c 814c functional derivatives 815c 816c dEn = Two*C1*s 817c dEd = Two*C2*s 818c dE = (dEn*Ed-En*dEd)/(Ed*Ed) 819c dFdW = at*(at1 + Two*at2*W1 + Three*at3*W2 820c & + Four*at4*W3 + Five*at5*W4 821c & + Six*at6*W5 + Seven*at7*W6 822c & + Eight*at8*W7 + Nine*at9*W8 823c & + F10*at10*W9 + F11*at11*W10) 824c dWdT = Two/((One + Tsig)**2) 825c dTdR = (0.5d0*P32*rho13*rho13)/tauu 826c dTdTau = -TauUEG/tauu**2 827c dGGAdR = F4o3*rho13*Ax*(fL+fNL*(E-s*dE)) 828c dFdR = dFdW*dWdT*dTdR 829c dFdTau=dFdW*dWdT*dTdTau 830c dGGAdG =(fNL*dE*s/(Two*Gamma2)) 831c Amat(n,D1_RA) = Amat(n,D1_RA) + dGGAdR*Fsig 832c & + (fL+fNL*E)*Ax*rho43*dFdR 833c Cmat(n,D1_GAA) = Cmat(n,D1_GAA) 834c & + Two*dGGAdG*Ax*rho43*Fsig 835c Mmat(n,D1_TA) = Mmat(n,D1_TA) 836c & + 0.5d0*rho43*Ax*(fL+fNL*E)*dFdTau 837 83810 continue 839 840 841c UUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUnrestricted 842 else ! ipol=2 843c 844c ======> SPIN-UNRESTRICTED <====== 845 846c 847c use spin density functional theory ie n-->2n 848c Ex=(1/2)Ex[2*alpha] + (1/2)Ex[2*beta] 849 850 do 20 n = 1, nq 851 if (rho(n,R_A)+rho(n,R_B).lt.DTol) goto 20 852c 853c Alpha ALPHA ALPHA 854c 855 if (rho(n,R_A).lt.0.5d0*DTol) goto 25 856 rhoo = Two*rho(n,R_A) 857 rho43 = rhoo**F4o3 858 rrho = 1.0d0/rhoo ! reciprocal of rho 859 rho13 = rho43*rrho 860 rho53 = rhoo**F53 861c 862 tauN = tau(n,T_A) 863 tauu = Two*tauN 864 TauUEG=0.3d0*P32*rho53 865c Tsig =TauUEG/tauu 866c Wsig =(Tsig-One)/(Tsig+One) 867 Wsig =(TauUEG-tauu)/(TauUEG+tauu) 868 W1=Wsig 869 W2=Wsig*W1 870 W3=Wsig*W2 871 W4=Wsig*W3 872 W5=Wsig*W4 873 W6=Wsig*W5 874 W7=Wsig*W6 875 W8=Wsig*W7 876 W9=Wsig*W8 877 W10=Wsig*W9 878 W11=Wsig*W10 879 Fsig =at*(at0 + at1*W1+ at2*W2 + at3*W3 880 & + at4*W4 + at5*W5 + at6*W6 + at7*W7 881 & + at8*W8 + at9*W9 + at10*W10 + at11*W11) 882 883 884c Gamma2 = delrho(n,1,1)*delrho(n,1,1) + 885c & delrho(n,2,1)*delrho(n,2,1) + 886c & delrho(n,3,1)*delrho(n,3,1) 887 Gamma2 = rgamma(n,G_AA) 888 Gamma2 = Four*Gamma2 889c Gamma = sqrt(Gamma2) 890c s = Cs*Gamma/rho43 891c s2 = s*s 892 s2 = Cs*Cs*Gamma2/(rho43*rho43) 893 En = C1*s2 894 Ed = One + C2*s2 895 E = En/Ed 896c Ex = Ex + rho43*Ax*(fL+fNL*E)*(Fsig)*0.5d0*qwght(n) 897 func(n)=func(n)+rho43*Ax*(fL+fNL*E)*(Fsig)*.5d0 898c 899c functional derivatives 900c 901c dEn = Two*C1*s 902c dEd = Two*C2*s 903c dE = (dEn*Ed-En*dEd)/(Ed*Ed) 904c dFdW = at*(at1 + Two*at2*W1 + Three*at3*W2 905c & + Four*at4*W3 + Five*at5*W4 906c & + Six*at6*W5 + Seven*at7*W6 907c & + Eight*at8*W7 + Nine*at9*W8 908c & + F10*at10*W9 + F11*at11*W10) 909c dWdT = Two/((One + Tsig)**2) 910c dTdR = Two*(0.5d0*P32*rho13*rho13)/tauu 911c dTdTau = -Two*TauUEG/tauu**2 912c dGGAdR = Two*F4o3*rho13*Ax*(fL+fNL*(E-s*dE)) 913c dFdR = dFdW*dWdT*dTdR 914c dFdTau=dFdW*dWdT*dTdTau 915c dGGAdG =Four*(fNL*dE*s/(Two*Gamma2)) 916 917c Amat(n,D1_RA) = Amat(n,D1_RA) + (dGGAdR*(Fsig) 918c & + (fL+fNL*E)*Ax*rho43*dFdR)*0.5d0 919c Cmat(n,D1_GAA) = Cmat(n,D1_GAA) 920c & + dGGAdG*Ax*rho43*(Fsig)*0.5d0 921c Mmat(n,D1_TA) = Mmat(n,D1_TA) 922c & + 0.5d0*0.5d0*rho43*Ax*(fL+fNL*E)*dFdTau 923 924c write (*,*) "Ex,Amat(n,1),Cmat(n,1)", 925c & Ex,Amat(n,1),Cmat(n,1) 926 927c 928c Beta BETA BETA 929c 930 93125 continue 932 933c 934c Beta 935c 936 if (rho(n,R_B).lt.0.5d0*DTol) goto 20 937 rhoo = Two*rho(n,R_B) 938 rho43 = rhoo**F4o3 939 rrho = 1.0d0/rhoo ! reciprocal of rho 940 rho13 = rho43*rrho 941 rho53 = rhoo**F53 942c 943 tauN = tau(n,T_B) 944 tauu = Two*tauN 945 TauUEG=0.3d0*P32*rho53 946c Tsig =TauUEG/tauu 947c Wsig =(Tsig-One)/(Tsig+One) 948 Wsig =(TauUEG-tauu)/(TauUEG+tauu) 949 W1=Wsig 950 W2=Wsig*W1 951 W3=Wsig*W2 952 W4=Wsig*W3 953 W5=Wsig*W4 954 W6=Wsig*W5 955 W7=Wsig*W6 956 W8=Wsig*W7 957 W9=Wsig*W8 958 W10=Wsig*W9 959 W11=Wsig*W10 960 Fsig =at*(at0+at1*W1+ at2*W2 + at3*W3 961 & + at4*W4 + at5*W5 + at6*W6 + at7*W7 962 & + at8*W8 + at9*W9 + at10*W10 + at11*W11) 963 964 965c Gamma2 = delrho(n,1,2)*delrho(n,1,2) + 966c & delrho(n,2,2)*delrho(n,2,2) + 967c & delrho(n,3,2)*delrho(n,3,2) 968 Gamma2 = rgamma(n,G_BB) 969 Gamma2 = Four*Gamma2 970c Gamma = sqrt(Gamma2) 971c s = Cs*Gamma/rho43 972c s2 = s*s 973 s2 = Cs*Cs*Gamma2/(rho43*rho43) 974 En = C1*s2 975 Ed = One + C2*s2 976 E = En/Ed 977c Ex = Ex + rho43*Ax*(fL+fNL*E)*(Fsig)*0.5d0*qwght(n) 978 func(n)= func(n)+rho43*Ax*(fL+fNL*E)*(Fsig)*.5d0 979c 980c functional derivatives 981c 982c dEn = Two*C1*s 983c dEd = Two*C2*s 984c dE = (dEn*Ed-En*dEd)/(Ed*Ed) 985c dFdW = at*(at1 + Two*at2*W1 + Three*at3*W2 986c & + Four*at4*W3 + Five*at5*W4 987c & + Six*at6*W5 + Seven*at7*W6 988c & + Eight*at8*W7 + Nine*at9*W8 989c & + F10*at10*W9 + F11*at11*W10) 990c dWdT = Two/((One + Tsig)**2) 991c dTdR = Two*(0.5d0*P32*rho13*rho13)/tauu 992c dTdTau = -Two*TauUEG/tauu**2 993c dGGAdR = Two*F4o3*rho13*Ax*(fL+fNL*(E-s*dE)) 994c dFdR = dFdW*dWdT*dTdR 995c dFdTau=dFdW*dWdT*dTdTau 996c dGGAdG =Four*(fNL*dE*s/(Two*Gamma2)) 997 998c Amat(n,D1_RB) = Amat(n,D1_RB) + (dGGAdR*(Fsig) 999c & + (fL+fNL*E)*Ax*rho43*dFdR)*0.5d0 1000c Cmat(n,D1_GBB) = Cmat(n,D1_GBB) 1001c & + dGGAdG*Ax*rho43*(Fsig)*0.5d0 1002c Mmat(n,D1_TB) = Mmat(n,D1_TB) 1003c & + 0.5d0*0.5d0*rho43*Ax*(fL+fNL*E)*dFdTau 1004c 100520 continue 1006 endif 1007c 1008 return 1009 end 1010#if !defined(NWAD_PRINT) 1011#define NWAD_PRINT 1012c 1013c Compile source again for Maxima 1014c 1015#include "nwxc_x_m06.F" 1016#endif 1017#if !defined(SECOND_DERIV) 1018#define SECOND_DERIV 1019c 1020c Compile source again for the 2nd derivative case 1021c 1022#include "nwxc_x_m06.F" 1023#endif 1024#if !defined(THIRD_DERIV) 1025#define THIRD_DERIV 1026c 1027c Compile source again for the 3rd derivative case 1028c 1029#include "nwxc_x_m06.F" 1030#endif 1031#undef NWAD_PRINT 1032C> @} 1033