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