1c Uniform-gas correlation of Perdew and Wang 1991 2c 3c This has the same form as VWN functional V, only with a different 4c form for the parameterized functionals of rs. The VWN V code is 5c reused. 6* 7* $Id$ 8* 9#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 10 Subroutine xc_pw91lda(tol_rho, fac, lfac, nlfac, rho, Amat, nq, 11 & ipol, Ec, qwght, ldew, func) 12#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV) 13#include "dft2drv.fh" 14 Subroutine xc_pw91lda_d2(tol_rho, fac, lfac, nlfac, rho, 15 & Amat, Amat2, nq, ipol, Ec, qwght, ldew, func) 16#else 17#include "dft2drv.fh" 18#include "dft3drv.fh" 19 Subroutine xc_pw91lda_d3(tol_rho, fac, lfac, nlfac, rho, 20 & Amat, Amat2, Amat3, nq, ipol, Ec, qwght, ldew, func) 21#endif 22 implicit none 23c 24 integer nq, ipol 25 double precision fac, Ec, tol_rho 26 logical ldew, lfac, nlfac 27 double precision func(*) ! value of the functional [output] 28c 29c Charge Density 30c 31 double precision rho(nq,ipol*(ipol+1)/2) 32c 33c Quadrature Weights 34c 35 double precision qwght(nq) 36c 37c Partial Derivatives of the Correlation Energy Functional 38c 39 double precision Amat(nq,ipol) 40#ifdef SECOND_DERIV 41 double precision Amat2(nq,*) 42#endif 43#ifdef THIRD_DERIV 44 double precision Amat3(nq,NCOL_AMAT3) 45#endif 46c 47 double precision onethird, fourthirds, twothirds, pi 48 double precision fivethirds, seventhirds, threehalf 49 Parameter (onethird = 1.D0/3.D0, fourthirds = 4.D0/3.D0) 50 Parameter (twothirds = 2.D0/3.D0) 51 Parameter (fivethirds = 5.0d0/3.0d0) 52 Parameter (seventhirds = 7.0d0/3.0d0) 53 Parameter (threehalf = 3.0d0/2.0d0) 54 Parameter (pi = 3.1415926535898D0) 55c 56c Functional Parameters 57c 58 double precision A(3), alp(3), b(4,3) 59 save A, alp, b 60c 61 double precision e(3), d1e(3), rhoval, rs, d1rs, x, d1x, 62 & h1, d1h1, h2, d1h2, 63 & d1zeta(2), d1ersz(2), d1edrho(2), zeta, fz, d1fz, eps, 64 & dec_rs1, dec_rsz, d1dec_rs1, d1dec_rsz(2) 65 double precision devwn_rsz, d1devwn_rsz(2), zeta2, zeta3, zeta4, 66 & d2fz0, beta_rs1, d1beta_rs1, t_vwn, d1t_vwn 67#ifdef SECOND_DERIV 68 double precision d2beta_rs1, d2t_vwn, d2devwn_rsz(3) 69 double precision d2rs, d2x, d2h1, d2h2, 70 & d2e(3), d2zeta(3), d2dec_rs1, d2dec_rsz(3), 71 & d2ersz(3), d2edrho(3), d2fz, rrho2 72#endif 73 double precision p0, p1, p2, p3 74 double precision p4 75c 76#ifdef THIRD_DERIV 77 double precision d3beta_rs1, d3t_vwn, d3devwn_rsz(4) 78 double precision d3rs, d3x, d3h1, d3h2, 79 & d3e(3), d3zeta(4), d3dec_rs1, d3dec_rsz(4), 80 & d3ersz(4), d3edrho(4), d3fz, rrho3 81#endif 82c 83 integer i, n, initial 84 save initial 85c Daniel (10-19-12): Parameters are taken from the paper: 86c Phys. Rev. B 1992, 45, 13244. 87 data A / 0.0310907d0, 0.01554535d0, 0.0168869d0 / 88 data alp / 0.21370d0, 0.20548d0, 0.11125d0 / 89 data b / 7.5957d0, 3.5876d0, 1.6382d0, 0.49294d0, 90 & 14.1189d0, 6.1977d0, 3.3662d0, 0.62517d0, 91 & 10.357d0, 3.6231d0, 0.88026d0, 0.49671d0 / 92 data initial /1/ 93c 94c Define miscellaneous parameters. 95c 96 p0 = (1.0d0/(fourthirds*pi))**onethird 97 p1 = 0.5D0/(2.d0**onethird - 1.d0) 98 p2 = fourthirds*p1 99 p3 = onethird*p2 100c For XC-third derivative 101 p4 = -twothirds*p3 102 d2fz0 = 2.d0*p3 103 if (initial.eq.1)then 104 initial = 0 105c For convenience store -2A as A and multiply betas by 2A 106 do i = 1, 3 107 A(i) = -2d0*A(i) 108 do n = 1, 4 109 b(n,i) = -A(i)*b(n,i) 110 enddo 111 enddo 112c Finally, change the sign on A for spin stiffness since 113c the negative of that is fitted in the PW'91 paper. We can't 114c just take the negative of A at the start since A also contributes 115c to the argument of the ln function. 116 A(3) = -A(3) 117 endif 118c 119c ======> BOTH SPIN-RESTRICTED AND UNRESTRICTED <====== 120c 121 do 200 n = 1, nq 122 if (rho(n,1).lt.tol_rho)goto 200 123c 124 rhoval = rho(n,1) 125 rs = p0*rhoval**(-onethird) 126 d1rs = -onethird*rs/rhoval 127 x = sqrt(rs) 128 d1x = 0.5d0/x 129#ifdef SECOND_DERIV 130 d2rs = -fourthirds*d1rs/rhoval 131 d2x = -0.5d0*d1x/rs 132#endif 133#ifdef THIRD_DERIV 134 d3rs = -seventhirds*d2rs/rhoval 135 d3x = -threehalf*d2x/rs 136#endif 137c 138c Evaluate the individual correlation energy formulas 139c 140c Note that the Monte Carlo form (p = 1) is used for h2. 141c 142 do i = 1, 3 143 h2 = x*(b(1,i) + x*(b(2,i) + x*(b(3,i) + x*b(4,i)))) 144 d1h2 = b(1,i) 145 & + x*(2d0*b(2,i) + x*(3d0*b(3,i) + 4d0*x*b(4,i))) 146#ifdef SECOND_DERIV 147 d2h2 = 2d0*b(2,i) + x*(6d0*b(3,i) + 12d0*x*b(4,i)) 148#endif 149#ifdef THIRD_DERIV 150 d3h2 = 6.0d0*b(3,i) + 24.0d0*x*b(4,i) 151#endif 152c 153 h1 = DLOG(1d0+1d0/h2) 154 d1h1 = -d1h2/(h2*(h2+1d0)) 155#ifdef SECOND_DERIV 156 d2h1 = d1h1*d1h1*(2d0*h2+1d0) - d2h2/(h2*(h2+1d0)) 157#endif 158#ifdef THIRD_DERIV 159 d3h1 = 2d0*d2h1*d1h1*(2d0*h2+1d0) 160 1 + 2d0*d1h1*d1h1*d1h2 161 2 - d3h2/(h2*(h2+1d0)) 162 3 - d2h2*d1h1*(2d0*h2+1d0)/(h2*(h2+1d0)) 163#endif 164c 165 e(i) = A(i)*(1d0+alp(i)*rs)*h1 166 d1e(i) = A(i)*(2d0*alp(i)*x*h1+(1d0+alp(i)*rs)*d1h1) 167#ifdef SECOND_DERIV 168 d2e(i) = A(i)*(2d0*alp(i)*h1+4d0*alp(i)*x*d1h1 169 & +(1d0+alp(i)*rs)*d2h1) 170#endif 171#ifdef THIRD_DERIV 172 d3e(i) = A(i)*( 6d0*alp(i)*d1h1 + 6d0*alp(i)*x*d2h1 173 1 + (1d0+alp(i)*rs)*d3h1 ) 174#endif 175c 176c Transform derivatives wrt x to derivatives wrt rs 177c 178#ifdef THIRD_DERIV 179 d3e(i) = d3e(i)*d1x*d1x*d1x + 3.0d0*d2e(i)*d1x*d2x 180 & + d1e(i)*d3x 181#endif 182#ifdef SECOND_DERIV 183c Do 2nd derivative first so the x first derivative in d1e 184c is not lost 185 d2e(i) = d2e(i)*d1x*d1x + d1e(i)*d2x 186#endif 187 d1e(i) = d1e(i)*d1x 188 enddo 189c 190c Compute the polarization function and its derivatives 191c 192 if (ipol.eq.1) then 193 zeta = 0.0d0 194 else 195 zeta = (rho(n,2) - rho(n,3))/rhoval 196 endif 197 if (zeta.gt.1.d0)then 198 zeta = 1.d0 199 elseif (zeta.lt.-1.d0)then 200 zeta =-1.d0 201 endif 202 fz = ((1.d0+zeta)**fourthirds + 203 & (1.d0-zeta)**fourthirds - 2.d0)*p1 204 d1fz = ((1.d0+zeta)**onethird - 205 & (1.d0-zeta)**onethird)*p2 206 d1zeta(1) = (1.d0-zeta)/rhoval 207 d1zeta(2) =-(1.d0+zeta)/rhoval 208#ifdef SECOND_DERIV 209 if(dabs(zeta).lt.tol_rho) then 210 d2fz = d2fz0 211 else 212 if (1.0d0+zeta.le.tol_rho) then 213 d2fz = ((1.d0-zeta)**(-twothirds))*p3 214 else if (1.0d0-zeta.le.tol_rho) then 215 d2fz = ((1.d0+zeta)**(-twothirds))*p3 216 else 217 d2fz = ((1.d0+zeta)**(-twothirds) + 218 & (1.d0-zeta)**(-twothirds))*p3 219 endif 220 endif 221 rrho2 = 2.d0/(rhoval*rhoval) 222c 1 = aa, 2 = ab, 3 = bb 223 d2zeta(1) =-rrho2*(1.d0-zeta) 224 d2zeta(2) = rrho2*zeta 225 d2zeta(3) = rrho2*(1.d0+zeta) 226#endif 227#ifdef THIRD_DERIV 228 if (dabs(zeta+1.0d0).le.tol_rho) then 229 d3fz = (-(1.0d0-zeta)**(-fivethirds))*p4 230 else if (dabs(zeta-1.0d0).le.tol_rho) then 231 d3fz = ((1.0d0+zeta)**(-fivethirds))*p4 232 else 233 d3fz = ((1.0d0+zeta)**(-fivethirds) - 234 & (1.0d0-zeta)**(-fivethirds))*p4 235 end if 236 rrho3 = 1.0d0/(rhoval*rhoval*rhoval) 237c 1 = aaa, 2 = aab, 3 = abb, 4 = bbb 238 d3zeta(1) = 6.0d0*(1.0d0-zeta)*rrho3 239 d3zeta(2) = 2.0d0*(1.0d0-3.0d0*zeta)*rrho3 240 d3zeta(3) = -2.0d0*(1.0d0+3.0d0*zeta)*rrho3 241 d3zeta(4) = -6.0d0*(1.0d0+3.0d0*zeta)*rrho3 242#endif 243c 244 dec_rs1 = e(2)-e(1) 245 d1dec_rs1 = d1e(2)-d1e(1) 246#ifdef SECOND_DERIV 247 d2dec_rs1 = d2e(2)-d2e(1) 248#endif 249#ifdef THIRD_DERIV 250 d3dec_rs1 = d3e(2)-d3e(1) 251#endif 252c 253 beta_rs1 = e(2)-e(1) 254 d1beta_rs1 = d1e(2)-d1e(1) 255 zeta2 = zeta*zeta 256 zeta3 = zeta2*zeta 257 zeta4 = zeta3*zeta 258 t_vwn = d2fz0*beta_rs1-e(3) 259 d1t_vwn = d2fz0*d1beta_rs1-d1e(3) 260 devwn_rsz = fz/d2fz0*(e(3)+t_vwn*zeta4) 261 d1devwn_rsz(1) = fz/d2fz0*(d1e(3)+d1t_vwn*zeta4) 262 d1devwn_rsz(2) = d1fz/d2fz0*(e(3)+t_vwn*zeta4) 263 & + fz/d2fz0*t_vwn*4.d0*zeta3 264#ifdef SECOND_DERIV 265 d2beta_rs1 = d2e(2)-d2e(1) 266 d2t_vwn = d2fz0*d2beta_rs1-d2e(3) 267 d2devwn_rsz(1) = fz/d2fz0*(d2e(3)+d2t_vwn*zeta4) 268 d2devwn_rsz(2) = d1fz/d2fz0*(d1e(3)+d1t_vwn*zeta4) 269 & + fz/d2fz0*d1t_vwn*4.d0*zeta3 270 d2devwn_rsz(3) = d2fz/d2fz0*(e(3)+t_vwn*zeta4) 271 & + d1fz/d2fz0*t_vwn*8.d0*zeta3 272 & + fz/d2fz0*t_vwn*12.d0*zeta2 273#endif 274#ifdef THIRD_DERIV 275 d3beta_rs1 = d3e(2)-d3e(1) 276 d3t_vwn = d2fz0*d3beta_rs1-d3e(3) 277c Derivatives: 1 = drsdrsdrs, 2 = drsdrsdzeta, 3 = drsdzetadzeta, 278c 4 = dzetadzetadzeta 279 d3devwn_rsz(1) = fz/d2fz0*(d3e(3)+d3t_vwn*zeta4) 280 d3devwn_rsz(2) = d1fz/d2fz0*(d2e(3)+d2t_vwn*zeta4) 281 & + fz/d2fz0*d2t_vwn*4.0d0*zeta3 282 d3devwn_rsz(3) = d2fz/d2fz0*(d1e(3)+d1t_vwn*zeta4) 283 & + d1fz/d2fz0*d1t_vwn*8.0d0*zeta3 284 & + fz/d2fz0*d1t_vwn*12.0d0*zeta2 285 d3devwn_rsz(4) = d3fz/d2fz0*(e(3)+t_vwn*zeta4) 286 & + d2fz/d2fz0*t_vwn*12.0d0*zeta3 287 & + d1fz/d2fz0*t_vwn*36.0d0*zeta2 288 & + fz/d2fz0*t_vwn*24.0d0*zeta 289#endif 290c 291c Compute the function deltaEc(rs,zeta) function and its derivatives 292c wrt rs and zeta for the spin-unrestricted case - the rest has the 293c same form for all VWN functionals and is handled in the header files. 294c 295 dec_rsz = devwn_rsz 296 d1dec_rsz(1) = d1devwn_rsz(1) 297 d1dec_rsz(2) = d1devwn_rsz(2) 298#ifdef SECOND_DERIV 299 d2dec_rsz(1) = d2devwn_rsz(1) 300 d2dec_rsz(2) = d2devwn_rsz(2) 301 d2dec_rsz(3) = d2devwn_rsz(3) 302#endif 303#ifdef THIRD_DERIV 304 d3dec_rsz(1) = d3devwn_rsz(1) 305 d3dec_rsz(2) = d3devwn_rsz(2) 306 d3dec_rsz(3) = d3devwn_rsz(3) 307 d3dec_rsz(4) = d3devwn_rsz(4) 308#endif 309c 310c Finish off the unrestricted case: 311c Assemble the entire functional and its derivatives given the 312c parameterization-dependent part deltaEc(rs,zeta) and its derivatives 313c 314 eps = e(1) + dec_rsz 315 d1ersz(1) = d1e(1) + d1dec_rsz(1) 316 d1ersz(2) = d1dec_rsz(2) 317 d1edrho(1) = d1ersz(1)*d1rs + d1ersz(2)*d1zeta(1) 318 d1edrho(2) = d1ersz(1)*d1rs + d1ersz(2)*d1zeta(2) 319 Ec = Ec + eps*qwght(n)*rhoval*fac 320 if (ldew) func(n) = func(n) + eps*rhoval*fac 321 Amat(n,1) = Amat(n,1) + (eps + rhoval*d1edrho(1))*fac 322 if (ipol.eq.2) 323 & Amat(n,2) = Amat(n,2) + (eps + rhoval*d1edrho(2))*fac 324#ifdef SECOND_DERIV 325c 1 = rsrs, 2 = rsz, 3 = zz 326 d2ersz(1) = d2e(1) + d2dec_rsz(1) 327 d2ersz(2) = d2dec_rsz(2) 328 d2ersz(3) = d2dec_rsz(3) 329c 1 = aa, 2 = ab, 3 = bb 330 d2edrho(1) = d2ersz(1)*d1rs*d1rs 331 & + d2ersz(2)*d1rs*d1zeta(1)*2.d0 332 & + d2ersz(3)*d1zeta(1)*d1zeta(1) 333 & + d1ersz(1)*d2rs 334 & + d1ersz(2)*d2zeta(1) 335 d2edrho(2) = d2ersz(1)*d1rs*d1rs 336 & + d2ersz(2)*d1rs*(d1zeta(1)+d1zeta(2)) 337 & + d2ersz(3)*d1zeta(1)*d1zeta(2) 338 & + d1ersz(1)*d2rs 339 & + d1ersz(2)*d2zeta(2) 340 d2edrho(3) = d2ersz(1)*d1rs*d1rs 341 & + d2ersz(2)*d1rs*d1zeta(2)*2.d0 342 & + d2ersz(3)*d1zeta(2)*d1zeta(2) 343 & + d1ersz(1)*d2rs 344 & + d1ersz(2)*d2zeta(3) 345 Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) 346 & + (2.d0*d1edrho(1) + rhoval*d2edrho(1))*fac 347 Amat2(n,D2_RA_RB) = Amat2(n,D2_RA_RB) 348 & + (d1edrho(1) + d1edrho(2) + rhoval*d2edrho(2))*fac 349 if (ipol.eq.2) 350 & Amat2(n,D2_RB_RB) = Amat2(n,D2_RB_RB) 351 & + (2.d0*d1edrho(2) + rhoval*d2edrho(3))*fac 352#endif 353#ifdef THIRD_DERIV 354c 1 = rsrsrs, 2 = rsrsz, 3 = rszz, 4 = zzz 355 d3ersz(1) = d3e(1) + d3dec_rsz(1) 356 d3ersz(2) = d3dec_rsz(2) 357 d3ersz(3) = d3dec_rsz(3) 358 d3ersz(4) = d3dec_rsz(4) 359c 1 = aaa, 2 = aab, 3 = abb, 4 = bbb 360 d3edrho(1) = d3ersz(1)*d1rs*d1rs*d1rs 361 & + d2ersz(1)*d1rs*d2rs*3.0d0 362 & + d3ersz(3)*d1rs*d1zeta(1)*d1zeta(1)*3.0d0 363 & + d2ersz(2)*d1rs*d2zeta(1)*3.0d0 364 & + d1ersz(1)*d3rs 365 & + d2ersz(2)*d1zeta(1)*d2rs*3.0d0 366 & + d3ersz(2)*d1zeta(1)*d1rs*d1rs*3.0d0 367 & + d3ersz(4)*d1zeta(1)*d1zeta(1)*d1zeta(1) 368 & + d2ersz(3)*d1zeta(1)*d2zeta(1)*3.0d0 369 & + d1ersz(2)*d3zeta(1) 370 d3edrho(2) = d3ersz(1)*d1rs*d1rs*d1rs 371 & + d2ersz(1)*d1rs*d2rs*3.0d0 372 & + d3ersz(3)*d1rs*(d1zeta(1)*d1zeta(1) 373 & + d1zeta(1)*d1zeta(2)*2.0d0) 374 & + d2ersz(2)*d1rs*(d2zeta(2)*2.0d0 375 & + d2zeta(1)) 376 & + d1ersz(1)*d3rs 377 & + d2ersz(2)*d2rs*(d1zeta(1)*2.0d0 378 & + d1zeta(2)) 379 & + d3ersz(2)*d1rs*d1rs*(d1zeta(2) 380 & + d1zeta(1)*2.0d0) 381 & + d3ersz(4)*d1zeta(2)*d1zeta(1)*d1zeta(1) 382 & + d2ersz(3)*(d1zeta(1)*d2zeta(2)*2.0d0 383 & + d1zeta(2)*d2zeta(1)) 384 & + d1ersz(2)*d3zeta(2) 385 d3edrho(3) = d3ersz(1)*d1rs*d1rs*d1rs 386 & + d2ersz(1)*d1rs*d2rs*3.0d0 387 & + d3ersz(3)*d1rs*(d1zeta(2)*d1zeta(2) 388 & + d1zeta(2)*d1zeta(1)*2.0d0) 389 & + d2ersz(2)*d1rs*(d2zeta(2)*2.0d0 390 & + d2zeta(3)) 391 & + d1ersz(1)*d3rs 392 & + d2ersz(2)*d2rs*(d1zeta(2)*2.0d0 393 & + d1zeta(1)) 394 & + d3ersz(2)*d1rs*d1rs*(d1zeta(1) 395 & + d1zeta(2)*2.0d0) 396 & + d3ersz(4)*d1zeta(1)*d1zeta(2)*d1zeta(2) 397 & + d2ersz(3)*(d1zeta(2)*d2zeta(2)*2.0d0 398 & + d1zeta(1)*d2zeta(3)) 399 & + d1ersz(2)*d3zeta(3) 400 d3edrho(4) = d3ersz(1)*d1rs*d1rs*d1rs 401 & + d2ersz(1)*d1rs*d2rs*3.0d0 402 & + d3ersz(3)*d1rs*d1zeta(2)*d1zeta(2)*3.0d0 403 & + d2ersz(2)*d1rs*d2zeta(3)*3.0d0 404 & + d1ersz(1)*d3rs 405 & + d2ersz(2)*d1zeta(2)*d2rs*3.0d0 406 & + d3ersz(2)*d1zeta(2)*d1rs*d1rs*3.0d0 407 & + d3ersz(4)*d1zeta(2)*d1zeta(2)*d1zeta(2) 408 & + d2ersz(3)*d1zeta(2)*d2zeta(3)*3.0d0 409 & + d1ersz(2)*d3zeta(4) 410c 411 Amat3(n,D3_RA_RA_RA) = Amat3(n,D3_RA_RA_RA) 412 & + (3.0d0*d2edrho(1) + rhoval*d3edrho(1))*fac 413 Amat3(n,D3_RA_RA_RB) = Amat3(n,D3_RA_RA_RB) 414 & + (d2edrho(1) + 2.0d0*d2edrho(2) + rhoval*d3edrho(2))*fac 415 Amat3(n,D3_RA_RB_RB) = Amat3(n,D3_RA_RB_RB) 416 & + (2.0d0*d2edrho(2) + d2edrho(3) + rhoval*d3edrho(3))*fac 417 if (ipol.eq.2) 418 & Amat3(n,D3_RB_RB_RB) = Amat3(n,D3_RB_RB_RB) 419 & + (3.0d0*d2edrho(3) + rhoval*d3edrho(4))*fac 420#endif 421 200 continue 422c 423 return 424 end 425c 426#ifndef SECOND_DERIV 427#define SECOND_DERIV 428c 429c Compile source again for the 2nd derivative case 430c 431#include "xc_pw91lda.F" 432#endif 433c 434#ifndef THIRD_DERIV 435#define THIRD_DERIV 436c 437c Compile source again for the 3rd derivative case 438c 439#include "xc_pw91lda.F" 440#endif 441