(git:374b731)
Loading...
Searching...
No Matches
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! **************************************************************************************************
12MODULE 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
29CONTAINS
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
489END MODULE whittaker
490
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::dfac
real(kind=dp), dimension(-1:2 *maxfac+1), parameter, public dfac
Definition mathconstants.F:61
mathconstants::rootpi
real(kind=dp), parameter, public rootpi
Definition mathconstants.F:114
mathconstants::fac
real(kind=dp), dimension(0:maxfac), parameter, public fac
Definition mathconstants.F:37
whittaker
Calculates special integrals.
Definition whittaker.F:12
whittaker::epsilon
real(kind=dp), parameter epsilon
Definition whittaker.F:23
whittaker::whittaker_c0a
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
whittaker::whittaker_ci
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