1* 2* $Id$ 3* 4 5*** a = 0.2 am=0.8 6**** a*HF + (1-a)*Dirac - 0.72*Dirac + 0.72*Becke88 + 0.81*LYP + (1-0.81)*Vosko 7**** a*HF + (1-a-0.72)*Dirac + 0.72*Becke88 + 0.81*LYP + 0.19*Vosko 8**** a*HF + (am/0.8)*((0.28-0.2)*Dirac + 0.72*Becke88) + 0.81*LYP + 0.19*Vosko 9**** a*HF + (am/0.8)*(0.08*Dirac + 0.72*Becke88) + 0.81*LYP + 0.19*Vosko 10**** a*HF + (am)*(0.1*Dirac + 0.9*Becke88) + 0.81*LYP + 0.19*Vosko 11**** a*HF + (1-a)*(0.1*Dirac + 0.9*Becke88) + 0.81*LYP + 0.19*Vosko 12 13* ************************************************** 14* * * 15* * gen_B3LYP_BW_unrestricted * 16* * * 17* ************************************************** 18 19 subroutine gen_B3LYP_BW_unrestricted(n2ft3d, 20 & dn_in,agr_in, 21 & x_parameter,c_parameter, 22 & xce,fn,fdn) 23 24 implicit none 25 26 integer n2ft3d 27 real*8 dn_in(n2ft3d,2) 28 real*8 agr_in(n2ft3d,3) 29 real*8 xce(n2ft3d) 30 real*8 fn(n2ft3d,2) 31 real*8 fdn(n2ft3d,3) 32 real*8 x_parameter, c_parameter 33 34 35c*---- parameters given by vosko et al -----------------* 36 real*8 ap,af,x0p,x0f,bp,bf,cpt,cft 37 parameter (ap = 3.109070d-02, af = 1.554530d-02) 38 parameter (x0p=-1.049800d-01, x0f=-3.250000d-01) 39 parameter (bp = 3.727440d+00, bf = 7.060420d+00) 40 parameter (cpt = 1.293520d+01, cft = 1.805780d+01) 41c*------------------------------------------------------* 42 43* constants calculated from vosko's parameters 44 real*8 xp,xf,qp,qf,xx0p,xx0f,fc1,fd1,crs 45 real*8 cp1,cp2,cp3,cp4,cp5,cp6,dp1,dp2,dp3,dp4,dp5,dp6,dp7 46 real*8 cf1,cf2,cf3,cf4,cf5,cf6,df1,df2,df3,df4,df5,df6,df7 47 parameter (xp = -4.581653d-01, xf = -5.772521d-01) 48 parameter (qp = 6.151991d+00, qf = 4.730927d+00) 49 parameter (xx0p= 1.255491d+01, xx0f= 1.586879d+01) 50 parameter (cp1 = 3.109070d-02, cf1 = 1.554530d-02) 51 parameter (cp2 = 9.690228d-04, cf2 = 2.247860d-03) 52 parameter (cp3 = 1.049800d-01, cf3 = 3.250000d-01) 53 parameter (cp4 = 3.878329d-02, cf4 = 5.249122d-02) 54 parameter (cp5 = 3.075995d+00, cf5 = 2.365463d+00) 55 parameter (cp6 = 1.863720d+00, cf6 = 3.530210d+00) 56 parameter (dp1 = 6.218140d-02, df1 = 3.109060d-02) 57 parameter (dp2 = 1.938045d-03, df2 = 4.495720d-03) 58 parameter (dp3 = 1.049800d-01, df3 = 3.250000d-01) 59 parameter (dp4 = -3.205972d-02, df4 = -1.779316d-02) 60 parameter (dp5 = -1.192972d-01, df5 = -1.241661d-01) 61 parameter (dp6 = 1.863720d+00, df6 = 3.530210d+00) 62 parameter (dp7 = 9.461748d+00, df7 = 5.595417d+00) 63 parameter (fc1 = 1.923661d+00, fd1 = 2.564881d+00) 64 parameter (crs = 7.876233d-01) 65 real*8 one6th 66 parameter (one6th=1.0d0/6.0d0) 67 68 69 real*8 ETA,ETA2,DNS_CUT 70 parameter (ETA = 0.0d0) 71 parameter (ETA2 = 1.0d-20) 72 parameter (DNS_CUT = 1.0d-18) 73 74 real*8 tolrho,minagr 75 parameter (tolrho=2.0e-11,minagr=1.0e-12) 76 !parameter (tolrho=2.0e-8,minagr=1.0e-10) 77 78* ***LYP parameters**************** 79 real*8 a,b,c,d 80 parameter (a = 0.04918d0) 81 parameter (b = 0.132d0) 82 parameter (c = 0.2533d0) 83 parameter (d = 0.349d0) 84* ***Becke88 parameters************* 85 real*8 beta 86 parameter (beta = 0.0042d0) 87 88 real*8 thrd,AA 89 parameter (thrd = 1.0d0/3.0d0) 90 91 integer j 92 real*8 twthrd,frthrd,fvthrd,snthrd,etthrd 93 real*8 pi, Cf,lda_c 94 95 real*8 nup,ndn,n,agrup,agrdn,agr 96 real*8 agr2, agrup2,agrdn2 97 real*8 n2,nup2,ndn2 98 real*8 gamma,gamma_u,gamma_uu 99 real*8 gamma_d,gamma_dd,gamma_du 100 real*8 F,F_u,F_uu 101 real*8 F_d,F_dd,F_du 102 real*8 n13d 103 real*8 enthrd 104 real*8 G,G_u,G_uu 105 real*8 G_d,G_dd,G_du 106 real*8 fc,fclda,H0,Hu,Hd 107 real*8 ce 108 real*8 Hu_u,Hu_d,Hd_d,Hd_u 109 real*8 fc_u,fc_d 110 real*8 H0_u,H0_d 111 real*8 fclda_u,fclda_d 112 real*8 fc_agr, fc_agrup,fc_agrdn 113 114 real*8 chiup,chidn 115 real*8 chiup2,chidn2,chiupSQ,chidnSQ 116 real*8 Kup,Kdn,F1up,F1dn,F2up,F2dn 117 real*8 xeup_pbe,fdnxup_pbe,fdagrxup_pbe 118 real*8 xedn_pbe,fdnxdn_pbe,fdagrxdn_pbe 119 real*8 xeup,xedn,xe 120 real*8 fdnxup,fdnxdn,fdagrxup,fdagrxdn 121 real*8 Q,Q1,P,P1,n_m,n2_m,n3_m,n4_m 122 real*8 n13d_m,n13d2_m,n13d3_m,d3 123 real*8 n_thrd,nupndn 124 real*8 n_mthrd,n_mfrthrd,n_mfvthrd 125 real*8 nup_etthrd,nup_fvthrd 126 real*8 ndn_etthrd,ndn_fvthrd 127 real*8 x,xi,ff,dff,xxp,dxxp,xxf,dxxf 128 real*8 ex_p,ux_p,ex_f,ux_f 129 real*8 ec_p,uc_p,ec_f,uc_f 130 real*8 xe_dirac,fnxup_dirac,fnxdn_dirac 131 real*8 ce_vosko,fncup_vosko,fncdn_vosko 132 133 134 twthrd = thrd*2.0d0 135 frthrd = thrd*4.0d0 136 fvthrd = thrd*5.0d0 137 snthrd = thrd*7.0d0 138 etthrd = thrd*8.0d0 139 140 pi = 4.0d0*datan(1.0d0) 141 Cf = dble(3.0d0*pi*pi) 142 Cf = dble(3.0d0*(Cf**(2.0d0/3.0d0))/10.0d0) 143 144 lda_c =(3.0d0/2.0d0)*(3.0d0/(4.0d0*pi))**(1.0d0/3.0d0) 145 AA = 2.0d0**thrd-1.0d0 146 147 148!$OMP DO 149 do j=1,n2ft3d 150 nup = dn_in(j,1) + 0.5d0*ETA2 151 ndn = dn_in(j,2) + 0.5d0*ETA2 152 153 agrup = agr_in(j,1) 154 agrdn = agr_in(j,2) 155 agrup2 = agrup*agrup 156 agrdn2 = agrdn*agrdn 157 158 n = nup + ndn 159 agr = agr_in(j,3) 160 agr2 = agr*agr 161 162 n2 = n*n 163 nup2 = nup*nup 164 ndn2 = ndn*ndn 165 166* ***** LSDA terms **** 167 x = crs/n**one6th 168 xi = (nup-ndn)/n 169 170 ff = ( (1.0d0+xi)**frthrd 171 > + (1.0d0-xi)**frthrd - 2.0d0)/(2.0d0*AA) 172 dff = frthrd*( (1.0d0+xi)**thrd 173 > - (1.0d0-xi)**thrd) /(2.0d0*AA) 174 175* **** dirac exchange **** 176 ex_p = (xp/x**2) 177 ux_p = frthrd*(xp/x**2) 178 ex_f = (xf/x**2) 179 ux_f = frthrd*(xf/x**2) 180 xe_dirac = ex_p + ff*(ex_f-ex_p) 181 fnxup_dirac = ux_p + ff*(ux_f-ux_p) + (+1.0d0-xi)*dff*(ex_f-ex_p) 182 fnxdn_dirac = ux_p + ff*(ux_f-ux_p) + (-1.0d0-xi)*dff*(ex_f-ex_p) 183 184* **** vosko correlation **** 185 xxp = x**2 + bp*x + cpt 186 dxxp = 2.0d0*x + bp 187 ec_p = cp1*dlog(x**2/xxp) + cp2*dlog( (x+cp3)*(x+cp3)/xxp) 188 > + cp4*datan(cp5/(x+cp6)) 189 uc_p = ec_p 190 > - one6th*x*( dp1/x + dp2/(x+cp3) + dp4*dxxp/xxp 191 > + dp5/( (x+dp6)*(x+dp6)+dp7)) 192 193 xxf = x**2 + bf*x + cft 194 dxxf = 2.0d0*x + bf 195 ec_f = cf1*dlog(x**2/xxf) + cf2*dlog( (x+cf3)*(x+cf3)/xxf) 196 > + cf4*datan(cf5/(x+cf6)) 197 uc_f = ec_f 198 > - one6th*x*( df1/x + df2/(x+cf3) + df4*dxxf/xxf 199 > + df5/( (x+df6)*(x+df6)+df7)) 200 ce_vosko = ec_p + ff*(ec_f - ec_p) 201 fncup_vosko = uc_p + ff*(uc_f-uc_p) + (+1.0d0-xi)*dff*(ec_f-ec_p) 202 fncdn_vosko = uc_p + ff*(uc_f-uc_p) + (-1.0d0-xi)*dff*(ec_f-ec_p) 203 204 if ((dn_in(j,1)+dn_in(j,2)).lt.DNS_CUT) then 205 xe = 0.0d0 206 fdnxup = 0.0d0 207 fdnxdn = 0.0d0 208 fdagrxup = 0.0d0 209 fdagrxdn = 0.0d0 210 else 211 212* *******exchange part*************** 213 214 215* **************UP******************* 216 if (dn_in(j,1).lt.DNS_CUT) then 217 xeup = 0.0d0 218 fdnxup = 0.0d0 219 fdagrxup = 0.0d0 220 else 221 chiup = agrup/nup**(4.0d0/3.0d0) 222 chiup2 = chiup*chiup 223 chiupSQ = dsqrt(1.0d0+chiup2) 224 225 Kup = 6.0d0*beta*dlog(chiup+chiupSQ) 226 F1up = chiup2/(1.0d0 + chiup*Kup) 227 xeup = -nup**thrd*(lda_c + beta*F1up) 228 F2up = (2.0d0 + chiup*Kup 229 & - 6.0d0*beta*chiup2/chiupSQ) 230 & /(1.0d0+chiup*Kup)**2.0d0 231 fdnxup = -nup**(thrd)*(4.0d0/3.0d0) 232 & *(lda_c+beta*(F1up-chiup2*F2up)) 233 fdagrxup = -beta*chiup*F2up 234 if ((fdnxup-xeup).gt.0.0d0) then 235 call gen_PBE96_x_unrestricted(nup,agrup, 236 > xeup_pbe,fdnxup_pbe,fdagrxup_pbe) 237 238 call Becke_smalln_correction(nup,nup**thrd,1.0d0, 239 > beta,lda_c,chiup,chiup2, 240 > chiupSQ,Kup,F1up,F2up, 241 > xeup_pbe,fdnxup_pbe,fdagrxup_pbe, 242 > xeup,fdnxup,fdagrxup) 243 244* call Becke_smalln_correction(nup,nup**thrd,agrup, 245* > beta,lda_c, 246* > chiup,chiup2,chiupSQ,Kup,F1up,F2up, 247* > xeup,fdnxup,fdagrxup) 248 end if 249 end if 250* ************END UP***************** 251 252 253* *************DOWN****************** 254 if (dn_in(j,2).lt.DNS_CUT) then 255 xedn = 0.0d0 256 fdnxdn = 0.0d0 257 fdagrxdn = 0.0d0 258 else 259 chidn = agrdn/ndn**(4.0d0/3.0d0) 260 chidn2 = chidn*chidn 261 chidnSQ = dsqrt(1.0d0+chidn2) 262 263 Kdn = 6.0d0*beta*dlog(chidn+chidnSQ) 264 F1dn = chidn2/(1.0d0 + chidn*Kdn) 265 xedn = -ndn**thrd*(lda_c + beta*F1dn) 266 F2dn = (2.0d0 + chidn*Kdn 267 & - 6.0d0*beta*chidn2/chidnSQ) 268 & /(1.0d0+chidn*Kdn)**2.0d0 269 fdnxdn = -ndn**(thrd)*(4.0d0/3.0d0) 270 & *(lda_c+beta*(F1dn-chidn2*F2dn)) 271 fdagrxdn = -beta*chidn*F2dn 272 if ((fdnxdn-xedn).gt.0.0d0) then 273 call gen_PBE96_x_unrestricted(ndn,agrdn, 274 > xedn_pbe,fdnxdn_pbe,fdagrxdn_pbe) 275 276 call Becke_smalln_correction(ndn,ndn**thrd,1.0d0, 277 > beta,lda_c,chidn,chidn2, 278 > chidnSQ,Kdn,F1dn,F2dn, 279 > xedn_pbe,fdnxdn_pbe,fdagrxdn_pbe, 280 > xedn,fdnxdn,fdagrxdn) 281 282* call Becke_smalln_correction(ndn,ndn**thrd,agrdn, 283* > beta,lda_c, 284* > chidn,chidn2,chidnSQ,Kdn,F1dn,F2dn, 285* > xedn,fdnxdn,fdagrxdn) 286 end if 287 end if 288* ***********END DOWN**************** 289 290 xe = (xeup*nup + xedn*ndn)/n 291 292* *******end excange part************ 293 end if 294 295 296* *******correlation part************ 297 n_thrd = n**thrd 298 n_m = 1.0d0/n 299 n2_m = n_m*n_m 300 n3_m = n2_m*n_m 301 n4_m = n3_m*n_m 302 nupndn = nup*ndn 303 n_mthrd = 1.0d0/n_thrd 304 n_mfrthrd = n_mthrd*n_m 305 n_mfvthrd = n_mfrthrd*n_mthrd 306 307 nup_etthrd = nup**etthrd 308 nup_fvthrd = nup**fvthrd 309 ndn_etthrd = ndn**etthrd 310 ndn_fvthrd = ndn**fvthrd 311 312 313 314 315 gamma = (4.0d0*nupndn)*n2_m 316 gamma_u = 4.0d0*ndn*n2_m - 8.0d0*nupndn*n3_m 317 gamma_d = 4.0d0*nup*n2_m - 8.0d0*nupndn*n3_m 318 319 gamma_uu = -16.0d0*ndn*n3_m + 24.0d0*nupndn*n4_m 320 gamma_dd = -16.0d0*nup*n3_m + 24.0d0*nupndn*n4_m 321 322 gamma_du = (6.0d0*gamma-4.0d0)*n2_m 323 324 325 326 d3 = d/3.0d0 327 n13d_m = 1.0d0/(1.0d0 + d*n_mthrd) 328 n13d2_m = n13d_m*n13d_m 329 n13d3_m = n13d2_m*n13d_m 330 F = gamma*n13d_m 331 F_u = gamma_u*n13d_m + d3*gamma*n_mfrthrd*n13d2_m 332 F_d = gamma_d*n13d_m + d3*gamma*n_mfrthrd*n13d2_m 333 334 F_uu = gamma_uu*n13d_m + d3*(gamma_u+gamma_u)*n_mfrthrd*n13d2_m 335 & - (4.0d0/9.0d0)*d*gamma*n_mfrthrd*n_m*n13d2_m 336 & + (2.0d0/9.0d0)*d*d*gamma*n_mfrthrd*n_mfrthrd*n13d3_m 337 338 F_dd = gamma_dd*n13d_m + d3*(gamma_d+gamma_d)*n_mfrthrd*n13d2_m 339 & - (4.0d0/9.0d0)*d*gamma*n_mfrthrd*n_m*n13d2_m 340 & + (2.0d0/9.0d0)*d*d*gamma*n_mfrthrd*n_mfrthrd*n13d3_m 341 342 F_du = gamma_du*n13d_m + d3*(gamma_u+gamma_d)*n_mfrthrd*n13d2_m 343 & - (4.0d0/9.0d0)*d*gamma*n_mfrthrd*n_m*n13d2_m 344 & + (2.0d0/9.0d0)*d*d*gamma*n_mfrthrd*n_mfrthrd*n13d3_m 345 346 347 enthrd = dexp(-c*n_mthrd) 348 Q = enthrd*n_mfvthrd 349c Q1 = (1.0d0/3.0d0)*c*n_mfrthrd *enthrd 350c & - (5.0d0/3.0d0)*n_mfvthrd*n_m*enthrd 351 Q1 = (Q/3.0d0)*(c*n_mfrthrd - 5.0*n_m) 352 G = F*Q 353 354 P = (1.0d0/3.0d0)*c*n_mfrthrd - (5.0d0/3.0d0)*n_m 355 P1 = ((-4.0d0/9.0d0)*c*n_mfrthrd + (5.d0/3.0d0)*n_m)*n_m 356 357 G_u = F_u*Q + G*P 358 G_d = F_d*Q + G*P 359 360 G_uu = F_uu*Q + F_u*Q1 + G_u*P + G*P1 361 G_dd = F_dd*Q + F_d*Q1 + G_d*P + G*P1 362 G_du = F_du*Q + F_d*Q1 + G_u*P + G*P1 363 364 365 366 fclda = -a*F*n - 2.0d0*a*b*G*Cf*(2.0d0**twthrd) 367 & *(nup_etthrd + ndn_etthrd) 368 369 fclda_u = -a*F_u*n - a*F 370 & - 2.0d0*a*b*G_u*Cf*(2.0d0**twthrd) 371 & *(nup_etthrd + ndn_etthrd) 372 & - 2.0d0*a*b*G*Cf*(2.0d0**twthrd) 373 & *(8.0d0/3.0d0)*nup_fvthrd 374 375 fclda_d = -a*F_d*n - a*F 376 & - 2.0d0*a*b*G_d*Cf*(2.0d0**twthrd) 377 & *(nup_etthrd + ndn_etthrd) 378 & - 2.0d0*a*b*G*Cf*(2.0d0**twthrd) 379 & *(8.0d0/3.0d0)*ndn_fvthrd 380 381 382 383 H0 = (a*b/2.0d0)*(G 384 & + (1.0d0/3.0d0)*(nup*G_d + ndn*G_u) 385 & + (1.0d0/4.0d0)*(nup*G_u + ndn*G_d)) 386 H0_u=(a*b/2.0d0)*(G_u 387 & + (1.0d0/3.0d0)*(G_d + nup*G_du + ndn*G_uu) 388 & + (1.0d0/4.0d0)*(G_u + nup*G_uu + ndn*G_du)) 389 H0_d=(a*b/2.0d0)*(G_d 390 & + (1.0d0/3.0d0)*(G_u + ndn*G_du + nup*G_dd) 391 & + (1.0d0/4.0d0)*(G_d + ndn*G_dd + nup*G_du)) 392 393 394 Hu = (a*b/18.0d0)*(G + (15.0d0/4.0d0)*nup*G_u 395 & - (9.0d0/4.0d0)*ndn*G_d 396 & - (3.0d0)*nup*G_d 397 & + (3.0d0/2.0d0)*ndn*G_u) 398 399 Hu_u = (a*b/18.0d0)*(G_u + (15.0d0/4.0d0)*(G_u + nup*G_uu) 400 & - (9.0d0/4.0d0)*(ndn*G_du) 401 & - (3.0d0)*(G_d + nup*G_du) 402 & + (3.0d0/2.0d0)*(ndn*G_uu)) 403 404 Hu_d = (a*b/18.0d0)*(G_d + (15.0d0/4.0d0)*(nup*G_du) 405 & - (9.0d0/4.0d0)*(G_d+ndn*G_dd) 406 & - (3.0d0)*(nup*G_dd) 407 & + (3.0d0/2.0d0)*(G_u + ndn*G_du)) 408 409 410 411 Hd = (a*b/18.0d0)*(G + (15.0d0/4.0d0)*ndn*G_d 412 & - (9.0d0/4.0d0)*nup*G_u 413 & - (3.0d0)*ndn*G_u 414 & + (3.0d0/2.0d0)*nup*G_d) 415 416 Hd_d = (a*b/18.0d0)*(G_d + (15.0d0/4.0d0)*(G_d + ndn*G_dd) 417 & - (9.0d0/4.0d0)*(nup*G_du) 418 & - (3.0d0)*(G_u + ndn*G_du) 419 & + (3.0d0/2.0d0)*(nup*G_dd)) 420 421 Hd_u = (a*b/18.0d0)*(G_u + (15.0d0/4.0d0)*(ndn*G_du) 422 & - (9.0d0/4.0d0)*(G_u+nup*G_uu) 423 & - (3.0d0)*(ndn*G_uu) 424 & + (3.0d0/2.0d0)*(G_d + nup*G_du)) 425 426 427 fc = fclda + H0*agr2 + Hu*agrup2 + Hd*agrdn2 428 429* ***calculate derivatives w.r.t up and down density 430 fc_u = fclda_u + H0_u*agr2 + Hu_u*agrup2 + Hd_u*agrdn2 431 fc_d = fclda_d + H0_d*agr2 + Hu_d*agrup2 + Hd_d*agrdn2 432 433* ***calculate derivatives w.r.t. up,down and total density gradients 434 fc_agr = 2.0d0*H0*agr 435 fc_agrup = 2.0d0*Hu*agrup 436 fc_agrdn = 2.0d0*Hd*agrdn 437 438 ce = fc/n 439 440 441* *******end correlation part******** 442 443c write(*,*) "n,nup,ndn,agrup,agrdn,:",j,n,nup,ndn,agrup,agrdn,agr 444c write(*,*) "xe,ce :",j,xe,ce 445c write(*,*) "fdnxup,fdncup :",j,fdnxup,fc_u 446c write(*,*) "fdnxdn,fdncdn :",j,fdnxdn,fc_d 447c write(*,*) "fdagrxup,fdargrcup:",j,fdagrxup, 448c > fc_agrup,fc_agr 449c write(*,*) "fdagrxdn,fdargrcdn:",j,fdagrxdn, 450c > fc_agrdn,fc_agr 451c 452c write(*,*) "restricted fdagrc",fc_agr+0.5d0*(fc_agrdn+fc_agrup) 453c write(*,*) "fc_lda,fdnc_lda:",fclda,fclda_u,fclda_d 454c write(*,*) "Hu,Hd,H0:",Hu,Hd,H0 455c write(*,*) "Ho:",H0+0.25d0*Hu+0.25d0*Hd 456c write(*,*) "Ho_n:",0.5d0*(H0_u+H0_d+0.25d0*(Hu_u+Hu_d+Hd_d+Hd_u)) 457c write(*,*) "F,G:",F,G 458c write(*,*) "F_u,G_u:",F_u,G_u 459c write(*,*) "F_d,G_d:",F_d,G_d 460c write(*,*) "G_uu,G_dd,G_du:",G_uu,G_dd,G_du 461c write(*,*) "F_uu,F_dd,F_du:",F_uu,F_dd,F_du 462c write(*,*) "Gamma's:",gamma,gamma_u,gamma_d, 463c > gamma_uu,gamma_dd,gamma_du 464c write(*,*) "n13d",n13d 465c write(*,*) "0.5/nup2,0.5/ndn2:",0.5d0/nup2,0.5d0/ndn2 466c write(*,*) "fc:",fc 467c write(*,*) "0.5d0*(F_uu+F_du)",0.5d0*(F_uu + F_du) 468c write(*,*) "0.5d0*(G_uu+G_du)",0.5d0*(G_uu + G_du) 469c write(*,*) "nup*G_u...",nup*G_u,ndn*G_d,nup*G_d,ndn*G_u 470c write(*,*) 471 472 473**** a*HF + (1-a)*(0.1*Dirac + 0.9*Becke88) + 0.81*LYP + 0.19*Vosko 474**** x_parameter = (1-a) 475 476* ***return blyp exchange correlation values*** 477 xce(j) = x_parameter*0.9d0*xe + c_parameter*0.81d0*ce 478 fn(j,1) = x_parameter*0.9d0*fdnxup + c_parameter*0.81d0*fc_u 479 fn(j,2) = x_parameter*0.9d0*fdnxdn + c_parameter*0.81d0*fc_d 480 481 fdn(j,1) = x_parameter*0.9d0*fdagrxup + c_parameter*0.81*fc_agrup 482 fdn(j,2) = x_parameter*0.9d0*fdagrxdn + c_parameter*0.81*fc_agrdn 483 fdn(j,3) = c_parameter*0.81*fc_agr 484 485 xce(j) = xce(j) + x_parameter*0.10d0*xe_dirac 486 > + c_parameter*0.19d0*ce_vosko 487 fn(j,1) = fn(j,1) + x_parameter*0.10d0*fnxup_dirac 488 > + c_parameter*0.19d0*fncup_vosko 489 fn(j,2) = fn(j,2) + x_parameter*0.10d0*fnxdn_dirac 490 > + c_parameter*0.19d0*fncdn_vosko 491 492 end do 493!$OMP END DO 494 return 495 end 496 497* ************************************************** 498* * * 499* * gen_B3LYP_BW_restricted * 500* * * 501* ************************************************** 502 503* blyp restricted calc. 504* 505* 506* 507 508* subroutine gen_B3LYP_BW_restricted(n2ft3d,rho_in,agr_in,xce,xcp,fn,fdn) 509* input: n2ft3d grid 510* rho_in density 511* agr_in absolute gradient of density 512* x_parameter: scale parameter for exchange 513* c_parameter: scale parameter for correlation 514* output: xce exchange correlation energy density 515* fn d(n*exc)/dn 516* fdn d(n*exc)/d(|grad n|) 517 518 subroutine gen_B3LYP_BW_restricted(n2ft3d,rho_in,agr_in, 519 & x_parameter, 520 & c_parameter,xce,fn,fdn) 521 implicit none 522 523 integer n2ft3d 524 real*8 rho_in(n2ft3d) 525 real*8 agr_in(n2ft3d) 526 real*8 xce(n2ft3d) 527 real*8 fn(n2ft3d) 528 real*8 fdn(n2ft3d) 529 real*8 x_parameter, c_parameter 530 531c*---- parameters given by vosko et al -----------------* 532 real*8 ap,af,x0p,x0f,bp,bf,cpt,cft 533 parameter (ap = 3.109070d-02, af = 1.554530d-02) 534 parameter (x0p=-1.049800d-01, x0f=-3.250000d-01) 535 parameter (bp = 3.727440d+00, bf = 7.060420d+00) 536 parameter (cpt = 1.293520d+01, cft = 1.805780d+01) 537c*------------------------------------------------------* 538 539* constants calculated from vosko's parameters 540 real*8 xp,xf,qp,qf,xx0p,xx0f,fc1,fd1,crs 541 real*8 cp1,cp2,cp3,cp4,cp5,cp6,dp1,dp2,dp3,dp4,dp5,dp6,dp7 542 real*8 cf1,cf2,cf3,cf4,cf5,cf6,df1,df2,df3,df4,df5,df6,df7 543 parameter (xp = -4.581653d-01, xf = -5.772521d-01) 544 parameter (qp = 6.151991d+00, qf = 4.730927d+00) 545 parameter (xx0p= 1.255491d+01, xx0f= 1.586879d+01) 546 parameter (cp1 = 3.109070d-02, cf1 = 1.554530d-02) 547 parameter (cp2 = 9.690228d-04, cf2 = 2.247860d-03) 548 parameter (cp3 = 1.049800d-01, cf3 = 3.250000d-01) 549 parameter (cp4 = 3.878329d-02, cf4 = 5.249122d-02) 550 parameter (cp5 = 3.075995d+00, cf5 = 2.365463d+00) 551 parameter (cp6 = 1.863720d+00, cf6 = 3.530210d+00) 552 parameter (dp1 = 6.218140d-02, df1 = 3.109060d-02) 553 parameter (dp2 = 1.938045d-03, df2 = 4.495720d-03) 554 parameter (dp3 = 1.049800d-01, df3 = 3.250000d-01) 555 parameter (dp4 = -3.205972d-02, df4 = -1.779316d-02) 556 parameter (dp5 = -1.192972d-01, df5 = -1.241661d-01) 557 parameter (dp6 = 1.863720d+00, df6 = 3.530210d+00) 558 parameter (dp7 = 9.461748d+00, df7 = 5.595417d+00) 559 parameter (fc1 = 1.923661d+00, fd1 = 2.564881d+00) 560 parameter (crs = 7.876233d-01) 561 real*8 for3rd,one6th 562 parameter (for3rd=4.0d0/3.0d0,one6th=1.0d0/6.0d0) 563 564 565 566*******local declarations*************************************** 567 integer i 568 real*8 Fc,Gc,C1,C2 569 real*8 n, n_thrd,n_fv,n_fr,n_tw 570 real*8 n_m,n_mthrd,n_mfrthrd,n_mfvthrd 571 real*8 agr,agr2, chi, chi2,chiSQ,sd 572 real*8 K 573 real*8 F1, F2 574 real*8 xe_pbe,fdnx_pbe,fdagrx_pbe 575 real*8 xe,fdnx,fdagrx,xe_dirac,fnx_dirac,x,xx 576 real*8 ce,fdnc, fdagrc,fdnc_lda,ce_vosko,fnc_vosko 577 578 real*8 fc_lda,Ho,Ho_n,Gc_n,Gc_nn,Fc_n,Fc_nn 579 real*8 P,P_n 580 581******* constants ********************************** 582 real*8 pi,thrd,two_thrd 583 parameter (pi=3.14159265358979311599d0) 584 parameter (thrd=1.0d0/3.0d0) 585 parameter (two_thrd=1.25992104989487319066d0) ! two_thrd = 2**thrd 586 587 588*******density cutoff parameters******************** 589 real*8 DNS_CUT, ETA 590 parameter (DNS_CUT = 1.0d-18) 591 parameter (ETA = 1.0d-20) 592 593 real*8 tolrho,minagr 594 parameter (tolrho=2.0e-11,minagr=1.0e-12) 595****** Becke constants ***************************** 596 real*8 lda_c,beta 597 parameter (beta = 0.0042d0) 598 parameter (lda_c = 0.93052573634910018540d0) ! lda_c = (3/2)*(3/(4*pi))**(1/3) 599*******LYP correlation parameters a, b, c, d******** 600 real*8 a,b,c,d,Cf 601 parameter (a = 0.04918d0) 602 parameter (b = 0.132d0) 603 parameter (c = 0.2533d0) 604 parameter (d = 0.349d0) 605 parameter (Cf =2.87123400018819108225d0) ! Cf = (3/10)*(3*pi*pi)**(2/3) 606 607****** collated LYP parameters ********************* 608 real*8 ho1,ho2,ho3,thrd_d,abCf,p2,p3,p4,p5,p6 609 parameter (ho1=(19.0d0/36.0d0)*a*b) 610 parameter (ho2=( 7.0d0/24.0d0)*a*b) 611 parameter (ho3=(59.0d0/72.0d0)*a*b) 612 parameter (thrd_d=thrd*d) 613 parameter (abCf=a*b*Cf) 614 parameter (p2=-5.0d0*thrd,p3=c*thrd) 615 parameter (p4=-4.0d0*thrd*thrd_d) 616 parameter (p5=5.0d0*thrd) 617 parameter (p6=-4.0d0*thrd*thrd*c) 618 619**************************************************** 620 621 622!$OMP DO 623 do i=1,n2ft3d 624 n = rho_in(i) + ETA 625 agr = agr_in(i) 626 n_thrd = n**thrd 627 628* ***** LDA terms **** 629 x = crs/n**one6th 630 631 !**** dirac exchange **** 632 xe_dirac = (xp/x**2) 633 fnx_dirac = for3rd*(xp/x**2) 634 635 !**** vosko correlation **** 636 xx=1.0d0/(x*(x+bp)+cpt) 637 ce_vosko=cp1*dlog(xx*x**2)+cp2*dlog(xx*(x+cp3)**2) 638 > +cp4*datan(cp5/(x+cp6)) 639 fnc_vosko = ce_vosko-one6th*x*( 640 > dp1/x+dp2/(x+dp3)+dp4*xx*(2.0d0*x+bp) 641 > +dp5/((x+dp6)**2+dp7) ) 642 643 644 if (rho_in(i).lt.DNS_CUT) then 645 xe = 0.0d0 646 fdnx = 0.0d0 647 fdagrx = 0.0d0 648 else 649 650****************************************************************** 651* *******calc. becke exchange energy density, fnx, fdnx******* 652***************************************************************** 653 sd = 1.0d0/(n_thrd*n) 654 chi = two_thrd*agr*sd 655 chi2 = chi*chi 656 chiSQ = dsqrt(1.0d0+chi2) 657 658 K = 6.0d0*beta*dlog(chi+chiSQ) 659 F1 = chi2/(1.0d0+chi*K) 660 xe = -n_thrd*(lda_c+beta*F1)/two_thrd 661 F2 = (2.0d0 + chi*K-(chi2)*6.0d0*beta 662 & /chiSQ) 663 & /((1.0d0+chi*K)*(1.0d0+chi*K)) 664 fdnx = -(n_thrd/two_thrd)*dble(4.0d0/3.0d0) 665 & *(lda_c+beta*(F1-chi2*F2)) 666 fdagrx = -beta*chi*F2 667 if ((fdnx-xe).gt.0.0d0) then 668 call gen_PBE96_x_restricted(n,agr, 669 > xe_pbe,fdnx_pbe,fdagrx_pbe) 670 call Becke_smalln_correction(n,n**thrd,2.0d0,beta, 671 > lda_c,chi,chi2,chiSQ,K,F1,F2, 672 > xe_pbe,fdnx_pbe,fdagrx_pbe, 673 > xe,fdnx,fdagrx) 674 675* call Becke_smalln_correction(n,n_thrd,agr,beta,lda_c, 676* > chi,chi2,chiSQ,K,F1,F2, 677* > xe,fdnx,fdagrx) 678 end if 679 end if 680 681 682*******final result for restricted LYP***************************** 683 agr2 = agr*agr 684 685 n_m = 1.0d0/n 686 n_mthrd = 1.0d0/n_thrd 687 n_mfrthrd = n_mthrd*n_m 688 n_mfvthrd = n_mfrthrd*n_mthrd 689 690 n_fv = n_thrd*n_thrd*n_thrd*n_thrd*n_thrd 691 n_fr = n_thrd*n_thrd*n_thrd*n_thrd 692 n_tw = n_thrd*n_thrd 693 694 695 Fc = (1.0d0/(1.0d0+d*n_mthrd)) 696 Fc_n = thrd_d*n_mfrthrd*Fc*Fc 697 Fc_nn = thrd_d*(-4.0d0*thrd)*n_mfrthrd*n_m*Fc*Fc 698 > + thrd_d*n_mfrthrd*2.0d0*Fc*Fc_n 699 700 Gc = Fc*dexp(-c*n_mthrd)*n_mfvthrd 701 702 P = (thrd_d*Fc*n_mfrthrd - 5.0d0*thrd*n_m + c*thrd*n_mfrthrd) 703 704 P_n = thrd_d*Fc_n*n_mfrthrd 705 > - thrd_d*Fc *n_mfrthrd*(4.0d0*thrd*n_m) 706 > + 5.0d0*thrd*n_m*n_m 707 > - 4.0d0*thrd*thrd*c*n_mfrthrd*n_m 708 709c P = (thrd_d*Fc*n_mfrthrd + p2*n_m + p3*n_mfrthrd) 710c P_n = thrd_d*Fc_n*n_mfrthrd 711c > + p4*Fc*n_mfrthrd*n_m 712c > + p5*n_m*n_m 713c > + p6*n_mfrthrd*n_m 714 715 Gc_n = Gc*P 716 Gc_nn = Gc*P*P + Gc*P_n 717 718 719 fc_lda = -a*Fc*n 720 > - abCf*Gc*(n_fr*n_fr) 721 fdnc_lda = -a*Fc_n*n 722 > - a*Fc 723 > - abCf*Gc_n*n_fr*n_fr 724 > - 8.0d0*thrd*abCf*Gc*n_fv 725 726c Ho = (19.0d0/36.0d0)*(a*b*Gc) + (7.0d0/24.0d0)*a*b*Gc_n*n 727c Ho_n = (59.0d0/72.0d0)*(a*b*Gc_n) + (7.0d0/24.0d0)*a*b*Gc_nn*n 728 Ho = ho1*Gc + ho2*Gc_n*n 729 Ho_n = ho3*Gc_n + ho2*Gc_nn*n 730 731 ce = (fc_lda + Ho*agr2)*n_m 732 fdnc = fdnc_lda + Ho_n*agr2 733 734 fdagrc = 2.0d0*Ho*agr 735 736c write(*,*) "n,agr :",i,n,agr 737c write(*,*) "xe,ce :",i,xe,ce 738c write(*,*) "fdnx,fdnc :",i,fdnx,fdnc 739c write(*,*) "fdagrx,fdargrc:",i,fdagrx,fdagrc 740c write(*,*) "fc_lda,fdnc_lda:",fc_lda,fdnc_lda 741c write(*,*) "Ho :",Ho 742c write(*,*) "Ho_n:",Ho_n 743c write(*,*) "F,G:",Fc,Gc 744c write(*,*) "F_n,G_n:",Fc_n,Gc_n 745c write(*,*) "G_nn:",Gc_nn 746c write(*,*) "F_nn:",Fc_nn 747c write(*,*) "n13d:",1.0d0+d*n_mthrd 748c write(*,*) "fc:",fc_lda + Ho*agr2 749c write(*,*) 750 751**** a*HF + (1-a)*(0.1*Dirac + 0.9*Becke88) + 0.81*LYP + 0.19*Vosko 752**** x_parameter = (1-a) 753 754 xce(i) = x_parameter*0.9d0*xe + c_parameter*0.81d0*ce 755 fn(i) = x_parameter*0.9d0*fdnx + c_parameter*0.81d0*fdnc 756 fdn(i) = x_parameter*0.9d0*fdagrx + c_parameter*0.81d0*fdagrc 757 758 xce(i) = xce(i) + x_parameter*0.10d0*xe_dirac 759 > + c_parameter*0.19d0*ce_vosko 760 fn(i) = fn(i) + x_parameter*0.10d0*fnx_dirac 761 > + c_parameter*0.19d0*fnc_vosko 762 end do 763!$OMP END DO 764 return 765 end 766