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 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 3927574 : 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 3927574 : l = l1 + l2
61 :
62 3927574 : IF (MOD(l, 2) /= 0) THEN
63 0 : CPABORT("Total angular momentum has to be even")
64 : END IF
65 3927574 : IF (l1 < l2) THEN
66 0 : CPABORT("l1 >= l2")
67 : END IF
68 :
69 216643514 : wc(:) = 0.0_dp
70 3927574 : t1 = SQRT(alpha)
71 3927574 : y = REAL(l, dp)
72 216643514 : DO i = 1, n
73 212715940 : x = r(i)
74 216643514 : 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 28271426 : alpha**5*x**12/120._dp/(13._dp + y))
80 : ELSE
81 184444514 : wc(i) = -rootpi*erfa(i)*alpha*dfac(l + 1)
82 683871668 : DO k = 0, l/2
83 : wc(i) = wc(i) + expa(i)*x**(2*k + 1)*t1**(2*k + 3)* &
84 683871668 : dfac(l + 1)/dfac(2*k + 1)*2**(k + 1)
85 : END DO
86 184444514 : wc(i) = -wc(i)/2._dp**(l/2 + 2)/t1**(l + 5)/x**(l2 + 1)
87 : END IF
88 : END DO
89 :
90 3927574 : 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 0 : 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 0 : IF (MOD(l, 2) /= 0) THEN
123 0 : CPABORT("Angular momentum has to be even")
124 : END IF
125 :
126 0 : wc(:) = 0.0_dp
127 :
128 0 : SELECT CASE (l)
129 :
130 : CASE DEFAULT
131 :
132 0 : t1 = SQRT(alpha)
133 0 : DO i = 1, n
134 0 : x = r(i)
135 0 : wc(i) = -rootpi*erfa(i)*alpha*dfac(l + 1)
136 0 : DO k = 0, l/2
137 : wc(i) = wc(i) + expa(i)*x**(2*k + 1)*t1**(2*k + 3)* &
138 0 : dfac(l + 1)/dfac(2*k + 1)*2**(k + 1)
139 : END DO
140 0 : wc(i) = -wc(i)/2._dp**(l/2 + 2)/t1**(l + 5)
141 : END DO
142 :
143 : CASE (0)
144 :
145 0 : t1 = SQRT(alpha)
146 0 : t2 = t1**2
147 0 : t11 = rootpi
148 0 : DO i = 1, n
149 0 : x = r(i)
150 0 : t5 = x**2
151 0 : t7 = expa(i)
152 0 : t13 = erfa(i)
153 0 : t18 = -1._dp/t2/t1*(2._dp*x*t7*t1 - t11*t13)/4._dp
154 0 : wc(i) = t18
155 : END DO
156 :
157 : CASE (2)
158 :
159 0 : t1 = SQRT(alpha)
160 0 : t2 = t1**2
161 0 : t3 = t2**2
162 0 : t17 = rootpi
163 0 : DO i = 1, n
164 0 : x = r(i)
165 0 : t6 = x**2
166 0 : t9 = expa(i)
167 0 : t19 = erfa(i)
168 0 : t25 = -1._dp/t3/t1*(4._dp*t6*x*t9*t2*t1 + 6._dp*x*t9*t1 - 3*t17*t19)/8._dp
169 0 : wc(i) = t25
170 : END DO
171 :
172 : CASE (4)
173 :
174 0 : t1 = SQRT(alpha)
175 0 : t2 = t1**2
176 0 : t3 = t2*t1
177 0 : t4 = t2**2
178 0 : t23 = rootpi
179 0 : DO i = 1, n
180 0 : x = r(i)
181 0 : t7 = x**2
182 0 : t8 = t7**2
183 0 : t11 = expa(i)
184 0 : 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 0 : 15._dp*t23*t25)/16._dp
187 0 : wc(i) = t31
188 : END DO
189 :
190 : CASE (6)
191 :
192 0 : t8 = SQRT(alpha)
193 0 : t9 = t8**2
194 0 : t10 = t9**2
195 0 : t11 = t10**2
196 0 : t17 = t9*t8
197 0 : t28 = rootpi
198 0 : DO i = 1, n
199 0 : x = r(i)
200 0 : t1 = x**2
201 0 : t2 = t1*x
202 0 : t3 = t1**2
203 0 : t6 = expa(i)
204 0 : 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 0 : 210._dp*x*t6*t17 - 105._dp*t28*t30*alpha)/t11/t17/32._dp
207 0 : wc(i) = t39
208 : END DO
209 :
210 : CASE (8)
211 :
212 0 : t8 = SQRT(alpha)
213 0 : t9 = t8**2
214 0 : t10 = t9*t8
215 0 : t11 = t9**2
216 0 : t12 = t11**2
217 0 : t34 = rootpi
218 0 : DO i = 1, n
219 0 : x = r(i)
220 0 : t1 = x**2
221 0 : t2 = t1**2
222 0 : t3 = t2**2
223 0 : t6 = expa(i)
224 0 : t16 = t1*x
225 0 : t28 = t11*t8
226 0 : 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 0 : 1260._dp*t16*t6*t28 + 1890._dp*x*t6*t10 - 945._dp*t34*t36*alpha)/t12/t28/64._dp
229 0 : wc(i) = t45
230 : END DO
231 :
232 : CASE (10)
233 :
234 0 : t9 = SQRT(alpha)
235 0 : t10 = t9**2
236 0 : t11 = t10**2
237 0 : t12 = t11*t9
238 0 : t13 = t11**2
239 0 : t19 = t10*t9
240 0 : t30 = t11*t19
241 0 : t39 = rootpi
242 0 : DO i = 1, n
243 0 : x = r(i)
244 0 : t1 = x**2
245 0 : t2 = t1*x
246 0 : t3 = t1**2
247 0 : t4 = t3**2
248 0 : t7 = expa(i)
249 0 : 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 0 : t13/t30/128._dp
254 0 : wc(i) = t50
255 : END DO
256 :
257 : CASE (12)
258 :
259 0 : t9 = SQRT(alpha)
260 0 : t10 = t9**2
261 0 : t11 = t10*t9
262 0 : t12 = t10**2
263 0 : t13 = t12*t11
264 0 : t14 = t12**2
265 0 : t44 = rootpi
266 0 : DO i = 1, n
267 0 : x = r(i)
268 0 : t1 = x**2
269 0 : t2 = t1**2
270 0 : t3 = t2*x
271 0 : t4 = t2**2
272 0 : t7 = expa(i)
273 0 : t18 = t1*x
274 0 : t21 = t12*t9
275 0 : t46 = erfa(i)
276 0 : 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 0 : t51/t9/256._dp
281 0 : wc(i) = t56
282 : END DO
283 :
284 : CASE (14)
285 :
286 0 : t10 = SQRT(alpha)
287 0 : t11 = t10**2
288 0 : t12 = t11**2
289 0 : t13 = t12**2
290 0 : t14 = t13**2
291 0 : t21 = t11*t10
292 0 : t22 = t12*t21
293 0 : t28 = t12*t10
294 0 : t50 = rootpi
295 0 : DO i = 1, n
296 0 : x = r(i)
297 0 : t1 = x**2
298 0 : t2 = t1*x
299 0 : t3 = t1**2
300 0 : t4 = t3*t2
301 0 : t5 = t3**2
302 0 : t8 = expa(i)
303 0 : t18 = t3*x
304 0 : 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 0 : 2027025._dp*t50*t52*alpha)/t14/t21/512._dp
310 0 : wc(i) = t61
311 : END DO
312 :
313 : END SELECT
314 :
315 0 : 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 3920130 : 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 3920130 : IF (MOD(l, 2) /= 0) THEN
351 0 : CPABORT("Angular momentum has to be even")
352 : END IF
353 :
354 216254910 : wc(:) = 0.0_dp
355 :
356 : SELECT CASE (l)
357 :
358 : CASE DEFAULT
359 :
360 0 : DO i = 1, n
361 0 : x = r(i)
362 0 : wc(i) = 0._dp
363 0 : DO k = 0, l/2
364 0 : wc(i) = wc(i) + alpha**k*x**(2*k)*fac(l/2)/fac(k)
365 : END DO
366 0 : wc(i) = 0.5_dp*wc(i)/alpha**(l/2 + 1)*expa(i)
367 : END DO
368 :
369 : CASE (0)
370 :
371 152070104 : DO i = 1, n
372 149335140 : x = r(i)
373 149335140 : t2 = x**2
374 149335140 : t4 = expa(i)
375 149335140 : t6 = 1._dp/alpha*t4/2._dp
376 152070104 : wc(i) = t6
377 : END DO
378 :
379 : CASE (2)
380 :
381 1064258 : t6 = alpha**2
382 57936298 : DO i = 1, n
383 56872040 : x = r(i)
384 56872040 : t1 = x**2
385 56872040 : t2 = alpha*t1
386 56872040 : t3 = expa(i)
387 56872040 : t9 = t3*(t2 + 1)/t6/2._dp
388 57936298 : wc(i) = t9
389 : END DO
390 :
391 : CASE (4)
392 :
393 114958 : t5 = alpha**2
394 5945118 : DO i = 1, n
395 5830160 : x = r(i)
396 5830160 : t1 = x**2
397 5830160 : t2 = alpha*t1
398 5830160 : t3 = expa(i)
399 5830160 : t4 = t1**2
400 5830160 : t13 = t3*(t4*t5 + 2._dp*t2 + 2._dp)/t5/alpha/2._dp
401 5945118 : wc(i) = t13
402 : END DO
403 :
404 : CASE (6)
405 :
406 5882 : t6 = alpha**2
407 5882 : t14 = t6**2
408 299922 : DO i = 1, n
409 294040 : x = r(i)
410 294040 : t1 = x**2
411 294040 : t2 = alpha*t1
412 294040 : t3 = expa(i)
413 294040 : t4 = t1**2
414 294040 : t17 = t3*(t4*t1*t6*alpha + 3._dp*t4*t6 + 6._dp*t2 + 6._dp)/t14/2._dp
415 299922 : wc(i) = t17
416 : END DO
417 :
418 : CASE (8)
419 :
420 60 : t6 = alpha**2
421 60 : t7 = t6**2
422 3060 : DO i = 1, n
423 3000 : x = r(i)
424 3000 : t1 = x**2
425 3000 : t2 = alpha*t1
426 3000 : t3 = expa(i)
427 3000 : t4 = t1**2
428 3000 : t5 = t4**2
429 3000 : t21 = t3*(t5*t7 + 4._dp*t4*t1*t6*alpha + 12._dp*t4*t6 + 24._dp*t2 + 24._dp)/t7/alpha/2._dp
430 3060 : wc(i) = t21
431 : END DO
432 :
433 : CASE (10)
434 :
435 8 : t7 = alpha**2
436 8 : t8 = t7**2
437 408 : DO i = 1, n
438 400 : x = r(i)
439 400 : t1 = x**2
440 400 : t2 = alpha*t1
441 400 : t3 = expa(i)
442 400 : t4 = t1**2
443 400 : 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 400 : 120._dp*t2 + 120._dp)/t8/t7/2._dp
446 408 : wc(i) = t25
447 : END DO
448 :
449 : CASE (12)
450 :
451 0 : t7 = alpha**2
452 0 : t8 = t7**2
453 0 : t18 = t7*alpha
454 0 : DO i = 1, n
455 0 : x = r(i)
456 0 : t1 = x**2
457 0 : t2 = alpha*t1
458 0 : t3 = expa(i)
459 0 : t4 = t1**2
460 0 : 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 0 : 360._dp*t4*t7 + 720._dp*t2 + 720._dp)/t8/t18/2._dp
463 0 : wc(i) = t29
464 : END DO
465 :
466 : CASE (14)
467 :
468 0 : t8 = alpha**2
469 0 : t9 = t8*alpha
470 0 : t10 = t8**2
471 0 : t30 = t10**2
472 3920130 : DO i = 1, n
473 0 : x = r(i)
474 0 : t1 = x**2
475 0 : t2 = alpha*t1
476 0 : t3 = expa(i)
477 0 : t4 = t1**2
478 0 : t5 = t4*t1
479 0 : t6 = t4**2
480 : t33 = t3*(t6*t5*t10*t9 + 7*t6*t4*t10*t8 + 42._dp*t6*t1*t10*alpha + &
481 0 : 210._dp*t6*t10 + 840._dp*t5*t9 + 2520._dp*t4*t8 + 5040._dp*t2 + 5040._dp)/t30/2._dp
482 0 : wc(i) = t33
483 : END DO
484 :
485 : END SELECT
486 :
487 3920130 : END SUBROUTINE whittaker_ci
488 :
489 : END MODULE whittaker
490 :
|
