1c M06 suite correlation functional 2c META GGA 3C utilizes ingredients: 4c rho - density 5c delrho - gradient of density 6c tau (tauN)- K.S kinetic energy density 7c ijzy - 1 M06-L 8c ijzy - 2 M06-HF 9c ijzy - 3 M06 10c ijzy - 4 M06-2X 11c ijzy - 5 revM06-L 12c ijzy - 6 revM06 13 14 Subroutine xc_cm06(tol_rho, cfac, lcfac, nlcfac, rho, delrho, 15 & nq, ipol, Ec, qwght, ldew, func, 16 & tau, Amat, Cmat, Mmat,ijzy) 17 18 19c 20c$Id$ 21c 22c 23c [a] Zhao, Y. and Truhlar, D. G. J. Chem. Phys. 2006, 125, 194101; 24c [b] Zhao, Y. and Truhlar, D. G. J. Phys. Chem. A (2006),110(49),13126-13130. 25c [c] Grafenstein, J., Izotov, D. and Cremer, D. J. Chem. Phys. 2007, 127, 214103. 26 27 28 implicit none 29c 30#include "errquit.fh" 31c 32c 33c Input and other parameters 34c 35 integer ipol, nq 36 37 double precision cfac 38 logical lcfac, nlcfac 39 40 logical lfac, nlfac 41 double precision fac 42 double precision tol_rho 43 44c 45c Threshold parameters 46c 47 double precision F1, F2, F3, F4,COpp 48 Data COpp/0.0031d0/,F1/1.0d0/,F2/2.0d0/, 49 & F3/3.0d0/,F4/4.0d0/ 50c 51c Correlation energy 52c 53 double precision Ec 54c 55c Charge Density 56c 57 double precision rho(nq,ipol*(ipol+1)/2) 58c 59c Charge Density Gradient 60c 61 double precision delrho(nq,3,ipol), gammaval, gam12 62 63c 64c Kinetic Energy Density 65c 66 double precision tau(nq,ipol), tauN 67 68c 69c Quadrature Weights 70c 71 double precision qwght(nq) 72c 73 logical ldew 74 double precision func(*) 75c 76c Sampling Matrices for the XC Potential 77c 78 double precision Amat(nq,ipol), Cmat(nq,*) 79 double precision Mmat(nq,*) 80 81 integer n, ijzy 82 83c call to the m05css subroutine 84 double precision PA,GAA,TA,FA,FPA,FGA,FTA,EUA,EUEGA,ChiA,EUPA 85 &,ChiAP,ChiAG 86 double precision PB,GBB,TB,FB,FPB,FGB,FTB,EUB,EUEGB,ChiB,EUPB 87 &,ChiBP,ChiBG 88c 89 double precision sop, sopp0, sopp1,sopp2, sopp3, sopp4 90 double precision Pi, F6, F43, Pi34, F13, 91 &RS,RSP,Zeta,dZdA,dZdB,PotLC,dLdS,dLdZ,d2LdSS,d2LdSZ,d2LdZZ 92 double precision P, EUEG, U, W 93 double precision dUdChiA,dUdChiB,dUdPA,dUdPB,dUdGA,dUdGB, 94 &dWdU,dWdPA,dWdPB, dWdGA,dWdGB,EUEGPA,EUEGPB 95 96c 97c ======> BOTH SPIN-RESTRICETED AND UNRESTRICTED <====== 98c 99 sop=1.0d0 100 sopp0= 0d0 101 sopp1= 0d0 102 sopp2= 0d0 103 sopp3= 0d0 104 sopp4= 0d0 105 if (ijzy.eq.1) then 106C Parameters for M06-L Correlation 107 sopp0= 6.042374D-01 108 sopp1= 1.776783D+02 109 sopp2= -2.513252D+02 110 sopp3= 7.635173D+01 111 sopp4= -1.255699D+01 112 elseif (ijzy.eq.2) then 113C Parameters for M06-HF Correlation 114 sopp0= 1.674634D+00 115 sopp1= 5.732017D+01 116 sopp2= 5.955416D+01 117 sopp3= -2.311007D+02 118 sopp4= 1.255199D+02 119 elseif (ijzy.eq.3) then 120C Parameters for M06 Correlation 121 sopp0= 3.741539D+00 122 sopp1= 2.187098D+02 123 sopp2= -4.531252D+02 124 sopp3= 2.936479D+02 125 sopp4= -6.287470D+01 126 elseif (ijzy.eq.4) then 127C Parameters for M06-2X Correlation 128 sopp0= 8.833596D-01 129 sopp1= 3.357972D+01 130 sopp2= -7.043548D+01 131 sopp3= 4.978271D+01 132 sopp4= -1.852891D+01 133 elseif (ijzy.eq.5) then 134C Parameters for revM06-L Correlation 135 sopp0= 3.44360696D-01 136 sopp1= -5.57080242D-01 137 sopp2= -2.009821162D+00 138 sopp3= -1.857641887D+00 139 sopp4= -1.076639864D+00 140 elseif (ijzy.eq.6) then 141C Parameters for revM06 Correlation 142 sopp0= 1.222401598D+00 143 sopp1= 6.613907336D-01 144 sopp2= -1.884581043D+00 145 sopp3= -2.780360568D+00 146 sopp4= -3.068579344D+00 147 else 148 call errquit("xc_cm06: illegal value of ijzy",ijzy,UERR) 149 endif 150 151 call xc_cvs98(tol_rho, 1.0d0, lfac, nlcfac, 152 & rho, delrho, nq, ipol, 153 & Ec, qwght, ldew,func,tau,Amat,Cmat,Mmat,ijzy+1) 154 155 156 157 Pi = F4*ATan(F1) 158 F6=6.0d0 159 F43 = F4 / F3 160 Pi34 = F3 / (F4*Pi) 161 F13 = F1 / F3 162 163 do 20 n = 1, nq 164 if (rho(n,1).lt.Tol_Rho) goto 20 165 if (ipol.eq.1) then 166c 167c get the density, gradient, and tau for the alpha spin from the total 168c 169 PA = rho(n,1)/F2 170 GAA = ( delrho(n,1,1)*delrho(n,1,1) + 171 & delrho(n,2,1)*delrho(n,2,1) + 172 & delrho(n,3,1)*delrho(n,3,1))/F4 173 if(sqrt(gaa).lt.tol_rho) goto 20 174c In the m05css subroutine, we use 2*TA as the tau, so we do not divide 175c the tau by 2 here 176 177 TA = tau(n,1) 178 if(ta.lt.tol_rho) goto 20 179 180 Call m06css(Tol_Rho,PA,GAA,TA,FA,FPA,FGA,FTA,EUA, 181 & ChiA,EUPA,ChiAP,ChiAG,ijzy) 182 PB = PA 183 GBB = GAA 184 TB = TA 185 FB = FA 186 FPB = FPA 187 FGB = FGA 188 FTB = FTA 189 EUB = EUA 190 ChiB = ChiA 191 EUPB = EUPA 192 ChiBP = ChiAP 193 ChiBG = ChiAG 194 195 Ec = Ec + 2.d0*FA*qwght(n) !factor of 2 account for both spin 196 if(ldew) func(n)=func(n)+ FA*2d0 197 Amat(n,1)=Amat(n,1)+ FPA 198 Cmat(n,1)= Cmat(n,1) + FGA 199 Mmat(n,1)= Mmat(n,1) + FTA 200c write (*,*) "PA,GAA,TA",PA,GAA,TA 201c write (*,*) "FPA,FGA,FTA",FPA,FGA,FTA 202c UUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUnrestricted 203 else ! ipol=2 204c 205c ======> SPIN-UNRESTRICTED <====== 206c 207 ChiA = 0.0d0 208 ChiAP = 0.0d0 209 ChiAG = 0.0d0 210 ChiB = 0.0d0 211 ChiBP = 0.0d0 212 ChiBG = 0.0d0 213 EUA = 0.0d0 214 EUB = 0.0d0 215 EUPA = 0.0d0 216 EUPB = 0.0d0 217c 218c alpha 219c 220 221 PA = rho(n,2) 222 if (PA.le.Tol_Rho) go to 25 223 GAA = delrho(n,1,1)*delrho(n,1,1) + 224 & delrho(n,2,1)*delrho(n,2,1) + 225 & delrho(n,3,1)*delrho(n,3,1) 226c 227c In the m05css subroutine, we use 2*TA as the tau 228c 229 TA = 2*tau(n,1) 230 231 if(sqrt(gaa).lt.tol_rho) goto 25 232 if(ta.lt.tol_rho) goto 25 233 234 Call m06css(Tol_Rho,PA,GAA,TA,FA,FPA,FGA,FTA,EUA, 235 & ChiA,EUPA,ChiAP,ChiAG,ijzy) 236 Ec = Ec + FA*qwght(n) 237 if(ldew) func(n)=func(n)+ FA 238 Amat(n,1)=Amat(n,1)+ FPA 239 Cmat(n,1)= Cmat(n,1) + FGA 240 Mmat(n,1)= Mmat(n,1) + FTA 241c 242c In the m05css subroutine, we use 2*TB as the tau, 243c 244c 245c Beta 246c 247 25 continue 248 PB = rho(n,3) 249 if (PB.le.Tol_Rho) go to 30 250 GBB = delrho(n,1,2)*delrho(n,1,2) + 251 & delrho(n,2,2)*delrho(n,2,2) + 252 & delrho(n,3,2)*delrho(n,3,2) 253 254 TB = 2*tau(n,2) 255 256 if(sqrt(gbb).lt.tol_rho) goto 30 257 if(tb.lt.tol_rho) goto 30 258 259 Call m06css(Tol_Rho,PB,GBB,TB,FB,FPB,FGB,FTB,EUB, 260 & ChiB,EUPB,ChiBP,ChiBG,ijzy) 261 Ec = Ec + FB*qwght(n) 262 if(ldew) func(n)=func(n)+ FB 263 Amat(n,2)= Amat(n,2)+ FPB 264 Cmat(n,3)= Cmat(n,3) + FGB 265 Mmat(n,2)= Mmat(n,2) + FTB 266 endif 267 268 30 continue 269 P = PA + PB 270 271 If((PA.gt.Tol_Rho).and.(PB.gt.Tol_Rho)) then 272 RS = (Pi34/P) ** F13 273 RSP = -RS/(F3*P) 274 Zeta = (PA-PB)/P 275 dZdA = (F1-Zeta)/P 276 dZdB = (-F1-Zeta)/P 277 Call lsdac(tol_rho, 278 R RS,Zeta,PotLC,dLdS,dLdZ,d2LdSS,d2LdSZ, 279 $ d2LdZZ) 280 EUEG = P*PotLC - EUA - EUB 281 U = COpp*(ChiA+ChiB)/(F1 + COpp*(ChiA+ChiB)) 282 W = sopp0+U*(sopp1+U*(sopp2+U*(sopp3+U*sopp4))) 283 Ec = Ec + sop*EUEG*W*qwght(n) 284 if(ldew) func(n)=func(n)+ sop*EUEG*W 285 dUdChiA =COpp/(F1 + COpp*(ChiA+ChiB))**2 286 dUdChiB =COpp/(F1 + COpp*(ChiA+ChiB))**2 287 dUdPA= dUdChiA*ChiAP 288 dUdPB= dUdChiB*ChiBP 289 dUdGA= dUdChiA*ChiAG 290 dUdGB= dUdChiB*ChiBG 291 dWdU =sopp1+U*(F2*sopp2+U*(F3*sopp3+U*F4*sopp4)) 292 dWdPA= dWdU*dUdPA 293 dWdPB= dWdU*dUdPB 294 dWdGA= dWdU*dUdGA 295 dWdGB= dWdU*dUdGB 296 EUEGPA = PotLC + P*dLdS*RSP + P*dLdZ*dZdA - EUPA 297 EUEGPB = PotLC + P*dLdS*RSP + P*dLdZ*dZdB - EUPB 298 if (ipol.eq.1) then 299 Amat(n,1) = Amat(n,1) + sop*(EUEGPA*W + EUEG*dWdPA) 300 Cmat(n,1)= Cmat(n,1) + sop*(EUEG*dWdGA) 301 else 302 Amat(n,1) = Amat(n,1) + sop*(EUEGPA*W + EUEG*dWdPA) 303 Amat(n,2) = Amat(n,2) + sop*(EUEGPB*W + EUEG*dWdPB) 304 Cmat(n,1) = Cmat(n,1) + sop*EUEG*dWdGA 305 Cmat(n,3) = Cmat(n,3) + sop*(EUEG*dWdGB) 306 endif 307 endIf 308c write (*,*) "PA, PB, GAA, GBB,ipol",PA, PB, GAA, GBB,ipol 309c write (*,*) "FA, FB,FGA, FGB",FA, FB,FGA, FGB 310c Stop 31120 continue 312 end 313 314 Subroutine xc_cm06_d2() 315 implicit none 316 call errquit(' xc06: d2 not coded ',0,0) 317 return 318 end 319 320 321 322 323 Subroutine m06css(Tol_Rho,PX,GX,TX,F,FP,FG,FT,EUEG,Chi,EUEGP, 324 & ChiP,ChiG,ijzy) 325 Implicit none 326c 327#include "errquit.fh" 328C 329C Compute the same-spin part of the m05 correlation functional for one grid 330C point and one spin-case. 331C 332C 333 integer ijzy 334 double precision PX, GX, TX, F, FP, FG, FT, Tol_Rho 335 double precision EUEG, Chi, EUEGP, ChiP, ChiG 336 double precision Zero, Pt25, F1, F2, F3, F4, F5, F6, F8, F11 337 double precision ss, sss0,sss1, sss2, sss3, sss4, Css 338 double precision Pi, Pi34, F13, F23, F43, F53, F83, F113 339 double precision RS, FDUEG, D, Fscc, RSP, dFsccP, dFsccG 340 double precision E, W, U, dFsccT, dUdChi, dWdU, dWdP, dWdG 341 double precision d2LdSS,d2LdSZ,d2LdZZ,PotLC,dLdS,dLdZ 342 343 344 345 Data Zero/0.0d0/, Pt25/0.25d0/, F1/1.0d0/, F2/2.0d0/, F3/3.0d0/, 346 $ F4/4.0d0/, F5/5.0d0/, F6/6.0d0/, F8/8.0d0/, F11/11.0d0/, 347 $ Css/0.06d0/ 348C 349c Tol_Rho=1.0D-7 350c write (*,*) Tol_Rho 351 sss0= 0d0 352 sss1= 0d0 353 sss2= 0d0 354 sss3= 0d0 355 sss4= 0d0 356 ss=1.0 357 if (ijzy.eq.1) then 358C Parameters for M06-L Correlation 359 sss0= 5.349466D-01 360 sss1= 5.396620D-01 361 sss2= -3.161217D+01 362 sss3= 5.149592D+01 363 sss4= -2.919613D+01 364 elseif (ijzy.eq.2) then 365C Parameters for M06-HF Correlation 366 sss0= 1.023254D-01 367 sss1= -2.453783D+00 368 sss2= 2.913180D+01 369 sss3= -3.494358D+01 370 sss4= 2.315955D+01 371 elseif (ijzy.eq.3) then 372C Parameters for M06 Correlation 373 sss0= 5.094055D-01 374 sss1= -1.491085D+00 375 sss2= 1.723922D+01 376 sss3= -3.859018D+01 377 sss4= 2.845044D+01 378 elseif (ijzy.eq.4) then 379C Parameters for M06-2X Correlation 380 sss0= 3.097855D-01 381 sss1= -5.528642D+00 382 sss2= 1.347420D+01 383 sss3= -3.213623D+01 384 sss4= 2.846742D+01 385 elseif (ijzy.eq.5) then 386C Parameters for revM06-L Correlation 387 sss0= 1.227659748D+00 388 sss1= 8.55201283D-01 389 sss2= -3.113346677D+00 390 sss3= -2.239678026D+00 391 sss4= 3.54638962D-01 392 elseif (ijzy.eq.6) then 393C Parameters for revM06 Correlation 394 sss0= 9.017224575D-01 395 sss1= 2.079991827D-01 396 sss2= -1.823747562D+00 397 sss3= -1.384430429D+00 398 sss4= -4.423253381D-01 399 else 400 call errquit("m06css: illegal value of ijzy",ijzy,UERR) 401 endif 402 403 If ((PX.le.Tol_Rho)) then 404 EUEG = Zero 405 Chi = Zero 406 EUEGP = Zero 407 ChiP = Zero 408 ChiG = Zero 409 PX = Zero 410 GX = Zero 411 TX = Zero 412 F = Zero 413 FP = Zero 414 FG = Zero 415 FT = Zero 416 else 417 Pi = F4*ATan(F1) 418 Pi34 = F3 / (F4*Pi) 419 F13 = F1 / F3 420 F23 = F2 / F3 421 F43 = F2 * F23 422 F53 = F5 / F3 423 F83 = F8 / F3 424 F113 = F11 / F3 425 FDUEG = (F3/F5)*(F6*Pi*Pi)**F23 426 RS = (Pi34/PX) ** F13 427 Call lsdac(tol_rho, 428 R RS,F1,PotLC,dLdS,dLdZ,d2LdSS,d2LdSZ,d2LdZZ) 429 EUEG = PX*PotLC 430 D = TX - Pt25*GX/PX 431C DUEG = FDUEG*PX**F53 432 Chi = GX/PX**F83 433 U = Css*Chi/(F1 + Css*Chi) 434 W = sss0+U*(sss1+U*(sss2+U*(sss3+U*sss4))) 435 Fscc=D/TX 436 E = Fscc*W*EUEG 437 F = E*ss 438 RSP = -RS/(F3*Px) 439 ChiG = F1/PX**F83 440 ChiP = -F83*Chi/PX 441 dFsccP=Pt25*GX/(TX*PX**2) 442 dFsccG=-Pt25/(TX*PX) 443 dFsccT=Pt25*GX/(PX*TX**2) 444 dUdChi=Css/((F1+Css*Chi)**2) 445 dWdU=sss1+U*(F2*sss2+U*(F3*sss3+U*F4*sss4)) 446 dWdP=dWdU*dUdChi*ChiP 447 dWdG=dWdU*dUdChi*ChiG 448 EUEGP = PotLC + PX*dLdS*RSP 449 FP = ss*(dFsccP*W*EUEG 450 $ + Fscc*dWdP*EUEG 451 $ + Fscc*W*EUEGP) 452 FG = ss*(dFsccG*W*EUEG 453 $ + Fscc*dWdG*EUEG) 454 455 FT = ss*(dFsccT*W*EUEG) 456 Endif 457 458 Return 459 End 460 461 462