1c 2c ----------------------------------------------------------------------- 3c Uniform electron gas exchange functional for the erfc(r)/r interaction 4c as implemented in the following paper: 5c "A well-tempered density functional theory of electrons in molecules" 6c Ester Livshits & Roi Baer, Phys. Chem. Chem. Phys., 9, 2932 (2007) 7c The other relevant publication is: 8c R. Baer, D. Neuhauser, Phys. Rev. Lett., 94, 043002 (2005) 9c ----------------------------------------------------------------------- 10c 11#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 12 subroutine xc_bnl(tol_rho, fac, lfac, nlfac, rho, Amat, nq, 13 & ipol, Ex, qwght, ldew, func) 14#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV) 15c For locations of 2nd derivatives of functionals in array 16#include "dft2drv.fh" 17 subroutine xc_bnl_d2(tol_rho, fac, lfac, nlfac, rho, Amat, 18 & Amat2, nq, ipol, Ex, qwght, ldew, func) 19#else 20#include "dft3drv.fh" 21 subroutine xc_bnl_d3(tol_rho, fac, lfac, nlfac, rho, Amat, 22 1 Amat2, Amat3, nq, ipol, Ex, qwght, ldew, func) 23#endif 24c 25 implicit none 26c 27#include "errquit.fh" 28#include "stdio.fh" 29c 30ccase...start 31#include "case.fh" 32ccase...end 33c 34 integer nq, ipol, n 35 double precision fac, Ex, total 36 logical ldew, lfac, nlfac 37 double precision func(*) ! value of the functional [output] 38 double precision tol_rho 39 double precision rho(nq,(ipol*(ipol+1))/2) ! charge density 40 double precision qwght(nq) ! quadrature weights 41 double precision Amat(nq,ipol) ! partial first derivatives 42 double precision F(nq),RA(nq),RB(nq) 43 double precision rhoA, rhoB, rhoTotal, rhoA1, rhoB1 44 double precision gamma 45c double precision fA, fB, fpA, fpB, fppA, fppB 46 double precision fA, fB, fpA, fpB 47 double precision EpsX 48 double precision EpsXprime 49c double precision EpsTwoXprime 50c 51#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 52c double precision Amat2(nq,*) ! partial second derivatives 53 double precision Amat2(nq,NCOL_AMAT2) ! partial second derivatives 54 double precision fppA, fppB 55 double precision EpsTwoXprime 56#endif 57#ifdef THIRD_DERIV 58 double precision Amat3(nq,NCOL_AMAT3) 59 double precision fpppA, fpppB 60 double precision EpsThreeXprime 61#endif 62c 63c ----------------------------------------------------------------------- 64c Preliminaries 65c ----------------------------------------------------------------------- 66c 67 gamma = cam_omega 68c 69 do n = 1,nq 70 if (ipol.eq.1) then ! spin-restricted 71 rA(n) = rho(n,1) 72 rB(n) = 0.d0 73 else ! spin-unrestricted 74 rA(n) = rho(n,2) 75 rB(n) = rho(n,3) 76 end if 77 end do 78c 79c ----------------------------------------------------------------------- 80c Calculate the first and second derivatives 81c ----------------------------------------------------------------------- 82c 83 total = 0.d0 84 do n = 1,nq 85 rhoA = rA(n) 86 rhoB = rB(n) 87 rhoTotal = rhoA + rhoB ! total density at point 88 if (rhoTotal.gt.tol_rho) then 89 90 if (ipol.eq.1) then ! spin-restricted 91 rhoA1 = rhoA 92 rhoB1 = rhoB 93 else ! spin-unrestricted 94 rhoA1 = rhoA*2.0d0 95 rhoB1 = rhoB*2.0d0 96 end if 97 98 fA = EpsX(rhoA1,gamma) 99 fB = EpsX(rhoB1,gamma) 100 fpA = EpsXprime(rhoA1,gamma) 101 fpB = EpsXprime(rhoB1,gamma) 102 103 f(n) = fA * rhoA + fB * rhoB 104 Amat(n,1) = Amat(n,1) + (fpA*rhoA1+fA)*fac 105 if (ipol.gt.1) then 106 Amat(n,2) = Amat(n,2) + (fpB*rhoB1+fB)*fac 107 end if 108 109#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 110 fppA = EpsTwoXprime(rhoA1,gamma) 111 fppB = EpsTwoXprime(rhoB1,gamma) 112c 113 if (ipol.eq.1) then 114 Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) + 115 & (fppA*rhoA+2.0d0*fpA)*fac*2.0d0 116 else 117 Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) + 118 & (fppA*rhoA+fpA)*fac*4.0d0 119c Guard against case of no beta electrons, e.g. H atom 120 if (rho(n,3).gt.tol_rho) then 121c Amat2(n,D2_RB_RB) = Amat2(n,D2_RB_RB) + 122c & ((fppB*rhoB+fpB)*4)*fac 123 Amat2(n,D2_RB_RB) = Amat2(n,D2_RB_RB) + 124 & (fppB*rhoB+fpB)*fac*4.0d0 125 end if 126 end if 127#endif 128#ifdef THIRD_DERIV 129 fpppA = EpsThreeXprime(rhoA1,gamma) 130 fpppB = EpsThreeXprime(rhoB1,gamma) 131c 132 if (ipol.eq.1) then 133 Amat3(n,D3_RA_RA_RA) = Amat3(n,D3_RA_RA_RA) 134 1 + ( fpppA*rhoA + 3.0d0*fppA )*fac*4.0d0 135 else 136c Amat3(n,D3_RA_RA_RA) = Amat3(n,D3_RA_RA_RA) 137c 1 + ( fpppA*rhoA + 3.0d0*fppA )*fac*4.0d0 138 Amat3(n,D3_RA_RA_RA) = Amat3(n,D3_RA_RA_RA) 139 1 + ( fpppA*rhoA + 1.5d0*fppA )*fac*8.0d0 140c Guard against case of no beta electrons, e.g. H atom 141 if (rho(n,3).gt.tol_rho) then 142c Amat3(n,D3_RB_RB_RB) = Amat3(n,D3_RB_RB_RB) 143c & + ( fpppB*rhoB + 3.0d0*fppB )*fac*4.0d0 144 Amat3(n,D3_RB_RB_RB) = Amat3(n,D3_RB_RB_RB) 145 & + ( fpppB*rhoB + 1.5d0*fppB )*fac*8.0d0 146 end if 147 end if 148#endif 149 if (ldew) func(n) = func(n) + f(n)*fac 150 total = total + f(n)*qwght(n) 151 end if 152 end do 153 154 Ex = Ex + total*fac 155 156 return 157 end 158c 159c 160#ifndef SECOND_DERIV 161#define SECOND_DERIV 162c 163#include "xc_bnl.F" 164c 165c --------------------------------------------------------------------------------------- 166c Utility functions 167c --------------------------------------------------------------------------------------- 168c 169c --------------------------------------------------------------------------------------- 170c Return the value of pi 171c --------------------------------------------------------------------------------------- 172c 173 double precision function ValueOfPi() 174c 175 implicit none 176#include "xc_params.fh" 177c 178 ValueOfPi = pi 179 180 return 181 end 182c 183c --------------------------------------------------------------------------------------- 184c Evaluates the actual function 185c --------------------------------------------------------------------------------------- 186c 187 double precision function HqBNL(q) 188 189 implicit none 190#include "xc_params.fh" 191 192 double precision q,TwoSqrtPi,OneOverQ,q2,DERF 193 194 if (q .lt. 1D-15) then 195 HqBNL=1.d0 196 return 197 end if 198 199 OneOverQ = 1.0d0/q 200 TwoSqrtPi = 2.0d0*dsqrt(pi) 201 q2 = q**2.0d0 202 203 if (q .lt. 0.1d0) then 204 HqBNL=1.0d0-q*2.0d0/3.0d0*(TwoSqrtPi-q+q*(q2-2.0d0)) 205 return 206 end if 207 208 HqBNL=1.0d0-q*2.0d0/3.0d0*(TwoSqrtPi*DERF(OneOverQ)-q+ 209 $ q*(q2-2.0d0)*(1.0d0-exp(-OneOverQ*OneOverQ))) 210 211 return 212 end 213c 214c --------------------------------------------------------------------------------------- 215c Calculate the local Fermi vector for the provided density 216c --------------------------------------------------------------------------------------- 217c 218 double precision function FermiK(den) 219 220 implicit none 221#include "xc_params.fh" 222 223 double precision F13, den 224 225 F13 = 1.0D0 / 3.0D0 226 FermiK = (3.d0*pi*pi*den)**F13 227 228 return 229 end 230c 231c --------------------------------------------------------------------------------------- 232c Calculate the function EpsX at the given density value and gamma 233c --------------------------------------------------------------------------------------- 234c 235 double precision function EpsX(Rho,gamma) 236 237 implicit none 238#include "xc_params.fh" 239 240 double precision kF,RHO,gamma,Cs 241 double precision HqBNL 242 double precision FermiK 243 244 if (RHO.le.0D0) then 245 EpsX = 0.0D0 246 return 247 end if 248 249 kF = FermiK(Rho) 250 Cs = -3.0D0/(4.0d0*pi) 251 EpsX = Cs * kF * HqBNL(gamma/kF) 252 253 return 254 end 255c 256c --------------------------------------------------------------------------------------- 257c Calculate the first derivative of the function 258c --------------------------------------------------------------------------------------- 259c 260 double precision function HqBNLPrime(q) 261 262 implicit none 263#include "xc_params.fh" 264 265 double precision q,OneOverQ,q2,q3,DERF 266 267 q2 = q**2.0d0 268 q3 = q**3.0d0 269 if (q .lt. 0.1d0) then 270 HqBNLPrime = -4.0d0/3.0d0*(dsqrt(Pi)+2.0d0*q3-3.0d0*q) 271 return 272 end if 273 274 OneOverQ = 1.0d0/q 275 276 HqBNLPrime = 4.0d0/3.0d0*(q*(exp(-OneOverQ*OneOverQ)*(2.0d0*q2 277 $ -1.0d0)+(3.0d0-2.0d0*q2))-dsqrt(Pi)*DERF(OneOverQ)) 278 279 return 280 end 281c 282c --------------------------------------------------------------------------------------- 283c Calculate the first derivative of the local Fermi vector (it depends on the density) 284c --------------------------------------------------------------------------------------- 285c 286 double precision function FermiKPrime(den) 287 288 implicit none 289#include "xc_params.fh" 290 291 double precision F23, den 292 293 F23 = 2.0D0 / 3.0D0 294 FermiKPrime = (Pi/(3.0d0*den))**F23 295 296 return 297 end 298c 299c --------------------------------------------------------------------------------------- 300c Calculate the first derivative of q (q=gamma/kf) (it implicitly depends on the density) 301c --------------------------------------------------------------------------------------- 302c 303 double precision function QPrime(gamma,kF) 304 305 implicit none 306 307 double precision kF, FermiK2, gamma 308 309 FermiK2 = kF**2.0d0 310 QPrime = -gamma/FermiK2 311 312 return 313 end 314c 315c --------------------------------------------------------------------------------------- 316c Calculate the first derivative of EpsX 317c --------------------------------------------------------------------------------------- 318c 319 double precision function EpsXprime(Rho,gamma) 320 321 implicit none 322#include "xc_params.fh" 323 324 double precision Rho,gamma 325 double precision Cs,kF,CsPrime 326 327 double precision HqBNL 328 double precision HqBNLPrime 329 double precision QPrime 330 double precision FermiK 331 double precision FermiKPrime 332 333 kF = FermiK(Rho) 334 CsPrime = -3.0D0/(4.0d0*Pi) 335 Cs = CsPrime*kF 336 337 if (Rho.le.0d0) then 338 EpsXprime = 0.0d0 339 return 340 end if 341 342 EpsXprime = FermiKPrime(Rho)*(CsPrime*HqBNL(gamma/kF)+ 343 $ QPrime(gamma,kF)*HqBNLPrime(gamma/kF)*Cs) 344 345 return 346 end 347c 348c --------------------------------------------------------------------------------------- 349c Calculate the second derivative of the main function that consititutes the functional 350c --------------------------------------------------------------------------------------- 351c 352 double precision function HqBNLTwoPrime(q) 353 354 implicit none 355#include "xc_params.fh" 356 357 double precision q,OneOverQ,q2 358 359 q2 = q**2.0d0 360 if (q .lt. 0.1d0) then 361 HqBNLTwoPrime = 4.0d0-8.0d0*q2 362 return 363 end if 364 365 OneOverQ = 1.0d0/q 366 367 HqBNLTwoPrime = exp(-OneOverQ*OneOverQ)*(4.0d0+8.0d0*q2) 368 $ -8.0d0*q2+4.0d0 369 370 return 371 end 372c 373c --------------------------------------------------------------------------------------- 374c Calculate the second derivative of the local Fermi vector 375c --------------------------------------------------------------------------------------- 376c 377 double precision function FermiKTwoPrime(den) 378 379 implicit none 380#include "xc_params.fh" 381 382 double precision F13, den 383 384 F13 = 1.0D0/3.0D0 385 FermiKTwoPrime = -(8.0d0*Pi**2.0d0/(243.0d0*den**5.0d0))**F13 386 387 return 388 end 389c 390c --------------------------------------------------------------------------------------- 391c Calculate the second derivative of q 392c --------------------------------------------------------------------------------------- 393c 394 double precision function QTwoPrime(gamma,kF) 395 396 implicit none 397 398 double precision gamma, kF, FermiK3 399 400 FermiK3 = kF**3.0d0 401 QTwoPrime = (2.0d0*gamma)/FermiK3 402 403 return 404 end 405c 406c --------------------------------------------------------------------------------------- 407c Calculate the second derivative of EpsX 408c --------------------------------------------------------------------------------------- 409c 410 double precision function EpsTwoXprime(Rho,gamma) 411 412 implicit none 413#include "xc_params.fh" 414 415 double precision Rho,gamma 416 double precision kF,kFPrim,kFprim2,kF2prim 417 double precision q,qprim,qprim2,q2prim 418 double precision g,gprim,g2prim 419 double precision Cs,CsPrim 420 421 double precision FermiK 422 double precision FermiKPrime 423 double precision FermiKTwoPrime 424 double precision QPrime 425 double precision QTwoPrime 426 double precision HqBNL 427 double precision HqBNLPrime 428 double precision HqBNLTwoPrime 429 430 if (Rho.le.0d0) then 431 EpsTwoXprime = 0.0d0 432 return 433 end if 434 435 kF = FermiK(Rho) 436 kFPrim = FermiKPrime(Rho) 437 kFPrim2=kFPrim**2.0d0 438 kF2prim = FermiKTwoPrime(Rho) 439 CsPrim = -3.0d0/(4.0d0*Pi) 440 Cs = CsPrim * kF 441 q = gamma / kF 442 qprim = QPrime(gamma,kF) 443 Qprim2=qprim**2.0d0 444 q2prim = QTwoPrime(gamma,kF) 445 g = HqBNL(q) 446 gprim = HqBNLPrime(q) 447 g2prim = HqBNLTwoPrime(q) 448 449 EpsTwoXprime = 450 $ kFPrim2*(2.0d0*CsPrim*gprim*qprim 451 $ +Cs*(QPrim2*g2prim+gprim*Q2Prim)) 452 $ +kF2Prim*(g*CsPrim+Cs*gprim*qprim) 453 454 return 455 end 456#endif 457#ifndef THIRD_DERIV 458#define THIRD_DERIV 459c 460#include "xc_bnl.F" 461c 462c 463c --------------------------------------------------------------------------------------- 464c Calculate the third derivative of the main function that consititutes the functional 465c --------------------------------------------------------------------------------------- 466c 467 double precision function HqBNLThreePrime(q) 468 469 implicit none 470#include "xc_params.fh" 471 472 double precision q,OneOverQ 473 double precision q2, q3, q4 474 475 if (q .lt. 0.1d0) then 476 HqBNLThreePrime = -16.0d0*q 477 return 478 end if 479 480 OneOverQ = 1.0d0/q 481 q2 = q*q 482 q3 = q2*q 483 q4 = q3*q 484 485 HqBNLThreePrime = 8.0d0*( exp(-OneOverQ*OneOverQ) 486 1 + 2.0d0*q2*exp(-OneOverQ*OneOverQ) 487 2 - 2.0d0*q4 488 3 + 2.0d0*q4*exp(-OneOverQ*OneOverQ) ) / q3 489 490 return 491 end 492c 493c --------------------------------------------------------------------------------------- 494c Calculate the third derivative of the local Fermi vector 495c --------------------------------------------------------------------------------------- 496c 497 double precision function FermiKThreePrime(den) 498 499 implicit none 500#include "xc_params.fh" 501 502 double precision F13, den 503 504 F13 = 1.0D0/3.0D0 505 FermiKThreePrime = (10.0d0/9.0d0)* 506 1 (Pi**2.0d0/(9.0d0*den**8.0d0))**F13 507 508 return 509 end 510c 511c --------------------------------------------------------------------------------------- 512c Calculate the third derivative of q 513c --------------------------------------------------------------------------------------- 514c 515 double precision function QThreePrime(gamma,kF) 516 517 implicit none 518 519 double precision gamma, kF, FermiK4 520 521 FermiK4 = kF**4.0d0 522 QThreePrime = -(6.0d0*gamma)/FermiK4 523 524 return 525 end 526c 527c --------------------------------------------------------------------------------------- 528c Calculate the third derivative of EpsX 529c --------------------------------------------------------------------------------------- 530c 531 double precision function EpsThreeXprime(Rho,gamma) 532 533 implicit none 534#include "xc_params.fh" 535 536 double precision Rho, gamma 537c Fermi wavevector stuff 538 double precision kF, kFPrim, kF2prim, kF3prim 539 double precision kFPrim2, kFPrim3 540c q stuff 541 double precision q, qprim, q2prim, q3prim 542 double precision qprim2, qprim3 543c H(q) stuff 544 double precision g, gprim, g2prim, g3prim 545 double precision Cs 546 547 double precision FermiK 548 double precision FermiKPrime 549 double precision FermiKTwoPrime 550 double precision FermiKThreePrime 551 double precision QPrime 552 double precision QTwoPrime 553 double precision QThreePrime 554 double precision HqBNL 555 double precision HqBNLPrime 556 double precision HqBNLTwoPrime 557 double precision HqBNLThreePrime 558 559 if (Rho.le.0d0) then 560 EpsThreeXprime = 0.0d0 561 return 562 end if 563 564 kF = FermiK(Rho) 565 kFPrim = FermiKPrime(Rho) 566 kF2Prim = FermiKTwoPrime(Rho) 567 kF3Prim = FermiKThreePrime(Rho) 568c 569 kFPrim2 = kFPrim**2.0d0 570 kFPrim3 = kFPrim**3.0d0 571c 572 Cs = -3.0d0/(4.0d0*Pi) 573c 574 q = gamma / kF 575 qprim = QPrime(gamma,kF) 576 q2prim = QTwoPrime(gamma,kF) 577 q3prim = QThreePrime(gamma,kF) 578c 579 qprim2 = qprim**2.0d0 580 qprim3 = qprim**3.0d0 581c 582 g = HqBNL(q) 583 gprim = HqBNLPrime(q) 584 g2prim = HqBNLTwoPrime(q) 585 g3prim = HqBNLThreePrime(q) 586 587 EpsThreeXprime = Cs*kFPrim3*( 3.0d0*qprim2*g2prim 588 1 + 3.0d0*kF*qprim*g2prim*q2prim 589 2 + kF*qprim3*g3prim 590 3 + gprim*( 3.0d0*q2prim 591 4 + kF*q3prim ) ) 592 5 + 3.0d0*Cs*kFPrim*kF2Prim*( kF*qprim2*g2prim 593 6 + gprim*( 2.0d0*qprim 594 7 + kF*q2prim ) ) 595 8 + Cs*kF3Prim*( g + kF*gprim*qprim ) 596 597 return 598 end 599#endif 600c 601c $Id$ 602