1#ifndef SECOND_DERIV 2 Subroutine xc_spbe96(tol_rho, cfac, lcfac, nlcfac, rho, delrho, 3 & Amat, Cmat, nq, ipol, Ec, qwght, 4 & ldew,ffunc) 5#else 6 Subroutine xc_spbe96_d2(tol_rho, cfac, lcfac, nlcfac, 7 & rho, delrho, Amat, Amat2, Cmat, Cmat2, 8 & nq, ipol, Ec, qwght, ldew, ffunc) 9#endif 10c 11c$Id$ 12c 13 Implicit none 14#include "dft2drv.fh" 15c 16c 17c Input and other parameters 18c 19 integer lnq 20 integer ipol, nq 21 double precision dummy(1) 22 double precision cfac(*) 23 logical lcfac(*), nlcfac(*), ldew 24 logical lfac, nlfac, lfacl, nlfacl 25 double precision ffunc(*) 26 double precision fac, facl 27 double precision lqwght 28 double precision tol_rho 29c 30c Constants in PBE functional 31c 32 double precision GAMMA, BETA, PI 33 parameter (GAMMA = 0.03109069086965489503494086371273d0) 34 parameter (BETA = 0.06672455060314922d0) 35 parameter (PI = 3.1415926535897932385d0) 36c 37c Threshold parameters 38c 39 double precision TOLL, EXPTOL 40 double precision EPS 41 parameter (TOLL = 1.0D-40, EXPTOL = 40.0d0) 42 parameter (EPS = 1.0e-8) 43c 44c Correlation energy 45c 46 double precision Ec 47c 48c Charge Density 49c 50 double precision rho(nq,ipol*(ipol+1)/2) 51 double precision rho_t(3) 52c 53c Charge Density Gradient 54c 55 double precision delrho(nq,3,ipol) 56 double precision dsqgamma 57c 58c Quadrature Weights 59c 60 double precision qwght(nq) 61c 62c Sampling Matrices for the XC Potential 63c 64 double precision Amat(nq,ipol), Cmat(nq,*) 65#ifdef SECOND_DERIV 66c 67c Sampling Matrices for the XC Kernel 68c 69 double precision Amat2(nq,NCOL_AMAT2), Cmat2(nq,NCOL_CMAT2) 70#endif 71c 72c Intermediate derivative results, etc. 73c 74 integer n 75 double precision rhoval, gammaval 76 double precision nepsc, dnepscdn(2) 77 double precision epsc, depscdna, depscdnb 78 double precision H0, dH0dna, dH0dnb, dH0dg 79 double precision phi, dphidna, dphidnb, dphidzeta 80 double precision zeta, dzetadna, dzetadnb 81 double precision arglog, darglogdna, darglogdnb, darglogdg 82 double precision fAt, dfAtdt, dfAtdA 83 double precision fAtnum 84 double precision fAtden, dfAtdendt, dfAtdendA 85 double precision dfAtdna, dfAtdnb, dfAtdg 86 double precision A, dAdna, dAdnb 87 double precision t, dtdna, dtdnb, dtdg 88 double precision ks, dksdna, dksdnb 89 double precision argexp, dargexpdna, dargexpdnb 90 double precision expinA 91#ifdef SECOND_DERIV 92 double precision d2nepscdn2(NCOL_AMAT2) 93 double precision d2epscdna2, d2epscdnadnb, d2epscdnb2 94 double precision d2H0dna2, d2H0dnadnb, d2H0dnb2 95 double precision d2H0dnadg, d2H0dnbdg, d2H0dg2 96 double precision d2phidzeta2, d2phidna2, d2phidnadnb, d2phidnb2 97 double precision d2zetadna2, d2zetadnadnb, d2zetadnb2 98 double precision d2arglogdna2, d2arglogdnb2, d2arglogdnadnb 99 double precision d2arglogdnadg, d2arglogdnbdg, d2arglogdg2 100 double precision d2fAtdt2, d2fAtdA2, d2fAtdtdA, d2fAtdg2 101 double precision d2fAtdendt2, d2fAtdendtdA, d2fAtdendA2 102 double precision d2fAtdna2, d2fAtdnb2, d2fAtdnadnb 103 double precision d2fAtdnadg, d2fAtdnbdg 104 double precision d2Adna2, d2Adnadnb, d2Adnb2 105 double precision d2tdna2, d2tdnb2, d2tdnadnb 106 double precision d2tdg2, d2tdnadg, d2tdnbdg 107 double precision d2ksdna2, d2ksdnb2, d2ksdnadnb 108 double precision d2argexpdna2, d2argexpdnb2, d2argexpdnadnb 109#endif 110c 111c References: 112c [a] J. P. Perdew, K. Burke, and M. Ernzerhof, 113c {\it Generalized gradient approximation made simple}, 114c Phys.\ Rev.\ Lett. {\bf 77,} 3865 (1996). 115c [b] J. P. Perdew, K. Burke, and Y. Wang, {\it Real-space cutoff 116c construction of a generalized gradient approximation: The PW91 117c density functional}, submitted to Phys.\ Rev.\ B, Feb. 1996. 118c [c] J. P. Perdew and Y. Wang, Phys.\ Rev.\ B {\bf 45}, 13244 (1992). 119c 120c E_c(PBE) = Int n (epsilon_c + H0) dxdydz 121c 122c n*epsilon_c <=== supplied by another subroutine 123c d(n*epsilon_c)/d(na) <=== supplied by another subroutine 124c d2(n*epsilon_c)/d(na)d(na) <=== supplied by another subroutine 125c d2(n*epsilon_c)/d(na)d(nb) <=== supplied by another subroutine 126c d2(n*epsilon_c)/d(nb)d(nb) <=== supplied by another subroutine 127c 128c H0 = GAMMA * phi**3 * log{ 1 + BETA/GAMMA * t**2 * [ ... ]} 129c 130c phi = (1/2)[(1+zeta)**(2/3)+(1-zeta)**(2/3)] 131c 132c zeta = (na - nb)/n 133c 134c [ ... ] = (1 + A * t**2)/(1 + A * t**2 + A**2 * t**4) 135c 136c A = BETA/GAMMA [exp{-epsilon_c/(GAMMA*phi**3)}-1]**(-1) 137c 138c t = |Nabla n|/(2*phi*ks*n) 139c 140c ks = 2 * (3 * PI**2 * n)**(1/6) / sqrt(PI) 141c 142c |Nabla n| = sqrt(g_aa + g_bb + 2*g_ab) 143c 144c Names of variables 145c 146c E_c(PBE) : Ec 147c n (alpha+beta density) : rhoval 148c na, nb : rho(*,2), rho(*,3) 149c epsilon_c : epsc 150c H0 : H0 151c n*epsilon_c : nepsc 152c phi : phi 153c zeta : zeta 154c { ... } : arglog 155c [ ... ] : fAt 156c (1 + A * t**2) : fAtnum 157c (1 + A * t**2 + A**2 * t**4) : fAtden 158c A : A 159c t : t 160c |Nabla n| : gammaval 161c ks : ks 162c {-epsilon_c ... } : argexp 163c g_aa, g_bb, g_ab : g 164c 165c Derivatives of these are named like d...dna, d2...dnadnb, 166c d2...dna2, etc. 167c 168 fac = cfac(46) 169 lfac = lcfac(46) 170 nlfac = nlcfac(46) 171c 172c ======> BOTH SPIN-RESTRICETED AND UNRESTRICTED <====== 173c 174 do 20 n = 1, nq 175c 176c n and zeta = (na - nb)/n 177c 178 rhoval = rho(n,1) 179 rho_t(1) = rho(n,1) 180 if (ipol.eq.2) then 181 rho_t(2) = rho(n,2) 182 rho_t(3) = rho(n,3) 183 endif 184 if (rhoval.le.tol_rho) goto 20 185 if (ipol.eq.1) then 186 gammaval = delrho(n,1,1)*delrho(n,1,1) + 187 & delrho(n,2,1)*delrho(n,2,1) + 188 & delrho(n,3,1)*delrho(n,3,1) 189 else 190 gammaval = delrho(n,1,1)*delrho(n,1,1) + 191 & delrho(n,1,2)*delrho(n,1,2) + 192 & delrho(n,2,1)*delrho(n,2,1) + 193 & delrho(n,2,2)*delrho(n,2,2) + 194 & delrho(n,3,1)*delrho(n,3,1) + 195 & delrho(n,3,2)*delrho(n,3,2) + 196 & 2.d0*(delrho(n,1,1)*delrho(n,1,2) + 197 & delrho(n,2,1)*delrho(n,2,2) + 198 & delrho(n,3,1)*delrho(n,3,2)) 199 endif 200 dsqgamma = max(dsqrt(gammaval),tol_rho) 201 nepsc = 0.0d0 202 dnepscdn(1) = 0.0d0 203 if (ipol.eq.2) dnepscdn(2) = 0.0d0 204#ifdef SECOND_DERIV 205 d2nepscdn2(D2_RA_RA)=0.0d0 206 d2nepscdn2(D2_RA_RB)=0.0d0 207 if (ipol.eq.2) d2nepscdn2(D2_RB_RB)=0.0d0 208#endif 209c 210c call for LDA bit 211c 212 lnq = 1 213 lqwght = 1.0d0 214c 215 if (abs(cfac(1)).gt.EPS) then 216 facl = cfac(1) 217#ifndef SECOND_DERIV 218 call xc_vwn_5(tol_rho,facl,rho_t, 219 & dnepscdn,lnq,ipol,nepsc,lqwght,.false.,dummy) 220#else 221 call xc_vwn_5_d2(tol_rho,facl,rho_t, 222 & dnepscdn,d2nepscdn2,lnq,ipol,nepsc,lqwght, 223 & .false.,dummy) 224#endif 225 endif 226c 227 if (abs(cfac(3)).gt.EPS) then 228 facl = cfac(3) 229 lfacl = lcfac(3) 230 nlfacl = nlcfac(3) 231#ifndef SECOND_DERIV 232 call xc_p81(tol_rho,facl,lfacl,nlfacl,rho_t,dnepscdn, 233 & lnq,ipol,nepsc,lqwght,.false.,dummy) 234#else 235 call xc_p81_d2(tol_rho,facl,lfacl,nlfacl,rho_t, 236 & dnepscdn,d2nepscdn2,lnq,ipol,nepsc,lqwght, 237 & .false.,dummy) 238#endif 239 endif 240c 241 if (abs(cfac(6)).gt.EPS.or.lfac) then 242 facl = cfac(6) 243 if (abs(facl).lt.EPS) facl = 1.0d0 244 lfacl = lcfac(6) 245 nlfacl = nlcfac(6) 246#ifndef SECOND_DERIV 247 call xc_pw91lda(tol_rho,facl,lfacl,nlfacl,rho_t, 248 & dnepscdn,lnq,ipol,nepsc,lqwght, 249 & .false.,dummy) 250#else 251 call xc_pw91lda_d2(tol_rho,facl,lfacl,nlfacl,rho_t, 252 & dnepscdn,d2nepscdn2,lnq,ipol,nepsc,lqwght, 253 & .false.,dummy) 254#endif 255 endif 256c 257 if (abs(cfac(7)).gt.EPS) then 258 facl = cfac(7) 259#ifndef SECOND_DERIV 260 call xc_vwn_1_rpa(tol_rho,facl,rho_t, 261 & dnepscdn,lnq,ipol,nepsc,lqwght, 262 & .false.,dummy) 263#else 264 call xc_vwn_1_rpa_d2(tol_rho,facl,rho_t, 265 & dnepscdn,d2nepscdn2,lnq,ipol,nepsc,lqwght, 266 & .false.,dummy) 267#endif 268 endif 269c 270 if (abs(cfac(8)).gt.EPS) then 271 facl = cfac(8) 272#ifndef SECOND_DERIV 273 call xc_vwn_1(tol_rho,facl,rho_t, 274 & dnepscdn,lnq,ipol,nepsc,lqwght, 275 & .false.,dummy) 276#else 277 call xc_vwn_1_d2(tol_rho,facl,rho_t, 278 & dnepscdn,d2nepscdn2,lnq,ipol,nepsc,lqwght, 279 & .false.,dummy) 280#endif 281 endif 282c 283 if (abs(cfac(9)).gt.EPS) then 284 facl = cfac(9) 285#ifndef SECOND_DERIV 286 call xc_vwn_2(tol_rho,facl,rho_t, 287 & dnepscdn,lnq,ipol,nepsc,lqwght, 288 & .false.,dummy) 289#else 290 call xc_vwn_2_d2(tol_rho,facl,rho_t, 291 & dnepscdn,d2nepscdn2,lnq,ipol,nepsc,lqwght, 292 & .false.,dummy) 293#endif 294 endif 295c 296 if (abs(cfac(10)).gt.EPS) then 297 facl = cfac(10) 298#ifndef SECOND_DERIV 299 call xc_vwn_3(tol_rho,facl,rho_t, 300 & dnepscdn,lnq,ipol,nepsc,lqwght, 301 & .false.,dummy) 302#else 303 call xc_vwn_3_d2(tol_rho,facl,rho_t, 304 & dnepscdn,d2nepscdn2,lnq,ipol,nepsc,lqwght, 305 & .false.,dummy) 306#endif 307 endif 308c 309 if (abs(cfac(11)).gt.EPS) then 310 facl = cfac(11) 311#ifndef SECOND_DERIV 312 call xc_vwn_4(tol_rho,facl,rho_t, 313 & dnepscdn,lnq,ipol,nepsc,lqwght, 314 & .false.,dummy) 315#else 316 call xc_vwn_4_d2(tol_rho,facl,rho_t, 317 & dnepscdn,d2nepscdn2,lnq,ipol,nepsc,lqwght, 318 & .false.,dummy) 319#endif 320 endif 321c ================== 322c PBE non-local part 323c ================== 324 if(abs(nepsc).lt.tol_rho*tol_rho) goto 20 325c 326c epsilon_c = n*epsilon_c / n 327c 328 epsc = nepsc/rhoval 329 if (ipol.eq.1) then 330 depscdna = dnepscdn(1)/rhoval-nepsc/(rhoval**2) 331 depscdnb = depscdna 332 else 333 depscdna = dnepscdn(1)/rhoval-nepsc/(rhoval**2) 334 depscdnb = dnepscdn(2)/rhoval-nepsc/(rhoval**2) 335 endif 336#ifdef SECOND_DERIV 337 if (ipol.eq.1) then 338 d2epscdna2 = d2nepscdn2(D2_RA_RA)/rhoval 339 & -dnepscdn(1)/(rhoval**2) 340 & -dnepscdn(1)/(rhoval**2) 341 & +2.0d0*nepsc/(rhoval**3) 342 d2epscdnadnb = d2nepscdn2(D2_RA_RB)/rhoval 343 & -dnepscdn(1)/(rhoval**2) 344 & -dnepscdn(1)/(rhoval**2) 345 & +2.0d0*nepsc/(rhoval**3) 346 d2epscdnb2 = d2epscdna2 347 else 348 d2epscdna2 = d2nepscdn2(D2_RA_RA)/rhoval 349 & -dnepscdn(1)/(rhoval**2) 350 & -dnepscdn(1)/(rhoval**2) 351 & +2.0d0*nepsc/(rhoval**3) 352 d2epscdnadnb = d2nepscdn2(D2_RA_RB)/rhoval 353 & -dnepscdn(1)/(rhoval**2) 354 & -dnepscdn(2)/(rhoval**2) 355 & +2.0d0*nepsc/(rhoval**3) 356 d2epscdnb2 = d2nepscdn2(D2_RB_RB)/rhoval 357 & -dnepscdn(2)/(rhoval**2) 358 & -dnepscdn(2)/(rhoval**2) 359 & +2.0d0*nepsc/(rhoval**3) 360 endif 361#endif 362c 363c ks = 2*(3*PI**2*n)**(1/6)/sqrt(PI) and its derivs 364c 365 ks = 2.0d0*(3.0d0*PI*PI*rhoval)**(1.0d0/6.0d0)/dsqrt(PI) 366 dksdna = (1.0d0/6.0d0)*ks/rhoval 367 dksdnb = dksdna 368#ifdef SECOND_DERIV 369 d2ksdna2 = (1.0d0/6.0d0)*dksdna/rhoval 370 & - (1.0d0/6.0d0)*ks/(rhoval**2) 371 d2ksdnadnb = d2ksdna2 372 d2ksdnb2 = d2ksdna2 373#endif 374c 375c zeta = (na-nb)/n and its derivs 376c 377 if (ipol.eq.1) then 378 zeta = 0.0d0 379 else 380 zeta = (rho(n,2)-rho(n,3))/rhoval 381 endif 382 if(zeta.lt.-1.0d0) zeta=-1.0d0 383 if(zeta.gt. 1.0d0) zeta= 1.0d0 384 if (ipol.eq.1) then 385 dzetadna = 1.0d0/rhoval 386 dzetadnb = -1.0d0/rhoval 387#ifdef SECOND_DERIV 388 d2zetadna2 = -2.0d0/(rhoval**2) 389 d2zetadnadnb = 0.0d0 390 d2zetadnb2 = -2.0d0/(rhoval**2) 391#endif 392 else 393 dzetadna = 2.0d0*rho(n,3)/(rhoval**2) 394 dzetadnb = -2.0d0*rho(n,2)/(rhoval**2) 395#ifdef SECOND_DERIV 396 d2zetadna2 = -4.0d0*rho(n,3)/(rhoval**3) 397 d2zetadnadnb = 2.0d0*(rho(n,2)-rho(n,3))/(rhoval**3) 398 d2zetadnb2 = -4.0d0*rho(n,2)/(rhoval**3) 399#endif 400 endif 401c 402c phi = (1/2)[(1+zeta)**(2/3)+(1-zeta)**(2/3)] and its derivs 403c 404 phi = 0.5d0*((1.0d0+zeta)**(2.0d0/3.0d0) 405 & +(1.0d0-zeta)**(2.0d0/3.0d0)) 406 if ((1.0d0-zeta).lt.tol_rho) then 407 dphidzeta = 0.5d0*(2.0d0/3.0d0)*( 408 & (1.0d0+zeta)**(2.0d0/3.0d0)/(1.0d0+zeta)) 409 else if ((1.0d0+zeta).lt.tol_rho) then 410 dphidzeta = 0.5d0*(2.0d0/3.0d0)*( 411 & -(1.0d0-zeta)**(2.0d0/3.0d0)/(1.0d0-zeta)) 412 else 413 dphidzeta = 0.5d0*(2.0d0/3.0d0)*( 414 & (1.0d0+zeta)**(2.0d0/3.0d0)/(1.0d0+zeta) 415 & -(1.0d0-zeta)**(2.0d0/3.0d0)/(1.0d0-zeta)) 416 endif 417 dphidna = dphidzeta*dzetadna 418 dphidnb = dphidzeta*dzetadnb 419#ifdef SECOND_DERIV 420 if ((1.0d0-zeta).lt.tol_rho) then 421 d2phidzeta2 = -(1.0d0/9.0d0)*( 422 & (1.0d0+zeta)**(2.0d0/3.0d0)/((1.0d0+zeta)**2)) 423 else if ((1.0d0+zeta).lt.tol_rho) then 424 d2phidzeta2 = -(1.0d0/9.0d0)*( 425 & (1.0d0-zeta)**(2.0d0/3.0d0)/((1.0d0-zeta)**2)) 426 else 427 d2phidzeta2 = -(1.0d0/9.0d0)*( 428 & (1.0d0+zeta)**(2.0d0/3.0d0)/((1.0d0+zeta)**2) 429 & +(1.0d0-zeta)**(2.0d0/3.0d0)/((1.0d0-zeta)**2)) 430 endif 431 d2phidna2 = d2phidzeta2*dzetadna*dzetadna 432 & + dphidzeta*d2zetadna2 433 d2phidnadnb = d2phidzeta2*dzetadna*dzetadnb 434 & + dphidzeta*d2zetadnadnb 435 d2phidnb2 = d2phidzeta2*dzetadnb*dzetadnb 436 & + dphidzeta*d2zetadnb2 437#endif 438c 439c t = |Nabla n|/(2*phi*ks*n) and its derivs 440c 441 t = dsqgamma/(2.0d0*phi*ks*rhoval) 442 dtdna = -t/rhoval-t/phi*dphidna-t/ks*dksdna 443 dtdnb = -t/rhoval-t/phi*dphidnb-t/ks*dksdnb 444#ifdef SECOND_DERIV 445 d2tdna2 = - dtdna/rhoval 446 & + t/(rhoval**2) 447 & - dtdna/phi*dphidna 448 & + t/(phi**2)*(dphidna**2) 449 & - t/phi*d2phidna2 450 & - dtdna/ks*dksdna 451 & + t/(ks**2)*(dksdna**2) 452 & - t/ks*d2ksdna2 453 d2tdnadnb = - dtdnb/rhoval 454 & + t/(rhoval**2) 455 & - dtdnb/phi*dphidna 456 & + t/(phi**2)*(dphidna*dphidnb) 457 & - t/phi*d2phidnadnb 458 & - dtdnb/ks*dksdna 459 & + t/(ks**2)*(dksdna*dksdnb) 460 & - t/ks*d2ksdnadnb 461 d2tdnb2 = - dtdna/rhoval 462 & + t/(rhoval**2) 463 & - dtdnb/phi*dphidnb 464 & + t/(phi**2)*(dphidnb**2) 465 & - t/phi*d2phidnb2 466 & - dtdnb/ks*dksdnb 467 & + t/(ks**2)*(dksdnb**2) 468 & - t/ks*d2ksdnb2 469#endif 470c 471c { ... } in A (see below) and its derivs 472c 473 argexp = -epsc/GAMMA/(phi**3) 474 dargexpdna = -depscdna/GAMMA/(phi**3) 475 & +3.0d0*epsc/GAMMA/(phi**4)*dphidna 476 dargexpdnb = -depscdnb/GAMMA/(phi**3) 477 & +3.0d0*epsc/GAMMA/(phi**4)*dphidnb 478#ifdef SECOND_DERIV 479 d2argexpdna2 = -d2epscdna2/GAMMA/(phi**3) 480 & +3.0d0*depscdna/GAMMA/(phi**4)*dphidna 481 & +3.0d0*depscdna/GAMMA/(phi**4)*dphidna 482 & -12.0d0*epsc/GAMMA/(phi**5)*dphidna**2 483 & +3.0d0*epsc/GAMMA/(phi**4)*d2phidna2 484 d2argexpdnadnb = -d2epscdnadnb/GAMMA/(phi**3) 485 & +3.0d0*depscdna/GAMMA/(phi**4)*dphidnb 486 & +3.0d0*depscdnb/GAMMA/(phi**4)*dphidna 487 & -12.0d0*epsc/GAMMA/(phi**5)*dphidna*dphidnb 488 & +3.0d0*epsc/GAMMA/(phi**4)*d2phidnadnb 489 d2argexpdnb2 = -d2epscdnb2/GAMMA/(phi**3) 490 & +3.0d0*depscdnb/GAMMA/(phi**4)*dphidnb 491 & +3.0d0*depscdnb/GAMMA/(phi**4)*dphidnb 492 & -12.0d0*epsc/GAMMA/(phi**5)*dphidnb**2 493 & +3.0d0*epsc/GAMMA/(phi**4)*d2phidnb2 494#endif 495c 496c A = BETA/GAMMA [exp{-epsilon_c/(GAMMA*phi**3)}-1]**(-1) 497c 498 if (dabs(argexp).lt.EXPTOL) then 499 expinA=dexp(argexp) 500 else 501 expinA=0.0d0 502 endif 503 A = BETA/GAMMA/(expinA-1.0d0) 504 dAdna = -BETA/GAMMA*dargexpdna*expinA/(expinA-1.0d0)**2 505 dAdnb = -BETA/GAMMA*dargexpdnb*expinA/(expinA-1.0d0)**2 506#ifdef SECOND_DERIV 507 d2Adna2 = -BETA/GAMMA*d2argexpdna2 508 & *expinA/(expinA-1.0d0)**2 509 & - BETA/GAMMA*dargexpdna 510 & *dargexpdna*expinA/(expinA-1.0d0)**2 511 & + 2.0d0*BETA/GAMMA*dargexpdna*dargexpdna 512 & *expinA*expinA/(expinA-1.0d0)**3 513 d2Adnadnb = -BETA/GAMMA*d2argexpdnadnb 514 & *expinA/(expinA-1.0d0)**2 515 & - BETA/GAMMA*dargexpdna 516 & *dargexpdnb*expinA/(expinA-1.0d0)**2 517 & + 2.0d0*BETA/GAMMA*dargexpdna*dargexpdnb 518 & *expinA*expinA/(expinA-1.0d0)**3 519 d2Adnb2 = -BETA/GAMMA*d2argexpdnb2 520 & *expinA/(expinA-1.0d0)**2 521 & - BETA/GAMMA*dargexpdnb 522 & *dargexpdnb*expinA/(expinA-1.0d0)**2 523 & + 2.0d0*BETA/GAMMA*dargexpdnb*dargexpdnb 524 & *expinA*expinA/(expinA-1.0d0)**3 525#endif 526c 527c fAt = (1 + A * t**2)/(1 + A * t**2 + A**2 * t**4) and its derivs 528c 529 fAtnum = 1.0d0 530 fAtden = 1.0d0+A*t**2 531 fAt = fAtnum/fAtden 532 dfAtdendt = 2.0d0*A*t 533 dfAtdendA = t**2 534 dfAtdt = (-fAtnum*dfAtdendt)/(fAtden**2) 535 dfAtdA = (-fAtnum*dfAtdendA)/(fAtden**2) 536 dfAtdna = dfAtdt * dtdna + dfAtdA * dAdna 537 dfAtdnb = dfAtdt * dtdnb + dfAtdA * dAdnb 538#ifdef SECOND_DERIV 539 d2fAtdendt2 = 2.0d0*A 540 d2fAtdendtdA = 2.0d0*t 541 d2fAtdendA2 = 0.0d0 542 d2fAtdt2 = (-fAtnum*d2fAtdendt2) 543 & /(fAtden**2) 544 & -2.0d0*(-fAtnum*dfAtdendt) 545 & /(fAtden**3)*dfAtdendt 546 d2fAtdtdA = (-fAtnum*d2fAtdendtdA) 547 & /(fAtden**2) 548 & -2.0d0*(-fAtnum*dfAtdendt) 549 & /(fAtden**3)*dfAtdendA 550 d2fAtdA2 = (-fAtnum*d2fAtdendA2) 551 & /(fAtden**2) 552 & -2.0d0*(-fAtnum*dfAtdendA) 553 & /(fAtden**3)*dfAtdendA 554 d2fAtdna2 = d2fAtdt2*dtdna*dtdna + d2fAtdtdA*dtdna*dAdna 555 & + dfAtdt*d2tdna2 + d2fAtdtdA*dAdna*dtdna 556 & + d2fAtdA2*dAdna*dAdna + dfAtdA*d2Adna2 557 d2fAtdnb2 = d2fAtdt2*dtdnb*dtdnb + d2fAtdtdA*dtdnb*dAdnb 558 & + dfAtdt*d2tdnb2 + d2fAtdtdA*dAdnb*dtdnb 559 & + d2fAtdA2*dAdnb*dAdnb + dfAtdA*d2Adnb2 560 d2fAtdnadnb = d2fAtdt2*dtdna*dtdnb + d2fAtdtdA*dtdna*dAdnb 561 & + dfAtdt*d2tdnadnb + d2fAtdtdA*dAdna*dtdnb 562 & + d2fAtdA2*dAdna*dAdnb + dfAtdA*d2Adnadnb 563#endif 564c 565c arglog = 1 + BETA/GAMMA * t**2 * fAt and its derivs 566c 567 arglog = 1.0d0 + BETA/GAMMA*t**2*fAt 568 darglogdna = BETA/GAMMA*(2.0d0*t*dtdna*fAt 569 & +t*t*dfAtdna) 570 darglogdnb = BETA/GAMMA*(2.0d0*t*dtdnb*fAt 571 & +t*t*dfAtdnb) 572#ifdef SECOND_DERIV 573 d2arglogdna2 = BETA/GAMMA*(2.0d0*dtdna*dtdna*fAt 574 & +2.0d0*t*d2tdna2*fAt 575 & +2.0d0*t*dtdna*dfAtdna 576 & +2.0d0*t*dtdna*dfAtdna 577 & +t*t*d2fAtdna2) 578 d2arglogdnb2 = BETA/GAMMA*(2.0d0*dtdnb*dtdnb*fAt 579 & +2.0d0*t*d2tdnb2*fAt 580 & +2.0d0*t*dtdnb*dfAtdnb 581 & +2.0d0*t*dtdnb*dfAtdnb 582 & +t*t*d2fAtdnb2) 583 d2arglogdnadnb = BETA/GAMMA*(2.0d0*dtdna*dtdnb*fAt 584 & +2.0d0*t*d2tdnadnb*fAt 585 & +2.0d0*t*dtdna*dfAtdnb 586 & +2.0d0*t*dtdnb*dfAtdna 587 & +t*t*d2fAtdnadnb) 588#endif 589c 590c H0 = GAMMA * phi**3 * log{arglog} and its derivs 591c 592 H0 = GAMMA*(phi**3)*dlog(arglog) 593 dH0dna = GAMMA*(3.0d0*(phi**2)*dphidna*dlog(arglog) 594 & +(phi**3)*darglogdna/arglog) 595 dH0dnb = GAMMA*(3.0d0*(phi**2)*dphidnb*dlog(arglog) 596 & +(phi**3)*darglogdnb/arglog) 597#ifdef SECOND_DERIV 598 d2H0dna2 = GAMMA*(6.0d0*phi*dphidna*dphidna*dlog(arglog) 599 & +3.0d0*(phi**2)*d2phidna2*dlog(arglog) 600 & +6.0d0*(phi**2)*dphidna*darglogdna/arglog 601 & +(phi**3)*d2arglogdna2/arglog 602 & -(phi**3)*darglogdna*darglogdna/arglog/arglog) 603 d2H0dnadnb = GAMMA*(6.0d0*phi*dphidna*dphidnb*dlog(arglog) 604 & +3.0d0*(phi**2)*d2phidnadnb*dlog(arglog) 605 & +3.0d0*(phi**2)*dphidna*darglogdnb/arglog 606 & +3.0d0*(phi**2)*dphidnb*darglogdna/arglog 607 & +(phi**3)*d2arglogdnadnb/arglog 608 & -(phi**3)*darglogdna*darglogdnb/arglog/arglog) 609 d2H0dnb2 = GAMMA*(6.0d0*phi*dphidnb*dphidnb*dlog(arglog) 610 & +3.0d0*(phi**2)*d2phidnb2*dlog(arglog) 611 & +6.0d0*(phi**2)*dphidnb*darglogdnb/arglog 612 & +(phi**3)*d2arglogdnb2/arglog 613 & -(phi**3)*darglogdnb*darglogdnb/arglog/arglog) 614#endif 615c 616c Now we update Ec, Amat, and Amat2 617c 618 if (lfac) then 619 Ec = Ec+nepsc*qwght(n)*fac 620 if (ldew) ffunc(n)=ffunc(n)+nepsc*fac 621 endif 622 if (nlfac) then 623 Ec = Ec+(H0*rhoval)*qwght(n)*fac 624 if (ldew) ffunc(n)=ffunc(n)+(H0*rhoval)*fac 625 endif 626 if (lfac) then 627 Amat(n,1) = Amat(n,1) + dnepscdn(1)*fac 628 if (ipol.eq.2) Amat(n,2) = Amat(n,2) + dnepscdn(2)*fac 629#ifdef SECOND_DERIV 630 Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) 631 & + d2nepscdn2(D2_RA_RA)*fac 632 Amat2(n,D2_RA_RB) = Amat2(n,D2_RA_RB) 633 & + d2nepscdn2(D2_RA_RB)*fac 634 if (ipol.eq.2) 635 & Amat2(n,D2_RB_RB) = Amat2(n,D2_RB_RB) 636 & + d2nepscdn2(D2_RB_RB)*fac 637#endif 638 endif 639 if (nlfac)then 640 Amat(n,1) = Amat(n,1) + (H0 + 641 & rhoval*dH0dna)*fac 642 if (ipol.eq.2) Amat(n,2) = Amat(n,2) + (H0 + 643 & rhoval*dH0dnb)*fac 644#ifdef SECOND_DERIV 645 Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) 646 & + (2.d0*dH0dna + rhoval*d2H0dna2)*fac 647 Amat2(n,D2_RA_RB) = Amat2(n,D2_RA_RB) 648 & + (dH0dna + dH0dnb + rhoval*d2H0dnadnb)*fac 649 if (ipol.eq.2) 650 & Amat2(n,D2_RB_RB) = Amat2(n,D2_RB_RB) 651 & + (2.d0*dH0dnb + rhoval*d2H0dnb2)*fac 652#endif 653 endif 654c 655c Now we go into gradient-correction parts 656c Note that the functional depends on |Nabla n| through "t" only 657c 658 if (dsqgamma.gt.TOLL)then 659 dtdg = 0.25d0/(phi*ks*rhoval)/dsqgamma 660 dfAtdg = dfAtdt*dtdg 661 darglogdg = BETA/GAMMA*(2.0d0*t*dtdg*fAt+t*t*dfAtdg) 662 dH0dg = GAMMA*(phi**3)*darglogdg/arglog 663 if (ipol.eq.1) then 664 Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + dH0dg*rhoval*fac 665 Cmat(n,D1_GAB) = Cmat(n,D1_GAB) + dH0dg*rhoval*fac*2.0d0 666 else 667 Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + dH0dg*rhoval*fac 668 Cmat(n,D1_GAB) = Cmat(n,D1_GAB) + dH0dg*rhoval*fac*2.0d0 669 Cmat(n,D1_GBB) = Cmat(n,D1_GBB) + dH0dg*rhoval*fac 670 endif 671#ifdef SECOND_DERIV 672 d2tdg2 = -0.125d0/(phi*ks*rhoval)/(dsqgamma**3) 673 d2tdnadg = -dtdg/rhoval-dtdg/phi*dphidna 674 & -dtdg/ks*dksdna 675 d2tdnbdg = -dtdg/rhoval-dtdg/phi*dphidnb 676 & -dtdg/ks*dksdnb 677 d2fAtdg2 = d2fAtdt2*(dtdg**2)+dfAtdt*d2tdg2 678 d2fAtdnadg = d2fAtdt2*dtdg*dtdna 679 & +d2fAtdtdA*dtdg*dAdna 680 & +dfAtdt*d2tdnadg 681 d2fAtdnbdg = d2fAtdt2*dtdg*dtdnb 682 & +d2fAtdtdA*dtdg*dAdnb 683 & +dfAtdt*d2tdnbdg 684 d2arglogdnadg = BETA/GAMMA*(2.0d0*dtdna*dtdg*fAt 685 & +2.0d0*t*d2tdnadg*fAt 686 & +2.0d0*t*dtdg*dfAtdna 687 & +2.0d0*t*dtdna*dfAtdg 688 & +t*t*d2fAtdnadg) 689 d2arglogdnbdg = BETA/GAMMA*(2.0d0*dtdnb*dtdg*fAt 690 & +2.0d0*t*d2tdnbdg*fAt 691 & +2.0d0*t*dtdg*dfAtdnb 692 & +2.0d0*t*dtdnb*dfAtdg 693 & +t*t*d2fAtdnbdg) 694 d2arglogdg2 = BETA/GAMMA*(2.0d0*dtdg*dtdg*fAt 695 & +2.0d0*t*d2tdg2*fAt 696 & +2.0d0*t*dtdg*dfAtdg 697 & +2.0d0*t*dtdg*dfAtdg 698 & +t*t*d2fAtdg2) 699 d2H0dnadg = GAMMA*3.0d0*phi**2*dphidna*darglogdg/arglog 700 & + GAMMA*phi**3*d2arglogdnadg/arglog 701 & - GAMMA*phi**3*darglogdg*darglogdna/arglog**2 702 d2H0dnbdg = GAMMA*3.0d0*phi**2*dphidnb*darglogdg/arglog 703 & + GAMMA*phi**3*d2arglogdnbdg/arglog 704 & - GAMMA*phi**3*darglogdg*darglogdnb/arglog**2 705 d2H0dg2 = GAMMA*phi**3*d2arglogdg2/arglog 706 & - GAMMA*phi**3*darglogdg*darglogdg/arglog**2 707 Cmat2(n,D2_RA_GAA) = Cmat2(n,D2_RA_GAA) 708 & + (dH0dg + d2H0dnadg*rhoval)*fac 709 Cmat2(n,D2_RA_GAB) = Cmat2(n,D2_RA_GAB) 710 & + 2.0d0*(dH0dg + d2H0dnadg*rhoval)*fac 711 Cmat2(n,D2_RA_GBB) = Cmat2(n,D2_RA_GBB) 712 & + (dH0dg + d2H0dnadg*rhoval)*fac 713 Cmat2(n,D2_RB_GAA) = Cmat2(n,D2_RB_GAA) 714 & + (dH0dg + d2H0dnbdg*rhoval)*fac 715 Cmat2(n,D2_RB_GAB) = Cmat2(n,D2_RB_GAB) 716 & + 2.0d0*(dH0dg + d2H0dnbdg*rhoval)*fac 717 Cmat2(n,D2_RB_GBB) = Cmat2(n,D2_RB_GBB) 718 & + (dH0dg + d2H0dnbdg*rhoval)*fac 719 Cmat2(n,D2_GAA_GAA) = Cmat2(n,D2_GAA_GAA) 720 & + d2H0dg2*rhoval*fac 721 Cmat2(n,D2_GAA_GAB) = Cmat2(n,D2_GAA_GAB) 722 & + 2.0d0*d2H0dg2*rhoval*fac 723 Cmat2(n,D2_GAA_GBB) = Cmat2(n,D2_GAA_GBB) 724 & + d2H0dg2*rhoval*fac 725 Cmat2(n,D2_GAB_GAB) = Cmat2(n,D2_GAB_GAB) 726 & + 4.0d0*d2H0dg2*rhoval*fac 727 Cmat2(n,D2_GAB_GBB) = Cmat2(n,D2_GAB_GBB) 728 & + 2.0d0*d2H0dg2*rhoval*fac 729 Cmat2(n,D2_GBB_GBB) = Cmat2(n,D2_GBB_GBB) 730 & + d2H0dg2*rhoval*fac 731#endif 732 endif 733 20 continue 734c 735 return 736 end 737c 738#ifndef SECOND_DERIV 739#define SECOND_DERIV 740c 741c Compile source again for the 2nd derivative case 742c 743#include "xc_spbe96.F" 744#endif 745