1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Calculation of the incomplete Gamma function F_n(t) for multi-center 8!> integrals over Cartesian Gaussian functions. 9!> \par History 10!> - restructured and cleaned (24.05.2004,MK) 11!> \author Matthias Krack (07.01.1999) 12! ************************************************************************************************** 13MODULE gamma 14 15 USE kinds, ONLY: dp 16 USE mathconstants, ONLY: ifac,& 17 pi 18#include "../base/base_uses.f90" 19 20 IMPLICIT NONE 21 22 PRIVATE 23 24! *** Global parameters *** 25 26 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'gamma' 27 REAL(KIND=dp), PARAMETER :: teps = 1.0E-13_dp 28 29! *** Maximum n value of the tabulated F_n(t) function values *** 30 31 INTEGER, SAVE :: current_nmax = -1 32 33! *** F_n(t) table *** 34 35 REAL(KIND=dp), DIMENSION(:, :), ALLOCATABLE, SAVE :: ftable 36 37! *** Public subroutines *** 38 39 PUBLIC :: deallocate_md_ftable, & 40 fgamma_0, & 41 fgamma_ref, & 42 init_md_ftable 43 44CONTAINS 45 46! ************************************************************************************************** 47!> \brief Build a table of F_n(t) values in the range tmin <= t <= tmax 48!> with a stepsize of tdelta up to n equal to nmax. 49!> \param nmax ... 50!> \param tmin ... 51!> \param tmax ... 52!> \param tdelta ... 53!> \date 11.01.1999 54!> \par Parameters 55!> - nmax : Maximum n value of F_n(t). 56!> - tdelta: Difference between two consecutive t abcissas (step size). 57!> - tmax : Maximum t value. 58!> - tmin : Minimum t value. 59!> \author MK 60!> \version 1.0 61! ************************************************************************************************** 62 SUBROUTINE create_md_ftable(nmax, tmin, tmax, tdelta) 63 64 INTEGER, INTENT(IN) :: nmax 65 REAL(KIND=dp), INTENT(IN) :: tmin, tmax, tdelta 66 67 CHARACTER(len=*), PARAMETER :: routineN = 'create_md_ftable', & 68 routineP = moduleN//':'//routineN 69 70 INTEGER :: itab, itabmax, itabmin, n 71 REAL(KIND=dp) :: t 72 73 IF (current_nmax > -1) THEN 74 CALL cp_abort(__LOCATION__, & 75 "An incomplete Gamma function table is already "// & 76 " allocated. Use the init routine for an update") 77 END IF 78 79 IF (nmax < 0) THEN 80 CALL cp_abort(__LOCATION__, & 81 "A negative n value for the initialization of the "// & 82 "incomplete Gamma function is invalid") 83 END IF 84 85! *** Check arguments *** 86 87 IF ((tmax < 0.0_dp) .OR. & 88 (tmin < 0.0_dp) .OR. & 89 (tdelta <= 0.0_dp) .OR. & 90 (tmin > tmax)) THEN 91 CPABORT("Invalid arguments") 92 END IF 93 94 n = nmax + 6 95 96 itabmin = FLOOR(tmin/tdelta) 97 itabmax = CEILING((tmax - tmin)/tdelta) 98 99 ALLOCATE (ftable(0:n, itabmin:itabmax)) 100 ftable = 0.0_dp 101 102! *** Fill table *** 103 104 DO itab = itabmin, itabmax 105 t = REAL(itab, dp)*tdelta 106 ftable(0:n, itab) = fgamma_ref(n, t) 107 END DO 108 109! *** Save initialization status *** 110 111 current_nmax = nmax 112 113 END SUBROUTINE create_md_ftable 114 115! ************************************************************************************************** 116!> \brief Deallocate the table of F_n(t) values. 117!> \date 24.05.2004 118!> \author MK 119!> \version 1.0 120! ************************************************************************************************** 121 SUBROUTINE deallocate_md_ftable() 122 123 CHARACTER(len=*), PARAMETER :: routineN = 'deallocate_md_ftable', & 124 routineP = moduleN//':'//routineN 125 126 IF (current_nmax > -1) THEN 127 128 DEALLOCATE (ftable) 129 130 current_nmax = -1 131 132 END IF 133 134 END SUBROUTINE deallocate_md_ftable 135 136! ************************************************************************************************** 137!> \brief Calculation of the incomplete Gamma function F(t) for multicenter 138!> integrals over Gaussian functions. f returns a vector with all 139!> F_n(t) values for 0 <= n <= nmax. 140!> \param nmax ... 141!> \param t ... 142!> \param f ... 143!> \date 08.01.1999, 144!> \par History 145!> 09.06.1999, MK : Changed from a FUNCTION to a SUBROUTINE 146!> \par Literature 147!> L. E. McMurchie, E. R. Davidson, J. Comp. Phys. 26, 218 (1978) 148!> \par Parameters 149!> - f : The incomplete Gamma function F_n(t). 150!> - nmax: Maximum n value of F_n(t). 151!> - t : Argument of the incomplete Gamma function. 152!> - kmax: Maximum number of iterations. 153!> - expt: Exponential term in the upward recursion of F_n(t). 154!> \author MK 155!> \version 1.0 156! ************************************************************************************************** 157 SUBROUTINE fgamma_0(nmax, t, f) 158 159 INTEGER, INTENT(IN) :: nmax 160 REAL(KIND=dp), INTENT(IN) :: t 161 REAL(KIND=dp), DIMENSION(0:nmax), INTENT(OUT) :: f 162 163 INTEGER :: itab, k, n 164 REAL(KIND=dp) :: expt, g, tdelta, tmp, ttab 165 166! *** Calculate F(t) *** 167 168 IF (t < teps) THEN 169 170! *** Special cases: t = 0 *** 171 172 DO n = 0, nmax 173 f(n) = 1.0_dp/REAL(2*n + 1, dp) 174 END DO 175 176 ELSE IF (t <= 12.0_dp) THEN 177 178! *** 0 < t < 12 -> Taylor expansion *** 179 180 tdelta = 0.1_dp 181 182! *** Pretabulation of the F_n(t) function *** 183! *** for the Taylor series expansion *** 184 185 IF (nmax > current_nmax) THEN 186 CALL init_md_ftable(nmax) 187 END IF 188 189 itab = NINT(t/tdelta) 190 ttab = REAL(itab, dp)*tdelta 191 192 f(nmax) = ftable(nmax, itab) 193 194 tmp = 1.0_dp 195 DO k = 1, 6 196 tmp = tmp*(ttab - t) 197 f(nmax) = f(nmax) + ftable(nmax + k, itab)*tmp*ifac(k) 198 END DO 199 200 expt = EXP(-t) 201 202! *** Use the downward recursion relation to *** 203! *** generate the remaining F_n(t) values *** 204 205 DO n = nmax - 1, 0, -1 206 f(n) = (2.0_dp*t*f(n + 1) + expt)/REAL(2*n + 1, dp) 207 END DO 208 209 ELSE 210 211! *** t > 12 *** 212 213 IF (t <= 15.0_dp) THEN 214 215! *** 12 < t <= 15 -> Four term polynom expansion *** 216 217 g = 0.4999489092_dp - 0.2473631686_dp/t + & 218 0.321180909_dp/t**2 - 0.3811559346_dp/t**3 219 f(0) = 0.5_dp*SQRT(pi/t) - g*EXP(-t)/t 220 221 ELSE IF (t <= 18.0_dp) THEN 222 223! *** 15 < t <= 18 -> Three term polynom expansion *** 224 225 g = 0.4998436875_dp - 0.24249438_dp/t + 0.24642845_dp/t**2 226 f(0) = 0.5_dp*SQRT(pi/t) - g*EXP(-t)/t 227 228 ELSE IF (t <= 24.0_dp) THEN 229 230! *** 18 < t <= 24 -> Two term polynom expansion *** 231 232 g = 0.499093162_dp - 0.2152832_dp/t 233 f(0) = 0.5_dp*SQRT(pi/t) - g*EXP(-t)/t 234 235 ELSE IF (t <= 30.0_dp) THEN 236 237! *** 24 < t <= 30 -> One term polynom expansion *** 238 239 g = 0.49_dp 240 f(0) = 0.5_dp*SQRT(pi/t) - g*EXP(-t)/t 241 242 ELSE 243 244! *** t > 30 -> Asymptotic formula *** 245 246 f(0) = 0.5_dp*SQRT(pi/t) 247 248 END IF 249 250 IF (t > REAL(2*nmax + 36, dp)) THEN 251 expt = 0.0_dp 252 ELSE 253 expt = EXP(-t) 254 END IF 255 256! *** Use the upward recursion relation to *** 257! *** generate the remaining F_n(t) values *** 258 259 DO n = 1, nmax 260 f(n) = 0.5_dp*(REAL(2*n - 1, dp)*f(n - 1) - expt)/t 261 END DO 262 263 END IF 264 265 END SUBROUTINE fgamma_0 266 267! ************************************************************************************************** 268!> \brief Calculation of the incomplete Gamma function F(t) for multicenter 269!> integrals over Gaussian functions. f returns a vector with all 270!> F_n(t) values for 0 <= n <= nmax. 271!> \param nmax ... 272!> \param t ... 273!> \param f ... 274!> \date 08.01.1999 275!> \par Literature 276!> L. E. McMurchie, E. R. Davidson, J. Comp. Phys. 26, 218 (1978) 277!> \par Parameters 278!> - f : The incomplete Gamma function F_n(t). 279!> - nmax: Maximum n value of F_n(t). 280!> - t : Argument of the incomplete Gamma function. 281!> \author MK 282!> \version 1.0 283! ************************************************************************************************** 284 SUBROUTINE fgamma_1(nmax, t, f) 285 286 INTEGER, INTENT(IN) :: nmax 287 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: t 288 REAL(KIND=dp), DIMENSION(SIZE(t, 1), 0:nmax), & 289 INTENT(OUT) :: f 290 291 INTEGER :: i, itab, k, n 292 REAL(KIND=dp) :: expt, g, tdelta, tmp, ttab 293 294 DO i = 1, SIZE(t, 1) 295 296! *** Calculate F(t) *** 297 298 IF (t(i) < teps) THEN 299 300! *** Special cases: t = 0 *** 301 302 DO n = 0, nmax 303 f(i, n) = 1.0_dp/REAL(2*n + 1, dp) 304 END DO 305 306 ELSE IF (t(i) <= 12.0_dp) THEN 307 308! *** 0 < t < 12 -> Taylor expansion *** 309 310 tdelta = 0.1_dp 311 312! *** Pretabulation of the F_n(t) function *** 313! *** for the Taylor series expansion *** 314 315 IF (nmax > current_nmax) THEN 316 CALL init_md_ftable(nmax) 317 END IF 318 319 itab = NINT(t(i)/tdelta) 320 ttab = REAL(itab, dp)*tdelta 321 322 f(i, nmax) = ftable(nmax, itab) 323 324 tmp = 1.0_dp 325 DO k = 1, 6 326 tmp = tmp*(ttab - t(i)) 327 f(i, nmax) = f(i, nmax) + ftable(nmax + k, itab)*tmp*ifac(k) 328 END DO 329 330 expt = EXP(-t(i)) 331 332! *** Use the downward recursion relation to *** 333! *** generate the remaining F_n(t) values *** 334 335 DO n = nmax - 1, 0, -1 336 f(i, n) = (2.0_dp*t(i)*f(i, n + 1) + expt)/REAL(2*n + 1, dp) 337 END DO 338 339 ELSE 340 341! *** t > 12 *** 342 343 IF (t(i) <= 15.0_dp) THEN 344 345! *** 12 < t <= 15 -> Four term polynom expansion *** 346 347 g = 0.4999489092_dp - 0.2473631686_dp/t(i) + & 348 0.321180909_dp/t(i)**2 - 0.3811559346_dp/t(i)**3 349 f(i, 0) = 0.5_dp*SQRT(pi/t(i)) - g*EXP(-t(i))/t(i) 350 351 ELSE IF (t(i) <= 18.0_dp) THEN 352 353! *** 15 < t <= 18 -> Three term polynom expansion *** 354 355 g = 0.4998436875_dp - 0.24249438_dp/t(i) + 0.24642845_dp/t(i)**2 356 f(i, 0) = 0.5_dp*SQRT(pi/t(i)) - g*EXP(-t(i))/t(i) 357 358 ELSE IF (t(i) <= 24.0_dp) THEN 359 360! *** 18 < t <= 24 -> Two term polynom expansion *** 361 362 g = 0.499093162_dp - 0.2152832_dp/t(i) 363 f(i, 0) = 0.5_dp*SQRT(pi/t(i)) - g*EXP(-t(i))/t(i) 364 365 ELSE IF (t(i) <= 30.0_dp) THEN 366 367! *** 24 < t <= 30 -> One term polynom expansion *** 368 369 g = 0.49_dp 370 f(i, 0) = 0.5_dp*SQRT(pi/t(i)) - g*EXP(-t(i))/t(i) 371 372 ELSE 373 374! *** t > 30 -> Asymptotic formula *** 375 376 f(i, 0) = 0.5_dp*SQRT(pi/t(i)) 377 378 END IF 379 380 IF (t(i) > REAL(2*nmax + 36, dp)) THEN 381 expt = 0.0_dp 382 ELSE 383 expt = EXP(-t(i)) 384 END IF 385 386! *** Use the upward recursion relation to *** 387! *** generate the remaining F_n(t) values *** 388 389 DO n = 1, nmax 390 f(i, n) = 0.5_dp*(REAL(2*n - 1, dp)*f(i, n - 1) - expt)/t(i) 391 END DO 392 393 END IF 394 395 END DO 396 397 END SUBROUTINE fgamma_1 398 399! ************************************************************************************************** 400!> \brief Calculation of the incomplete Gamma function F_n(t) using a 401!> spherical Bessel function expansion. fgamma_ref returns a 402!> vector with all F_n(t) values for 0 <= n <= nmax. 403!> For t values greater than 50 an asymptotic formula is used. 404!> This function is expected to return accurate F_n(t) values 405!> for any combination of n and t, but the calculation is slow 406!> and therefore the function may only be used for a pretabulation 407!> of F_n(t) values or for reference calculations. 408!> \param nmax ... 409!> \param t ... 410!> \return ... 411!> \date 07.01.1999 412!> \par Literature 413!> F. E. Harris, Int. J. Quant. Chem. 23, 1469 (1983) 414!> \par Parameters 415!> - expt : Exponential term in the downward recursion of F_n(t). 416!> - factor : Prefactor of the Bessel function expansion. 417!> - nmax : Maximum n value of F_n(t). 418!> - p : Product of the Bessel function quotients. 419!> - r : Quotients of the Bessel functions. 420!> - sumterm: One term in the sum over products of Bessel functions. 421!> - t : Argument of the incomplete Gamma function. 422!> \author MK 423!> \version 1.0 424! ************************************************************************************************** 425 FUNCTION fgamma_ref(nmax, t) RESULT(f) 426 427 INTEGER, INTENT(IN) :: nmax 428 REAL(KIND=dp), INTENT(IN) :: t 429 REAL(KIND=dp), DIMENSION(0:nmax) :: f 430 431 CHARACTER(len=*), PARAMETER :: routineN = 'fgamma_ref', routineP = moduleN//':'//routineN 432 INTEGER, PARAMETER :: kmax = 50 433 REAL(KIND=dp), PARAMETER :: eps = EPSILON(0.0_dp) 434 435 INTEGER :: j, k, n 436 REAL(KIND=dp) :: expt, factor, p, sumterm, sumtot, term 437 REAL(KIND=dp), DIMENSION(kmax+10) :: r 438 439! ------------------------------------------------------------------ 440! *** Initialization *** 441 442 f(:) = 0.0_dp 443 444 IF (t < teps) THEN 445 446! *** Special case: t = 0 => analytic expression *** 447 448 DO n = 0, nmax 449 f(n) = 1.0_dp/REAL(2*n + 1, dp) 450 END DO 451 452 ELSE IF (t <= 50.0_dp) THEN 453 454! *** Initialize ratios of Bessel functions *** 455 456 r(kmax + 10) = 0.0_dp 457 458 DO j = kmax + 9, 1, -1 459 r(j) = -t/(REAL(4*j + 2, dp) - t*r(j + 1)) 460 END DO 461 462 factor = 2.0_dp*SINH(0.5_dp*t)*EXP(-0.5_dp*t)/t 463 464 DO n = 0, nmax 465 466! *** Initialize iteration *** 467 468 sumtot = factor/REAL(2*n + 1, dp) 469 term = 1.0_dp 470 471! *** Begin the summation and recursion *** 472 473 DO k = 1, kmax 474 475 term = term*REAL(2*n - 2*k + 1, dp)/REAL(2*n + 2*k + 1, dp) 476 477! *** Product of Bessel function quotients *** 478 479 p = 1.0_dp 480 481 DO j = 1, k 482 p = p*r(j) 483 END DO 484 485 sumterm = factor*term*p*REAL(2*k + 1, dp)/REAL(2*n + 1, dp) 486 487 IF (ABS(sumterm) < eps) THEN 488 489! *** Iteration converged *** 490 491 EXIT 492 493 ELSE IF (k == kmax) THEN 494 495! *** No convergence with kmax iterations *** 496 497 CPABORT("Maximum number of iterations reached") 498 499 ELSE 500 501! *** Add the current term to the sum and continue the iteration *** 502 503 sumtot = sumtot + sumterm 504 505 END IF 506 507 END DO 508 509 f(n) = sumtot 510 511 END DO 512 513 ELSE 514 515! *** Use asymptotic formula for t > 50 *** 516 517 f(0) = 0.5_dp*SQRT(pi/t) 518 519! *** Use the upward recursion relation to *** 520! *** generate the remaining F_n(t) values *** 521 522 expt = EXP(-t) 523 524 DO n = 1, nmax 525 f(n) = 0.5_dp*(REAL(2*n - 1, dp)*f(n - 1) - expt)/t 526 END DO 527 528 END IF 529 530 END FUNCTION fgamma_ref 531 532! ************************************************************************************************** 533!> \brief Initalize a table of F_n(t) values in the range 0 <= t <= 12 with 534!> a stepsize of 0.1 up to n equal to nmax for the Taylor series 535!> expansion used by McMurchie-Davidson (MD). 536!> \param nmax ... 537!> \date 10.06.1999 538!> \par Parameters 539!> - nmax : Maximum n value of F_n(t). 540!> \author MK 541!> \version 1.0 542! ************************************************************************************************** 543 SUBROUTINE init_md_ftable(nmax) 544 INTEGER, INTENT(IN) :: nmax 545 546 CHARACTER(len=*), PARAMETER :: routineN = 'init_md_ftable', routineP = moduleN//':'//routineN 547 548 IF (nmax < 0) THEN 549 CALL cp_abort(__LOCATION__, & 550 "A negative n value for the initialization of the "// & 551 "incomplete Gamma function is invalid") 552 END IF 553 554! *** Check, if the current initialization is sufficient *** 555 556 IF (nmax > current_nmax) THEN 557 558 CALL deallocate_md_ftable() 559 560! *** Pretabulation of the F_n(t) function *** 561! *** for the Taylor series expansion *** 562 563 CALL create_md_ftable(nmax, 0.0_dp, 12.0_dp, 0.1_dp) 564 565 END IF 566 567 END SUBROUTINE init_md_ftable 568 569END MODULE gamma 570