1* ************************************************ 2* * * 3* * nwpw_m06_x * 4* * * 5* ************************************************ 6 subroutine nwpw_m06_x(pi,thrd,fvthrd,etthrd, 7 > a0,a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11, 8 > C1,C2,Cx,P23, 9 > n,agr,tau, 10 > xe,dfdnx,dfdagrx,dfdtaux) 11 implicit none 12* ***** input ***** 13 real*8 pi,thrd,fvthrd,etthrd 14 real*8 a0,a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11 15 real*8 C1,C2,Cx,P23 16 real*8 n,agr,tau 17* ***** output ***** 18 real*8 xe,dfdnx,dfdagrx,dfdtaux 19* ***** local declarations ***** 20 real*8 n_13,n_23,n_53,n_83,agr2 21 real*8 tauU,t,dt_dn,dt_dtau 22 real*8 w,dw_dt,w1,w2,w3,w4,w5,w6,w7,w8,w9,w10,w11 23 real*8 fw,dfw_dw,dfw_dn,dfw_dtau 24 real*8 x,dx_dn,dx_dagr 25 real*8 enum,eden,eden2,etmp,detmp_dx 26 real*8 Fpbe,dFpbe_dn,dFpbe_dagr 27 28 n_13 = n**thrd 29 n_23 = n_13*n_13 30 n_53 = n_23*n 31 n_83 = n_53*n 32 agr2 = agr*agr 33 34 tauU = P23*n_53 35 t = tauU/tau 36 dt_dn = fvthrd*t/n 37 dt_dtau = -t/tau 38 39 w = (t - 1.0d0)/(t + 1.0d0) 40 dw_dt = 2.0d0/((1.0d0 + t)**2.0d0) 41 42 w1 = w 43 w2 = w1*w 44 w3 = w2*w 45 w4 = w3*w 46 w5 = w4*w 47 w6 = w5*w 48 w7 = w6*w 49 w8 = w7*w 50 w9 = w8*w 51 w10 = w9*w 52 w11 = w10*w 53 54 fw = a0 + a1*w1 + a2*w2 + a3*w3 + a4*w4 + a5*w5 55 & + a6*w6 + a7*w7 + a8*w8 + a9*w9 + a10*w10 + a11*w11 56 dfw_dw = a1 + 2.0d0*a2*w1 + 3.0d0*a3*w2 57 & + 4.0d0*a4*w3 + 5.0d0*a5*w4 + 6.0d0*a6*w5 58 & + 7.0d0*a7*w6 + 8.0d0*a8*w7 + 9.0d0*a9*w8 59 & + 10.0d0*a10*w9 + 11.0d0*a11*w10 60 dfw_dn = dfw_dw*dw_dt*dt_dn 61 dfw_dtau = dfw_dw*dw_dt*dt_dtau 62 63 x = agr2/n_83 64 dx_dn = -etthrd*x/n 65 dx_dagr = 2.0d0*x/agr 66 67 enum = C1*x 68 eden = 1.0d0 + C2*x 69 eden2 = eden*eden 70 etmp = -enum/eden 71 detmp_dx = -(C1*eden - C2*enum)/eden2 72 73 Fpbe = n_13*(Cx + etmp) 74 dFpbe_dn = n_13*detmp_dx*dx_dn + thrd*Fpbe/n 75 dFpbe_dagr = n_13*detmp_dx*dx_dagr 76 77 xe = Fpbe*fw 78 dfdnx = n*(dFpbe_dn*fw + Fpbe*dfw_dn) + xe 79 dfdagrx = n*dFpbe_dagr*fw 80 dfdtaux = n*Fpbe*dfw_dtau 81 82 return 83 end 84* ************************************************ 85* * * 86* * nwpw_m06_c_ss * 87* * * 88* ************************************************ 89 subroutine nwpw_m06_c_ss(cgss,css0,css1,css2,css3,css4, 90 > n,xs,dxs_dn,dxs_dagr,zs,dzs_dn,dzs_dtau, 91 > ess0c,dess0c_dn, 92 > cess,dfss_dn,dfss_dagr,dfss_dtau) 93 implicit none 94* ***** input ***** 95 real*8 cgss,css0,css1,css2,css3,css4 96 real*8 n,xs,dxs_dn,dxs_dagr,zs,dzs_dn,dzs_dtau 97* ***** output ***** 98 real*8 ess0c,dess0c_dn 99 real*8 cess,dfss_dn,dfss_dagr,dfss_dtau 100* ***** local declarations ***** 101 real*8 D,dD_dx,dD_dz,dD_dn,dD_dagr,dD_dtau 102 real*8 U,dU_dx,U2,U3,U4,gss,dgss_dU 103 real*8 rs,drs_dn 104 real*8 pwc,dpwc_drs,dummy,ness0c 105* ***** constants ***** 106 real*8 pi,thrd 107 parameter (pi = 3.14159265358979311599d0) 108 parameter (thrd = 1.0d0/3.0d0) 109* ***** density cutoff parameters ***** 110 real*8 tol 111 parameter (tol = 1.0d-18) 112* ***** M06 constants ***** 113 real*8 cf 114 parameter (cf = 9.115599720d0) 115 116 rs = (0.75d0/(pi*n))**thrd 117 drs_dn = -thrd*rs/n 118 119 call gen_PW91_c_rz(tol,rs,1.0d0,pwc,dpwc_drs,dummy) 120 121 ess0c = pwc 122 ness0c = n*ess0c 123 dess0c_dn = dpwc_drs*drs_dn 124 125 U = cgss*xs/(1.0d0 + cgss*xs) 126 dU_dx = cgss/((1.0d0 + cgss*xs)**2.0d0) 127 128 U2 = U*U 129 U3 = U2*U 130 U4 = U3*U 131 132 gss = css0 + css1*U + css2*U2 + css3*U3 + css4*U4 133 dgss_dU = css1 + 2.0d0*css2*U + 3.0d0*css3*U2 + 4.0d0*css4*U3 134 135 D = 1.0d0 - 0.25d0*xs/(zs+cf) 136 if (D .lt. 0.0d0) then 137 D = 0.0d0 138 dD_dn = 0.0d0 139 dD_dagr = 0.0d0 140 dD_dtau = 0.0d0 141 else 142 dD_dx = -1.0d0/(4.0d0*(zs + cf)) 143 dD_dz = xs/(4.0d0*(zs + cf)*(zs + cf)) 144 dD_dn = dD_dx*dxs_dn + dD_dz*dzs_dn 145 dD_dagr = dD_dx*dxs_dagr 146 dD_dtau = dD_dz*dzs_dtau 147 end if 148 149 cess = ess0c*gss*D 150 dfss_dn = cess + n*dess0c_dn*gss*D 151 & + ness0c*dgss_dU*dU_dx*dxs_dn*D + ness0c*gss*dD_dn 152 dfss_dagr = ness0c*(dgss_dU*dU_dx*dxs_dagr*D + gss*dD_dagr) 153 dfss_dtau = ness0c*gss*dD_dtau 154 155 return 156 end 157 158* ************************************************ 159* * * 160* * nwpw_m06_c_restricted * 161* * * 162* ************************************************ 163 subroutine nwpw_m06_c_restricted(cgss,css0,css1,css2,css3,css4, 164 > cgopp,copp0,copp1,copp2,copp3,copp4, 165 > n,nup,xup,dxup_dnup,dxup_dagrup, 166 > zup,dzup_dnup,dzup_dtauup,rs,drs_dn, 167 > ce,fnc,fdnc,fdtauc) 168 implicit none 169* ***** input ***** 170 real*8 cgss,css0,css1,css2,css3,css4 171 real*8 cgopp,copp0,copp1,copp2,copp3,copp4 172 real*8 n,nup,xup,dxup_dnup,dxup_dagrup 173 real*8 zup,dzup_dnup,dzup_dtauup 174 real*8 rs,drs_dn 175* ***** output ***** 176 real*8 ce,fnc,fdnc,fdtauc 177* ***** local declarations ***** 178 real*8 euu0c,deuu0c_dn 179 real*8 euuc,dfuuc_dn,dfuuc_dagr,dfuuc_dtau 180 real*8 eud0c,deud0c_drs,dummy,eud1c,fud1c,dfud1c_dn 181 real*8 x,U,dU_dx,U2,U3,U4,gopp,dgopp_dU 182 real*8 eudc,dfudc_dn,dfudc_dagr,dfudc_dtau 183* ***** density cutoff parameters ***** 184 real*8 tol 185 parameter (tol = 1.0d-18) 186 187* ***** same spin (alpha-alpha or beta-beta) ***** 188 call nwpw_m06_c_ss(cgss,css0,css1,css2,css3,css4, 189 > nup,xup,dxup_dnup,dxup_dagrup, 190 > zup,dzup_dnup,dzup_dtauup, 191 > euu0c,deuu0c_dn, 192 > euuc,dfuuc_dn,dfuuc_dagr,dfuuc_dtau) 193 194* ***** opposite spin (alpha-beta) ***** 195 call gen_PW91_c_rz(tol,rs,0.0d0,eud0c,deud0c_drs,dummy) 196 eud1c = eud0c - euu0c 197 fud1c = n*eud1c 198 dfud1c_dn = eud0c + n*deud0c_drs*drs_dn 199 & - euu0c - nup*deuu0c_dn 200 201 x = 2.0d0*xup 202 U = cgopp*x/(1.0d0 + cgopp*x) 203 dU_dx = cgopp/((1.0d0 + cgopp*x)**2.0d0) 204 205 U2 = U*U 206 U3 = U2*U 207 U4 = U3*U 208 209 gopp = copp0 + copp1*U + copp2*U2 + copp3*U3 + copp4*U4 210 dgopp_dU = copp1 + 2.0d0*copp2*U 211 & + 3.0d0*copp3*U2 + 4.0d0*copp4*U3 212 213 eudc = eud1c*gopp 214 dfudc_dn = dfud1c_dn*gopp + fud1c*dgopp_dU*dU_dx*dxup_dnup 215 dfudc_dagr = fud1c*dgopp_dU*dU_dx*dxup_dagrup 216 dfudc_dtau = 0.0d0 217 218 ce = eudc + euuc 219 fnc = dfudc_dn + dfuuc_dn 220 fdnc = dfudc_dagr + dfuuc_dagr 221 fdtauc = dfuuc_dtau 222 223 return 224 end 225 226* ************************************************ 227* * * 228* * nwpw_m06_c_unrestricted * 229* * * 230* ************************************************ 231 subroutine nwpw_m06_c_unrestricted(cgss,css0,css1,css2,css3,css4, 232 > cgopp,copp0,copp1,copp2,copp3,copp4, 233 > n,nup,agrup,tauup,ndn,agrdn,taudn, 234 > xup,dxup_dnup,dxup_dagrup,zup,dzup_dnup,dzup_dtauup, 235 > xdn,dxdn_dndn,dxdn_dagrdn,zdn,dzdn_dndn,dzdn_dtaudn, 236 > rs,drs_dn,zeta,dzeta_dnup,dzeta_dndn, 237 > ce,fnupc,fndnc,fdnupc,fdndnc,fdtauupc,fdtaudnc) 238 implicit none 239* ***** input ***** 240 real*8 cgss,css0,css1,css2,css3,css4 241 real*8 cgopp,copp0,copp1,copp2,copp3,copp4 242 real*8 n,nup,agrup,tauup,ndn,agrdn,taudn 243 real*8 xup,dxup_dnup,dxup_dagrup,zup,dzup_dnup,dzup_dtauup 244 real*8 xdn,dxdn_dndn,dxdn_dagrdn,zdn,dzdn_dndn,dzdn_dtaudn 245 real*8 rs,drs_dn,zeta,dzeta_dnup,dzeta_dndn 246* ***** output ***** 247 real*8 ce,fnupc,fndnc,fdnupc,fdndnc,fdtauupc,fdtaudnc 248* ***** local declarations ***** 249 real*8 euu0c,deuu0c_dnup 250 real*8 euuc,dfuuc_dnup,dfuuc_dagrup,dfuuc_dtauup 251 real*8 edd0c,dedd0c_dndn 252 real*8 eddc,dfddc_dndn,dfddc_dagrdn,dfddc_dtaudn 253 real*8 eud0c,deud0c_drs,deud0c_dzeta 254 real*8 x,U,dU_dx,U2,U3,U4,gopp,dgopp_dU 255 real*8 eudc,eud1c,fud1c,dfud1c_dnup,dfud1c_dndn 256 real*8 dfudc_dnup,dfudc_dagrup,dfudc_dtauup 257 real*8 dfudc_dndn,dfudc_dagrdn,dfudc_dtaudn 258* ***** density cutoff parameters ***** 259 real*8 tol 260 parameter (tol = 1.0d-18) 261 262* ***** same spin ***** 263* ***** alpha-alpha ***** 264 call nwpw_m06_c_ss(cgss,css0,css1,css2,css3,css4, 265 > nup,xup,dxup_dnup,dxup_dagrup, 266 > zup,dzup_dnup,dzup_dtauup, 267 > euu0c,deuu0c_dnup, 268 > euuc,dfuuc_dnup,dfuuc_dagrup,dfuuc_dtauup) 269* ***** beta-beta ***** 270 call nwpw_m06_c_ss(cgss,css0,css1,css2,css3,css4, 271 > ndn,xdn,dxdn_dndn,dxdn_dagrdn, 272 > zdn,dzdn_dndn,dzdn_dtaudn, 273 > edd0c,dedd0c_dndn, 274 > eddc,dfddc_dndn,dfddc_dagrdn,dfddc_dtaudn) 275 276* ***** opposite spin (alpha-beta) ***** 277 call gen_PW91_c_rz(tol,rs,zeta,eud0c,deud0c_drs,deud0c_dzeta) 278 eud1c = eud0c - (euu0c*nup + edd0c*ndn)/n 279 fud1c = n*eud1c 280 281 x = xup + xdn 282 U = cgopp*x/(1.0d0 + cgopp*x) 283 dU_dx = cgopp/((1.0d0 + cgopp*x)**2.0d0) 284 285 U2 = U*U 286 U3 = U2*U 287 U4 = U3*U 288 289 gopp = copp0 + copp1*U + copp2*U2 + copp3*U3 + copp4*U4 290 dgopp_dU = copp1 + 2.0d0*copp2*U 291 & + 3.0d0*copp3*U2 + 4.0d0*copp4*U3 292 293 eudc = eud1c*gopp 294 dfud1c_dnup = eud0c + n*(deud0c_drs*drs_dn 295 & + deud0c_dzeta*dzeta_dnup) - euu0c 296 & - nup*deuu0c_dnup 297 dfud1c_dndn = eud0c + n*(deud0c_drs*drs_dn 298 & + deud0c_dzeta*dzeta_dndn) - edd0c 299 & - ndn*dedd0c_dndn 300 301 dfudc_dnup = dfud1c_dnup*gopp + fud1c*dgopp_dU*dU_dx*dxup_dnup 302 dfudc_dndn = dfud1c_dndn*gopp + fud1c*dgopp_dU*dU_dx*dxdn_dndn 303 dfudc_dagrup = fud1c*dgopp_dU*dU_dx*dxup_dagrup 304 dfudc_dagrdn = fud1c*dgopp_dU*dU_dx*dxdn_dagrdn 305 dfudc_dtauup = 0.0d0 306 dfudc_dtaudn = 0.0d0 307 308 ce = eudc + (euuc*nup + eddc*ndn)/n 309 fnupc = dfudc_dnup + dfuuc_dnup 310 fndnc = dfudc_dndn + dfddc_dndn 311 fdnupc = dfudc_dagrup + dfuuc_dagrup 312 fdndnc = dfudc_dagrdn + dfddc_dagrdn 313 fdtauupc = dfuuc_dtauup 314 fdtaudnc = dfddc_dtaudn 315 316 return 317 end 318