1#if !defined SECOND_DERIV && !defined THIRD_DERIV 2 Subroutine xc_cpbe96(tol_rho, cfac, lcfac, nlcfac, rho, delrho, 3 & Amat, Cmat, nq, ipol, Ec, qwght, 4 & ldew,ffunc) 5#elif defined(SECOND_DERIV) && !defined THIRD_DERIV 6 Subroutine xc_cpbe96_d2(tol_rho, cfac, lcfac, nlcfac, 7 & rho, delrho, Amat, Amat2, Cmat, Cmat2, 8 & nq, ipol, Ec, qwght, ldew, ffunc) 9#else 10 Subroutine xc_cpbe96_d3(tol_rho, cfac, lcfac, nlcfac, 11 & rho, delrho, Amat, Amat2, Amat3, 12 & Cmat, Cmat2, Cmat3, 13 & nq, ipol, Ec, qwght, ldew, ffunc) 14#endif 15c 16 Implicit none 17#include "dft2drv.fh" 18#include "dft3drv.fh" 19c 20c 21c Input and other parameters 22c 23 integer lnq 24 integer ipol, nq, id_cpbe 25 double precision dummy(1) 26 double precision cfac(*) 27 logical lcfac(*), nlcfac(*), ldew 28 logical lfac, nlfac, lfacl, nlfacl 29 double precision ffunc(*) 30 double precision fac, facl 31 double precision lqwght 32 double precision tol_rho 33c 34 double precision epsMS, p14f, p14a, p14b 35 double precision rho13, rho23, rho43, rho53 36 double precision BETA, dBETAdr, d2BETAdr2, d3BETAdr3 37 parameter (epsMS=1.e-8) 38 parameter (p14a=0.1d0) 39 parameter (p14b=0.1778d0) 40c 41c Constants in PBE functional 42c 43 double precision GAMMA, BETAzero, PI 44 parameter (GAMMA = 0.03109069086965489503494086371273d0) 45 parameter (BETAzero = 0.06672455060314922d0) 46 parameter (PI = 3.1415926535897932385d0) 47c 48c Threshold parameters 49c 50 double precision TOLL 51 double precision EXPTOL 52 double precision EPS 53 double precision eps_argexp 54c 55 parameter (TOLL = 1.0D-40) 56 parameter (EXPTOL = 40.0d0) 57 parameter (EPS = 1.0e-8) 58 parameter (eps_argexp = 1.0e-14) 59 60c 61c Correlation energy 62c 63 double precision Ec 64c 65c Charge Density 66c 67 double precision rho(nq,ipol*(ipol+1)/2) 68 double precision rho_t(3) 69c 70c Charge Density Gradient 71c 72 double precision delrho(nq,3,ipol) 73 double precision dsqgamma 74c 75c Quadrature Weights 76c 77 double precision qwght(nq) 78c 79c Sampling Matrices for the XC Potential 80c 81 double precision Amat(nq,ipol), Cmat(nq,3) 82c 83#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 84c 85c Sampling Matrices for the XC Kernel 86c 87 double precision Amat2(nq,NCOL_AMAT2), Cmat2(nq,NCOL_CMAT2) 88#endif 89#ifdef THIRD_DERIV 90c 91c Sampling Matrices for the XC third derivatives 92c 93 double precision Amat3(nq,NCOL_AMAT3), Cmat3(nq,NCOL_CMAT3) 94#endif 95c 96c Intermediate derivative results, etc. 97c 98 integer n 99 double precision rhoval, gammaval 100 double precision nepsc, dnepscdn(2) 101 double precision epsc, depscdna, depscdnb 102 double precision H0, dH0dna, dH0dnb, dH0dg 103 double precision phi, dphidna, dphidnb, dphidzeta 104 double precision zeta, dzetadna, dzetadnb 105 double precision arglog, darglogdna, darglogdnb, darglogdg 106 double precision fAt, dfAtdt, dfAtdA 107 double precision fAtnum, dfAtnumdt, dfAtnumdA 108 double precision fAtden, dfAtdendt, dfAtdendA 109 double precision dfAtdna, dfAtdnb, dfAtdg 110 double precision A, dAdna, dAdnb 111 double precision t, dtdna, dtdnb, dtdg 112 double precision ks, dksdna, dksdnb 113 double precision argexp, dargexpdna, dargexpdnb 114 double precision expinA 115 double precision Mz, dMzdna, dMzdnb 116 double precision Nzrt, dNzrtdna, dNzrtdnb 117 double precision dNzrtdg 118 double precision rrho, rrho2 119 double precision opz, omz 120 double precision rks 121 double precision rphi, rphi2, rphi3, rphi4 122 double precision expmone, expmone2 123 double precision t2, t3, t4 124 double precision fAtden2 125 double precision phi2, phi3 126c 127#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 128 double precision d2nepscdn2(NCOL_AMAT2) 129 double precision d2epscdna2, d2epscdnadnb, d2epscdnb2 130 double precision d2H0dna2, d2H0dnadnb, d2H0dnb2 131 double precision d2H0dnadg, d2H0dnbdg, d2H0dg2 132 double precision d2phidzeta2, d2phidna2, d2phidnadnb, d2phidnb2 133 double precision d2zetadna2, d2zetadnadnb, d2zetadnb2 134 double precision d2arglogdna2, d2arglogdnb2, d2arglogdnadnb 135 double precision d2arglogdnadg, d2arglogdnbdg, d2arglogdg2 136 double precision d2fAtdt2, d2fAtdA2, d2fAtdtdA, d2fAtdg2 137 double precision d2fAtnumdt2, d2fAtnumdtdA, d2fAtnumdA2 138 double precision d2fAtdendt2, d2fAtdendtdA, d2fAtdendA2 139 double precision d2fAtdna2, d2fAtdnb2, d2fAtdnadnb 140 double precision d2fAtdnadg, d2fAtdnbdg 141 double precision d2Adna2, d2Adnadnb, d2Adnb2 142 double precision d2tdna2, d2tdnb2, d2tdnadnb 143 double precision d2tdg2, d2tdnadg, d2tdnbdg 144 double precision d2ksdna2, d2ksdnb2, d2ksdnadnb 145 double precision d2argexpdna2, d2argexpdnb2, d2argexpdnadnb 146 double precision d2Mzdna2, d2Mzdnadnb, d2Mzdnb2 147 double precision d2Nzrtdna2, d2Nzrtdnadnb, d2Nzrtdnb2 148 double precision d2Nzrtdnadg, d2Nzrtdnbdg, d2Nzrtdg2 149 double precision rrho3 150 double precision opz2, omz2 151 double precision rks2 152 double precision rphi5 153 double precision expmone3 154 double precision fAtden3 155#endif 156#ifdef THIRD_DERIV 157 double precision d3nepscdn3(NCOL_AMAT3) 158 double precision d3epscdna3, d3epscdna2dnb, d3epscdnadnb2, 159 1 d3epscdnb3 160 double precision d3ksdna3, d3ksdna2dnb, d3ksdnadnb2, d3ksdnb3 161 double precision d3zetadna3, d3zetadna2dnb, d3zetadnadnb2, 162 1 d3zetadnb3 163 double precision d3phidzeta3, d3phidna3, d3phidna2dnb, 164 1 d3phidnadnb2, d3phidnb3 165 double precision d3tdna3, d3tdna2dnb, d3tdnadnb2, d3tdnb3 166 double precision d3argexpdna3, d3argexpdna2dnb, d3argexpdnadnb2, 167 1 d3argexpdnb3 168 double precision d3Adna3, d3Adna2dnb, d3Adnadnb2, d3Adnb3 169 double precision d3fAtnumdt3, d3fAtdendt3, d3fAtnumdt2dA, 170 1 d3fAtdendt2dA, d3fAtnumdtdA2, d3fAtdendtdA2, 171 2 d3fAtnumdA3, d3fAtdendA3 172 double precision d3fAtdt3, d3fAtdt2dA, d3fAtdtdA2, d3fAtdA3 173 double precision d3fAtdna3, d3fAtdna2dnb, d3fAtdnadnb2, 174 1 d3fAtdnb3 175 double precision d3arglogdna3, d3arglogdna2dnb, d3arglogdnadnb2, 176 1 d3arglogdnb3 177 double precision d3H0dna3, d3H0dna2dnb, d3H0dnadnb2, d3H0dnb3 178 double precision d3tdg3, d3tdna2dg, d3tdnadnbdg, d3tdnbdnadg, 179 1 d3tdnb2dg, d3tdnadg2, d3tdnbdg2 180 double precision d3fAtdg3, d3fAtdna2dg, d3fAtdnadnbdg, 181 1 d3fAtdnbdnadg, 182 2 d3fAtdnb2dg, d3fAtdnadg2, d3fAtdnbdg2 183 double precision d3arglogdg3, d3arglogdna2dg, d3arglogdnadnbdg, 184 1 d3arglogdnbdnadg, 185 2 d3arglogdnb2dg, d3arglogdnadg2, d3arglogdnbdg2 186 double precision d3H0dg3, d3H0dna2dg, d3H0dnadnbdg, d3H0dnbdnadg, 187 1 d3H0dnb2dg, d3H0dnadg2, d3H0dnbdg2 188 double precision d3Mzdna3, d3Mzdna2dnb, d3Mzdnadnb2, d3Mzdnb3 189 double precision d3Nzrtdna3, d3Nzrtdna2dnb, d3Nzrtdnadnb2, 190 1 d3Nzrtdnb3 191 double precision d3Nzrtdg3, d3Nzrtdnadg2, d3Nzrtdnbdg2, 192 1 d3Nzrtdna2dg, d3Nzrtdnadnbdg, d3Nzrtdnb2dg 193 double precision rrho4 194 double precision opz3, omz3 195 double precision rks3 196 double precision rphi6 197 double precision expmone4 198 double precision fAtden4 199#endif 200c 201c References: 202c [a] J. P. Perdew, K. Burke, and M. Ernzerhof, 203c {\it Generalized gradient approximation made simple}, 204c Phys.\ Rev.\ Lett. {\bf 77,} 3865 (1996). 205c [b] J. P. Perdew, K. Burke, and Y. Wang, {\it Real-space cutoff 206c construction of a generalized gradient approximation: The PW91 207c density functional}, submitted to Phys.\ Rev.\ B, Feb. 1996. 208c [c] J. P. Perdew and Y. Wang, Phys.\ Rev.\ B {\bf 45}, 13244 (1992). 209c 210c E_c(PBE) = Int n (epsilon_c + H0) dxdydz 211c 212c n*epsilon_c <=== supplied by another subroutine 213c d(n*epsilon_c)/d(na) <=== supplied by another subroutine 214c d2(n*epsilon_c)/d(na)d(na) <=== supplied by another subroutine 215c d2(n*epsilon_c)/d(na)d(nb) <=== supplied by another subroutine 216c d2(n*epsilon_c)/d(nb)d(nb) <=== supplied by another subroutine 217c 218c H0 = GAMMA * phi**3 * log{ 1 + BETA/GAMMA * t**2 * [ ... ]} 219c 220c phi = (1/2)[(1+zeta)**(2/3)+(1-zeta)**(2/3)] 221c 222c zeta = (na - nb)/n 223c 224c [ ... ] = (1 + A * t**2)/(1 + A * t**2 + A**2 * t**4) 225c 226c A = BETA/GAMMA [exp{-epsilon_c/(GAMMA*phi**3)}-1]**(-1) 227c 228c t = |Nabla n|/(2*phi*ks*n) 229c 230c ks = 2 * (3 * PI**2 * n)**(1/6) / sqrt(PI) 231c 232c |Nabla n| = sqrt(g_aa + g_bb + 2*g_ab) 233c 234c Names of variables 235c 236c E_c(PBE) : Ec 237c n (alpha+beta density) : rhoval 238c na, nb : rho(*,2), rho(*,3) 239c epsilon_c : epsc 240c H0 : H0 241c n*epsilon_c : nepsc 242c phi : phi 243c zeta : zeta 244c { ... } : arglog 245c [ ... ] : fAt 246c (1 + A * t**2) : fAtnum 247c (1 + A * t**2 + A**2 * t**4) : fAtden 248c A : A 249c t : t 250c |Nabla n| : gammaval 251c ks : ks 252c {-epsilon_c ... } : argexp 253c g_aa, g_bb, g_ab : g 254c 255c Derivatives of these are named like d...dna, d2...dnadnb, 256c d2...dna2, etc. 257c 258c 259 p14f=(3d0/(4d0*acos(-1d0)))**(1d0/3d0) ! gfortran error if this is defined as a parameter 260 261 if (abs(cfac(12)).gt.epsMS) then 262c original PBE correlation functional 263 id_cpbe = 1 264 fac = cfac(12) 265 lfac = lcfac(12) 266 nlfac = nlcfac(12) 267 elseif (abs(cfac(46)).gt.epsMS) then 268c simplified PBE correlation functional 269 id_cpbe = 2 270 fac = cfac(46) 271 lfac = lcfac(46) 272 nlfac = nlcfac(46) 273 elseif (abs(cfac(64)).gt.epsMS) then 274c simplified PBE correlation functional 275 id_cpbe = 3 276 fac = cfac(64) 277 lfac = lcfac(64) 278 nlfac = nlcfac(64) 279 endif 280c 281c ======> BOTH SPIN-RESTRICTED AND UNRESTRICTED <====== 282c 283 do 20 n = 1, nq 284c 285c n and zeta = (na - nb)/n 286c 287 rhoval = rho(n,1) 288 rho_t(1) = rho(n,1) 289 if (ipol.eq.2) then 290 rho_t(2) = rho(n,2) 291 rho_t(3) = rho(n,3) 292 endif 293 if (rhoval.le.tol_rho) goto 20 294c Daniel (7-24-12): gammaval is gamma^2 in the exchange part of the 295c calculation. 296 if (ipol.eq.1) then 297 gammaval = delrho(n,1,1)*delrho(n,1,1) + 298 & delrho(n,2,1)*delrho(n,2,1) + 299 & delrho(n,3,1)*delrho(n,3,1) 300 else 301 gammaval = delrho(n,1,1)*delrho(n,1,1) + 302 & delrho(n,1,2)*delrho(n,1,2) + 303 & delrho(n,2,1)*delrho(n,2,1) + 304 & delrho(n,2,2)*delrho(n,2,2) + 305 & delrho(n,3,1)*delrho(n,3,1) + 306 & delrho(n,3,2)*delrho(n,3,2) + 307 & 2.d0*(delrho(n,1,1)*delrho(n,1,2) + 308 & delrho(n,2,1)*delrho(n,2,2) + 309 & delrho(n,3,1)*delrho(n,3,2)) 310 endif 311c Daniel (7-26-12): Is this due to numerical instabilities? 312 dsqgamma = max(dsqrt(gammaval),tol_rho) 313c Daniel (7-24-12): variable for storing the product of the density 314c and the correlation energy per particle, and its derivatives. This 315c is later used for storing the correlation energy from the LDA part 316c of the calculation. 317 nepsc = 0.0d0 318 dnepscdn(1) = 0.0d0 319 if (ipol.eq.2) dnepscdn(2) = 0.0d0 320c Daniel (7-29-12): Probably needed for XC third derivatives also 321#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 322 d2nepscdn2(D2_RA_RA)=0.0d0 323 d2nepscdn2(D2_RA_RB)=0.0d0 324 if (ipol.eq.2) d2nepscdn2(D2_RB_RB)=0.0d0 325#endif 326#ifdef THIRD_DERIV 327 d3nepscdn3(D3_RA_RA_RA)=0.0d0 328 d3nepscdn3(D3_RA_RA_RB)=0.0d0 329 d3nepscdn3(D3_RA_RB_RB)=0.0d0 330 if (ipol.eq.2) d3nepscdn3(D3_RB_RB_RB)=0.0d0 331#endif 332c 333c call for LDA bit 334c 335 lnq = 1 336 lqwght = 1.0d0 337c 338 if (abs(cfac(1)).gt.EPS) then 339 facl = cfac(1) 340#if !defined SECOND_DERIV && !defined THIRD_DERIV 341 call xc_vwn_5(tol_rho,facl,rho_t, 342 & dnepscdn,lnq,ipol,nepsc,lqwght,.false.,dummy) 343#elif defined(SECOND_DERIV) && !defined THIRD_DERIV 344 call xc_vwn_5_d2(tol_rho,facl,rho_t, 345 & dnepscdn,d2nepscdn2,lnq,ipol,nepsc,lqwght, 346 & .false.,dummy) 347#else 348 call xc_vwn_5_d3(tol_rho,facl,rho_t, 349 & dnepscdn,d2nepscdn2,d3nepscdn3,lnq,ipol,nepsc,lqwght, 350 & .false.,dummy) 351#endif 352 endif 353c 354 if (abs(cfac(3)).gt.EPS) then 355 facl = cfac(3) 356 lfacl = lcfac(3) 357 nlfacl = nlcfac(3) 358#if !defined SECOND_DERIV && !defined THIRD_DERIV 359 call xc_p81(tol_rho,facl,lfacl,nlfacl,rho_t,dnepscdn, 360 & lnq,ipol,nepsc,lqwght,.false.,dummy) 361#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV) 362 call xc_p81_d2(tol_rho,facl,lfacl,nlfacl,rho_t, 363 & dnepscdn,d2nepscdn2,lnq,ipol,nepsc,lqwght, 364 & .false.,dummy) 365#else 366 call xc_p81_d3(tol_rho,facl,lfacl,nlfacl,rho_t, 367 & dnepscdn,d2nepscdn2,d3nepscdn3,lnq,ipol,nepsc,lqwght, 368 & .false.,dummy) 369#endif 370 endif 371c 372c Daniel (2-19-13): This is the default choice when PBE96 correlation 373c is used. 374 if (abs(cfac(6)).gt.EPS.or.lfac) then 375 facl = cfac(6) 376 if (abs(facl).lt.EPS) facl = 1.0d0 377 lfacl = lcfac(6) 378 nlfacl = nlcfac(6) 379#if !defined SECOND_DERIV && !defined THIRD_DERIV 380 call xc_pw91lda(tol_rho,facl,lfacl,nlfacl,rho_t, 381 & dnepscdn,lnq,ipol,nepsc,lqwght, 382 & .false.,dummy) 383#elif defined(SECOND_DERIV) && !defined THIRD_DERIV 384 call xc_pw91lda_d2(tol_rho,facl,lfacl,nlfacl,rho_t, 385 & dnepscdn,d2nepscdn2,lnq,ipol,nepsc,lqwght, 386 & .false.,dummy) 387#else 388 call xc_pw91lda_d3(tol_rho,facl,lfacl,nlfacl,rho_t, 389 & dnepscdn,d2nepscdn2,d3nepscdn3,lnq,ipol,nepsc,lqwght, 390 & .false.,dummy) 391#endif 392 endif 393c 394 if (abs(cfac(7)).gt.EPS) then 395 facl = cfac(7) 396#if !defined SECOND_DERIV && !defined THIRD_DERIV 397 call xc_vwn_1_rpa(tol_rho,facl,rho_t, 398 & dnepscdn,lnq,ipol,nepsc,lqwght, 399 & .false.,dummy) 400#elif defined(SECOND_DERIV) && !defined THIRD_DERIV 401 call xc_vwn_1_rpa_d2(tol_rho,facl,rho_t, 402 & dnepscdn,d2nepscdn2,lnq,ipol,nepsc,lqwght, 403 & .false.,dummy) 404#else 405 call xc_vwn_1_rpa_d3(tol_rho,facl,rho_t, 406 & dnepscdn,d2nepscdn2,d3nepscdn3,lnq,ipol,nepsc,lqwght, 407 & .false.,dummy) 408#endif 409 endif 410c 411 if (abs(cfac(8)).gt.EPS) then 412 facl = cfac(8) 413#if !defined SECOND_DERIV && !defined THIRD_DERIV 414 call xc_vwn_1(tol_rho,facl,rho_t, 415 & dnepscdn,lnq,ipol,nepsc,lqwght, 416 & .false.,dummy) 417#elif defined(SECOND_DERIV) && !defined THIRD_DERIV 418 call xc_vwn_1_d2(tol_rho,facl,rho_t, 419 & dnepscdn,d2nepscdn2,lnq,ipol,nepsc,lqwght, 420 & .false.,dummy) 421#else 422 call xc_vwn_1_d3(tol_rho,facl,rho_t, 423 & dnepscdn,d2nepscdn2,d3nepscdn3,lnq,ipol,nepsc,lqwght, 424 & .false.,dummy) 425#endif 426 endif 427c 428 if (abs(cfac(9)).gt.EPS) then 429 facl = cfac(9) 430#if !defined SECOND_DERIV && !defined THIRD_DERIV 431 call xc_vwn_2(tol_rho,facl,rho_t, 432 & dnepscdn,lnq,ipol,nepsc,lqwght, 433 & .false.,dummy) 434#elif defined(SECOND_DERIV) && !defined THIRD_DERIV 435 call xc_vwn_2_d2(tol_rho,facl,rho_t, 436 & dnepscdn,d2nepscdn2,lnq,ipol,nepsc,lqwght, 437 & .false.,dummy) 438#else 439 call xc_vwn_2_d3(tol_rho,facl,rho_t, 440 & dnepscdn,d2nepscdn2,d3nepscdn3,lnq,ipol,nepsc,lqwght, 441 & .false.,dummy) 442#endif 443 endif 444c 445 if (abs(cfac(10)).gt.EPS) then 446 facl = cfac(10) 447#if !defined SECOND_DERIV && !defined THIRD_DERIV 448 call xc_vwn_3(tol_rho,facl,rho_t, 449 & dnepscdn,lnq,ipol,nepsc,lqwght, 450 & .false.,dummy) 451#elif defined(SECOND_DERIV) && !defined THIRD_DERIV 452 call xc_vwn_3_d2(tol_rho,facl,rho_t, 453 & dnepscdn,d2nepscdn2,lnq,ipol,nepsc,lqwght, 454 & .false.,dummy) 455#else 456 call xc_vwn_3_d3(tol_rho,facl,rho_t, 457 & dnepscdn,d2nepscdn2,d3nepscdn3,lnq,ipol,nepsc,lqwght, 458 & .false.,dummy) 459#endif 460 endif 461c 462 if (abs(cfac(11)).gt.EPS) then 463 facl = cfac(11) 464#if !defined SECOND_DERIV && !defined THIRD_DERIV 465 call xc_vwn_4(tol_rho,facl,rho_t, 466 & dnepscdn,lnq,ipol,nepsc,lqwght, 467 & .false.,dummy) 468#elif defined(SECOND_DERIV) && !defined THIRD_DERIV 469 call xc_vwn_4_d2(tol_rho,facl,rho_t, 470 & dnepscdn,d2nepscdn2,lnq,ipol,nepsc,lqwght, 471 & .false.,dummy) 472#else 473 call xc_vwn_4_d3(tol_rho,facl,rho_t, 474 & dnepscdn,d2nepscdn2,d3nepscdn3,lnq,ipol,nepsc,lqwght, 475 & .false.,dummy) 476#endif 477 endif 478c ================== 479c PBE non-local part 480c ================== 481 if(abs(nepsc).lt.tol_rho*tol_rho) goto 20 482c 483c epsilon_c = n*epsilon_c / n 484c 485c Daniel (7-24-12): Regardless of the form of the HEG correlation used 486c above, nepsc = eps*qwght*rhoval*cfac. We need to gather the 487c contributions from different spins. The derivatives given are 488c just the derivative of the epsilon from the local (LDA) part. Since 489c epsilon and its derivatives are multiplied by the density in the 490c routines above, we end up taking derivatives of n*epsilon_c / n. 491 rrho = 1.0d0/rhoval 492 rrho2 = rrho*rrho 493c 494 epsc = nepsc*rrho 495 if (ipol.eq.1) then 496 depscdna = dnepscdn(1)*rrho-nepsc*rrho2 497 depscdnb = depscdna 498 else 499 depscdna = dnepscdn(1)*rrho-nepsc*rrho2 500 depscdnb = dnepscdn(2)*rrho-nepsc*rrho2 501 endif 502c Daniel (7-29-12): Probably needed for XC third derivatives also 503#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 504 rrho3 = rrho2*rrho 505 if (ipol.eq.1) then 506 d2epscdna2 = d2nepscdn2(D2_RA_RA)*rrho 507 & -dnepscdn(1)*rrho2 508 & -dnepscdn(1)*rrho2 509 & +2.0d0*nepsc*rrho3 510 d2epscdnadnb = d2nepscdn2(D2_RA_RB)*rrho 511 & -dnepscdn(1)*rrho2 512 & -dnepscdn(1)*rrho2 513 & +2.0d0*nepsc*rrho3 514 d2epscdnb2 = d2epscdna2 515 else 516 d2epscdna2 = d2nepscdn2(D2_RA_RA)*rrho 517 & -dnepscdn(1)*rrho2 518 & -dnepscdn(1)*rrho2 519 & +2.0d0*nepsc*rrho3 520 d2epscdnadnb = d2nepscdn2(D2_RA_RB)*rrho 521 & -dnepscdn(1)*rrho2 522 & -dnepscdn(2)*rrho2 523 & +2.0d0*nepsc*rrho3 524 d2epscdnb2 = d2nepscdn2(D2_RB_RB)*rrho 525 & -dnepscdn(2)*rrho2 526 & -dnepscdn(2)*rrho2 527 & +2.0d0*nepsc*rrho3 528 endif 529#endif 530c 531#ifdef THIRD_DERIV 532c Subtle point here! We don't actually calculate the first derivative 533c of the LDA correlation energy with respect to the beta spin density 534c for a restricted calculation, so we need to use that na = nb 535 rrho4 = rrho3*rrho 536 if (ipol.eq.1) then 537 d3epscdna3 = d3nepscdn3(D3_RA_RA_RA)*rrho 538 1 - d2nepscdn2(D2_RA_RA)*rrho2 539 2 - 2.0d0*d2nepscdn2(D2_RA_RA)*rrho2 540 3 + 6.0d0*dnepscdn(1)*rrho3 541 4 - 6.0d0*nepsc*rrho4 542 d3epscdna2dnb = d3nepscdn3(D3_RA_RA_RB)*rrho 543 1 - d2nepscdn2(D2_RA_RA)*rrho2 544 2 - 2.0d0*d2nepscdn2(D2_RA_RB)*rrho2 545 3 + 4.0d0*dnepscdn(1)*rrho3 546 4 + 2.0d0*dnepscdn(1)*rrho3 547 5 - 6.0d0*nepsc*rrho4 548 d3epscdnadnb2 = d3epscdna2dnb 549 d3epscdnb3 = d3epscdna3 550 else 551 d3epscdna3 = d3nepscdn3(D3_RA_RA_RA)*rrho 552 1 - 3.0d0*d2nepscdn2(D2_RA_RA)*rrho2 553 2 + 6.0d0*dnepscdn(1)*rrho3 554 3 - 6.0d0*nepsc*rrho4 555 d3epscdna2dnb = d3nepscdn3(D3_RA_RA_RB)*rrho 556 1 - d2nepscdn2(D2_RA_RA)*rrho2 557 2 - 2.0d0*d2nepscdn2(D2_RA_RB)*rrho2 558 3 + 4.0d0*dnepscdn(1)*rrho3 559 4 + 2.0d0*dnepscdn(2)*rrho3 560 5 - 6.0d0*nepsc*rrho4 561 d3epscdnadnb2 = d3nepscdn3(D3_RA_RB_RB)*rrho 562 1 - d2nepscdn2(D2_RB_RB)*rrho2 563 2 - 2.0d0*d2nepscdn2(D2_RA_RB)*rrho2 564 3 + 2.0d0*dnepscdn(1)*rrho3 565 4 + 4.0d0*dnepscdn(2)*rrho3 566 5 - 6.0d0*nepsc*rrho4 567 d3epscdnb3 = d3nepscdn3(D3_RB_RB_RB)*rrho 568 1 - 3.0d0*d2nepscdn2(D2_RB_RB)*rrho2 569 2 + 6.0d0*dnepscdn(2)*rrho3 570 3 - 6.0d0*nepsc*rrho4 571 endif 572#endif 573c 574c ks = 2*(3*PI**2*n)**(1/6)/sqrt(PI) and its derivs 575c 576c Daniel (7-24-12): Thomas-Fermi screening vector 577 ks = 2.0d0*(3.0d0*PI*PI*rhoval)**(1.0d0/6.0d0)/dsqrt(PI) 578 dksdna = (1.0d0/6.0d0)*ks*rrho 579 dksdnb = dksdna 580#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 581 d2ksdna2 = (1.0d0/6.0d0)*dksdna*rrho 582 & - (1.0d0/6.0d0)*ks*rrho2 583 d2ksdnadnb = d2ksdna2 584 d2ksdnb2 = d2ksdna2 585#endif 586#ifdef THIRD_DERIV 587 d3ksdna3 = (1.0d0/6.0d0)*d2ksdna2*rrho 588 1 - (1.0d0/3.0d0)*dksdna*rrho2 589 2 + (1.0d0/3.0d0)*ks*rrho3 590 d3ksdna2dnb = d3ksdna3 591 d3ksdnadnb2 = d3ksdna3 592 d3ksdnb3 = d3ksdna3 593#endif 594c 595c zeta = (na-nb)/n and its derivs 596c 597c Daniel (7-24-12): Spin-polarization 598 if (ipol.eq.1) then 599 zeta = 0.0d0 600 else 601 zeta = (rho(n,2)-rho(n,3))*rrho 602 endif 603 if(zeta.lt.-1.0d0) zeta=-1.0d0 604 if(zeta.gt. 1.0d0) zeta= 1.0d0 605 if (ipol.eq.1) then 606 dzetadna = 1.0d0*rrho 607 dzetadnb = -1.0d0*rrho 608#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 609 d2zetadna2 = -2.0d0*rrho2 610 d2zetadnadnb = 0.0d0 611c Daniel (8-27-12): This should be a positive sign. 612c OLD 613c d2zetadnb2 = -2.0d0/(rhoval**2) 614c OLD 615 d2zetadnb2 = 2.0d0*rrho2 616#endif 617c Daniel (7-29-12): Third derivative stuff 618#ifdef THIRD_DERIV 619 d3zetadna3 = 6.0d0*rrho3 620 d3zetadna2dnb = 2.0d0*rrho3 621 d3zetadnadnb2 = -2.0d0*rrho3 622 d3zetadnb3 = -6.0d0*rrho3 623#endif 624 else 625 dzetadna = 2.0d0*rho(n,3)*rrho2 626 dzetadnb = -2.0d0*rho(n,2)*rrho2 627#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 628 d2zetadna2 = -4.0d0*rho(n,3)*rrho3 629 d2zetadnadnb = 2.0d0*(rho(n,2)-rho(n,3))*rrho3 630c Daniel (8-28-12): This should be a positive sign 631c OLD 632c d2zetadnb2 = -4.0d0*rho(n,2)/(rhoval**3) 633c OLD 634 d2zetadnb2 = 4.0d0*rho(n,2)*rrho3 635#endif 636c Daniel (7-29-12): Third derivative stuff 637#ifdef THIRD_DERIV 638 d3zetadna3 = 12.0d0*rho(n,3)*rrho4 639 d3zetadna2dnb = -4.0d0*(rho(n,2)-2.0d0*rho(n,3))* 640 1 rrho4 641 d3zetadnadnb2 = 4.0d0*(rho(n,2)-2.0d0*rho(n,3))* 642 1 rrho4 643 d3zetadnb3 = -12.0d0*rho(n,2)*rrho4 644#endif 645 endif 646c 647c phi = (1/2)[(1+zeta)**(2/3)+(1-zeta)**(2/3)] and its derivs 648c 649c Daniel (7-24-12): Spin-scaling factor 650 opz = 1.0d0 + zeta 651 omz = 1.0d0 - zeta 652c 653 phi = 0.5d0*(opz**(2.0d0/3.0d0) 654 & +omz**(2.0d0/3.0d0)) 655 if (omz.lt.tol_rho) then 656 dphidzeta = 0.5d0*(2.0d0/3.0d0)*( 657 & opz**(2.0d0/3.0d0)/opz) 658 else if (opz.lt.tol_rho) then 659 dphidzeta = 0.5d0*(2.0d0/3.0d0)*( 660 & -omz**(2.0d0/3.0d0)/omz) 661 else 662 dphidzeta = 0.5d0*(2.0d0/3.0d0)*( 663 & opz**(2.0d0/3.0d0)/opz 664 & -omz**(2.0d0/3.0d0)/omz) 665 endif 666 dphidna = dphidzeta*dzetadna 667 dphidnb = dphidzeta*dzetadnb 668#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 669 opz2 = opz*opz 670 omz2 = omz*omz 671c 672 if (omz.lt.tol_rho) then 673 d2phidzeta2 = -(1.0d0/9.0d0)*( 674 & opz**(2.0d0/3.0d0)/opz2) 675 else if (opz.lt.tol_rho) then 676 d2phidzeta2 = -(1.0d0/9.0d0)*( 677 & omz**(2.0d0/3.0d0)/omz2) 678 else 679 d2phidzeta2 = -(1.0d0/9.0d0)*( 680 & opz**(2.0d0/3.0d0)/opz2 681 & +omz**(2.0d0/3.0d0)/omz2) 682 endif 683 d2phidna2 = d2phidzeta2*dzetadna*dzetadna 684 & + dphidzeta*d2zetadna2 685 d2phidnadnb = d2phidzeta2*dzetadna*dzetadnb 686 & + dphidzeta*d2zetadnadnb 687 d2phidnb2 = d2phidzeta2*dzetadnb*dzetadnb 688 & + dphidzeta*d2zetadnb2 689#endif 690#ifdef THIRD_DERIV 691c Daniel (7-30-12): Testing is done here to prevent numerical problems 692 opz3 = opz2*opz 693 omz3 = omz2*omz 694c 695 if (omz.lt.tol_rho) then 696 d3phidzeta3 = (4.0d0/27.0d0)*( 697 1 opz**(2.0d0/3.0d0)/opz3) 698 else if (opz.lt.tol_rho) then 699 d3phidzeta3 = (4.0d0/27.0d0)*( 700 1 -omz**(2.0d0/3.0d0)/omz3) 701 else 702 d3phidzeta3 = (4.0d0/27.0d0)*( 703 1 opz**(2.0d0/3.0d0)/opz3 704 2 -omz**(2.0d0/3.0d0)/omz3) 705 endif 706 d3phidna3 = d3phidzeta3*dzetadna*dzetadna*dzetadna 707 1 + 3.0d0*d2phidzeta2*d2zetadna2*dzetadna 708 2 + dphidzeta*d3zetadna3 709 d3phidna2dnb = d3phidzeta3*dzetadna*dzetadna*dzetadnb 710 1 + 2.0d0*d2phidzeta2*d2zetadnadnb*dzetadna 711 2 + d2phidzeta2*d2zetadna2*dzetadnb 712 3 + dphidzeta*d3zetadna2dnb 713 d3phidnadnb2 = d3phidzeta3*dzetadna*dzetadnb*dzetadnb 714 1 + 2.0d0*d2phidzeta2*d2zetadnadnb*dzetadnb 715 2 + d2phidzeta2*d2zetadnb2*dzetadna 716 3 + dphidzeta*d3zetadnadnb2 717 d3phidnb3 = d3phidzeta3*dzetadnb*dzetadnb*dzetadnb 718 1 + 3.0d0*d2phidzeta2*d2zetadnb2*dzetadnb 719 2 + dphidzeta*d3zetadnb3 720#endif 721c 722c t = |Nabla n|/(2*phi*ks*n) and its derivs 723c 724 rks = 1.0d0/ks 725 rphi = 1.0d0/phi 726 rphi2 = rphi*rphi 727 rphi3 = rphi2*rphi 728c 729 t = dsqgamma*rphi*rks*rrho/2.0d0 730 dtdna = -t*rrho-t*rphi*dphidna-t*rks*dksdna 731 dtdnb = -t*rrho-t*rphi*dphidnb-t*rks*dksdnb 732#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 733 rks2 = rks*rks 734c 735 d2tdna2 = - dtdna*rrho 736 & + t*rrho2 737 & - dtdna*rphi*dphidna 738 & + t*rphi2*dphidna*dphidna 739 & - t*rphi*d2phidna2 740 & - dtdna*rks*dksdna 741 & + t*rks2*dksdna*dksdna 742 & - t*rks*d2ksdna2 743 d2tdnadnb = - dtdnb*rrho 744 & + t*rrho2 745 & - dtdnb*rphi*dphidna 746 & + t*rphi2*dphidna*dphidnb 747 & - t*rphi*d2phidnadnb 748 & - dtdnb*rks*dksdna 749 & + t*rks2*dksdna*dksdnb 750 & - t*rks*d2ksdnadnb 751c 752 d2tdnb2 = - dtdnb*rrho 753 & + t*rrho2 754 & - dtdnb*rphi*dphidnb 755 & + t*rphi2*dphidnb*dphidnb 756 & - t*rphi*d2phidnb2 757 & - dtdnb*rks*dksdnb 758 & + t*rks2*dksdnb*dksdnb 759 & - t*rks*d2ksdnb2 760#endif 761c 762#ifdef THIRD_DERIV 763 rks3 = rks2*rks 764c 765 d3tdna3 = -2.0d0*t*rphi3*dphidna*dphidna*dphidna 766 1 - t*rks*d3ksdna3 767 2 - 2.0d0*t*rrho3 768 3 + 3.0d0*t*rphi2*d2phidna2*dphidna 769 4 - t*rphi*d3phidna3 770 5 - 2.0d0*t*rks3*dksdna*dksdna*dksdna 771 6 + 3.0d0*t*rks2*d2ksdna2*dksdna 772 7 + (dphidna*dphidna*rphi2 773 8 - d2phidna2*rphi 774 9 + dksdna*dksdna*rks2 775 A - d2ksdna2*rks 776 B + rrho2)*dtdna*2.0d0 777 C - (dphidna*rphi + dksdna*rks 778 D + rrho)*d2tdna2 779 d3tdna2dnb = -2.0d0*t*rphi3*dphidna*dphidna*dphidnb 780 1 - t*rks*d3ksdna2dnb 781 2 - 2.0d0*t*rrho3 782 3 + 2.0d0*t*rphi2*d2phidnadnb*dphidna 783 4 + t*rphi2*d2phidna2*dphidnb 784 5 - t*rphi*d3phidna2dnb 785 6 - 2.0d0*t*rks3*dksdna*dksdna*dksdnb 786 7 + 2.0d0*t*rks2*d2ksdnadnb*dksdna 787 8 + t*rks2*d2ksdna2*dksdnb 788 9 + (dphidna*dphidna*rphi2 789 A - d2phidna2*rphi 790 B + dksdna*dksdna*rks2 791 C - d2ksdna2*rks 792 D + rrho2)*dtdnb 793 E + (dphidna*dphidnb*rphi2 794 F - d2phidnadnb*rphi 795 G + dksdna*dksdnb*rks2 796 H - d2ksdnadnb*rks 797 I + rrho2)*dtdna 798 J - (dphidna*rphi + dksdna*rks 799 K + rrho)*d2tdnadnb 800 d3tdnadnb2 = -2.0d0*t*rphi3*dphidna*dphidnb*dphidnb 801 1 - t*rks*d3ksdnadnb2 802 2 - 2.0d0*t*rrho3 803 3 + 2.0d0*t*rphi2*d2phidnadnb*dphidnb 804 4 + t*rphi2*d2phidnb2*dphidna 805 5 - t*rphi*d3phidnadnb2 806 6 - 2.0d0*t*rks3*dksdna*dksdnb*dksdnb 807 7 + 2.0d0*t*rks2*d2ksdnadnb*dksdnb 808 8 + t*rks2*d2ksdnb2*dksdna 809 9 + (dphidna*dphidnb*rphi2 810 A - d2phidnadnb*rphi 811 B + dksdna*dksdnb*rks2 812 C - d2ksdnadnb*rks 813 D + rrho2)*dtdnb 814 E + (dphidna*dphidnb*rphi2 815 F - d2phidnadnb*rphi 816 G + dksdna*dksdnb*rks2 817 H - d2ksdnadnb*rks 818 I + rrho2)*dtdnb 819 J - (dphidna*rphi + dksdna*rks 820 K + rrho)*d2tdnb2 821 d3tdnb3 = -2.0d0*t*rphi3*dphidnb*dphidnb*dphidnb 822 1 - t*rks*d3ksdnb3 823 2 - 2.0d0*t*rrho3 824 3 + 3.0d0*t*rphi2*d2phidnb2*dphidnb 825 4 - t*rphi*d3phidnb3 826 5 - 2.0d0*t*rks3*dksdnb*dksdnb*dksdnb 827 6 + 3.0d0*t*rks2*d2ksdnb2*dksdnb 828 7 + (dphidnb*dphidnb*rphi2 829 8 - d2phidnb2*rphi 830 9 + dksdnb*dksdnb*rks2 831 A - d2ksdnb2*rks 832 B + rrho2)*dtdnb*2.0d0 833 C - (dphidnb*rphi + dksdnb*rks 834 D + rrho)*d2tdnb2 835#endif 836c 837c { ... } in A (see below) and its derivs 838c 839 rphi4 = rphi3*rphi 840c 841 argexp = -epsc*rphi3/GAMMA 842 dargexpdna = -depscdna*rphi3/GAMMA 843 & +3.0d0*epsc*rphi4*dphidna/GAMMA 844 dargexpdnb = -depscdnb*rphi3/GAMMA 845 & +3.0d0*epsc*rphi4*dphidnb/GAMMA 846#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 847 rphi5 = rphi4*rphi 848 d2argexpdna2 = -d2epscdna2*rphi3/GAMMA 849 & +3.0d0*depscdna*rphi4*dphidna/GAMMA 850 & +3.0d0*depscdna*rphi4*dphidna/GAMMA 851 & -12.0d0*epsc*rphi5*dphidna*dphidna/GAMMA 852 & +3.0d0*epsc*rphi4*d2phidna2/GAMMA 853 d2argexpdnadnb = -d2epscdnadnb*rphi3/GAMMA 854 & +3.0d0*depscdna*rphi4*dphidnb/GAMMA 855 & +3.0d0*depscdnb*rphi4*dphidna/GAMMA 856 & -12.0d0*epsc*rphi5*dphidna*dphidnb/GAMMA 857 & +3.0d0*epsc*rphi4*d2phidnadnb/GAMMA 858 d2argexpdnb2 = -d2epscdnb2*rphi3/GAMMA 859 & +3.0d0*depscdnb*rphi4*dphidnb/GAMMA 860 & +3.0d0*depscdnb*rphi4*dphidnb/GAMMA 861 & -12.0d0*epsc*rphi5*dphidnb*dphidnb/GAMMA 862 & +3.0d0*epsc*rphi4*d2phidnb2/GAMMA 863#endif 864c 865#ifdef THIRD_DERIV 866 rphi6 = rphi5*rphi 867c 868 d3argexpdna3 = -12.0d0*rphi5* 869 1 ( 3.0d0*dphidna*dphidna*depscdna 870 2 + 3.0d0*dphidna*d2phidna2*epsc)/GAMMA 871 3 + 60.0d0*rphi6* 872 4 dphidna*dphidna*dphidna*epsc/GAMMA 873 5 + 3.0d0*rphi4* 874 6 ( 3.0d0*d2phidna2*depscdna 875 7 + 3.0d0*dphidna*d2epscdna2 876 8 + d3phidna3*epsc)/GAMMA 877 9 - d3epscdna3*rphi3/GAMMA 878 d3argexpdna2dnb = -12.0d0*rphi5* 879 1 ( 2.0d0*dphidna*dphidnb*depscdna 880 2 + dphidna*dphidna*depscdnb 881 3 + 2.0d0*dphidna*d2phidnadnb*epsc 882 4 + dphidnb*d2phidna2*epsc)/GAMMA 883 5 + 60.0d0*rphi6* 884 6 dphidna*dphidna*dphidnb*epsc/GAMMA 885 7 + 3.0d0*rphi4* 886 8 ( 2.0d0*d2phidnadnb*depscdna 887 9 + d2phidna2*depscdnb 888 A + 2.0d0*dphidna*d2epscdnadnb 889 B + dphidnb*d2epscdna2 890 C + d3phidna2dnb*epsc)/GAMMA 891 D - d3epscdna2dnb*rphi3/GAMMA 892 d3argexpdnadnb2 = -12.0d0*rphi5* 893 1 ( 2.0d0*dphidna*dphidnb*depscdnb 894 2 + dphidnb*dphidnb*depscdna 895 3 + 2.0d0*dphidnb*d2phidnadnb*epsc 896 4 + dphidna*d2phidnb2*epsc)/GAMMA 897 5 + 60.0d0*rphi6* 898 6 dphidna*dphidnb*dphidnb*epsc/GAMMA 899 7 + 3.0d0*rphi4* 900 8 ( 2.0d0*d2phidnadnb*depscdnb 901 9 + d2phidnb2*depscdna 902 A + 2.0d0*dphidnb*d2epscdnadnb 903 B + dphidna*d2epscdnb2 904 C + d3phidnadnb2*epsc)/GAMMA 905 D - d3epscdnadnb2*rphi3/GAMMA 906 d3argexpdnb3 = -12.0d0*rphi5* 907 1 ( 3.0d0*dphidnb*dphidnb*depscdnb 908 2 + 3.0d0*dphidnb*d2phidnb2*epsc)/GAMMA 909 3 + 60.0d0*rphi6* 910 4 dphidnb*dphidnb*dphidnb*epsc/GAMMA 911 5 + 3.0d0*rphi4* 912 6 ( 3.0d0*d2phidnb2*depscdnb 913 7 + 3.0d0*dphidnb*d2epscdnb2 914 8 + d3phidnb3*epsc)/GAMMA 915 9 - d3epscdnb3*rphi3/GAMMA 916#endif 917c 918c A = BETA/GAMMA [exp{-epsilon_c/(GAMMA*phi**3)}-1]**(-1) 919c 920c Daniel (7-24-12): "A" function 921c this avoids expmone to get a zero value that will result in 1./0. 922 if (dabs(argexp).lt.eps_argexp) then 923 expinA=dexp(eps_argexp) 924 elseif (dabs(argexp).lt.EXPTOL) then 925 expinA=dexp(argexp) 926 else 927 expinA=0.0d0 928 endif 929 expmone = expinA - 1.0d0 930 expmone2 = expmone*expmone 931c 932c revTPSS variant of cPBE functional 933c 934 if (id_cpbe.eq.3) then 935c 936 rho13 = rhoval**(1d0/3d0) 937 rho23 = rho13*rho13 938 rho43 = rho23*rho23 939 rho53 = rho43*rho13 940c 941 BETA = BETAzero*(rho13+p14a*p14f)/(rho13+p14b*p14f) 942 dBETAdr = BETAzero*(1d0/3d0)*(p14b*p14f-p14a*p14f) / 943 & (rho23*(rho13+p14b*p14f)**2) 944 d2BETAdr2 = BETAzero*(-2d0/9d0)*(p14b*p14f-p14a*p14f)* 945 & (2d0/rho43 + p14b*p14f/rho53) / (rho13+p14b*p14f)**3 946#ifdef THIRD_DERIV 947 stop 'not yet implemented 3rd derivative of BETA vs r' 948#endif 949c 950 else 951 BETA = BETAzero 952 dBETAdr = 0.0d0 953 d2BETAdr2 = 0.0d0 954 d3BETAdr3 = 0.0d0 955 endif 956c 957 A = BETA/GAMMA/expmone 958 dAdna = -BETA/GAMMA*dargexpdna*expinA/expmone2 959 dAdnb = -BETA/GAMMA*dargexpdnb*expinA/expmone2 960c 961 if (id_cpbe.eq.3) then 962 dAdna = dAdna + dBETAdr/GAMMA/expmone 963 dAdnb = dAdnb + dBETAdr/GAMMA/expmone 964 endif 965c 966#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 967 expmone3 = expmone2*expmone 968c 969 d2Adna2 = -BETA/GAMMA*d2argexpdna2 970 & *expinA/expmone2 971 & - BETA/GAMMA*dargexpdna 972 & *dargexpdna*expinA/expmone2 973 & + 2.0d0*BETA/GAMMA*dargexpdna*dargexpdna 974 & *expinA*expinA/expmone3 975 d2Adnadnb = -BETA/GAMMA*d2argexpdnadnb 976 & *expinA/expmone2 977 & - BETA/GAMMA*dargexpdna 978 & *dargexpdnb*expinA/expmone2 979 & + 2.0d0*BETA/GAMMA*dargexpdna*dargexpdnb 980 & *expinA*expinA/expmone3 981 d2Adnb2 = -BETA/GAMMA*d2argexpdnb2 982 & *expinA/expmone2 983 & - BETA/GAMMA*dargexpdnb 984 & *dargexpdnb*expinA/expmone2 985 & + 2.0d0*BETA/GAMMA*dargexpdnb*dargexpdnb 986 & *expinA*expinA/expmone3 987c 988 if (id_cpbe.eq.3) then 989 d2Adna2 = d2Adna2 + d2BETAdr2/GAMMA/expmone 990 & - dBETAdr/GAMMA*dargexpdna*expinA/expmone2 991 d2Adnadnb = d2Adnadnb + d2BETAdr2/GAMMA/expmone 992 & - dBETAdr/GAMMA*dargexpdna*expinA/expmone2 993 & - dBETAdr/GAMMA*dargexpdnb*expinA/expmone2 994 d2Adnb2 = d2Adnb2 + d2BETAdr2/GAMMA/expmone 995 & - dBETAdr/GAMMA*dargexpdnb*expinA/expmone2 996 endif 997c 998#endif 999#ifdef THIRD_DERIV 1000 expmone4 = expmone3*expmone 1001c 1002 d3Adna3 = -6.0d0*BETA/GAMMA*dargexpdna*dargexpdna* 1003 1 dargexpdna*expinA*expinA*expinA/ 1004 2 expmone4 1005 3 + 6.0d0*BETA/GAMMA* 1006 4 ( dargexpdna*d2argexpdna2 1007 5 + dargexpdna*dargexpdna*dargexpdna)* 1008 6 expinA*expinA/expmone3 1009 7 - BETA/GAMMA* 1010 8 ( d3argexpdna3 + 3.0d0*d2argexpdna2*dargexpdna 1011 9 + dargexpdna*dargexpdna*dargexpdna)* 1012 A expinA/expmone2 1013 d3Adna2dnb = -6.0d0*BETA/GAMMA*dargexpdna*dargexpdna* 1014 1 dargexpdnb*expinA*expinA*expinA/ 1015 2 expmone4 1016 3 + 2.0d0*BETA/GAMMA* 1017 4 ( 2.0d0*dargexpdna*d2argexpdnadnb 1018 5 + 2.0d0*dargexpdna*dargexpdna*dargexpdnb 1019 6 + dargexpdnb*d2argexpdna2 1020 7 + dargexpdna*dargexpdna*dargexpdnb)* 1021 8 expinA*expinA/expmone3 1022 9 - BETA/GAMMA* 1023 A ( d3argexpdna2dnb 1024 B + 2.0d0*d2argexpdnadnb*dargexpdna 1025 C + d2argexpdna2*dargexpdnb 1026 D + dargexpdna*dargexpdna*dargexpdnb)* 1027 E expinA/expmone2 1028 d3Adnadnb2 = -6.0d0*BETA/GAMMA*dargexpdna*dargexpdnb* 1029 1 dargexpdnb*expinA*expinA*expinA/ 1030 2 expmone4 1031 3 + 2.0d0*BETA/GAMMA* 1032 4 ( 2.0d0*dargexpdnb*d2argexpdnadnb 1033 5 + 2.0d0*dargexpdna*dargexpdnb*dargexpdnb 1034 6 + dargexpdna*d2argexpdnb2 1035 7 + dargexpdna*dargexpdnb*dargexpdnb)* 1036 8 expinA*expinA/expmone3 1037 9 - BETA/GAMMA* 1038 A ( d3argexpdnadnb2 1039 B + 2.0d0*d2argexpdnadnb*dargexpdnb 1040 C + d2argexpdnb2*dargexpdna 1041 D + dargexpdna*dargexpdnb*dargexpdnb)* 1042 E expinA/expmone2 1043 d3Adnb3 = -6.0d0*BETA/GAMMA*dargexpdnb*dargexpdnb* 1044 1 dargexpdnb*expinA*expinA*expinA/ 1045 2 expmone4 1046 3 + 6.0d0*BETA/GAMMA* 1047 4 ( dargexpdnb*d2argexpdnb2 1048 5 + dargexpdnb*dargexpdnb*dargexpdnb)* 1049 6 expinA*expinA/expmone3 1050 7 - BETA/GAMMA* 1051 8 ( d3argexpdnb3 + 3.0d0*d2argexpdnb2*dargexpdnb 1052 9 + dargexpdnb*dargexpdnb*dargexpdnb)* 1053 A expinA/expmone2 1054#endif 1055c 1056c fAt = (1 + A * t**2)/(1 + A * t**2 + A**2 * t**4) and its derivs 1057c 1058 t2 = t*t 1059 t3 = t2*t 1060 t4 = t3*t 1061c 1062 if (id_cpbe.eq.1.or.id_cpbe.eq.3) then 1063 fAtnum = 1.0d0+A*t2 1064 fAtden = 1.0d0+A*t2+A*A*t4 1065 fAtden2 = fAtden*fAtden 1066 fAt = fAtnum/fAtden 1067 dfAtnumdt = 2.0d0*A*t 1068 dfAtnumdA = t2 1069 dfAtdendt = 2.0d0*A*t+4.0d0*A*A*t3 1070 dfAtdendA = t2+2.0d0*A*t4 1071 dfAtdt = (dfAtnumdt*fAtden-fAtnum*dfAtdendt)/fAtden2 1072 dfAtdA = (dfAtnumdA*fAtden-fAtnum*dfAtdendA)/fAtden2 1073 else 1074 fAtnum = 1.0d0 1075 fAtden = 1.0d0+A*t**2 1076 fAt = fAtnum/fAtden 1077 dfAtdendt = 2.0d0*A*t 1078 dfAtdendA = t**2 1079 dfAtdt = (-fAtnum*dfAtdendt)/(fAtden**2) 1080 dfAtdA = (-fAtnum*dfAtdendA)/(fAtden**2) 1081 endif 1082 dfAtdna = dfAtdt * dtdna + dfAtdA * dAdna 1083 dfAtdnb = dfAtdt * dtdnb + dfAtdA * dAdnb 1084#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 1085 fAtden3 = fAtden2*fAtden 1086c 1087 if (id_cpbe.eq.1.or.id_cpbe.eq.3) then 1088 d2fAtnumdt2 = 2.0d0*A 1089 d2fAtdendt2 = 2.0d0*A+12.0d0*A**2*t2 1090 d2fAtnumdtdA = 2.0d0*t 1091 d2fAtdendtdA = 2.0d0*t+8.0d0*A*t3 1092 d2fAtnumdA2 = 0.0d0 1093 d2fAtdendA2 = 2.0d0*t4 1094 d2fAtdt2 = (d2fAtnumdt2*fAtden-fAtnum*d2fAtdendt2) 1095 & /(fAtden2) 1096 & -2.0d0*(dfAtnumdt*fAtden-fAtnum*dfAtdendt) 1097 & /(fAtden3)*dfAtdendt 1098 d2fAtdtdA = (d2fAtnumdtdA*fAtden+dfAtnumdt*dfAtdendA 1099 & -dfAtnumdA*dfAtdendt-fAtnum*d2fAtdendtdA) 1100 & /(fAtden2) 1101 & -2.0d0*(dfAtnumdt*fAtden-fAtnum*dfAtdendt) 1102 & /(fAtden3)*dfAtdendA 1103 d2fAtdA2 = (d2fAtnumdA2*fAtden-fAtnum*d2fAtdendA2) 1104 & /(fAtden2) 1105 & -2.0d0*(dfAtnumdA*fAtden-fAtnum*dfAtdendA) 1106 & /(fAtden3)*dfAtdendA 1107 else 1108 d2fAtdendt2 = 2.0d0*A 1109 d2fAtdendtdA = 2.0d0*t 1110 d2fAtdendA2 = 0.0d0 1111 d2fAtdt2 = (-fAtnum*d2fAtdendt2) 1112 & /(fAtden**2) 1113 & -2.0d0*(-fAtnum*dfAtdendt) 1114 & /(fAtden**3)*dfAtdendt 1115 d2fAtdtdA = (-fAtnum*d2fAtdendtdA) 1116 & /(fAtden**2) 1117 & -2.0d0*(-fAtnum*dfAtdendt) 1118 & /(fAtden**3)*dfAtdendA 1119 d2fAtdA2 = (-fAtnum*d2fAtdendA2) 1120 & /(fAtden**2) 1121 & -2.0d0*(-fAtnum*dfAtdendA) 1122 & /(fAtden**3)*dfAtdendA 1123 endif 1124c 1125 d2fAtdna2 = d2fAtdt2*dtdna*dtdna + d2fAtdtdA*dtdna*dAdna 1126 & + dfAtdt*d2tdna2 + d2fAtdtdA*dAdna*dtdna 1127 & + d2fAtdA2*dAdna*dAdna + dfAtdA*d2Adna2 1128 d2fAtdnb2 = d2fAtdt2*dtdnb*dtdnb + d2fAtdtdA*dtdnb*dAdnb 1129 & + dfAtdt*d2tdnb2 + d2fAtdtdA*dAdnb*dtdnb 1130 & + d2fAtdA2*dAdnb*dAdnb + dfAtdA*d2Adnb2 1131 d2fAtdnadnb = d2fAtdt2*dtdna*dtdnb + d2fAtdtdA*dtdna*dAdnb 1132 & + dfAtdt*d2tdnadnb + d2fAtdtdA*dAdna*dtdnb 1133 & + d2fAtdA2*dAdna*dAdnb + dfAtdA*d2Adnadnb 1134#endif 1135#ifdef THIRD_DERIV 1136 fAtden4 = fAtden3*fAtden 1137c 1138 d3fAtnumdt3 = 0.0d0 1139 d3fAtdendt3 = 24.0d0*A**2*t 1140 d3fAtnumdt2dA = 2.0d0 1141 d3fAtdendt2dA = 2.0d0+24.0d0*A*t2 1142 d3fAtnumdtdA2 = 0.0d0 1143 d3fAtdendtdA2 = 8.0d0*t3 1144 d3fAtnumdA3 = 0.0d0 1145 d3fAtdendA3 = 0.0d0 1146c Here we pay attention to which derivatives are zero. 1147 d3fAtdt3 = -( fAtnum*d3fAtdendt3 1148 1 + 3.0d0*dfAtnumdt*d2fAtdendt2 1149 2 + 3.0d0*d2fAtnumdt2*dfAtdendt)/ 1150 3 (fAtden2) 1151 4 + 6.0d0*( fAtnum*d2fAtdendt2*dfAtdendt 1152 5 + dfAtnumdt*dfAtdendt*dfAtdendt)/ 1153 6 (fAtden3) 1154 7 - 6.0d0*fAtnum*dfAtdendt*dfAtdendt*dfAtdendt/ 1155 8 (fAtden4) 1156 d3fAtdt2dA = -( fAtnum*d3fAtdendt2dA 1157 1 + 2.0d0*dfAtnumdt*d2fAtdendtdA 1158 2 + dfAtnumdA*d2fAtdendt2 1159 3 + 2.0d0*d2fAtnumdtdA*dfAtdendt 1160 4 + d2fAtnumdt2*dfAtdendA)/ 1161 5 (fAtden2) 1162 6 + 2.0d0*( 2.0d0*fAtnum*d2fAtdendtdA*dfAtdendt 1163 7 + fAtnum*d2fAtdendt2*dfAtdendA 1164 8 + 2.0d0*dfAtnumdt*dfAtdendA*dfAtdendt 1165 9 + dfAtnumdA*dfAtdendt*dfAtdendt)/ 1166 A (fAtden3) 1167 B - 6.0d0*fAtnum*dfAtdendt*dfAtdendt*dfAtdendA/ 1168 C (fAtden4) 1169 D + d3fAtnumdt2dA/fAtden 1170 d3fAtdtdA2 = -( fAtnum*d3fAtdendtdA2 1171 1 + 2.0d0*dfAtnumdA*d2fAtdendtdA 1172 2 + dfAtnumdt*d2fAtdendA2 1173 3 + 2.0d0*d2fAtnumdtdA*dfAtdendA)/ 1174 4 (fAtden2) 1175 5 + 2.0d0*( 2.0d0*fAtnum*d2fAtdendtdA*dfAtdendA 1176 6 + fAtnum*d2fAtdendA2*dfAtdendt 1177 7 + 2.0d0*dfAtnumdA*dfAtdendA*dfAtdendt 1178 8 + dfAtnumdt*dfAtdendA*dfAtdendA)/ 1179 9 (fAtden3) 1180 A - 6.0d0*fAtnum*dfAtdendt*dfAtdendA*dfAtdendA/ 1181 B (fAtden4) 1182 d3fAtdA3 = -3.0d0*dfAtnumdA*d2fAtdendA2/ 1183 1 (fAtden2) 1184 2 + 6.0d0*( fAtnum*d2fAtdendA2*dfAtdendA 1185 3 + dfAtnumdA*dfAtdendA*dfAtdendA)/ 1186 3 (fAtden3) 1187 4 - 6.0d0*fAtnum*dfAtdendA*dfAtdendA*dfAtdendA/ 1188 5 (fAtden4) 1189c 1190 d3fAtdna3 = d3fAtdA3*dAdna*dAdna*dAdna 1191 1 + 3.0d0*d3fAtdtdA2*dAdna*dAdna*dtdna 1192 2 + 3.0d0*d2fAtdA2*d2Adna2*dAdna 1193 3 + 3.0d0*d3fAtdt2dA*dAdna*dtdna*dtdna 1194 4 + 3.0d0*d2fAtdtdA*(d2Adna2*dtdna + dAdna*d2tdna2) 1195 5 + dfAtdA*d3Adna3 + d3fAtdt3*dtdna*dtdna*dtdna 1196 6 + dfAtdt*d3tdna3 1197 7 + 3.0d0*d2fAtdt2*d2tdna2*dtdna 1198 d3fAtdna2dnb = d3fAtdA3*dAdna*dAdna*dAdnb 1199 1 + 2.0d0*d3fAtdtdA2*dAdna*dAdnb*dtdna 1200 2 + d3fAtdtdA2*dAdna*dAdna*dtdnb 1201 3 + 2.0d0*d2fAtdA2*d2Adnadnb*dAdna 1202 4 + d2fAtdA2*d2Adna2*dAdnb 1203 5 + 2.0d0*d3fAtdt2dA*dAdna*dtdna*dtdnb 1204 6 + d3fAtdt2dA*dAdnb*dtdna*dtdna 1205 7 + 2.0d0*d2fAtdtdA*(d2Adnadnb*dtdna 1206 8 + dAdna*d2tdnadnb) 1207 9 + d2fAtdtdA*(d2Adna2*dtdnb + dAdnb*d2tdna2) 1208 A + dfAtdA*d3Adna2dnb + d3fAtdt3*dtdna*dtdna*dtdnb 1209 B + dfAtdt*d3tdna2dnb 1210 C + 2.0d0*d2fAtdt2*d2tdnadnb*dtdna 1211 D + d2fAtdt2*d2tdna2*dtdnb 1212 d3fAtdnadnb2 = d3fAtdA3*dAdna*dAdnb*dAdnb 1213 1 + 2.0d0*d3fAtdtdA2*dAdna*dAdnb*dtdnb 1214 2 + d3fAtdtdA2*dAdnb*dAdnb*dtdna 1215 3 + 2.0d0*d2fAtdA2*d2Adnadnb*dAdnb 1216 4 + d2fAtdA2*d2Adnb2*dAdna 1217 5 + 2.0d0*d3fAtdt2dA*dAdnb*dtdna*dtdnb 1218 6 + d3fAtdt2dA*dAdna*dtdnb*dtdnb 1219 7 + 2.0d0*d2fAtdtdA*(d2Adnadnb*dtdnb 1220 8 + dAdnb*d2tdnadnb) 1221 9 + d2fAtdtdA*(d2Adnb2*dtdna + dAdna*d2tdnb2) 1222 A + dfAtdA*d3Adnadnb2 + d3fAtdt3*dtdna*dtdnb*dtdnb 1223 B + dfAtdt*d3tdnadnb2 1224 C + 2.0d0*d2fAtdt2*d2tdnadnb*dtdnb 1225 D + d2fAtdt2*d2tdnb2*dtdna 1226 d3fAtdnb3 = d3fAtdA3*dAdnb*dAdnb*dAdnb 1227 1 + 3.0d0*d3fAtdtdA2*dAdnb*dAdnb*dtdnb 1228 2 + 3.0d0*d2fAtdA2*d2Adnb2*dAdnb 1229 3 + 3.0d0*d3fAtdt2dA*dAdnb*dtdnb*dtdnb 1230 4 + 3.0d0*d2fAtdtdA*(d2Adnb2*dtdnb + dAdnb*d2tdnb2) 1231 5 + dfAtdA*d3Adnb3 + d3fAtdt3*dtdnb*dtdnb*dtdnb 1232 6 + dfAtdt*d3tdnb3 1233 7 + 3.0d0*d2fAtdt2*d2tdnb2*dtdnb 1234#endif 1235c 1236c arglog = 1 + BETA/GAMMA * t**2 * fAt and its derivs 1237c 1238 arglog = 1.0d0 + BETA/GAMMA*t2*fAt 1239 darglogdna = BETA/GAMMA*(2.0d0*t*dtdna*fAt 1240 & +t2*dfAtdna) 1241 darglogdnb = BETA/GAMMA*(2.0d0*t*dtdnb*fAt 1242 & +t2*dfAtdnb) 1243c 1244 if (id_cpbe.eq.3) then 1245 darglogdna = darglogdna + dBETAdr/GAMMA*t**2*fAt 1246 darglogdnb = darglogdnb + dBETAdr/GAMMA*t**2*fAt 1247 endif 1248c 1249#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 1250 d2arglogdna2 = BETA/GAMMA*(2.0d0*dtdna*dtdna*fAt 1251 & +2.0d0*t*d2tdna2*fAt 1252 & +2.0d0*t*dtdna*dfAtdna 1253 & +2.0d0*t*dtdna*dfAtdna 1254 & +t2*d2fAtdna2) 1255 d2arglogdnb2 = BETA/GAMMA*(2.0d0*dtdnb*dtdnb*fAt 1256 & +2.0d0*t*d2tdnb2*fAt 1257 & +2.0d0*t*dtdnb*dfAtdnb 1258 & +2.0d0*t*dtdnb*dfAtdnb 1259 & +t2*d2fAtdnb2) 1260 d2arglogdnadnb = BETA/GAMMA*(2.0d0*dtdna*dtdnb*fAt 1261 & +2.0d0*t*d2tdnadnb*fAt 1262 & +2.0d0*t*dtdna*dfAtdnb 1263 & +2.0d0*t*dtdnb*dfAtdna 1264 & +t2*d2fAtdnadnb) 1265c 1266 if (id_cpbe.eq.3) then 1267 d2arglogdna2 = d2arglogdna2 1268 & + d2BETAdr2/GAMMA*t**2*fAt 1269 & + dBETAdr/GAMMA*(2.0d0*t*dtdna*fAt+t*t*dfAtdna) 1270 d2arglogdnb2 = d2arglogdnb2 1271 & + d2BETAdr2/GAMMA*t**2*fAt 1272 & + dBETAdr/GAMMA*(2.0d0*t*dtdnb*fAt+t*t*dfAtdnb) 1273 d2arglogdnadnb = d2arglogdnadnb 1274 & + d2BETAdr2/GAMMA*t**2*fAt 1275 & + dBETAdr/GAMMA*(2.0d0*t*dtdna*fAt+t*t*dfAtdna) 1276 & + dBETAdr/GAMMA*(2.0d0*t*dtdnb*fAt+t*t*dfAtdnb) 1277 endif 1278c 1279#endif 1280#ifdef THIRD_DERIV 1281 d3arglogdna3 = BETA/GAMMA*(6.0d0*dfAtdna*dtdna*dtdna 1282 1 +6.0d0*dfAtdna*d2tdna2*t 1283 2 +6.0d0*d2fAtdna2*dtdna*t 1284 3 +6.0d0*fAt*dtdna*d2tdna2 1285 4 +d3fAtdna3*t2 1286 5 +2.0d0*fAt*d3tdna3*t) 1287 d3arglogdna2dnb = BETA/GAMMA*(4.0d0*dfAtdna*dtdna*dtdnb 1288 1 +2.0d0*dfAtdnb*dtdna*dtdna 1289 2 +4.0d0*dfAtdna*d2tdnadnb*t 1290 3 +2.0d0*dfAtdnb*d2tdna2*t 1291 3 +4.0d0*d2fAtdnadnb*dtdna*t 1292 4 +2.0d0*d2fAtdna2*dtdnb*t 1293 5 +4.0d0*fAt*dtdna*d2tdnadnb 1294 6 +2.0d0*fAt*dtdnb*d2tdna2 1295 7 +d3fAtdna2dnb*t2 1296 8 +2.0d0*fAt*d3tdna2dnb*t) 1297 d3arglogdnadnb2 = BETA/GAMMA*(4.0d0*dfAtdnb*dtdna*dtdnb 1298 1 +2.0d0*dfAtdna*dtdnb*dtdnb 1299 2 +4.0d0*dfAtdnb*d2tdnadnb*t 1300 3 +2.0d0*dfAtdna*d2tdnb2*t 1301 3 +4.0d0*d2fAtdnadnb*dtdnb*t 1302 4 +2.0d0*d2fAtdnb2*dtdna*t 1303 5 +4.0d0*fAt*dtdnb*d2tdnadnb 1304 6 +2.0d0*fAt*dtdna*d2tdnb2 1305 7 +d3fAtdnadnb2*t2 1306 8 +2.0d0*fAt*d3tdnadnb2*t) 1307 d3arglogdnb3 = BETA/GAMMA*(6.0d0*dfAtdnb*dtdnb*dtdnb 1308 1 +6.0d0*dfAtdnb*d2tdnb2*t 1309 2 +6.0d0*d2fAtdnb2*dtdnb*t 1310 3 +6.0d0*fAt*dtdnb*d2tdnb2 1311 4 +d3fAtdnb3*t2 1312 5 +2.0d0*fAt*d3tdnb3*t) 1313#endif 1314c 1315c H0 = GAMMA * phi**3 * log{arglog} and its derivs 1316c 1317c Daniel - I redefine the enhancement factor as a product of two 1318c functions: 1319c 1320c Mz = phi**3 1321c Nzrt = log{arglog} 1322c H0 = GAMMA * Mz * Nzrt 1323c 1324c This makes the third derivatives substantially easier to read. 1325 phi2 = phi*phi 1326 phi3 = phi2*phi 1327 Mz = phi3 1328 dMzdna = 3.0d0*phi2*dphidna 1329 dMzdnb = 3.0d0*phi2*dphidnb 1330c 1331 Nzrt = dlog(arglog) 1332 dNzrtdna = darglogdna/arglog 1333 dNzrtdnb = darglogdnb/arglog 1334c 1335 H0 = GAMMA*Mz*Nzrt 1336 dH0dna = GAMMA*(dMzdna*Nzrt + Mz*dNzrtdna) 1337 dH0dnb = GAMMA*(dMzdnb*Nzrt + Mz*dNzrtdnb) 1338#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 1339 d2Mzdna2 = 3.0d0*phi2*d2phidna2 + 6.0d0*phi*dphidna*dphidna 1340 d2Mzdnadnb = 3.0d0*phi2*d2phidnadnb 1341 1 + 6.0d0*phi*dphidna*dphidnb 1342 d2Mzdnb2 = 3.0d0*phi2*d2phidnb2 + 6.0d0*phi*dphidnb*dphidnb 1343c 1344 d2Nzrtdna2 = d2arglogdna2/arglog - 1345 1 darglogdna*darglogdna/arglog/arglog 1346 d2Nzrtdnadnb = d2arglogdnadnb/arglog - 1347 1 darglogdna*darglogdnb/arglog/arglog 1348 d2Nzrtdnb2 = d2arglogdnb2/arglog - 1349 1 darglogdnb*darglogdnb/arglog/arglog 1350c 1351 d2H0dna2 = GAMMA*( d2Mzdna2*Nzrt + 2.0d0*dMzdna*dNzrtdna 1352 & + Mz*d2Nzrtdna2) 1353 d2H0dnadnb = GAMMA*( d2Mzdnadnb*Nzrt + dMzdna*dNzrtdnb 1354 & + dMzdnb*dNzrtdna + Mz*d2Nzrtdnadnb) 1355 d2H0dnb2 = GAMMA*( d2Mzdnb2*Nzrt + 2.0d0*dMzdnb*dNzrtdnb 1356 & + Mz*d2Nzrtdnb2) 1357#endif 1358#ifdef THIRD_DERIV 1359 d3Mzdna3 = 3.0d0*phi2*d3phidna3 1360 1 + 18.0d0*phi*d2phidna2*dphidna 1361 2 + 6.0d0*dphidna*dphidna*dphidna 1362 d3Mzdna2dnb = 3.0d0*phi2*d3phidna2dnb 1363 1 + 12.0d0*phi*d2phidnadnb*dphidna 1364 2 + 6.0d0*phi*d2phidna2*dphidnb 1365 3 + 6.0d0*dphidna*dphidna*dphidnb 1366 d3Mzdnadnb2 = 3.0d0*phi2*d3phidnadnb2 1367 1 + 12.0d0*phi*d2phidnadnb*dphidnb 1368 2 + 6.0d0*phi*d2phidnb2*dphidna 1369 3 + 6.0d0*dphidna*dphidnb*dphidnb 1370 d3Mzdnb3 = 3.0d0*phi2*d3phidnb3 1371 1 + 18.0d0*phi*d2phidnb2*dphidnb 1372 2 + 6.0d0*dphidnb*dphidnb*dphidnb 1373c 1374 d3Nzrtdna3 = d3arglogdna3/arglog 1375 1 - 3.0d0*d2arglogdna2*darglogdna/arglog/arglog 1376 2 + 2.0d0*darglogdna*darglogdna*darglogdna/ 1377 3 arglog/arglog/arglog 1378 d3Nzrtdna2dnb = d3arglogdna2dnb/arglog 1379 1 - 2.0d0*d2arglogdnadnb*darglogdna/arglog/arglog 1380 2 - d2arglogdna2*darglogdnb/arglog/arglog 1381 3 + 2.0d0*darglogdna*darglogdna*darglogdnb/ 1382 4 arglog/arglog/arglog 1383 d3Nzrtdnadnb2 = d3arglogdnadnb2/arglog 1384 1 - 2.0d0*d2arglogdnadnb*darglogdnb/arglog/arglog 1385 2 - d2arglogdnb2*darglogdna/arglog/arglog 1386 3 + 2.0d0*darglogdna*darglogdnb*darglogdnb/ 1387 4 arglog/arglog/arglog 1388 d3Nzrtdnb3 = d3arglogdnb3/arglog 1389 1 - 3.0d0*d2arglogdnb2*darglogdnb/arglog/arglog 1390 2 + 2.0d0*darglogdnb*darglogdnb*darglogdnb/ 1391 3 arglog/arglog/arglog 1392c 1393 d3H0dna3 = GAMMA*( d3Mzdna3*Nzrt 1394 1 + 3.0d0*d2Mzdna2*dNzrtdna 1395 2 + 3.0d0*dMzdna*d2Nzrtdna2 1396 3 + Mz*d3Nzrtdna3) 1397 d3H0dna2dnb = GAMMA*( d3Mzdna2dnb*Nzrt 1398 1 + 2.0d0*d2Mzdnadnb*dNzrtdna 1399 2 + d2Mzdna2*dNzrtdnb 1400 3 + 2.0d0*dMzdna*d2Nzrtdnadnb 1401 4 + dMzdnb*d2Nzrtdna2 1402 5 + Mz*d3Nzrtdna2dnb) 1403 d3H0dnadnb2 = GAMMA*( d3Mzdnadnb2*Nzrt 1404 1 + 2.0d0*d2Mzdnadnb*dNzrtdnb 1405 2 + d2Mzdnb2*dNzrtdna 1406 3 + 2.0d0*dMzdnb*d2Nzrtdnadnb 1407 4 + dMzdna*d2Nzrtdnb2 1408 5 + Mz*d3Nzrtdnadnb2) 1409 d3H0dnb3 = GAMMA*( d3Mzdnb3*Nzrt 1410 1 + 3.0d0*d2Mzdnb2*dNzrtdnb 1411 2 + 3.0d0*dMzdnb*d2Nzrtdnb2 1412 3 + Mz*d3Nzrtdnb3) 1413#endif 1414c 1415c Now we update Ec, Amat, and Amat2 1416c 1417c Daniel (3-5-13): I believe that lfac is set to false for all 1418c calculations that I do. 1419 if (lfac) then 1420 Ec = Ec+nepsc*qwght(n)*fac 1421 if (ldew) ffunc(n)=ffunc(n)+nepsc*fac 1422 endif 1423 if (nlfac) then 1424 Ec = Ec+(H0*rhoval)*qwght(n)*fac 1425 if (ldew) ffunc(n)=ffunc(n)+(H0*rhoval)*fac 1426 endif 1427 if (lfac) then 1428 Amat(n,1) = Amat(n,1) + dnepscdn(1)*fac 1429 if (ipol.eq.2) Amat(n,2) = Amat(n,2) + dnepscdn(2)*fac 1430#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 1431 Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) 1432 & + d2nepscdn2(D2_RA_RA)*fac 1433 Amat2(n,D2_RA_RB) = Amat2(n,D2_RA_RB) 1434 & + d2nepscdn2(D2_RA_RB)*fac 1435 if (ipol.eq.2) 1436 & Amat2(n,D2_RB_RB) = Amat2(n,D2_RB_RB) 1437 & + d2nepscdn2(D2_RB_RB)*fac 1438#endif 1439c Daniel (7-31-12): XC third derivatives (LDA part) 1440#ifdef THIRD_DERIV 1441 Amat3(n,D3_RA_RA_RA) = Amat3(n,D3_RA_RA_RA) 1442 1 + d3nepscdn3(D3_RA_RA_RA)*fac 1443 Amat3(n,D3_RA_RA_RB) = Amat3(n,D3_RA_RA_RB) 1444 1 + d3nepscdn3(D3_RA_RA_RB)*fac 1445 Amat3(n,D3_RA_RB_RB) = Amat3(n,D3_RA_RB_RB) 1446 1 + d3nepscdn3(D3_RA_RB_RB)*fac 1447 if (ipol.eq.2) 1448 1 Amat3(n,D3_RB_RB_RB) = Amat3(n,D3_RB_RB_RB) 1449 2 + d3nepscdn3(D3_RB_RB_RB)*fac 1450#endif 1451 endif 1452 if (nlfac)then 1453 Amat(n,1) = Amat(n,1) + (H0 + 1454 & rhoval*dH0dna)*fac 1455 if (ipol.eq.2) Amat(n,2) = Amat(n,2) + (H0 + 1456 & rhoval*dH0dnb)*fac 1457#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 1458 Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) 1459 & + (2.d0*dH0dna + rhoval*d2H0dna2)*fac 1460 Amat2(n,D2_RA_RB) = Amat2(n,D2_RA_RB) 1461 & + (dH0dna + dH0dnb + rhoval*d2H0dnadnb)*fac 1462 if (ipol.eq.2) 1463 & Amat2(n,D2_RB_RB) = Amat2(n,D2_RB_RB) 1464 & + (2.d0*dH0dnb + rhoval*d2H0dnb2)*fac 1465#endif 1466c Daniel (7-31-12): XC third derivatives (GGA part) 1467#ifdef THIRD_DERIV 1468 Amat3(n,D3_RA_RA_RA) = Amat3(n,D3_RA_RA_RA) 1469 1 + (3.0d0*d2H0dna2 + rhoval*d3H0dna3)*fac 1470 Amat3(n,D3_RA_RA_RB) = Amat3(n,D3_RA_RA_RB) 1471 1 + (d2H0dna2 + 2.0d0*d2H0dnadnb + rhoval*d3H0dna2dnb)*fac 1472 Amat3(n,D3_RA_RB_RB) = Amat3(n,D3_RA_RB_RB) 1473 1 + (d2H0dnb2 + 2.0d0*d2H0dnadnb + rhoval*d3H0dnadnb2)*fac 1474 if (ipol.eq.2) 1475 1 Amat3(n,D3_RB_RB_RB) = Amat3(n,D3_RB_RB_RB) 1476 2 + (3.0d0*d2H0dnb2 + rhoval*d3H0dnb3)*fac 1477#endif 1478 endif 1479c 1480c Now we go into gradient-correction parts 1481c Note that the functional depends on |Nabla n| through "t" only 1482c 1483 if (dsqgamma.gt.TOLL)then 1484 dtdg = 0.25d0*rphi*rks*rrho/dsqgamma 1485 dfAtdg = dfAtdt*dtdg 1486 darglogdg = BETA/GAMMA*(2.0d0*t*dtdg*fAt+t2*dfAtdg) 1487c 1488 dNzrtdg = darglogdg/arglog 1489c 1490 dH0dg = GAMMA*Mz*dNzrtdg 1491 if (ipol.eq.1) then 1492 Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + dH0dg*rhoval*fac 1493 Cmat(n,D1_GAB) = Cmat(n,D1_GAB) + dH0dg*rhoval*fac*2.0d0 1494 else 1495 Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + dH0dg*rhoval*fac 1496 Cmat(n,D1_GAB) = Cmat(n,D1_GAB) + dH0dg*rhoval*fac*2.0d0 1497 Cmat(n,D1_GBB) = Cmat(n,D1_GBB) + dH0dg*rhoval*fac 1498 endif 1499#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 1500 d2tdg2 = -0.125d0*rphi*rks*rrho/ 1501 & (dsqgamma*dsqgamma*dsqgamma) 1502 d2tdnadg = -dtdg*rrho-dtdg*rphi*dphidna 1503 & -dtdg*rks*dksdna 1504 d2tdnbdg = -dtdg*rrho-dtdg*rphi*dphidnb 1505 & -dtdg*rks*dksdnb 1506 d2fAtdg2 = d2fAtdt2*dtdg*dtdg+dfAtdt*d2tdg2 1507 d2fAtdnadg = d2fAtdt2*dtdg*dtdna 1508 & +d2fAtdtdA*dtdg*dAdna 1509 & +dfAtdt*d2tdnadg 1510 d2fAtdnbdg = d2fAtdt2*dtdg*dtdnb 1511 & +d2fAtdtdA*dtdg*dAdnb 1512 & +dfAtdt*d2tdnbdg 1513 d2arglogdnadg = BETA/GAMMA*(2.0d0*dtdna*dtdg*fAt 1514 & +2.0d0*t*d2tdnadg*fAt 1515 & +2.0d0*t*dtdg*dfAtdna 1516 & +2.0d0*t*dtdna*dfAtdg 1517 & +t2*d2fAtdnadg) 1518 d2arglogdnbdg = BETA/GAMMA*(2.0d0*dtdnb*dtdg*fAt 1519 & +2.0d0*t*d2tdnbdg*fAt 1520 & +2.0d0*t*dtdg*dfAtdnb 1521 & +2.0d0*t*dtdnb*dfAtdg 1522 & +t2*d2fAtdnbdg) 1523c 1524 if (id_cpbe.eq.3) then 1525 d2arglogdnadg = d2arglogdnadg 1526 & + dBETAdr/GAMMA*(2.0d0*t*dtdg*fAt+t*t*dfAtdg) 1527 d2arglogdnbdg = d2arglogdnbdg 1528 & + dBETAdr/GAMMA*(2.0d0*t*dtdg*fAt+t*t*dfAtdg) 1529 endif 1530c 1531 d2arglogdg2 = BETA/GAMMA*(2.0d0*dtdg*dtdg*fAt 1532 & +2.0d0*t*d2tdg2*fAt 1533 & +2.0d0*t*dtdg*dfAtdg 1534 & +2.0d0*t*dtdg*dfAtdg 1535 & +t2*d2fAtdg2) 1536c 1537 d2Nzrtdg2 = d2arglogdg2/arglog - 1538 1 darglogdg*darglogdg/arglog/arglog 1539 d2Nzrtdnadg = d2arglogdnadg/arglog - 1540 1 darglogdna*darglogdg/arglog/arglog 1541 d2Nzrtdnbdg = d2arglogdnbdg/arglog - 1542 1 darglogdnb*darglogdg/arglog/arglog 1543c 1544 d2H0dnadg = GAMMA*dMzdna*dNzrtdg 1545 & + GAMMA*Mz*d2Nzrtdnadg 1546 d2H0dnbdg = GAMMA*dMzdnb*dNzrtdg 1547 & + GAMMA*Mz*d2Nzrtdnbdg 1548 d2H0dg2 = GAMMA*Mz*d2Nzrtdg2 1549c 1550c ORIGINAL CODE 1551c Cmat2(n,D2_RA_GAA) = Cmat2(n,D2_RA_GAA) 1552c & + (dH0dg + d2H0dnadg*rhoval)*fac 1553c Cmat2(n,D2_RA_GAB) = Cmat2(n,D2_RA_GAB) 1554c & + 2.0d0*(dH0dg + d2H0dnadg*rhoval)*fac 1555c Cmat2(n,D2_RA_GBB) = Cmat2(n,D2_RA_GBB) 1556c & + (dH0dg + d2H0dnadg*rhoval)*fac 1557c Cmat2(n,D2_RB_GAA) = Cmat2(n,D2_RB_GAA) 1558c & + (dH0dg + d2H0dnbdg*rhoval)*fac 1559c Cmat2(n,D2_RB_GAB) = Cmat2(n,D2_RB_GAB) 1560c & + 2.0d0*(dH0dg + d2H0dnbdg*rhoval)*fac 1561c Cmat2(n,D2_RB_GBB) = Cmat2(n,D2_RB_GBB) 1562c & + (dH0dg + d2H0dnbdg*rhoval)*fac 1563c Cmat2(n,D2_GAA_GAA) = Cmat2(n,D2_GAA_GAA) 1564c & + d2H0dg2*rhoval*fac 1565c Cmat2(n,D2_GAA_GAB) = Cmat2(n,D2_GAA_GAB) 1566c & + 2.0d0*d2H0dg2*rhoval*fac 1567c Cmat2(n,D2_GAA_GBB) = Cmat2(n,D2_GAA_GBB) 1568c & + d2H0dg2*rhoval*fac 1569c Cmat2(n,D2_GAB_GAB) = Cmat2(n,D2_GAB_GAB) 1570c & + 4.0d0*d2H0dg2*rhoval*fac 1571c Cmat2(n,D2_GAB_GBB) = Cmat2(n,D2_GAB_GBB) 1572c & + 2.0d0*d2H0dg2*rhoval*fac 1573c Cmat2(n,D2_GBB_GBB) = Cmat2(n,D2_GBB_GBB) 1574c & + d2H0dg2*rhoval*fac 1575c ORIGINAL CODE 1576c 1577c Daniel (3-11-13): Only do beta terms if they are needed 1578c since the extra calculations are unnecessary for restricted 1579c calculations. 1580 Cmat2(n,D2_RA_GAA) = Cmat2(n,D2_RA_GAA) 1581 & + (dH0dg + d2H0dnadg*rhoval)*fac 1582 Cmat2(n,D2_RA_GAB) = Cmat2(n,D2_RA_GAB) 1583 & + 2.0d0*(dH0dg + d2H0dnadg*rhoval)*fac 1584 Cmat2(n,D2_RA_GBB) = Cmat2(n,D2_RA_GBB) 1585 & + (dH0dg + d2H0dnadg*rhoval)*fac 1586c 1587 Cmat2(n,D2_GAA_GAA) = Cmat2(n,D2_GAA_GAA) 1588 & + d2H0dg2*rhoval*fac 1589 Cmat2(n,D2_GAA_GAB) = Cmat2(n,D2_GAA_GAB) 1590 & + 2.0d0*d2H0dg2*rhoval*fac 1591 Cmat2(n,D2_GAA_GBB) = Cmat2(n,D2_GAA_GBB) 1592 & + d2H0dg2*rhoval*fac 1593 Cmat2(n,D2_GAB_GAB) = Cmat2(n,D2_GAB_GAB) 1594 & + 4.0d0*d2H0dg2*rhoval*fac 1595 if (ipol.eq.2) then 1596 Cmat2(n,D2_RB_GAA) = Cmat2(n,D2_RB_GAA) 1597 & + (dH0dg + d2H0dnbdg*rhoval)*fac 1598 Cmat2(n,D2_RB_GAB) = Cmat2(n,D2_RB_GAB) 1599 & + 2.0d0*(dH0dg + d2H0dnbdg*rhoval)*fac 1600 Cmat2(n,D2_RB_GBB) = Cmat2(n,D2_RB_GBB) 1601 & + (dH0dg + d2H0dnbdg*rhoval)*fac 1602c 1603 Cmat2(n,D2_GAB_GBB) = Cmat2(n,D2_GAB_GBB) 1604 & + 2.0d0*d2H0dg2*rhoval*fac 1605 Cmat2(n,D2_GBB_GBB) = Cmat2(n,D2_GBB_GBB) 1606 & + d2H0dg2*rhoval*fac 1607 endif 1608#endif 1609#ifdef THIRD_DERIV 1610c Derivatives of t, the reduced density gradient 1611 d3tdg3 = 0.1875d0*rphi*rks*rrho/dsqgamma/dsqgamma/ 1612 1 dsqgamma/dsqgamma/dsqgamma 1613 d3tdna2dg = dtdg*( dphidna*dphidna*rphi2 1614 1 - d2phidna2*rphi 1615 2 + dksdna*dksdna*rks2 1616 3 - d2ksdna2*rks 1617 4 + rrho2) 1618 5 - d2tdnadg*( dphidna*rphi 1619 6 + dksdna*rks 1620 7 + rrho) 1621 d3tdnadnbdg = dtdg*( dphidna*dphidnb*rphi2 1622 1 - d2phidnadnb*rphi 1623 2 + dksdna*dksdnb*rks2 1624 3 - d2ksdnadnb*rks 1625 4 + rrho2) 1626 5 - d2tdnbdg*( dphidna*rphi 1627 6 + dksdna*rks 1628 7 + rrho) 1629 d3tdnbdnadg = dtdg*( dphidna*dphidnb*rphi2 1630 1 - d2phidnadnb*rphi 1631 2 + dksdna*dksdnb*rks2 1632 3 - d2ksdnadnb*rks 1633 4 + rrho2) 1634 5 - d2tdnadg*( dphidnb*rphi 1635 6 + dksdnb*rks 1636 7 + rrho) 1637 d3tdnb2dg = dtdg*( dphidnb*dphidnb*rphi2 1638 1 - d2phidnb2*rphi 1639 2 + dksdnb*dksdnb*rks2 1640 3 - d2ksdnb2*rks 1641 4 + rrho2) 1642 5 - d2tdnbdg*( dphidnb*rphi 1643 6 + dksdnb*rks 1644 7 + rrho) 1645 d3tdnadg2 = -d2tdg2*( dphidna*rphi + dksdna*rks 1646 1 + rrho) 1647 d3tdnbdg2 = -d2tdg2*( dphidnb*rphi + dksdnb*rks 1648 1 + rrho) 1649c Derivatives of the rational function fAt 1650 d3fAtdg3 = d3fAtdt3*dtdg*dtdg*dtdg 1651 1 + 3.0d0*d2fAtdt2*d2tdg2*dtdg 1652 2 + dfAtdt*d3tdg3 1653 d3fAtdna2dg = d3fAtdtdA2*dAdna*dAdna*dtdg 1654 1 + 2.0d0*d3fAtdt2dA*dAdna*dtdna*dtdg 1655 2 + d2fAtdtdA*( d2Adna2*dtdg 1656 3 + 2.0d0*dAdna*d2tdnadg) 1657 4 + d3fAtdt3*dtdna*dtdna*dtdg 1658 5 + d2fAtdt2*( d2tdna2*dtdg 1659 6 + 2.0d0*d2tdnadg*dtdna) 1660 7 + dfAtdt*d3tdna2dg 1661 d3fAtdnadnbdg = d3fAtdtdA2*dAdna*dAdnb*dtdg 1662 1 + d3fAtdt2dA*( dAdna*dtdnb 1663 2 + dAdnb*dtdna)*dtdg 1664 3 + d2fAtdtdA*( d2Adnadnb*dtdg 1665 4 + dAdna*d2tdnbdg 1666 5 + dAdnb*d2tdnadg) 1667 6 + d3fAtdt3*dtdna*dtdnb*dtdg 1668 7 + d2fAtdt2*( d2tdnadnb*dtdg 1669 8 + d2tdnadg*dtdnb 1670 9 + d2tdnbdg*dtdna) 1671 A + dfAtdt*d3tdnadnbdg 1672 d3fAtdnbdnadg = d3fAtdnadnbdg 1673 d3fAtdnb2dg = d3fAtdtdA2*dAdnb*dAdnb*dtdg 1674 1 + 2.0d0*d3fAtdt2dA*dAdnb*dtdnb*dtdg 1675 2 + d2fAtdtdA*( d2Adnb2*dtdg 1676 3 + 2.0d0*dAdnb*d2tdnbdg) 1677 4 + d3fAtdt3*dtdnb*dtdnb*dtdg 1678 5 + d2fAtdt2*( d2tdnb2*dtdg 1679 6 + 2.0d0*d2tdnbdg*dtdnb) 1680 7 + dfAtdt*d3tdnb2dg 1681 d3fAtdnadg2 = d3fAtdt2dA*dAdna*dtdg*dtdg 1682 1 + d2fAtdtdA*dAdna*d2tdg2 1683 2 + d3fAtdt3*dtdna*dtdg*dtdg 1684 3 + d2fAtdt2*( 2.0d0*d2tdnadg*dtdg 1685 4 + dtdna*d2tdg2) 1686 5 + dfAtdt*d3tdnadg2 1687 d3fAtdnbdg2 = d3fAtdt2dA*dAdnb*dtdg*dtdg 1688 1 + d2fAtdtdA*dAdnb*d2tdg2 1689 2 + d3fAtdt3*dtdnb*dtdg*dtdg 1690 3 + d2fAtdt2*( 2.0d0*d2tdnbdg*dtdg 1691 4 + dtdnb*d2tdg2) 1692 5 + dfAtdt*d3tdnbdg2 1693c Derivatives of the logarithm argument arglog 1694 d3arglogdg3 = BETA/GAMMA*(d3fAtdg3*t2 1695 1 + 6.0d0*d2fAtdg2*dtdg*t 1696 2 + 6.0d0*dfAtdg*( dtdg*dtdg 1697 3 + d2tdg2*t) 1698 4 + 6.0d0*fAt*dtdg*d2tdg2 1699 5 + 2.0d0*fAt*d3tdg3*t) 1700 d3arglogdna2dg = 2.0d0*BETA/GAMMA*( 1701 1 2.0d0*d2fAtdnadg*dtdna*t 1702 2 + 2.0d0*dfAtdna*d2tdnadg*t 1703 3 + 2.0d0*dfAtdna*dtdna*dtdg 1704 4 + dfAtdg*dtdna*dtdna 1705 5 + 2.0d0*fAt*d2tdnadg*dtdna 1706 6 + d2fAtdna2*dtdg*t 1707 7 + dfAtdg*d2tdna2*t 1708 8 + fAt*d3tdna2dg*t 1709 9 + fAt*d2tdna2*dtdg 1710 A + 0.50d0*d3fAtdna2dg*t2) 1711 d3arglogdnadnbdg = 2.0d0*BETA/GAMMA*( 1712 1 d2fAtdnadg*dtdnb*t 1713 2 + d2fAtdnbdg*dtdna*t 1714 3 + dfAtdna*d2tdnbdg*t 1715 4 + dfAtdnb*d2tdnadg*t 1716 5 + dfAtdna*dtdnb*dtdg 1717 6 + dfAtdnb*dtdna*dtdg 1718 7 + dfAtdg*dtdna*dtdnb 1719 8 + fAt*d2tdnadg*dtdnb 1720 9 + fAt*d2tdnbdg*dtdna 1721 A + d2fAtdnadnb*dtdg*t 1722 B + dfAtdg*d2tdnadnb*t 1723 C + fAt*d3tdnadnbdg*t 1724 D + fAt*d2tdnadnb*dtdg 1725 E + 0.50d0*d3fAtdnadnbdg*t2) 1726 d3arglogdnbdnadg = 2.0d0*BETA/GAMMA*( 1727 1 d2fAtdnbdg*dtdna*t 1728 2 + d2fAtdnadg*dtdnb*t 1729 3 + dfAtdnb*d2tdnadg*t 1730 4 + dfAtdna*d2tdnbdg*t 1731 5 + dfAtdnb*dtdna*dtdg 1732 6 + dfAtdna*dtdnb*dtdg 1733 7 + dfAtdg*dtdnb*dtdna 1734 8 + fAt*d2tdnbdg*dtdna 1735 9 + fAt*d2tdnadg*dtdnb 1736 A + d2fAtdnadnb*dtdg*t 1737 B + dfAtdg*d2tdnadnb*t 1738 C + fAt*d3tdnbdnadg*t 1739 D + fAt*d2tdnadnb*dtdg 1740 E + 0.50d0*d3fAtdnbdnadg*t2) 1741 d3arglogdnb2dg = 2.0d0*BETA/GAMMA*( 1742 1 2.0d0*d2fAtdnbdg*dtdnb*t 1743 2 + 2.0d0*dfAtdnb*d2tdnbdg*t 1744 3 + 2.0d0*dfAtdnb*dtdnb*dtdg 1745 4 + dfAtdg*dtdnb*dtdnb 1746 5 + 2.0d0*fAt*d2tdnbdg*dtdnb 1747 6 + d2fAtdnb2*dtdg*t 1748 7 + dfAtdg*d2tdnb2*t 1749 8 + fAt*d3tdnb2dg*t 1750 9 + fAt*d2tdnb2*dtdg 1751 A + 0.50d0*d3fAtdnb2dg*t2) 1752 d3arglogdnadg2 = 2.0d0*BETA/GAMMA*( 1753 1 2.0d0*d2fAtdnadg*dtdg*t 1754 2 + dfAtdna*d2tdg2*t 1755 3 + dfAtdna*dtdg*dtdg 1756 4 + d2fAtdg2*dtdna*t 1757 5 + 2.0d0*dfAtdg*d2tdnadg*t 1758 6 + 2.0d0*dfAtdg*dtdna*dtdg 1759 7 + 2.0d0*fAt*d2tdnadg*dtdg 1760 8 + fAt*dtdna*d2tdg2 1761 9 + fAt*d3tdnadg2*t 1762 A + 0.50d0*d3fAtdnadg2*t2) 1763 d3arglogdnbdg2 = 2.0d0*BETA/GAMMA*( 1764 1 2.0d0*d2fAtdnbdg*dtdg*t 1765 2 + dfAtdnb*d2tdg2*t 1766 3 + dfAtdnb*dtdg*dtdg 1767 4 + d2fAtdg2*dtdnb*t 1768 5 + 2.0d0*dfAtdg*d2tdnbdg*t 1769 6 + 2.0d0*dfAtdg*dtdnb*dtdg 1770 7 + 2.0d0*fAt*d2tdnbdg*dtdg 1771 8 + fAt*dtdnb*d2tdg2 1772 9 + fAt*d3tdnbdg2*t 1773 A + 0.50d0*d3fAtdnbdg2*t2) 1774c Derivatives of the correlation enhancement factor 1775 d3Nzrtdg3 = d3arglogdg3/arglog 1776 1 - 3.0d0*d2arglogdg2*darglogdg/arglog/arglog 1777 2 + 2.0d0*darglogdg*darglogdg*darglogdg/arglog/ 1778 3 arglog/arglog 1779 d3Nzrtdnadg2 = d3arglogdnadg2/arglog 1780 1 - (2.0d0*d2arglogdnadg*darglogdg + 1781 2 darglogdna*d2arglogdg2)/arglog/arglog 1782 3 + 2.0d0*darglogdna*darglogdg*darglogdg/ 1783 4 arglog/arglog/arglog 1784 d3Nzrtdnbdg2 = d3arglogdnbdg2/arglog 1785 1 - (2.0d0*d2arglogdnbdg*darglogdg + 1786 2 darglogdnb*d2arglogdg2)/arglog/arglog 1787 3 + 2.0d0*darglogdnb*darglogdg*darglogdg/ 1788 4 arglog/arglog/arglog 1789 d3Nzrtdna2dg = d3arglogdna2dg/arglog 1790 1 - (2.0d0*d2arglogdnadg*darglogdna + 1791 2 d2arglogdna2*darglogdg)/arglog/arglog 1792 3 + 2.0d0*darglogdna*darglogdna*darglogdg/ 1793 4 arglog/arglog/arglog 1794 d3Nzrtdnadnbdg = d3arglogdnadnbdg/arglog 1795 1 - (d2arglogdnadg*darglogdnb + 1796 2 d2arglogdnbdg*darglogdna + 1797 3 d2arglogdnadnb*darglogdg)/arglog/arglog 1798 4 + 2.0d0*darglogdna*darglogdnb*darglogdg/ 1799 5 arglog/arglog/arglog 1800 d3Nzrtdnb2dg = d3arglogdnb2dg/arglog 1801 1 - (2.0d0*d2arglogdnbdg*darglogdnb + 1802 2 d2arglogdnb2*darglogdg)/arglog/arglog 1803 3 + 2.0d0*darglogdnb*darglogdnb*darglogdg/ 1804 4 arglog/arglog/arglog 1805c 1806 d3H0dg3 = GAMMA*Mz*d3Nzrtdg3 1807 d3H0dna2dg = GAMMA*( d2Mzdna2*dNzrtdg 1808 1 + 2.0d0*dMzdna*d2Nzrtdnadg 1809 2 + Mz*d3Nzrtdna2dg) 1810 d3H0dnadnbdg = GAMMA*( d2Mzdnadnb*dNzrtdg 1811 1 + dMzdna*d2Nzrtdnbdg 1812 2 + dMzdnb*d2Nzrtdnadg 1813 3 + Mz*d3Nzrtdnadnbdg) 1814 d3H0dnbdnadg = d3H0dnadnbdg 1815 d3H0dnb2dg = GAMMA*( d2Mzdnb2*dNzrtdg 1816 1 + 2.0d0*dMzdnb*d2Nzrtdnbdg 1817 2 + Mz*d3Nzrtdnb2dg) 1818 d3H0dnadg2 = GAMMA*( dMzdna*d2Nzrtdg2 1819 1 + Mz*d3Nzrtdnadg2) 1820 d3H0dnbdg2 = GAMMA*( dMzdnb*d2Nzrtdg2 1821 1 + Mz*d3Nzrtdnbdg2) 1822c There are 31 unique 3rd order functional derivatives involving the 1823c gradient of the electron density. Note that many permutations of 1824c indices are identical so only one permutation is stored. 1825c Keep in mind that gamma = gaa + 2*gab + gbb for a spin polarized 1826c functional, so dgamma/dgaa = 1 and dgamma/dgab = 2 1827c 1828c Notes so far (may reduce storage requirements, worry about this when 1829c we know the code works) 1830c 1. D3_RA_RA_GAA = D3_RA_RA_GBB 1831c 2. D3_RA_RB_GAA = D3_RA_RB_GBB 1832c 3. D3_RB_RB_GAA = D3_RB_RB_GBB 1833c 4. D3_RA_GAA_GAA = D3_RA_GAA_GBB = D3_RA_GBB_GBB 1834c 5. D3_RA_GAA_GAB = D3_RA_GAB_GBB 1835c 6. D3_RB_GAA_GAA = D3_RB_GAA_GBB = D3_RB_GBB_GBB 1836c 7. D3_RB_GAA_GAB = D3_RB_GAB_GBB 1837c 8. D3_GAA_GAA_GAA = D3_GAA_GAA_GBB = D3_GAA_GBB_GBB = D3_GBB_GBB_GBB 1838c 9. D3_GAA_GAA_GAB = D3_GAA_GAB_GBB = D3_GAB_GBB_GBB 1839c 10. D3_GAA_GAB_GAB = D3_GAB_GAB_GBB 1840c 1841c It looks like we can remove 15 of the 31 elements, which could be a 1842c huge benefit here. 1843c 1844c Mixed derivatives dradradg 1845 Cmat3(n,D3_RA_RA_GAA) = Cmat3(n,D3_RA_RA_GAA) 1846 1 + (2.0d0*d2H0dnadg + d3H0dna2dg*rhoval)*fac 1847 Cmat3(n,D3_RA_RA_GAB) = Cmat3(n,D3_RA_RA_GAB) 1848 1 + 2.0d0*(2.0d0*d2H0dnadg + d3H0dna2dg*rhoval)*fac 1849 Cmat3(n,D3_RA_RA_GBB) = Cmat3(n,D3_RA_RA_GBB) 1850 1 + (2.0d0*d2H0dnadg + d3H0dna2dg*rhoval)*fac 1851c Mixed derivatives dradrbdg 1852 Cmat3(n,D3_RA_RB_GAA) = Cmat3(n,D3_RA_RB_GAA) 1853 1 + ( d2H0dnadg + d2H0dnbdg 1854 2 + d3H0dnadnbdg*rhoval)*fac 1855 Cmat3(n,D3_RA_RB_GAB) = Cmat3(n,D3_RA_RB_GAB) 1856 1 + 2.0d0*( d2H0dnadg + d2H0dnbdg 1857 2 + d3H0dnadnbdg*rhoval)*fac 1858 Cmat3(n,D3_RA_RB_GBB) = Cmat3(n,D3_RA_RB_GBB) 1859 1 + ( d2H0dnadg + d2H0dnbdg 1860 2 + d3H0dnadnbdg*rhoval)*fac 1861c Mixed derivatives drbdrbdg 1862 if (ipol.eq.2) then 1863 Cmat3(n,D3_RB_RB_GAA) = Cmat3(n,D3_RB_RB_GAA) 1864 1 + (2.0d0*d2H0dnbdg + d3H0dnb2dg*rhoval)*fac 1865 Cmat3(n,D3_RB_RB_GAB) = Cmat3(n,D3_RB_RB_GAB) 1866 1 + 2.0d0*(2.0d0*d2H0dnbdg + d3H0dnb2dg*rhoval)*fac 1867 Cmat3(n,D3_RB_RB_GBB) = Cmat3(n,D3_RB_RB_GBB) 1868 1 + (2.0d0*d2H0dnbdg + d3H0dnb2dg*rhoval)*fac 1869 endif 1870c Mixed derivatives dradgdg 1871 Cmat3(n,D3_RA_GAA_GAA) = Cmat3(n,D3_RA_GAA_GAA) 1872 1 + (d2H0dg2 + d3H0dnadg2*rhoval)*fac 1873 Cmat3(n,D3_RA_GAA_GAB) = Cmat3(n,D3_RA_GAA_GAB) 1874 1 + 2.0d0*(d2H0dg2 + d3H0dnadg2*rhoval)*fac 1875 Cmat3(n,D3_RA_GAA_GBB) = Cmat3(n,D3_RA_GAA_GBB) 1876 1 + (d2H0dg2 + d3H0dnadg2*rhoval)*fac 1877 Cmat3(n,D3_RA_GAB_GAB) = Cmat3(n,D3_RA_GAB_GAB) 1878 1 + 4.0d0*(d2H0dg2 + d3H0dnadg2*rhoval)*fac 1879 Cmat3(n,D3_RA_GAB_GBB) = Cmat3(n,D3_RA_GAB_GBB) 1880 1 + 2.0d0*(d2H0dg2 + d3H0dnadg2*rhoval)*fac 1881 Cmat3(n,D3_RA_GBB_GBB) = Cmat3(n,D3_RA_GBB_GBB) 1882 1 + (d2H0dg2 + d3H0dnadg2*rhoval)*fac 1883c Mixed derivatives drbdgdg 1884 if (ipol.eq.2) then 1885 Cmat3(n,D3_RB_GAA_GAA) = Cmat3(n,D3_RB_GAA_GAA) 1886 1 + (d2H0dg2 + d3H0dnbdg2*rhoval)*fac 1887 Cmat3(n,D3_RB_GAA_GAB) = Cmat3(n,D3_RB_GAA_GAB) 1888 1 + 2.0d0*(d2H0dg2 + d3H0dnbdg2*rhoval)*fac 1889 Cmat3(n,D3_RB_GAA_GBB) = Cmat3(n,D3_RB_GAA_GBB) 1890 1 + (d2H0dg2 + d3H0dnbdg2*rhoval)*fac 1891 Cmat3(n,D3_RB_GAB_GAB) = Cmat3(n,D3_RB_GAB_GAB) 1892 1 + 4.0d0*(d2H0dg2 + d3H0dnbdg2*rhoval)*fac 1893 Cmat3(n,D3_RB_GAB_GBB) = Cmat3(n,D3_RB_GAB_GBB) 1894 1 + 2.0d0*(d2H0dg2 + d3H0dnbdg2*rhoval)*fac 1895 Cmat3(n,D3_RB_GBB_GBB) = Cmat3(n,D3_RB_GBB_GBB) 1896 1 + (d2H0dg2 + d3H0dnbdg2*rhoval)*fac 1897 endif 1898c Derivatives dgaadgdg 1899 Cmat3(n,D3_GAA_GAA_GAA) = Cmat3(n,D3_GAA_GAA_GAA) 1900 1 + (d3H0dg3*rhoval)*fac 1901 Cmat3(n,D3_GAA_GAA_GAB) = Cmat3(n,D3_GAA_GAA_GAB) 1902 1 + 2.0d0*(d3H0dg3*rhoval)*fac 1903 Cmat3(n,D3_GAA_GAA_GBB) = Cmat3(n,D3_GAA_GAA_GBB) 1904 1 + (d3H0dg3*rhoval)*fac 1905 Cmat3(n,D3_GAA_GAB_GAB) = Cmat3(n,D3_GAA_GAB_GAB) 1906 1 + 4.0d0*(d3H0dg3*rhoval)*fac 1907 Cmat3(n,D3_GAA_GAB_GBB) = Cmat3(n,D3_GAA_GAB_GBB) 1908 1 + 2.0d0*(d3H0dg3*rhoval)*fac 1909 Cmat3(n,D3_GAA_GBB_GBB) = Cmat3(n,D3_GAA_GBB_GBB) 1910 1 + (d3H0dg3*rhoval)*fac 1911c Derivatives dgabdgdg 1912 Cmat3(n,D3_GAB_GAB_GAB) = Cmat3(n,D3_GAB_GAB_GAB) 1913 1 + 8.0d0*(d3H0dg3*rhoval)*fac 1914 if (ipol.eq.2) then 1915 Cmat3(n,D3_GAB_GAB_GBB) = Cmat3(n,D3_GAB_GAB_GBB) 1916 1 + 4.0d0*(d3H0dg3*rhoval)*fac 1917 Cmat3(n,D3_GAB_GBB_GBB) = Cmat3(n,D3_GAB_GBB_GBB) 1918 1 + 2.0d0*(d3H0dg3*rhoval)*fac 1919c Derivatives dgbbdgdg 1920 Cmat3(n,D3_GBB_GBB_GBB) = Cmat3(n,D3_GBB_GBB_GBB) 1921 1 + (d3H0dg3*rhoval)*fac 1922 endif 1923#endif 1924 endif 1925 20 continue 1926c 1927 return 1928 end 1929c 1930#ifndef SECOND_DERIV 1931#define SECOND_DERIV 1932c 1933c Compile source again for the 2nd derivative case 1934c 1935#include "xc_pbe96.F" 1936#endif 1937#ifndef THIRD_DERIV 1938#define THIRD_DERIV 1939c 1940c Compile source again for the 3rd derivative case 1941c 1942#include "xc_pbe96.F" 1943#endif 1944c $Id$ 1945