1 subroutine x_rks_s(r,f,dfdra) 2 implicit none 3c 4c This subroutine evaluates the spin polarised exchange functional 5c in the Local Density Approximation [1], and the corresponding 6c potential. Often this functional is referred to as the Dirac 7c functional [2] or Slater functional. 8c 9c [1] F. Bloch, Zeitschrift fuer Physik, Vol. 57 (1929) 545. 10c 11c [2] P.A.M. Dirac, Proceedings of the Cambridge Philosophical 12c Society, Vol. 26 (1930) 376. 13c 14c Parameters: 15c 16c r the total electron density 17c f On return the functional value 18c dfdra On return the derivative of f with respect to alpha electron 19c density 20c 21c 22c Ax = -3/4*(6/pi)**(1/3) 23c Bx = -(6/pi)**(1/3) 24c C = (1/2)**(1/3) 25 real*8 Ax, Bx, C 26 parameter(Ax = -0.930525736349100025002010218071667d0) 27 parameter(Bx = -1.24070098179880003333601362409556d0) 28 parameter(C = 0.793700525984099737375852819636154d0) 29c 30 real*8 third 31 parameter(third = 0.333333333333333333333333333333333d0) 32c 33 real*8 r 34 real*8 f, dfdra 35c 36 real*8 ra13 37c 38 ra13 = C*r**third 39 f = Ax*r*ra13 40 dfdra = Bx*ra13 41c 42 end 43c----------------------------------------------------------------------- 44 subroutine x_uks_s(ra,rb,f,dfdra,dfdrb) 45 implicit none 46c 47c This subroutine evaluates the spin polarised exchange functional 48c in the Local Density Approximation [1], and the corresponding 49c potential. Often this functional is referred to as the Dirac 50c functional [2] or Slater functional. 51c 52c [1] F. Bloch, Zeitschrift fuer Physik, Vol. 57 (1929) 545. 53c 54c [2] P.A.M. Dirac, Proceedings of the Cambridge Philosophical 55c Society, Vol. 26 (1930) 376. 56c 57c Parameters: 58c 59c ra the alpha electron density 60c rb the beta electron density 61c f On return the functional value 62c dfdra On return the derivative of f with respect to ra 63c dfdrb On return the derivative of f with respect to rb 64c 65c 66c Ax = -3/4*(6/pi)**(1/3) 67c Bx = -(6/pi)**(1/3) 68 real*8 Ax, Bx 69 parameter(Ax = -0.930525736349100025002010218071667d0) 70 parameter(Bx = -1.24070098179880003333601362409556d0) 71c 72 real*8 third 73 parameter(third = 0.333333333333333333333333333333333d0) 74c 75 real*8 ra, rb 76 real*8 f, dfdra, dfdrb 77c 78 real*8 ra13, rb13 79c 80 ra13 = ra**third 81 rb13 = rb**third 82 f = Ax*(ra*ra13+rb*rb13) 83 dfdra = Bx*ra13 84 dfdrb = Bx*rb13 85c 86 end 87 subroutine c_rks_vwn5(r,f,dfdra) 88 implicit none 89c 90c This subroutine evaluates the Vosko, Wilk and Nusair correlation 91c functional number 5 [1] for the closed shell case, with the 92c parametrisation as given in table 5. 93c 94c The original code was obtained from Dr. Phillip Young, 95c with corrections from Dr. Paul Sherwood. 96c 97c [1] S.H. Vosko, L. Wilk, and M. Nusair 98c "Accurate spin-dependent electron liquid correlation energies 99c for local spin density calculations: a critical analysis", 100c Can.J.Phys, Vol. 58 (1980) 1200-1211. 101c 102c Parameters: 103c 104c r the total electron density 105c f On return the functional value 106c dfdra On return the derivative of f with respect to the alpha 107c electron density 108c 109 real*8 r 110 real*8 f, dfdra 111c 112 real*8 t4,t5,t6,t7 113 real*8 a2,b2,c2,d2 114 real*8 P1,P2,P3,P4 115 real*8 alpha_rho13 116 real*8 srho,srho13 117 real*8 iv,iv2,inv,i1,i2,i3 118 real*8 pp1,pp2 119c 120 real*8 n13, n43, n16, one, two, four 121 parameter(n13 = 0.333333333333333333333333333333d0) 122 parameter(n43 = 1.333333333333333333333333333333d0) 123 parameter(n16 = 0.166666666666666666666666666666d0) 124 parameter(one = 1.0d0) 125 parameter(two = 2.0d0) 126 parameter(four= 4.0d0) 127C 128C VWN interpolation parameters 129C 130C paramagnetic 131 a2 = 0.0621814d0 132 b2 = 3.72744d0 133 c2 = 12.9352d0 134 d2 = -0.10498d0 135C 136C t4 = (1/(4/3)*pi)**(1/3) 137 t4=0.620350490899399531d0 138C 139C t5 = 0.5/(2**(1/3)-1) 140 t5 = 1.92366105093153617d0 141C 142C t6 = 2/(3*(2**(1/3)-1)) 143 t6 = 2.56488140124204822d0 144C 145C t7 = 2.25*(2**(1/3)-1) 146 t7 = 0.584822362263464735d0 147C 148C Paramagnetic interpolation constants 149C 150 P1 = 6.15199081975907980d0 151 P2 = a2 *0.5d0 152 P3 = 0.193804554230887489d-02 *0.5d0 153 P4 = 0.775665897562260176d-01 *0.5d0 154C 155C closed shell case 156 srho = r 157 srho13 = srho**n13 158 alpha_rho13 = (0.5d0**n13)*srho 159 iv2 = T4/srho13 160 iv = sqrt(iv2) 161C 162C paramagnetic 163 inv = 1.0d0/(iv2+b2*iv+c2) 164 i1 = log(iv2*inv) 165 i2 = log((iv-d2)*(iv-d2)*inv) 166c corrected b1->b2 ps Apr98 167 i3 = atan(P1/(2.0d0*iv+b2)) 168 pp1 = P2*i1 + P3*i2 + P4*i3 169 pp2 = a2*(1.0d0/iv-iv*inv*(1.0d0+b2/(iv-d2))) 170C 171 f = pp1*srho 172 dfdra = pp1 - n16*iv*pp2 173c 174 end 175c----------------------------------------------------------------------- 176 subroutine c_uks_vwn5(ra,rb,f,dfdra,dfdrb) 177 implicit none 178c 179c This subroutine evaluates the Vosko, Wilk and Nusair correlation 180c functional number 5 [1], with the parametrisation as given in 181c table 5. 182c 183c The original code was obtained from Dr. Phillip Young, 184c with corrections from Dr. Paul Sherwood. 185c 186c [1] S.H. Vosko, L. Wilk, and M. Nusair 187c "Accurate spin-dependent electron liquid correlation energies 188c for local spin density calculations: a critical analysis", 189c Can.J.Phys, Vol. 58 (1980) 1200-1211. 190c 191c Parameters: 192c 193c ra the alpha-electron density 194c rb the beta-electron density 195c f On return the functional value 196c dfdra On return the derivative of f with respect to ra 197c dfdrb On return the derivative of f with respect to rb 198c 199 real*8 ra, rb 200 real*8 f, dfdra, dfdrb 201c 202 real*8 t1,t2,t4,t5,t6,t7 203 real*8 a1,b1,c1,d1 204 real*8 a2,b2,c2,d2 205 real*8 a3,b3,c3,d3 206 real*8 S1,S2,S3,S4 207 real*8 P1,P2,P3,P4 208 real*8 F1,F2,F3,F4 209 real*8 inter1,inter2 210 real*8 alpha_rho13,beta_rho13 211 real*8 srho,srho13 212 real*8 iv,iv2,inv,i1,i2,i3 213 real*8 vwn1,vwn2 214 real*8 ss1,ss2,pp1,pp2,ff1,ff2 215 real*8 zeta,zeta3,zeta4 216 real*8 tau,dtau,v 217c 218 real*8 n13, n43, n16, one, two, four, tn13, tn43 219c tn13 = 2**(1/3) 220c tn43 = 2**(4/3) 221 parameter(n13 = 0.333333333333333333333333333333d0) 222 parameter(n43 = 1.333333333333333333333333333333d0) 223 parameter(n16 = 0.166666666666666666666666666667d0) 224 parameter(tn13= 1.25992104989487316476721060727823d0) 225 parameter(tn43= 2.51984209978974632953442121455646d0) 226 parameter(one = 1.0d0) 227 parameter(two = 2.0d0) 228 parameter(four= 4.0d0) 229C 230C VWN interpolation parameters 231C 232C spin stiffness 233 a1 = -0.337737278807791058d-01 234 b1 = 1.13107d0 235 c1 = 13.0045d0 236 d1 = -0.00475840d0 237C paramagnetic 238 a2 = 0.0621814d0 239 b2 = 3.72744d0 240 c2 = 12.9352d0 241 d2 = -0.10498d0 242C ferromagnetic 243c try cadpac/nwchem value (.5*a2) 244 a3 = 0.0310907d0 245 b3 = 7.06042d0 246 c3 = 18.0578d0 247 d3 = -0.32500d0 248C 249C t4 = (1/(4/3)*pi)**(1/3) 250 t4=0.620350490899399531d0 251C 252C t5 = 0.5/(2**(1/3)-1) 253 t5 = 1.92366105093153617d0 254C 255C t6 = 2/(3*(2**(1/3)-1)) 256 t6 = 2.56488140124204822d0 257C 258C t7 = 2.25*(2**(1/3)-1) 259 t7 = 0.584822362263464735d0 260C 261C Spin stiffness interpolation constants 262C 263 S1 = 7.12310891781811772d0 264 S2 = a1 *0.5d0 265 S3 = -0.139834647015288626d-04 *0.5d0 266 S4 = -0.107301836977671539d-01 *0.5d0 267C 268C Paramagnetic interpolation constants 269C 270 P1 = 6.15199081975907980d0 271 P2 = a2 *0.5d0 272 P3 = 0.193804554230887489d-02 *0.5d0 273 P4 = 0.775665897562260176d-01 *0.5d0 274C 275C Ferromagnetic interpolation constants 276C 277 F1 = 4.73092690956011364d0 278 F2 = a3 *0.5d0 279c 280c F3 = -0.244185082989490298d-02 *0.5d0 281c F4 = -0.570212323620622186d-01 *0.5d0 282c 283c try nwchem values 284c 285 F3 = 0.224786709554261133d-02 286 F4 = 0.524913931697809227d-01 287C 288C Interpolation intervals 289C 290 inter1 = 1.0d0-1.0d-10 291 inter2 = -1.0d0+1.0d-10 292C 293C open shell case 294 alpha_rho13 = ra**n13 295 beta_rho13 = rb**n13 296 srho = ra+rb 297 srho13 = srho**n13 298 iv2 = T4/srho13 299 iv = sqrt(iv2) 300C 301C spin-stiffness 302 inv = 1.0d0/(iv2+b1*iv+c1) 303 i1 = log(iv2*inv) 304 i2 = log((iv-d1)*(iv-d1)*inv) 305 i3 = atan(S1/(2.0d0*iv+b1)) 306 ss1 = S2*i1 + S3*i2 + S4*i3 307 ss2 = a1*(1.0d0/iv-iv*inv*(1.0d0+b1/(iv-d1))) 308C 309C paramagnetic 310 inv = 1.0d0/(iv2+b2*iv+c2) 311 i1 = log(iv2*inv) 312 i2 = log((iv-d2)*(iv-d2)*inv) 313c corrected b1->b2 ps Apr98 314 i3 = atan(P1/(2.0d0*iv+b2)) 315 pp1 = P2*i1 + P3*i2 + P4*i3 316 pp2 = a2*(1.0d0/iv-iv*inv*(1.0d0+b2/(iv-d2))) 317C 318C ferromagnetic 319 inv = 1.0d0/(iv2+b3*iv+c3) 320 i1 = log(iv2*inv) 321 i2 = log((iv-d3)*(iv-d3)*inv) 322 i3 = atan(F1/(2.0d0*iv+b3)) 323 ff1 = F2*i1 + F3*i2 + F4*i3 324 ff2 = a3*(1.0d0/iv-iv*inv*(1.0d0+b3/(iv-d3))) 325C 326C polarisation function 327c 328 zeta = (ra-rb)/srho 329 zeta3 = zeta*zeta*zeta 330 zeta4 = zeta3*zeta 331 if(zeta.gt.inter1) then 332 vwn1 = (tn43-two)*t5 333 vwn2 = (tn13)*t6 334 elseif(zeta.lt.inter2) then 335 vwn1 = (tn43-two)*t5 336 vwn2 = -(tn13)*t6 337 else 338 vwn1 = ((one+zeta)**n43+(one-zeta)**n43-two)*t5 339 vwn2 = ((one+zeta)**n13-(one-zeta)**n13)*t6 340 endif 341 ss1 = ss1*t7 342 ss2 = ss2*t7 343 tau = ff1-pp1-ss1 344 dtau = ff2-pp2-ss2 345c 346 v = pp1+vwn1*(ss1+tau*zeta4) 347 f = v*srho 348c 349 t1 = v - n16*iv*(pp2+vwn1*(ss2+dtau*zeta4)) 350 t2 = vwn2*(ss1+tau*zeta4)+vwn1*four*tau*zeta3 351 dfdra = t1+t2*(one-zeta) 352 dfdrb = t1-t2*(one+zeta) 353c 354 end 355