1#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 2#if !defined(NWAD_PRINT) 3C> \ingroup nwxc 4C> @{ 5C> 6C> \file nwxc_c_m06.F 7C> Implementation of the M06 correlation functional 8C> 9C> @} 10#endif 11#endif 12C> 13C> \ingroup nwxc_priv 14C> @{ 15C> 16C> \brief The M06 correlation functional 17C> 18C> The M06 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, DG Truhlar, 28C> "A new local density functional for main-group thermochemistry, 29C> transition metal bonding, thermochemical kinetics, and noncovalent 30C> interactions", 31C> J. Chem. Phys. <b>125</b>, 194101 (2006), DOI: 32C> <a href="https://doi.org/10.1063/1.2370993"> 33C> 10.1063/1.2370993</a>. 34C> 35C> [2] Y Zhao, DG Truhlar, 36C> "Density functional for spectroscopy: No long-range self-interaction 37C> error, good performance for Rydberg and charge-transfer states, 38C> and better performance on average than B3LYP for ground states", 39C> J. Phys. Chem. A <b>110</b>, 13126-13130 (2006), DOI: 40C> <a href="https://doi.org/10.1021/jp066479k"> 41C> 10.1021/jp066479k</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 M06 suite correlation functional 51c META GGA 52C utilizes ingredients: 53c rho - density 54c delrho - gradient of density 55c tau (tauN)- K.S kinetic energy density 56c ijzy - 1 M06-L 57c ijzy - 2 M06-HF 58c ijzy - 3 M06 59c ijzy - 4 M06-2X 60 61#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 62#if defined(NWAD_PRINT) 63 Subroutine nwxc_c_m06_p(param, tol_rho, ipol, nq, wght, rho, 64 & rgamma, tau, func) 65#else 66 Subroutine nwxc_c_m06(param, tol_rho, ipol, nq, wght, rho, rgamma, 67 & tau, func) 68#endif 69#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV) 70 Subroutine nwxc_c_m06_d2(param, tol_rho, ipol, nq, wght, rho, 71 & rgamma, tau, func) 72#else 73 Subroutine nwxc_c_m06_d3(param, tol_rho, ipol, nq, wght, rho, 74 & rgamma, tau, func) 75#endif 76c 77c$Id$ 78c 79c [a] Zhao, Y. and Truhlar, D. G. J. Chem. Phys. 2006, 125, 194101; 80c [b] Zhao, Y. and Truhlar, D. G. J. Phys. Chem. A (2006),110(49),13126-13130. 81c 82#include "nwad.fh" 83c 84 implicit none 85c 86#include "intf_nwxc_c_lsda.fh" 87#include "intf_nwxc_c_vs98.fh" 88#include "intf_nwxc_m06css.fh" 89c 90#include "nwxc_param.fh" 91c 92c Input and other parameters 93c 94#if defined(NWAD_PRINT) 95#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 96 type(nwad_dble)::param(*) 97 type(nwad_dble):: sopp0, sopp1,sopp2, sopp3, sopp4 98#else 99 double precision param(*) 100 double precision sopp0, sopp1,sopp2, sopp3, sopp4 101#endif 102#else 103 double precision param(*) !< [Input] Parameters of functional 104 !< - param(1): \f$ d_{C\alpha\beta,0} \f$ 105 !< - param(2): \f$ d_{C\alpha\beta,1} \f$ 106 !< - param(3): \f$ d_{C\alpha\beta,2} \f$ 107 !< - param(4): \f$ d_{C\alpha\beta,3} \f$ 108 !< - param(5): \f$ d_{C\alpha\beta,4} \f$ 109 !< - param(6): \f$ d_{C\alpha\beta,5} \f$ 110 !< - param(7): \f$ d_{C\sigma\sigma,0} \f$ 111 !< - param(8): \f$ d_{C\sigma\sigma,1} \f$ 112 !< - param(9): \f$ d_{C\sigma\sigma,2} \f$ 113 !< - param(10): \f$ d_{C\sigma\sigma,3} \f$ 114 !< - param(11): \f$ d_{C\sigma\sigma,4} \f$ 115 !< - param(12): \f$ d_{C\sigma\sigma,5} \f$ 116 !< - param(13): \f$ c_{C\alpha\beta,0} \f$ 117 !< - param(14): \f$ c_{C\alpha\beta,1} \f$ 118 !< - param(15): \f$ c_{C\alpha\beta,2} \f$ 119 !< - param(16): \f$ c_{C\alpha\beta,3} \f$ 120 !< - param(17): \f$ c_{C\alpha\beta,4} \f$ 121 !< - param(18): \f$ c_{C\sigma\sigma,0} \f$ 122 !< - param(19): \f$ c_{C\sigma\sigma,1} \f$ 123 !< - param(20): \f$ c_{C\sigma\sigma,2} \f$ 124 !< - param(21): \f$ c_{C\sigma\sigma,3} \f$ 125 !< - param(22): \f$ c_{C\sigma\sigma,4} \f$ 126 double precision sopp0, sopp1,sopp2, sopp3, sopp4 127#endif 128 double precision tol_rho !< [Input] The lower limit on the density 129 integer nq !< [Input] The number of points 130 integer ipol !< [Input] The number of spin channels 131 double precision wght !< [Input] The weight of the functional 132c 133c Charge Density 134c 135 type(nwad_dble)::rho(nq,*) !< [Input] The density 136c 137c Charge Density Gradient Norm 138c 139 type(nwad_dble)::rgamma(nq,*) !< [Input] The density gradient norm 140c 141c Kinetic Energy Density 142c 143 type(nwad_dble)::tau(nq,*) !< [Input] The kinetic energy density 144c 145c Functional values 146c 147 type(nwad_dble)::func(*) !< [Output] The functional value 148c 149c Sampling Matrices for the XC Potential 150c 151c double precision Amat(nq,*) !< [Output] Derivative wrt density 152c double precision Cmat(nq,*) !< [Output] Derivative wrt rgamma 153c double precision Mmat(nq,*) !< [Output] Derivative wrt tau 154c 155c Threshold parameters 156c 157 double precision F1, F2, F3, F4,COpp 158 Data COpp/0.0031d0/,F1/1.0d0/,F2/2.0d0/, 159 & F3/3.0d0/,F4/4.0d0/ 160 161 integer n 162 163c call to the m06css subroutine 164 type(nwad_dble)::PA,GAA,TA,FA,EUA,ChiA 165 double precision FPA,FGA,FTA,EUEGA,EUPA,ChiAP,ChiAG 166 type(nwad_dble)::PB,GBB,TB,FB,EUB,ChiB 167 double precision FPB,FGB,FTB,EUEGB,EUPB,ChiBP,ChiBG 168c 169 type(nwad_dble)::RS,Zeta,PotLC,P,U,W,EUEG 170 double precision sop 171 double precision Pi, F6, F43, Pi34, F13, 172 &RSP,dZdA,dZdB,dLdS,dLdZ,d2LdSS,d2LdSZ,d2LdZZ 173 double precision dUdChiA,dUdChiB,dUdPA,dUdPB,dUdGA,dUdGB, 174 &dWdU,dWdPA,dWdPB, dWdGA,dWdGB,EUEGPA,EUEGPB 175 176c 177c ======> BOTH SPIN-RESTRICETED AND UNRESTRICTED <====== 178c 179 sop=1.0d0 180 sopp0= param(13) 181 sopp1= param(14) 182 sopp2= param(15) 183 sopp3= param(16) 184 sopp4= param(17) 185c if (ijzy.eq.1) then 186C Parameters for M06-L Correlation 187c sopp0= 6.042374D-01 188c sopp1= 1.776783D+02 189c sopp2= -2.513252D+02 190c sopp3= 7.635173D+01 191c sopp4= -1.255699D+01 192c elseif (ijzy.eq.2) then 193c Parameters for M06-HF Correlation 194c sopp0= 1.674634D+00 195c sopp1= 5.732017D+01 196c sopp2= 5.955416D+01 197c sopp3= -2.311007D+02 198c sopp4= 1.255199D+02 199c elseif (ijzy.eq.3) then 200c Parameters for M06 Correlation 201c sopp0= 3.741539D+00 202c sopp1= 2.187098D+02 203c sopp2= -4.531252D+02 204c sopp3= 2.936479D+02 205c sopp4= -6.287470D+01 206c elseif (ijzy.eq.4) then 207c Parameters for M06-2X Correlation 208c sopp0= 8.833596D-01 209c sopp1= 3.357972D+01 210c sopp2= -7.043548D+01 211c sopp3= 4.978271D+01 212c sopp4= -1.852891D+01 213c endif 214 215#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 216#if defined(NWAD_PRINT) 217 call nwxc_c_vs98_p(param, tol_rho, ipol, nq, wght, rho, 218 & rgamma, tau, func) 219#else 220 call nwxc_c_vs98(param, tol_rho, ipol, nq, wght, rho, 221 & rgamma, tau, func) 222#endif 223#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV) 224 call nwxc_c_vs98_d2(param, tol_rho, ipol, nq, wght, rho, 225 & rgamma, tau, func) 226#else 227 call nwxc_c_vs98_d3(param, tol_rho, ipol, nq, wght, rho, 228 & rgamma, tau, func) 229#endif 230 231 Pi = F4*ATan(F1) 232 F6=6.0d0 233 F43 = F4 / F3 234 Pi34 = F3 / (F4*Pi) 235 F13 = F1 / F3 236 237 do 20 n = 1, nq 238 EUA = 0.0d0 239 EUB = 0.0d0 240 ChiA = 0.0d0 241 ChiB = 0.0d0 242 if (ipol.eq.1) then 243 if (rho(n,R_T).lt.Tol_Rho) goto 20 244c 245c get the density, gradient, and tau for the alpha spin from the total 246c 247 PA = rho(n,R_T)/F2 248c GAA = ( delrho(n,1,1)*delrho(n,1,1) + 249c & delrho(n,2,1)*delrho(n,2,1) + 250c & delrho(n,3,1)*delrho(n,3,1))/F4 251 PB = PA 252 GAA = rgamma(n,G_TT)/F4 253c if(sqrt(gaa).lt.tol_rho) goto 20 254c In the m06css subroutine, we use 2*TA as the tau, so we do not divide 255c the tau by 2 here 256 257 TA = tau(n,T_T) 258 if(ta.lt.tol_rho) goto 30 259 260#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 261#if defined(NWAD_PRINT) 262 Call nwxc_m06css_p(param,Tol_Rho,PA,GAA,TA,FA,FPA,FGA,FTA,EUA, 263 & ChiA,EUPA,ChiAP,ChiAG) 264#else 265 Call nwxc_m06css(param,Tol_Rho,PA,GAA,TA,FA,FPA,FGA,FTA,EUA, 266 & ChiA,EUPA,ChiAP,ChiAG) 267#endif 268#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV) 269 Call nwxc_m06css_d2(param,Tol_Rho,PA,GAA,TA,FA,FPA,FGA,FTA,EUA, 270 & ChiA,EUPA,ChiAP,ChiAG) 271#else 272 Call nwxc_m06css_d3(param,Tol_Rho,PA,GAA,TA,FA,FPA,FGA,FTA,EUA, 273 & ChiA,EUPA,ChiAP,ChiAG) 274#endif 275 GBB = GAA 276 TB = TA 277 FB = FA 278 FPB = FPA 279 FGB = FGA 280 FTB = FTA 281 EUB = EUA 282 ChiB = ChiA 283 EUPB = EUPA 284 ChiBP = ChiAP 285 ChiBG = ChiAG 286 287c Ec = Ec + 2.d0*FA*qwght(n) !factor of 2 account for both spin 288 func(n)=func(n)+ FA*2d0*wght 289c Amat(n,D1_RA) = Amat(n,D1_RA)+ FPA*wght 290c Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + FGA*wght 291c Mmat(n,D1_TA) = Mmat(n,D1_TA) + FTA*wght 292c write (*,*) "PA,GAA,TA",PA,GAA,TA 293c write (*,*) "FPA,FGA,FTA",FPA,FGA,FTA 294c UUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUnrestricted 295 else ! ipol=2 296c 297c ======> SPIN-UNRESTRICTED <====== 298c 299 PA = 0.0d0 300 PB = 0.0d0 301 GAA = 0.0d0 302 GBB = 0.0d0 303 TA = 0.0d0 304 TB = 0.0d0 305c 306c alpha 307c 308 if (rho(n,R_A).le.0.5d0*Tol_Rho) go to 25 309 PA = rho(n,R_A) 310c GAA = delrho(n,1,1)*delrho(n,1,1) + 311c & delrho(n,2,1)*delrho(n,2,1) + 312c & delrho(n,3,1)*delrho(n,3,1) 313c 314c In the m06css subroutine, we use 2*TA as the tau 315c 316 if (tau(n,T_A).le.0.5d0*Tol_Rho) go to 25 317 GAA = rgamma(n,G_AA) 318 TA = 2.0d0*tau(n,T_A) 319#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 320#if defined(NWAD_PRINT) 321 Call nwxc_m06css_p(param,Tol_Rho,PA,GAA,TA,FA,FPA,FGA,FTA,EUA, 322 & ChiA,EUPA,ChiAP,ChiAG) 323#else 324 Call nwxc_m06css(param,Tol_Rho,PA,GAA,TA,FA,FPA,FGA,FTA,EUA, 325 & ChiA,EUPA,ChiAP,ChiAG) 326#endif 327#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV) 328 Call nwxc_m06css_d2(param,Tol_Rho,PA,GAA,TA,FA,FPA,FGA,FTA,EUA, 329 & ChiA,EUPA,ChiAP,ChiAG) 330#else 331 Call nwxc_m06css_d3(param,Tol_Rho,PA,GAA,TA,FA,FPA,FGA,FTA,EUA, 332 & ChiA,EUPA,ChiAP,ChiAG) 333#endif 334c Ec = Ec + FA*qwght(n) 335 func(n)=func(n)+ FA*wght 336c Amat(n,D1_RA) = Amat(n,D1_RA)+ FPA*wght 337c Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + FGA*wght 338c Mmat(n,D1_TA) = Mmat(n,D1_TA) + FTA*wght 339c 340c In the m06css subroutine, we use 2*TB as the tau, 341c 342c 343c Beta 344c 345 25 continue 346 if (rho(n,R_B).le.0.5d0*Tol_Rho) go to 30 347 PB = rho(n,R_B) 348c GBB = delrho(n,1,2)*delrho(n,1,2) + 349c & delrho(n,2,2)*delrho(n,2,2) + 350c & delrho(n,3,2)*delrho(n,3,2) 351 352 if (tau(n,T_B).le.0.5d0*Tol_Rho) go to 30 353 GBB = rgamma(n,G_BB) 354 TB = 2.0d0*tau(n,T_B) 355#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 356#if defined(NWAD_PRINT) 357 Call nwxc_m06css_p(param,Tol_Rho,PB,GBB,TB,FB,FPB,FGB,FTB,EUB, 358 & ChiB,EUPB,ChiBP,ChiBG) 359#else 360 Call nwxc_m06css(param,Tol_Rho,PB,GBB,TB,FB,FPB,FGB,FTB,EUB, 361 & ChiB,EUPB,ChiBP,ChiBG) 362#endif 363#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV) 364 Call nwxc_m06css_d2(param,Tol_Rho,PB,GBB,TB,FB,FPB,FGB,FTB,EUB, 365 & ChiB,EUPB,ChiBP,ChiBG) 366#else 367 Call nwxc_m06css_d3(param,Tol_Rho,PB,GBB,TB,FB,FPB,FGB,FTB,EUB, 368 & ChiB,EUPB,ChiBP,ChiBG) 369#endif 370c Ec = Ec + FB*qwght(n) 371 func(n)=func(n)+ FB*wght 372c Amat(n,D1_RB) = Amat(n,D1_RB)+ FPB 373c Cmat(n,D1_GBB) = Cmat(n,D1_GBB) + FGB 374c Mmat(n,D1_TB) = Mmat(n,D1_TB) + FTB 375 endif 376 377 30 continue 378 P = PA + PB 379 380 If((PA.gt.0.5d0*Tol_Rho).or.(PB.gt.0.5d0*Tol_Rho)) then 381 RS = (Pi34/P) ** F13 382c RSP = -RS/(F3*P) 383 Zeta = (PA-PB)/P 384c dZdA = (F1-Zeta)/P 385c dZdB = (-F1-Zeta)/P 386#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 387#if defined(NWAD_PRINT) 388 Call nwxc_c_lsda_p(tol_rho, 389 R RS,Zeta,PotLC,dLdS,dLdZ,d2LdSS,d2LdSZ, 390 $ d2LdZZ) 391#else 392 Call nwxc_c_lsda(tol_rho, 393 R RS,Zeta,PotLC,dLdS,dLdZ,d2LdSS,d2LdSZ, 394 $ d2LdZZ) 395#endif 396#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV) 397 Call nwxc_c_lsda_d2(tol_rho, 398 R RS,Zeta,PotLC,dLdS,dLdZ,d2LdSS,d2LdSZ, 399 $ d2LdZZ) 400#else 401 Call nwxc_c_lsda_d3(tol_rho, 402 R RS,Zeta,PotLC,dLdS,dLdZ,d2LdSS,d2LdSZ, 403 $ d2LdZZ) 404#endif 405 EUEG = P*PotLC - EUA - EUB 406 U = COpp*(ChiA+ChiB)/(F1 + COpp*(ChiA+ChiB)) 407 W = sopp0+U*(sopp1+U*(sopp2+U*(sopp3+U*sopp4))) 408c Ec = Ec + sop*EUEG*W*qwght(n) 409 func(n)=func(n)+ sop*EUEG*W*wght 410c dUdChiA =COpp/(F1 + COpp*(ChiA+ChiB))**2 411c dUdChiB =COpp/(F1 + COpp*(ChiA+ChiB))**2 412c dUdPA= dUdChiA*ChiAP 413c dUdPB= dUdChiB*ChiBP 414c dUdGA= dUdChiA*ChiAG 415c dUdGB= dUdChiB*ChiBG 416c dWdU =sopp1+U*(F2*sopp2+U*(F3*sopp3+U*F4*sopp4)) 417c dWdPA= dWdU*dUdPA 418c dWdPB= dWdU*dUdPB 419c dWdGA= dWdU*dUdGA 420c dWdGB= dWdU*dUdGB 421c EUEGPA = PotLC + P*dLdS*RSP + P*dLdZ*dZdA - EUPA 422c EUEGPB = PotLC + P*dLdS*RSP + P*dLdZ*dZdB - EUPB 423c if (ipol.eq.1) then 424c Amat(n,D1_RA) = Amat(n,D1_RA) 425c + + sop*(EUEGPA*W + EUEG*dWdPA)*wght 426c Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + sop*(EUEG*dWdGA)*wght 427c else 428c Amat(n,D1_RA) = Amat(n,D1_RA) 429c + + sop*(EUEGPA*W + EUEG*dWdPA)*wght 430c Amat(n,D1_RB) = Amat(n,D1_RB) 431c + + sop*(EUEGPB*W + EUEG*dWdPB)*wght 432c Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + sop*EUEG*dWdGA*wght 433c Cmat(n,D1_GBB) = Cmat(n,D1_GBB) + sop*(EUEG*dWdGB)*wght 434c endif 435 endIf 436c write (*,*) "PA, PB, GAA, GBB,ipol",PA, PB, GAA, GBB,ipol 437c write (*,*) "FA, FB,FGA, FGB",FA, FB,FGA, FGB 438c Stop 43920 continue 440 end 441 442 443#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 444#if defined(NWAD_PRINT) 445 Subroutine nwxc_m06css_p(param,Tol_Rho,PX,GX,TX,F,FP,FG,FT,EUEG, 446 & Chi,EUEGP,ChiP,ChiG) 447#else 448 Subroutine nwxc_m06css(param,Tol_Rho,PX,GX,TX,F,FP,FG,FT,EUEG, 449 & Chi,EUEGP,ChiP,ChiG) 450#endif 451#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV) 452 Subroutine nwxc_m06css_d2(param,Tol_Rho,PX,GX,TX,F,FP,FG,FT,EUEG, 453 & Chi,EUEGP,ChiP,ChiG) 454#else 455 Subroutine nwxc_m06css_d3(param,Tol_Rho,PX,GX,TX,F,FP,FG,FT,EUEG, 456 & Chi,EUEGP,ChiP,ChiG) 457#endif 458#include "nwad.fh" 459 Implicit none 460c 461#include "intf_nwxc_c_lsda.fh" 462c 463C 464C Compute the same-spin part of the m06 correlation functional for one grid 465C point and one spin-case. 466C 467C 468#if defined(NWAD_PRINT) 469#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 470 type(nwad_dble)::param(22) 471 type(nwad_dble)::sss0, sss1, sss2, sss3, sss4 472#else 473 double precision param(22) 474 double precision sss0, sss1, sss2, sss3, sss4 475#endif 476#else 477 double precision param(22) 478 double precision sss0, sss1, sss2, sss3, sss4 479#endif 480 type(nwad_dble)::PX, GX, TX, F, EUEG, Fscc, Chi 481 type(nwad_dble)::Rs, D, E, W, U, Zeta, PotLC 482 double precision FP, FG, FT, Tol_Rho 483 double precision EUEGP, ChiP, ChiG 484 double precision Zero, Pt25, F1, F2, F3, F4, F5, F6, F8, F11 485 double precision ss, Css 486 double precision Pi, Pi34, F13, F23, F43, F53, F83, F113 487 double precision FDUEG, RSP, dFsccP, dFsccG 488 double precision dFsccT, dUdChi, dWdU, dWdP, dWdG 489 double precision d2LdSS,d2LdSZ,d2LdZZ,dLdS,dLdZ 490 491 492 493 Data Zero/0.0d0/, Pt25/0.25d0/, F1/1.0d0/, F2/2.0d0/, F3/3.0d0/, 494 $ F4/4.0d0/, F5/5.0d0/, F6/6.0d0/, F8/8.0d0/, F11/11.0d0/, 495 $ Css/0.06d0/ 496C 497c Tol_Rho=1.0D-7 498c write (*,*) Tol_Rho 499 ss=1.0 500 sss0= param(18) 501 sss1= param(19) 502 sss2= param(20) 503 sss3= param(21) 504 sss4= param(22) 505c if (ijzy.eq.1) then 506C Parameters for M06-L Correlation 507c sss0= 5.349466D-01 508c sss1= 5.396620D-01 509c sss2= -3.161217D+01 510c sss3= 5.149592D+01 511c sss4= -2.919613D+01 512c elseif (ijzy.eq.2) then 513c Parameters for M06-HF Correlation 514c sss0= 1.023254D-01 515c sss1= -2.453783D+00 516c sss2= 2.913180D+01 517c sss3= -3.494358D+01 518c sss4= 2.315955D+01 519c elseif (ijzy.eq.3) then 520c Parameters for M06 Correlation 521c sss0= 5.094055D-01 522c sss1= -1.491085D+00 523c sss2= 1.723922D+01 524c sss3= -3.859018D+01 525c sss4= 2.845044D+01 526c elseif (ijzy.eq.4) then 527c Parameters for M06-2X Correlation 528c sss0= 3.097855D-01 529c sss1= -5.528642D+00 530c sss2= 1.347420D+01 531c sss3= -3.213623D+01 532c sss4= 2.846742D+01 533c endif 534 535 If ((PX.le.Tol_Rho)) then 536 EUEG = Zero 537 Chi = Zero 538 EUEGP = Zero 539 ChiP = Zero 540 ChiG = Zero 541 PX = Zero 542 GX = Zero 543 TX = Zero 544 F = Zero 545 FP = Zero 546 FG = Zero 547 FT = Zero 548 else 549 Pi = F4*ATan(F1) 550 Pi34 = F3 / (F4*Pi) 551 F13 = F1 / F3 552 F23 = F2 / F3 553 F43 = F2 * F23 554 F53 = F5 / F3 555 F83 = F8 / F3 556 F113 = F11 / F3 557 FDUEG = (F3/F5)*(F6*Pi*Pi)**F23 558 RS = (Pi34/PX) ** F13 559 Zeta = F1 560#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 561#if defined(NWAD_PRINT) 562 Call nwxc_c_lsda_p(tol_rho, 563 R RS,Zeta,PotLC,dLdS,dLdZ,d2LdSS,d2LdSZ,d2LdZZ) 564#else 565 Call nwxc_c_lsda(tol_rho, 566 R RS,Zeta,PotLC,dLdS,dLdZ,d2LdSS,d2LdSZ,d2LdZZ) 567#endif 568#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV) 569 Call nwxc_c_lsda_d2(tol_rho, 570 R RS,Zeta,PotLC,dLdS,dLdZ,d2LdSS,d2LdSZ,d2LdZZ) 571#else 572 Call nwxc_c_lsda_d3(tol_rho, 573 R RS,Zeta,PotLC,dLdS,dLdZ,d2LdSS,d2LdSZ,d2LdZZ) 574#endif 575 EUEG = PX*PotLC 576 D = TX - Pt25*GX/PX 577C DUEG = FDUEG*PX**F53 578 Chi = GX/(PX**F83) 579 U = Css*Chi/(F1 + Css*Chi) 580 W = sss0+U*(sss1+U*(sss2+U*(sss3+U*sss4))) 581 Fscc=D/TX 582 E = Fscc*W*EUEG 583 F = E*ss 584c RSP = -RS/(F3*Px) 585c ChiG = F1/PX**F83 586c ChiP = -F83*Chi/PX 587c dFsccP=Pt25*GX/(TX*PX**2) 588c dFsccG=-Pt25/(TX*PX) 589c dFsccT=Pt25*GX/(PX*TX**2) 590c dUdChi=Css/((F1+Css*Chi)**2) 591c dWdU=sss1+U*(F2*sss2+U*(F3*sss3+U*F4*sss4)) 592c dWdP=dWdU*dUdChi*ChiP 593c dWdG=dWdU*dUdChi*ChiG 594c EUEGP = PotLC + PX*dLdS*RSP 595c FP = ss*(dFsccP*W*EUEG 596c $ + Fscc*dWdP*EUEG 597c $ + Fscc*W*EUEGP) 598c FG = ss*(dFsccG*W*EUEG 599c $ + Fscc*dWdG*EUEG) 600 601c FT = ss*(dFsccT*W*EUEG) 602 Endif 603 604 Return 605 End 606#ifndef NWAD_PRINT 607#define NWAD_PRINT 608c 609c Compile source again for the 2nd derivative case 610c 611#include "nwxc_c_m06.F" 612#endif 613#ifndef SECOND_DERIV 614#define SECOND_DERIV 615c 616c Compile source again for the 2nd derivative case 617c 618#include "nwxc_c_m06.F" 619#endif 620#ifndef THIRD_DERIV 621#define THIRD_DERIV 622c 623c Compile source again for the 3rd derivative case 624c 625#include "nwxc_c_m06.F" 626#endif 627#undef NWAD_PRINT 628C> @} 629