1c PBE exchange functional 2c 3c References: 4c [a] J.P. Perdew, K. Burke, and M. Ernzerhof, PRL 77, 3865 (1996). 5c [b] J.P. Perdew and Y. Wang, Phys. Rev. B 33, 8800 (1986).; 6c 40, 3399 (1989) (E). 7c Hammer, Hansen and Norskov, PRB 59, 7413 (1999) [RPBE] 8c Zhang and Yang, PRL 80, 890 (1998) [RevPBE] 9c 10#if !defined SECOND_DERIV && !defined THIRD_DERIV 11 Subroutine xc_xpbe96(whichf, 12 W tol_rho, fac, lfac, nlfac, rho, delrho, 13 & Amat, Cmat, nq, ipol, Ex, qwght,ldew,func) 14#elif defined(SECOND_DERIV) && !defined THIRD_DERIV 15 Subroutine xc_xpbe96_d2(whichf, 16 W tol_rho, fac, lfac, nlfac, rho, delrho, 17 & Amat, Amat2, Cmat, Cmat2, nq, ipol, Ex, 18 & qwght,ldew,func) 19#else 20 Subroutine xc_xpbe96_d3(whichf, 21 W tol_rho, fac, lfac, nlfac, rho, delrho, 22 & Amat, Amat2, Amat3, Cmat, Cmat2, Cmat3, 23 & nq, ipol, Ex, qwght,ldew,func) 24#endif 25c 26 implicit none 27c 28#include "dft2drv.fh" 29#include "dft3drv.fh" 30c 31 character*4 whichf 32 double precision fac, Ex 33 integer nq, ipol 34 logical lfac, nlfac,ldew 35 double precision func(*) ! value of the functional [output] 36c 37c Charge Density & Its Cube Root 38c 39 double precision rho(nq,ipol*(ipol+1)/2) 40c 41c Charge Density Gradient 42c 43 double precision delrho(nq,3,ipol) 44c 45c Quadrature Weights 46c 47 double precision qwght(nq) 48c 49c Sampling Matrices for the XC Potential & Energy 50c 51 double precision amat(nq,ipol), cmat(nq,*) 52c 53#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 54 double precision Amat2(nq,NCOL_AMAT2), Cmat2(nq,NCOL_CMAT2) 55#endif 56#ifdef THIRD_DERIV 57 double precision Amat3(nq,NCOL_AMAT3), Cmat3(nq,NCOL_CMAT3) 58#endif 59c 60 double precision tol_rho, pi, um, uk, umk,ukrev,umkrev 61 double precision C, Cs 62 double precision F43, F13, F23 63c 64#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 65 double precision F73 66#endif 67#ifdef THIRD_DERIV 68 double precision F10d3 69#endif 70 parameter(um=0.2195149727645171d0, uk=0.8040d0, umk=um/uk) 71 parameter(ukrev=1.245d0, umkrev=um/ukrev) 72 parameter (F43=4.d0/3.d0, F13=1.d0/3.d0, F23=2.0d0/3.0d0) 73#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 74 parameter (F73=7.d0/3.d0) 75#endif 76#ifdef THIRD_DERIV 77 parameter (F10d3=10.0d0/3.0d0) 78#endif 79c 80 integer n 81 double precision rrho, rho43, rho13, gamma, gam12, s, d1s(2), 82 & d, g, gp, d1g(2) 83#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 84 double precision rhom23, d2s(3), gpp, d2g(3) 85#endif 86#ifdef THIRD_DERIV 87 double precision d3s(4), d3g(4), rhom53, gppp 88#endif 89 double precision gpbe0, gpbe1, gpbe2, gpbe3 90 double precision grpbe0, grpbe1, grpbe2, grpbe3 91 double precision grevpbe0, grevpbe1, grevpbe2, grevpbe3 92c Original PBE 93 gpbe0(s)= uk*(1d0 - 1d0/(1d0+umk*s*s)) 94 gpbe1(s)= 2d0*um*s/(1d0+umk*s*s)**2 95 gpbe2(s)= 2d0*um*(1d0-4d0*umk*s*s/(1d0+umk*s*s))/(1d0+umk*s*s)**2 96 gpbe3(s)= 24.0d0*umk*um*s* 97 1 (2.0d0*umk*s*s/(1.0d0+umk*s*s)-1.0d0)/(1.0d0+umk*s*s)**3 98c revPBE by Zhang et al. 99 grevpbe0(s)= ukrev*(1d0 - 1d0/(1d0+umkrev*s*s)) 100 grevpbe1(s)= 2d0*um*s/(1d0+umkrev*s*s)**2 101 grevpbe2(s)= 2d0*um*(1d0-4d0*umkrev*s*s/(1d0+umkrev*s*s))/ 102 / (1d0+umkrev*s*s)**2 103 grevpbe3(s)= 24.0d0*umkrev*um*s* 104 1 (2.0d0*umkrev*s*s/(1.0d0+umkrev*s*s)-1.0d0)/ 105 2 (1.0d0+umkrev*s*s)**3 106c RPBE by Hammer et al. 107 grpbe0(s)= uk*(1d0 - exp(-umk*s*s)) 108 grpbe1(s)= 2d0*um*s*exp(-umk*s*s) 109 grpbe2(s)= 2d0*um*exp(-umk*s*s)*(1d0-2d0*umk*s*s) 110 grpbe3(s)= -4.0d0*umk*um*s*exp(-umk*s*s)*(3d0-2d0*umk*s*s) 111c 112 pi = acos(-1.d0) 113 C = -3d0/(4d0*pi)*(3d0*pi*pi)**F13 114 Cs = 0.5d0/(3d0*pi*pi)**F13 115 Cs = Cs * C ! account for including C in rho43 116c 117 if (ipol.eq.1 )then 118c 119c ======> SPIN-RESTRICTED <====== 120c 121c Daniel (9-28-12): There are somewhat mysterious coefficients involved 122c in the evaluation of the functional and its derivatives. We must 123c recall that the exchange energy is always written based on the 124c spin-scaling relationship for exchange: 125c 126c Ex[rho] = Ex[rho_a,rho_b] = 0.5*( Ex[2*rho_a] + Ex[2*rho_b] ) 127c 128c Thus, the electron density is always written: 129c rho -> 2*rho_s 130c gamma -> 4*gamma_ss 131c 132c Rationalization for the coefficients is mathematically justified below: 133c 134c ---------------------------- 135c Amat -> 0.5*2 = 1 136c Cmat -> 0.5*4 = 2 137c ---------------------------- 138c Amat2 -> 0.5*2*2 = 2 139c Cmat2(rg) -> 0.5*2*4 = 4 140c Cmat2(gg) -> 0.5*4*4 = 8 141c ---------------------------- 142c Amat3 -> 0.5*2*2*2 = 4 143c Cmat3(rrg) -> 0.5*2*2*4 = 8 144c Cmat3(rgg) -> 0.5*2*4*4 = 16 145c Cmat3(ggg) -> 0.5*4*4*4 = 32 146c ---------------------------- 147c 148c If, instead, the author of this code had decided to divide the total 149c density (rho(n,1)) by 2 in constructing the density and gamma, those 150c coefficients would be unnecessary. 151c 152#ifdef IFCV81 153CDEC$ NOSWP 154#endif 155 do 10 n = 1, nq 156 if (rho(n,1).lt.tol_rho) goto 10 157c#ifdef THIRD_DERIV 158c write(6,*) 'rho', rho(n,1) 159c#endif 160 rho43 = C*rho(n,1)**F43 161 rrho = 1d0/rho(n,1) 162 rho13 = F43*rho43*rrho 163#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 164 rhom23 = F13*rho13*rrho 165#endif 166#ifdef THIRD_DERIV 167 rhom53 = F23*rhom23*rrho 168#endif 169 if (lfac) then 170 Ex = Ex + rho43*qwght(n)*fac 171 if(ldew)func(n) = func(n) + rho43*fac 172 Amat(n,1) = Amat(n,1) + rho13*fac 173#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 174 Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) + 2d0*rhom23*fac 175#endif 176#ifdef THIRD_DERIV 177 Amat3(n,D3_RA_RA_RA) = Amat3(n,D3_RA_RA_RA) 178 1 - 4.0d0*rhom53*fac 179#endif 180 endif 181c 182 gamma = delrho(n,1,1)*delrho(n,1,1) + 183 & delrho(n,2,1)*delrho(n,2,1) + 184 & delrho(n,3,1)*delrho(n,3,1) 185 gam12 = dsqrt(gamma) 186 if (.not.(nlfac.and.gam12.gt.tol_rho**2)) goto 10 187c 188 189 s = Cs*gam12/rho43 190 d1s(1) = -F43*s*rrho 191 d1s(2) = 0.5d0*s/gamma 192c 193c Evaluate the GC part of F(s), i.e. g(s) = F(s) - 1 194c 195 if(whichf.eq.'revp') then 196 g=grevpbe0(s) 197 gp=grevpbe1(s) 198 elseif(whichf.eq.'rpbe') then 199 g=grpbe0(s) 200 gp=grpbe1(s) 201 else 202 g=gpbe0(s) 203 gp=gpbe1(s) 204 endif 205c 206c Daniel (7-27-12): gp is the derivative of the rational function, 207c or whatever the function in the revision is. 208c First derivatives of the enhancement factor 209 d1g(1) = gp*d1s(1) 210 d1g(2) = gp*d1s(2) 211 Ex = Ex + rho43*g*qwght(n)*fac 212 if(ldew)func(n) = func(n) + rho43*g*fac 213 Amat(n,1) = Amat(n,1) + (rho13*g+rho43*d1g(1))*fac 214 Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + 2d0*rho43*d1g(2)*fac 215#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 216 d2s(1) = -F73*d1s(1)*rrho 217 d2s(2) = -F43*d1s(2)*rrho 218 d2s(3) = -0.5d0*d1s(2)/gamma 219 if(whichf.eq.'revp') then 220 gpp=grevpbe2(s) 221 elseif(whichf.eq.'rpbe') then 222 gpp=grpbe2(s) 223 else 224 gpp=gpbe2(s) 225 endif 226c Second derivatives of the enhancement factor 227 d2g(1) = gp*d2s(1) + gpp*d1s(1)*d1s(1) 228 d2g(2) = gp*d2s(2) + gpp*d1s(1)*d1s(2) 229 d2g(3) = gp*d2s(3) + gpp*d1s(2)*d1s(2) 230 Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) 231 & +(rhom23*g 232 & + 2.d0*rho13*d1g(1) 233 & + rho43*d2g(1))*fac*2d0 234 Cmat2(n,D2_RA_GAA) = Cmat2(n,D2_RA_GAA) 235 & +(rho13*d1g(2) 236 & + rho43*d2g(2))*fac*4d0 237 Cmat2(n,D2_GAA_GAA) = Cmat2(n,D2_GAA_GAA) 238 & + rho43*d2g(3)*fac*8d0 239#endif 240#ifdef THIRD_DERIV 241c 1 = drdrdr, 2 = drdrdg, 3 = drdgdg, 4 = dgdgdg 242 d3s(1) = -F10d3*d2s(1)*rrho 243 d3s(2) = 0.5d0*d2s(1)/gamma 244 d3s(3) = -F43*d2s(3)*rrho 245 d3s(4) = -1.5d0*d2s(3)/gamma 246 if(whichf.eq.'revp') then 247 gppp = grevpbe3(s) 248 elseif(whichf.eq.'rpbe') then 249 gppp = grpbe3(s) 250 else 251 gppp = gpbe3(s) 252 endif 253c Third derivatives of the enhancement factor 254 d3g(1) = gp*d3s(1) + 3.0d0*gpp*d1s(1)*d2s(1) 255 1 + gppp*d1s(1)*d1s(1)*d1s(1) 256 d3g(2) = gp*d3s(2) 257 1 + gpp*d1s(2)*d2s(1) 258 2 + 2.0d0*gpp*d1s(1)*d2s(2) 259 3 + gppp*d1s(1)*d1s(1)*d1s(2) 260 d3g(3) = gp*d3s(3) 261 1 + gpp*d1s(1)*d2s(3) 262 2 + 2.0d0*gpp*d1s(2)*d2s(2) 263 3 + gppp*d1s(1)*d1s(2)*d1s(2) 264 d3g(4) = gp*d3s(4) + 3.0d0*gpp*d1s(2)*d2s(3) 265 1 + gppp*d1s(2)*d1s(2)*d1s(2) 266c 267 Amat3(n,D3_RA_RA_RA) = Amat3(n,D3_RA_RA_RA) 268 1 + (-rhom53*g 269 2 + 3.0d0*rhom23*d1g(1) 270 3 + 3.0d0*rho13*d2g(1) 271 4 + rho43*d3g(1))*fac*4.0d0 272 Cmat3(n,D3_RA_RA_GAA) = Cmat3(n,D3_RA_RA_GAA) 273 1 + (rhom23*d1g(2) 274 2 + 2.0d0*rho13*d2g(2) 275 3 + rho43*d3g(2))*fac*8.0d0 276 Cmat3(n,D3_RA_GAA_GAA) = Cmat3(n,D3_RA_GAA_GAA) 277 1 + (rho13*d2g(3) 278 2 + rho43*d3g(3))*fac*16.0d0 279 Cmat3(n,D3_GAA_GAA_GAA) = Cmat3(n,D3_GAA_GAA_GAA) 280 1 + (rho43*d3g(4))*fac*32.0d0 281#endif 282 10 continue 283c 284 else 285c 286c ======> SPIN-UNRESTRICTED <====== 287c 288c Daniel (9-28-12): There are somewhat mysterious coefficients involved 289c in the evaluation of the functional and its derivatives. We must 290c recall that the exchange energy is always written based on the 291c spin-scaling relationship for exchange: 292c 293c Ex[rho] = Ex[rho_a,rho_b] = 0.5*( Ex[2*rho_a] + Ex[2*rho_b] ) 294c 295c Thus, the electron density is always written: 296c rho -> 2*rho_s 297c gamma -> 4*gamma_ss 298c 299c It seems like the derivatives should be correctly balanced by the 300c following coefficients: 301c 302c ----------------------------- 303c Amat -> 0.5*2 = 1 304c Cmat -> 0.5*1 = 0.5 305c ----------------------------- 306c Amat2 -> 0.5*2*2 = 2 307c Cmat2(rg) -> 0.5*2*1 = 1 308c Cmat2(gg) -> 0.5*1*1 = 0.5 309c ----------------------------- 310c Amat3 -> 0.5*2*2*2 = 4 311c Cmat3(rrg) -> 0.5*2*2*1 = 2 312c Cmat3(rgg) -> 0.5*2*1*1 = 1 313c Cmat3(ggg) -> 0.5*1*1*1 = 0.5 314c ----------------------------- 315c 316#ifdef IFCV81 317CDEC$ NOSWP 318#endif 319 do 20 n = 1, nq 320 if (rho(n,1).lt.tol_rho) goto 20 321c 322c Alpha 323c 324 if (rho(n,2).lt.tol_rho) goto 25 325 rho43 = C*(2d0*rho(n,2))**F43 326 rrho = 0.5d0/rho(n,2) 327 rho13 = F43*rho43*rrho 328#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 329 rhom23 = F13*rho13*rrho 330#endif 331#ifdef THIRD_DERIV 332 rhom53 = F23*rhom23*rrho 333#endif 334 if (lfac) then 335 Ex = Ex + rho43*qwght(n)*fac*0.5d0 336 if(ldew)func(n) = func(n) + rho43*fac*0.5d0 337 Amat(n,1) = Amat(n,1) + rho13*fac 338#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 339 Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) + 2d0*rhom23*fac 340#endif 341#ifdef THIRD_DERIV 342 Amat3(n,D3_RA_RA_RA) = Amat3(n,D3_RA_RA_RA) 343 1 - 4.0d0*rhom53*fac 344#endif 345 endif 346c 347 gamma = delrho(n,1,1)*delrho(n,1,1) + 348 & delrho(n,2,1)*delrho(n,2,1) + 349 & delrho(n,3,1)*delrho(n,3,1) 350 gam12 = 2d0*dsqrt(gamma) 351 if (.not.(nlfac.and.gam12.gt.tol_rho**2)) goto 25 352c 353 s = Cs*gam12/rho43 354 d1s(1) = -F43*s*rrho 355 d1s(2) = 0.5d0*s/gamma 356c 357c Evaluate the GC part of F(s), i.e. g(s) = F(s) - 1 358c 359 360 if(whichf.eq.'revp') then 361 g=grevpbe0(s) 362 gp=grevpbe1(s) 363 elseif(whichf.eq.'rpbe') then 364 g=grpbe0(s) 365 gp=grpbe1(s) 366 else 367 g=gpbe0(s) 368 gp=gpbe1(s) 369 endif 370c Daniel (9-28-12): Factors of 2 are inconsistent with the restricted 371c calculations because a gam12 is doubled above and 372 d1g(1) = gp*d1s(1) 373 d1g(2) = gp*d1s(2) 374 Ex = Ex + rho43*g*qwght(n)*fac*0.5d0 375 if(ldew)func(n) = func(n) + rho43*g*fac*0.5d0 376 Amat(n,1) = Amat(n,1) + (rho13*g+rho43*d1g(1))*fac 377 Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + 0.5d0*rho43*d1g(2)*fac 378#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 379 d2s(1) = -F73*d1s(1)*rrho 380 d2s(2) = -F43*d1s(2)*rrho 381 d2s(3) = -0.5d0*d1s(2)/gamma 382 if(whichf.eq.'revp') then 383 gpp=grevpbe2(s) 384 elseif(whichf.eq.'rpbe') then 385 gpp=grpbe2(s) 386 else 387 gpp=gpbe2(s) 388 endif 389 d2g(1) = gp*d2s(1) + gpp*d1s(1)*d1s(1) 390 d2g(2) = gp*d2s(2) + gpp*d1s(1)*d1s(2) 391 d2g(3) = gp*d2s(3) + gpp*d1s(2)*d1s(2) 392c 393 Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) 394 & +(rhom23*g 395 & + 2.d0*rho13*d1g(1) 396 & + rho43*d2g(1))*fac*2d0 397 Cmat2(n,D2_RA_GAA) = Cmat2(n,D2_RA_GAA) 398 & +(rho13*d1g(2) 399 & + rho43*d2g(2))*fac 400 Cmat2(n,D2_GAA_GAA) = Cmat2(n,D2_GAA_GAA) 401 & + rho43*d2g(3)*fac*0.5d0 402#endif 403#ifdef THIRD_DERIV 404c 1 = drdrdr, 2 = drdrdg, 3 = drdgdg, 4 = dgdgdg 405 d3s(1) = -F10d3*d2s(1)*rrho 406 d3s(2) = 0.5d0*d2s(1)/gamma 407 d3s(3) = -F43*d2s(3)*rrho 408 d3s(4) = -1.5d0*d2s(3)/gamma 409 if(whichf.eq.'revp') then 410 gppp = grevpbe3(s) 411 elseif(whichf.eq.'rpbe') then 412 gppp = grpbe3(s) 413 else 414 gppp = gpbe3(s) 415 endif 416c Third derivatives of the enhancement factor 417 d3g(1) = gp*d3s(1) + 3.0d0*gpp*d1s(1)*d2s(1) 418 1 + gppp*d1s(1)*d1s(1)*d1s(1) 419 d3g(2) = gp*d3s(2) 420 1 + gpp*d1s(2)*d2s(1) 421 2 + 2.0d0*gpp*d1s(1)*d2s(2) 422 3 + gppp*d1s(1)*d1s(1)*d1s(2) 423 d3g(3) = gp*d3s(3) 424 1 + gpp*d1s(1)*d2s(3) 425 2 + 2.0d0*gpp*d1s(2)*d2s(2) 426 3 + gppp*d1s(1)*d1s(2)*d1s(2) 427 d3g(4) = gp*d3s(4) + 3.0d0*gpp*d1s(2)*d2s(3) 428 1 + gppp*d1s(2)*d1s(2)*d1s(2) 429c 430 Amat3(n,D3_RA_RA_RA) = Amat3(n,D3_RA_RA_RA) 431 1 + (-rhom53*g 432 2 + 3.0d0*rhom23*d1g(1) 433 3 + 3.0d0*rho13*d2g(1) 434 4 + rho43*d3g(1))*fac*4.0d0 435 Cmat3(n,D3_RA_RA_GAA) = Cmat3(n,D3_RA_RA_GAA) 436 1 + (rhom23*d1g(2) 437 2 + 2.0d0*rho13*d2g(2) 438 3 + rho43*d3g(2))*fac*2.0d0 439 Cmat3(n,D3_RA_GAA_GAA) = Cmat3(n,D3_RA_GAA_GAA) 440 1 + (rho13*d2g(3) 441 2 + rho43*d3g(3))*fac 442 Cmat3(n,D3_GAA_GAA_GAA) = Cmat3(n,D3_GAA_GAA_GAA) 443 1 + (rho43*d3g(4))*fac*0.5d0 444#endif 445c 446c Beta 447c 448 25 continue 449 if (rho(n,3).lt.tol_rho) goto 20 450 rho43 = C*(2d0*rho(n,3))**F43 451 rrho = 0.5d0/rho(n,3) 452 rho13 = F43*rho43*rrho 453#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 454 rhom23 = F13*rho13*rrho 455#endif 456#ifdef THIRD_DERIV 457 rhom53 = F23*rhom23*rrho 458#endif 459 if (lfac) then 460 Ex = Ex + rho43*qwght(n)*fac*0.5d0 461 if(ldew)func(n) = func(n) + rho43*fac*0.5d0 462 Amat(n,2) = Amat(n,2) + rho13*fac 463#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 464 Amat2(n,D2_RB_RB) = Amat2(n,D2_RB_RB) + 2d0*rhom23*fac 465#endif 466#ifdef THIRD_DERIV 467 Amat3(n,D3_RB_RB_RB) = Amat3(n,D3_RB_RB_RB) 468 1 - 4.0d0*rhom53*fac 469#endif 470 endif 471c 472 gamma = delrho(n,1,2)*delrho(n,1,2) + 473 & delrho(n,2,2)*delrho(n,2,2) + 474 & delrho(n,3,2)*delrho(n,3,2) 475 gam12 = 2d0*dsqrt(gamma) 476 if (.not.(nlfac.and.gam12.gt.tol_rho**2)) goto 20 477c 478 s = Cs*gam12/rho43 479 d1s(1) = -F43*s*rrho 480 d1s(2) = 0.5d0*s/gamma 481c 482c Evaluate the GC part of F(s), i.e. g(s) = F(s) - 1 483c 484 if(whichf.eq.'revp') then 485 g=grevpbe0(s) 486 gp=grevpbe1(s) 487 elseif(whichf.eq.'rpbe') then 488 g=grpbe0(s) 489 gp=grpbe1(s) 490 else 491 g=gpbe0(s) 492 gp=gpbe1(s) 493 endif 494c 495 d1g(1) = gp*d1s(1) 496 d1g(2) = gp*d1s(2) 497 Ex = Ex + rho43*g*qwght(n)*fac*0.5d0 498 if(ldew)func(n) = func(n) + rho43*g*fac*0.5d0 499 Amat(n,2) = Amat(n,2) + (rho13*g+rho43*d1g(1))*fac 500 Cmat(n,D1_GBB) = Cmat(n,D1_GBB) + 0.5d0*rho43*d1g(2)*fac 501#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 502 d2s(1) = -F73*d1s(1)*rrho 503 d2s(2) = -F43*d1s(2)*rrho 504 d2s(3) = -0.5d0*d1s(2)/gamma 505 if(whichf.eq.'revp') then 506 gpp=grevpbe2(s) 507 elseif(whichf.eq.'rpbe') then 508 gpp=grpbe2(s) 509 else 510 gpp=gpbe2(s) 511 endif 512 d2g(1) = gp*d2s(1) + gpp*d1s(1)*d1s(1) 513 d2g(2) = gp*d2s(2) + gpp*d1s(1)*d1s(2) 514 d2g(3) = gp*d2s(3) + gpp*d1s(2)*d1s(2) 515 Amat2(n,D2_RB_RB) = Amat2(n,D2_RB_RB) 516 & +(rhom23*g 517 & + 2.d0*rho13*d1g(1) 518 & + rho43*d2g(1))*fac*2d0 519 Cmat2(n,D2_RB_GBB) = Cmat2(n,D2_RB_GBB) 520 & +(rho13*d1g(2) 521 & + rho43*d2g(2))*fac 522 Cmat2(n,D2_GBB_GBB) = Cmat2(n,D2_GBB_GBB) 523 & + rho43*d2g(3)*fac*0.5d0 524#endif 525#ifdef THIRD_DERIV 526c 1 = drdrdr, 2 = drdrdg, 3 = drdgdg, 4 = dgdgdg 527 d3s(1) = -F10d3*d2s(1)*rrho 528 d3s(2) = 0.5d0*d2s(1)/gamma 529 d3s(3) = -F43*d2s(3)*rrho 530 d3s(4) = -1.5d0*d2s(3)/gamma 531 if(whichf.eq.'revp') then 532 gppp = grevpbe3(s) 533 elseif(whichf.eq.'rpbe') then 534 gppp = grpbe3(s) 535 else 536 gppp = gpbe3(s) 537 endif 538c Third derivatives of the enhancement factor 539 d3g(1) = gp*d3s(1) + 3.0d0*gpp*d1s(1)*d2s(1) 540 1 + gppp*d1s(1)*d1s(1)*d1s(1) 541 d3g(2) = gp*d3s(2) 542 1 + gpp*d1s(2)*d2s(1) 543 2 + 2.0d0*gpp*d1s(1)*d2s(2) 544 3 + gppp*d1s(1)*d1s(1)*d1s(2) 545 d3g(3) = gp*d3s(3) 546 1 + gpp*d1s(1)*d2s(3) 547 2 + 2.0d0*gpp*d1s(2)*d2s(2) 548 3 + gppp*d1s(1)*d1s(2)*d1s(2) 549 d3g(4) = gp*d3s(4) + 3.0d0*gpp*d1s(2)*d2s(3) 550 1 + gppp*d1s(2)*d1s(2)*d1s(2) 551c 552 Amat3(n,D3_RB_RB_RB) = Amat3(n,D3_RB_RB_RB) 553 1 + (-rhom53*g 554 2 + 3.0d0*rhom23*d1g(1) 555 3 + 3.0d0*rho13*d2g(1) 556 4 + rho43*d3g(1))*fac*4.0d0 557 Cmat3(n,D3_RB_RB_GBB) = Cmat3(n,D3_RB_RB_GBB) 558 1 + (rhom23*d1g(2) 559 2 + 2.0d0*rho13*d2g(2) 560 3 + rho43*d3g(2))*fac*2.0d0 561 Cmat3(n,D3_RB_GBB_GBB) = Cmat3(n,D3_RB_GBB_GBB) 562 1 + (rho13*d2g(3) 563 2 + rho43*d3g(3))*fac 564 Cmat3(n,D3_GBB_GBB_GBB) = Cmat3(n,D3_GBB_GBB_GBB) 565 1 + (rho43*d3g(4))*fac*0.5d0 566#endif 567c 568 20 continue 569 endif 570c 571 return 572 end 573#ifndef SECOND_DERIV 574#define SECOND_DERIV 575c 576c Compile source again for the 2nd derivative case 577c 578#include "xc_xpbe96.F" 579#endif 580#ifndef THIRD_DERIV 581#define THIRD_DERIV 582c 583c Compile source again for the 3rd derivative case 584c 585#include "xc_xpbe96.F" 586#endif 587c $Id$ 588