1#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 2 Subroutine xc_perdew86(tol_rho, fac, lfac, nlfac, rho, delrho, 3 & Amat, Cmat, nq, ipol, Ec, qwght, 4 , ldew, ffunc) 5#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV) 6 Subroutine xc_perdew86_d2(tol_rho, fac, lfac, nlfac, rho, delrho, 7 & Amat, Amat2, Cmat, Cmat2, nq, ipol, Ec, 8 , qwght, ldew, ffunc) 9#else 10 Subroutine xc_perdew86_d3(tol_rho, fac, lfac, nlfac, rho, delrho, 11 1 Amat, Amat2, Amat3, Cmat, Cmat2, Cmat3, nq, ipol, Ec, 12 2 qwght, ldew, ffunc) 13#endif 14c 15c$Id$ 16c 17 implicit none 18c 19#include "dft2drv.fh" 20#include "dft3drv.fh" 21c 22 double precision tol_rho, fac ! [input] 23 integer nq, ipol ! [input] 24 double precision Ec ! [input/output] 25 logical lfac, nlfac, ldew 26 double precision ffunc(*) ! value of the functional [output] 27c 28c Charge Density 29c 30 double precision rho(nq,ipol*(ipol+1)/2) 31c 32c Charge Density Gradient 33c 34 double precision delrho(nq,3,ipol) 35c 36c Quadrature Weights 37c 38 double precision qwght(nq) 39c 40c Sampling Matrices for the XC Potential & Energy 41c 42 double precision Amat(nq,ipol), Cmat(nq,*) 43#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 44 double precision Amat2(nq,NCOL_AMAT2), Cmat2(nq,NCOL_CMAT2) 45#endif 46#ifdef THIRD_DERIV 47 double precision Amat3(nq,NCOL_AMAT3), Cmat3(nq,NCOL_CMAT3) 48#endif 49 double precision TOLL, EXPTOL, alpha, beta, pgamma, delta, 50 & beta10, ftilde, zzz, fff, pfff, CINF, ONE, 51 & ONE3, THREE, FOUR3, SEV6, FIVE3, 52 & TWO3, FIVE6, pi 53 double precision SEVEN3, EIGHT3 54 Parameter (TOLL = 1.D-40, EXPTOL = 80.d0) 55 Parameter (alpha = 0.023266D0, beta = 7.389D-6, 56 & pgamma = 8.723d0, delta = 0.472d0, beta10 = 10000.d0*beta) 57 parameter (ftilde = 0.11d0, zzz = 0.001667d0, fff = 0.002568d0) 58 parameter(pfff = 1.745d0, CINF = zzz+fff) 59 Parameter (ONE = 1.D0, ONE3 = 1.d0/3.d0, THREE = 3.d0) 60 Parameter (FOUR3 = 4.D0/3.D0, SEV6 = 7.d0/6.d0) 61 parameter (FIVE3 = 5.d0/3.d0, TWO3 = 2.d0/3.d0, FIVE6 = 5.d0/6.d0) 62 parameter (pi = 3.1415926535897932385d0) 63 parameter (SEVEN3 = 7.0d0/3.0d0, EIGHT3 = 8.0d0/3.0d0) 64c 65c Mlynarski Salahub PRB 43, 1399 (1991) 66c 67 integer n 68 double precision rsfact, rs, rs2, rs3 69 double precision rhoval, rho13, rho43, rho76, arho 70 double precision d1rs 71#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 72 double precision d2rs 73#endif 74#ifdef THIRD_DERIV 75 double precision d3rs 76#endif 77 double precision gamma, gam12 78 double precision anum, aden, d1anum, d1aden, Cn, d1Cn, 79 & expfac, phi, d1phi(2), dlnphi, func, d1f(3), 80 & dlnfrho(2), dlnfgam 81 double precision zeta, d1z(2), d, dm1, adp, d1d(2), t, 82 & dt12, d1dt12 83 double precision aden2 84#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 85 double precision d2anum, d2aden, rrho2, d2z(3), dpp, d2d(3), 86 & d2phi(3), d2dt12, d2Cn 87 double precision aden3 88 double precision arho2 89 double precision d2lnphi 90 double precision d2f(3) 91 double precision d2lnfrho(3), d2lnfrg(2), d2lnfgam 92#endif 93#ifdef THIRD_DERIV 94 double precision d3lnphi 95 double precision d3anum, d3aden, d3Cn, d3phi(4) 96 double precision d3lnfrho(4), d3lnfgam 97 double precision d3f(3) 98 double precision aden4 99 double precision arho3 100#endif 101c 102 rsfact = (0.75d0/pi)**ONE3 103c 104 if (ipol.eq.1 )then 105c 106c ======> SPIN-RESTRICTED <====== 107c 108 do 10 n = 1, nq 109 rhoval = rho(n,1) 110 if (rhoval.lt.tol_rho) goto 10 111 arho=1.d0/rhoval 112 rho13 = abs(rhoval)**ONE3 113 rho43 = rhoval*rho13 114 rho76 = abs(rhoval)**SEV6 115 rs = rsfact/rho13 116 rs2 = rs*rs 117 rs3 = rs2*rs 118 d1rs = -ONE3*rs*arho 119#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 120 d2rs = -FOUR3*d1rs*arho 121#endif 122#ifdef THIRD_DERIV 123 d3rs = -SEVEN3*d2rs*arho 124#endif 125 gamma = delrho(n,1,1)*delrho(n,1,1) + 126 & delrho(n,2,1)*delrho(n,2,1) + 127 & delrho(n,3,1)*delrho(n,3,1) 128 gam12 = sqrt(abs(gamma)) 129c 130c C(n) 131c 132 anum = fff+alpha*rs+beta*rs2 133 aden = 1.d0+pgamma*rs+delta*rs2+beta10*rs3 134 Cn = zzz + anum/aden 135 d1anum = alpha + 2d0*beta*rs 136 d1aden = pgamma + 2d0*delta*rs + 3d0*beta10*rs2 137#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 138 d2anum = 2d0*beta 139 d2aden = 2d0*delta + 6d0*beta10*rs 140#endif 141#ifdef THIRD_DERIV 142 d3anum = 0.0d0 143 d3aden = 6.0d0*beta10 144#endif 145c First compute rs derivative 146 aden2 = aden*aden 147c 148 d1Cn = d1anum/aden - anum*d1aden/aden2 149#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 150 aden3 = aden2*aden 151c 152c d2Cn = d2anum/aden - (2d0*d1anum*d1aden+anum*d2aden)/aden**2 153c & + 2d0*anum*d1aden**2/aden**3 154 d2Cn = d2anum/aden - (2d0*d1anum*d1aden+anum*d2aden)/aden2 155 & + 2d0*anum*d1aden**2/aden3 156#endif 157#ifdef THIRD_DERIV 158 aden4 = aden3*aden 159c 160 d3Cn = -( 3.0d0*d2anum*d1aden + 3.0d0*d1anum*d2aden 161 1 + anum*d3aden )/aden2 162 2 + 6.0d0*( d1anum*d1aden**2 163 3 + anum*d2aden*d1aden )/aden3 164 4 - 6.0d0*anum*d1aden**3/aden4 165#endif 166c Convert to rho derivative 167#ifdef THIRD_DERIV 168 d3Cn = d3Cn*d1rs*d1rs*d1rs 169 1 + 3.0d0*d2Cn*d2rs*d1rs 170 2 + d1Cn*d3rs 171#endif 172c 173#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 174 d2Cn = d2Cn*d1rs*d1rs + d1Cn*d2rs 175#endif 176 d1Cn = d1Cn*d1rs 177c 178c phi(n,gradn) 179c 180 expfac = 0.d0 181 phi = (pfff*ftilde)*(CINF/Cn)*gam12/rho76 182 if (phi.lt.EXPTOL) expfac = exp(-phi) 183 dlnphi = -(d1Cn/Cn + SEV6/rhoval) 184 d1phi(1) = phi*dlnphi 185#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 186 arho2 = arho*arho 187c 188 d2lnphi = (d1Cn/Cn)**2 - d2Cn/Cn + SEV6*arho2 189c 190c d2phi(1) = d1phi(1)*dlnphi 191c & + phi*((d1Cn/Cn)**2 - d2Cn/Cn + SEV6/rhoval**2) 192 d2phi(1) = d1phi(1)*dlnphi + phi*d2lnphi 193#endif 194c 195#ifdef THIRD_DERIV 196 arho3 = arho2*arho 197c 198 d3lnphi = -2.0d0*(d1Cn/Cn)**3 199 1 + 3.0d0*(d2Cn/Cn)*(d1Cn/Cn) 200 2 - d3Cn/Cn 201 3 - SEVEN3*arho3 202 d3phi(1) = d2phi(1)*dlnphi 203 1 + 2.0d0*d1phi(1)*d2lnphi 204 2 + phi*d3lnphi 205#endif 206c 207c functional 208c 209 func = expfac*Cn*gamma/rho43 210 dlnfrho(1) = d1Cn/Cn - (d1phi(1) + FOUR3/rhoval) 211 d1f(1) = dlnfrho(1)*func 212 Amat(n,1) = Amat(n,1) + d1f(1)*fac 213 if (gam12.gt.TOLL)then 214 d1phi(2) = phi / (2d0*gamma) 215 dlnfgam = 1d0/gamma - d1phi(2) 216 d1f(3) = func*dlnfgam 217 Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + d1f(3)*fac 218 Cmat(n,D1_GAB) = Cmat(n,D1_GAB) + d1f(3)*2D0*fac 219#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 220 d2phi(2) = d1phi(2)*dlnphi 221 d2phi(3) =-d1phi(2)/(2d0*gamma) 222c rr terms 223c!!! Which of the following are actually needed for restricted? 224c!!! Should treat derivatives of d as zero? d is a constant? 225c Daniel (11-19-12): d is a constant (it equals 1) for a restricted 226c calculation, since there is no spin-polarization. Thus, the 227c derivatives are zero. 228 d2lnfrho(1) = -d2phi(1) - (d1Cn/Cn)**2 + d2Cn/Cn 229 1 + FOUR3*arho2 230c 231 d2f(1) = d1f(1)*dlnfrho(1) 232 1 + func*d2lnfrho(1) 233c 234 t = d2f(1)*fac 235c 236 Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) + t 237 Amat2(n,D2_RA_RB) = Amat2(n,D2_RA_RB) + t 238c OLD CODE 239c t = -d2phi(1) - (d1Cn/Cn)**2 + d2Cn/Cn + FOUR3/rhoval**2 240c Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) 241c & + (d1f(1)*dlnfrho(1) 242c & + func*t)*fac 243c Amat2(n,D2_RA_RB) = Amat2(n,D2_RA_RB) 244c & + (d1f(1)*dlnfrho(1) 245c & + func*t)*fac 246c OLD CODE 247#if 0 248 Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) 249 & + (d1f(1)*dlnfrho(1) 250 & + func*(d1d(1)*d1d(1)*dm1**2-d2d(1)*dm1+t))*fac 251 Amat2(n,D2_RA_RB) = Amat2(n,D2_RA_RB) 252 & + (d1f(1)*dlnfrho(2) 253 & + func*(d1d(1)*d1d(2)*dm1**2-d2d(2)*dm1+t))*fac 254#endif 255c rg terms 256 d2lnfrg(1) = -d2phi(2) 257 d2f(2) = (d1f(1)*dlnfgam + func*d2lnfrg(1)) 258 t = d2f(2)*fac 259c 260 Cmat2(n,D2_RA_GAA) = Cmat2(n,D2_RA_GAA) + t 261 Cmat2(n,D2_RA_GAB) = Cmat2(n,D2_RA_GAB) + t*2d0 262 Cmat2(n,D2_RA_GBB) = Cmat2(n,D2_RA_GBB) + t 263c OLD CODE 264c t = (d1f(1)*dlnfgam - func*d2phi(2))*fac 265c Cmat2(n,D2_RA_GAA) = Cmat2(n,D2_RA_GAA) + t 266c Cmat2(n,D2_RA_GAB) = Cmat2(n,D2_RA_GAB) + t*2d0 267c Cmat2(n,D2_RA_GBB) = Cmat2(n,D2_RA_GBB) + t 268c OLD CODE 269c gg terms 270 d2lnfgam = -1.0d0/gamma**2 - d2phi(3) 271 d2f(3) = d1f(3)*dlnfgam + func*d2lnfgam 272 t = d2f(3)*fac 273c 274 Cmat2(n,D2_GAA_GAA) = Cmat2(n,D2_GAA_GAA) + t 275 Cmat2(n,D2_GAA_GBB) = Cmat2(n,D2_GAA_GBB) + t 276 Cmat2(n,D2_GAA_GAB) = Cmat2(n,D2_GAA_GAB) + t*2d0 277 Cmat2(n,D2_GAB_GAB) = Cmat2(n,D2_GAB_GAB) + t*4d0 278c OLD CODE 279c t = (d1f(3)*dlnfgam - func*(1d0/gamma**2+d2phi(3)))*fac 280c Cmat2(n,D2_GAA_GAA) = Cmat2(n,D2_GAA_GAA) + t 281c Cmat2(n,D2_GAA_GBB) = Cmat2(n,D2_GAA_GBB) + t 282c Cmat2(n,D2_GAA_GAB) = Cmat2(n,D2_GAA_GAB) + t*2d0 283c Cmat2(n,D2_GAB_GAB) = Cmat2(n,D2_GAB_GAB) + t*4d0 284c OLD CODE 285#endif 286c 287#ifdef THIRD_DERIV 288c rrr terms 289 d3lnfrho(1) = -d3phi(1) 290 1 + 2.0d0*(d1Cn/Cn)**3 291 2 - 3.0d0*(d2Cn/Cn)*(d1Cn/Cn) 292 3 + d3Cn/Cn 293 4 - EIGHT3*arho3 294c 295 d3f(1) = d2f(1)*dlnfrho(1) 296 1 + 2.0d0*d1f(1)*d2lnfrho(1) 297 2 + func*d3lnfrho(1) 298c 299 t = d3f(1)*fac 300c 301 Amat3(n,D3_RA_RA_RA) = Amat3(n,D3_RA_RA_RA) + t 302 Amat3(n,D3_RA_RA_RB) = Amat3(n,D3_RA_RA_RB) + t 303 Amat3(n,D3_RA_RB_RB) = Amat3(n,D3_RA_RB_RB) + t 304c rrg terms 305 d3phi(2) = d2phi(2)*dlnphi + d1phi(2)*d2lnphi 306c 307 t = ( d2f(2)*dlnfrho(1) 308 1 - d1f(1)*d2phi(2) 309 2 + d1f(3)*d2lnfrho(1) 310 3 - func*d3phi(2) )*fac 311c 312 Cmat3(n,D3_RA_RA_GAA) = Cmat3(n,D3_RA_RA_GAA) + t 313 Cmat3(n,D3_RA_RA_GAB) = Cmat3(n,D3_RA_RA_GAB) + t*2.0d0 314 Cmat3(n,D3_RA_RA_GBB) = Cmat3(n,D3_RA_RA_GBB) + t 315 Cmat3(n,D3_RA_RB_GAA) = Cmat3(n,D3_RA_RB_GAA) + t 316 Cmat3(n,D3_RA_RB_GAB) = Cmat3(n,D3_RA_RB_GAB) + t*2.0d0 317 Cmat3(n,D3_RA_RB_GBB) = Cmat3(n,D3_RA_RB_GBB) + t 318c rgg terms 319 d3phi(3) = -d2phi(3)*dlnphi 320c 321 t = ( d2f(2)*dlnfgam 322 1 + d1f(1)*d2lnfgam 323 2 + d1f(3)*d2lnfrg(1) 324 3 + func*d3phi(3) )*fac 325c 326 Cmat3(n,D3_RA_GAA_GAA) = Cmat3(n,D3_RA_GAA_GAA) + t 327 Cmat3(n,D3_RA_GAA_GAB) = Cmat3(n,D3_RA_GAA_GAB) + t*2.0d0 328 Cmat3(n,D3_RA_GAA_GBB) = Cmat3(n,D3_RA_GAA_GBB) + t 329 Cmat3(n,D3_RA_GAB_GAB) = Cmat3(n,D3_RA_GAB_GAB) + t*4.0d0 330 Cmat3(n,D3_RA_GAB_GBB) = Cmat3(n,D3_RA_GAB_GBB) + t*2.0d0 331 Cmat3(n,D3_RA_GBB_GBB) = Cmat3(n,D3_RA_GBB_GBB) + t 332c ggg terms 333 d3phi(4) = -3.0d0*d2phi(3)/(2.0d0*gamma) 334 d3lnfgam = 2.0d0/gamma**3 - d3phi(4) 335c 336 t = ( d2f(3)*dlnfgam 337 1 + 2.0d0*d1f(3)*d2lnfgam 338 2 + func*d3lnfgam )*fac 339c 340 Cmat3(n,D3_GAA_GAA_GAA) = Cmat3(n,D3_GAA_GAA_GAA) + t 341 Cmat3(n,D3_GAA_GAA_GAB) = Cmat3(n,D3_GAA_GAA_GAB) 342 1 + t*2.0d0 343 Cmat3(n,D3_GAA_GAA_GBB) = Cmat3(n,D3_GAA_GAA_GBB) + t 344 Cmat3(n,D3_GAA_GAB_GAB) = Cmat3(n,D3_GAA_GAB_GAB) 345 1 + t*4.0d0 346 Cmat3(n,D3_GAA_GAB_GBB) = Cmat3(n,D3_GAA_GAB_GBB) 347 1 + t*2.0d0 348 Cmat3(n,D3_GAA_GBB_GBB) = Cmat3(n,D3_GAA_GBB_GBB) + t 349 Cmat3(n,D3_GAB_GAB_GAB) = Cmat3(n,D3_GAB_GAB_GAB) 350 1 + t*8.0d0 351#endif 352 endif 353 Ec = Ec + func*qwght(n)*fac 354 if (ldew) ffunc(n)=ffunc(n)+func*fac 355 10 continue 356 else 357c 358c ======> SPIN-UNRESTRICTED <====== 359c 360 do 20 n = 1, nq 361 rhoval = rho(n,1) 362 if (rhoval.lt.tol_rho) goto 20 363 arho=1.d0/rhoval 364 rho13 = abs(rhoval)**ONE3 365 rho43 = rhoval*rho13 366 rho76 = abs(rhoval)**SEV6 367 rs = rsfact/rho13 368 rs2 = rs*rs 369 rs3 = rs2*rs 370 d1rs = -ONE3*rs*arho 371#ifdef SECOND_DERIV 372 d2rs = -FOUR3*d1rs*arho 373#endif 374 gamma = delrho(n,1,1)*delrho(n,1,1) + 375 & delrho(n,2,1)*delrho(n,2,1) + 376 & delrho(n,3,1)*delrho(n,3,1) + 377 & delrho(n,1,2)*delrho(n,1,2) + 378 & delrho(n,2,2)*delrho(n,2,2) + 379 & delrho(n,3,2)*delrho(n,3,2) + 380 & 2.d0*(delrho(n,1,1)*delrho(n,1,2) + 381 & delrho(n,2,1)*delrho(n,2,2) + 382 & delrho(n,3,1)*delrho(n,3,2)) 383 gam12 = sqrt(abs(gamma)) 384 zeta = (rho(n,2) - rho(n,3))*arho 385 if(zeta.lt.-1d0) zeta=-1d0 386 if(zeta.gt.1d0) zeta=1d0 387 d1z(1) = (1.d0 - zeta)*arho 388 d1z(2) = -(1.d0 + zeta)*arho 389#ifdef SECOND_DERIV 390 rrho2 = 2.d0*arho*arho 391c 1 = aa, 2 = ab, 3 = bb 392 d2z(1) =-rrho2*(1.d0-zeta) 393 d2z(2) = rrho2*zeta 394 d2z(3) = rrho2*(1.d0+zeta) 395#endif 396c 397c d(zeta) 398c 399 dt12 = ((ONE+zeta)*.5d0)**FIVE3 + ((ONE-zeta)*.5d0)**FIVE3 400 d1dt12 = FIVE3*0.5d0*( 401 & ((ONE+zeta)*.5d0)**TWO3 - ((ONE-zeta)*.5d0)**TWO3 ) 402 d = 2.d0**ONE3*dsqrt(dt12) 403 dm1 = 1.d0/d 404 adp = 0.5d0*d/dt12*d1dt12 405 d1d(1) = adp*d1z(1) 406 d1d(2) = adp*d1z(2) 407#ifdef SECOND_DERIV 408 if ((1.d0-zeta).lt.tol_rho) then 409 d2dt12 = FIVE3*TWO3*0.25d0*(((ONE+zeta)*.5d0)**(-ONE3)) 410 else if ((1.d0+zeta).lt.tol_rho) then 411 d2dt12 = FIVE3*TWO3*0.25d0*(((ONE-zeta)*.5d0)**(-ONE3)) 412 else 413 d2dt12 = FIVE3*TWO3*0.25d0*( 414 & ((ONE+zeta)*.5d0)**(-ONE3) + ((ONE-zeta)*.5d0)**(-ONE3) ) 415 end if 416c 417 dpp =-0.5d0*adp/dt12*d1dt12 418 & + 2.d0**(-TWO3)*d2dt12/dsqrt(dt12) 419 d2d(1) = dpp*d1z(1)*d1z(1) + adp*d2z(1) 420 d2d(2) = dpp*d1z(1)*d1z(2) + adp*d2z(2) 421 d2d(3) = dpp*d1z(2)*d1z(2) + adp*d2z(3) 422#endif 423c 424c C(n) 425c 426 anum = fff+alpha*rs+beta*rs2 427 aden = 1.d0+pgamma*rs+delta*rs2+beta10*rs3 428 Cn = zzz + anum/aden 429 d1anum = alpha + 2d0*beta*rs 430 d1aden = pgamma + 2d0*delta*rs + 3d0*beta10*rs2 431#ifdef SECOND_DERIV 432 d2anum = 2d0*beta 433 d2aden = 2d0*delta + 6d0*beta10*rs 434#endif 435c First compute rs derivative 436 d1Cn = d1anum/aden - anum*d1aden/aden**2 437#ifdef SECOND_DERIV 438 d2Cn = d2anum/aden - (2d0*d1anum*d1aden+anum*d2aden)/aden**2 439 & + 2d0*anum*d1aden**2/aden**3 440#endif 441c Convert to rho derivative 442#ifdef SECOND_DERIV 443 d2Cn = d2Cn*d1rs*d1rs + d1Cn*d2rs 444#endif 445 d1Cn = d1Cn*d1rs 446c 447c phi(n,gradn) 448c 449 expfac = 0.d0 450 phi = (pfff*ftilde)*(CINF/Cn)*gam12/rho76 451 if (phi.lt.EXPTOL) expfac = exp(-phi) 452 dlnphi = -(d1Cn/Cn + SEV6/rhoval) 453 d1phi(1) = phi*dlnphi 454#ifdef SECOND_DERIV 455 d2phi(1) = d1phi(1)*dlnphi 456 & + phi*((d1Cn/Cn)**2 - d2Cn/Cn + SEV6/rhoval**2) 457#endif 458c 459c functional 460c 461 func = expfac*Cn*gamma/rho43*dm1 462 t = d1Cn/Cn - (d1phi(1) + FOUR3/rhoval) 463 dlnfrho(1) = t - dm1*d1d(1) 464 dlnfrho(2) = t - dm1*d1d(2) 465 d1f(1) = dlnfrho(1)*func 466 d1f(2) = dlnfrho(2)*func 467 Amat(n,1) = Amat(n,1) + d1f(1)*fac 468 Amat(n,2) = Amat(n,2) + d1f(2)*fac 469 if (gam12.gt.TOLL)then 470 d1phi(2) = phi / (2d0*gamma) 471 dlnfgam = 1d0/gamma - d1phi(2) 472 d1f(3) = func*dlnfgam 473 Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + d1f(3)*fac 474 Cmat(n,D1_GAB) = Cmat(n,D1_GAB) + d1f(3)*2D0*fac 475 Cmat(n,D1_GBB) = Cmat(n,D1_GBB) + d1f(3)*fac 476#ifdef SECOND_DERIV 477 d2phi(2) = d1phi(2)*dlnphi 478 d2phi(3) =-d1phi(2)/(2d0*gamma) 479c 480 t = -d2phi(1) - (d1Cn/Cn)**2 + d2Cn/Cn + FOUR3/rhoval**2 481 Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) 482 & + (d1f(1)*dlnfrho(1) 483 & + func*(d1d(1)*d1d(1)*dm1**2-d2d(1)*dm1+t))*fac 484 Amat2(n,D2_RA_RB) = Amat2(n,D2_RA_RB) 485 & + (d1f(1)*dlnfrho(2) 486 & + func*(d1d(1)*d1d(2)*dm1**2-d2d(2)*dm1+t))*fac 487 Amat2(n,D2_RB_RB) = Amat2(n,D2_RB_RB) 488 & + (d1f(2)*dlnfrho(2) 489 & + func*(d1d(2)*d1d(2)*dm1**2-d2d(3)*dm1+t))*fac 490c 491 t = (d1f(1)*dlnfgam - func*d2phi(2))*fac 492 Cmat2(n,D2_RA_GAA) = Cmat2(n,D2_RA_GAA) + t 493 Cmat2(n,D2_RA_GAB) = Cmat2(n,D2_RA_GAB) + t*2d0 494 Cmat2(n,D2_RA_GBB) = Cmat2(n,D2_RA_GBB) + t 495 t = (d1f(2)*dlnfgam - func*d2phi(2))*fac 496 Cmat2(n,D2_RB_GAA) = Cmat2(n,D2_RB_GAA) + t 497 Cmat2(n,D2_RB_GAB) = Cmat2(n,D2_RB_GAB) + t*2d0 498 Cmat2(n,D2_RB_GBB) = Cmat2(n,D2_RB_GBB) + t 499c 500 t = (d1f(3)*dlnfgam - func*(1d0/gamma**2+d2phi(3)))*fac 501 Cmat2(n,D2_GAA_GAA) = Cmat2(n,D2_GAA_GAA) + t 502 Cmat2(n,D2_GAA_GBB) = Cmat2(n,D2_GAA_GBB) + t 503 Cmat2(n,D2_GBB_GBB) = Cmat2(n,D2_GBB_GBB) + t 504 Cmat2(n,D2_GAA_GAB) = Cmat2(n,D2_GAA_GAB) + t*2d0 505 Cmat2(n,D2_GAB_GBB) = Cmat2(n,D2_GAB_GBB) + t*2d0 506 Cmat2(n,D2_GAB_GAB) = Cmat2(n,D2_GAB_GAB) + t*4d0 507#endif 508 endif 509 Ec = Ec + func*qwght(n)*fac 510 if (ldew) ffunc(n)=ffunc(n)+func*fac 511 20 continue 512 endif 513 return 514 end 515 516#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 517 Subroutine xc_p81(tol_rho, fac, lfac, nlfac, rho, Amat, nq, ipol, 518 & Ec, qwght, ldew, func) 519#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV) 520#include "dft2drv.fh" 521 Subroutine xc_p81_d2(tol_rho, fac, lfac, nlfac, rho, Amat, Amat2, 522 & nq, ipol, Ec, qwght, ldew, func) 523#else 524#include "dft3drv.fh" 525 Subroutine xc_p81_d3(tol_rho, fac, lfac, nlfac, rho, Amat, Amat2, 526 & Amat3, nq, ipol, Ec, qwght, ldew, func) 527#endif 528c Daniel (4-2-13): Third derivatives aren't implemented for Perdew 81 529c LDA yet. The preprocessor stuff is here to keep the compiler from 530c complaining. We need this functional for Perdew 86 to work. 531c 532c Ceperley Alder LDA from Perdew Zunger PRB 23, 5048 (1981) 533c 534 implicit none 535c 536 integer nq, ipol 537 logical lfac, nlfac, ldew 538 double precision func(*) ! value of the functional [output] 539 double precision Ec, fac 540c 541c Charge Density 542c 543 double precision rho(nq,ipol*(ipol+1)/2) 544c 545c Quadrature Weights 546c 547 double precision qwght(nq) 548c 549c Sampling Matrices for the XC Potential & Energy 550c 551 double precision Amat(nq,ipol) 552#ifdef SECOND_DERIV 553c double precision Amat2(nq,*) 554 double precision Amat2(nq,NCOL_AMAT2) 555#endif 556c 557#ifdef THIRD_DERIV 558 double precision Amat3(nq,NCOL_AMAT3) 559#endif 560 double precision A(2), B(2), C(2), D(2), G(2), B1(2), B2(2), 561 & pi, tol_rho, ONE3, FOUR3, TWO3 562 double precision FIVE3, SEVEN3 563 save A, B, C, D, G, B1, B2 564 parameter (pi = 3.1415926535897932385d0) 565 Parameter (ONE3 = 1.d0/3.d0, FOUR3 = 4.D0/3.D0) 566 Parameter (TWO3 = 2.d0/3.d0) 567 Parameter (FIVE3 = 5.0d0/3.0d0, SEVEN3 = 7.0d0/3.0d0) 568 integer n, i 569 double precision rhoval, rs, alnrs, d1rs, e(2), d1e(2), rden(2), 570 & d1den(2), d1zeta(2), d1ersz(2), d1edrho(2), eps, 571 & sqrtrs, fz, d1fz, zeta 572#ifdef SECOND_DERIV 573 double precision d2rs, d2e(2), d2den(2), d2zeta(3), d2ersz(3), 574 & d2edrho(3), d2fzeta, d2fz, rrho2 575#endif 576#ifdef THIRD_DERIV 577 double precision d3rs, d3fz, rrho3, d3zeta(4), d3den(2), d3e(2), 578 1 d3ersz(4), d3edrho(4) 579#endif 580 double precision x, fzeta, d1fzeta, rsfact 581 fzeta(x) = ((1.d0+x)**FOUR3 + 582 & (1.d0-x)**FOUR3 - 2.d0) / (2.d0**FOUR3-2.d0) 583 d1fzeta(x) = FOUR3*((1.d0+x)**ONE3 - 584 & (1.d0-x)**ONE3) / (2.d0**FOUR3-2.d0) 585#ifdef SECOND_DERIV 586 d2fzeta(x) = ONE3*FOUR3*((1.d0+x)**(-TWO3) + 587 & (1.d0-x)**(-TWO3)) / (2.d0**FOUR3-2.d0) 588#endif 589 data A / 0.0311d0, 0.01555d0 / 590 data B / -0.048d0, -0.0269d0 / 591 data C / 0.0020d0, 0.0007d0 / 592 data D / -0.0116d0, -0.0048d0 / 593 data G / -.1423d0, -.0843d0 / 594 data B1 / 1.0529d0, 1.3981d0 / 595 data B2 / 0.3334d0, 0.2611d0 / 596c 597 rsfact = (0.75d0/pi)**ONE3 598c 599c ======> BOTH SPIN-RESTRICTED AND UNRESTRICTED <====== 600c 601 do n = 1, nq 602 if (rho(n,1).gt.tol_rho)then 603 rhoval = rho(n,1) 604 if (ipol.eq.1) then 605 zeta = 0.0d0 606 d1zeta(1) = 1.d0/rhoval 607 d1zeta(2) =-1.d0/rhoval 608 fz = 0d0 609 d1fz = 0d0 610 else 611 zeta = (rho(n,2)-rho(n,3))/rhoval 612 if(zeta.lt.-1d0) zeta=-1d0 613 if(zeta.gt.1d0) zeta=1d0 614 fz = fzeta(zeta) 615 d1fz = d1fzeta(zeta) 616 d1zeta(1) = (1.d0-zeta)/rhoval 617 d1zeta(2) =-(1.d0+zeta)/rhoval 618 endif 619 rs = rsfact/abs(rhoval)**ONE3 620 d1rs = -ONE3*rs/rhoval 621#ifdef SECOND_DERIV 622 d2rs = -FOUR3*d1rs/rhoval 623 if ((1.d0-zeta).lt.tol_rho) then 624 d2fz = (1.d0+zeta)**(-TWO3) 625 else if ((1.d0+zeta).lt.tol_rho) then 626 d2fz = (1.d0-zeta)**(-TWO3) 627 else 628 d2fz = (1.d0+zeta)**(-TWO3) + (1.d0-zeta)**(-TWO3) 629 end if 630 d2fz = d2fz*ONE3*FOUR3/(2.d0**FOUR3-2.d0) 631c 632 rrho2 = 2.d0/(rhoval*rhoval) 633c 1 = aa, 2 = ab, 3 = bb 634 d2zeta(1) =-rrho2*(1.d0-zeta) 635 d2zeta(2) = rrho2*zeta 636 d2zeta(3) = rrho2*(1.d0+zeta) 637#endif 638c 639#ifdef THIRD_DERIV 640 d3rs = -SEVEN3*d2rs/rhoval 641 if ((1.d0-zeta).lt.tol_rho) then 642 d3fz = (1.d0+zeta)**(-FIVE3) 643 else if ((1.d0+zeta).lt.tol_rho) then 644 d3fz = (1.d0-zeta)**(-FIVE3) 645 else 646 d3fz = (1.d0+zeta)**(-FIVE3) + (1.d0-zeta)**(-FIVE3) 647 end if 648 d3fz = -d3fz*TWO3*ONE3*FOUR3/(2.d0**FOUR3-2.d0) 649c 650 rrho3 = 2.0d0/(rhoval*rhoval*rhoval) 651c 652c 1 = aaa, 2 = aab, 3 = abb, 4 = bbb 653 d3zeta(1) = 3.0d0*rrho3*(1.0d0 - zeta) 654 d3zeta(2) = rrho3*(1.0d0 - 3.0d0*zeta) 655 d3zeta(3) = -rrho3*(1.0d0 + 3.0d0*zeta) 656 d3zeta(4) = -3.0d0*rrho3*(1.0d0 + zeta) 657#endif 658 if (rs.lt.1.d0)then 659 alnrs = log(rs) 660 do i = 1, 2 661 e(i) = A(i)*alnrs+B(i)+C(i)*rs*alnrs+D(i)*rs 662 d1e(i) = A(i)/rs+C(i)*(alnrs+1d0)+D(i) 663#ifdef SECOND_DERIV 664 d2e(i) = (C(i)-A(i)/rs)/rs 665#endif 666#ifdef THIRD_DERIV 667 d3e(i) = 2.0d0*A(i)/(rs*rs*rs) 668 1 - C(i)/(rs*rs) 669#endif 670 enddo 671 else 672 sqrtrs = sqrt(rs) 673 do i = 1, 2 674 rden(i) = 1.d0/(1.d0+B1(i)*sqrtrs+B2(i)*rs) 675 d1den(i) = B1(i)/(2.d0*sqrtrs)+B2(i) 676 e(i) = G(i)*rden(i) 677 d1e(i) = -G(i)*d1den(i)*rden(i)**2 678#ifdef SECOND_DERIV 679 d2den(i) = -B1(i)/(4.d0*rs*sqrtrs) 680 d2e(i) = G(i)*rden(i)**2 681 & *(2.d0*d1den(i)**2*rden(i)-d2den(i)) 682#endif 683#ifdef THIRD_DERIV 684 d3den(i) = 3.0d0*B1(i)/(8.0d0*rs*rs*sqrtrs) 685 d3e(i) = G(i)*rden(i)*rden(i)* 686 1 ( 6.0d0*( d1den(i)*d2den(i)*rden(i) 687 2 - d1den(i)*d1den(i)*d1den(i)* 688 3 rden(i)*rden(i) ) 689 4 - d3den(i) ) 690#endif 691 enddo 692 endif 693 eps = e(1) + fz*(e(2)-e(1)) 694 d1ersz(1) = d1e(1) + fz*(d1e(2)-d1e(1)) 695 d1ersz(2) = d1fz*(e(2)-e(1)) 696 d1edrho(1) = d1ersz(1)*d1rs + d1ersz(2)*d1zeta(1) 697 d1edrho(2) = d1ersz(1)*d1rs + d1ersz(2)*d1zeta(2) 698 Ec = Ec + eps*qwght(n)*rhoval*fac 699 if (ldew) func(n) = func(n) + eps*rhoval*fac 700 Amat(n,1) = Amat(n,1) + (eps + rhoval*d1edrho(1))*fac 701 if (ipol.eq.2) 702 & Amat(n,2) = Amat(n,2) + (eps + rhoval*d1edrho(2))*fac 703#ifdef SECOND_DERIV 704c 1 = rsrs, 2 = rsz, 3 = zz 705 d2ersz(1) = d2e(1) + fz*(d2e(2)-d2e(1)) 706 d2ersz(2) = d1fz*(d1e(2)-d1e(1)) 707 d2ersz(3) = d2fz*(e(2)-e(1)) 708c 1 = aa, 2 = ab, 3 = bb 709 d2edrho(1) = d2ersz(1)*d1rs*d1rs 710 & + d2ersz(2)*d1rs*d1zeta(1)*2.d0 711 & + d2ersz(3)*d1zeta(1)*d1zeta(1) 712 & + d1ersz(1)*d2rs 713 & + d1ersz(2)*d2zeta(1) 714 d2edrho(2) = d2ersz(1)*d1rs*d1rs 715 & + d2ersz(2)*d1rs*(d1zeta(1)+d1zeta(2)) 716 & + d2ersz(3)*d1zeta(1)*d1zeta(2) 717 & + d1ersz(1)*d2rs 718 & + d1ersz(2)*d2zeta(2) 719 d2edrho(3) = d2ersz(1)*d1rs*d1rs 720 & + d2ersz(2)*d1rs*d1zeta(2)*2.d0 721 & + d2ersz(3)*d1zeta(2)*d1zeta(2) 722 & + d1ersz(1)*d2rs 723 & + d1ersz(2)*d2zeta(3) 724 Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) 725 & + (2.d0*d1edrho(1) + rhoval*d2edrho(1))*fac 726 Amat2(n,D2_RA_RB) = Amat2(n,D2_RA_RB) 727 & + (d1edrho(1) + d1edrho(2) + rhoval*d2edrho(2))*fac 728 if (ipol.eq.2) 729 & Amat2(n,D2_RB_RB) = Amat2(n,D2_RB_RB) 730 & + (2.d0*d1edrho(2) + rhoval*d2edrho(3))*fac 731#endif 732#ifdef THIRD_DERIV 733c 1 = rsrsrs, 2 = rsrsz, 3 = rszz, 4 = zzz 734 d3ersz(1) = d3e(1) + fz*(d3e(2)-d3e(1)) 735 d3ersz(2) = d1fz*(d2e(2)-d2e(1)) 736 d3ersz(3) = d2fz*(d1e(2)-d1e(1)) 737 d3ersz(4) = d3fz*(e(2)-e(1)) 738c 739c 1 = aaa, 2 = aab, 3 = abb, 4 = bbb 740 d3edrho(1) = d3ersz(1)*d1rs*d1rs*d1rs 741 1 + d2ersz(1)*d1rs*d2rs*3.0d0 742 2 + d3ersz(3)*d1rs*d1zeta(1)*d1zeta(1)*3.0d0 743 3 + d2ersz(2)*d1rs*d2zeta(1)*3.0d0 744 4 + d1ersz(1)*d3rs 745 5 + d2ersz(2)*d1zeta(1)*d2rs*3.0d0 746 6 + d3ersz(2)*d1zeta(1)*d1rs*d1rs*3.0d0 747 7 + d3ersz(4)*d1zeta(1)*d1zeta(1)*d1zeta(1) 748 8 + d2ersz(3)*d1zeta(1)*d2zeta(1)*3.0d0 749 9 + d1ersz(2)*d3zeta(1) 750 d3edrho(2) = d3ersz(1)*d1rs*d1rs*d1rs 751 1 + d2ersz(1)*d1rs*d2rs*3.0d0 752 2 + d3ersz(3)*d1rs*(d1zeta(1)*d1zeta(1) 753 3 + d1zeta(1)*d1zeta(2)*2.0d0) 754 4 + d2ersz(2)*d1rs*(d2zeta(2)*2.0d0 755 5 + d2zeta(1)) 756 6 + d1ersz(1)*d3rs 757 7 + d2ersz(2)*d2rs*(d1zeta(1)*2.0d0 758 8 + d1zeta(2)) 759 9 + d3ersz(2)*d1rs*d1rs*(d1zeta(2) 760 A + d1zeta(1)*2.0d0) 761 B + d3ersz(4)*d1zeta(2)*d1zeta(1)*d1zeta(1) 762 C + d2ersz(3)*(d1zeta(1)*d2zeta(2)*2.0d0 763 D + d1zeta(2)*d2zeta(1)) 764 E + d1ersz(2)*d3zeta(2) 765 d3edrho(3) = d3ersz(1)*d1rs*d1rs*d1rs 766 1 + d2ersz(1)*d1rs*d2rs*3.0d0 767 2 + d3ersz(3)*d1rs*(d1zeta(2)*d1zeta(2) 768 3 + d1zeta(2)*d1zeta(1)*2.0d0) 769 4 + d2ersz(2)*d1rs*(d2zeta(2)*2.0d0 770 5 + d2zeta(3)) 771 6 + d1ersz(1)*d3rs 772 7 + d2ersz(2)*d2rs*(d1zeta(2)*2.0d0 773 8 + d1zeta(1)) 774 9 + d3ersz(2)*d1rs*d1rs*(d1zeta(1) 775 A + d1zeta(2)*2.0d0) 776 B + d3ersz(4)*d1zeta(1)*d1zeta(2)*d1zeta(2) 777 C + d2ersz(3)*(d1zeta(2)*d2zeta(2)*2.0d0 778 D + d1zeta(1)*d2zeta(3)) 779 E + d1ersz(2)*d3zeta(3) 780 d3edrho(4) = d3ersz(1)*d1rs*d1rs*d1rs 781 1 + d2ersz(1)*d1rs*d2rs*3.0d0 782 2 + d3ersz(3)*d1rs*d1zeta(2)*d1zeta(2)*3.0d0 783 3 + d2ersz(2)*d1rs*d2zeta(3)*3.0d0 784 4 + d1ersz(1)*d3rs 785 5 + d2ersz(2)*d1zeta(2)*d2rs*3.0d0 786 6 + d3ersz(2)*d1zeta(2)*d1rs*d1rs*3.0d0 787 7 + d3ersz(4)*d1zeta(2)*d1zeta(2)*d1zeta(2) 788 8 + d2ersz(3)*d1zeta(2)*d2zeta(3)*3.0d0 789 9 + d1ersz(2)*d3zeta(4) 790c 791 Amat3(n,D3_RA_RA_RA) = Amat3(n,D3_RA_RA_RA) 792 1 + ( 3.0d0*d2edrho(1) + rhoval*d3edrho(1) )*fac 793 Amat3(n,D3_RA_RA_RB) = Amat3(n,D3_RA_RA_RB) 794 1 + ( d2edrho(1) + 2.0d0*d2edrho(2) 795 2 + rhoval*d3edrho(2) )*fac 796 Amat3(n,D3_RA_RB_RB) = Amat3(n,D3_RA_RB_RB) 797 1 + ( 2.0d0*d2edrho(2) + d2edrho(3) 798 2 + rhoval*d3edrho(3) )*fac 799 if (ipol.eq.2) 800 1 Amat3(n,D3_RB_RB_RB) = Amat3(n,D3_RB_RB_RB) 801 2 + ( 3.0d0*d2edrho(3) + rhoval*d3edrho(4) )*fac 802#endif 803 endif 804 enddo 805 return 806 end 807c 808#ifndef SECOND_DERIV 809#define SECOND_DERIV 810c 811c Compile source again for the 2nd derivative case 812c 813#include "xc_perdew86.F" 814#endif 815c 816#ifndef THIRD_DERIV 817#define THIRD_DERIV 818c 819c Compile source again for the 3rd derivative case 820c 821#include "xc_perdew86.F" 822#endif 823