1* ************************************************ 2* * * 3* * gen_M06L_restricted * 4* * * 5* ************************************************ 6 7* This function returns the M06-L exchange-correlation 8* energy density, xce, and its derivatives with respect 9* to n, |grad n|, tau. 10 11* 12* Entry - n2ft3d : number of grid points 13* rho_in(*) : density (nup+ndn) 14* agr_in(*): |grad rho_in| 15* tau_in(*): tau 16* x_parameter: scale parameter for exchange 17* c_parameter: scale parameter for correlation 18* 19* Exit - xce(n2ft3d) : M06-L exchange correlation energy density 20* fn(n2ft3d) : d(n*xce)/dn 21* fdn(n2ft3d) : d(n*xce)/d|grad n| 22* fdtau(n2ft3d) : d(n*xce)/dtau 23* 24 subroutine gen_M06L_restricted(n2ft3d,rho_in,agr_in,tau_in, 25 > x_parameter,c_parameter, 26 > xce,fn,fdn,fdtau) 27 implicit none 28* ***** input ***** 29 integer n2ft3d 30 real*8 rho_in(*),agr_in(*),tau_in(*) 31 real*8 x_parameter,c_parameter 32* ***** output ***** 33 real*8 xce(*),fn(*),fdn(*),fdtau(*) 34* ***** local declarations ***** 35 integer i 36 real*8 Cx,P23 37 real*8 nup,agrup,tauup 38 real*8 n,agr,tau 39 real*8 xr,zr,gamma 40 real*8 inv_n,n_13,n_43,n_53,n_83,agr2 41 real*8 x,dx_dn,dx_dagr,z,dz_dn,dz_dtau 42 real*8 GG,dGdx,dGdz,dG_dn,dG_dagr,dG_dtau 43 real*8 n_onethird,kf,ks,ss,P0,FF,Fs,ex_lda,fdnx_const 44 real*8 pbex,pbefnx,pbefdnx 45 real*8 tauU,t,dt_dn,dt_dtau 46 real*8 w,dw_dt,w1,w2,w3,w4,w5,w6,w7,w8,w9,w10,w11 47 real*8 fw,dfw_dw,dfw_dn,dfw_dtau 48 real*8 ex,fnx,fdnx,fdtaux 49 real*8 ess0c,dess0c_drs,dess0c_dn 50 real*8 U,dU_dx,U2,U3,U4 51 real*8 gss,dgss_dU,dgss_dx,gopp,dgopp_dU,dgopp_dx 52 real*8 hss,dhss_dx,dhss_dz,hopp,dhopp_dx,dhopp_dz 53 real*8 ghss,ghopp,dgh_dx,dgh_dn,dgh_dagr,dgh_dtau 54 real*8 eud0c,deud0c_drs 55 real*8 eudc,eud1c,fud1c,dfud1c_dn 56 real*8 dfudc_dn,dfudc_dagr,dfudc_dtau 57 real*8 rs,drs_dn,dummy 58 real*8 wt1,dwt1,aw,bw,cw,dw 59 real*8 D,dD_dx,dD_dz,dD_dn,dD_dagr,dD_dtau 60 real*8 Dt,dD_dDt,dDt_dx,dDt_dz 61 real*8 ec,fnc,fdnc,fdtauc 62* ***** constants ***** 63 real*8 pi,cf,thrd,twthrd,frthrd,fvthrd,etthrd 64 parameter (pi = 3.14159265358979311599d0) 65 parameter (cf = 9.115599720d0) 66c cf = 3.0d0/5.0d0*(6.0d0*pi*pi)**(2.0d0/3.0d0)) 67 parameter (thrd = 1.0d0/3.0d0) 68 parameter (twthrd = 2.0d0/3.0d0) 69 parameter (frthrd = 4.0d0/3.0d0) 70 parameter (fvthrd = 5.0d0/3.0d0) 71 parameter (etthrd = 8.0d0/3.0d0) 72* ***** density cutoff parametersi ***** 73 real*8 tol,ETA,t1,t2,thresd,thresx 74 parameter (tol = 1.0d-10) 75 parameter (ETA = 1.0d-20) 76 parameter (t1 = 1.0d7) 77 parameter (t2 = 5.0d7) 78 parameter (thresd = 5.0d-3) 79 parameter (thresx = 1.0d8) 80c ***** PBE96 GGA exchange constants ****** 81 real*8 MU,KAPPA 82 parameter (MU = 0.2195149727645171d0) 83 parameter (KAPPA = 0.8040000000000000d0) 84* ***** VS98 constants ***** 85 real*8 aax,bbx,ccx,ddx,eex,ffx,alphax 86 real*8 aass,bbss,ccss,ddss,eess,ffss,alphass 87 real*8 aaopp,bbopp,ccopp,ddopp,eeopp,ffopp,alphaopp 88 real*8 clda 89 parameter (clda = 0.9305257363491d0) 90 parameter (aax = 6.012244d-1) 91 parameter (bbx = 4.748822d-3) 92 parameter (ccx = -8.635108d-3) 93 parameter (ddx = -9.308062d-6) 94 parameter (eex = 4.482811d-5) 95 parameter (ffx = 0.000000d0) 96 parameter (alphax = 0.00186726d0) 97 parameter (aass = 4.650534d-1) 98 parameter (bbss = 1.617589d-1) 99 parameter (ccss = 1.833657d-1) 100 parameter (ddss = 4.692100d-4) 101 parameter (eess = -4.990573d-3) 102 parameter (ffss = 0.000000d0) 103 parameter (alphass = 0.00515088d0) 104 parameter (aaopp = 3.957626d-1) 105 parameter (bbopp = -5.614546d-1) 106 parameter (ccopp = 1.403963d-2) 107 parameter (ddopp = 9.831442d-4) 108 parameter (eeopp = -3.577176d-3) 109 parameter (ffopp = 0.000000d0) 110 parameter (alphaopp = 0.00304966d0) 111* ***** M06 constants ***** 112 real*8 ma0,ma1,ma2,ma3,ma4,ma5,ma6,ma7,ma8,ma9,ma10,ma11 113 real*8 C1,C2 114 real*8 css,ccss0,ccss1,ccss2,ccss3,ccss4 115 real*8 copp,ccopp0,ccopp1,ccopp2,ccopp3,ccopp4 116 parameter (ma0 = 3.987756d-1) 117 parameter (ma1 = 2.548219d-1) 118 parameter (ma2 = 3.923994d-1) 119 parameter (ma3 = -2.103655d0) 120 parameter (ma4 = -6.302147d0) 121 parameter (ma5 = 1.097615d1) 122 parameter (ma6 = 3.097273d1) 123 parameter (ma7 = -2.318489d1) 124 parameter (ma8 = -5.673480d1) 125 parameter (ma9 = 2.160364d1) 126 parameter (ma10 = 3.421814d1) 127 parameter (ma11 = -9.049762d0) 128 parameter (C1 = 3.36116d-3) 129 parameter (C2 = 4.49267d-3) 130 parameter (css = 0.06d0) 131 parameter (ccss0 = 5.349466d-1) 132 parameter (ccss1 = 5.396620d-1) 133 parameter (ccss2 = -3.161217d1) 134 parameter (ccss3 = 5.149592d1) 135 parameter (ccss4 = -2.919613d1) 136 parameter (copp = 0.0031d0) 137 parameter (ccopp0 = 6.042374d-1) 138 parameter (ccopp1 = 1.776783d2) 139 parameter (ccopp2 = -2.513252d2) 140 parameter (ccopp3 = 7.635173d1) 141 parameter (ccopp4 = -1.255699d1) 142 143 Cx = -1.50d0*(0.75d0/pi)**thrd 144 P23 = 0.60d0*(6.0d0*pi*pi)**twthrd 145 fdnx_const = -3.0d0/(8.0d0*pi) 146 147 do i=1,n2ft3d 148 n = rho_in(i) + ETA 149 agr = agr_in(i) + ETA 150 tau = 2.0d0*tau_in(i) + ETA 151 152 n = 0.50d0*n 153 agr = 0.50d0*agr 154 tau = 0.50d0*tau 155 156 agr2 = agr*agr 157 inv_n = 1.0d0/n 158 n_13 = n**thrd 159 n_43 = n_13*n 160 n_53 = n_43*n_13 161 n_83 = n_53*n 162 163 x = agr2/n_83 164 dx_dn = -etthrd*x*inv_n 165 dx_dagr = 2.0d0*agr/n_83 166 z = tau/n_53 - cf 167 dz_dn = -fvthrd*tau/n_83 168 dz_dtau = 1.0d0/n_53 169 170* ***** VS98 Exchange ***** 171 gamma = 1.0d0 + alphax*(x + z) 172 xr = x/gamma 173 zr = z/gamma 174 175 call nwpw_GVT4(aax,bbx,ccx,ddx,eex,ffx,alphax,alphax, 176 > xr,zr,gamma,GG,dGdx,dGdz) 177 178 dG_dn = dGdx*dx_dn + dGdz*dz_dn 179 dG_dagr = dGdx*dx_dagr 180 dG_dtau = dGdz*dz_dtau 181 182 ex = -clda*n_13*GG 183 fnx = frthrd*ex - clda*n_43*dG_dn 184 fdnx = -clda*n_43*dG_dagr 185 fdtaux = -clda*n_43*dG_dtau 186 187* ***** M06-L Enhancement Factor ***** 188 tauU = P23*n_53 189 t = tauU/tau 190 dt_dn = fvthrd*t/n 191 dt_dtau = -t/tau 192 193 w = 1.0d0 194 dw_dt = 0.0d0 195 196 if(t .le. t1) then 197 w = (t - 1.0d0)/(t + 1.0d0) 198 dw_dt = 2.0d0/((1.0d0 + t)**2.0d0) 199 else if(t .gt. t1 .and. t .lt. t2) then 200 wt1 = (t1 - 1.0d0)/(t1 + 1.0d0) 201 dwt1 = 2.0d0/((1.0d0 + t1)**2.0d0) 202 aw = (3.0d0*t1*t2*t2 - t2**3.0d0)*wt1/(t1 - t2) 203 & + (t1**3.0d0 - 3.0d0*t1*t1*t2)/(t1 - t2) 204 & - t1*t2*t2*dwt1 205 aw = aw/((t1 - t2)*(t1 - t2)) 206 bw = -6.0d0*t1*t2*wt1/(t1 - t2) + 6.0d0*t1*t2/(t1 - t2) 207 & + (t2*t2 + 2.0d0*t1*t2)*dwt1 208 bw = bw/((t1 - t2)*(t1 - t2)) 209 cw = 3.0d0*(t1 + t2)*wt1/(t1 - t2) 210 & - 3.0d0*(t1 + t2)/(t1 - t2) - (t1 + 2.0d0*t2)*dwt1 211 cw = cw/((t1 - t2)*(t1 - t2)) 212 dw = -2.0d0*wt1/(t1 - t2) + 2.0d0/(t1 - t2) + dwt1 213 dw = dw/((t1 - t2)*(t1 - t2)) 214 215 w = aw + bw*t + cw*t*t + dw*t*t*t 216 dw_dt = bw + 2.0d0*cw*t + 3.0d0*dw*t*t 217 else if(t .ge. t2) then 218 w = 1.0d0 219 dw_dt = 0.0d0 220 end if 221 222 w1 = w 223 w2 = w1*w 224 w3 = w2*w 225 w4 = w3*w 226 w5 = w4*w 227 w6 = w5*w 228 w7 = w6*w 229 w8 = w7*w 230 w9 = w8*w 231 w10 = w9*w 232 w11 = w10*w 233 234 fw = ma0 + ma1*w1 + ma2*w2 + ma3*w3 + ma4*w4 235 & + ma5*w5 + ma6*w6 + ma7*w7 + ma8*w8 + ma9*w9 236 & + ma10*w10 + ma11*w11 237 dfw_dw = ma1 + 2.0d0*ma2*w1 + 3.0d0*ma3*w2 238 & + 4.0d0*ma4*w3 + 5.0d0*ma5*w4 + 6.0d0*ma6*w5 239 & + 7.0d0*ma7*w6 + 8.0d0*ma8*w7 + 9.0d0*ma9*w8 240 & + 10.0d0*ma10*w9 + 11.0d0*ma11*w10 241 dfw_dn = dfw_dw*dw_dt*dt_dn 242 dfw_dtau = dfw_dw*dw_dt*dt_dtau 243 244* ***** PBE96 Exchange ***** 245 n_onethird = (3.0d0*n/pi)**thrd 246 ex_lda = -0.75d0*n_onethird 247 248 kf = (3.0d0*pi*pi*n)**thrd 249 ss = agr/(2.0d0*kf*n) 250 P0 = 1.0d0 + (MU/KAPPA)*ss*ss 251 252 FF = (1.0d0 + KAPPA - KAPPA/P0) 253 Fs = 2.0d0*MU/(P0*P0)*ss 254 255 pbex = ex_lda*FF 256 pbefnx = frthrd*(pbex - ex_lda*Fs*ss) 257 pbefdnx = fdnx_const*Fs 258 259 ex = ex + pbex*fw 260 fnx = fnx + pbefnx*fw + n*pbex*dfw_dn 261 fdnx = fdnx + pbefdnx*fw 262 fdtaux = fdtaux + n*pbex*dfw_dtau 263 264 265* ***** VS98 Correlation ***** 266* ***** Same-Spin (alpha-alpha/beta-beta) ***** 267 rs = (0.75d0/(pi*n))**thrd 268 drs_dn = -thrd*rs/n 269 270 call gen_PW91_c_rz(tol,rs,1.0d0,ess0c,dess0c_drs,dummy) 271 dess0c_dn = dess0c_drs*drs_dn 272 273 Dt = 1.0d0 - 0.25d0*x/(z + cf) 274 D = 0.0d0 275 dD_dn = 0.0d0 276 dD_dagr = 0.0d0 277 dD_dtau = 0.0d0 278 279 if (Dt .le. 0.0d0) then 280 D = 0.0d0 281 dD_dn = 0.0d0 282 dD_dagr = 0.0d0 283 dD_dtau = 0.0d0 284 else if (Dt .gt. 0.0d0 .and. Dt .lt. thresd) then 285 D = 2.0d0*Dt*Dt/thresd - Dt*Dt*Dt/(thresd*thresd) 286 dD_dDt = 4.0d0*Dt/thresd - 3.0d0*Dt*Dt/(thresd*thresd) 287 dDt_dx = -0.25d0/(z + cf) 288 dDt_dz = 0.25d0*x/((z + cf)*(z + cf)) 289 dD_dn = dD_dDt*(dDt_dx*dx_dn + dDt_dz*dz_dn) 290 dD_dagr = dD_dDt*dDt_dx*dx_dagr 291 dD_dtau = dD_dDt*dDt_dz*dz_dtau 292 else if(Dt .ge. thresd) then 293 D = Dt 294 dD_dx = -0.25d0/(z + cf) 295 dD_dz = 0.25d0*x/((z + cf)*(z + cf)) 296 dD_dn = dD_dx*dx_dn + dD_dz*dz_dn 297 dD_dagr = dD_dx*dx_dagr 298 dD_dtau = dD_dz*dz_dtau 299 end if 300 301 gamma = 1.0d0 + alphass*(x + z) 302 xr = x/gamma 303 zr = z/gamma 304 305 if (x .ge. thresx) then 306 U = 1.0d0 307 dU_dx = 0.0d0 308 else 309 U = css*x/(1.0d0 + css*x) 310 dU_dx = css/((1.0d0 + css*x)*(1.0d0 + css*x)) 311 end if 312 313 U2 = U*U 314 U3 = U2*U 315 U4 = U3*U 316 317 gss = ccss0 + ccss1*U + ccss2*U2 + ccss3*U3 + ccss4*U4 318 dgss_dU = ccss1 + 2.0d0*ccss2*U + 3.0d0*ccss3*U2 319 & + 4.0d0*ccss4*U3 320 dgss_dx = dgss_dU*dU_dx 321 322 call nwpw_GVT4(aass,bbss,ccss,ddss,eess,ffss,alphass,alphass, 323 > xr,zr,gamma,hss,dhss_dx,dhss_dz) 324 325 ghss = gss + hss 326 dgh_dx = dgss_dx + dhss_dx 327 dgh_dn = dgh_dx*dx_dn + dhss_dz*dz_dn 328 dgh_dagr = dgh_dx*dx_dagr 329 dgh_dtau = dhss_dz*dz_dtau 330 331 ec = ess0c*ghss*D 332 fnc = ec + n*dess0c_drs*drs_dn*ghss*D 333 & + n*ess0c*dgh_dn*D 334 & + n*ess0c*ghss*dD_dn 335 fdnc = n*ess0c*(dgh_dagr*D + ghss*dD_dagr) 336 fdtauc = n*ess0c*(dgh_dtau*D + ghss*dD_dtau) 337 338* ***** Opposite-Spin (alpha-beta) ***** 339 n = 2.0d0*n 340 341 rs = (0.75d0/(pi*n))**thrd 342 drs_dn = -thrd*rs/n 343 344 call gen_PW91_c_rz(tol,rs,0.0d0,eud0c,deud0c_drs,dummy) 345 eud1c = eud0c - ess0c 346 fud1c = n*eud1c 347 dfud1c_dn = eud0c + n*deud0c_drs*drs_dn 348 & - ess0c - 0.50d0*n*dess0c_dn 349 350 x = 2.0d0*x 351 z = 2.0d0*z 352 353 gamma = 1.0d0 + alphaopp*(x + z) 354 xr = x/gamma 355 zr = z/gamma 356 357 if (x .ge. thresx) then 358 U = 1.0d0 359 dU_dx = 0.0d0 360 else 361 U = copp*x/(1.0d0 + copp*x) 362 dU_dx = copp/((1.0d0 + copp*x)*(1.0d0 + copp*x)) 363 end if 364 365 U2 = U*U 366 U3 = U2*U 367 U4 = U3*U 368 369 gopp = ccopp0 + ccopp1*U + ccopp2*U2 + ccopp3*U3 + ccopp4*U4 370 dgopp_dU = ccopp1 + 2.0d0*ccopp2*U 371 & + 3.0d0*ccopp3*U2 + 4.0d0*ccopp4*U3 372 dgopp_dx = dgopp_dU*dU_dx 373 374 call nwpw_GVT4(aaopp,bbopp,ccopp,ddopp,eeopp,ffopp, 375 > alphaopp,alphaopp, 376 > xr,zr,gamma,hopp,dhopp_dx,dhopp_dz) 377 378 ghopp = gopp + hopp 379 dgh_dx = dgopp_dx + dhopp_dx 380 dgh_dn = dgh_dx*dx_dn + dhopp_dz*dz_dn 381 dgh_dagr = dgh_dx*dx_dagr 382 dgh_dtau = dhopp_dz*dz_dtau 383 384 eudc = eud1c*ghopp 385 dfudc_dn = dfud1c_dn*ghopp + fud1c*dgh_dn 386 dfudc_dagr = fud1c*dgh_dagr 387 dfudc_dtau = fud1c*dgh_dtau 388 389 ec = ec + eudc 390 fnc = fnc + dfudc_dn 391 fdnc = fdnc + dfudc_dagr 392 fdtauc = fdtauc + dfudc_dtau 393 394 xce(i) = x_parameter*ex + c_parameter*ec 395 fn(i) = x_parameter*fnx + c_parameter*fnc 396 fdn(i) = x_parameter*fdnx + c_parameter*fdnc 397 fdtau(i) = x_parameter*fdtaux + c_parameter*fdtauc 398 end do 399 400 return 401 end 402* ************************************************ 403* * * 404* * gen_M06L_unrestricted * 405* * * 406* ************************************************ 407 408* This function returns the M06-L exchange-correlation 409* energy density, xce, and its derivatives with respect 410* to nup, ndn, |grad nup|, |grad ndn|, tauup, taudn. 411 412* 413* Entry - n2ft3d : number of grid points 414* rho_in(*,2) : density (nup and ndn) 415* agr_in(*,3): |grad rho_in| (nup, ndn and n) 416* tau_in(*,2): tau (nup and ndn) 417* x_parameter: scale parameter for exchange 418* c_parameter: scale parameter for correlation 419* 420* Exit - xce(n2ft3d) : M06-L exchange correlation energy density 421* fn(n2ft3d,2) : d(n*xce)/dnup, d(n*xce)/dndn 422* fdn(n2ft3d,3) : d(n*xce)/d|grad nup|, d(n*xce)/d|grad ndn|, d(n*xce)/d|grad n| 423* fdtau(n2ft3d,2) : d(n*xce)/dtauup, d(n*xce)/dtaudn 424* 425 subroutine gen_M06L_unrestricted(n2ft3d,rho_in,agr_in,tau_in, 426 > x_parameter,c_parameter, 427 > xce,fn,fdn,fdtau) 428 implicit none 429* ***** input ***** 430 integer n2ft3d 431 real*8 rho_in(n2ft3d,2),agr_in(n2ft3d,3),tau_in(n2ft3d,2) 432 real*8 x_parameter,c_parameter 433* ***** output ***** 434 real*8 xce(n2ft3d),fn(n2ft3d,2),fdn(n2ft3d,3),fdtau(n2ft3d,2) 435* ***** local declarations ***** 436 integer i 437 real*8 Cx,P23 438 real*8 nup,agrup,tauup 439 real*8 ndn,agrdn,taudn 440 real*8 n,agr,tau 441 real*8 x,dxdn,dxdagr,z,dzdn,dzdtau 442 real*8 xr,zr,gamma 443 real*8 inv_nup,nup_13,nup_43,nup_53,nup_83,agrup2 444 real*8 inv_ndn,ndn_13,ndn_43,ndn_53,ndn_83,agrdn2 445 real*8 xup,dxup_dnup,dxup_dagrup,zup,dzup_dnup,dzup_dtauup 446 real*8 xdn,dxdn_dndn,dxdn_dagrdn,zdn,dzdn_dndn,dzdn_dtaudn 447 real*8 GG,dGdx,dGdz 448 real*8 n_onethird,kf,ks,ss,P0,FF,Fs,ex_lda,fdnx_const 449 real*8 pbex,pbefnx,pbefdnx 450 real*8 tauU,t,dt_dn,dt_dtau 451 real*8 w,dw_dt,w1,w2,w3,w4,w5,w6,w7,w8,w9,w10,w11 452 real*8 fw,dfw_dw,dfw_dn,dfw_dtau 453 real*8 eupx,fnupx,fdnupx,fdtauupx 454 real*8 ednx,fndnx,fdndnx,fdtaudnx 455 real*8 euu0c,deuu0c_dnup 456 real*8 U,dU_dx,U2,U3,U4,gss,dgss_dU 457 real*8 hss,dhss_dx,dhss_dz,ghss,dgh_dx,dgh_dn,dgh_dagr,dgh_dtau 458 real*8 euuc,dfuuc_dnup,dfuuc_dagrup,dfuuc_dtauup 459 real*8 edd0c,dedd0c_dndn 460 real*8 eddc,dfddc_dndn,dfddc_dagrdn,dfddc_dtaudn 461 real*8 eud0c,deud0c_drs,deud0c_dzeta 462 real*8 eudc,eud1c,fud1c,dfud1c_dnup,dfud1c_dndn 463 real*8 gopp,dgopp_dU,hopp,dhopp_dx,dhopp_dz 464 real*8 ghopp,dgh_dnup,dgh_dndn 465 real*8 dgh_dagrup,dgh_dagrdn,dgh_dtauup,dgh_dtaudn 466 real*8 dfudc_dnup,dfudc_dagrup,dfudc_dtauup 467 real*8 dfudc_dndn,dfudc_dagrdn,dfudc_dtaudn 468 real*8 rs,drs_dn,zeta,dzeta_dnup,dzeta_dndn,deuu0c_drs,dedd0c_drs 469 real*8 wt1,dwt1,aw,bw,cw,dw 470 real*8 D,dD_dx,dD_dz,dD_dn,dD_dagr,dD_dtau,dummy 471 real*8 Dt,dD_dDt,dDt_dx,dDt_dz 472 real*8 eupc,fnupc,fdnupc,fdtauupc 473 real*8 ednc,fndnc,fdndnc,fdtaudnc 474 real*8 ex,ec 475 real*8 xbar,xb2,xb3,xb4,xb5,dgss_dx,dgopp_dx 476* ***** constants ***** 477 real*8 pi,cf,thrd,twthrd,frthrd,fvthrd,etthrd 478 parameter (pi = 3.14159265358979311599d0) 479 parameter (cf = 9.115599720d0) 480c cf = 3.0d0/5.0d0*(6.0d0*pi*pi)**(2.0d0/3.0d0)) 481 parameter (thrd = 1.0d0/3.0d0) 482 parameter (twthrd = 2.0d0/3.0d0) 483 parameter (frthrd = 4.0d0/3.0d0) 484 parameter (fvthrd = 5.0d0/3.0d0) 485 parameter (etthrd = 8.0d0/3.0d0) 486* ***** density cutoff parametersi ***** 487 real*8 tol,ETA,t1,t2,thresd,thresx 488 parameter (tol = 1.0d-10) 489 parameter (ETA = 1.0d-20) 490 parameter (t1 = 1.0d7) 491 parameter (t2 = 5.0d7) 492 parameter (thresd = 5.0d-3) 493 parameter (thresx = 1.0d8) 494c ***** PBE96 GGA exchange constants ****** 495 real*8 MU,KAPPA 496 parameter (MU = 0.2195149727645171d0) 497 parameter (KAPPA = 0.8040000000000000d0) 498* ***** VS98 constants ***** 499 real*8 aax,bbx,ccx,ddx,eex,ffx,alphax 500 real*8 aass,bbss,ccss,ddss,eess,ffss,alphass 501 real*8 aaopp,bbopp,ccopp,ddopp,eeopp,ffopp,alphaopp 502 real*8 clda 503 parameter (clda = 0.9305257363491d0) 504 parameter (aax = 6.012244d-1) 505 parameter (bbx = 4.748822d-3) 506 parameter (ccx = -8.635108d-3) 507 parameter (ddx = -9.308062d-6) 508 parameter (eex = 4.482811d-5) 509 parameter (ffx = 0.000000d0) 510 parameter (alphax = 0.00186726d0) 511 parameter (aass = 4.650534d-1) 512 parameter (bbss = 1.617589d-1) 513 parameter (ccss = 1.833657d-1) 514 parameter (ddss = 4.692100d-4) 515 parameter (eess = -4.990573d-3) 516 parameter (ffss = 0.000000d0) 517 parameter (alphass = 0.00515088d0) 518 parameter (aaopp = 3.957626d-1) 519 parameter (bbopp = -5.614546d-1) 520 parameter (ccopp = 1.403963d-2) 521 parameter (ddopp = 9.831442d-4) 522 parameter (eeopp = -3.577176d-3) 523 parameter (ffopp = 0.000000d0) 524 parameter (alphaopp = 0.00304966d0) 525* ***** M06 constants ***** 526 real*8 ma0,ma1,ma2,ma3,ma4,ma5,ma6,ma7,ma8,ma9,ma10,ma11 527 real*8 C1,C2 528 real*8 css,ccss0,ccss1,ccss2,ccss3,ccss4 529 real*8 copp,ccopp0,ccopp1,ccopp2,ccopp3,ccopp4 530 parameter (ma0 = 3.987756d-1) 531 parameter (ma1 = 2.548219d-1) 532 parameter (ma2 = 3.923994d-1) 533 parameter (ma3 = -2.103655d0) 534 parameter (ma4 = -6.302147d0) 535 parameter (ma5 = 1.097615d1) 536 parameter (ma6 = 3.097273d1) 537 parameter (ma7 = -2.318489d1) 538 parameter (ma8 = -5.673480d1) 539 parameter (ma9 = 2.160364d1) 540 parameter (ma10 = 3.421814d1) 541 parameter (ma11 = -9.049762d0) 542 parameter (C1 = 3.36116d-3) 543 parameter (C2 = 4.49267d-3) 544 parameter (css = 0.06d0) 545 parameter (ccss0 = 5.349466d-1) 546 parameter (ccss1 = 5.396620d-1) 547 parameter (ccss2 = -3.161217d1) 548 parameter (ccss3 = 5.149592d1) 549 parameter (ccss4 = -2.919613d1) 550 parameter (copp = 0.0031d0) 551 parameter (ccopp0 = 6.042374d-1) 552 parameter (ccopp1 = 1.776783d2) 553 parameter (ccopp2 = -2.513252d2) 554 parameter (ccopp3 = 7.635173d1) 555 parameter (ccopp4 = -1.255699d1) 556 557 Cx = -1.50d0*(0.75d0/pi)**thrd 558 P23 = 0.60d0*(6.0d0*pi*pi)**twthrd 559 fdnx_const = -3.0d0/(8.0d0*pi) 560 561 do i=1,n2ft3d 562 nup = rho_in(i,1) + ETA 563 agrup = agr_in(i,1) + ETA 564 tauup = tau_in(i,1) + ETA 565 ndn = rho_in(i,2) + ETA 566 agrdn = agr_in(i,2) + ETA 567 taudn = tau_in(i,2) + ETA 568 569 n = nup + ndn 570 571* ***** M06-L Exchange ***** 572* ***** Up ***** 573 agrup2 = agrup*agrup 574 inv_nup = 1.0d0/nup 575 nup_13 = nup**thrd 576 nup_43 = nup_13*nup 577 nup_53 = nup_43*nup_13 578 nup_83 = nup_53*nup 579 580 xup = agrup2/nup_83 581 dxup_dnup = -etthrd*xup*inv_nup 582 dxup_dagrup = 2.0d0*agrup/nup_83 583 zup = tauup/nup_53 - cf 584 dzup_dnup = -fvthrd*tauup/nup_83 585 dzup_dtauup = 1.0d0/nup_53 586 587* ***** VS98 Exchange ***** 588 gamma = 1.0d0 + alphax*(xup + zup) 589 xr = xup/gamma 590 zr = zup/gamma 591 592 call nwpw_GVT4(aax,bbx,ccx,ddx,eex,ffx,alphax,alphax, 593 > xr,zr,gamma,GG,dGdx,dGdz) 594 595 eupx = -clda*nup_13*GG 596 fnupx = frthrd*eupx - clda*nup_43*(dGdx*dxup_dnup 597 & + dGdz*dzup_dnup) 598 fdnupx = -clda*nup_43*(dGdx*dxup_dagrup) 599 fdtauupx = -clda*nup_43*(dGdz*dzup_dtauup) 600 601* ***** M06-L Enhancement Factor ***** 602 tauU = P23*nup_53 603 t = tauU/tauup 604 dt_dn = fvthrd*t/nup 605 dt_dtau = -t/tauup 606 607 w = 1.0d0 608 dw_dt = 0.0d0 609 610 if(t .le. t1) then 611 w = (t - 1.0d0)/(t + 1.0d0) 612 dw_dt = 2.0d0/((1.0d0 + t)**2.0d0) 613 else if(t .gt. t1 .and. t .lt. t2) then 614 wt1 = (t1 - 1.0d0)/(t1 + 1.0d0) 615 dwt1 = 2.0d0/((1.0d0 + t1)**2.0d0) 616 aw = (3.0d0*t1*t2*t2 - t2**3.0d0)*wt1/(t1 - t2) 617 & + (t1**3.0d0 - 3.0d0*t1*t1*t2)/(t1 - t2) 618 & - t1*t2*t2*dwt1 619 aw = aw/((t1 - t2)*(t1 - t2)) 620 bw = -6.0d0*t1*t2*wt1/(t1 - t2) + 6.0d0*t1*t2/(t1 - t2) 621 & + (t2*t2 + 2.0d0*t1*t2)*dwt1 622 bw = bw/((t1 - t2)*(t1 - t2)) 623 cw = 3.0d0*(t1 + t2)*wt1/(t1 - t2) 624 & - 3.0d0*(t1 + t2)/(t1 - t2) - (t1 + 2.0d0*t2)*dwt1 625 cw = cw/((t1 - t2)*(t1 - t2)) 626 dw = -2.0d0*wt1/(t1 - t2) + 2.0d0/(t1 - t2) + dwt1 627 dw = dw/((t1 - t2)*(t1 - t2)) 628 629 w = aw + bw*t + cw*t*t + dw*t*t*t 630 dw_dt = bw + 2.0d0*cw*t + 3.0d0*dw*t*t 631 else if(t .ge. t2) then 632 w = 1.0d0 633 dw_dt = 0.0d0 634 end if 635 636 w1 = w 637 w2 = w1*w 638 w3 = w2*w 639 w4 = w3*w 640 w5 = w4*w 641 w6 = w5*w 642 w7 = w6*w 643 w8 = w7*w 644 w9 = w8*w 645 w10 = w9*w 646 w11 = w10*w 647 648 fw = ma0 + ma1*w1 + ma2*w2 + ma3*w3 + ma4*w4 649 & + ma5*w5 + ma6*w6 + ma7*w7 + ma8*w8 + ma9*w9 650 & + ma10*w10 + ma11*w11 651 dfw_dw = ma1 + 2.0d0*ma2*w1 + 3.0d0*ma3*w2 652 & + 4.0d0*ma4*w3 + 5.0d0*ma5*w4 + 6.0d0*ma6*w5 653 & + 7.0d0*ma7*w6 + 8.0d0*ma8*w7 + 9.0d0*ma9*w8 654 & + 10.0d0*ma10*w9 + 11.0d0*ma11*w10 655 dfw_dn = dfw_dw*dw_dt*dt_dn 656 dfw_dtau = dfw_dw*dw_dt*dt_dtau 657 658* ***** PBE96 Exchange ***** 659 n_onethird = (3.0d0*nup/pi)**thrd 660 ex_lda = -0.75d0*n_onethird 661 662 kf = (3.0d0*pi*pi*nup)**thrd 663 ss = agrup/(2.0d0*kf*nup) 664 P0 = 1.0d0 + (MU/KAPPA)*ss*ss 665 666 FF = (1.0d0 + KAPPA - KAPPA/P0) 667 Fs = 2.0d0*MU/(P0*P0)*ss 668 669 pbex = ex_lda*FF 670 pbefnx = frthrd*(pbex - ex_lda*Fs*ss) 671 pbefdnx = fdnx_const*Fs 672 673 eupx = eupx + pbex*fw 674 fnupx = fnupx + pbefnx*fw + nup*pbex*dfw_dn 675 fdnupx = fdnupx + pbefdnx*fw 676 fdtauupx = fdtauupx + nup*pbex*dfw_dtau 677 678 679* ***** Down ***** 680 agrdn2 = agrdn*agrdn 681 inv_ndn = 1.0d0/ndn 682 ndn_13 = ndn**thrd 683 ndn_43 = ndn_13*ndn 684 ndn_53 = ndn_43*ndn_13 685 ndn_83 = ndn_53*ndn 686 687 xdn = agrdn2/ndn_83 688 dxdn_dndn = -etthrd*xdn*inv_ndn 689 dxdn_dagrdn = 2.0d0*agrdn/ndn_83 690 zdn = taudn/ndn_53 - cf 691 dzdn_dndn = -fvthrd*taudn/ndn_83 692 dzdn_dtaudn = 1.0d0/ndn_53 693 694 gamma = 1.0d0 + alphax*(xdn + zdn) 695 xr = xdn/gamma 696 zr = zdn/gamma 697 698 call nwpw_GVT4(aax,bbx,ccx,ddx,eex,ffx,alphax,alphax, 699 > xr,zr,gamma,GG,dGdx,dGdz) 700 701 ednx = -clda*ndn_13*GG 702 fndnx = frthrd*ednx - clda*ndn_43*(dGdx*dxdn_dndn 703 & + dGdz*dzdn_dndn) 704 fdndnx = -clda*ndn_43*(dGdx*dxdn_dagrdn) 705 fdtaudnx = -clda*ndn_43*(dGdz*dzdn_dtaudn) 706 707* ***** M06-L Enhancement Factor ***** 708 tauU = P23*ndn_53 709 t = tauU/taudn 710 dt_dn = fvthrd*t/ndn 711 dt_dtau = -t/taudn 712 713 w = 1.0d0 714 dw_dt = 0.0d0 715 716 if(t .le. t1) then 717 w = (t - 1.0d0)/(t + 1.0d0) 718 dw_dt = 2.0d0/((1.0d0 + t)**2.0d0) 719 else if(t .gt. t1 .and. t .lt. t2) then 720 wt1 = (t1 - 1.0d0)/(t1 + 1.0d0) 721 dwt1 = 2.0d0/((1.0d0 + t1)**2.0d0) 722 aw = (3.0d0*t1*t2*t2 - t2**3.0d0)*wt1/(t1 - t2) 723 & + (t1**3.0d0 - 3.0d0*t1*t1*t2)/(t1 - t2) 724 & - t1*t2*t2*dwt1 725 aw = aw/((t1 - t2)*(t1 - t2)) 726 bw = -6.0d0*t1*t2*wt1/(t1 - t2) + 6.0d0*t1*t2/(t1 - t2) 727 & + (t2*t2 + 2.0d0*t1*t2)*dwt1 728 bw = bw/((t1 - t2)*(t1 - t2)) 729 cw = 3.0d0*(t1 + t2)*wt1/(t1 - t2) 730 & - 3.0d0*(t1 + t2)/(t1 - t2) - (t1 + 2.0d0*t2)*dwt1 731 cw = cw/((t1 - t2)*(t1 - t2)) 732 dw = -2.0d0*wt1/(t1 - t2) + 2.0d0/(t1 - t2) + dwt1 733 dw = dw/((t1 - t2)*(t1 - t2)) 734 735 w = aw + bw*t + cw*t*t + dw*t*t*t 736 dw_dt = bw + 2.0d0*cw*t + 3.0d0*dw*t*t 737 else if(t .ge. t2) then 738 w = 1.0d0 739 dw_dt = 0.0d0 740 end if 741 742 w1 = w 743 w2 = w1*w 744 w3 = w2*w 745 w4 = w3*w 746 w5 = w4*w 747 w6 = w5*w 748 w7 = w6*w 749 w8 = w7*w 750 w9 = w8*w 751 w10 = w9*w 752 w11 = w10*w 753 754 fw = ma0 + ma1*w1 + ma2*w2 + ma3*w3 + ma4*w4 755 & + ma5*w5 + ma6*w6 + ma7*w7 + ma8*w8 + ma9*w9 756 & + ma10*w10 + ma11*w11 757 dfw_dw = ma1 + 2.0d0*ma2*w1 + 3.0d0*ma3*w2 758 & + 4.0d0*ma4*w3 + 5.0d0*ma5*w4 + 6.0d0*ma6*w5 759 & + 7.0d0*ma7*w6 + 8.0d0*ma8*w7 + 9.0d0*ma9*w8 760 & + 10.0d0*ma10*w9 + 11.0d0*ma11*w10 761 dfw_dn = dfw_dw*dw_dt*dt_dn 762 dfw_dtau = dfw_dw*dw_dt*dt_dtau 763 764* ***** PBE96 Exchange ***** 765 n_onethird = (3.0d0*ndn/pi)**thrd 766 ex_lda = -0.75d0*n_onethird 767 768 kf = (3.0d0*pi*pi*ndn)**thrd 769 ss = agrdn/(2.0d0*kf*ndn) 770 P0 = 1.0d0 + (MU/KAPPA)*ss*ss 771 772 FF = (1.0d0 + KAPPA - KAPPA/P0) 773 Fs = 2.0d0*MU/(P0*P0)*ss 774 775 pbex = ex_lda*FF 776 pbefnx = frthrd*(pbex - ex_lda*Fs*ss) 777 pbefdnx = fdnx_const*Fs 778 779 ednx = ednx + pbex*fw 780 fndnx = fndnx + pbefnx*fw + ndn*pbex*dfw_dn 781 fdndnx = fdndnx + pbefdnx*fw 782 fdtaudnx = fdtaudnx + ndn*pbex*dfw_dtau 783 784 ex = (eupx*nup + ednx*ndn)/n 785 786* ***** M06-L Correlation ***** 787* ***** Same-Spin ***** 788* ***** alpha-alpha ***** 789 rs = (0.75d0/(pi*nup))**thrd 790 drs_dn = -thrd*rs/nup 791 792 call gen_PW91_c_rz(tol,rs,1.0d0,euu0c,deuu0c_drs,dummy) 793 deuu0c_dnup = deuu0c_drs*drs_dn 794 795 Dt = 1.0d0 - 0.25d0*xup/(zup + cf) 796 D = 0.0d0 797 dD_dn = 0.0d0 798 dD_dagr = 0.0d0 799 dD_dtau = 0.0d0 800 801 if (Dt .le. 0.0d0) then 802 D = 0.0d0 803 dD_dn = 0.0d0 804 dD_dagr = 0.0d0 805 dD_dtau = 0.0d0 806 else if (Dt .gt. 0.0d0 .and. Dt .lt. thresd) then 807 D = 2.0d0*Dt*Dt/thresd - Dt*Dt*Dt/(thresd*thresd) 808 dD_dDt = 4.0d0*Dt/thresd - 3.0d0*Dt*Dt/(thresd*thresd) 809 dDt_dx = -0.25d0/(zup + cf) 810 dDt_dz = 0.25d0*xup/((zup + cf)*(zup + cf)) 811 dD_dn = dD_dDt*(dDt_dx*dxup_dnup + dDt_dz*dzup_dnup) 812 dD_dagr = dD_dDt*dDt_dx*dxup_dagrup 813 dD_dtau = dD_dDt*dDt_dz*dzup_dtauup 814 else if(Dt .ge. thresd) then 815 D = Dt 816 dD_dx = -0.25d0/(zup + cf) 817 dD_dz = 0.25d0*xup/((zup + cf)*(zup + cf)) 818 dD_dn = dD_dx*dxup_dnup + dD_dz*dzup_dnup 819 dD_dagr = dD_dx*dxup_dagrup 820 dD_dtau = dD_dz*dzup_dtauup 821 end if 822 823 gamma = 1.0d0 + alphass*(xup + zup) 824 xr = xup/gamma 825 zr = zup/gamma 826 827 if (xup .ge. thresx) then 828 U = 1.0d0 829 dU_dx = 0.0d0 830 else 831 U = css*xup/(1.0d0 + css*xup) 832 dU_dx = css/((1.0d0 + css*xup)*(1.0d0 + css*xup)) 833 end if 834 835 U2 = U*U 836 U3 = U2*U 837 U4 = U3*U 838 839 gss = ccss0 + ccss1*U + ccss2*U2 + ccss3*U3 + ccss4*U4 840 dgss_dU = ccss1 + 2.0d0*ccss2*U + 3.0d0*ccss3*U2 841 & + 4.0d0*ccss4*U3 842 dgss_dx = dgss_dU*dU_dx 843 844 call nwpw_GVT4(aass,bbss,ccss,ddss,eess,ffss,alphass,alphass, 845 > xr,zr,gamma,hss,dhss_dx,dhss_dz) 846 847 ghss = gss + hss 848 dgh_dx = dgss_dx + dhss_dx 849 dgh_dn = dgh_dx*dxup_dnup + dhss_dz*dzup_dnup 850 dgh_dagr = dgh_dx*dxup_dagrup 851 dgh_dtau = dhss_dz*dzup_dtauup 852 853 eupc = euu0c*ghss*D 854 fnupc = eupc + nup*deuu0c_drs*drs_dn*ghss*D 855 & + nup*euu0c*dgh_dn*D 856 & + nup*euu0c*ghss*dD_dn 857 fdnupc = nup*euu0c*(dgh_dagr*D + ghss*dD_dagr) 858 fdtauupc = nup*euu0c*(dgh_dtau*D + ghss*dD_dtau) 859 860* ***** beta-beta ***** 861 rs = (0.75d0/(pi*ndn))**thrd 862 drs_dn = -thrd*rs/ndn 863 864 call gen_PW91_c_rz(tol,rs,1.0d0,edd0c,dedd0c_drs,dummy) 865 dedd0c_dndn = dedd0c_drs*drs_dn 866 867 Dt = 1.0d0 - 0.25d0*xdn/(zdn + cf) 868 D = 0.0d0 869 dD_dn = 0.0d0 870 dD_dagr = 0.0d0 871 dD_dtau = 0.0d0 872 873 if (Dt .le. 0.0d0) then 874 D = 0.0d0 875 dD_dn = 0.0d0 876 dD_dagr = 0.0d0 877 dD_dtau = 0.0d0 878 else if (Dt .gt. 0.0d0 .and. Dt .lt. thresd) then 879 D = 2.0d0*Dt*Dt/thresd - Dt*Dt*Dt/(thresd*thresd) 880 dD_dDt = 4.0d0*Dt/thresd - 3.0d0*Dt*Dt/(thresd*thresd) 881 dDt_dx = -0.25d0/(zdn + cf) 882 dDt_dz = 0.25d0*xdn/((zdn + cf)*(zdn + cf)) 883 dD_dn = dD_dDt*(dDt_dx*dxdn_dndn + dDt_dz*dzdn_dndn) 884 dD_dagr = dD_dDt*dDt_dx*dxdn_dagrdn 885 dD_dtau = dD_dDt*dDt_dz*dzdn_dtaudn 886 else if(Dt .ge. thresd) then 887 D = Dt 888 dD_dx = -0.25d0/(zdn + cf) 889 dD_dz = 0.25d0*xdn/((zdn + cf)*(zdn + cf)) 890 dD_dn = dD_dx*dxdn_dndn + dD_dz*dzdn_dndn 891 dD_dagr = dD_dx*dxdn_dagrdn 892 dD_dtau = dD_dz*dzdn_dtaudn 893 end if 894 895 gamma = 1.0d0 + alphass*(xdn + zdn) 896 xr = xdn/gamma 897 zr = zdn/gamma 898 899 if (xdn .ge. thresx) then 900 U = 1.0d0 901 dU_dx = 0.0d0 902 else 903 U = css*xdn/(1.0d0 + css*xdn) 904 dU_dx = css/((1.0d0 + css*xdn)*(1.0d0 + css*xdn)) 905 end if 906 907 U2 = U*U 908 U3 = U2*U 909 U4 = U3*U 910 911 gss = ccss0 + ccss1*U + ccss2*U2 + ccss3*U3 + ccss4*U4 912 dgss_dU = ccss1 + 2.0d0*ccss2*U + 3.0d0*ccss3*U2 913 & + 4.0d0*ccss4*U3 914 dgss_dx = dgss_dU*dU_dx 915 916 call nwpw_GVT4(aass,bbss,ccss,ddss,eess,ffss,alphass,alphass, 917 > xr,zr,gamma,hss,dhss_dx,dhss_dz) 918 919 ghss = gss + hss 920 dgh_dx = dgss_dx + dhss_dx 921 dgh_dn = dgh_dx*dxdn_dndn + dhss_dz*dzdn_dndn 922 dgh_dagr = dgh_dx*dxdn_dagrdn 923 dgh_dtau = dhss_dz*dzdn_dtaudn 924 925 ednc = edd0c*ghss*D 926 fndnc = ednc + ndn*dedd0c_drs*drs_dn*ghss*D 927 & + ndn*edd0c*dgh_dn*D 928 & + ndn*edd0c*ghss*dD_dn 929 fdndnc = ndn*edd0c*(dgh_dagr*D + ghss*dD_dagr) 930 fdtaudnc = ndn*edd0c*(dgh_dtau*D + ghss*dD_dtau) 931 932 ec = (eupc*nup + ednc*ndn)/n 933 934* ***** Opposite-Spin (alpha-beta) ***** 935 rs = (0.75d0/(pi*n))**thrd 936 drs_dn = -thrd*rs/n 937 zeta = (nup - ndn)/n 938 dzeta_dnup = ( 1.0d0 - zeta)/n 939 dzeta_dndn = (-1.0d0 - zeta)/n 940 941 call gen_PW91_c_rz(tol,rs,zeta,eud0c,deud0c_drs,deud0c_dzeta) 942 eud1c = eud0c - (euu0c*nup + edd0c*ndn)/n 943 fud1c = n*eud1c 944 dfud1c_dnup = eud0c + n*(deud0c_drs*drs_dn 945 & + deud0c_dzeta*dzeta_dnup) - euu0c 946 & - nup*deuu0c_dnup 947 dfud1c_dndn = eud0c + n*(deud0c_drs*drs_dn 948 & + deud0c_dzeta*dzeta_dndn) - edd0c 949 & - ndn*dedd0c_dndn 950 951 x = xup + xdn 952 z = zup + zdn 953 954 gamma = 1.0d0 + alphaopp*(x + z) 955 xr = x/gamma 956 zr = z/gamma 957 958 if (x .ge. thresx) then 959 U = 1.0d0 960 dU_dx = 0.0d0 961 else 962 U = copp*x/(1.0d0 + copp*x) 963 dU_dx = copp/((1.0d0 + copp*x)*(1.0d0 + copp*x)) 964 end if 965 966 U2 = U*U 967 U3 = U2*U 968 U4 = U3*U 969 970 gopp = ccopp0 + ccopp1*U + ccopp2*U2 + ccopp3*U3 + ccopp4*U4 971 dgopp_dU = ccopp1 + 2.0d0*ccopp2*U 972 & + 3.0d0*ccopp3*U2 + 4.0d0*ccopp4*U3 973 dgopp_dx = dgopp_dU*dU_dx 974 975 call nwpw_GVT4(aaopp,bbopp,ccopp,ddopp,eeopp,ffopp, 976 > alphaopp,alphaopp, 977 > xr,zr,gamma,hopp,dhopp_dx,dhopp_dz) 978 979 ghopp = gopp + hopp 980 dgh_dx = dgopp_dx + dhopp_dx 981 dgh_dnup = dgh_dx*dxup_dnup + dhopp_dz*dzup_dnup 982 dgh_dndn = dgh_dx*dxdn_dndn + dhopp_dz*dzdn_dndn 983 dgh_dagrup = dgh_dx*dxup_dagrup 984 dgh_dagrdn = dgh_dx*dxdn_dagrdn 985 dgh_dtauup = dhopp_dz*dzup_dtauup 986 dgh_dtaudn = dhopp_dz*dzdn_dtaudn 987 988 eudc = eud1c*ghopp 989 dfudc_dnup = dfud1c_dnup*ghopp + fud1c*dgh_dnup 990 dfudc_dndn = dfud1c_dndn*ghopp + fud1c*dgh_dndn 991 dfudc_dagrup = fud1c*dgh_dagrup 992 dfudc_dagrdn = fud1c*dgh_dagrdn 993 dfudc_dtauup = fud1c*dgh_dtauup 994 dfudc_dtaudn = fud1c*dgh_dtaudn 995 996 ec = ec + eudc 997 fnupc = fnupc + dfudc_dnup 998 fndnc = fndnc + dfudc_dndn 999 fdnupc = fdnupc + dfudc_dagrup 1000 fdndnc = fdndnc + dfudc_dagrdn 1001 fdtauupc = fdtauupc + dfudc_dtauup 1002 fdtaudnc = fdtaudnc + dfudc_dtaudn 1003 1004 xce(i) = x_parameter*ex + c_parameter*ec 1005 fn(i,1) = x_parameter*fnupx + c_parameter*fnupc 1006 fn(i,2) = x_parameter*fndnx + c_parameter*fndnc 1007 fdn(i,1) = x_parameter*fdnupx + c_parameter*fdnupc 1008 fdn(i,2) = x_parameter*fdndnx + c_parameter*fdndnc 1009 fdn(i,3) = 0.0d0 1010 fdtau(i,1) = x_parameter*fdtauupx + c_parameter*fdtauupc 1011 fdtau(i,2) = x_parameter*fdtaudnx + c_parameter*fdtaudnc 1012 end do 1013 1014 return 1015 end 1016 1017