1#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 2 Subroutine xc_lyp88(tol_rho, fac, rho, delrho, 3 & Amat, Cmat, nq, ipol, Ec, qwght, ldew, func) 4#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV) 5 Subroutine xc_lyp88_d2(tol_rho, fac, rho, delrho, 6 & Amat, Amat2, Cmat, Cmat2, nq, ipol, Ec, 7 & qwght, ldew, func) 8#else 9 Subroutine xc_lyp88_d3(tol_rho, fac, rho, delrho, 10 & Amat, Amat2, Amat3, Cmat, Cmat2, Cmat3, 11 & nq, ipol, Ec, qwght, ldew, func) 12#endif 13c 14C$Id$ 15c 16 implicit none 17c 18#include "dft2drv.fh" 19#include "dft3drv.fh" 20c 21 double precision fac ! [input] 22 integer nq 23 integer ipol 24 double precision Ec 25 logical ldew 26 double precision func(*) ! value of the functional [output] 27c 28c Charge Density & Its Cube Root 29c 30 double precision rho(nq,(ipol*(ipol+1))/2) 31c 32c Charge Density Gradient 33c 34 double precision delrho(nq,3,ipol) 35c 36c Quadrature Weights 37c 38 double precision qwght(nq) 39c 40c Sampling Matrices for the XC Potential & Energy 41c 42 double precision Amat(nq,ipol), Cmat(nq,*) 43#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 44 double precision Amat2(nq,NCOL_AMAT2), Cmat2(nq,NCOL_CMAT2) 45#endif 46#ifdef THIRD_DERIV 47 double precision Amat3(nq,NCOL_AMAT3), Cmat3(nq,NCOL_CMAT3) 48#endif 49 double precision F13, F43, F113, F83, F53, F19, F79, P1, tol_rho, 50 & A, B, C, D 51c 52 Parameter (F13 = 1.D0/3.D0, F43 = 4.D0*F13, F113 = 11.D0*F13, 53 & F83 = 8.D0*F13, F53 = 5.D0*F13, F19 = 1.D0/9.D0, 54 & F79 = 7.D0*F19) 55#ifdef THIRD_DERIV 56 double precision F23, F73, F223 57 Parameter (F23 = 2.0D0*F13, F73 = 7.0D0*F13, F223 = 22.0*F13) 58#endif 59c 60c P1 = 2**(11/3)*(3/10)*(3*PI**2)**(2/3) 61c 62 Parameter (P1 = 0.3646239897876487D+02) 63c 64c Colle-Salvetti Empirical Parameters 65c 66 Parameter (A = 0.04918D0) 67 Parameter (B = 0.13200D0) 68 Parameter (C = 0.25330D0) 69 Parameter (D = 0.34900D0) 70c 71c Compute the partial derivatives of the correlation functional of 72c Lee, Yang and Parr. 73c 74c References: 75c 76c Colle & Salvetti, Theor. Chim. Acta 37, 329 (1975) 77c Lee, Yang & Parr, Phys. Rev. B 37, 785 (1988) 78c Miehlich, Savin, Stoll & Preuss, Chem. Phys. Lett. 157, 200 (1989) 79c Johnson, Gill & Pople, J. Chem. Phys. 98, 5612 (1993) 80c 81 integer n 82 double precision c1, c2, ab, ratcd 83 double precision rrho, rhoa, rhob, rrho2, rhoa2, rhob2, rhoab, 84 & rhoa53, rhob53, rhoa83, rhob83, rho2, 85 & h1, h2, h3, d1h1, d1h2, d1h3, om, d1om, de, d1de, de11, de47, 86 & t, d1t(2) 87c 88 double precision xrarb, d1xrarb(2) 89 double precision tm_in, d1tm_in(2) 90c 91 double precision gaa, gab, gbb 92 double precision f1, f2, d1f1(2), d1f2(2), f, d1f(5), 93 & d2fgaa(2), d2fgab(2), d2fgbb(2) 94c 95#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 96 double precision d2h1, d2h2, d2h3, d2om, d2de, rrho3, d2f1(3), 97 & d2f2(3), rhoa113, rhob113, d3fgaa(3), d3fgab(3), d3fgbb(3), 98 & d2t(3), d2f(3) 99 double precision dt 100 double precision d2xrarb(3) 101 double precision d2tm_in(3) 102#endif 103#ifdef THIRD_DERIV 104 double precision rrho4, rhoa23, d3h1, d3h2, d3h3, d3om, d3de, 105 1 d3f1(4), d3f2(4), d4fgaa(4), d4fgab(4), d4fgbb(4), d3t(4), 106 2 d3f(4) 107 double precision rhob23 108 double precision ddt 109 double precision d3xrarb(4) 110 double precision d3tm_in(4) 111#endif 112c 113c Coefficients of first two terms in LYP functional and other 114c commonly occurring factors 115c 116 c1 = -4.0d0*a 117 c2 = -P1*a*b 118 ab = a*b 119 ratcd = c/d 120c 121 if (ipol.eq.1)then 122c 123c ======> SPIN-RESTRICTED <====== 124c 125 do 10 n = 1, nq 126 if (rho(n,1).lt.tol_rho)goto 10 127c rrho = 1d0/rho(n,1) 128 rrho = 1.0d0/rho(n,1) 129 rhoa = 0.5d0*rho(n,1) 130 rrho2 = rrho*rrho 131c rho2 = 1d0/rrho2 132 rho2 = 1.0d0/rrho2 133 rhoa2 = rhoa*rhoa 134 rhoab = rhoa2 135c rhoa53 = abs(rhoa)**F53*sign(1d0,rhoa) 136 rhoa53 = abs(rhoa)**F53*sign(1.0d0,rhoa) 137c rhoa83 = abs(rhoa)**F83*sign(1d0,rhoa) 138 rhoa83 = abs(rhoa)**F83*sign(1.0d0,rhoa) 139#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 140 rrho3 = rrho*rrho2 141 rhoa113 = rhoa*rhoa83 142#endif 143c 144#ifdef THIRD_DERIV 145 rrho4 = rrho*rrho3 146 rhoa23 = abs(rhoa)**F23*sign(1.0d0,rhoa) 147#endif 148c 149 h2 = d*abs(rrho)**F13 150 d1h2 = -F13*h2*rrho 151#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 152 d2h2 = -F43*d1h2*rrho 153#endif 154c 155#ifdef THIRD_DERIV 156 d3h2 = -F73*d2h2*rrho 157#endif 158c 159 h3 = ratcd*h2 160 d1h3 = ratcd*d1h2 161#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 162 d2h3 = ratcd*d2h2 163#endif 164c 165#ifdef THIRD_DERIV 166 d3h3 = ratcd*d3h2 167#endif 168c 169c h1 = 1d0/(1d0+h2) 170 h1 = 1.0d0/(1.0d0+h2) 171 d1h1 = -h1*h1*d1h2 172#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 173c d2h1 = -(2d0*h1*d1h1*d1h2 + h1*h1*d2h2) 174 d2h1 = -(2.0d0*h1*d1h1*d1h2 + h1*h1*d2h2) 175#endif 176c 177#ifdef THIRD_DERIV 178 d3h1 = -6.0d0*d1h1*d1h1*d1h2 179 1 - 6.0d0*h1*d2h2*d1h1 180 2 - h1*h1*d3h2 181#endif 182c 183! om = exp(-h3)*h1*rrho**F113 184 om = exp(-h3)*h1*abs(rrho)**F113 185 t = d1h3+h1*d1h2+F113*rrho 186 d1om = -om*t 187#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 188 dt = d2h3+d1h1*d1h2+h1*d2h2-F113*rrho2 189c d2om = -(d1om*t+om*(d2h3+d1h1*d1h2+h1*d2h2-F113*rrho2)) 190 d2om = -(d1om*t+om*dt) 191#endif 192c 193#ifdef THIRD_DERIV 194 ddt = d3h3 + d2h1*d1h2 + 2.0d0*d1h1*d2h2 195 1 + h1*d3h2 + F223*rrho3 196 d3om = -(ddt*om + 2.0d0*d1om*dt + d2om*t) 197#endif 198c 199 de = h3+h1*h2 200 d1de = d1h3 + d1h1*h2 + h1*d1h2 201#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 202c d2de = d2h3 + d2h1*h2 + h1*d2h2 + 2d0*d1h1*d1h2 203 d2de = d2h3 + d2h1*h2 + h1*d2h2 + 2.0d0*d1h1*d1h2 204#endif 205c 206#ifdef THIRD_DERIV 207 d3de = d3h3 + d3h1*h2 + 3.0d0*d2h1*d1h2 208 1 + 3.0d0*d1h1*d2h2 + h1*d3h2 209#endif 210c 211 f1 = h1*rhoab*rrho 212 d1f1(1) = d1h1*rhoab*rrho - h1*rhoab*rrho2 213 d1f1(1) = d1f1(1) + h1*rhoa*rrho 214#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 215c Daniel (10-30-12): I'd suggest rewriting the second derivatives of 216c term 1 in a much simpler format. This relies on the fact that 217c 2*rhoa (or 2*rhob) is rho for a restricted calculation. 218c d2f1(1) = d2h1*rhoab*rrho + 2d0*d1h1*(rhoa*rrho-rhoab*rrho2) 219c & + 2d0*h1*(-rhoa*rrho2+rhoab*rrho3) 220c d2f1(2) = d2h1*rhoab*rrho + d1h1*(1d0-2d0*rhoab*rrho2) 221c & + 2d0*h1*rhoab*rrho3 222 d2f1(1) = d2h1*rhoab*rrho 223 1 + d1h1*(1.0d0 - rhoa*rrho) 224 2 + h1*(-rrho + rhoa*rrho2) 225 d2f1(2) = d2h1*rhoab*rrho 226 1 + d1h1*(1.0d0-rhoa*rrho) 227 2 + h1*rhoa*rrho2 228#endif 229c 230#ifdef THIRD_DERIV 231 d3f1(1) = d3h1*rhoab*rrho 232 1 + d2h1*( 1.0d0 - rhoab*rrho2) 233 2 + d1h1*(-1.5d0*rrho) 234 3 + h1*(1.5d0*rrho2) 235 d3f1(2) = d3h1*rhoab*rrho 236 1 + d2h1*(1.0d0 - rhoab*rrho2) 237 2 + d1h1*(0.50d0*rrho) 238 3 + h1*(-0.50d0*rrho2) 239 d3f1(3) = d3f1(2) 240 d3f1(4) = d3f1(1) 241#endif 242c 243 xrarb = rhoab*(rhoa83+rhoa83) 244 f2 = om*xrarb 245c 246 d1xrarb(1) = rhoa*(F113*rhoa83+rhoa83) 247c 248c d1f2(1) = d1om*rhoab*(rhoa83+rhoa83) 249c d1f2(1) = d1f2(1) + om*rhoa*(F113*rhoa83+rhoa83) 250 d1f2(1) = d1om*xrarb + om*d1xrarb(1) 251#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 252c d2f2(1) = d2om*rhoab*(rhoa83+rhoa83) 253c & + 2d0*d1om*rhoa*(F113*rhoa83+rhoa83) 254c & + om*rhoa*F113*F83*rhoa53 255c d2f2(2) = d2om*rhoab*(rhoa83+rhoa83) 256c & + d1om*(rhoa113+rhoa113+F113*(rhoa*rhoa83+rhoa*rhoa83)) 257c & + om*F113*(rhoa83+rhoa83) 258c 259 d2xrarb(1) = rhoa*F113*F83*rhoa53 260 d2xrarb(2) = F113*(rhoa83+rhoa83) 261c 262 d2f2(1) = d2om*xrarb 263 & + 2.0d0*d1om*d1xrarb(1) 264 & + om*d2xrarb(1) 265 d2f2(2) = d2om*xrarb 266 & + 2.0d0*d1om*d1xrarb(1) 267 & + om*d2xrarb(2) 268#endif 269c 270#ifdef THIRD_DERIV 271 d3xrarb(1) = rhoa*F113*F83*F53*rhoa23 272 d3xrarb(2) = F113*F83*rhoa53 273 d3xrarb(3) = F113*F83*rhoa53 274c 275 d3f2(1) = d3om*xrarb 276 1 + 3.0d0*d2om*d1xrarb(1) 277 2 + 3.0d0*d1om*d2xrarb(1) 278 3 + om*d3xrarb(1) 279 d3f2(2) = d3om*xrarb 280 1 + 3.0d0*d2om*d1xrarb(1) 281 2 + d1om*(d2xrarb(1) + 2.0d0*d2xrarb(2)) 282 3 + om*d3xrarb(2) 283 d3f2(3) = d3f2(2) 284 d3f2(4) = d3f2(1) 285#endif 286c 287 gaa =(delrho(n,1,1)*delrho(n,1,1) + 288 & delrho(n,2,1)*delrho(n,2,1) + 289 & delrho(n,3,1)*delrho(n,3,1))*0.25d0 290c 291c de11 = de - 11d0 292c de47 = 47d0 - 7d0*de 293 de11 = de - 11.0d0 294 de47 = 47.0d0 - 7.0d0*de 295c 296c Daniel (10-30-12): tm_in is what I call Qi (the inside term) 297 tm_in = 1.0d0 - 3.0d0*de - de11*rhoa*rrho 298c Daniel (10-23-12): "t" is what I call Q or S. 299c t = F19*rhoab*(1d0-3d0*de-de11*rhoa*rrho)-rhoa2 300 t = F19*rhoab*tm_in - rhoa2 301c Daniel (10-30-12): Derivatives of the inside term 302 d1tm_in(1) = -(3.0d0+rhoa*rrho)*d1de-de11*rhoa*rrho2 303 d1tm_in(2) = -(3.0d0+rhoa*rrho)*d1de+de11*rhoa*rrho2 304c Daniel (10-23-12): d1t(1) is the derivative with respect to rhoa, 305c and d1t(2) is the derivative with respect to rhob. 306c d1t(1) = F19*(rhoa*(1d0-3d0*de-de11*rhoa*rrho) 307c & - rhoab*((3d0+rhoa*rrho)*d1de+de11*rhoa*rrho2)) 308c d1t(2) = F19*(rhoa*(1d0-3d0*de-de11*rhoa*rrho) 309c & + rhoab*(-(3d0+rhoa*rrho)*d1de+de11*rhoa*rrho2)) 310c & - 2d0*rhoa 311 d1t(1) = F19*( rhoa*tm_in + rhoab*d1tm_in(1) ) 312 d1t(2) = F19*( rhoa*tm_in + rhoab*d1tm_in(2) ) 313 & - 2.0d0*rhoa 314c 315 d1f(3) = -ab*om*t 316c 317 d2fgaa(1) = -ab*( d1om*t + om*d1t(1) ) 318 d2fgaa(2) = -ab*( d1om*t + om*d1t(2) ) 319#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 320c Daniel (10-30-12): Derivatives of the inside term, Qi. 321 d2tm_in(1) = -(3.0d0+rhoa*rrho)*d2de 322 1 - 2.0d0*d1de*rhoa*rrho2 323 2 + 2.0d0*de11*rhoa*rrho3 324 d2tm_in(2) = -(3.0d0+rhoa*rrho) ! Written without d2de 325 d2tm_in(3) = -(3.0d0+rhoa*rrho)*d2de 326 1 + 2.0d0*d1de*rhoa*rrho2 327 2 - 2.0d0*de11*rhoa*rrho3 328c 329c d2t(1) = -F19*( 330c & 2d0*rhoa*((3d0+rhoa*rrho)*d1de+de11*rhoa*rrho2) 331c & + rhoab*((3d0+rhoa*rrho)*d2de+2d0*d1de*rhoa*rrho2 332c & -2d0*de11*rhoa*rrho3)) 333c d2t(2) = F19*(1d0-3d0*de-de11*rhoa*rrho 334c & - rho(n,1)*(3d0+rhoa*rrho)*d1de 335c & - rhoab*((3d0+rhoa*rrho)*d2de)) 336c d2t(3) = -F19*( 337c & 2d0*rhoa*((3d0+rhoa*rrho)*d1de-de11*rhoa*rrho2) 338c & + rhoab*((3d0+rhoa*rrho)*d2de-2d0*d1de*rhoa*rrho2 339c & +2d0*de11*rhoa*rrho3)) 340c & - 2d0 341 d2t(1) = F19*( 2.0d0*rhoa*d1tm_in(1) 342 1 + rhoab*d2tm_in(1) ) 343 d2t(2) = F19*( tm_in 344 1 + rho(n,1)*d2tm_in(2)*d1de 345 2 + rhoab*d2tm_in(2)*d2de ) 346 d2t(3) = F19*( 2.0d0*rhoa*d1tm_in(2) 347 1 + rhoab*d2tm_in(3) ) 348 2 - 2.0d0 349c 350c d3fgaa(1) = -ab*(d2om*t+2d0*d1om*d1t(1)+om*d2t(1)) 351c d3fgaa(2) = -ab*(d2om*t+d1om*(d1t(1)+d1t(2))+om*d2t(2)) 352c d3fgaa(3) = -ab*(d2om*t+2d0*d1om*d1t(2)+om*d2t(3)) 353 d3fgaa(1) = -ab*( d2om*t + 2.0d0*d1om*d1t(1) + om*d2t(1) ) 354 d3fgaa(2) = -ab*( d2om*t + d1om*( d1t(1) + d1t(2) ) 355 1 + om*d2t(2) ) 356 d3fgaa(3) = -ab*( d2om*t + 2.0d0*d1om*d1t(2) + om*d2t(3) ) 357#endif 358c 359#ifdef THIRD_DERIV 360 d3tm_in(1) = -( 3.0d0 + rhoa*rrho )*d3de 361 1 - 3.0d0*d2de*rhoa*rrho2 362 2 + 6.0d0*d1de*rhoa*rrho3 363 3 - 6.0d0*de11*rhoa*rrho4 364 d3tm_in(2) = -( 3.0d0 + rhoa*rrho )*d3de 365 1 - d2de*rhoa*rrho2 366 2 + 2.0d0*d1de*rhoa*rrho3 367 3 - 2.0d0*de11*rhoa*rrho4 368 d3tm_in(3) = -( 3.0d0 + rhoa*rrho )*d3de 369 1 + d2de*rhoa*rrho2 370 2 - 2.0d0*d1de*rhoa*rrho3 371 3 + 2.0d0*de11*rhoa*rrho4 372 d3tm_in(4) = -( 3.0d0 + rhoa*rrho )*d3de 373 1 + 3.0d0*d2de*rhoa*rrho2 374 2 - 6.0d0*d1de*rhoa*rrho3 375 3 + 6.0d0*de11*rhoa*rrho4 376c 377 d3t(1) = F19*( 3.0d0*rhoa*d2tm_in(1) 378 1 + rhoab*d3tm_in(1) ) 379 d3t(2) = F19*( 2.0d0*d1tm_in(1) 380 1 + rhoa*d2tm_in(1) 381 2 + 2.0d0*rhoa*d2tm_in(2)*d2de 382 3 + rhoab*d3tm_in(2) ) 383 d3t(3) = F19*( 2.0d0*d1tm_in(2) 384 1 + 2.0d0*rhoa*d2tm_in(2)*d2de 385 2 + rhoa*d2tm_in(3) 386 3 + rhoab*d3tm_in(3) ) 387 d3t(4) = F19*( 3.0d0*rhoa*d2tm_in(3) 388 1 + rhoab*d3tm_in(4) ) 389c 390 d4fgaa(1) = -ab*( d3om*t + 3.0d0*d2om*d1t(1) 391 1 + 3.0d0*d1om*d2t(1) + om*d3t(1) ) 392 d4fgaa(2) = -ab*( d3om*t + d2om*(2.0d0*d1t(1) + d1t(2)) 393 1 + d1om*(d2t(1) + 2.0d0*d2t(2)) + om*d3t(2) ) 394 d4fgaa(3) = -ab*( d3om*t + d2om*(d1t(1) + 2.0d0*d1t(2)) 395 1 + d1om*(2.0d0*d2t(2) + d2t(3)) + om*d3t(3) ) 396 d4fgaa(4) = -ab*( d3om*t + 3.0d0*d2om*d1t(2) 397 1 + 3.0d0*d1om*d2t(3) + om*d3t(4) ) 398#endif 399c 400c Daniel (10-23-12): "t" is what I call R. 401 t = F19*rhoab*de47-F43*rho2 402c 403 d1t(1) = F19*rhoa*de47 - F79*rhoab*d1de - F83*rho(n,1) 404c 405 d1f(4) = -ab*om*t 406c 407 d2fgab(1) = -ab*( d1om*t + om*d1t(1) ) 408#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 409c d2t(1) = -F79*(2d0*rhoa*d1de+rhoab*d2de) - F83 410c d2t(2) = F19*de47 - F79*(rho(n,1)*d1de+rhoab*d2de) - F83 411 d2t(1) = -F79*( 2.0d0*rhoa*d1de + rhoab*d2de ) - F83 412 d2t(2) = F19*de47 - F79*( rho(n,1)*d1de + rhoab*d2de ) 413 1 - F83 414c 415c d3fgab(1) = -ab*(d2om*t+2d0*d1om*d1t(1)+om*d2t(1)) 416c d3fgab(2) = -ab*(d2om*t+2d0*d1om*d1t(1)+om*d2t(2)) 417 d3fgab(1) = -ab*( d2om*t + 2.0d0*d1om*d1t(1) + om*d2t(1) ) 418 d3fgab(2) = -ab*( d2om*t + 2.0d0*d1om*d1t(1) + om*d2t(2) ) 419c To keep compilers quiet (WE MAY NEED TO FIX THIS) 420 d3fgab(3) = d3fgab(1) 421#endif 422c 423#ifdef THIRD_DERIV 424 d3t(1) = -F79*( 3.0d0*rhoa*d2de + rhoab*d3de ) 425 d3t(2) = -F79*( 2.0d0*d1de + 3.0d0*rhoa*d2de 426 1 + rhoab*d3de ) 427 d3t(3) = d3t(2) 428 d3t(4) = d3t(1) 429c 430 d4fgab(1) = -ab*( d3om*t + 3.0d0*d2om*d1t(1) 431 1 + 3.0d0*d1om*d2t(1) + om*d3t(1) ) 432 d4fgab(2) = -ab*( d3om*t + 3.0d0*d2om*d1t(1) 433 1 + d1om*(d2t(1) + 2.0d0*d2t(2)) 434 2 + om*d3t(2) ) 435 d4fgab(3) = -ab*( d3om*t + 3.0d0*d2om*d1t(1) 436 1 + d1om*(d2t(1) + 2.0d0*d2t(2)) 437 2 + om*d3t(3) ) 438 d4fgab(4) = -ab*( d3om*t + 3.0d0*d2om*d1t(1) 439 1 + 3.0d0*d1om*d2t(1) + om*d3t(4) ) 440#endif 441c 442 d2fgbb(1) = d2fgaa(2) 443#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 444 d3fgbb(1) = d3fgaa(3) 445 d3fgbb(2) = d3fgaa(2) 446 d3fgbb(3) = d3fgaa(1) 447#endif 448c 449#ifdef THIRD_DERIV 450 d4fgbb(1) = d4fgaa(4) 451 d4fgbb(2) = d4fgaa(3) 452 d4fgbb(3) = d4fgaa(2) 453 d4fgbb(4) = d4fgaa(1) 454#endif 455c 456 f = c1*f1 + c2*f2 + gaa*(2d0*d1f(3) + d1f(4)) 457 d1f(1) = c1*d1f1(1) + c2*d1f2(1) 458 & + gaa*(d2fgaa(1) + d2fgab(1) + d2fgbb(1)) 459#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 460 d2f(1) = c1*d2f1(1) + c2*d2f2(1) 461 & + gaa*(d3fgaa(1) + d3fgab(1) + d3fgbb(1)) 462 d2f(2) = c1*d2f1(2) + c2*d2f2(2) 463 & + gaa*(d3fgaa(2) + d3fgab(2) + d3fgbb(2)) 464#endif 465c 466#ifdef THIRD_DERIV 467 d3f(1) = c1*d3f1(1) + c2*d3f2(1) 468 1 + gaa*(d4fgaa(1) + d4fgab(1) + d4fgbb(1)) 469 d3f(2) = c1*d3f1(2) + c2*d3f2(2) 470 1 + gaa*(d4fgaa(2) + d4fgab(2) + d4fgbb(2)) 471 d3f(3) = c1*d3f1(3) + c2*d3f2(3) 472 1 + gaa*(d4fgaa(3) + d4fgab(3) + d4fgbb(3)) 473 d3f(4) = c1*d3f1(4) + c2*d3f2(4) 474 1 + gaa*(d4fgaa(4) + d4fgab(4) + d4fgbb(4)) 475#endif 476c 477 Ec = Ec + f*fac*qwght(n) 478 if (ldew) func(n) = func(n) + f*fac 479 Amat(n,1) = Amat(n,1) + d1f(1)*fac 480 Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + d1f(3)*fac 481 Cmat(n,D1_GAB) = Cmat(n,D1_GAB) + d1f(4)*fac 482#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 483 Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) + d2f(1)*fac 484 Amat2(n,D2_RA_RB) = Amat2(n,D2_RA_RB) + d2f(2)*fac 485 Cmat2(n,D2_RA_GAA) = Cmat2(n,D2_RA_GAA) + d2fgaa(1)*fac 486 Cmat2(n,D2_RA_GAB) = Cmat2(n,D2_RA_GAB) + d2fgab(1)*fac 487 Cmat2(n,D2_RA_GBB) = Cmat2(n,D2_RA_GBB) + d2fgbb(1)*fac 488#endif 489c 490#ifdef THIRD_DERIV 491 Amat3(n,D3_RA_RA_RA) = Amat3(n,D3_RA_RA_RA) + d3f(1)*fac 492 Amat3(n,D3_RA_RA_RB) = Amat3(n,D3_RA_RA_RB) + d3f(2)*fac 493 Amat3(n,D3_RA_RB_RB) = Amat3(n,D3_RA_RB_RB) + d3f(3)*fac 494c 495 Cmat3(n,D3_RA_RA_GAA) = Cmat3(n,D3_RA_RA_GAA) 496 1 + d3fgaa(1)*fac 497 Cmat3(n,D3_RA_RA_GAB) = Cmat3(n,D3_RA_RA_GAB) 498 1 + d3fgab(1)*fac 499 Cmat3(n,D3_RA_RA_GBB) = Cmat3(n,D3_RA_RA_GBB) 500 1 + d3fgbb(1)*fac 501c 502 Cmat3(n,D3_RA_RB_GAA) = Cmat3(n,D3_RA_RB_GAA) 503 1 + d3fgaa(2)*fac 504 Cmat3(n,D3_RA_RB_GAB) = Cmat3(n,D3_RA_RB_GAB) 505 1 + d3fgab(2)*fac 506 Cmat3(n,D3_RA_RB_GBB) = Cmat3(n,D3_RA_RB_GBB) 507 1 + d3fgbb(2)*fac 508#endif 509c 510 10 continue 511 else 512c 513c ======> SPIN-UNRESTRICTED <====== 514c 515 do 20 n = 1,nq 516 if (rho(n,1).lt.tol_rho)goto 20 517 rrho = 1d0/rho(n,1) 518 rhoa = max(0.0d0,rho(n,2)) 519 rhob = max(0.0d0,rho(n,3)) 520 rrho2 = rrho*rrho 521 rho2 = 1d0/rrho2 522 rhoa2 = rhoa*rhoa 523 rhob2 = rhob*rhob 524 rhoab = rhoa*rhob 525 rhoa53 = rhoa**F53 526 rhob53 = rhob**F53 527csedo rhoa53 = abs(rhoa)**F53 528csedo rhob53 = abs(rhob)**F53 529 rhoa83 = rhoa*rhoa53 530 rhob83 = rhob*rhob53 531#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 532 rrho3 = rrho*rrho2 533 rhoa113 = rhoa*rhoa83 534 rhob113 = rhob*rhob83 535#endif 536c 537#ifdef THIRD_DERIV 538 rrho4 = rrho*rrho3 539 rhoa23 = rhoa**F23 540 rhob23 = rhob**F23 541#endif 542c 543cedo h2 = d*abs(rrho)**F13*sign(1d0,rrho) 544 h2 = d*rrho**F13 545 d1h2 = -F13*h2*rrho 546#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 547 d2h2 = -F43*d1h2*rrho 548#endif 549c 550#ifdef THIRD_DERIV 551 d3h2 = -F73*d2h2*rrho 552#endif 553c 554 h3 = ratcd*h2 555 d1h3 = ratcd*d1h2 556#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 557 d2h3 = ratcd*d2h2 558#endif 559c 560#ifdef THIRD_DERIV 561 d3h3 = ratcd*d3h2 562#endif 563c 564 h1 = 1d0/(1d0+h2) 565 d1h1 = -h1*h1*d1h2 566#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 567 d2h1 = -(2d0*h1*d1h1*d1h2 + h1*h1*d2h2) 568#endif 569c 570#ifdef THIRD_DERIV 571 d3h1 = -6.0d0*d1h1*d1h1*d1h2 572 1 - 6.0d0*h1*d2h2*d1h1 573 2 - h1*h1*d3h2 574#endif 575c 576 om = exp(-h3)*h1*rrho**F113 577cedo om = exp(-h3)*h1*abs(rrho)**F113*sign(1d0,rrho) 578 t = d1h3+h1*d1h2+F113*rrho 579 d1om = -om*t 580#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 581 dt = d2h3+d1h1*d1h2+h1*d2h2-F113*rrho2 582c d2om = -(d1om*t+om*(d2h3+d1h1*d1h2+h1*d2h2-F113*rrho2)) 583 d2om = -(d1om*t+om*dt) 584#endif 585c 586#ifdef THIRD_DERIV 587 ddt = d3h3 + d2h1*d1h2 + 2.0d0*d1h1*d2h2 588 1 + h1*d3h2 + F223*rrho3 589 d3om = -(ddt*om + 2.0d0*d1om*dt + d2om*t) 590#endif 591c 592 de = h3+h1*h2 593 d1de = d1h3 + d1h1*h2 + h1*d1h2 594#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 595 d2de = d2h3 + d2h1*h2 + h1*d2h2 + 2d0*d1h1*d1h2 596#endif 597c 598#ifdef THIRD_DERIV 599 d3de = d3h3 + d3h1*h2 + 3.0d0*d2h1*d1h2 600 1 + 3.0d0*d1h1*d2h2 + h1*d3h2 601#endif 602c 603c Daniel (3-21-13): f1 is J/(-4*a) in my notes. 604 f1 = h1*rhoab*rrho 605 d1f1(1) = d1h1*rhoab*rrho - h1*rhoab*rrho2 606 d1f1(2) = d1f1(1) 607 d1f1(1) = d1f1(1) + h1*rhob*rrho 608 d1f1(2) = d1f1(2) + h1*rhoa*rrho 609#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 610 d2f1(1) = d2h1*rhoab*rrho + 2d0*d1h1*(rhob*rrho-rhoab*rrho2) 611 & + 2d0*h1*(-rhob*rrho2+rhoab*rrho3) 612 d2f1(2) = d2h1*rhoab*rrho + d1h1*(1d0-2d0*rhoab*rrho2) 613 & + 2d0*h1*rhoab*rrho3 614 d2f1(3) = d2h1*rhoab*rrho + 2d0*d1h1*(rhoa*rrho-rhoab*rrho2) 615 & + 2d0*h1*(-rhoa*rrho2+rhoab*rrho3) 616#endif 617c 618#ifdef THIRD_DERIV 619 d3f1(1) = d3h1*rhoab*rrho 620 1 + 3.0d0*d2h1*(rhob*rrho-rhoab*rrho2) 621 2 + 6.0d0*d1h1*(-rhob*rrho2+rhoab*rrho3) 622 3 + 6.0d0*h1*(rhob*rrho3-rhoab*rrho4) 623 d3f1(2) = d3h1*rhoab*rrho 624 1 + d2h1*( (2.0d0*rhob + rhoa)*rrho 625 2 - 3.0d0*rhoab*rrho2 ) 626 3 + 2.0d0*d1h1*( rrho 627 4 - (2.0d0*rhob + rhoa)*rrho2 628 5 + 3.0d0*rhoab*rrho3 ) 629 6 + 2.0d0*h1*( -rrho2 630 7 + (2.0d0*rhob + rhoa)*rrho3 631 8 - 3.0d0*rhoab*rrho4 ) 632 d3f1(3) = d3h1*rhoab*rrho 633 1 + d2h1*( (rhob + 2.0d0*rhoa)*rrho 634 2 - 3.0d0*rhoab*rrho2 ) 635 3 + 2.0d0*d1h1*( rrho 636 4 - (rhob + 2.0d0*rhoa)*rrho2 637 5 + 3.0d0*rhoab*rrho3 ) 638 6 + 2.0d0*h1*( -rrho2 639 7 + (rhob + 2.0d0*rhoa)*rrho3 640 8 - 3.0d0*rhoab*rrho4 ) 641 d3f1(4) = d3h1*rhoab*rrho 642 1 + 3.0d0*d2h1*(rhoa*rrho-rhoab*rrho2) 643 2 + 6.0d0*d1h1*(-rhoa*rrho2+rhoab*rrho3) 644 3 + 6.0d0*h1*(rhoa*rrho3-rhoab*rrho4) 645#endif 646c 647c Daniel (10-30-12): Define xrarb here 648 xrarb = rhoab*(rhoa83+rhob83) 649c 650c f2 =om*rhoab*(rhoa83+rhob83) 651 f2 =om*xrarb 652c 653 d1xrarb(1) = rhob*(F113*rhoa83+rhob83) 654 d1xrarb(2) = rhoa*(F113*rhob83+rhoa83) 655c 656c d1f2(1) = d1om*rhoab*(rhoa83+rhob83) 657 d1f2(1) = d1om*xrarb 658 d1f2(2) = d1f2(1) 659c d1f2(1) = d1f2(1) + om*rhob*(F113*rhoa83+rhob83) 660 d1f2(1) = d1f2(1) + om*d1xrarb(1) 661c d1f2(2) = d1f2(2) + om*rhoa*(F113*rhob83+rhoa83) 662 d1f2(2) = d1f2(2) + om*d1xrarb(2) 663#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 664 d2xrarb(1) = rhob*F113*F83*rhoa53 665 d2xrarb(2) = F113*(rhoa83+rhob83) 666 d2xrarb(3) = rhoa*F113*F83*rhob53 667c 668c d2f2(1) = d2om*rhoab*(rhoa83+rhob83) 669c & + 2d0*d1om*rhob*(F113*rhoa83+rhob83) 670c & + om*rhob*F113*F83*rhoa53 671c d2f2(2) = d2om*rhoab*(rhoa83+rhob83) 672c & + d1om*(rhoa113+rhob113+F113*(rhob*rhoa83+rhoa*rhob83)) 673c & + om*F113*(rhoa83+rhob83) 674c d2f2(3) = d2om*rhoab*(rhoa83+rhob83) 675c & + 2d0*d1om*rhoa*(F113*rhob83+rhoa83) 676c & + om*rhoa*F113*F83*rhob53 677 d2f2(1) = d2om*xrarb 678 & + 2d0*d1om*d1xrarb(1) 679 & + om*d2xrarb(1) 680 d2f2(2) = d2om*xrarb 681 & + d1om*(d1xrarb(1) + d1xrarb(2)) 682 & + om*d2xrarb(2) 683 d2f2(3) = d2om*xrarb 684 & + 2d0*d1om*d1xrarb(2) 685 & + om*d2xrarb(3) 686#endif 687c 688#ifdef THIRD_DERIV 689 d3xrarb(1) = rhob*F113*F83*F53*rhoa23 690 d3xrarb(2) = F113*F83*rhoa53 691 d3xrarb(3) = F113*F83*rhob53 692 d3xrarb(4) = rhoa*F113*F83*F53*rhob23 693c 694 d3f2(1) = d3om*xrarb 695 1 + 3.0d0*d2om*d1xrarb(1) 696 2 + 3.0d0*d1om*d2xrarb(1) 697 3 + om*d3xrarb(1) 698 d3f2(2) = d3om*xrarb 699 1 + d2om*( d1xrarb(2) + 2.0d0*d1xrarb(1) ) 700 2 + d1om*( 2.0d0*d2xrarb(2) + d2xrarb(1) ) 701 3 + om*d3xrarb(2) 702 d3f2(3) = d3om*xrarb 703 1 + d2om*( 2.0d0*d1xrarb(2) + d1xrarb(1) ) 704 2 + d1om*( d2xrarb(3) + 2.0d0*d2xrarb(2) ) 705 3 + om*d3xrarb(3) 706 d3f2(4) = d3om*xrarb 707 1 + 3.0d0*d2om*d1xrarb(2) 708 2 + 3.0d0*d1om*d2xrarb(3) 709 3 + om*d3xrarb(4) 710#endif 711c 712 gaa = delrho(n,1,1)*delrho(n,1,1) + 713 & delrho(n,2,1)*delrho(n,2,1) + 714 & delrho(n,3,1)*delrho(n,3,1) 715 gab = delrho(n,1,1)*delrho(n,1,2) + 716 & delrho(n,2,1)*delrho(n,2,2) + 717 & delrho(n,3,1)*delrho(n,3,2) 718 gbb = delrho(n,1,2)*delrho(n,1,2) + 719 & delrho(n,2,2)*delrho(n,2,2) + 720 & delrho(n,3,2)*delrho(n,3,2) 721c 722 de11 = de - 11d0 723 de47 = 47d0 - 7d0*de 724c Daniel (10-30-12): tm_in is what I call Qi (the inside term) 725 tm_in = 1.0d0 - 3.0d0*de - de11*rhoa*rrho 726c 727c t = F19*rhoab*(1d0-3d0*de-de11*rhoa*rrho)-rhob2 728 t = F19*rhoab*tm_in-rhob2 729c Daniel (10-30-12): Derivatives of the inside term 730 d1tm_in(1) = -(3.0d0+rhoa*rrho)*d1de-de11*rhob*rrho2 731 d1tm_in(2) = -(3.0d0+rhoa*rrho)*d1de+de11*rhoa*rrho2 732c 733c d1t(1) = F19*(rhob*(1d0-3d0*de-de11*rhoa*rrho) 734c & - rhoab*((3d0+rhoa*rrho)*d1de+de11*rhob*rrho2)) 735c d1t(2) = F19*(rhoa*(1d0-3d0*de-de11*rhoa*rrho) 736c & + rhoab*(-(3d0+rhoa*rrho)*d1de+de11*rhoa*rrho2)) 737c & - 2d0*rhob 738 d1t(1) = F19*(rhob*tm_in + rhoab*d1tm_in(1)) 739 d1t(2) = F19*(rhoa*tm_in + rhoab*d1tm_in(2)) - 2.0d0*rhob 740c 741 d1f(3) = -ab*om*t 742c 743 d2fgaa(1) = -ab*( d1om*t + om*d1t(1) ) 744 d2fgaa(2) = -ab*( d1om*t + om*d1t(2) ) 745#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 746c Daniel (10-30-12): Derivatives of the inside term, Qi. 747 d2tm_in(1) = -(3.0d0+rhoa*rrho)*d2de 748 1 - 2.0d0*d1de*rhob*rrho2 749 2 + 2.0d0*de11*rhob*rrho3 750 d2tm_in(2) = -(3.0d0+rhoa*rrho)*d2de 751 1 + ( rhoa - rhob)*d1de*rrho2 752 2 + ( rhob - rhoa)*de11*rrho3 753 d2tm_in(3) = -(3.0d0+rhoa*rrho)*d2de 754 1 + 2.0d0*d1de*rhoa*rrho2 755 2 - 2.0d0*de11*rhoa*rrho3 756c 757c d2t(1) = -F19*( 758c & 2d0*rhob*((3d0+rhoa*rrho)*d1de+de11*rhob*rrho2) 759c & + rhoab*((3d0+rhoa*rrho)*d2de+2d0*d1de*rhob*rrho2 760c & -2d0*de11*rhob*rrho3)) 761c d2t(2) = F19*(1d0-3d0*de-de11*rhoa*rrho 762c & - rho(n,1)*(3d0+rhoa*rrho)*d1de 763c & - rhoab*((3d0+rhoa*rrho)*d2de-d1de*(rhoa-rhob)*rrho2 764c & +de11*(rhoa-rhob)*rrho3)) 765c d2t(3) = -F19*( 766c & 2d0*rhoa*((3d0+rhoa*rrho)*d1de-de11*rhoa*rrho2) 767c & + rhoab*((3d0+rhoa*rrho)*d2de-2d0*d1de*rhoa*rrho2 768c & +2d0*de11*rhoa*rrho3)) 769c & - 2d0 770 d2t(1) = F19*( 2.0d0*rhob*d1tm_in(1) 771 & + rhoab*d2tm_in(1) ) 772 d2t(2) = F19*( tm_in 773 & - rho(n,1)*(3.0d0+rhoa*rrho)*d1de 774 & + rhoab*d2tm_in(2) ) 775 d2t(3) = F19*( 2.0d0*rhoa*d1tm_in(2) 776 & + rhoab*d2tm_in(3) ) 777 & - 2.0d0 778c 779 d3fgaa(1) = -ab*( d2om*t + 2d0*d1om*d1t(1) + om*d2t(1) ) 780 d3fgaa(2) = -ab*( d2om*t + d1om*( d1t(1) + d1t(2) ) 781 1 + om*d2t(2) ) 782 d3fgaa(3) = -ab*( d2om*t + 2d0*d1om*d1t(2) + om*d2t(3) ) 783#endif 784c 785#ifdef THIRD_DERIV 786 d3tm_in(1) = -( 3.0d0 + rhoa*rrho )*d3de 787 1 - 3.0d0*d2de*rhob*rrho2 788 2 + 6.0d0*d1de*rhob*rrho3 789 3 - 6.0d0*de11*rhob*rrho4 790 d3tm_in(2) = -( 3.0d0 + rhoa*rrho )*d3de 791 1 + d2de*rhoa*rrho2 792 2 - 2.0d0*d2de*rhob*rrho2 793 3 - 2.0d0*d1de*rhoa*rrho3 794 4 + 4.0d0*d1de*rhob*rrho3 795 5 - 2.0d0*de11*(2.0d0*rhob-rhoa)*rrho4 796 d3tm_in(3) = -( 3.0d0 + rhoa*rrho )*d3de 797 1 + 2.0d0*d2de*rhoa*rrho2 798 2 - d2de*rhob*rrho2 799 3 - 4.0d0*d1de*rhoa*rrho3 800 4 + 2.0d0*d1de*rhob*rrho3 801 5 - 2.0d0*de11*(rhob-2.0d0*rhoa)*rrho4 802 d3tm_in(4) = -( 3.0d0 + rhoa*rrho )*d3de 803 1 + 3.0d0*d2de*rhoa*rrho2 804 2 - 6.0d0*d1de*rhoa*rrho3 805 3 + 6.0d0*de11*rhoa*rrho4 806c 807 d3t(1) = F19*( 3.0d0*rhob*d2tm_in(1) 808 1 + rhoab*d3tm_in(1) ) 809 d3t(2) = F19*( 2.0d0*d1tm_in(1) 810 1 + rhoa*d2tm_in(1) 811 2 + 2.0d0*rhob*d2tm_in(2) 812 3 + rhoab*d3tm_in(2) ) 813 d3t(3) = F19*( 2.0d0*d1tm_in(2) 814 1 + 2.0d0*rhoa*d2tm_in(2) 815 2 + rhob*d2tm_in(3) 816 3 + rhoab*d3tm_in(3) ) 817 d3t(4) = F19*( 3.0d0*rhoa*d2tm_in(3) 818 1 + rhoab*d3tm_in(4) ) 819c 820 d4fgaa(1) = -ab*( d3om*t + 3.0d0*d2om*d1t(1) 821 1 + 3.0d0*d1om*d2t(1) + om*d3t(1) ) 822 d4fgaa(2) = -ab*( d3om*t + d2om*(2.0d0*d1t(1) + d1t(2)) 823 1 + d1om*(d2t(1) + 2.0d0*d2t(2)) + om*d3t(2) ) 824 d4fgaa(3) = -ab*( d3om*t + d2om*(d1t(1) + 2.0d0*d1t(2)) 825 1 + d1om*(2.0d0*d2t(2) + d2t(3)) + om*d3t(3) ) 826 d4fgaa(4) = -ab*( d3om*t + 3.0d0*d2om*d1t(2) 827 1 + 3.0d0*d1om*d2t(3) + om*d3t(4) ) 828#endif 829c 830 t = F19*rhoab*de47 - F43*rho2 831 d1t(1) = F19*rhob*de47 - F79*rhoab*d1de - F83*rho(n,1) 832 d1t(2) = F19*rhoa*de47 - F79*rhoab*d1de - F83*rho(n,1) 833 d1f(4) = -ab*om*t 834 d2fgab(1) = -ab*( d1om*t + om*d1t(1) ) 835 d2fgab(2) = -ab*( d1om*t + om*d1t(2) ) 836#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 837 d2t(1) = -F79*( 2d0*rhob*d1de + rhoab*d2de ) - F83 838 d2t(2) = F19*de47 - F79*( rho(n,1)*d1de + rhoab*d2de ) 839 1 - F83 840 d2t(3) = -F79*( 2d0*rhoa*d1de + rhoab*d2de ) - F83 841 d3fgab(1) = -ab*( d2om*t + 2d0*d1om*d1t(1) + om*d2t(1) ) 842 d3fgab(2) = -ab*( d2om*t + d1om*( d1t(1) + d1t(2) ) 843 1 + om*d2t(2) ) 844 d3fgab(3) = -ab*( d2om*t + 2d0*d1om*d1t(2) + om*d2t(3) ) 845#endif 846c 847#ifdef THIRD_DERIV 848 d3t(1) = -F79*( 3.0d0*rhob*d2de + rhoab*d3de ) 849 d3t(2) = -F79*( 2.0d0*d1de + 2.0d0*rhob*d2de 850 1 + rhoa*d2de + rhoab*d3de ) 851 d3t(3) = -F79*( 2.0d0*d1de + rhob*d2de 852 1 + 2.0d0*rhoa*d2de + rhoab*d3de ) 853 d3t(4) = -F79*( 3.0d0*rhoa*d2de + rhoab*d3de ) 854c 855 d4fgab(1) = -ab*( d3om*t + 3.0d0*d2om*d1t(1) 856 1 + 3.0d0*d1om*d2t(1) + om*d3t(1) ) 857 d4fgab(2) = -ab*( d3om*t + 3.0d0*d2om*d1t(1) 858 1 + d1om*(d2t(1) + 2.0d0*d2t(2)) 859 2 + om*d3t(2) ) 860 d4fgab(3) = -ab*( d3om*t + 3.0d0*d2om*d1t(1) 861 1 + d1om*(d2t(1) + 2.0d0*d2t(2)) 862 2 + om*d3t(3) ) 863 d4fgab(4) = -ab*( d3om*t + 3.0d0*d2om*d1t(1) 864 1 + 3.0d0*d1om*d2t(1) + om*d3t(4) ) 865#endif 866c Daniel (3-21-13): tm_in is what I call Si (the inside term) 867 tm_in = 1.0d0 - 3.0d0*de - de11*rhob*rrho 868c 869c t = F19*rhoab*( 1d0 - 3d0*de - de11*rhob*rrho ) - rhoa2 870 t = F19*rhoab*tm_in - rhoa2 871c Daniel (10-30-12): Derivatives of the inside term 872 d1tm_in(1) = -(3.0d0+rhob*rrho)*d1de + de11*rhob*rrho2 873 d1tm_in(2) = -(3.0d0+rhob*rrho)*d1de - de11*rhoa*rrho2 874c 875c d1t(1) = F19*( rhob*( 1d0 - 3d0*de - de11*rhob*rrho ) 876c 1 + rhoab*( -( 3d0 + rhob*rrho )*d1de 877c 2 + de11*rhob*rrho2 ) ) 878c 3 - 2d0*rhoa 879c d1t(2) = F19*( rhoa*( 1d0 - 3d0*de - de11*rhob*rrho ) 880c 1 - rhoab*( ( 3d0 + rhob*rrho )*d1de 881c 2 + de11*rhoa*rrho2 ) ) 882 d1t(1) = F19*( rhob*tm_in + rhoab*d1tm_in(1) ) 883 1 - 2.0d0*rhoa 884 d1t(2) = F19*( rhoa*tm_in + rhoab*d1tm_in(2) ) 885 d1f(5) = -ab*om*t 886 d2fgbb(1) = -ab*( d1om*t + om*d1t(1) ) 887 d2fgbb(2) = -ab*( d1om*t + om*d1t(2) ) 888#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 889c Daniel (10-30-12): Derivatives of the inside term, Ri. 890 d2tm_in(1) = -(3.0d0+rhob*rrho)*d2de 891 1 + 2.0d0*d1de*rhob*rrho2 892 2 - 2.0d0*de11*rhob*rrho3 893 d2tm_in(2) = -(3.0d0+rhoa*rrho)*d2de 894 1 + (rhob - rhoa)*d1de*rrho2 895 2 + (rhoa - rhob)*de11*rrho3 896 d2tm_in(3) = -(3.0d0+rhoa*rrho)*d2de 897 1 - 2.0d0*d1de*rhoa*rrho2 898 2 + 2.0d0*de11*rhoa*rrho3 899c 900c d2t(1) = -F19*( 901c & 2d0*rhob*((3d0+rhob*rrho)*d1de-de11*rhob*rrho2) 902c & + rhoab*((3d0+rhob*rrho)*d2de-2d0*d1de*rhob*rrho2 903c & +2d0*de11*rhob*rrho3)) 904c & - 2d0 905c d2t(2) = F19*(1d0-3d0*de-de11*rhob*rrho 906c & - rho(n,1)*(3d0+rhob*rrho)*d1de 907c & - rhoab*((3d0+rhob*rrho)*d2de+d1de*(rhoa-rhob)*rrho2 908c & -de11*(rhoa-rhob)*rrho3)) 909c d2t(3) = -F19*( 910c & 2d0*rhoa*((3d0+rhob*rrho)*d1de+de11*rhoa*rrho2) 911c & + rhoab*((3d0+rhob*rrho)*d2de+2d0*d1de*rhoa*rrho2 912c & -2d0*de11*rhoa*rrho3)) 913 d2t(1) = F19*( 2.0d0*rhob*d1tm_in(1) + rhoab*d2tm_in(1) ) 914 1 - 2.0d0 915 d2t(2) = F19*( tm_in 916 1 - rho(n,1)*( 3.0d0 + rhob*rrho )*d1de 917 2 + rhoab*d2tm_in(2) ) 918 d2t(3) = F19*( 2.0d0*rhoa*d1tm_in(2) + rhoab*d2tm_in(3) ) 919c 920 d3fgbb(1) = -ab*( d2om*t + 2.0d0*d1om*d1t(1) + om*d2t(1)) 921 d3fgbb(2) = -ab*( d2om*t + d1om*( d1t(1) + d1t(2) ) 922 1 +om*d2t(2) ) 923 d3fgbb(3) = -ab*( d2om*t + 2.0d0*d1om*d1t(2) + om*d2t(3) ) 924#endif 925c 926#ifdef THIRD_DERIV 927 d3tm_in(1) = -( 3.0d0 + rhoa*rrho )*d3de 928 1 + 3.0d0*d2de*rhob*rrho2 929 2 - 6.0d0*d1de*rhob*rrho3 930 3 + 6.0d0*de11*rhob*rrho4 931 d3tm_in(2) = -( 3.0d0 + rhoa*rrho )*d3de 932 1 - d2de*rhoa*rrho2 933 2 + 2.0d0*d2de*rhob*rrho2 934 3 + 2.0d0*d1de*rhoa*rrho3 935 4 - 4.0d0*d1de*rhob*rrho3 936 5 - 2.0d0*de11*(rhoa-2.0d0*rhob)*rrho4 937 d3tm_in(3) = -( 3.0d0 + rhoa*rrho )*d3de 938 1 - 2.0d0*d2de*rhoa*rrho2 939 2 + d2de*rhob*rrho2 940 3 + 4.0d0*d1de*rhoa*rrho3 941 4 - 2.0d0*d1de*rhob*rrho3 942 5 - 2.0d0*de11*(2.0d0*rhoa-rhoa)*rrho4 943 d3tm_in(4) = -( 3.0d0 + rhoa*rrho )*d3de 944 1 - 3.0d0*d2de*rhoa*rrho2 945 2 + 6.0d0*d1de*rhoa*rrho3 946 3 - 6.0d0*de11*rhoa*rrho4 947c 948 d3t(1) = F19*( 3.0d0*rhob*d2tm_in(1) 949 1 + rhoab*d3tm_in(1) ) 950 d3t(2) = F19*( 2.0d0*d1tm_in(1) 951 1 + rhoa*d2tm_in(1) 952 2 + 2.0d0*rhob*d2tm_in(2) 953 3 + rhoab*d3tm_in(2) ) 954 d3t(3) = F19*( 2.0d0*d1tm_in(2) 955 1 + 2.0d0*rhoa*d2tm_in(2) 956 2 + rhob*d2tm_in(3) 957 3 + rhoab*d3tm_in(3) ) 958 d3t(4) = F19*( 3.0d0*rhoa*d2tm_in(3) 959 1 + rhoab*d3tm_in(4) ) 960c 961 d4fgbb(1) = -ab*( d3om*t + 3.0d0*d2om*d1t(1) 962 1 + 3.0d0*d1om*d2t(1) + om*d3t(1) ) 963 d4fgbb(2) = -ab*( d3om*t + d2om*(2.0d0*d1t(1) + d1t(2)) 964 1 + d1om*(d2t(1) + 2.0d0*d2t(2)) + om*d3t(2) ) 965 d4fgbb(3) = -ab*( d3om*t + d2om*(d1t(1) + 2.0d0*d1t(2)) 966 1 + d1om*(2.0d0*d2t(2) + d2t(3)) + om*d3t(3) ) 967 d4fgbb(4) = -ab*( d3om*t + 3.0d0*d2om*d1t(2) 968 1 + 3.0d0*d1om*d2t(3) + om*d3t(4) ) 969#endif 970c 971 f = c1*f1 + c2*f2 + gaa*d1f(3) + gab*d1f(4) + gbb*d1f(5) 972 d1f(1) = c1*d1f1(1) + c2*d1f2(1) 973 & + gaa*d2fgaa(1) + gab*d2fgab(1) + gbb*d2fgbb(1) 974 d1f(2) = c1*d1f1(2) + c2*d1f2(2) 975 & + gaa*d2fgaa(2) + gab*d2fgab(2) + gbb*d2fgbb(2) 976#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 977 d2f(1) = c1*d2f1(1) + c2*d2f2(1) 978 & + gaa*d3fgaa(1) + gab*d3fgab(1) + gbb*d3fgbb(1) 979 d2f(2) = c1*d2f1(2) + c2*d2f2(2) 980 & + gaa*d3fgaa(2) + gab*d3fgab(2) + gbb*d3fgbb(2) 981 d2f(3) = c1*d2f1(3) + c2*d2f2(3) 982 & + gaa*d3fgaa(3) + gab*d3fgab(3) + gbb*d3fgbb(3) 983#endif 984c 985#ifdef THIRD_DERIV 986 d3f(1) = c1*d3f1(1) + c2*d3f2(1) 987 1 + gaa*d4fgaa(1) + gab*d4fgab(1) + gbb*d4fgbb(1) 988 d3f(2) = c1*d3f1(2) + c2*d3f2(2) 989 1 + gaa*d4fgaa(2) + gab*d4fgab(2) + gbb*d4fgbb(2) 990 d3f(3) = c1*d3f1(3) + c2*d3f2(3) 991 1 + gaa*d4fgaa(3) + gab*d4fgab(3) + gbb*d4fgbb(3) 992 d3f(4) = c1*d3f1(4) + c2*d3f2(4) 993 1 + gaa*d4fgaa(4) + gab*d4fgab(4) + gbb*d4fgbb(4) 994#endif 995c 996 Ec = Ec + f*fac*qwght(n) 997 if (ldew) func(n) = func(n) + f*fac 998 Amat(n,1) = Amat(n,1) + d1f(1)*fac 999 Amat(n,2) = Amat(n,2) + d1f(2)*fac 1000 Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + d1f(3)*fac 1001 Cmat(n,D1_GAB) = Cmat(n,D1_GAB) + d1f(4)*fac 1002 Cmat(n,D1_GBB) = Cmat(n,D1_GBB) + d1f(5)*fac 1003#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 1004 Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) + d2f(1)*fac 1005 Amat2(n,D2_RA_RB) = Amat2(n,D2_RA_RB) + d2f(2)*fac 1006 Amat2(n,D2_RB_RB) = Amat2(n,D2_RB_RB) + d2f(3)*fac 1007 Cmat2(n,D2_RA_GAA) = Cmat2(n,D2_RA_GAA) + d2fgaa(1)*fac 1008 Cmat2(n,D2_RA_GAB) = Cmat2(n,D2_RA_GAB) + d2fgab(1)*fac 1009 Cmat2(n,D2_RA_GBB) = Cmat2(n,D2_RA_GBB) + d2fgbb(1)*fac 1010 Cmat2(n,D2_RB_GAA) = Cmat2(n,D2_RB_GAA) + d2fgaa(2)*fac 1011 Cmat2(n,D2_RB_GAB) = Cmat2(n,D2_RB_GAB) + d2fgab(2)*fac 1012 Cmat2(n,D2_RB_GBB) = Cmat2(n,D2_RB_GBB) + d2fgbb(2)*fac 1013#endif 1014c 1015#ifdef THIRD_DERIV 1016 Amat3(n,D3_RA_RA_RA) = Amat3(n,D3_RA_RA_RA) + d3f(1)*fac 1017 Amat3(n,D3_RA_RA_RB) = Amat3(n,D3_RA_RA_RB) + d3f(2)*fac 1018 Amat3(n,D3_RA_RB_RB) = Amat3(n,D3_RA_RB_RB) + d3f(3)*fac 1019 Amat3(n,D3_RB_RB_RB) = Amat3(n,D3_RB_RB_RB) + d3f(4)*fac 1020c 1021 Cmat3(n,D3_RA_RA_GAA) = Cmat3(n,D3_RA_RA_GAA) 1022 1 + d3fgaa(1)*fac 1023 Cmat3(n,D3_RA_RA_GAB) = Cmat3(n,D3_RA_RA_GAB) 1024 1 + d3fgab(1)*fac 1025 Cmat3(n,D3_RA_RA_GBB) = Cmat3(n,D3_RA_RA_GBB) 1026 1 + d3fgbb(1)*fac 1027c 1028 Cmat3(n,D3_RA_RB_GAA) = Cmat3(n,D3_RA_RB_GAA) 1029 1 + d3fgaa(2)*fac 1030 Cmat3(n,D3_RA_RB_GAB) = Cmat3(n,D3_RA_RB_GAB) 1031 1 + d3fgab(2)*fac 1032 Cmat3(n,D3_RA_RB_GBB) = Cmat3(n,D3_RA_RB_GBB) 1033 1 + d3fgbb(2)*fac 1034c 1035 Cmat3(n,D3_RB_RB_GAA) = Cmat3(n,D3_RB_RB_GAA) 1036 1 + d3fgaa(3)*fac 1037 Cmat3(n,D3_RB_RB_GAB) = Cmat3(n,D3_RB_RB_GAB) 1038 1 + d3fgab(3)*fac 1039 Cmat3(n,D3_RB_RB_GBB) = Cmat3(n,D3_RB_RB_GBB) 1040 1 + d3fgbb(3)*fac 1041#endif 1042c 1043 20 continue 1044 endif 1045 return 1046 end 1047#ifndef SECOND_DERIV 1048#define SECOND_DERIV 1049c 1050c Compile source again for the 2nd derivative case 1051c 1052#include "xc_lyp88.F" 1053#endif 1054c 1055#ifndef THIRD_DERIV 1056#define THIRD_DERIV 1057c 1058c Compile source again for the 3rd derivative case 1059c 1060#include "xc_lyp88.F" 1061#endif 1062