1#if defined(FUJITSU_VPP) 2!ocl scalar 3#endif 4#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 5 Subroutine xc_camb88(tol_rho, fac, lfac, nlfac, rho, delrho, 6 & Amat, Cmat, nq, ipol, Ex, qwght,ldew,func) 7#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV) 8 Subroutine xc_camb88_d2(tol_rho, fac, lfac, nlfac, rho, delrho, 9 & Amat, Amat2, Cmat, Cmat2, nq, ipol, Ex, 10 & qwght,ldew,func) 11#else 12 Subroutine xc_camb88_d3(tol_rho, fac, lfac, nlfac, rho, delrho, 13 1 Amat, Amat2, Amat3, Cmat, Cmat2, Cmat3, 14 2 nq, ipol, Ex, qwght,ldew,func) 15#endif 16c 17C$Id$ 18c 19c Coulomb attenuated Becke88 functional 20c 21 implicit none 22c 23#include "dft2drv.fh" 24#include "dft3drv.fh" 25c 26 double precision tol_rho, fac, Ex 27 integer nq, ipol 28 logical lfac, nlfac,ldew 29 double precision func(*) ! value of the functional [output] 30c 31c Charge Density 32c 33 double precision rho(nq,ipol*(ipol+1)/2) 34c 35c Charge Density Gradient 36c 37 double precision delrho(nq,3,ipol) 38c 39c Quadrature Weights 40c 41 double precision qwght(nq) 42c 43c Sampling Matrices for the XC Potential 44c 45 double precision Amat(nq,ipol), Cmat(nq,*) 46c 47#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 48c 49c Second Derivatives of the Exchange Energy Functional 50c 51 double precision Amat2(nq,NCOL_AMAT2), Cmat2(nq,NCOL_CMAT2) 52 double precision A2tmp, C2tmp, C3tmp 53#endif 54#ifdef THIRD_DERIV 55c 56c Third Derivatives of the Exchange Energy Functional 57c 58 double precision Amat3(nq,NCOL_AMAT3), Cmat3(nq,NCOL_CMAT3) 59 double precision A3tmp, C4tmp, C5tmp, C6tmp 60#endif 61c 62 double precision BETA 63 Parameter (BETA = 0.0042D0) 64c 65c References: 66c 67c Becke, Phys. Rev. A 38, 3098 (1988) 68c Johnson, Gill & Pople, J. Chem. Phys. 98, 5612 (1993) 69c 70c*************************************************************************** 71c 72 integer n 73 double precision arcsinh, darcsinh, d2arcsinh 74 double precision C, rho13, rho43, gamma, x, g, gdenom, dg, 75 & dgdenom, t, Etmp, Atmp, Ctmp 76 double precision gdenom2 77#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 78 double precision rhom23, d2g, d2gdenom 79 double precision gdenom3 80#endif 81c 82#ifdef THIRD_DERIV 83 double precision rhom53, d3g, d3gdenom 84 double precision gdenom4 85#endif 86c 87 arcsinh(x)=log(x+dsqrt(1d0+x*x)) 88 darcsinh(x)=1d0/dsqrt(1d0+x*x) 89 d2arcsinh(x) = -x/dsqrt(1d0+x*x)**3 90c 91c Uniform electron gas constant 92c 93 C = -(1.5d0)*(0.75d0/acos(-1d0))**(1d0/3d0) 94c 95 if (ipol.eq.1) then 96c 97c ======> SPIN-RESTRICTED <====== 98c 99 do 10 n = 1, nq 100 if (rho(n,1).lt.tol_rho) goto 10 101c 102c Spin alpha: 103c 104 rho13 = (0.5d0*rho(n,1))**(1.d0/3.d0) 105 rho43 = rho13**4 106 gamma = delrho(n,1,1)*delrho(n,1,1) + 107 & delrho(n,2,1)*delrho(n,2,1) + 108 & delrho(n,3,1)*delrho(n,3,1) 109 if (dsqrt(gamma).gt.tol_rho)then 110 gamma = 0.25d0 * gamma 111 x = dsqrt(gamma) / rho43 112 else 113 x = 0d0 114 endif 115c 116 gdenom = 1d0 + 6d0*BETA*x*arcsinh(x) 117 gdenom2 = gdenom*gdenom 118 g = -BETA*x*x / gdenom 119 dgdenom = 6d0*BETA*(arcsinh(x) + x*darcsinh(x)) 120c dg = BETA*x*(x*dgdenom - 2d0*gdenom) / gdenom**2 121 dg = BETA*x*(x*dgdenom - 2d0*gdenom) / gdenom2 122c 123 Etmp = 0.d0 124 Atmp = 0.d0 125 Ctmp = 0.d0 126 if (lfac) then 127 Etmp = 2d0*rho43*C 128 Atmp = (4d0/3d0)*rho13*C 129 endif 130c 131 if (nlfac) then 132 Etmp = Etmp + 2d0*rho43*g 133 Atmp = Atmp + (4d0/3d0)*rho13*(g-x*dg) 134 endif 135c 136 if (x.gt.tol_rho) then 137 Ctmp = 0.5d0 * dg / sqrt(gamma) 138 endif 139c 140#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 141 A2tmp = 0d0 142 C2tmp = 0d0 143 C3tmp = 0d0 144 if(lfac) g = g + C ! Add local contribution back to g 145 rhom23 = rho13 / (0.5d0*rho(n,1)) 146 d2gdenom = 6d0*BETA*darcsinh(x)*(2d0 - x*x/(x*x+1d0)) 147 gdenom3 = gdenom2*gdenom 148 d2g = -2d0*BETA/gdenom + 4d0*BETA*x*dgdenom/gdenom2 149 & + BETA*x*x*d2gdenom/gdenom2 150 & - 2d0*BETA*x*x*(dgdenom)**2/gdenom3 151c 152c rr 153 A2tmp = (4d0/9d0)*rhom23*(g-x*dg+4d0*x*x*d2g) 154c rg 155 C2tmp = - (4d0/3d0)*(rhom23**2/rho(n,1))*d2g 156c gg 157 if (x.gt.tol_rho) then 158 C3tmp = - 0.25d0*gamma**(-1.5d0)*(dg-x*d2g) 159 endif 160#endif 161#ifdef THIRD_DERIV 162 A3tmp = 0.0d0 163 C4tmp = 0.0d0 164 C5tmp = 0.0d0 165 C6tmp = 0.0d0 166c 167 rhom53 = rhom23 / (0.5d0*rho(n,1)) 168c 169 d3gdenom = 6.0d0*BETA* 170 1 d2arcsinh(x)*( 3.0d0 171 2 - (2.0d0*x*x - 1.0d0)/(1.0d0 + x*x)) 172c 173 gdenom4 = gdenom3*gdenom 174c 175 d3g = 6.0d0*BETA*dgdenom/gdenom2 176 1 - 12.0d0*BETA*x*dgdenom*dgdenom/gdenom3 177 2 + 6.0d0*BETA*x*d2gdenom/gdenom2 178 3 + 6.0d0*BETA*x*x*dgdenom*dgdenom*dgdenom/gdenom4 179 4 - 6.0d0*BETA*x*x*dgdenom*d2gdenom/gdenom3 180 5 + BETA*x*x*d3gdenom/gdenom2 181c 182c rrr 183 A3tmp = (8.0d0/27.0d0)*rhom53*(-g + x*dg 184 1 - 18.0d0*x*x*d2g 185 2 - 8.0d0*x*x*x*d3g) 186c 187c rrg 188 C4tmp = (2.0d0/9.0d0)*(rhom23/gamma)*( 7.0d0*x*x*d2g 189 1 + 4.0d0*x*x*x*d3g) 190c 191c rgg 192 C5tmp = -(8.0d0/3.0d0)*(rhom23/rho(n,1)**3)/dsqrt(gamma) 193 1 *d3g 194c 195c ggg 196 if (x.gt.tol_rho) then 197 C6tmp = (1.0d0/8.0d0)*gamma**(-2.5d0)*( 3.0d0*dg 198 1 - 3.0d0*x*d2g 199 2 + x*x*d3g) 200 endif 201#endif 202c 203#ifdef THIRD_DERIV 204 call xc_att_xc_d3(rho(n,1),ipol,Etmp,Atmp,Ctmp,A2tmp, 205 & C2tmp,C3tmp,A3tmp,C4tmp,C5tmp,C6tmp) 206c 207 Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) + A2tmp*fac 208 Cmat2(n,D2_RA_GAA) = Cmat2(n,D2_RA_GAA) + C2tmp*fac 209 Cmat2(n,D2_GAA_GAA) = Cmat2(n,D2_GAA_GAA) + C3tmp*fac 210c 211 Amat3(n,D3_RA_RA_RA) = Amat3(n,D3_RA_RA_RA) + A3tmp*fac 212 Cmat3(n,D3_RA_RA_GAA) = Cmat3(n,D3_RA_RA_GAA) + C4tmp*fac 213 Cmat3(n,D3_RA_GAA_GAA) = Cmat3(n,D3_RA_GAA_GAA) + C5tmp* 214 * fac 215 Cmat3(n,D3_GAA_GAA_GAA) = Cmat3(n,D3_GAA_GAA_GAA) + C6tmp* 216 * fac 217#elif defined(SECOND_DERIV) 218 call xc_att_xc_d2(rho(n,1),ipol,Etmp,Atmp,Ctmp,A2tmp, 219 & C2tmp,C3tmp) 220c 221 Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) + A2tmp*fac 222 Cmat2(n,D2_RA_GAA) = Cmat2(n,D2_RA_GAA) + C2tmp*fac 223 Cmat2(n,D2_GAA_GAA) = Cmat2(n,D2_GAA_GAA) + C3tmp*fac 224#else 225 call xc_att_xc(rho(n,1),ipol,Etmp,Atmp,Ctmp) 226#endif 227 Ex = Ex + qwght(n)*Etmp*fac 228 if(ldew) func(n) = func(n) + Etmp*fac 229 Amat(n,1) = Amat(n,1) + Atmp*fac 230 Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + Ctmp*fac 231 10 continue 232c 233 else 234c 235c ======> SPIN-UNRESTRICTED <====== 236c 237 do 20 n = 1, nq 238 if (rho(n,1).lt.tol_rho) goto 20 239 if (rho(n,2).lt.tol_rho) goto 25 240c 241c Spin alpha: 242c 243 rho13 = rho(n,2)**(1.d0/3.d0) 244 rho43 = rho13*rho(n,2) 245 gamma = delrho(n,1,1)*delrho(n,1,1) + 246 & delrho(n,2,1)*delrho(n,2,1) + 247 & delrho(n,3,1)*delrho(n,3,1) 248 if (dsqrt(gamma).gt.tol_rho)then 249 x = dsqrt(gamma) / rho43 250 else 251 x = 0d0 252 endif 253c 254 gdenom = 1d0 + 6d0*BETA*x*arcsinh(x) 255 g = -BETA*x*x / gdenom 256 dgdenom = 6d0*BETA*(arcsinh(x) + x*darcsinh(x)) 257 gdenom2 = gdenom*gdenom 258 dg = BETA*x*(x*dgdenom - 2d0*gdenom) / gdenom2 259c 260 Etmp = 0.d0 261 Atmp = 0.d0 262 Ctmp = 0.d0 263 if (lfac) then 264 Etmp = rho43*C 265 Atmp = (4d0/3d0)*rho13*C 266 endif 267c 268 if (nlfac) then 269 Etmp = Etmp + rho43*g 270 Atmp = Atmp + (4d0/3d0)*rho13*(g-x*dg) 271 endif 272c 273 if (x.gt.tol_rho) then 274 Ctmp = 0.5d0*dg / sqrt(gamma) 275 endif 276c 277#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 278 if (lfac) g = g + C ! Add local contribution back to g 279 rhom23 = rho13 / rho(n,2) 280 d2gdenom = 6d0*BETA*darcsinh(x)*(2d0 - x*x/(x*x+1d0)) 281 gdenom3 = gdenom2*gdenom 282 d2g = -2d0*BETA/gdenom + 4d0*BETA*x*dgdenom/gdenom2 283 & + BETA*x*x*d2gdenom/gdenom2 284 & - 2d0*BETA*x*x*(dgdenom)**2/gdenom3 285c d2g = -2d0*BETA/gdenom + 4d0*BETA*x*dgdenom/gdenom**2 286c & + BETA*x*x*d2gdenom/gdenom**2 287c & - 2d0*BETA*x*x*(dgdenom)**2/gdenom**3 288c 289 A2tmp = (4d0/9d0)*rhom23*(g-x*dg+4d0*x*x*d2g) 290c C2tmp = (2d0/3d0)*(rhom23**2/rho(n,2))*d2g 291 C2tmp = -(2d0/3d0)*(rhom23**2/rho(n,2))*d2g 292 if (x.gt.tol_rho) then 293 C3tmp = -0.25d0*gamma**(-1.5d0)*(dg-x*d2g) 294 endif 295#endif 296#ifdef THIRD_DERIV 297 rhom53 = rhom23 / rho(n,2) 298c 299 d3gdenom = 6.0d0*BETA* 300 1 d2arcsinh(x)*( 3.0d0 301 2 - (2.0d0*x*x - 1.0d0)/(1.0d0 + x*x)) 302c 303 gdenom4 = gdenom3*gdenom 304c 305 d3g = 6.0d0*BETA*dgdenom/gdenom2 306 1 - 12.0d0*BETA*x*dgdenom*dgdenom/gdenom3 307 2 + 6.0d0*BETA*x*d2gdenom/gdenom2 308 3 + 6.0d0*BETA*x*x*dgdenom*dgdenom*dgdenom/gdenom4 309 4 - 6.0d0*BETA*x*x*dgdenom*d2gdenom/gdenom3 310 5 + BETA*x*x*d3gdenom/gdenom2 311c 312 A3tmp = (8.0d0/27.0d0)*rhom53*(-g + x*dg 313 1 - 18.0d0*x*x*d2g 314 2 - 8.0d0*x*x*x*d3g) 315c 316 C4tmp = (2.0d0/9.0d0)*(rhom23/gamma)*( 7.0d0*x*x*d2g 317 1 + 4.0d0*x*x*x*d3g) 318c 319 C5tmp = -(1.0d0/3.0d0)*(rhom23/rho(n,2)**3)/dsqrt(gamma) 320 1 *d3g 321c 322 if (x.gt.tol_rho) then 323 C6tmp = (1.0d0/8.0d0)*gamma**(-2.5d0)*( 3.0d0*dg 324 1 - 3.0d0*x*d2g 325 2 + x*x*d3g) 326 endif 327#endif 328c 329#ifdef THIRD_DERIV 330 call xc_att_xc_d3(rho(n,2),ipol,Etmp,Atmp,Ctmp,A2tmp, 331 & C2tmp,C3tmp,A3tmp,C4tmp,C5tmp,C6tmp) 332c 333 Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) + A2tmp*fac 334 Cmat2(n,D2_RA_GAA) = Cmat2(n,D2_RA_GAA) + C2tmp*fac 335 Cmat2(n,D2_GAA_GAA) = Cmat2(n,D2_GAA_GAA) + C3tmp*fac 336c 337 Amat3(n,D3_RA_RA_RA) = Amat3(n,D3_RA_RA_RA) + A3tmp*fac 338 Cmat3(n,D3_RA_RA_GAA) = Cmat3(n,D3_RA_RA_GAA) + C4tmp*fac 339 Cmat3(n,D3_RA_GAA_GAA) = Cmat3(n,D3_RA_GAA_GAA) + C5tmp* 340 * fac 341 Cmat3(n,D3_GAA_GAA_GAA) = Cmat3(n,D3_GAA_GAA_GAA) + C6tmp* 342 * fac 343#elif defined(SECOND_DERIV) 344 call xc_att_xc_d2(rho(n,2),ipol,Etmp,Atmp,Ctmp,A2tmp,C2tmp, 345 & C3tmp) 346 Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) + A2tmp*fac 347 Cmat2(n,D2_RA_GAA) = Cmat2(n,D2_RA_GAA) + C2tmp*fac 348 Cmat2(n,D2_GAA_GAA) = Cmat2(n,D2_GAA_GAA) + C3tmp*fac 349#else 350 call xc_att_xc(rho(n,2),ipol,Etmp,Atmp,Ctmp) 351#endif 352 Ex = Ex + qwght(n)*Etmp*fac 353 if(ldew) func(n) = func(n) + Etmp*fac 354 Amat(n,1) = Amat(n,1) + Atmp*fac 355 Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + Ctmp*fac 356c 357 25 continue 358c 359c Spin beta: 360c 361 if (rho(n,3).lt.tol_rho) goto 20 362c 363 rho13 = rho(n,3)**(1.d0/3.d0) 364 rho43 = rho13*rho(n,3) 365 gamma = delrho(n,1,2)*delrho(n,1,2) + 366 & delrho(n,2,2)*delrho(n,2,2) + 367 & delrho(n,3,2)*delrho(n,3,2) 368 if (dsqrt(gamma).gt.tol_rho)then 369 x = dsqrt(gamma) / rho43 370 else 371 x = 0d0 372 endif 373c 374 gdenom = 1d0 + 6d0*BETA*x*arcsinh(x) 375 g = -BETA*x*x / gdenom 376 dgdenom = 6d0*BETA*(arcsinh(x) + x*darcsinh(x)) 377 gdenom2 = gdenom*gdenom 378 dg = BETA*x*(x*dgdenom - 2d0*gdenom) / gdenom2 379c 380 Etmp = 0.d0 381 Atmp = 0.d0 382 Ctmp = 0.d0 383 if (lfac) then 384 Etmp = rho43*C 385 Atmp = (4d0/3d0)*rho13*C 386 endif 387c 388 if (nlfac) then 389 Etmp = Etmp + rho43*g 390 Atmp = Atmp + (4d0/3d0)*rho13*(g-x*dg) 391 endif 392c 393 if (x.gt.tol_rho) then 394 Ctmp = 0.5d0*dg / sqrt(gamma) 395 endif 396c 397#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 398 A2tmp = 0d0 399 C2tmp = 0d0 400 C3tmp = 0d0 401 if(lfac) g = g + C ! Add local contribution back to g 402 rhom23 = rho13 / rho(n,3) 403 d2gdenom = 6d0*BETA*darcsinh(x)*(2d0 - x*x/(x*x+1d0)) 404 gdenom3 = gdenom2*gdenom 405 d2g = -2d0*BETA/gdenom + 4d0*BETA*x*dgdenom/gdenom2 406 & + BETA*x*x*d2gdenom/gdenom2 407 & - 2d0*BETA*x*x*(dgdenom)**2/gdenom3 408 A2tmp = (4d0/9d0)*rhom23*(g-x*dg+4d0*x*x*d2g) 409 C2tmp = -(2d0/3d0)*(rhom23**2/rho(n,3))*d2g 410 if (x.gt.tol_rho) then 411 C3tmp = - 0.25d0*gamma**(-1.5d0)*(dg-x*d2g) 412 endif 413#endif 414#ifdef THIRD_DERIV 415 rhom53 = rhom23 / rho(n,3) 416c 417 d3gdenom = 6.0d0*BETA* 418 1 d2arcsinh(x)*( 3.0d0 419 2 - (2.0d0*x*x - 1.0d0)/(1.0d0 + x*x)) 420c 421 gdenom4 = gdenom3*gdenom 422c 423 d3g = 6.0d0*BETA*dgdenom/gdenom2 424 1 - 12.0d0*BETA*x*dgdenom*dgdenom/gdenom3 425 2 + 6.0d0*BETA*x*d2gdenom/gdenom2 426 3 + 6.0d0*BETA*x*x*dgdenom*dgdenom*dgdenom/gdenom4 427 4 - 6.0d0*BETA*x*x*dgdenom*d2gdenom/gdenom3 428 5 + BETA*x*x*d3gdenom/gdenom2 429c 430 A3tmp = (8.0d0/27.0d0)*rhom53*(-g + x*dg 431 1 - 18.0d0*x*x*d2g 432 2 - 8.0d0*x*x*x*d3g) 433c 434 C4tmp = (2.0d0/9.0d0)*(rhom23/gamma)*( 7.0d0*x*x*d2g 435 1 + 4.0d0*x*x*x*d3g) 436c 437 C5tmp = -(1.0d0/3.0d0)*(rhom23/rho(n,3)**3)/dsqrt(gamma) 438 1 *d3g 439c 440 if (x.gt.tol_rho) then 441 C6tmp = (1.0d0/8.0d0)*gamma**(-2.5d0)*( 3.0d0*dg 442 1 - 3.0d0*x*d2g 443 2 + x*x*d3g) 444 endif 445#endif 446c 447#ifdef THIRD_DERIV 448 call xc_att_xc_d3(rho(n,3),ipol,Etmp,Atmp,Ctmp,A2tmp, 449 & C2tmp,C3tmp,A3tmp,C4tmp,C5tmp,C6tmp) 450c 451 Amat2(n,D2_RB_RB) = Amat2(n,D2_RB_RB) + A2tmp*fac 452 Cmat2(n,D2_RB_GBB) = Cmat2(n,D2_RB_GBB) + C2tmp*fac 453 Cmat2(n,D2_GBB_GBB) = Cmat2(n,D2_GBB_GBB) + C3tmp*fac 454c 455 Amat3(n,D3_RB_RB_RB) = Amat3(n,D3_RB_RB_RB) + A3tmp*fac 456 Cmat3(n,D3_RB_RB_GBB) = Cmat3(n,D3_RB_RB_GBB) + C4tmp*fac 457 Cmat3(n,D3_RB_GBB_GBB) = Cmat3(n,D3_RB_GBB_GBB) + C5tmp* 458 * fac 459 Cmat3(n,D3_GBB_GBB_GBB) = Cmat3(n,D3_GBB_GBB_GBB) + C6tmp* 460 * fac 461#elif defined(SECOND_DERIV) 462 call xc_att_xc_d2(rho(n,3),ipol,Etmp,Atmp,Ctmp,A2tmp,C2tmp, 463 & C3tmp) 464 Amat2(n,D2_RB_RB) = Amat2(n,D2_RB_RB) + A2tmp*fac 465 Cmat2(n,D2_RB_GBB) = Cmat2(n,D2_RB_GBB) + C2tmp*fac 466 Cmat2(n,D2_GBB_GBB) = Cmat2(n,D2_GBB_GBB) + C3tmp*fac 467#else 468 call xc_att_xc(rho(n,3),ipol,Etmp,Atmp,Ctmp) 469#endif 470 Ex = Ex + qwght(n)*Etmp*fac 471 if(ldew) func(n) = func(n) + Etmp*fac 472 Amat(n,2) = Amat(n,2) + Atmp*fac 473 Cmat(n,D1_GBB) = Cmat(n,D1_GBB) + Ctmp*fac 474 20 continue 475c 476 endif 477c 478 return 479 end 480#ifndef SECOND_DERIV 481#define SECOND_DERIV 482c 483c Compile source again for the 2nd derivative case 484c 485#include "xc_camb88.F" 486#endif 487c 488#ifndef THIRD_DERIV 489#define THIRD_DERIV 490c 491c Compile source again for the 3rd derivative case 492c 493#include "xc_camb88.F" 494#endif 495