(git:0de0cc2)
whittaker.F
Go to the documentation of this file.
1 !--------------------------------------------------------------------------------------------------!
2 ! CP2K: A general program to perform molecular dynamics simulations !
3 ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 ! !
5 ! SPDX-License-Identifier: GPL-2.0-or-later !
6 !--------------------------------------------------------------------------------------------------!
7 
8 ! **************************************************************************************************
9 !> \brief Calculates special integrals
10 !> \author JGH 10-08-2004
11 ! **************************************************************************************************
12 MODULE whittaker
13 
14  USE kinds, ONLY: dp
15  USE mathconstants, ONLY: dfac,&
16  fac,&
17  rootpi
18 #include "../base/base_uses.f90"
19 
20  IMPLICIT NONE
21 
22  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'whittaker'
23  REAL(kind=dp), PARAMETER :: epsilon = 1.e-2_dp
24 
25  PRIVATE
26 
27  PUBLIC :: whittaker_c0a, whittaker_ci
28 
29 CONTAINS
30 
31 ! **************************************************************************************************
32 !> \brief int(y^(2+l1+l2) * exp(-alpha*y*y),y=0..x) / x^(l2+1);
33 !> wc(:) :: output
34 !> r(:) :: coordinate
35 !> expa(:) :: exp(-alpha*r(:)**2)
36 !> erfa(:) :: erf(sqrt(alpha)*r(:))
37 !> alpha :: exponent
38 !> l1, l2 :: L-quantum number
39 !> n :: number of points
40 !>
41 !> \param wc ...
42 !> \param r ...
43 !> \param expa ...
44 !> \param erfa ...
45 !> \param alpha ...
46 !> \param l1 ...
47 !> \param l2 ...
48 !> \param n ...
49 !> \author JGH 10-08-2004
50 ! **************************************************************************************************
51  SUBROUTINE whittaker_c0a(wc, r, expa, erfa, alpha, l1, l2, n)
52  INTEGER, INTENT(IN) :: n, l2, l1
53  REAL(kind=dp), INTENT(IN) :: alpha
54  REAL(kind=dp), DIMENSION(n), INTENT(IN) :: erfa, expa, r
55  REAL(kind=dp), DIMENSION(n), INTENT(OUT) :: wc
56 
57  INTEGER :: i, k, l
58  REAL(dp) :: t1, x, y
59 
60  l = l1 + l2
61 
62  IF (mod(l, 2) /= 0) THEN
63  cpabort("Total angular momentum has to be even")
64  END IF
65  IF (l1 < l2) THEN
66  cpabort("l1 >= l2")
67  END IF
68 
69  wc(:) = 0.0_dp
70  t1 = sqrt(alpha)
71  y = real(l, dp)
72  DO i = 1, n
73  x = r(i)
74  IF (t1*x < epsilon) THEN
75  wc(i) = x**l1*(x**2/(3._dp + y) - alpha*x**4/(5._dp + y) + &
76  alpha**2*x**6/(14._dp + 2._dp*y) - &
77  alpha**3*x**8/(54._dp + 6._dp*y) + &
78  alpha**4*x**10/(256._dp + 24._dp*y) - &
79  alpha**5*x**12/120._dp/(13._dp + y))
80  ELSE
81  wc(i) = -rootpi*erfa(i)*alpha*dfac(l + 1)
82  DO k = 0, l/2
83  wc(i) = wc(i) + expa(i)*x**(2*k + 1)*t1**(2*k + 3)* &
84  dfac(l + 1)/dfac(2*k + 1)*2**(k + 1)
85  END DO
86  wc(i) = -wc(i)/2._dp**(l/2 + 2)/t1**(l + 5)/x**(l2 + 1)
87  END IF
88  END DO
89 
90  END SUBROUTINE whittaker_c0a
91 
92 ! **************************************************************************************************
93 !> \brief int(y^(2+l) * exp(-alpha*y*y),y=0..x);
94 !> wc(:) :: output
95 !> r(:) :: coordinate
96 !> expa(:) :: exp(-alpha*r(:)**2)
97 !> erfa(:) :: erf(sqrt(alpha)*r(:))
98 !> alpha :: exponent
99 !> l :: L-quantum number
100 !> n :: number of points
101 !>
102 !> \param wc ...
103 !> \param r ...
104 !> \param expa ...
105 !> \param erfa ...
106 !> \param alpha ...
107 !> \param l ...
108 !> \param n ...
109 !> \author JGH 10-08-2004
110 ! **************************************************************************************************
111  SUBROUTINE whittaker_c0(wc, r, expa, erfa, alpha, l, n)
112  INTEGER, INTENT(IN) :: n, l
113  REAL(kind=dp), INTENT(IN) :: alpha
114  REAL(kind=dp), DIMENSION(n), INTENT(IN) :: erfa, expa, r
115  REAL(kind=dp), DIMENSION(n), INTENT(OUT) :: wc
116 
117  INTEGER :: i, k
118  REAL(dp) :: t1, t10, t11, t12, t13, t14, t16, t17, t18, t19, t2, t21, t22, t23, t25, t28, &
119  t3, t30, t31, t34, t36, t39, t4, t41, t44, t45, t46, t5, t50, t51, t52, t56, t6, t61, t7, &
120  t8, t9, x
121 
122  IF (mod(l, 2) /= 0) THEN
123  cpabort("Angular momentum has to be even")
124  END IF
125 
126  wc(:) = 0.0_dp
127 
128  SELECT CASE (l)
129 
130  CASE DEFAULT
131 
132  t1 = sqrt(alpha)
133  DO i = 1, n
134  x = r(i)
135  wc(i) = -rootpi*erfa(i)*alpha*dfac(l + 1)
136  DO k = 0, l/2
137  wc(i) = wc(i) + expa(i)*x**(2*k + 1)*t1**(2*k + 3)* &
138  dfac(l + 1)/dfac(2*k + 1)*2**(k + 1)
139  END DO
140  wc(i) = -wc(i)/2._dp**(l/2 + 2)/t1**(l + 5)
141  END DO
142 
143  CASE (0)
144 
145  t1 = sqrt(alpha)
146  t2 = t1**2
147  t11 = rootpi
148  DO i = 1, n
149  x = r(i)
150  t5 = x**2
151  t7 = expa(i)
152  t13 = erfa(i)
153  t18 = -1._dp/t2/t1*(2._dp*x*t7*t1 - t11*t13)/4._dp
154  wc(i) = t18
155  END DO
156 
157  CASE (2)
158 
159  t1 = sqrt(alpha)
160  t2 = t1**2
161  t3 = t2**2
162  t17 = rootpi
163  DO i = 1, n
164  x = r(i)
165  t6 = x**2
166  t9 = expa(i)
167  t19 = erfa(i)
168  t25 = -1._dp/t3/t1*(4._dp*t6*x*t9*t2*t1 + 6._dp*x*t9*t1 - 3*t17*t19)/8._dp
169  wc(i) = t25
170  END DO
171 
172  CASE (4)
173 
174  t1 = sqrt(alpha)
175  t2 = t1**2
176  t3 = t2*t1
177  t4 = t2**2
178  t23 = rootpi
179  DO i = 1, n
180  x = r(i)
181  t7 = x**2
182  t8 = t7**2
183  t11 = expa(i)
184  t25 = erfa(i)
185  t31 = -1._dp/t4/t3*(8._dp*t8*x*t11*t4*t1 + 20._dp*t7*x*t11*t3 + 30._dp*x*t11*t1 - &
186  15._dp*t23*t25)/16._dp
187  wc(i) = t31
188  END DO
189 
190  CASE (6)
191 
192  t8 = sqrt(alpha)
193  t9 = t8**2
194  t10 = t9**2
195  t11 = t10**2
196  t17 = t9*t8
197  t28 = rootpi
198  DO i = 1, n
199  x = r(i)
200  t1 = x**2
201  t2 = t1*x
202  t3 = t1**2
203  t6 = expa(i)
204  t30 = erfa(i)
205  t39 = -(16._dp*t3*t2*t6*t11*t8 + 56._dp*t3*x*t6*t10*t17 + 140._dp*t2*t6*t10*t8 + &
206  210._dp*x*t6*t17 - 105._dp*t28*t30*alpha)/t11/t17/32._dp
207  wc(i) = t39
208  END DO
209 
210  CASE (8)
211 
212  t8 = sqrt(alpha)
213  t9 = t8**2
214  t10 = t9*t8
215  t11 = t9**2
216  t12 = t11**2
217  t34 = rootpi
218  DO i = 1, n
219  x = r(i)
220  t1 = x**2
221  t2 = t1**2
222  t3 = t2**2
223  t6 = expa(i)
224  t16 = t1*x
225  t28 = t11*t8
226  t36 = erfa(i)
227  t45 = -(32._dp*t3*x*t6*t12*t10 + 144._dp*t2*t16*t6*t12*t8 + 504._dp*t2*x*t6*t11*t10 + &
228  1260._dp*t16*t6*t28 + 1890._dp*x*t6*t10 - 945._dp*t34*t36*alpha)/t12/t28/64._dp
229  wc(i) = t45
230  END DO
231 
232  CASE (10)
233 
234  t9 = sqrt(alpha)
235  t10 = t9**2
236  t11 = t10**2
237  t12 = t11*t9
238  t13 = t11**2
239  t19 = t10*t9
240  t30 = t11*t19
241  t39 = rootpi
242  DO i = 1, n
243  x = r(i)
244  t1 = x**2
245  t2 = t1*x
246  t3 = t1**2
247  t4 = t3**2
248  t7 = expa(i)
249  t41 = erfa(i)
250  t50 = -(64._dp*t4*t2*t7*t13*t12 + 352._dp*t4*x*t7*t13*t19 + &
251  1584._dp*t3*t2*t7*t13*t9 + 5544._dp*t3*x*t7*t30 + &
252  13860._dp*t2*t7*t12 + 20790._dp*x*t7*t19 - 10395._dp*t39*t41*alpha)/ &
253  t13/t30/128._dp
254  wc(i) = t50
255  END DO
256 
257  CASE (12)
258 
259  t9 = sqrt(alpha)
260  t10 = t9**2
261  t11 = t10*t9
262  t12 = t10**2
263  t13 = t12*t11
264  t14 = t12**2
265  t44 = rootpi
266  DO i = 1, n
267  x = r(i)
268  t1 = x**2
269  t2 = t1**2
270  t3 = t2*x
271  t4 = t2**2
272  t7 = expa(i)
273  t18 = t1*x
274  t21 = t12*t9
275  t46 = erfa(i)
276  t51 = t14**2
277  t56 = -(128._dp*t4*t3*t7*t13*t14 + 832._dp*t4*t18*t7*t14*t21 + &
278  4576._dp*t4*x*t7*t14*t11 + 20592._dp*t2*t18*t7*t14*t9 + 72072._dp*t3*t7*t13 + &
279  180180._dp*t18*t7*t21 + 270270._dp*x*t7*t11 - 135135._dp*t44*t46*alpha)/ &
280  t51/t9/256._dp
281  wc(i) = t56
282  END DO
283 
284  CASE (14)
285 
286  t10 = sqrt(alpha)
287  t11 = t10**2
288  t12 = t11**2
289  t13 = t12**2
290  t14 = t13**2
291  t21 = t11*t10
292  t22 = t12*t21
293  t28 = t12*t10
294  t50 = rootpi
295  DO i = 1, n
296  x = r(i)
297  t1 = x**2
298  t2 = t1*x
299  t3 = t1**2
300  t4 = t3*t2
301  t5 = t3**2
302  t8 = expa(i)
303  t18 = t3*x
304  t52 = erfa(i)
305  t61 = -(256._dp*t5*t4*t8*t14*t10 + 1920._dp*t5*t18*t8*t13*t22 + &
306  12480._dp*t5*t2*t8*t13*t28 + 68640._dp*t5*x*t8*t13*t21 + &
307  308880._dp*t4*t8*t13*t10 + 1081080._dp*t18*t8*t22 + &
308  2702700._dp*t2*t8*t28 + 4054050._dp*x*t8*t21 - &
309  2027025._dp*t50*t52*alpha)/t14/t21/512._dp
310  wc(i) = t61
311  END DO
312 
313  END SELECT
314 
315  END SUBROUTINE whittaker_c0
316 
317 ! **************************************************************************************************
318 !> \brief int(y^(l+1) * exp(-alpha*y*y),y=x..infinity);
319 !>
320 !> (-1 - 1/2 l~) 2
321 !> 1/2 alpha GAMMA(1/2 l + 1, alpha x )
322 !>
323 !>
324 !> wc(:) :: output
325 !> r(:) :: coordinate
326 !> expa(:) :: exp(-alpha*r(:)**2)
327 !> alpha :: exponent
328 !> l :: L-quantum number
329 !> n :: number of points
330 !>
331 !> \param wc ...
332 !> \param r ...
333 !> \param expa ...
334 !> \param alpha ...
335 !> \param l ...
336 !> \param n ...
337 !> \author JGH 10-08-2004
338 ! **************************************************************************************************
339  SUBROUTINE whittaker_ci(wc, r, expa, alpha, l, n)
340  INTEGER, INTENT(IN) :: n, l
341  REAL(kind=dp), INTENT(IN) :: alpha
342  REAL(kind=dp), DIMENSION(n), INTENT(IN) :: expa, r
343  REAL(kind=dp), DIMENSION(n), INTENT(OUT) :: wc
344 
345  INTEGER :: i, k
346  REAL(dp) :: t1, t10, t13, t14, t17, t18, t2, t21, &
347  t25, t29, t3, t30, t33, t4, t5, t6, &
348  t7, t8, t9, x
349 
350  IF (mod(l, 2) /= 0) THEN
351  cpabort("Angular momentum has to be even")
352  END IF
353 
354  wc(:) = 0.0_dp
355 
356  SELECT CASE (l)
357 
358  CASE DEFAULT
359 
360  DO i = 1, n
361  x = r(i)
362  wc(i) = 0._dp
363  DO k = 0, l/2
364  wc(i) = wc(i) + alpha**k*x**(2*k)*fac(l/2)/fac(k)
365  END DO
366  wc(i) = 0.5_dp*wc(i)/alpha**(l/2 + 1)*expa(i)
367  END DO
368 
369  CASE (0)
370 
371  DO i = 1, n
372  x = r(i)
373  t2 = x**2
374  t4 = expa(i)
375  t6 = 1._dp/alpha*t4/2._dp
376  wc(i) = t6
377  END DO
378 
379  CASE (2)
380 
381  t6 = alpha**2
382  DO i = 1, n
383  x = r(i)
384  t1 = x**2
385  t2 = alpha*t1
386  t3 = expa(i)
387  t9 = t3*(t2 + 1)/t6/2._dp
388  wc(i) = t9
389  END DO
390 
391  CASE (4)
392 
393  t5 = alpha**2
394  DO i = 1, n
395  x = r(i)
396  t1 = x**2
397  t2 = alpha*t1
398  t3 = expa(i)
399  t4 = t1**2
400  t13 = t3*(t4*t5 + 2._dp*t2 + 2._dp)/t5/alpha/2._dp
401  wc(i) = t13
402  END DO
403 
404  CASE (6)
405 
406  t6 = alpha**2
407  t14 = t6**2
408  DO i = 1, n
409  x = r(i)
410  t1 = x**2
411  t2 = alpha*t1
412  t3 = expa(i)
413  t4 = t1**2
414  t17 = t3*(t4*t1*t6*alpha + 3._dp*t4*t6 + 6._dp*t2 + 6._dp)/t14/2._dp
415  wc(i) = t17
416  END DO
417 
418  CASE (8)
419 
420  t6 = alpha**2
421  t7 = t6**2
422  DO i = 1, n
423  x = r(i)
424  t1 = x**2
425  t2 = alpha*t1
426  t3 = expa(i)
427  t4 = t1**2
428  t5 = t4**2
429  t21 = t3*(t5*t7 + 4._dp*t4*t1*t6*alpha + 12._dp*t4*t6 + 24._dp*t2 + 24._dp)/t7/alpha/2._dp
430  wc(i) = t21
431  END DO
432 
433  CASE (10)
434 
435  t7 = alpha**2
436  t8 = t7**2
437  DO i = 1, n
438  x = r(i)
439  t1 = x**2
440  t2 = alpha*t1
441  t3 = expa(i)
442  t4 = t1**2
443  t5 = t4**2
444  t25 = t3*(t5*t1*t8*alpha + 5._dp*t5*t8 + 20._dp*t4*t1*t7*alpha + 60._dp*t4*t7 + &
445  120._dp*t2 + 120._dp)/t8/t7/2._dp
446  wc(i) = t25
447  END DO
448 
449  CASE (12)
450 
451  t7 = alpha**2
452  t8 = t7**2
453  t18 = t7*alpha
454  DO i = 1, n
455  x = r(i)
456  t1 = x**2
457  t2 = alpha*t1
458  t3 = expa(i)
459  t4 = t1**2
460  t5 = t4**2
461  t29 = t3*(t5*t4*t8*t7 + 6._dp*t5*t1*t8*alpha + 30._dp*t5*t8 + 120._dp*t4*t1*t18 + &
462  360._dp*t4*t7 + 720._dp*t2 + 720._dp)/t8/t18/2._dp
463  wc(i) = t29
464  END DO
465 
466  CASE (14)
467 
468  t8 = alpha**2
469  t9 = t8*alpha
470  t10 = t8**2
471  t30 = t10**2
472  DO i = 1, n
473  x = r(i)
474  t1 = x**2
475  t2 = alpha*t1
476  t3 = expa(i)
477  t4 = t1**2
478  t5 = t4*t1
479  t6 = t4**2
480  t33 = t3*(t6*t5*t10*t9 + 7*t6*t4*t10*t8 + 42._dp*t6*t1*t10*alpha + &
481  210._dp*t6*t10 + 840._dp*t5*t9 + 2520._dp*t4*t8 + 5040._dp*t2 + 5040._dp)/t30/2._dp
482  wc(i) = t33
483  END DO
484 
485  END SELECT
486 
487  END SUBROUTINE whittaker_ci
488 
489 END MODULE whittaker
490 
Defines the basic variable types.
Definition: kinds.F:23
integer, parameter, public dp
Definition: kinds.F:34
Definition of mathematical constants and functions.
Definition: mathconstants.F:16
real(kind=dp), dimension(-1:2 *maxfac+1), parameter, public dfac
Definition: mathconstants.F:61
real(kind=dp), parameter, public rootpi
real(kind=dp), dimension(0:maxfac), parameter, public fac
Definition: mathconstants.F:37
Calculates special integrals.
Definition: whittaker.F:12
real(kind=dp), parameter epsilon
Definition: whittaker.F:23
subroutine, public whittaker_c0a(wc, r, expa, erfa, alpha, l1, l2, n)
int(y^(2+l1+l2) * exp(-alpha*y*y),y=0..x) / x^(l2+1); wc(:) :: output r(:) :: coordinate expa(:) :: e...
Definition: whittaker.F:52
subroutine, public whittaker_ci(wc, r, expa, alpha, l, n)
int(y^(l+1) * exp(-alpha*y*y),y=x..infinity);
Definition: whittaker.F:340