1#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 2#if !defined(NWAD_PRINT) 3C> \ingroup nwxc 4C> @{ 5C> 6C> \file nwxc_c_lyp.F 7C> The LYP correlation functional 8C> 9C> @} 10#endif 11#endif 12C> 13C> \ingroup nwxc_priv 14C> @{ 15C> 16C> \brief Evaluate the LYP correlation functional 17C> 18C> Evaluate the LYP correlation functional [1-4]. 19C> 20C> ### References ### 21C> 22C> [1] C. Lee, W. Yang, R.G. Parr, 23C> "Development of the Colle-Salvetti correlation-energy formula into 24C> a functional of the electron density", Phys. Rev. B <b>37</b>, 25C> 785-789 (1988), DOI: 26C> <a href="https://doi.org/10.1103/PhysRevB.37.785"> 27C> 10.1103/PhysRevB.37.785</a>. 28C> 29C> [2] R. Colle, O. Salvetti, 30C> "Approximate calculation of the correlation energy for the closed 31C> shells", Theor. Chim. Acta <b>37</b>, 329-334 (1975), DOI: 32C> <a href="https://doi.org/10.1007/BF01028401"> 33C> 10.1007/BF01028401</a>. 34C> 35C> [3] B. Miehlich, A. Savin, H. Stoll, H. Preuss, 36C> "Results obtained with the correlation energy density functionals of 37C> Becke, and Lee, Yang and Parr", Chem. Phys. Lett. <b>157</b>, 200-206 38C> (1989), DOI: <a href="https://doi.org/10.1016/0009-2614(89)87234-3"> 39C> 10.1016/0009-2614(89)87234-3</a>. 40C> 41C> [4] B.G. Johnson, P.M.W. Gill, J.A. Pople, 42C> "The performance of a family of density functional methods", 43C> J. Chem. Phys. <b>98</b>, 5612-5626 (1993), DOI: 44C> <a href="https://doi.org/10.1063/1.464906"> 45C> 10.1063/1.464906</a>. 46C> 47#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 48#if defined(NWAD_PRINT) 49 Subroutine nwxc_c_lyp_p(tol_rho, ipol, nq, wght, rho, rgamma, 50 & func) 51#else 52 Subroutine nwxc_c_lyp(tol_rho, ipol, nq, wght, rho, rgamma, 53 & func) 54#endif 55#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV) 56 Subroutine nwxc_c_lyp_d2(tol_rho, ipol, nq, wght, rho, 57 & rgamma, func) 58#else 59 Subroutine nwxc_c_lyp_d3(tol_rho, ipol, nq, wght, rho, 60 & rgamma, func) 61#endif 62c 63C$Id$ 64c 65#include "nwad.fh" 66 implicit none 67c 68#include "nwxc_param.fh" 69c 70c Input and other parameters 71c 72 double precision tol_rho !< [Input] The lower limit on the density 73 integer ipol !< [Input] The number of spin channels 74 integer nq !< [Input] The number of points 75 double precision wght !< [Input] The weight of the functional 76c 77c Charge Density 78c 79 type(nwad_dble)::rho(nq,*) !< [Input] The density 80c 81c Charge Density Gradient 82c 83 type(nwad_dble)::rgamma(nq,*) !< [Input] The norm of the density gradients 84c 85c Sampling Matrices for the XC Potential 86c 87 type(nwad_dble)::func(nq) !< [Output] The value of the functional 88c double precision Amat(nq,*) !< [Output] The derivative wrt rho 89c double precision Cmat(nq,*) !< [Output] The derivative wrt rgamma 90#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 91c 92c Sampling Matrices for the XC Kernel 93c 94c double precision Amat2(nq,*) !< [Output] The 2nd derivative wrt rho 95c double precision Cmat2(nq,*) !< [Output] The 2nd derivative wrt rgamma 96 !< and possibly rho 97#endif 98#if defined(THIRD_DERIV) 99c 100c Sampling Matrices for the XC Kernel 101c 102c double precision Amat3(nq,*) !< [Output] The 3rd derivative wrt rho 103c double precision Cmat3(nq,*) !< [Output] The 3rd derivative wrt rgamma 104 !< and possibly rho 105#endif 106 double precision F13, F43, F113, F83, F53, F19, F79, P1, 107 & A, B, C, D 108c 109 Parameter (F13 = 1.D0/3.D0, F43 = 4.D0*F13, F113 = 11.D0*F13, 110 & F83 = 8.D0*F13, F53 = 5.D0*F13, F19 = 1.D0/9.D0, 111 & F79 = 7.D0*F19) 112c 113c P1 = 2**(11/3)*(3/10)*(3*PI**2)**(2/3) 114c 115 Parameter (P1 = 0.3646239897876487D+02) 116c 117c Colle-Salvetti Empirical Parameters 118c 119 Parameter (A = 0.04918D0) 120 Parameter (B = 0.13200D0) 121 Parameter (C = 0.25330D0) 122 Parameter (D = 0.34900D0) 123c 124c Compute the partial derivatives of the correlation functional of 125c Lee, Yang and Parr. 126c 127c References: 128c 129c Colle & Salvetti, Theor. Chim. Acta 37, 329 (1975) 130c Lee, Yang & Parr, Phys. Rev. B 37, 785 (1988) 131c Miehlich, Savin, Stoll & Preuss, Chem. Phys. Lett. 157, 200 (1989) 132c Johnson, Gill & Pople, J. Chem. Phys. 98, 5612 (1993) 133c 134 integer n 135 double precision c1, c2, ab, ratcd 136 type(nwad_dble)::rrho, rhoa, rhob, d1t(2) 137 type(nwad_dble)::rhoa83, rhob83 138 type(nwad_dble)::rrho2, rhoa2, rhob2, rhoab, rho2 139 type(nwad_dble)::h1, h2, h3, om, t, tm_in 140 type(nwad_dble)::de, de11, de47 141 type(nwad_dble)::zero 142 type(nwad_dble)::xrarb 143c 144 double precision d1xrarb(2) 145 double precision d1tm_in(2) 146c 147 type(nwad_dble)::gaa, gab, gbb 148 type(nwad_dble)::f1, f2, d1f1(2), d1f2(2), f, d1f(5), 149 & d2fgaa(2), d2fgab(2), d2fgbb(2) 150c 151c Coefficients of first two terms in LYP functional and other 152c commonly occurring factors 153c 154 zero = 0.0d0 155 c1 = -4.0d0*a 156 c2 = -P1*a*b 157 ab = a*b 158 ratcd = c/d 159c 160 if (ipol.eq.1)then 161c 162c ======> SPIN-RESTRICTED <====== 163c 164 do 10 n = 1, nq 165 if (rho(n,R_T).lt.tol_rho)goto 10 166 rrho = 1.0d0/rho(n,R_T) 167 rhoa = 0.5d0*rho(n,R_T) 168 rrho2 = rrho*rrho 169 rho2 = 1.0d0/rrho2 170 rhoa2 = rhoa*rhoa 171 rhoab = rhoa2 172c rhoa53 = rhoa**F53 173 rhoa83 = rhoa**F83 174c 175 h2 = d*rrho**F13 176c d1h2 = -F13*h2*rrho 177c 178 h3 = ratcd*h2 179c d1h3 = ratcd*d1h2 180c 181 h1 = 1d0/(1d0+h2) 182c d1h1 = -h1*h1*d1h2 183c 184 om = exp(-h3)*h1*rrho**F113 185c t = d1h3+h1*d1h2+F113*rrho 186c d1om = -om*t 187c 188 de = h3+h1*h2 189c d1de = d1h3 + d1h1*h2 + h1*d1h2 190c 191 f1 = h1*rhoab*rrho 192c d1f1(1) = d1h1*rhoab*rrho - h1*rhoab*rrho2 193c d1f1(1) = d1f1(1) + h1*rhoa*rrho 194c 195 xrarb = rhoab*(rhoa83+rhoa83) 196 f2 = om*xrarb 197c d1xrarb(1) = rhoa*(F113*rhoa83+rhoa83) 198c d1f2(1) = d1om*xrarb + om*d1xrarb(1) 199c 200c gaa =(delrho(n,1,1)*delrho(n,1,1) + 201c & delrho(n,2,1)*delrho(n,2,1) + 202c & delrho(n,3,1)*delrho(n,3,1))*0.25d0 203 gaa = rgamma(n,G_TT)*0.25d0 204c 205 de11 = de - 11.0d0 206 de47 = 47.0d0 - 7.0d0*de 207c 208c Daniel (10-30-12): tm_in is what I call Qi (the inside term) 209 tm_in = 1.0d0 - 3.0d0*de - de11*rhoa*rrho 210c Daniel (10-23-12): "t" is what I call Q or S. 211 t = F19*rhoab*tm_in - rhoa2 212c Daniel (10-30-12): Derivatives of the inside term 213c d1tm_in(1) = -(3.0d0+rhoa*rrho)*d1de-de11*rhoa*rrho2 214c d1tm_in(2) = -(3.0d0+rhoa*rrho)*d1de+de11*rhoa*rrho2 215c Daniel (10-23-12): d1t(1) is the derivative with respect to rhoa, 216c and d1t(2) is the derivative with respect to rhob. 217c d1t(1) = F19*( rhoa*tm_in + rhoab*d1tm_in(1) ) 218c d1t(2) = F19*( rhoa*tm_in + rhoab*d1tm_in(2) ) 219c & - 2.0d0*rhoa 220c 221 d1f(3) = -ab*om*t 222c 223c d2fgaa(1) = -ab*( d1om*t + om*d1t(1) ) 224c d2fgaa(2) = -ab*( d1om*t + om*d1t(2) ) 225c 226c Daniel (10-23-12): "t" is what I call R. 227 t = F19*rhoab*de47-F43*rho2 228c d1t(1) = F19*rhoa*de47 - F79*rhoab*d1de - F83*rho(n,R_T) 229 d1f(4) = -ab*om*t 230c d2fgab(1) = -ab*(d1om*t+om*d1t(1)) 231c 232c d2fgbb(1) = d2fgaa(2) 233c 234 f = c1*f1 + c2*f2 + gaa*(2d0*d1f(3) + d1f(4)) 235c d1f(1) = c1*d1f1(1) + c2*d1f2(1) 236c & + gaa*(d2fgaa(1) + d2fgab(1) + d2fgbb(1)) 237c 238 func(n) = func(n) + f*wght 239c Amat(n,D1_RA) = Amat(n,D1_RA) + d1f(1)*wght 240c Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + d1f(3)*wght 241c Cmat(n,D1_GAB) = Cmat(n,D1_GAB) + d1f(4)*wght 242c 243 10 continue 244 else 245c 246c ======> SPIN-UNRESTRICTED <====== 247c 248 do 20 n = 1,nq 249 if (rho(n,R_A)+rho(n,R_B).lt.tol_rho)goto 20 250 rrho = 1.0d0/(rho(n,R_A)+rho(n,R_B)) 251c rhoa = max(zero,rho(n,R_A)) 252c rhob = max(zero,rho(n,R_B)) 253 if (rho(n,R_A).ge.0.5d0*tol_rho) then 254 rhoa = rho(n,R_A) 255 else 256 rhoa = 0.0d0 257 rrho = 1.0d0/rho(n,R_B) 258 endif 259 if (rho(n,R_B).ge.0.5d0*tol_rho) then 260 rhob = rho(n,R_B) 261 else 262 rhob = 0.0d0 263 rrho = 1.0d0/rho(n,R_A) 264 endif 265 rrho2 = rrho*rrho 266 rho2 = 1d0/rrho2 267 rhoa2 = rhoa*rhoa 268 rhob2 = rhob*rhob 269 rhoab = rhoa*rhob 270c rhoa53 = rhoa**F53 271c rhob53 = rhob**F53 272 rhoa83 = rhoa**F83 273 rhob83 = rhob**F83 274c 275 h2 = d*rrho**F13 276c d1h2 = -F13*h2*rrho 277c 278 h3 = ratcd*h2 279c d1h3 = ratcd*d1h2 280c 281 h1 = 1d0/(1d0+h2) 282c d1h1 = -h1*h1*d1h2 283c 284 om = exp(-h3)*h1*rrho**F113 285c t = d1h3+h1*d1h2+F113*rrho 286c d1om = -om*t 287c 288 de = h3+h1*h2 289c d1de = d1h3 + d1h1*h2 + h1*d1h2 290c 291c Daniel (3-21-13): f1 is J/(-4*a) in my notes. 292 f1 = h1*rhoab*rrho 293c d1f1(1) = d1h1*rhoab*rrho - h1*rhoab*rrho2 294c d1f1(2) = d1f1(1) 295c d1f1(1) = d1f1(1) + h1*rhob*rrho 296c d1f1(2) = d1f1(2) + h1*rhoa*rrho 297c 298c Daniel (10-30-12): Define xrarb here 299 xrarb = rhoab*(rhoa83+rhob83) 300c 301 f2 =om*xrarb 302c 303c d1xrarb(1) = rhob*(F113*rhoa83+rhob83) 304c d1xrarb(2) = rhoa*(F113*rhob83+rhoa83) 305c 306c d1f2(1) = d1om*xrarb 307c d1f2(2) = d1f2(1) 308c d1f2(1) = d1f2(1) + om*d1xrarb(1) 309c d1f2(2) = d1f2(2) + om*d1xrarb(2) 310c 311 gaa = rgamma(n,G_AA) 312 gbb = rgamma(n,G_BB) 313 gab = rgamma(n,G_AB) 314 if (rho(n,R_A).lt.0.5d0*tol_rho) then 315 gaa = 0.0d0 316 gab = 0.0d0 317 endif 318 if (rho(n,R_B).lt.0.5d0*tol_rho) then 319 gbb = 0.0d0 320 gab = 0.0d0 321 endif 322c 323 de11 = de - 11d0 324 de47 = 47d0 - 7d0*de 325c Daniel (10-30-12): tm_in is what I call Qi (the inside term) 326 tm_in = 1.0d0 - 3.0d0*de - de11*rhoa*rrho 327 t = F19*rhoab*tm_in-rhob2 328c Daniel (10-30-12): Derivatives of the inside term 329c d1tm_in(1) = -(3.0d0+rhoa*rrho)*d1de-de11*rhob*rrho2 330c d1tm_in(2) = -(3.0d0+rhoa*rrho)*d1de+de11*rhoa*rrho2 331c d1t(1) = F19*(rhob*tm_in + rhoab*d1tm_in(1)) 332c d1t(2) = F19*(rhoa*tm_in + rhoab*d1tm_in(2)) - 2.0d0*rhob 333c 334 d1f(3) = -ab*om*t 335c d2fgaa(1) = -ab*(d1om*t+om*d1t(1)) 336c d2fgaa(2) = -ab*(d1om*t+om*d1t(2)) 337c 338 t = F19*rhoab*de47-F43*rho2 339c d1t(1) = F19*rhob*de47 - F79*rhoab*d1de 340c & - F83*(rho(n,R_A)+rho(n,R_B)) 341c d1t(2) = F19*rhoa*de47 - F79*rhoab*d1de 342c & - F83*(rho(n,R_A)+rho(n,R_B)) 343 d1f(4) = -ab*om*t 344c d2fgab(1) = -ab*(d1om*t+om*d1t(1)) 345c d2fgab(2) = -ab*(d1om*t+om*d1t(2)) 346c Daniel (3-21-13): tm_in is what I call Si (the inside term) 347 tm_in = 1.0d0 - 3.0d0*de - de11*rhob*rrho 348 t = F19*rhoab*tm_in - rhoa2 349c Daniel (10-30-12): Derivatives of the inside term 350c d1tm_in(1) = -(3.0d0+rhob*rrho)*d1de + de11*rhob*rrho2 351c d1tm_in(2) = -(3.0d0+rhob*rrho)*d1de - de11*rhoa*rrho2 352c d1t(1) = F19*( rhob*tm_in + rhoab*d1tm_in(1) ) 353c 1 - 2.0d0*rhoa 354c d1t(2) = F19*( rhoa*tm_in + rhoab*d1tm_in(2) ) 355 d1f(5) = -ab*om*t 356c d2fgbb(1) = -ab*( d1om*t + om*d1t(1) ) 357c d2fgbb(2) = -ab*( d1om*t + om*d1t(2) ) 358c 359 f = c1*f1 + c2*f2 + gaa*d1f(3) + gab*d1f(4) + gbb*d1f(5) 360c d1f(1) = c1*d1f1(1) + c2*d1f2(1) 361c & + gaa*d2fgaa(1) + gab*d2fgab(1) + gbb*d2fgbb(1) 362c d1f(2) = c1*d1f1(2) + c2*d1f2(2) 363c & + gaa*d2fgaa(2) + gab*d2fgab(2) + gbb*d2fgbb(2) 364c 365 func(n) = func(n) + f*wght 366c Amat(n,D1_RA) = Amat(n,D1_RA) + d1f(1)*wght 367c Amat(n,D1_RB) = Amat(n,D1_RB) + d1f(2)*wght 368c Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + d1f(3)*wght 369c Cmat(n,D1_GAB) = Cmat(n,D1_GAB) + d1f(4)*wght 370c Cmat(n,D1_GBB) = Cmat(n,D1_GBB) + d1f(5)*wght 371c 372 20 continue 373 endif 374 return 375 end 376#ifndef NWAD_PRINT 377#define NWAD_PRINT 378c 379c Compile source again for Maxima 380c 381#include "nwxc_c_lyp.F" 382#endif 383#ifndef SECOND_DERIV 384#define SECOND_DERIV 385c 386c Compile source again for the 2nd derivative case 387c 388#include "nwxc_c_lyp.F" 389#endif 390#ifndef THIRD_DERIV 391#define THIRD_DERIV 392c 393c Compile source again for the 3rd derivative case 394c 395#include "nwxc_c_lyp.F" 396#endif 397#undef NWAD_PRINT 398C> @} 399