1#ifndef SECOND_DERIV 2 Subroutine xc_perdew91(tol_rho, cfac, lcfac, nlcfac, rho, delrho, 3 & Amat, Cmat, nq, ipol, Ec, qwght, 4 & ldew,ffunc) 5#else 6 Subroutine xc_perdew91_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 Input and other parameters 17c 18 integer lnq 19 integer ipol, nq 20 double precision dummy(1) 21 double precision cfac(*) 22 logical lcfac(*), nlcfac(*), ldew 23 logical lfac, nlfac, lfacl, nlfacl 24 double precision ffunc(*) 25 double precision fac, facl 26 double precision lqwght 27 double precision tol_rho 28c 29c Constants in PBE functional 30c 31 double precision ALPHA, BETA, NU, CCZERO, CX, PI 32 double precision cnoname, ca, cb, cc, cd 33 parameter (ALPHA = 0.09d0) 34 parameter (CCZERO = 0.004235d0) 35 parameter (CX = -0.001667d0) 36 parameter (PI = 3.1415926535897932385d0) 37 parameter (cnoname = 2.568d0) 38 parameter (ca = 23.266d0) 39 parameter (cb = 7.389d-3) 40 parameter (cc = 8.723d0) 41 parameter (cd = 0.472d0) 42c 43c Threshold parameters 44c 45 double precision TOLL, EXPTOL 46 double precision EPS 47 parameter (TOLL = 1.0D-40, EXPTOL = 40.0d0) 48 parameter (EPS = 1.0e-8) 49c 50c Correlation energy 51c 52 double precision Ec 53c 54c Charge Density 55c 56 double precision rho(nq,ipol*(ipol+1)/2) 57 double precision rho_t(3) 58c 59c Charge Density Gradient 60c 61 double precision delrho(nq,3,ipol) 62 double precision dsqgamma 63c 64c Quadrature Weights 65c 66 double precision qwght(nq) 67c 68c Sampling Matrices for the XC Potential 69c 70 double precision Amat(nq,ipol), Cmat(nq,*) 71#ifdef SECOND_DERIV 72c 73c Sampling Matrices for the XC Kernel 74c 75 double precision Amat2(nq,NCOL_AMAT2), Cmat2(nq,NCOL_CMAT2) 76#endif 77c 78c Intermediate derivative results, etc. 79c 80 integer n 81 double precision rhoval, gammaval 82 double precision nepsc, dnepscdn(2) 83 double precision epsc, depscdna, depscdnb 84 double precision H0, dH0dna, dH0dnb, dH0dg 85 double precision phi, dphidna, dphidnb, dphidzeta 86 double precision zeta, dzetadna, dzetadnb 87 double precision arglog, darglogdna, darglogdnb, darglogdg 88 double precision fAt, dfAtdt, dfAtdA 89 double precision fAtnum, dfAtnumdt, dfAtnumdA 90 double precision fAtden, dfAtdendt, dfAtdendA 91 double precision dfAtdna, dfAtdnb, dfAtdg 92 double precision A, dAdna, dAdnb 93 double precision t, dtdna, dtdnb, dtdg 94 double precision ks, dksdna, dksdnb 95 double precision argexp, dargexpdna, dargexpdnb 96 double precision expinA 97 double precision rs, drsdna, drsdnb 98 double precision ccrs, dccrsdna, dccrsdnb 99 double precision cnum, dcnumdna, dcnumdnb 100 double precision cden, dcdendna, dcdendnb 101 double precision kf, dkfdna, dkfdnb 102 double precision H1argexp, dH1argexpdna, dH1argexpdnb 103 double precision dH1argexpdg 104 double precision H1prefac, dH1prefacdna, dH1prefacdnb 105 double precision dH1prefacdg 106 double precision expinH1 107 double precision H1, dH1dna, dH1dnb, dH1dg 108#ifdef SECOND_DERIV 109 double precision d2nepscdn2(NCOL_AMAT2) 110 double precision d2epscdna2, d2epscdnadnb, d2epscdnb2 111 double precision d2H0dna2, d2H0dnadnb, d2H0dnb2 112 double precision d2H0dnadg, d2H0dnbdg, d2H0dg2 113 double precision d2phidzeta2, d2phidna2, d2phidnadnb, d2phidnb2 114 double precision d2zetadna2, d2zetadnadnb, d2zetadnb2 115 double precision d2arglogdna2, d2arglogdnb2, d2arglogdnadnb 116 double precision d2arglogdnadg, d2arglogdnbdg, d2arglogdg2 117 double precision d2fAtdt2, d2fAtdA2, d2fAtdtdA, d2fAtdg2 118 double precision d2fAtnumdt2, d2fAtnumdtdA, d2fAtnumdA2 119 double precision d2fAtdendt2, d2fAtdendtdA, d2fAtdendA2 120 double precision d2fAtdna2, d2fAtdnb2, d2fAtdnadnb 121 double precision d2fAtdnadg, d2fAtdnbdg 122 double precision d2Adna2, d2Adnadnb, d2Adnb2 123 double precision d2tdna2, d2tdnb2, d2tdnadnb 124 double precision d2tdg2, d2tdnadg, d2tdnbdg 125 double precision d2ksdna2, d2ksdnb2, d2ksdnadnb 126 double precision d2argexpdna2, d2argexpdnb2, d2argexpdnadnb 127 double precision d2rsdna2, d2rsdnb2, d2rsdnadnb 128 double precision d2ccrsdna2, d2ccrsdnb2, d2ccrsdnadnb 129 double precision d2cnumdna2, d2cnumdnb2, d2cnumdnadnb 130 double precision d2cdendna2, d2cdendnb2, d2cdendnadnb 131 double precision d2kfdna2, d2kfdnb2, d2kfdnadnb 132 double precision d2H1argexpdna2, d2H1argexpdnadnb, d2H1argexpdnb2 133 double precision d2H1argexpdnadg, d2H1argexpdnbdg, d2H1argexpdg2 134 double precision d2H1prefacdna2, d2H1prefacdnadnb, d2H1prefacdnb2 135 double precision d2H1prefacdnadg, d2H1prefacdnbdg, d2H1prefacdg2 136 double precision d2H1dna2, d2H1dnb2, d2H1dnadnb 137 double precision d2H1dnadg, d2H1dnbdg, d2H1dg2 138#endif 139c 140c References: 141c [a] J. P. Perdew et al., Phys. Rev. B 46, 6671 (1992). 142c [b] M. Rasolt and D. J. W. Geldart, Phys. Rev. B 34, 1325 (1986). 143c 144c E_c(PW91) = Int n (epsilon_c + H0 + H1) dxdydz 145c 146c n*epsilon_c <=== supplied by another subroutine 147c d(n*epsilon_c)/d(na) <=== supplied by another subroutine 148c d2(n*epsilon_c)/d(na)d(na) <=== supplied by another subroutine 149c d2(n*epsilon_c)/d(na)d(nb) <=== supplied by another subroutine 150c d2(n*epsilon_c)/d(nb)d(nb) <=== supplied by another subroutine 151c 152c H0 = BETA**2/2/ALPHA * phi**3 * log{ 1 + 2*ALPHA/BETA * t**2 * [ ... ]} 153c 154c BETA = NU * CCZERO 155c 156c NU = (16/PI)(3*PI**2)**(1/3) 157c 158c phi = (1/2)[(1+zeta)**(2/3)+(1-zeta)**(2/3)] 159c 160c zeta = (na - nb)/n 161c 162c [ ... ] = (1 + A * t**2)/(1 + A * t**2 + A**2 * t**4) 163c 164c A = 2*ALPHA/BETA [exp{-2*ALPHA*epsilon_c/(BETA**2*phi**3)}-1]**(-1) 165c 166c t = |Nabla n|/(2*phi*ks*n) 167c 168c ks = 2 * (3 * PI**2 * n)**(1/6) / sqrt(PI) 169c 170c |Nabla n| = sqrt(g_aa + g_bb + 2*g_ab) 171c 172c H1 = NU * [CC(rs) - CCZERO - 3/7*CX] * Phi**3 * t**2 173c * exp[-100 * phi**4 * (ks**2/kf**2) * t**2] 174c 175c CC(rs) = 0.001 * (cnoname + ca*rs + cb*rs**2)/(1+cc*rs+cd*rs**2+10*cb*rs**3) 176c - CX 177c 178c rs = (3/4 / PI / n)**(1/3) 179c 180c Names of variables 181c 182c E_c(PW91) : Ec 183c n (alpha+beta density) : rhoval 184c na, nb : rho(*,2), rho(*,3) 185c epsilon_c : epsc 186c H0 : H0 187c n*epsilon_c : nepsc 188c phi (also called "g") : phi 189c zeta : zeta 190c { ... } : arglog 191c [ ... ] : fAt 192c (1 + A * t**2) : fAtnum 193c (1 + A * t**2 + A**2 * t**4) : fAtden 194c A : A 195c t : t 196c |Nabla n| : gammaval 197c ks : ks 198c {-epsilon_c ... } : argexp 199c g_aa, g_bb, g_ab : g 200c H1 : H1 201c CC(rs) : ccrs 202c CCZERO : CCZERO 203c CX : CX 204c kf : kf 205c rs : rs 206c 207c Derivatives of these are named like d...dna, d2...dnadnb, 208c d2...dna2, etc. 209c 210 fac = cfac(5) 211 lfac = lcfac(5) 212 nlfac = nlcfac(5) 213c 214 NU = (16.0d0/PI)*(3.0d0*PI**2)**(1.0d0/3.0d0) 215 BETA = NU * CCZERO 216c 217c ======> BOTH SPIN-RESTRICETED AND UNRESTRICTED <====== 218c 219 do 20 n = 1, nq 220c 221c n and zeta = (na - nb)/n 222c 223 rhoval = rho(n,1) 224 rho_t(1) = rho(n,1) 225 if (ipol.eq.2) then 226 rho_t(2) = rho(n,2) 227 rho_t(3) = rho(n,3) 228 endif 229 if (rhoval.le.tol_rho) goto 20 230 if (ipol.eq.1) then 231 gammaval = delrho(n,1,1)*delrho(n,1,1) + 232 & delrho(n,2,1)*delrho(n,2,1) + 233 & delrho(n,3,1)*delrho(n,3,1) 234 else 235 gammaval = delrho(n,1,1)*delrho(n,1,1) + 236 & delrho(n,1,2)*delrho(n,1,2) + 237 & delrho(n,2,1)*delrho(n,2,1) + 238 & delrho(n,2,2)*delrho(n,2,2) + 239 & delrho(n,3,1)*delrho(n,3,1) + 240 & delrho(n,3,2)*delrho(n,3,2) + 241 & 2.d0*(delrho(n,1,1)*delrho(n,1,2) + 242 & delrho(n,2,1)*delrho(n,2,2) + 243 & delrho(n,3,1)*delrho(n,3,2)) 244 endif 245 dsqgamma = max(dsqrt(gammaval),tol_rho) 246 nepsc = 0.0d0 247 dnepscdn(1) = 0.0d0 248 if (ipol.eq.2) dnepscdn(2) = 0.0d0 249#ifdef SECOND_DERIV 250 d2nepscdn2(D2_RA_RA)=0.0d0 251 d2nepscdn2(D2_RA_RB)=0.0d0 252 if (ipol.eq.2) d2nepscdn2(D2_RB_RB)=0.0d0 253#endif 254c 255c call for LDA bit 256c 257 lnq = 1 258 lqwght = 1.0d0 259c 260 if (abs(cfac(1)).gt.EPS) then 261 facl = cfac(1) 262#ifndef SECOND_DERIV 263 call xc_vwn_5(tol_rho,facl,rho_t, 264 & dnepscdn,lnq,ipol,nepsc,lqwght,.false.,dummy) 265#else 266 call xc_vwn_5_d2(tol_rho,facl,rho_t, 267 & dnepscdn,d2nepscdn2,lnq,ipol,nepsc,lqwght, 268 & .false.,dummy) 269#endif 270 endif 271c 272 if (abs(cfac(3)).gt.EPS) then 273 facl = cfac(3) 274 lfacl = lcfac(3) 275 nlfacl = nlcfac(3) 276#ifndef SECOND_DERIV 277 call xc_p81(tol_rho,facl,lfacl,nlfacl,rho_t,dnepscdn, 278 & lnq,ipol,nepsc,lqwght,.false.,dummy) 279#else 280 call xc_p81_d2(tol_rho,facl,lfacl,nlfacl,rho_t, 281 & dnepscdn,d2nepscdn2,lnq,ipol,nepsc,lqwght, 282 & .false.,dummy) 283#endif 284 endif 285c 286 if (abs(cfac(6)).gt.EPS.or.lfac) then 287 facl = cfac(6) 288 if (abs(facl).lt.EPS) facl = 1.0d0 289 lfacl = lcfac(6) 290 nlfacl = nlcfac(6) 291#ifndef SECOND_DERIV 292 call xc_pw91lda(tol_rho,facl,lfacl,nlfacl,rho_t, 293 & dnepscdn,lnq,ipol,nepsc,lqwght, 294 & .false.,dummy) 295#else 296 call xc_pw91lda_d2(tol_rho,facl,lfacl,nlfacl,rho_t, 297 & dnepscdn,d2nepscdn2,lnq,ipol,nepsc,lqwght, 298 & .false.,dummy) 299#endif 300 endif 301c 302 if (abs(cfac(7)).gt.EPS) then 303 facl = cfac(7) 304#ifndef SECOND_DERIV 305 call xc_vwn_1_rpa(tol_rho,facl,rho_t, 306 & dnepscdn,lnq,ipol,nepsc,lqwght, 307 & .false.,dummy) 308#else 309 call xc_vwn_1_rpa_d2(tol_rho,facl,rho_t, 310 & dnepscdn,d2nepscdn2,lnq,ipol,nepsc,lqwght, 311 & .false.,dummy) 312#endif 313 endif 314c 315 if (abs(cfac(8)).gt.EPS) then 316 facl = cfac(8) 317#ifndef SECOND_DERIV 318 call xc_vwn_1(tol_rho,facl,rho_t, 319 & dnepscdn,lnq,ipol,nepsc,lqwght, 320 & .false.,dummy) 321#else 322 call xc_vwn_1_d2(tol_rho,facl,rho_t, 323 & dnepscdn,d2nepscdn2,lnq,ipol,nepsc,lqwght, 324 & .false.,dummy) 325#endif 326 endif 327c 328 if (abs(cfac(9)).gt.EPS) then 329 facl = cfac(9) 330#ifndef SECOND_DERIV 331 call xc_vwn_2(tol_rho,facl,rho_t, 332 & dnepscdn,lnq,ipol,nepsc,lqwght, 333 & .false.,dummy) 334#else 335 call xc_vwn_2_d2(tol_rho,facl,rho_t, 336 & dnepscdn,d2nepscdn2,lnq,ipol,nepsc,lqwght, 337 & .false.,dummy) 338#endif 339 endif 340c 341 if (abs(cfac(10)).gt.EPS) then 342 facl = cfac(10) 343#ifndef SECOND_DERIV 344 call xc_vwn_3(tol_rho,facl,rho_t, 345 & dnepscdn,lnq,ipol,nepsc,lqwght, 346 & .false.,dummy) 347#else 348 call xc_vwn_3_d2(tol_rho,facl,rho_t, 349 & dnepscdn,d2nepscdn2,lnq,ipol,nepsc,lqwght, 350 & .false.,dummy) 351#endif 352 endif 353c 354 if (abs(cfac(11)).gt.EPS) then 355 facl = cfac(11) 356#ifndef SECOND_DERIV 357 call xc_vwn_4(tol_rho,facl,rho_t, 358 & dnepscdn,lnq,ipol,nepsc,lqwght, 359 & .false.,dummy) 360#else 361 call xc_vwn_4_d2(tol_rho,facl,rho_t, 362 & dnepscdn,d2nepscdn2,lnq,ipol,nepsc,lqwght, 363 & .false.,dummy) 364#endif 365 endif 366c ============ 367c PW91 H0 part 368c ============ 369 if(abs(nepsc).lt.tol_rho*tol_rho) goto 20 370c 371c epsilon_c = n*epsilon_c / n 372c 373 epsc = nepsc/rhoval 374 if (ipol.eq.1) then 375 depscdna = dnepscdn(1)/rhoval-nepsc/(rhoval**2) 376 depscdnb = depscdna 377 else 378 depscdna = dnepscdn(1)/rhoval-nepsc/(rhoval**2) 379 depscdnb = dnepscdn(2)/rhoval-nepsc/(rhoval**2) 380 endif 381#ifdef SECOND_DERIV 382 if (ipol.eq.1) then 383 d2epscdna2 = d2nepscdn2(D2_RA_RA)/rhoval 384 & -dnepscdn(1)/(rhoval**2) 385 & -dnepscdn(1)/(rhoval**2) 386 & +2.0d0*nepsc/(rhoval**3) 387 d2epscdnadnb = d2nepscdn2(D2_RA_RB)/rhoval 388 & -dnepscdn(1)/(rhoval**2) 389 & -dnepscdn(1)/(rhoval**2) 390 & +2.0d0*nepsc/(rhoval**3) 391 d2epscdnb2 = d2epscdna2 392 else 393 d2epscdna2 = d2nepscdn2(D2_RA_RA)/rhoval 394 & -dnepscdn(1)/(rhoval**2) 395 & -dnepscdn(1)/(rhoval**2) 396 & +2.0d0*nepsc/(rhoval**3) 397 d2epscdnadnb = d2nepscdn2(D2_RA_RB)/rhoval 398 & -dnepscdn(1)/(rhoval**2) 399 & -dnepscdn(2)/(rhoval**2) 400 & +2.0d0*nepsc/(rhoval**3) 401 d2epscdnb2 = d2nepscdn2(D2_RB_RB)/rhoval 402 & -dnepscdn(2)/(rhoval**2) 403 & -dnepscdn(2)/(rhoval**2) 404 & +2.0d0*nepsc/(rhoval**3) 405 endif 406#endif 407c 408c ks = 2*(3*PI**2*n)**(1/6)/sqrt(PI) and its derivs 409c 410 ks = 2.0d0*(3.0d0*PI*PI*rhoval)**(1.0d0/6.0d0)/dsqrt(PI) 411 dksdna = (1.0d0/6.0d0)*ks/rhoval 412 dksdnb = dksdna 413#ifdef SECOND_DERIV 414 d2ksdna2 = (1.0d0/6.0d0)*dksdna/rhoval 415 & - (1.0d0/6.0d0)*ks/(rhoval**2) 416 d2ksdnadnb = d2ksdna2 417 d2ksdnb2 = d2ksdna2 418#endif 419c 420c zeta = (na-nb)/n and its derivs 421c 422 if (ipol.eq.1) then 423 zeta = 0.0d0 424 else 425 zeta = (rho(n,2)-rho(n,3))/rhoval 426 endif 427 if(zeta.lt.-1.0d0) zeta=-1.0d0 428 if(zeta.gt. 1.0d0) zeta= 1.0d0 429 if (ipol.eq.1) then 430 dzetadna = 1.0d0/rhoval 431 dzetadnb = -1.0d0/rhoval 432#ifdef SECOND_DERIV 433 d2zetadna2 = -2.0d0/(rhoval**2) 434 d2zetadnadnb = 0.0d0 435 d2zetadnb2 = -2.0d0/(rhoval**2) 436#endif 437 else 438 dzetadna = 2.0d0*rho(n,3)/(rhoval**2) 439 dzetadnb = -2.0d0*rho(n,2)/(rhoval**2) 440#ifdef SECOND_DERIV 441 d2zetadna2 = -4.0d0*rho(n,3)/(rhoval**3) 442 d2zetadnadnb = 2.0d0*(rho(n,2)-rho(n,3))/(rhoval**3) 443 d2zetadnb2 = -4.0d0*rho(n,2)/(rhoval**3) 444#endif 445 endif 446c 447c phi = (1/2)[(1+zeta)**(2/3)+(1-zeta)**(2/3)] and its derivs 448c 449 phi = 0.5d0*((1.0d0+zeta)**(2.0d0/3.0d0) 450 & +(1.0d0-zeta)**(2.0d0/3.0d0)) 451 if ((1.0d0-zeta).lt.tol_rho) then 452 dphidzeta = 0.5d0*(2.0d0/3.0d0)*( 453 & (1.0d0+zeta)**(2.0d0/3.0d0)/(1.0d0+zeta)) 454 else if ((1.0d0+zeta).lt.tol_rho) then 455 dphidzeta = 0.5d0*(2.0d0/3.0d0)*( 456 & -(1.0d0-zeta)**(2.0d0/3.0d0)/(1.0d0-zeta)) 457 else 458 dphidzeta = 0.5d0*(2.0d0/3.0d0)*( 459 & (1.0d0+zeta)**(2.0d0/3.0d0)/(1.0d0+zeta) 460 & -(1.0d0-zeta)**(2.0d0/3.0d0)/(1.0d0-zeta)) 461 endif 462 dphidna = dphidzeta*dzetadna 463 dphidnb = dphidzeta*dzetadnb 464#ifdef SECOND_DERIV 465 if ((1.0d0-zeta).lt.tol_rho) then 466 d2phidzeta2 = -(1.0d0/9.0d0)*( 467 & (1.0d0+zeta)**(2.0d0/3.0d0)/((1.0d0+zeta)**2)) 468 else if ((1.0d0+zeta).lt.tol_rho) then 469 d2phidzeta2 = -(1.0d0/9.0d0)*( 470 & (1.0d0-zeta)**(2.0d0/3.0d0)/((1.0d0-zeta)**2)) 471 else 472 d2phidzeta2 = -(1.0d0/9.0d0)*( 473 & (1.0d0+zeta)**(2.0d0/3.0d0)/((1.0d0+zeta)**2) 474 & +(1.0d0-zeta)**(2.0d0/3.0d0)/((1.0d0-zeta)**2)) 475 endif 476 d2phidna2 = d2phidzeta2*dzetadna*dzetadna 477 & + dphidzeta*d2zetadna2 478 d2phidnadnb = d2phidzeta2*dzetadna*dzetadnb 479 & + dphidzeta*d2zetadnadnb 480 d2phidnb2 = d2phidzeta2*dzetadnb*dzetadnb 481 & + dphidzeta*d2zetadnb2 482#endif 483c 484c t = |Nabla n|/(2*phi*ks*n) and its derivs 485c 486 t = dsqgamma/(2.0d0*phi*ks*rhoval) 487 dtdna = -t/rhoval-t/phi*dphidna-t/ks*dksdna 488 dtdnb = -t/rhoval-t/phi*dphidnb-t/ks*dksdnb 489#ifdef SECOND_DERIV 490 d2tdna2 = - dtdna/rhoval 491 & + t/(rhoval**2) 492 & - dtdna/phi*dphidna 493 & + t/(phi**2)*(dphidna**2) 494 & - t/phi*d2phidna2 495 & - dtdna/ks*dksdna 496 & + t/(ks**2)*(dksdna**2) 497 & - t/ks*d2ksdna2 498 d2tdnadnb = - dtdnb/rhoval 499 & + t/(rhoval**2) 500 & - dtdnb/phi*dphidna 501 & + t/(phi**2)*(dphidna*dphidnb) 502 & - t/phi*d2phidnadnb 503 & - dtdnb/ks*dksdna 504 & + t/(ks**2)*(dksdna*dksdnb) 505 & - t/ks*d2ksdnadnb 506 d2tdnb2 = - dtdna/rhoval 507 & + t/(rhoval**2) 508 & - dtdnb/phi*dphidnb 509 & + t/(phi**2)*(dphidnb**2) 510 & - t/phi*d2phidnb2 511 & - dtdnb/ks*dksdnb 512 & + t/(ks**2)*(dksdnb**2) 513 & - t/ks*d2ksdnb2 514#endif 515c 516c { ... } in A (see below) and its derivs 517c 518 argexp = -2.0d0*ALPHA*epsc/BETA**2/(phi**3) 519 dargexpdna = -2.0d0*ALPHA*depscdna/BETA**2/(phi**3) 520 & +6.0d0*ALPHA*epsc/BETA**2/(phi**4)*dphidna 521 dargexpdnb = -2.0d0*ALPHA*depscdnb/BETA**2/(phi**3) 522 & +6.0d0*ALPHA*epsc/BETA**2/(phi**4)*dphidnb 523#ifdef SECOND_DERIV 524 d2argexpdna2 = -2.0d0*ALPHA*d2epscdna2/BETA**2/(phi**3) 525 & +6.0d0*ALPHA*depscdna/BETA**2/(phi**4)*dphidna 526 & +6.0d0*ALPHA*depscdna/BETA**2/(phi**4)*dphidna 527 & -24.0d0*ALPHA*epsc/BETA**2/(phi**5)*dphidna**2 528 & +6.0d0*ALPHA*epsc/BETA**2/(phi**4)*d2phidna2 529 d2argexpdnadnb = -2.0d0*ALPHA*d2epscdnadnb/BETA**2/(phi**3) 530 & +6.0d0*ALPHA*depscdna/BETA**2/(phi**4)*dphidnb 531 & +6.0d0*ALPHA*depscdnb/BETA**2/(phi**4)*dphidna 532 & -24.0d0*ALPHA*epsc/BETA**2/(phi**5)*dphidna*dphidnb 533 & +6.0d0*ALPHA*epsc/BETA**2/(phi**4)*d2phidnadnb 534 d2argexpdnb2 = -2.0d0*ALPHA*d2epscdnb2/BETA**2/(phi**3) 535 & +6.0d0*ALPHA*depscdnb/BETA**2/(phi**4)*dphidnb 536 & +6.0d0*ALPHA*depscdnb/BETA**2/(phi**4)*dphidnb 537 & -24.0d0*ALPHA*epsc/BETA**2/(phi**5)*dphidnb**2 538 & +6.0d0*ALPHA*epsc/BETA**2/(phi**4)*d2phidnb2 539#endif 540c 541c A = 2*ALPHA/BETA [exp{-2*ALPHA*epsilon_c/(BETA**2*phi**3)}-1]**(-1) 542c 543 if (dabs(argexp).lt.EXPTOL) then 544 expinA=dexp(argexp) 545 else 546 expinA=0.0d0 547 endif 548 A = 2.0d0*ALPHA/BETA/(expinA-1.0d0) 549 dAdna = -2.0d0*ALPHA/BETA*dargexpdna*expinA/(expinA-1.0d0)**2 550 dAdnb = -2.0d0*ALPHA/BETA*dargexpdnb*expinA/(expinA-1.0d0)**2 551#ifdef SECOND_DERIV 552 d2Adna2 = -2.0d0*ALPHA/BETA*d2argexpdna2 553 & *expinA/(expinA-1.0d0)**2 554 & - 2.0d0*ALPHA/BETA*dargexpdna 555 & *dargexpdna*expinA/(expinA-1.0d0)**2 556 & + 4.0d0*ALPHA/BETA*dargexpdna*dargexpdna 557 & *expinA*expinA/(expinA-1.0d0)**3 558 d2Adnadnb = -2.0d0*ALPHA/BETA*d2argexpdnadnb 559 & *expinA/(expinA-1.0d0)**2 560 & - 2.0d0*ALPHA/BETA*dargexpdna 561 & *dargexpdnb*expinA/(expinA-1.0d0)**2 562 & + 4.0d0*ALPHA/BETA*dargexpdna*dargexpdnb 563 & *expinA*expinA/(expinA-1.0d0)**3 564 d2Adnb2 = -2.0d0*ALPHA/BETA*d2argexpdnb2 565 & *expinA/(expinA-1.0d0)**2 566 & - 2.0d0*ALPHA/BETA*dargexpdnb 567 & *dargexpdnb*expinA/(expinA-1.0d0)**2 568 & + 4.0d0*ALPHA/BETA*dargexpdnb*dargexpdnb 569 & *expinA*expinA/(expinA-1.0d0)**3 570#endif 571c 572c fAt = (1 + A * t**2)/(1 + A * t**2 + A**2 * t**4) and its derivs 573c 574 fAtnum = 1.0d0+A*t**2 575 fAtden = 1.0d0+A*t**2+A**2*t**4 576 fAt = fAtnum/fAtden 577 dfAtnumdt = 2.0d0*A*t 578 dfAtnumdA = t**2 579 dfAtdendt = 2.0d0*A*t+4.0d0*A**2*t**3 580 dfAtdendA = t**2+2.0d0*A*t**4 581 dfAtdt = (dfAtnumdt*fAtden-fAtnum*dfAtdendt)/(fAtden**2) 582 dfAtdA = (dfAtnumdA*fAtden-fAtnum*dfAtdendA)/(fAtden**2) 583 dfAtdna = dfAtdt * dtdna + dfAtdA * dAdna 584 dfAtdnb = dfAtdt * dtdnb + dfAtdA * dAdnb 585#ifdef SECOND_DERIV 586 d2fAtnumdt2 = 2.0d0*A 587 d2fAtdendt2 = 2.0d0*A+12.0d0*A**2*t**2 588 d2fAtnumdtdA = 2.0d0*t 589 d2fAtdendtdA = 2.0d0*t+8.0d0*A*t**3 590 d2fAtnumdA2 = 0.0d0 591 d2fAtdendA2 = 2.0d0*t**4 592 d2fAtdt2 = (d2fAtnumdt2*fAtden-fAtnum*d2fAtdendt2) 593 & /(fAtden**2) 594 & -2.0d0*(dfAtnumdt*fAtden-fAtnum*dfAtdendt) 595 & /(fAtden**3)*dfAtdendt 596 d2fAtdtdA = (d2fAtnumdtdA*fAtden+dfAtnumdt*dfAtdendA 597 & -dfAtnumdA*dfAtdendt-fAtnum*d2fAtdendtdA) 598 & /(fAtden**2) 599 & -2.0d0*(dfAtnumdt*fAtden-fAtnum*dfAtdendt) 600 & /(fAtden**3)*dfAtdendA 601 d2fAtdA2 = (d2fAtnumdA2*fAtden-fAtnum*d2fAtdendA2) 602 & /(fAtden**2) 603 & -2.0d0*(dfAtnumdA*fAtden-fAtnum*dfAtdendA) 604 & /(fAtden**3)*dfAtdendA 605 d2fAtdna2 = d2fAtdt2*dtdna*dtdna + d2fAtdtdA*dtdna*dAdna 606 & + dfAtdt*d2tdna2 + d2fAtdtdA*dAdna*dtdna 607 & + d2fAtdA2*dAdna*dAdna + dfAtdA*d2Adna2 608 d2fAtdnb2 = d2fAtdt2*dtdnb*dtdnb + d2fAtdtdA*dtdnb*dAdnb 609 & + dfAtdt*d2tdnb2 + d2fAtdtdA*dAdnb*dtdnb 610 & + d2fAtdA2*dAdnb*dAdnb + dfAtdA*d2Adnb2 611 d2fAtdnadnb = d2fAtdt2*dtdna*dtdnb + d2fAtdtdA*dtdna*dAdnb 612 & + dfAtdt*d2tdnadnb + d2fAtdtdA*dAdna*dtdnb 613 & + d2fAtdA2*dAdna*dAdnb + dfAtdA*d2Adnadnb 614#endif 615c 616c arglog = 1 + 2*ALPHA/BETA * t**2 * fAt and its derivs 617c 618 arglog = 1.0d0 + 2.0d0*ALPHA/BETA*t**2*fAt 619 darglogdna = 2.0d0*ALPHA/BETA*(2.0d0*t*dtdna*fAt 620 & +t*t*dfAtdna) 621 darglogdnb = 2.0d0*ALPHA/BETA*(2.0d0*t*dtdnb*fAt 622 & +t*t*dfAtdnb) 623#ifdef SECOND_DERIV 624 d2arglogdna2 = 2.0d0*ALPHA/BETA*(2.0d0*dtdna*dtdna*fAt 625 & +2.0d0*t*d2tdna2*fAt 626 & +2.0d0*t*dtdna*dfAtdna 627 & +2.0d0*t*dtdna*dfAtdna 628 & +t*t*d2fAtdna2) 629 d2arglogdnb2 = 2.0d0*ALPHA/BETA*(2.0d0*dtdnb*dtdnb*fAt 630 & +2.0d0*t*d2tdnb2*fAt 631 & +2.0d0*t*dtdnb*dfAtdnb 632 & +2.0d0*t*dtdnb*dfAtdnb 633 & +t*t*d2fAtdnb2) 634 d2arglogdnadnb = 2.0d0*ALPHA/BETA*(2.0d0*dtdna*dtdnb*fAt 635 & +2.0d0*t*d2tdnadnb*fAt 636 & +2.0d0*t*dtdna*dfAtdnb 637 & +2.0d0*t*dtdnb*dfAtdna 638 & +t*t*d2fAtdnadnb) 639#endif 640c 641c H0 = BETA**2/2/ALPHA * phi**3 * log{arglog} and its derivs 642c 643 H0 = 0.5d0*BETA**2/ALPHA*(phi**3)*dlog(arglog) 644 dH0dna = 0.5d0*BETA**2/ALPHA*(3.0d0*(phi**2) 645 & *dphidna*dlog(arglog) 646 & +(phi**3)*darglogdna/arglog) 647 dH0dnb = 0.5d0*BETA**2/ALPHA*(3.0d0*(phi**2) 648 & *dphidnb*dlog(arglog) 649 & +(phi**3)*darglogdnb/arglog) 650#ifdef SECOND_DERIV 651 d2H0dna2 = 0.5d0*BETA**2/ALPHA 652 & *(6.0d0*phi*dphidna*dphidna*dlog(arglog) 653 & +3.0d0*(phi**2)*d2phidna2*dlog(arglog) 654 & +6.0d0*(phi**2)*dphidna*darglogdna/arglog 655 & +(phi**3)*d2arglogdna2/arglog 656 & -(phi**3)*darglogdna*darglogdna/arglog/arglog) 657 d2H0dnadnb = 0.5d0*BETA**2/ALPHA 658 & *(6.0d0*phi*dphidna*dphidnb*dlog(arglog) 659 & +3.0d0*(phi**2)*d2phidnadnb*dlog(arglog) 660 & +3.0d0*(phi**2)*dphidna*darglogdnb/arglog 661 & +3.0d0*(phi**2)*dphidnb*darglogdna/arglog 662 & +(phi**3)*d2arglogdnadnb/arglog 663 & -(phi**3)*darglogdna*darglogdnb/arglog/arglog) 664 d2H0dnb2 = 0.5d0*BETA**2/ALPHA 665 & *(6.0d0*phi*dphidnb*dphidnb*dlog(arglog) 666 & +3.0d0*(phi**2)*d2phidnb2*dlog(arglog) 667 & +6.0d0*(phi**2)*dphidnb*darglogdnb/arglog 668 & +(phi**3)*d2arglogdnb2/arglog 669 & -(phi**3)*darglogdnb*darglogdnb/arglog/arglog) 670#endif 671c ============ 672c PW91 H1 part 673c ============ 674c 675c rs = (3/4 / PI / n)**(1/3) and its derivs 676c 677 rs = (0.75d0 / PI / dabs(rhoval))**(1.0d0/3.0d0) 678 drsdna = (-1.0d0/3.0d0)*rs/rhoval 679 drsdnb = drsdna 680#ifdef SECOND_DERIV 681 d2rsdna2 = (-4.0d0/3.0d0)*drsdna/rhoval 682 d2rsdnb2 = d2rsdna2 683 d2rsdnadnb = d2rsdna2 684#endif 685c 686c CC(rs) and its derivs 687c 688 cnum = cnoname + ca*rs + cb*rs**2 689 cden = 1.0d0 + cc*rs + cd*rs**2 + 10.0d0*cb*rs**3 690 ccrs = 1.0d-3*cnum/cden - CX 691 dcnumdna = (ca + 2.0d0*cb*rs)*drsdna 692 dcnumdnb = (ca + 2.0d0*cb*rs)*drsdnb 693 dcdendna = (cc + 2.0d0*cd*rs + 30.0d0*cb*rs**2)*drsdna 694 dcdendnb = (cc + 2.0d0*cd*rs + 30.0d0*cb*rs**2)*drsdnb 695 dccrsdna = 1.0d-3*(dcnumdna*cden-cnum*dcdendna)/cden**2 696 dccrsdnb = 1.0d-3*(dcnumdnb*cden-cnum*dcdendnb)/cden**2 697#ifdef SECOND_DERIV 698 d2cnumdna2 = (ca + 2.0d0*cb*rs)*d2rsdna2 699 & + (2.0d0*cb)*drsdna*drsdna 700 d2cnumdnb2 = (ca + 2.0d0*cb*rs)*d2rsdnb2 701 & + (2.0d0*cb)*drsdnb*drsdnb 702 d2cnumdnadnb = (ca + 2.0d0*cb*rs)*d2rsdnadnb 703 & + (2.0d0*cb)*drsdna*drsdnb 704 d2cdendna2 = (cc + 2.0d0*cd*rs + 30.0d0*cb*rs**2)*d2rsdna2 705 & + (2.0d0*cd+60.0d0*cb*rs)*drsdna*drsdna 706 d2cdendnb2 = (cc + 2.0d0*cd*rs + 30.0d0*cb*rs**2)*d2rsdnb2 707 & + (2.0d0*cd+60.0d0*cb*rs)*drsdnb*drsdnb 708 d2cdendnadnb = (cc + 2.0d0*cd*rs + 30.0d0*cb*rs**2)*d2rsdnadnb 709 & + (2.0d0*cd+60.0d0*cb*rs)*drsdna*drsdnb 710 d2ccrsdna2 = 1.0d-3*(d2cnumdna2*cden-cnum*d2cdendna2)/cden**2 711 & - 1.0d-3*2.0d0*(dcnumdna*cden-cnum*dcdendna) 712 & /cden**3*dcdendna 713 d2ccrsdnb2 = 1.0d-3*(d2cnumdnb2*cden-cnum*d2cdendnb2)/cden**2 714 & - 1.0d-3*2.0d0*(dcnumdnb*cden-cnum*dcdendnb) 715 & /cden**3*dcdendnb 716 d2ccrsdnadnb 717 & = 1.0d-3*(d2cnumdnadnb*cden-cnum*d2cdendnadnb)/cden**2 718 & + 1.0d-3*(dcnumdna*dcdendnb-dcnumdnb*dcdendna)/cden**2 719 & - 1.0d-3*2.0d0*(dcnumdna*cden-cnum*dcdendna) 720 & /cden**3*dcdendnb 721#endif 722c 723c kf = (3 * PI**2 * n)**(1/3) 724c 725 kf = (3.0d0*PI**2*rhoval)**(1.0d0/3.0d0) 726 dkfdna = (1.0d0/3.0d0)*kf/rhoval 727 dkfdnb = dkfdna 728#ifdef SECOND_DERIV 729 d2kfdna2 = (-2.0d0/3.0d0)*dkfdna/rhoval 730 d2kfdnadnb = (-2.0d0/3.0d0)*dkfdna/rhoval 731 d2kfdnb2 = (-2.0d0/3.0d0)*dkfdnb/rhoval 732#endif 733c 734c H1argexp = -100 * phi**4 * (ks**2/kf**2) * t**2 and its derivs 735c 736 H1argexp = -100.0d0*phi**4*(ks**2/kf**2)*t**2 737 dH1argexpdna = 4.0d0*H1argexp/phi*dphidna 738 & +2.0d0*H1argexp/ks*dksdna 739 & -2.0d0*H1argexp/kf*dkfdna 740 & +2.0d0*H1argexp/t*dtdna 741 dH1argexpdnb = 4.0d0*H1argexp/phi*dphidnb 742 & +2.0d0*H1argexp/ks*dksdnb 743 & -2.0d0*H1argexp/kf*dkfdnb 744 & +2.0d0*H1argexp/t*dtdnb 745#ifdef SECOND_DERIV 746 d2H1argexpdna2 = 4.0d0*dH1argexpdna/phi*dphidna 747 & -4.0d0*H1argexp/phi**2*dphidna**2 748 & +4.0d0*H1argexp/phi*d2phidna2 749 & +2.0d0*dH1argexpdna/ks*dksdna 750 & -2.0d0*H1argexp/ks**2*dksdna**2 751 & +2.0d0*H1argexp/ks*d2ksdna2 752 & -2.0d0*dH1argexpdna/kf*dkfdna 753 & +2.0d0*H1argexp/kf**2*dkfdna**2 754 & -2.0d0*H1argexp/kf*d2kfdna2 755 & +2.0d0*dH1argexpdna/t*dtdna 756 & -2.0d0*H1argexp/t**2*dtdna**2 757 & +2.0d0*H1argexp/t*d2tdna2 758 d2H1argexpdnb2 = 4.0d0*dH1argexpdnb/phi*dphidnb 759 & -4.0d0*H1argexp/phi**2*dphidnb**2 760 & +4.0d0*H1argexp/phi*d2phidnb2 761 & +2.0d0*dH1argexpdnb/ks*dksdnb 762 & -2.0d0*H1argexp/ks**2*dksdnb**2 763 & +2.0d0*H1argexp/ks*d2ksdnb2 764 & -2.0d0*dH1argexpdnb/kf*dkfdnb 765 & +2.0d0*H1argexp/kf**2*dkfdnb**2 766 & -2.0d0*H1argexp/kf*d2kfdnb2 767 & +2.0d0*dH1argexpdnb/t*dtdnb 768 & -2.0d0*H1argexp/t**2*dtdnb**2 769 & +2.0d0*H1argexp/t*d2tdnb2 770 d2H1argexpdnadnb = 4.0d0*dH1argexpdna/phi*dphidnb 771 & -4.0d0*H1argexp/phi**2*dphidna*dphidnb 772 & +4.0d0*H1argexp/phi*d2phidnadnb 773 & +2.0d0*dH1argexpdna/ks*dksdnb 774 & -2.0d0*H1argexp/ks**2*dksdna*dksdnb 775 & +2.0d0*H1argexp/ks*d2ksdnadnb 776 & -2.0d0*dH1argexpdna/kf*dkfdnb 777 & +2.0d0*H1argexp/kf**2*dkfdna*dkfdnb 778 & -2.0d0*H1argexp/kf*d2kfdnadnb 779 & +2.0d0*dH1argexpdna/t*dtdnb 780 & -2.0d0*H1argexp/t**2*dtdna*dtdnb 781 & +2.0d0*H1argexp/t*d2tdnadnb 782#endif 783c 784c H1prefac = NU*[CC(rs) - CC(0)-3/7*CX] * g**3 * t**2 785c 786 H1prefac = NU*(ccrs - CCZERO - (3.0d0/7.0d0)*CX) 787 & * phi**3 * t**2 788 dH1prefacdna = NU*dccrsdna * phi**3 * t**2 789 & + 3.0d0*H1prefac/phi*dphidna 790 & + 2.0d0*H1prefac/t*dtdna 791 dH1prefacdnb = NU*dccrsdnb * phi**3 * t**2 792 & + 3.0d0*H1prefac/phi*dphidnb 793 & + 2.0d0*H1prefac/t*dtdnb 794#ifdef SECOND_DERIV 795 d2H1prefacdna2 = NU*d2ccrsdna2 * phi**3 * t**2 796 & + 3.0d0* NU*dccrsdna * phi**2*dphidna * t**2 797 & + 2.0d0* NU*dccrsdna * phi**3 * t*dtdna 798 & + 3.0d0*dH1prefacdna/phi*dphidna 799 & - 3.0d0*H1prefac/phi**2*dphidna**2 800 & + 3.0d0*H1prefac/phi*d2phidna2 801 & + 2.0d0*dH1prefacdna/t*dtdna 802 & - 2.0d0*H1prefac/t**2*dtdna**2 803 & + 2.0d0*H1prefac/t*d2tdna2 804 d2H1prefacdnb2 = NU*d2ccrsdnb2 * phi**3 * t**2 805 & + 3.0d0* NU*dccrsdnb * phi**2*dphidnb * t**2 806 & + 2.0d0* NU*dccrsdnb * phi**3 * t*dtdnb 807 & + 3.0d0*dH1prefacdnb/phi*dphidnb 808 & - 3.0d0*H1prefac/phi**2*dphidnb**2 809 & + 3.0d0*H1prefac/phi*d2phidnb2 810 & + 2.0d0*dH1prefacdnb/t*dtdnb 811 & - 2.0d0*H1prefac/t**2*dtdnb**2 812 & + 2.0d0*H1prefac/t*d2tdnb2 813 d2H1prefacdnadnb = NU*d2ccrsdnadnb * phi**3 * t**2 814 & + 3.0d0* NU*dccrsdna * phi**2*dphidnb * t**2 815 & + 2.0d0* NU*dccrsdna * phi**3 * t*dtdnb 816 & + 3.0d0*dH1prefacdna/phi*dphidnb 817 & - 3.0d0*H1prefac/phi**2*dphidna*dphidnb 818 & + 3.0d0*H1prefac/phi*d2phidnadnb 819 & + 2.0d0*dH1prefacdna/t*dtdnb 820 & - 2.0d0*H1prefac/t**2*dtdna*dtdnb 821 & + 2.0d0*H1prefac/t*d2tdnadnb 822#endif 823c 824c H1 = H1prefac * exp(H1argexp) 825c 826 if (dabs(H1argexp).lt.EXPTOL) then 827 expinH1 = dexp(H1argexp) 828 else 829 expinH1 = 0.0d0 830 endif 831 H1 = H1prefac * expinH1 832 dH1dna = dH1prefacdna * expinH1 833 & + H1prefac * dH1argexpdna * expinH1 834 dH1dnb = dH1prefacdnb * expinH1 835 & + H1prefac * dH1argexpdnb * expinH1 836#ifdef SECOND_DERIV 837 d2H1dna2 = d2H1prefacdna2 * expinH1 838 & + dH1prefacdna * dH1argexpdna * expinH1 839 & + dH1prefacdna * dH1argexpdna * expinH1 840 & + H1prefac * d2H1argexpdna2 * expinH1 841 & + H1prefac * dH1argexpdna * dH1argexpdna * expinH1 842 d2H1dnb2 = d2H1prefacdnb2 * expinH1 843 & + dH1prefacdnb * dH1argexpdnb * expinH1 844 & + dH1prefacdnb * dH1argexpdnb * expinH1 845 & + H1prefac * d2H1argexpdnb2 * expinH1 846 & + H1prefac * dH1argexpdnb * dH1argexpdnb * expinH1 847 d2H1dnadnb = d2H1prefacdnadnb * expinH1 848 & + dH1prefacdna * dH1argexpdnb * expinH1 849 & + dH1prefacdnb * dH1argexpdna * expinH1 850 & + H1prefac * d2H1argexpdnadnb * expinH1 851 & + H1prefac * dH1argexpdna * dH1argexpdnb * expinH1 852#endif 853c 854c Now we update Ec, Amat, and Amat2 855c 856 if (lfac) then 857 Ec = Ec+nepsc*qwght(n)*fac 858 if (ldew) ffunc(n)=ffunc(n)+nepsc*fac 859 endif 860 if (nlfac) then 861 Ec = Ec+(H0*rhoval+H1*rhoval)*qwght(n)*fac 862 if (ldew) ffunc(n)=ffunc(n)+(H0*rhoval+H1*rhoval)*fac 863 endif 864 if (lfac) then 865 Amat(n,1) = Amat(n,1) + dnepscdn(1)*fac 866 if (ipol.eq.2) Amat(n,2) = Amat(n,2) + dnepscdn(2)*fac 867#ifdef SECOND_DERIV 868 Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) 869 & + d2nepscdn2(D2_RA_RA)*fac 870 Amat2(n,D2_RA_RB) = Amat2(n,D2_RA_RB) 871 & + d2nepscdn2(D2_RA_RB)*fac 872 if (ipol.eq.2) 873 & Amat2(n,D2_RB_RB) = Amat2(n,D2_RB_RB) 874 & + d2nepscdn2(D2_RB_RB)*fac 875#endif 876 endif 877 if (nlfac)then 878 Amat(n,1) = Amat(n,1) + (H0 + H1 + 879 & rhoval*dH0dna + rhoval*dH1dna)*fac 880 if (ipol.eq.2) Amat(n,2) = Amat(n,2) + (H0 + H1 + 881 & rhoval*dH0dnb + rhoval*dH1dnb)*fac 882#ifdef SECOND_DERIV 883 Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) 884 & + (2.0d0*dH0dna + rhoval*d2H0dna2 885 & + 2.0d0*dH1dna + rhoval*d2H1dna2)*fac 886 Amat2(n,D2_RA_RB) = Amat2(n,D2_RA_RB) 887 & + (dH0dna + dH0dnb + rhoval*d2H0dnadnb 888 & + dH1dna + dH1dnb + rhoval*d2H1dnadnb)*fac 889 if (ipol.eq.2) 890 & Amat2(n,D2_RB_RB) = Amat2(n,D2_RB_RB) 891 & + (2.0d0*dH0dnb + rhoval*d2H0dnb2 892 & + 2.0d0*dH1dnb + rhoval*d2H1dnb2)*fac 893#endif 894 endif 895c 896c Now we go into gradient-correction parts 897c Note that the functional depends on |Nabla n| through "t" only 898c 899cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 900c goto 20 901cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 902 if (dsqgamma.gt.TOLL)then 903c 904c H0 part 905c 906 dtdg = 0.25d0/(phi*ks*rhoval)/dsqgamma 907 dfAtdg = dfAtdt*dtdg 908 darglogdg = 2.0d0*ALPHA/BETA*(2.0d0*t*dtdg*fAt+t*t*dfAtdg) 909 dH0dg = 0.5d0*BETA**2/ALPHA*(phi**3)*darglogdg/arglog 910c 911c H1 part 912c 913 dH1argexpdg = 2.0d0*H1argexp/t*dtdg 914 dH1prefacdg = 2.0d0*H1prefac/t*dtdg 915 dH1dg = dH1prefacdg * expinH1 916 & + H1prefac * dH1argexpdg * expinH1 917c 918c Now form Cmat 919c 920 if (ipol.eq.1) then 921 Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + dH0dg*rhoval*fac 922 & + dH1dg*rhoval*fac 923 Cmat(n,D1_GAB) = Cmat(n,D1_GAB) + dH0dg*rhoval*fac*2.0d0 924 & + dH1dg*rhoval*fac*2.0d0 925 else 926 Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + dH0dg*rhoval*fac 927 & + dH1dg*rhoval*fac 928 Cmat(n,D1_GAB) = Cmat(n,D1_GAB) + dH0dg*rhoval*fac*2.0d0 929 & + dH1dg*rhoval*fac*2.0d0 930 Cmat(n,D1_GBB) = Cmat(n,D1_GBB) + dH0dg*rhoval*fac 931 & + dH1dg*rhoval*fac 932 endif 933#ifdef SECOND_DERIV 934c 935c H0 part 936c 937 d2tdg2 = -0.125d0/(phi*ks*rhoval)/(dsqgamma**3) 938 d2tdnadg = -dtdg/rhoval-dtdg/phi*dphidna 939 & -dtdg/ks*dksdna 940 d2tdnbdg = -dtdg/rhoval-dtdg/phi*dphidnb 941 & -dtdg/ks*dksdnb 942 d2fAtdg2 = d2fAtdt2*(dtdg**2)+dfAtdt*d2tdg2 943 d2fAtdnadg = d2fAtdt2*dtdg*dtdna 944 & +d2fAtdtdA*dtdg*dAdna 945 & +dfAtdt*d2tdnadg 946 d2fAtdnbdg = d2fAtdt2*dtdg*dtdnb 947 & +d2fAtdtdA*dtdg*dAdnb 948 & +dfAtdt*d2tdnbdg 949 d2arglogdnadg = 2.0d0*ALPHA/BETA*(2.0d0*dtdna*dtdg*fAt 950 & +2.0d0*t*d2tdnadg*fAt 951 & +2.0d0*t*dtdg*dfAtdna 952 & +2.0d0*t*dtdna*dfAtdg 953 & +t*t*d2fAtdnadg) 954 d2arglogdnbdg = 2.0d0*ALPHA/BETA*(2.0d0*dtdnb*dtdg*fAt 955 & +2.0d0*t*d2tdnbdg*fAt 956 & +2.0d0*t*dtdg*dfAtdnb 957 & +2.0d0*t*dtdnb*dfAtdg 958 & +t*t*d2fAtdnbdg) 959 d2arglogdg2 = 2.0d0*ALPHA/BETA*(2.0d0*dtdg*dtdg*fAt 960 & +2.0d0*t*d2tdg2*fAt 961 & +2.0d0*t*dtdg*dfAtdg 962 & +2.0d0*t*dtdg*dfAtdg 963 & +t*t*d2fAtdg2) 964 d2H0dnadg = 0.5d0*BETA**2/ALPHA*3.0d0*phi**2 965 & *dphidna*darglogdg/arglog 966 & + 0.5d0*BETA**2/ALPHA*phi**3 967 & *d2arglogdnadg/arglog 968 & - 0.5d0*BETA**2/ALPHA*phi**3 969 & *darglogdg*darglogdna/arglog**2 970 d2H0dnbdg = 0.5d0*BETA**2/ALPHA*3.0d0 971 & *phi**2*dphidnb*darglogdg/arglog 972 & + 0.5d0*BETA**2/ALPHA*phi**3 973 & *d2arglogdnbdg/arglog 974 & - 0.5d0*BETA**2/ALPHA*phi**3 975 & *darglogdg*darglogdnb/arglog**2 976 d2H0dg2 = 0.5d0*BETA**2/ALPHA*phi**3 977 & *d2arglogdg2/arglog 978 & - 0.5d0*BETA**2/ALPHA*phi**3 979 & *darglogdg*darglogdg/arglog**2 980c 981c H1 part 982c 983 d2H1argexpdnadg = 2.0d0*dH1argexpdna/t*dtdg 984 & - 2.0d0*H1argexp/t**2*dtdg*dtdna 985 & + 2.0d0*H1argexp/t*d2tdnadg 986 d2H1argexpdnbdg = 2.0d0*dH1argexpdnb/t*dtdg 987 & - 2.0d0*H1argexp/t**2*dtdg*dtdnb 988 & + 2.0d0*H1argexp/t*d2tdnbdg 989 d2H1argexpdg2 = 2.0d0*dH1argexpdg/t*dtdg 990 & - 2.0d0*H1argexp/t**2*dtdg*dtdg 991 & + 2.0d0*H1argexp/t*d2tdg2 992 d2H1prefacdnadg = 2.0d0*dH1prefacdna/t*dtdg 993 & - 2.0d0*H1prefac/t**2*dtdna*dtdg 994 & + 2.0d0*H1prefac/t*d2tdnadg 995 d2H1prefacdnbdg = 2.0d0*dH1prefacdnb/t*dtdg 996 & - 2.0d0*H1prefac/t**2*dtdnb*dtdg 997 & + 2.0d0*H1prefac/t*d2tdnbdg 998 d2H1prefacdg2 = 2.0d0*dH1prefacdg/t*dtdg 999 & - 2.0d0*H1prefac/t**2*dtdg*dtdg 1000 & + 2.0d0*H1prefac/t*d2tdg2 1001 d2H1dnadg = d2H1prefacdnadg * expinH1 1002 & + dH1prefacdna * dH1argexpdg * expinH1 1003 & + dH1prefacdg * dH1argexpdna * expinH1 1004 & + H1prefac * d2H1argexpdnadg * expinH1 1005 & + H1prefac * dH1argexpdna * dH1argexpdg * expinH1 1006 d2H1dnbdg = d2H1prefacdnbdg * expinH1 1007 & + dH1prefacdnb * dH1argexpdg * expinH1 1008 & + dH1prefacdg * dH1argexpdnb * expinH1 1009 & + H1prefac * d2H1argexpdnbdg * expinH1 1010 & + H1prefac * dH1argexpdnb * dH1argexpdg * expinH1 1011 d2H1dg2 = d2H1prefacdg2 * expinH1 1012 & + dH1prefacdg * dH1argexpdg * expinH1 1013 & + dH1prefacdg * dH1argexpdg * expinH1 1014 & + H1prefac * d2H1argexpdg2 * expinH1 1015 & + H1prefac * dH1argexpdg * dH1argexpdg * expinH1 1016c 1017c Now form Cmat2 1018c 1019 Cmat2(n,D2_RA_GAA) = Cmat2(n,D2_RA_GAA) 1020 & + (dH0dg + d2H0dnadg*rhoval)*fac 1021 & + (dH1dg + d2H1dnadg*rhoval)*fac 1022 Cmat2(n,D2_RA_GAB) = Cmat2(n,D2_RA_GAB) 1023 & + 2.0d0*(dH0dg + d2H0dnadg*rhoval)*fac 1024 & + 2.0d0*(dH1dg + d2H1dnadg*rhoval)*fac 1025 Cmat2(n,D2_RA_GBB) = Cmat2(n,D2_RA_GBB) 1026 & + (dH0dg + d2H0dnadg*rhoval)*fac 1027 & + (dH1dg + d2H1dnadg*rhoval)*fac 1028 Cmat2(n,D2_RB_GAA) = Cmat2(n,D2_RB_GAA) 1029 & + (dH0dg + d2H0dnbdg*rhoval)*fac 1030 & + (dH1dg + d2H1dnbdg*rhoval)*fac 1031 Cmat2(n,D2_RB_GAB) = Cmat2(n,D2_RB_GAB) 1032 & + 2.0d0*(dH0dg + d2H0dnbdg*rhoval)*fac 1033 & + 2.0d0*(dH1dg + d2H1dnbdg*rhoval)*fac 1034 Cmat2(n,D2_RB_GBB) = Cmat2(n,D2_RB_GBB) 1035 & + (dH0dg + d2H0dnbdg*rhoval)*fac 1036 & + (dH1dg + d2H1dnbdg*rhoval)*fac 1037 Cmat2(n,D2_GAA_GAA) = Cmat2(n,D2_GAA_GAA) 1038 & + d2H0dg2*rhoval*fac 1039 & + d2H1dg2*rhoval*fac 1040 Cmat2(n,D2_GAA_GAB) = Cmat2(n,D2_GAA_GAB) 1041 & + 2.0d0*d2H0dg2*rhoval*fac 1042 & + 2.0d0*d2H1dg2*rhoval*fac 1043 Cmat2(n,D2_GAA_GBB) = Cmat2(n,D2_GAA_GBB) 1044 & + d2H0dg2*rhoval*fac 1045 & + d2H1dg2*rhoval*fac 1046 Cmat2(n,D2_GAB_GAB) = Cmat2(n,D2_GAB_GAB) 1047 & + 4.0d0*d2H0dg2*rhoval*fac 1048 & + 4.0d0*d2H1dg2*rhoval*fac 1049 Cmat2(n,D2_GAB_GBB) = Cmat2(n,D2_GAB_GBB) 1050 & + 2.0d0*d2H0dg2*rhoval*fac 1051 & + 2.0d0*d2H1dg2*rhoval*fac 1052 Cmat2(n,D2_GBB_GBB) = Cmat2(n,D2_GBB_GBB) 1053 & + d2H0dg2*rhoval*fac 1054 & + d2H1dg2*rhoval*fac 1055#endif 1056 endif 1057 20 continue 1058c 1059 return 1060 end 1061c 1062#ifndef SECOND_DERIV 1063#define SECOND_DERIV 1064c 1065c Compile source again for the 2nd derivative case 1066c 1067#include "xc_perdew91.F" 1068#endif 1069