1#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 2#if !defined(NWAD_PRINT) 3C> \ingroup nwxc 4C> @{ 5C> 6C> \file nwxc_c_perdew86.F 7C> The Perdew correlation functional of 1986 8C> 9C> @} 10#endif 11#endif 12C> \ingroup nwxc_priv 13C> @{ 14C> 15C> \brief Evaluate the Perdew 1986 correlation functional 16C> 17C> Evaluates the Perdew 1986 GGA correlation functional [1,2,3]. 18C> 19C> ### References ### 20C> 21C> [1] J.P. Perdew, 22C> "Density-functional approximation for the correlation energy of 23C> the inhomogeneous electron gas", Phys. Rev. B <b>33</b>, 24C> 8822–8824 (1986), DOI: 25C> <a href="https://doi.org/10.1103/PhysRevB.33.8822"> 26C> 10.1103/PhysRevB.33.8822</a>. 27C> 28C> [2] P. Mlynarski, D.R. Salahub, 29C> "Self-consistent implementation of nonlocal exchange and 30C> correlation in a Gaussian density-functional method", 31C> Phys. Rev. B <b>43</b>, 1399–1410 (1991), DOI: 32C> <a href="https://doi.org/10.1103/PhysRevB.43.1399"> 33C> 10.1103/PhysRevB.43.1399</a>. 34C> 35C> [3] J.P. Perdew, 36C> "Erratum: Density-functional approximation for the correlation 37C> energy of the inhomogeneous electron gas", Phys. Rev. B 38C> <b>34</b>, 7406–7406 (1986), DOI: 39C> <a href="https://doi.org/10.1103/PhysRevB.34.7406"> 40C> 10.1103/PhysRevB.34.7406</a>. 41C> 42#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 43#if defined(NWAD_PRINT) 44 Subroutine nwxc_c_perdew86_p(tol_rho, ipol, nq, wght, rho, rgamma, 45 & ffunc) 46#else 47 Subroutine nwxc_c_perdew86(tol_rho, ipol, nq, wght, rho, rgamma, 48 & ffunc) 49#endif 50#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV) 51 Subroutine nwxc_c_perdew86_d2(tol_rho, ipol, nq, wght, 52 & rho, rgamma, ffunc) 53#else 54 Subroutine nwxc_c_perdew86_d3(tol_rho, ipol, nq, wght, 55 & rho, rgamma, ffunc) 56#endif 57c 58c$Id$ 59c 60#include "nwad.fh" 61c 62 implicit none 63c 64#include "nwxc_param.fh" 65c 66c Input and other parameters 67c 68 double precision tol_rho !< [Input] The lower limit on the density 69 integer ipol !< [Input] The number of spin channels 70 integer nq !< [Input] The number of points 71 double precision wght !< [Input] The weight of the functional 72c 73c Charge Density 74c 75 type(nwad_dble)::rho(nq,*) !< [Input] The density 76c 77c Charge Density Gradient 78c 79 type(nwad_dble)::rgamma(nq,*) !< [Input] The norm of the density gradients 80c 81c Sampling Matrices for the XC Potential 82c 83 type(nwad_dble)::ffunc(nq) !< [Output] The value of the functional 84c double precision Amat(nq,*) !< [Output] The derivative wrt rho 85c double precision Cmat(nq,*) !< [Output] The derivative wrt rgamma 86#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 87c 88c Sampling Matrices for the XC Kernel 89c 90c double precision Amat2(nq,*) !< [Output] The 2nd derivative wrt rho 91c double precision Cmat2(nq,*) !< [Output] The 2nd derivative wrt rgamma 92c !< and possibly rho 93#endif 94#if defined(THIRD_DERIV) 95c 96c Sampling Matrices for the XC Kernel 97c 98c double precision Amat3(nq,*) !< [Output] The 3rd derivative wrt rho 99c double precision Cmat3(nq,*) !< [Output] The 3rd derivative wrt rgamma 100c !< and possibly rho 101#endif 102 double precision TOLL, EXPTOL, alpha, beta, pgamma, delta, 103 & beta10, ftilde, zzz, fff, pfff, CINF, ONE, 104 & ONE3, THREE, FOUR3, SEV6, FIVE3, 105 & TWO3, FIVE6, pi 106 double precision SEVEN3, EIGHT3 107 Parameter (TOLL = 1.D-40, EXPTOL = 80.d0) 108 Parameter (alpha = 0.023266D0, beta = 7.389D-6, 109 & pgamma = 8.723d0, delta = 0.472d0, beta10 = 10000.d0*beta) 110 parameter (ftilde = 0.11d0, zzz = 0.001667d0, fff = 0.002568d0) 111 parameter(pfff = 1.745d0, CINF = zzz+fff) 112 Parameter (ONE = 1.D0, ONE3 = 1.d0/3.d0, THREE = 3.d0) 113 Parameter (FOUR3 = 4.D0/3.D0, SEV6 = 7.d0/6.d0) 114 parameter (FIVE3 = 5.d0/3.d0, TWO3 = 2.d0/3.d0, FIVE6 = 5.d0/6.d0) 115 parameter (SEVEN3 = 7.0d0/3.0d0, EIGHT3 = 8.0d0/3.0d0) 116c parameter (pi = 3.1415926535897932385d0) 117c 118c Mlynarski Salahub PRB 43, 1399 (1991) 119c 120 integer n 121 type(nwad_dble)::rhoval 122 double precision rsfact 123 type(nwad_dble)::rs, rs2, rs3 124 type(nwad_dble)::rho13, rho43, rho76, arho 125 double precision d1rs 126#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 127c double precision d2rs 128#endif 129#if defined(THIRD_DERIV) 130c double precision d3rs 131#endif 132 type(nwad_dble)::gamma,gam12,zeta,func,dt12,phi,d,expfac 133 type(nwad_dble)::anum,aden,Cn,dm1 134 double precision d1anum, d1aden, d1Cn, 135 & d1phi(2), dlnphi, d1f(3), 136 & dlnfrho(2), dlnfgam 137 double precision d1z(2), adp, d1d(2), t, 138 & d1dt12 139 double precision aden2 140#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 141c double precision d2anum, d2aden, rrho2, d2z(3), dpp, d2d(3), 142c & d2phi(3), d2dt12, d2Cn 143c double precision aden3 144c double precision arho2 145c double precision d2lnphi 146c double precision d2f(3) 147c double precision d2lnfrho(3), d2lnfrg(2), d2lnfgam 148#endif 149#if defined(THIRD_DERIV) 150c double precision d3lnphi 151c double precision d3anum, d3aden, d3Cn, d3phi(4) 152c double precision d3lnfrho(4), d3lnfgam 153c double precision d3f(3) 154c double precision aden4 155c double precision arho3 156#endif 157c 158 pi = acos(-1.0d0) 159 rsfact = (0.75d0/pi)**ONE3 160c 161 if (ipol.eq.1 )then 162c 163c ======> SPIN-RESTRICTED <====== 164c 165 do 10 n = 1, nq 166 rhoval = rho(n,R_T) 167 if (rhoval.lt.tol_rho) goto 10 168 arho=1.d0/rhoval 169 rho13 = rhoval**ONE3 170 rho43 = rhoval*rho13 171 rho76 = rhoval**SEV6 172 rs = rsfact/rho13 173 rs2 = rs*rs 174 rs3 = rs2*rs 175c d1rs = -ONE3*rs*arho 176#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 177c d2rs = -FOUR3*d1rs*arho 178#endif 179#if defined(THIRD_DERIV) 180c d3rs = -SEVEN3*d2rs*arho 181#endif 182 gamma = rgamma(n,G_TT) 183c gamma = delrho(n,1,1)*delrho(n,1,1) + 184c & delrho(n,2,1)*delrho(n,2,1) + 185c & delrho(n,3,1)*delrho(n,3,1) 186 if (gamma.gt.tol_rho*tol_rho) then 187 gam12 = sqrt(gamma) 188 else 189 gam12 = 0.0d0 190 endif 191c 192c C(n) 193c 194 anum = fff+alpha*rs+beta*rs2 195 aden = 1.d0+pgamma*rs+delta*rs2+beta10*rs3 196 Cn = zzz + anum/aden 197c d1anum = alpha + 2d0*beta*rs 198c d1aden = pgamma + 2d0*delta*rs + 3d0*beta10*rs2 199#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 200c d2anum = 2d0*beta 201c d2aden = 2d0*delta + 6d0*beta10*rs 202#endif 203#if defined(THIRD_DERIV) 204c d3anum = 0.0d0 205c d3aden = 6.0d0*beta10 206#endif 207c First compute rs derivative 208c aden2 = aden*aden 209c d1Cn = d1anum/aden - anum*d1aden/aden2 210#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 211c aden3 = aden2*aden 212c d2Cn = d2anum/aden - (2d0*d1anum*d1aden+anum*d2aden)/aden2 213c & + 2d0*anum*d1aden**2/aden3 214#endif 215#if defined(THIRD_DERIV) 216c aden4 = aden3*aden 217c 218c d3Cn = -( 3.0d0*d2anum*d1aden + 3.0d0*d1anum*d2aden 219c 1 + anum*d3aden )/aden2 220c 2 + 6.0d0*( d1anum*d1aden**2 221c 3 + anum*d2aden*d1aden )/aden3 222c 4 - 6.0d0*anum*d1aden**3/aden4 223#endif 224c Convert to rho derivative 225#if defined(THIRD_DERIV) 226c d3Cn = d3Cn*d1rs*d1rs*d1rs 227c 1 + 3.0d0*d2Cn*d2rs*d1rs 228c 2 + d1Cn*d3rs 229#endif 230#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 231c d2Cn = d2Cn*d1rs*d1rs + d1Cn*d2rs 232#endif 233c d1Cn = d1Cn*d1rs 234c 235c phi(n,gradn) 236c 237 expfac = 0.d0 238 phi = (pfff*ftilde)*(CINF/Cn)*gam12/rho76 239 if (phi.lt.EXPTOL) expfac = exp(-phi) 240c dlnphi = -(d1Cn/Cn + SEV6/rhoval) 241c d1phi(1) = phi*dlnphi 242#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 243c arho2 = arho*arho 244c d2lnphi = (d1Cn/Cn)**2 - d2Cn/Cn + SEV6*arho2 245c d2phi(1) = d1phi(1)*dlnphi + phi*d2lnphi 246c d2phi(1) = d1phi(1)*dlnphi 247c & + phi*((d1Cn/Cn)**2 - d2Cn/Cn + SEV6/rhoval**2) 248#endif 249#if defined(THIRD_DERIV) 250c arho3 = arho2*arho 251c 252c d3lnphi = -2.0d0*(d1Cn/Cn)**3 253c 1 + 3.0d0*(d2Cn/Cn)*(d1Cn/Cn) 254c 2 - d3Cn/Cn 255c 3 - SEVEN3*arho3 256c d3phi(1) = d2phi(1)*dlnphi 257c 1 + 2.0d0*d1phi(1)*d2lnphi 258c 2 + phi*d3lnphi 259#endif 260c 261c functional 262c 263 func = expfac*Cn*gamma/rho43 264c dlnfrho(1) = d1Cn/Cn - (d1phi(1) + FOUR3/rhoval) 265c d1f(1) = dlnfrho(1)*func 266c Amat(n,D1_RA) = Amat(n,D1_RA) + d1f(1)*wght 267c if (gam12.gt.TOLL)then 268c d1phi(2) = phi / (2d0*gamma) 269c dlnfgam = 1d0/gamma - d1phi(2) 270c d1f(3) = func*dlnfgam 271c Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + d1f(3)*wght 272c Cmat(n,D1_GAB) = Cmat(n,D1_GAB) + d1f(3)*2D0*wght 273#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 274c d2phi(2) = d1phi(2)*dlnphi 275c d2phi(3) =-d1phi(2)/(2d0*gamma) 276c!!! Which of the following are actually needed for restricted? 277c!!! Should treat derivatives of d as zero? d is a constant? 278c Daniel (11-19-12): d is a constant (it equals 1) for a restricted 279c calculation, since there is no spin-polarization. Thus, the 280c derivatives are zero. 281c d2lnfrho(1) = -d2phi(1) - (d1Cn/Cn)**2 + d2Cn/Cn 282c 1 + FOUR3*arho2 283c 284c d2f(1) = d1f(1)*dlnfrho(1) 285c 1 + func*d2lnfrho(1) 286c 287c t = d2f(1)*wght 288c 289c Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) + t 290c Amat2(n,D2_RA_RB) = Amat2(n,D2_RA_RB) + t 291c & + (d1f(1)*dlnfrho(1) 292c & + func*t)*wght 293#if 0 294c Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) 295c & + (d1f(1)*dlnfrho(1) 296c & + func*(d1d(1)*d1d(1)*dm1**2-d2d(1)*dm1+t))*wght 297c Amat2(n,D2_RA_RB) = Amat2(n,D2_RA_RB) 298c & + (d1f(1)*dlnfrho(2) 299c & + func*(d1d(1)*d1d(2)*dm1**2-d2d(2)*dm1+t))*wght 300#endif 301c rg terms 302c d2lnfrg(1) = -d2phi(2) 303c d2f(2) = (d1f(1)*dlnfgam + func*d2lnfrg(1)) 304c t = d2f(2)*wght 305c 306c Cmat2(n,D2_RA_GAA) = Cmat2(n,D2_RA_GAA) + t 307c Cmat2(n,D2_RA_GAB) = Cmat2(n,D2_RA_GAB) + t*2d0 308c Cmat2(n,D2_RA_GBB) = Cmat2(n,D2_RA_GBB) + t 309c gg terms 310c d2lnfgam = -1.0d0/gamma**2 - d2phi(3) 311c d2f(3) = d1f(3)*dlnfgam + func*d2lnfgam 312c t = d2f(3)*wght 313c 314c Cmat2(n,D2_GAA_GAA) = Cmat2(n,D2_GAA_GAA) + t 315c Cmat2(n,D2_GAA_GBB) = Cmat2(n,D2_GAA_GBB) + t 316c Cmat2(n,D2_GAA_GAB) = Cmat2(n,D2_GAA_GAB) + t*2d0 317c Cmat2(n,D2_GAB_GAB) = Cmat2(n,D2_GAB_GAB) + t*4d0 318#endif 319#if defined(THIRD_DERIV) 320c rrr terms 321c d3lnfrho(1) = -d3phi(1) 322c 1 + 2.0d0*(d1Cn/Cn)**3 323c 2 - 3.0d0*(d2Cn/Cn)*(d1Cn/Cn) 324c 3 + d3Cn/Cn 325c 4 - EIGHT3*arho3 326c 327c d3f(1) = d2f(1)*dlnfrho(1) 328c 1 + 2.0d0*d1f(1)*d2lnfrho(1) 329c 2 + func*d3lnfrho(1) 330c 331c t = d3f(1)*wght 332c 333c Amat3(n,D3_RA_RA_RA) = Amat3(n,D3_RA_RA_RA) + t 334c Amat3(n,D3_RA_RA_RB) = Amat3(n,D3_RA_RA_RB) + t 335c Amat3(n,D3_RA_RB_RB) = Amat3(n,D3_RA_RB_RB) + t 336c rrg terms 337c d3phi(2) = d2phi(2)*dlnphi + d1phi(2)*d2lnphi 338c 339c t = ( d2f(2)*dlnfrho(1) 340c 1 - d1f(1)*d2phi(2) 341c 2 + d1f(3)*d2lnfrho(1) 342c 3 - func*d3phi(2) )*wght 343c 344c Cmat3(n,D3_RA_RA_GAA) = Cmat3(n,D3_RA_RA_GAA) + t 345c Cmat3(n,D3_RA_RA_GAB) = Cmat3(n,D3_RA_RA_GAB) + t*2.0d0 346c Cmat3(n,D3_RA_RA_GBB) = Cmat3(n,D3_RA_RA_GBB) + t 347c Cmat3(n,D3_RA_RB_GAA) = Cmat3(n,D3_RA_RB_GAA) + t 348c Cmat3(n,D3_RA_RB_GAB) = Cmat3(n,D3_RA_RB_GAB) + t*2.0d0 349c Cmat3(n,D3_RA_RB_GBB) = Cmat3(n,D3_RA_RB_GBB) + t 350c rgg terms 351c d3phi(3) = -d2phi(3)*dlnphi 352c 353c t = ( d2f(2)*dlnfgam 354c 1 + d1f(1)*d2lnfgam 355c 2 + d1f(3)*d2lnfrg(1) 356c 3 + func*d3phi(3) )*wght 357c 358c Cmat3(n,D3_RA_GAA_GAA) = Cmat3(n,D3_RA_GAA_GAA) + t 359c Cmat3(n,D3_RA_GAA_GAB) = Cmat3(n,D3_RA_GAA_GAB) + t*2.0d0 360c Cmat3(n,D3_RA_GAA_GBB) = Cmat3(n,D3_RA_GAA_GBB) + t 361c Cmat3(n,D3_RA_GAB_GAB) = Cmat3(n,D3_RA_GAB_GAB) + t*4.0d0 362c Cmat3(n,D3_RA_GAB_GBB) = Cmat3(n,D3_RA_GAB_GBB) + t*2.0d0 363c Cmat3(n,D3_RA_GBB_GBB) = Cmat3(n,D3_RA_GBB_GBB) + t 364c ggg terms 365c d3phi(4) = -3.0d0*d2phi(3)/(2.0d0*gamma) 366c d3lnfgam = 2.0d0/gamma**3 - d3phi(4) 367c 368c t = ( d2f(3)*dlnfgam 369c 1 + 2.0d0*d1f(3)*d2lnfgam 370c 2 + func*d3lnfgam )*wght 371c 372c Cmat3(n,D3_GAA_GAA_GAA) = Cmat3(n,D3_GAA_GAA_GAA) + t 373c Cmat3(n,D3_GAA_GAA_GAB) = Cmat3(n,D3_GAA_GAA_GAB) 374c 1 + t*2.0d0 375c Cmat3(n,D3_GAA_GAA_GBB) = Cmat3(n,D3_GAA_GAA_GBB) + t 376c Cmat3(n,D3_GAA_GAB_GAB) = Cmat3(n,D3_GAA_GAB_GAB) 377c 1 + t*4.0d0 378c Cmat3(n,D3_GAA_GAB_GBB) = Cmat3(n,D3_GAA_GAB_GBB) 379c 1 + t*2.0d0 380c Cmat3(n,D3_GAA_GBB_GBB) = Cmat3(n,D3_GAA_GBB_GBB) + t 381c Cmat3(n,D3_GAB_GAB_GAB) = Cmat3(n,D3_GAB_GAB_GAB) 382c 1 + t*8.0d0 383#endif 384c endif 385 ffunc(n)=ffunc(n)+func*wght 386 10 continue 387 else 388c 389c ======> SPIN-UNRESTRICTED <====== 390c 391 do 20 n = 1, nq 392 rhoval = 0.0d0 393 gamma = 0.0d0 394 if (rho(n,R_A).ge.0.5d0*tol_rho) then 395 rhoval = rhoval + rho(n,R_A) 396 gamma = gamma + rgamma(n,G_AA) 397 endif 398 if (rho(n,R_B).ge.0.5d0*tol_rho) then 399 rhoval = rhoval + rho(n,R_B) 400 gamma = gamma + rgamma(n,G_BB) 401 if (rho(n,R_A).ge.0.5d0*tol_rho) then 402 gamma = gamma + 2.0d0*rgamma(n,G_AB) 403 endif 404 endif 405 if (rhoval.lt.tol_rho) goto 20 406 arho=1.d0/rhoval 407 rho13 = abs(rhoval)**ONE3 408 rho43 = rhoval*rho13 409 rho76 = abs(rhoval)**SEV6 410 rs = rsfact/rho13 411 rs2 = rs*rs 412 rs3 = rs2*rs 413c d1rs = -ONE3*rs*arho 414#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 415c d2rs = -FOUR3*d1rs*arho 416#endif 417#if defined(THIRD_DERIV) 418c d3rs = -SEVEN3*d2rs*arho 419#endif 420c gamma = rgamma(n,G_AA)+rgamma(n,G_BB)+2.0d0*rgamma(n,G_AB) 421c gamma = delrho(n,1,1)*delrho(n,1,1) + 422c & delrho(n,2,1)*delrho(n,2,1) + 423c & delrho(n,3,1)*delrho(n,3,1) + 424c & delrho(n,1,2)*delrho(n,1,2) + 425c & delrho(n,2,2)*delrho(n,2,2) + 426c & delrho(n,3,2)*delrho(n,3,2) + 427c & 2.d0*(delrho(n,1,1)*delrho(n,1,2) + 428c & delrho(n,2,1)*delrho(n,2,2) + 429c & delrho(n,3,1)*delrho(n,3,2)) 430 if (gamma.gt.tol_rho*tol_rho) then 431 gam12 = sqrt(gamma) 432 else 433 gam12 = 0.0d0 434 endif 435 zeta = (rho(n,R_A) - rho(n,R_B))*arho 436 if(zeta.le.-1d0) zeta=-1d0 437 if(zeta.ge.1d0) zeta=1d0 438c d1z(1) = (1.d0 - zeta)*arho 439c d1z(2) = -(1.d0 + zeta)*arho 440#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 441c rrho2 = 2.d0*arho*arho 442c 1 = aa, 2 = ab, 3 = bb 443c d2z(1) =-rrho2*(1.d0-zeta) 444c d2z(2) = rrho2*zeta 445c d2z(3) = rrho2*(1.d0+zeta) 446#endif 447#if defined(THIRD_DERIV) 448c d3rs = -SEVEN3*d2rs*arho 449c if ((1.d0-zeta).lt.tol_rho) then 450c d3fz = (1.d0+zeta)**(-FIVE3) 451c else if ((1.d0+zeta).lt.tol_rho) then 452c d3fz = (1.d0-zeta)**(-FIVE3) 453c else 454c d3fz = (1.d0+zeta)**(-FIVE3) + (1.d0-zeta)**(-FIVE3) 455c end if 456c d3fz = -d3fz*TWO3*ONE3*FOUR3/(2.d0**FOUR3-2.d0) 457c 458c rrho3 = rrho2*arho 459c 460c 1 = aaa, 2 = aab, 3 = abb, 4 = bbb 461c d3z(1) = 3.0d0*rrho3*(1.0d0 - zeta) 462c d3z(2) = rrho3*(1.0d0 - 3.0d0*zeta) 463c d3z(3) = -rrho3*(1.0d0 + 3.0d0*zeta) 464c d3z(4) = -3.0d0*rrho3*(1.0d0 + zeta) 465#endif 466c 467c d(zeta) 468c 469 dt12 = 0.0d0 470 if (ONE+zeta.gt.1.0d-10) then 471 dt12 = dt12 + (0.5d0*(ONE+zeta))**FIVE3 472 endif 473 if (ONE-zeta.gt.1.0d-10) then 474 dt12 = dt12 + (0.5d0*(ONE-zeta))**FIVE3 475 endif 476c d1dt12 = FIVE3*0.5d0*( 477c & ((ONE+zeta)*.5d0)**TWO3 - ((ONE-zeta)*.5d0)**TWO3 ) 478 d = 2.d0**ONE3*sqrt(dt12) 479 dm1 = 1.d0/d 480c adp = 0.5d0*d/dt12*d1dt12 481c d1d(1) = adp*d1z(1) 482c d1d(2) = adp*d1z(2) 483#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 484c if ((1.d0-zeta).lt.tol_rho) then 485c d2dt12 = FIVE3*TWO3*0.25d0*(((ONE+zeta)*.5d0)**(-ONE3)) 486c else if ((1.d0+zeta).lt.tol_rho) then 487c d2dt12 = FIVE3*TWO3*0.25d0*(((ONE-zeta)*.5d0)**(-ONE3)) 488c else 489c d2dt12 = FIVE3*TWO3*0.25d0*( 490c & ((ONE+zeta)*.5d0)**(-ONE3) + ((ONE-zeta)*.5d0)**(-ONE3) ) 491c end if 492c 493c dpp =-0.5d0*adp/dt12*d1dt12 494c & + 2.d0**(-TWO3)*d2dt12/dsqrt(dt12) 495c d2d(1) = dpp*d1z(1)*d1z(1) + adp*d2z(1) 496c d2d(2) = dpp*d1z(1)*d1z(2) + adp*d2z(2) 497c d2d(3) = dpp*d1z(2)*d1z(2) + adp*d2z(3) 498#endif 499#if defined(THIRD_DERIV) 500c call errquit("nwxc_c_perdew86: no 3rd derivatives",0,0) 501#endif 502c 503c C(n) 504c 505 anum = fff+alpha*rs+beta*rs2 506 aden = 1.d0+pgamma*rs+delta*rs2+beta10*rs3 507 Cn = zzz + anum/aden 508c d1anum = alpha + 2d0*beta*rs 509c d1aden = pgamma + 2d0*delta*rs + 3d0*beta10*rs2 510#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 511c d2anum = 2d0*beta 512c d2aden = 2d0*delta + 6d0*beta10*rs 513#endif 514c First compute rs derivative 515c d1Cn = d1anum/aden - anum*d1aden/aden**2 516#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 517c d2Cn = d2anum/aden - (2d0*d1anum*d1aden+anum*d2aden)/aden**2 518c & + 2d0*anum*d1aden**2/aden**3 519#endif 520c Convert to rho derivative 521#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 522c d2Cn = d2Cn*d1rs*d1rs + d1Cn*d2rs 523#endif 524c d1Cn = d1Cn*d1rs 525c 526c phi(n,gradn) 527c 528 expfac = 0.d0 529 phi = (pfff*ftilde)*(CINF/Cn)*gam12/rho76 530 if (phi.lt.EXPTOL) expfac = exp(-phi) 531c dlnphi = -(d1Cn/Cn + SEV6/rhoval) 532c d1phi(1) = phi*dlnphi 533#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 534c d2phi(1) = d1phi(1)*dlnphi 535c & + phi*((d1Cn/Cn)**2 - d2Cn/Cn + SEV6/rhoval**2) 536#endif 537c 538c functional 539c 540 func = expfac*Cn*gamma/rho43*dm1 541c t = d1Cn/Cn - (d1phi(1) + FOUR3/rhoval) 542c dlnfrho(1) = t - dm1*d1d(1) 543c dlnfrho(2) = t - dm1*d1d(2) 544c d1f(1) = dlnfrho(1)*func 545c d1f(2) = dlnfrho(2)*func 546c Amat(n,D1_RA) = Amat(n,D1_RA) + d1f(1)*wght 547c Amat(n,D1_RB) = Amat(n,D1_RB) + d1f(2)*wght 548c if (gam12.gt.TOLL)then 549c d1phi(2) = phi / (2d0*gamma) 550c dlnfgam = 1d0/gamma - d1phi(2) 551c d1f(3) = func*dlnfgam 552c Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + d1f(3)*wght 553c Cmat(n,D1_GAB) = Cmat(n,D1_GAB) + d1f(3)*2D0*wght 554c Cmat(n,D1_GBB) = Cmat(n,D1_GBB) + d1f(3)*wght 555#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 556c d2phi(2) = d1phi(2)*dlnphi 557c d2phi(3) =-d1phi(2)/(2d0*gamma) 558c 559c t = -d2phi(1) - (d1Cn/Cn)**2 + d2Cn/Cn + FOUR3/rhoval**2 560c Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) 561c & + (d1f(1)*dlnfrho(1) 562c & + func*(d1d(1)*d1d(1)*dm1**2-d2d(1)*dm1+t))*wght 563c Amat2(n,D2_RA_RB) = Amat2(n,D2_RA_RB) 564c & + (d1f(1)*dlnfrho(2) 565c & + func*(d1d(1)*d1d(2)*dm1**2-d2d(2)*dm1+t))*wght 566c Amat2(n,D2_RB_RB) = Amat2(n,D2_RB_RB) 567c & + (d1f(2)*dlnfrho(2) 568c & + func*(d1d(2)*d1d(2)*dm1**2-d2d(3)*dm1+t))*wght 569c 570c t = (d1f(1)*dlnfgam - func*d2phi(2))*wght 571c Cmat2(n,D2_RA_GAA) = Cmat2(n,D2_RA_GAA) + t 572c Cmat2(n,D2_RA_GAB) = Cmat2(n,D2_RA_GAB) + t*2d0 573c Cmat2(n,D2_RA_GBB) = Cmat2(n,D2_RA_GBB) + t 574c t = (d1f(2)*dlnfgam - func*d2phi(2))*wght 575c Cmat2(n,D2_RB_GAA) = Cmat2(n,D2_RB_GAA) + t 576c Cmat2(n,D2_RB_GAB) = Cmat2(n,D2_RB_GAB) + t*2d0 577c Cmat2(n,D2_RB_GBB) = Cmat2(n,D2_RB_GBB) + t 578c 579c t = (d1f(3)*dlnfgam - func*(1d0/gamma**2+d2phi(3)))*wght 580c Cmat2(n,D2_GAA_GAA) = Cmat2(n,D2_GAA_GAA) + t 581c Cmat2(n,D2_GAA_GBB) = Cmat2(n,D2_GAA_GBB) + t 582c Cmat2(n,D2_GBB_GBB) = Cmat2(n,D2_GBB_GBB) + t 583c Cmat2(n,D2_GAA_GAB) = Cmat2(n,D2_GAA_GAB) + t*2d0 584c Cmat2(n,D2_GAB_GBB) = Cmat2(n,D2_GAB_GBB) + t*2d0 585c Cmat2(n,D2_GAB_GAB) = Cmat2(n,D2_GAB_GAB) + t*4d0 586#endif 587c endif 588 ffunc(n)=ffunc(n)+func*wght 589 20 continue 590 endif 591 return 592 end 593#ifndef NWAD_PRINT 594#define NWAD_PRINT 595c 596c Compile source again for Maxima 597c 598#include "nwxc_c_perdew86.F" 599#endif 600#ifndef SECOND_DERIV 601#define SECOND_DERIV 602c 603c Compile source again for the 2nd derivative case 604c 605#include "nwxc_c_perdew86.F" 606#endif 607#ifndef THIRD_DERIV 608#define THIRD_DERIV 609c 610c Compile source again for the 3rd derivative case 611c 612#include "nwxc_c_perdew86.F" 613#endif 614#undef NWAD_PRINT 615C> 616C> @} 617