Multidip 1.0
Multi-photon matrix elements
Loading...
Searching...
No Matches
multidip_tests.f90
Go to the documentation of this file.
1! Copyright 2020
2!
3! For a comprehensive list of the developers that contributed to these codes
4! see the UK-AMOR website.
5!
6! This file is part of UKRmol-out (UKRmol+ suite).
7!
8! UKRmol-out is free software: you can redistribute it and/or modify
9! it under the terms of the GNU General Public License as published by
10! the Free Software Foundation, either version 3 of the License, or
11! (at your option) any later version.
12!
13! UKRmol-out is distributed in the hope that it will be useful,
14! but WITHOUT ANY WARRANTY; without even the implied warranty of
15! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16! GNU General Public License for more details.
17!
18! You should have received a copy of the GNU General Public License
19! along with UKRmol-out (in source/COPYING). Alternatively, you can also visit
20! <https://www.gnu.org/licenses/>.
21!
22!> \brief MULTIDIP unit tests
23!> \author J Benda
24!> \date 2020
25!>
26!> This module contains a few sanity unit self-tests that can be executed using "multidip --test".
27!>
28module multidip_tests
29
30 use precisn_gbl, only: wp
31
32 implicit none
33
34contains
35
36 !> \brief Numerical unit tests
37 !> \author J Benda
38 !> \date 2020
39 !>
40 !> Run the low-level computational routines and compare their outputs to reference results.
41 !>
42 subroutine run_tests
43
44 print '(/,A)', 'Running all tests'
45
46 call test_sph
47 call test_permutations
48 call test_coul
49 call test_levin
50 call test_1p_eint
51 call test_1p_cint
52 call test_1p_gint
53 call test_2p_cint
54 call test_2p_gint
55 call test_ang_dist
56
57 end subroutine run_tests
58
59
60 !> \brief Test spherical harmonics
61 !> \author J Benda
62 !> \date 2025
63 !>
64 !> Verifies DIPELM functions for evaluations of real and complex spherical harmonics.
65 !>
66 subroutine test_sph
67
68 use precisn_gbl, only: cfp
69 use special_functions_gbl, only: cfp_resh
70 use dipelm_special_functions, only: a_legendre_p, a_sp_harm, a_re_sp_harm
71 use multidip_params, only: imu, pi
72
73 integer, parameter :: lmax = 3
74 real(wp), parameter :: thetas(5) = [0*pi, pi/4, pi/2, 3*pi/4, pi]
75 real(wp), parameter :: phis(5) = [0*pi, pi/4, pi/2, 3*pi/4, pi]
76
77 real(wp), allocatable :: Plm(:)
78 complex(wp), allocatable :: Ylm(:)
79
80 integer :: itheta, iphi, l, m, lm
81 real(wp) :: Plm_ref((lmax + 1)**2), Xlm, x, fac(0:2*lmax)
82 real(cfp) :: ex, ey, ez, SH(-lmax:lmax, 0:lmax)
83 complex(wp) :: Ylm_ref
84
85 fac(0) = 1
86 do l = 1, ubound(fac, 1)
87 fac(l) = l*fac(l - 1)
88 end do
89
90 print '(/,a,/)', 'Testing normalized associated Legendre polynomials...'
91 print '(2x,a,2x,a,2x,a,4x,a,22x,a)', 'theta', 'l', 'm', 'Plm', 'reference'
92
93 do itheta = 1, size(thetas)
94 call a_legendre_p(lmax, thetas(itheta), plm)
95 x = cos(thetas(itheta))
96 ! Plm from Wikipedia, Condon-Shortley factor included
97 plm_ref = [1._wp, sqrt(1-x*x)/2, x, -sqrt(1-x*x), (1-x*x)/8, x*sqrt(1-x*x)/2, (3*x*x-1)/2, -3*x*sqrt(1-x*x), &
98 3*(1-x*x), (1-x*x)**1.5/48, x*(1-x*x)/8, -(1-5*x*x)*sqrt(1-x*x)/8, (5*x*x-3)*x/2, 3*(1-5*x*x)*sqrt(1-x*x)/2, &
99 15*x*(1-x*x), -15*(1-x*x)**1.5]
100 do l = 0, lmax
101 do m = -l, l
102 lm = l*l + l + m + 1
103 plm_ref(lm) = plm_ref(lm) * sqrt(fac(l - m)/fac(l + m))
104 print '(f7.3,2i3,2e25.15,1x,a)', thetas(itheta), l, m, plm(lm), plm_ref(lm), &
105 merge(' ok', 'WRONG!', abs(plm(lm) - plm_ref(lm)) < 1e-10)
106 end do
107 end do
108 deallocate (plm)
109 end do
110
111 print '(/,a,/)', 'Testing complex spherical harmonics...'
112 print '(2x,a,2x,a,4x,a,2x,a,4x,a,47x,a)', 'theta', 'phi', 'l', 'm', 'Ylm', 'reference'
113
114 do itheta = 1, size(thetas)
115 x = cos(thetas(itheta))
116 ! Plm from Wikipedia, Condon-Shortley factor included
117 plm_ref = [1._wp, sqrt(1-x*x)/2, x, -sqrt(1-x*x), (1-x*x)/8, x*sqrt(1-x*x)/2, (3*x*x-1)/2, -3*x*sqrt(1-x*x), &
118 3*(1-x*x), (1-x*x)**1.5/48, x*(1-x*x)/8, -(1-5*x*x)*sqrt(1-x*x)/8, (5*x*x-3)*x/2, 3*(1-5*x*x)*sqrt(1-x*x)/2, &
119 15*x*(1-x*x), -15*(1-x*x)**1.5]
120 do iphi = 1, size(phis)
121 call a_sp_harm(lmax, thetas(itheta), phis(iphi), ylm)
122 do l = 0, lmax
123 do m = -l, l
124 lm = l*l + l + m + 1
125 ylm_ref = sqrt((2*l+1)/(4*pi)*fac(l - m)/fac(l + m)) * plm_ref(lm) * exp(imu*m*phis(iphi))
126 print '(2f7.3,2i3,4e25.15,1x,a)', thetas(itheta), phis(iphi), l, m, ylm(lm), ylm_ref, &
127 merge(' ok', 'WRONG!', abs(ylm(lm) - ylm_ref) <= 1e-10)
128 end do
129 end do
130 deallocate (ylm)
131 end do
132 end do
133
134 print '(/,a,/)', 'Testing real spherical harmonics...'
135 print '(2x,a,2x,a,4x,a,2x,a,4x,a,22x,a)', 'theta', 'phi', 'l', 'm', 'Xlm', 'reference'
136
137 do itheta = 1, size(thetas)
138 do iphi = 1, size(phis)
139 ex = sin(thetas(itheta))*cos(phis(iphi))
140 ey = sin(thetas(itheta))*sin(phis(iphi))
141 ez = cos(thetas(itheta))
142 call a_re_sp_harm(lmax, thetas(itheta), phis(iphi), ylm)
143 call cfp_resh(sh, ex, ey, ez, lmax)
144 do l = 0, lmax
145 do m = -l, l
146 lm = l*l + l + m + 1
147 xlm = real(ylm(lm))
148 print '(2f7.3,2i3,2e25.15,1x,a)', thetas(itheta), phis(iphi), l, m, xlm, sh(m, l), &
149 merge(' ok', 'WRONG!', abs(xlm - sh(m, l)) <= 1e-10)
150 end do
151 end do
152 deallocate (ylm)
153 end do
154 end do
155
156 end subroutine test_sph
157
158
159 !> \brief Permutations unit tests
160 !> \author J Benda
161 !> \date 2020
162 !>
163 !> Calculate a few permutations using the custom implementation of next_permutation to test
164 !> for correctness.
165 !>
166 subroutine test_permutations
167
168 use multidip_special, only: next_permutation
169
170 integer :: i, order(4)
171
172 order = -1
173
174 print '(/,A,/)', 'Testing permutations of 4 elements...'
175
176 i = 0
177 do while (next_permutation(order))
178 i = i + 1
179 print '(A,I2,A,*(1x,I0))', ' ', i, ':', order
180 end do
181
182 end subroutine test_permutations
183
184
185 !> \brief Negative-energy Coulomb function tests
186 !> \author J Benda
187 !> \date 2020
188 !>
189 !> Compare negative-energy Coulomb (Whittaker) functions calculated by UKRmol-out and by GSL with
190 !> the expected asymptotic behaviour (DLMF §13.19.3).
191 !>
192 subroutine test_coul
193
194 use multidip_params, only: ntermsasy
195 use multidip_special, only: coul_gsl, coul_ukrmol
196
197 real(wp) :: r = 100, z = 1, f, fp, g1, g1p, g2, g2p, g3, g3p, k
198 real(wp) :: Ek(4) = [ -1.000, -0.100, -0.010, -0.001 ]
199 integer :: ls(3) = [ 0, 1, 2 ], ipw, ie, n, l
200
201 print '(/,A)', 'Testing Whittaker functions'
202 print '(/,2x,A,3x,A,7x,A,20x,A,19x,A,17x,A,16x,A,20x,A)', 'l', 'Ek', 'GSL W', 'GSL W''', 'UKRmol W', 'UKRmol W''', &
203 'asy W', 'asy W'''
204
205 do ie = 1, size(ek)
206 k = sqrt(-2*ek(ie))
207 do ipw = 1, size(ls)
208 l = ls(ipw)
209 call coul_gsl(z, l, ek(ie), r, f, fp, g1, g1p)
210 call coul_ukrmol(z, l, ek(ie), r, f, fp, g2, g2p)
211 f = 1; g3 = 0; g3p = 0
212 do n = 0, ntermsasy
213 if (n > 0) f = f * (l + 1 - 1/k + n - 1) * (-l - 1/k + n - 1) / n / (-2*k*r)
214 g3 = g3 + f
215 g3p = g3p + (1/k - n)/r * f
216 end do
217 g3p = exp(-k*r) * (2*k*r)**(1/k) * (-k*g3 + g3p)
218 g3 = exp(-k*r) * (2*k*r)**(1/k) * g3
219 print '(I3,F8.3,6E25.15)', l, ek(ie), g1, g1p, g2, g2p, g3, g3p
220 end do
221 end do
222
223 end subroutine test_coul
224
225
226 !> \brief Tests of the Levin quadrature routine
227 !> \author J Benda
228 !> \date 2021
229 !>
230 subroutine test_levin
231
232 use multidip_romberg, only: nested_cgreen_correct_romberg
233 use multidip_levin, only: nested_cgreen_correct_levin
234
235 real :: t0, t1, t2
236 real(wp) :: ra = 15, rb = 100, c = 0, z = 1
237 integer :: m, s1 = +1, s2, l1, l2
238 complex(wp) :: k1, k2, romb, levn
239
240 print '(/,A)', 'Testing Levin quadrature'
241 print '(/,1x,A,2x,A,6x,A,4x,A,4x,A,4x,A,4x,A,2x,A,2x,A,2x,A,2x,A,4x,A,15x,A,15x,A,17x,A,15x,A,4x,A)', &
242 'Ra', 'Rb', 'Re k1', 'Im k1', 'Re k2', 'Im k2', 'm', 's1', 's2', 'l1', 'l2', &
243 're romberg', 'im romberg', 're levin', 'im levin', 'μs romberg', 'μs levin'
244
245 k1 = (0.005, 0.000)
246 k2 = (0.500, 0.000)
247
248 do l1 = 0, 1
249 do m = 0, 1
250 do s2 = -1, +1, 2
251
252 l2 = l1 + 1
253
254 call cpu_time(t0)
255 romb = 0; call nested_cgreen_correct_romberg(z, ra, rb, c, 1, s1, s2, [ m ], [ l1, l2 ], [ k1, k2 ], romb)
256 call cpu_time(t1)
257 levn = 0; call nested_cgreen_correct_levin (z, ra, rb, c, 1, s1, s2, [ m ], [ l1, l2 ], [ k1, k2 ], levn)
258 call cpu_time(t2)
259
260 print '(2(F4.0,1x),4(F8.3,1x),I4,SP,2I4,SS,2I4,4E25.15,2I12)', &
261 ra, rb, k1, k2, m, s1, s2, l1, l2, romb, levn, int(1e6*(t1 - t0)), int(1e6*(t2 - t1))
262
263 end do
264 end do
265 end do
266
267 k1 = (0.000, 0.005)
268 k2 = (0.500, 0.000)
269
270 do l1 = 0, 1
271 do m = 0, 1
272 do s2 = -1, +1, 2
273
274 l2 = l1 + 1
275
276 call cpu_time(t0)
277 romb = 0; call nested_cgreen_correct_romberg(z, ra, rb, c, 1, s1, s2, [ m ], [ l1, l2 ], [ k1, k2 ], romb)
278 call cpu_time(t1)
279 levn = 0; call nested_cgreen_correct_levin (z, ra, rb, c, 1, s1, s2, [ m ], [ l1, l2 ], [ k1, k2 ], levn)
280 call cpu_time(t2)
281
282 print '(2(F4.0,1x),4(F8.3,1x),I4,SP,2I4,SS,2I4,4E25.15,2I12)', &
283 ra, rb, k1, k2, m, s1, s2, l1, l2, romb, levn, int(1e6*(t1 - t0)), int(1e6*(t2 - t1))
284
285 end do
286 end do
287 end do
288
289 end subroutine test_levin
290
291
292 !> \brief One-photon exponential integral unit tests
293 !> \author J Benda
294 !> \date 2020
295 !>
296 !> Tests the function `nested_exp_integ` for several arguments.
297 !>
298 subroutine test_1p_eint
299
300 use multidip_integ, only: nested_exp_integ
301
302 real(wp) :: Z = 1, a = 150, c = 0
303 real(wp) :: k1(2) = [ 0.011080185297_wp, 0.680747949322_wp ]
304 real(wp) :: k2(2) = [ 1.069105413537_wp, 1.438602233160_wp ]
305
306 integer :: ie, m, s2, ms(1), ss(2)
307
308 complex(wp) :: ks(2), integ
309 complex(wp) :: reference(2,3,2) = reshape([(-4.772929e-01_wp, +8.885647e+01_wp), (-8.424650e-04_wp, +5.923864e-01_wp), &
310 (+9.979744e-06_wp, +3.949186e-03_wp), (+1.307372e+00_wp, +3.240230e+02_wp), &
311 (+3.965361e-02_wp, +2.158463e+00_wp), (+4.697322e-04_wp, +1.437268e-02_wp), &
312 (-2.211202e-01_wp, +7.029786e+01_wp), (-9.899674e-06_wp, +4.686525e-01_wp), &
313 (+9.695004e-06_wp, +3.124290e-03_wp), (-1.752969e+00_wp, +1.992839e+02_wp), &
314 (+8.061387e-05_wp, +1.328557e+00_wp), (+7.894722e-05_wp, +8.855647e-03_wp)], &
315 [2, 3, 2])
316
317 print '(/,A)', 'Testing one-photon exponential integrals'
318 print '(/,2x,A,1x,A,1x,A,2x,A,12x,A,13x,A,2x,A,3x,A,4x,A)', 'm', 's2', 's1', 'k2', 'k1', &
319 're calculated', 'im calculated', 're reference', 'im reference'
320
321 do ie = 1, size(k1)
322 do m = 0, 2
323 do s2 = 1, 2
324
325 ms = [ 1 - m ]
326 ss = [ (-1)**(s2+1), +1 ]
327 ks = [ k1(ie), k2(ie) ]
328
329 integ = nested_exp_integ(z, a, c, 1, ms, ss, ks)
330
331 print '(SP,3I3,SS,2F14.10,SP,4E16.7)', ms(1), ss(2), ss(1), real(ks(2),wp), real(ks(1),wp), &
332 real(integ, wp), aimag(integ), real(reference(s2,m+1,ie), wp), aimag(reference(s2,m+1,ie))
333 end do
334 end do
335 end do
336
337 end subroutine test_1p_eint
338
339
340 !> \brief One-photon Coulomb integral unit tests
341 !> \author J Benda
342 !> \date 2020 - 2023
343 !>
344 !> Tests the function `nested_coul_integ` for several arguments.
345 !>
346 subroutine test_1p_cint
347
348 use multidip_integ, only: nested_coul_integ
349
350 real(wp) :: a = 150, c = 0, z = 1
351 real(wp) :: k1(2) = [ 0.011080185297_wp, 0.680747949322_wp ]
352 real(wp) :: k2(2) = [ 1.069105413537_wp, 1.438602233160_wp ]
353
354 integer :: ie, ipw, s2, ms(1), ss(2), ls(2)
355 integer :: l1(8) = [ 1, 1, 3, 1, 3, 3, 3, 3 ]
356 integer :: l2(8) = [ 0, 2, 2, 2, 2, 4, 2, 4 ]
357
358 complex(wp) :: ks(2), integ
359 complex(wp) :: reference(2,8,2) = reshape([(-3.454030e+13_wp, +1.290212e+14_wp), (+4.518674e+14_wp, +1.770164e+14_wp), &
360 (-1.054015e+14_wp, -8.206518e+13_wp), (-3.390989e+14_wp, +3.473927e+14_wp), &
361 (-1.012103e+14_wp, -8.818320e+13_wp), (-3.209952e+14_wp, +3.673575e+14_wp), &
362 (-1.054015e+14_wp, -8.206518e+13_wp), (-3.390989e+14_wp, +3.473927e+14_wp), &
363 (-1.012103e+14_wp, -8.818320e+13_wp), (-3.209952e+14_wp, +3.673575e+14_wp), &
364 (+1.307948e+14_wp, +3.038993e+13_wp), (+1.111367e+14_wp, -4.753884e+14_wp), &
365 (-1.012103e+14_wp, -8.818320e+13_wp), (-3.209952e+14_wp, +3.673575e+14_wp), &
366 (+1.307948e+14_wp, +3.038993e+13_wp), (+1.111367e+14_wp, -4.753884e+14_wp), &
367 (-6.382646e+01_wp, +2.799838e+01_wp), (+7.071896e+01_wp, +1.844497e+02_wp), &
368 (+1.584632e+01_wp, -6.787748e+01_wp), (-1.900621e+02_wp, -5.395706e+01_wp), &
369 (+5.054597e+01_wp, +4.803432e+01_wp), (+4.962959e+01_wp, +1.911980e+02_wp), &
370 (+1.584632e+01_wp, -6.787748e+01_wp), (-1.900621e+02_wp, -5.395706e+01_wp), &
371 (+5.054597e+01_wp, +4.803432e+01_wp), (+4.962959e+01_wp, +1.911980e+02_wp), &
372 (-6.443654e+01_wp, -2.668020e+01_wp), (-1.150300e+02_wp, -1.606738e+02_wp), &
373 (+5.054597e+01_wp, +4.803432e+01_wp), (+4.962959e+01_wp, +1.911980e+02_wp), &
374 (-6.443654e+01_wp, -2.668020e+01_wp), (-1.150300e+02_wp, -1.606738e+02_wp)], &
375 [2, 8, 2])
376
377 print '(/,A)', 'Testing one-photon Coulomb integrals'
378 print '(/,1x,A,1x,A,1x,A,2x,A,12x,A,13x,A,2x,A,3x,A,4x,A)', 's2', 'l2', 'l1', 'k2', 'k1', &
379 're calculated', 'im calculated', 're reference', 'im reference'
380
381 do ie = 1, size(k1)
382 do ipw = 1, size(l1)
383 do s2 = 1, 2
384
385 ms = [ 1 ]
386 ss = [ (-1)**(s2+1), +1 ]
387 ls = [ l1(ipw), l2(ipw) ]
388 ks = [ k1(ie), k2(ie) ]
389
390 integ = nested_coul_integ(z, a, c, 1, ms, ss, ls, ks)
391
392 print '(SP,I3,SS,2I3,2F14.10,SP,4E16.7)', (-1)**(s2+1), l2(ipw), l1(ipw), k2(ie), k1(ie), &
393 real(integ, wp), aimag(integ), real(reference(s2,ipw,ie), wp), aimag(reference(s2,ipw,ie))
394 end do
395 end do
396 end do
397
398 end subroutine test_1p_cint
399
400
401 !> \brief One-photon Green integral unit tests
402 !> \author J Benda
403 !> \date 2021 - 2023
404 !>
405 !> Tests the function `nested_cgreen_integ` for several arguments. Compares results obtained from asymptotic integration
406 !> over interval (+50,+inf) to results obtained using a combination of numerical integration of exact Coulomb functions
407 !> at interval (+50,+150) and asymptotic formulas on (+150,+inf).
408 !>
409 subroutine test_1p_gint
410
411 use multidip_integ, only: nested_cgreen_integ
412
413 real(wp) :: a = 50, b = 150, c = 0, z = 1
414
415 real(wp) :: k1(2) = [ 0.011080185297_wp, 0.680747949322_wp ]
416 real(wp) :: k2(2) = [ 1.069105413537_wp, 1.438602233160_wp ]
417
418 integer :: l1(8) = [ 1, 1, 3, 1, 3, 3, 3, 3 ]
419 integer :: l2(8) = [ 0, 2, 2, 2, 2, 4, 2, 4 ]
420
421 integer :: ie, ipw, s1, s2, sa, sb, ms(1), ls(2)
422
423 complex(wp) :: ks(2), integ_ai, integ_bi, integ_ab
424
425 print '(/,A)', 'Testing one-photon Green integrals (numerical integration)'
426 print '(/,4(1x,A),2x,A,12x,A,12x,A,3x,A,3x,A,3x,A)', 's2', 's1', 'l2', 'l1', 'k2', 'k1', &
427 're numerical ', 'im numerical ', 're asymptotic', 'im asymptotic'
428
429 do ie = 1, size(k1)
430 do ipw = 1, size(l1)
431 do s2 = 1, 2
432 do s1 = 1, 2
433
434 ms = [ 1 ]
435 ls = [ l1(ipw), l2(ipw) ]
436 ks = [ k1(ie), k2(ie) ]
437
438 sa = (-1)**(s1 + 1)
439 sb = (-1)**(s2 + 1)
440
441 integ_ai = nested_cgreen_integ(z, a, a, c, 1, sa, sb, ms, ls, ks)
442 integ_bi = nested_cgreen_integ(z, b, b, c, 1, sa, sb, ms, ls, ks)
443 integ_ab = nested_cgreen_integ(z, a, b, c, 1, sa, sb, ms, ls, ks)
444
445 print '(SP,2I3,SS,2I3,2F14.10,SP,8E16.7)', sb, sa, l2(ipw), l1(ipw), k2(ie), k1(ie), &
446 real(integ_ab - integ_bi, wp), aimag(integ_ab - integ_bi), &
447 real(integ_ai - integ_bi, wp), aimag(integ_ai - integ_bi)
448
449 end do
450 end do
451 end do
452 end do
453
454 end subroutine test_1p_gint
455
456
457 !> \brief Two-photon Coulomb integral unit tests
458 !> \author J Benda
459 !> \date 2020 - 2023
460 !>
461 !> Tests the function `nested_coul_integ` for several arguments.
462 !>
463 subroutine test_2p_cint
464
465 use multidip_integ, only: nested_coul_integ
466
467 real(wp) :: a = 150, c = 0, z = 1
468
469 real(wp) :: k1(2) = [ 0.011080185297_wp, 0.680747949322_wp ]
470 real(wp) :: k2(2) = [ 1.069105413537_wp, 1.438602233160_wp ]
471 real(wp) :: k3(2) = [ 1.511902774652_wp, 1.917220644579_wp ]
472
473 integer :: l1(8) = [ 1, 1, 3, 1, 3, 3, 3, 3 ]
474 integer :: l2(8) = [ 0, 2, 2, 2, 2, 4, 2, 4 ]
475 integer :: l3(8) = [ 1, 1, 1, 3, 3, 3, 3, 3 ]
476
477 integer :: ie, ipw, s4, ms(2), ss(4), ls(4)
478
479 complex(wp) :: ks(4), integ1, integ2, integ
480 complex(wp) :: reference(2,8,2) = reshape([(-1.465883e+16_wp, -4.950701e+15_wp), (-3.842765e+16_wp, 8.620168e+16_wp), &
481 (-1.466146e+16_wp, -4.951551e+15_wp), (-3.843793e+16_wp, 8.623263e+16_wp), &
482 (-1.443664e+16_wp, -5.780723e+15_wp), (-3.378946e+16_wp, 8.865409e+16_wp), &
483 ( 1.520251e+16_wp, -2.902864e+15_wp), (-8.964058e+15_wp, -9.400575e+16_wp), &
484 ( 1.541476e+16_wp, -2.070261e+15_wp), (-1.420349e+16_wp, -9.382680e+16_wp), &
485 ( 1.542117e+16_wp, -2.071208e+15_wp), (-1.422236e+16_wp, -9.390101e+16_wp), &
486 ( 1.541476e+16_wp, -2.070261e+15_wp), (-1.420349e+16_wp, -9.382680e+16_wp), &
487 ( 1.542117e+16_wp, -2.071208e+15_wp), (-1.422236e+16_wp, -9.390101e+16_wp), &
488 (-4.322904e+03_wp, -4.890496e+03_wp), (-2.301503e+04_wp, 1.829823e+04_wp), &
489 (-4.323345e+03_wp, -4.890988e+03_wp), (-2.301886e+04_wp, 1.830140e+04_wp), &
490 ( 6.409175e+03_wp, -1.254518e+03_wp), ( 2.742317e+04_wp, 1.058896e+04_wp), &
491 ( 5.914821e+03_wp, 2.763232e+03_wp), ( 1.382079e+04_wp, -2.596191e+04_wp), &
492 (-5.378829e+03_wp, 3.705000e+03_wp), (-2.937519e+04_wp, 1.216512e+03_wp), &
493 (-5.380092e+03_wp, 3.705885e+03_wp), (-2.938678e+04_wp, 1.217206e+03_wp), &
494 (-5.378829e+03_wp, 3.705000e+03_wp), (-2.937519e+04_wp, 1.216512e+03_wp), &
495 (-5.380092e+03_wp, 3.705885e+03_wp), (-2.938678e+04_wp, 1.217206e+03_wp)], &
496 [2, 8, 2])
497
498 print '(/,A)', 'Testing two-photon Coulomb integrals'
499 print '(/,4(1x,A),2x,A,12x,A,12x,A,12x,A,3x,A,3x,A,4x,A)', 's4', 'l3', 'l2', 'l1', 'k3', 'k2', 'k1', &
500 're calculated', 'im calculated', 're reference', 'im reference'
501
502 do ie = 1, size(k1)
503 do ipw = 1, size(l1)
504 do s4 = 1, 2
505
506 ms = [ 1, 1 ]
507
508 ss = [ (-1)**(s4+1), +1, -1, +1 ]
509 ls = [ l1(ipw), l2(ipw), l2(ipw), l3(ipw) ]
510 ks = [ k1(ie), k2(ie), k2(ie), k3(ie) ]
511
512 integ1 = nested_coul_integ(z, a, c, 2, ms, ss, ls, ks)
513
514 ss = [ +1, +1, -1, (-1)**(s4+1) ]
515 ls = [ l3(ipw), l2(ipw), l2(ipw), l1(ipw) ]
516 ks = [ k3(ie), k2(ie), k2(ie), k1(ie) ]
517
518 integ2 = nested_coul_integ(z, a, c, 2, ms, ss, ls, ks)
519
520 integ = integ1 + integ2
521
522 print '(SP,I3,SS,3I3,3F14.10,SP,4E16.7)', (-1)**(s4+1), l3(ie), l2(ipw), l1(ipw), k3(ie), k2(ie), k1(ie), &
523 real(integ, wp), aimag(integ), real(reference(s4,ipw,ie), wp), aimag(reference(s4,ipw,ie))
524 end do
525 end do
526 end do
527
528 end subroutine test_2p_cint
529
530
531 !> \brief Two-photon Green integral unit tests
532 !> \author J Benda
533 !> \date 2020 - 2023
534 !>
535 !> Tests the function `nested_cgreen_integ` for several arguments.
536 !>
537 subroutine test_2p_gint
538
539 use multidip_integ, only: nested_cgreen_integ
540
541 real(wp) :: a = 150, c = 0, z = 1
542
543 real(wp) :: k1(2) = [ 0.011080185297_wp, 0.680747949322_wp ]
544 real(wp) :: k2(2) = [ 1.069105413537_wp, 1.438602233160_wp ]
545 real(wp) :: k3(2) = [ 1.511902774652_wp, 1.917220644579_wp ]
546
547 integer :: l1(8) = [ 1, 1, 3, 1, 3, 3, 3, 3 ]
548 integer :: l2(8) = [ 0, 2, 2, 2, 2, 4, 2, 4 ]
549 integer :: l3(8) = [ 1, 1, 1, 3, 3, 3, 3, 3 ]
550
551 integer :: ie, ipw, s4, ms(2), ls(3)
552
553 complex(wp) :: ks(3), integ1, integ2
554 complex(wp) :: reference(2,2,8,2) = reshape([(-8.310927e+15_wp, 7.523384e+15_wp), (-5.909915e+16_wp, 6.931934e+16_wp), &
555 ( 5.679523e+16_wp, 4.672527e+16_wp), ( 8.225897e+24_wp, 1.848188e+24_wp), &
556 (-6.506509e+15_wp, 2.066669e+16_wp), (-1.618385e+16_wp, 4.224605e+15_wp), &
557 ( 1.050096e+17_wp, 4.554189e+16_wp), (-4.312479e+24_wp, -7.240598e+24_wp), &
558 (-7.673705e+15_wp, 2.037601e+16_wp), (-1.647265e+16_wp, 3.342092e+15_wp), &
559 ( 1.078882e+17_wp, 3.987695e+16_wp), (-4.728274e+24_wp, -7.026109e+24_wp), &
560 (-4.503787e+15_wp, -2.119630e+16_wp), ( 1.619053e+16_wp, 4.287343e+15_wp), &
561 (-1.138495e+17_wp, 1.201762e+16_wp), ( 2.556995e+23_wp, 8.432496e+24_wp), &
562 (-3.344463e+15_wp, -2.151751e+16_wp), ( 1.600768e+16_wp, 5.199000e+15_wp), &
563 (-1.135684e+17_wp, 1.836697e+16_wp), ( 7.238276e+23_wp, 8.446832e+24_wp), &
564 (-8.612323e+15_wp, -1.723199e+16_wp), ( 3.535435e+16_wp, 4.008848e+16_wp), &
565 (-9.807095e+16_wp, 3.755844e+16_wp), ( 7.275593e+24_wp, 4.341652e+24_wp), &
566 (-3.344463e+15_wp, -2.151751e+16_wp), ( 1.600768e+16_wp, 5.199000e+15_wp), &
567 (-1.135684e+17_wp, 1.836697e+16_wp), ( 7.238276e+23_wp, 8.446832e+24_wp), &
568 (-8.612323e+15_wp, -1.723199e+16_wp), ( 3.535435e+16_wp, 4.008848e+16_wp), &
569 (-9.807095e+16_wp, 3.755844e+16_wp), ( 7.275593e+24_wp, 4.341652e+24_wp), &
570 (-1.247516e+03_wp, 3.118370e+03_wp), (-2.763140e+04_wp, 1.784564e+04_wp), &
571 ( 1.332880e+04_wp, 9.920923e+03_wp), ( 6.345060e+03_wp, 4.578023e+04_wp), &
572 (-3.896753e+03_wp, 9.080720e+02_wp), (-6.662204e+03_wp, 3.014237e+04_wp), &
573 ( 6.717303e+03_wp, 1.712697e+04_wp), ( 4.397613e+04_wp, -1.193480e+04_wp), &
574 ( 1.189054e+03_wp, -3.822275e+03_wp), (-2.263210e+04_wp, -2.099869e+04_wp), &
575 ( 1.136951e+04_wp, -1.445428e+04_wp), (-3.253996e+04_wp, -3.189015e+04_wp), &
576 ( 3.212357e+03_wp, -2.385969e+03_wp), ( 1.813020e+04_wp, -2.499891e+04_wp), &
577 (-1.298806e+04_wp, -1.303291e+04_wp), (-4.510261e+04_wp, -6.582287e+03_wp), &
578 ( 4.327121e+02_wp, 3.979850e+03_wp), ( 1.239242e+04_wp, 2.828936e+04_wp), &
579 (-4.667920e+03_wp, 1.779027e+04_wp), ( 1.713954e+04_wp, 4.222912e+04_wp), &
580 ( 1.147174e+03_wp, 5.355515e+03_wp), ( 1.429665e+04_wp, 1.753491e+04_wp), &
581 (-5.010378e+03_wp, 2.217364e+04_wp), ( 3.781810e+04_wp, 1.923663e+04_wp), &
582 ( 4.327121e+02_wp, 3.979850e+03_wp), ( 1.239242e+04_wp, 2.828936e+04_wp), &
583 (-4.667920e+03_wp, 1.779027e+04_wp), ( 1.713954e+04_wp, 4.222912e+04_wp), &
584 ( 1.147174e+03_wp, 5.355515e+03_wp), ( 1.429665e+04_wp, 1.753491e+04_wp), &
585 (-5.010378e+03_wp, 2.217364e+04_wp), ( 3.781810e+04_wp, 1.923663e+04_wp)], &
586 [2,2,8,2])
587
588 print '(/,A)', 'Testing two-photon Green integrals'
589 print '(/,5(1x,A),2x,A,12x,A,12x,A,12x,A,3x,A,3x,A,4x,A)', 's4', 's1', 'l3', 'l2', 'l1', 'k3', 'k2', 'k1', &
590 're calculated', 'im calculated', 're reference', 'im reference'
591
592 do ie = 1, size(k1)
593 do ipw = 1, size(l1)
594 do s4 = 1, 2
595
596 ms = [ 1, 1 ]
597 ls = [ l1(ipw), l2(ipw), l3(ipw) ]
598 ks = [ k1(ie), k2(ie), k3(ie) ]
599
600 integ1 = nested_cgreen_integ(z, a, a, c, 2, (-1)**(s4+1), +1, ms, ls, ks)
601 integ2 = nested_cgreen_integ(z, a, a, c, 2, (-1)**(s4+1), -1, ms, ls, ks)
602
603 print '(SP,2I3,SS,3I3,3F14.10,SP,8E16.7)', (-1)**(s4+1), +1, l3(ipw), l2(ipw), l1(ipw), k3(ie), k2(ie), k1(ie),&
604 real(integ1, wp), aimag(integ1), real(reference(1,s4,ipw,ie), wp), aimag(reference(1,s4,ipw,ie))
605 print '(SP,2I3,SS,3I3,3F14.10,SP,8E16.7)', (-1)**(s4+1), -1, l3(ipw), l2(ipw), l1(ipw), k3(ie), k2(ie), k1(ie),&
606 real(integ2, wp), aimag(integ2), real(reference(2,s4,ipw,ie), wp), aimag(reference(2,s4,ipw,ie))
607 end do
608 end do
609 end do
610 end subroutine test_2p_gint
611
612
613 !> \brief Angular algebra unit tests
614 !> \author J Benda
615 !> \date 2020
616 !>
617 !> Tests the function `beta_contraction_tensor` for several arguments.
618 !>
619 subroutine test_ang_dist
620
621 use dipelm_special_functions, only: threej
622 use multidip_params, only: rone
623 use multidip_special, only: beta_contraction_tensor, beta_2p_demekhin, beta_2p_arb_pol_sph_harm
624
625 integer :: J = 0, li = 1, lj = 1, mi = 0, mj = 0, ki = 0, kj = 0, qi = 0, qj = 0, p = 0
626 real(wp) :: beta, ref, beta_ap
627
628 print '(/,A)', 'Testing one-photon angular integrals'
629 print '(/,8A4,2A12)', 'J', 'p', 'li', 'lj', 'mi', 'mj', 'qi', 'qj', 'TJij', 'reference'
630
631 do j = 0, 2
632 do p = 0, +1
633 do mi = -1, +1
634 do mj = -1, +1
635 do qi = -1, +1
636 do qj = -1, +1
637 beta = beta_contraction_tensor(j, 1, [ p ], li, mi, [ qi ], lj, mj, [ qj ])
638 ref = (-1)**(p + mi + qi) * (2*j + 1) * sqrt((2*li + rone)*(2*lj + rone)) &
639 * threej(2*li, 0, 2*lj, 0, 2*j, 0) &
640 * threej(2*1, 2*p, 2*1, -2*p, 2*j, 0) &
641 * threej(2*li, 2*mi, 2*lj, -2*mj, 2*j, -2*(mi-mj)) &
642 * threej(2*1, 2*qi, 2*1, -2*qj, 2*j, -2*(mi-mj))
643 if (abs(beta) > 1e-10 .or. abs(ref) > 1e-10) then
644 print '(I4,SP,I4,SS,2I4,SP,4I4,2F12.7,A8)', &
645 j, p, li, lj, mi, mj, qi, qj, beta, ref, &
646 merge(' ok', 'WRONG!', abs(beta - ref) / (abs(beta) + abs(ref)) < 1e-5)
647 end if
648 end do
649 end do
650 end do
651 end do
652 end do
653 end do
654
655 print '(/,A)', 'Testing two-photon angular integrals'
656 print '(/,10A4,3A12)', 'J', 'p', 'li', 'lj', 'mi', 'mj', 'ki', 'kj', 'qi', 'qj', 'TJij', 'TJij(ap)', 'reference'
657
658 do j = 0, 2
659 do p = 0, +1
660 do mi = -1, +1
661 do mj = -1, +1
662 do ki = -1, +1
663 do kj = -1, +1
664 do qi = -1, +1
665 do qj = -1, +1
666 beta = beta_contraction_tensor(j, 2, [ p, p ], li, mi, [ ki, qi ], lj, mj, [ kj, qj ])
667 beta_ap = beta_2p_arb_pol_sph_harm(j, 0, 2, [p,p], li, mi, [ki,qi], [p,p], lj, mj, [kj,qj])
668 ref = beta_2p_demekhin(j, p, p, li, mi, ki, qi, lj, mj, kj, qj)
669 if (abs(beta) > 1e-10 .or. abs(beta_ap) > 1e-10 .or. abs(ref) > 1e-10) then
670 print '(I4,SP,I4,SS,2I4,SP,6I4,3F12.7,2A8)', &
671 j, p, li, lj, mi, mj, ki, kj, qi, qj, beta, beta_ap, ref, &
672 merge(' ok', 'WRONG!', abs(beta - ref) / (abs(beta) + abs(ref)) < 1e-5), &
673 merge(' ok', 'WRONG!', abs(beta_ap - ref) / (abs(beta_ap) + abs(ref)) < 1e-5)
674 end if
675 end do
676 end do
677 end do
678 end do
679 end do
680 end do
681 end do
682 end do
683
684 end subroutine test_ang_dist
685
686end module multidip_tests
multidip_integ
Special integrals needed by MULTIDIP.
Definition multidip_integ.f90:29
multidip_integ::nested_cgreen_integ
complex(wp) function nested_cgreen_integ(z, a, r0, c, n, sa, sb, m, l, k)
Multi-dimensional integral of Coulomb-Hankel and Coulomb-Green functions and powers.
Definition multidip_integ.f90:356
multidip_integ::nested_coul_integ
complex(wp) function nested_coul_integ(z, a, c, n, m, s, l, k)
Multi-dimensional triangular integral of Coulomb-Hankel functions and powers.
Definition multidip_integ.f90:228
multidip_integ::nested_exp_integ
complex(wp) function nested_exp_integ(z, a, c, n, m, s, k)
Multi-dimensional triangular integral of exponentials and powers.
Definition multidip_integ.f90:69
multidip_levin
Levin quadrature for numerical integration.
Definition multidip_levin.f90:37
multidip_levin::nested_cgreen_correct_levin
subroutine nested_cgreen_correct_levin(z, a, b, c, n, sa, sb, m, l, k, v)
Numerically correct asymptotic approximation of the Coulomb-Green integral (Levin integration)
Definition multidip_levin.f90:70
multidip_params
Hard-coded parameters of MULTIDIP.
Definition multidip_params.F90:29
multidip_params::pi
real(wp), parameter pi
Definition multidip_params.F90:60
multidip_params::imu
complex(wp), parameter imu
Definition multidip_params.F90:64
multidip_params::rone
real(wp), parameter rone
Definition multidip_params.F90:58
multidip_params::ntermsasy
integer, parameter ntermsasy
Number of terms in expansion of the Coulomb-Hankel functions.
Definition multidip_params.F90:35
multidip_romberg
Romberg quadrature for numerical integration.
Definition multidip_romberg.f90:26
multidip_romberg::nested_cgreen_correct_romberg
subroutine nested_cgreen_correct_romberg(z, a, b, c, n, sa, sb, m, l, k, v)
Numerically correct asymptotic approximation of the Coulomb-Green integral (Romberg integration)
Definition multidip_romberg.f90:59
multidip_special
Special functions and objects used by MULTIDIP.
Definition multidip_special.F90:26
multidip_special::beta_contraction_tensor
real(wp) function beta_contraction_tensor(j, n, p, li, mi, qi, lj, mj, qj)
Angular part of the beta parameters.
Definition multidip_special.F90:994
multidip_special::coul_ukrmol
subroutine coul_ukrmol(z, l, ek, r, f, fp, g, gp)
Coulomb functions (UKRmol)
Definition multidip_special.F90:664
multidip_special::coul_gsl
subroutine coul_gsl(z, l, ek, r, f, fp, g, gp)
Coulomb functions (GSL)
Definition multidip_special.F90:612
multidip_special::beta_2p_arb_pol_sph_harm
real(wp) function beta_2p_arb_pol_sph_harm(lbig, mbig, n, pb, lp, mp, qb, pa, l, m, qa)
Two-photon asymmetry parameter for the arbitrary polarized case.
Definition multidip_special.F90:1126
multidip_special::beta_2p_demekhin
real(wp) function beta_2p_demekhin(l, p1, p2, li, mi, q1i, q2i, lj, mj, q1j, q2j)
Two-photon asymmetry parameter.
Definition multidip_special.F90:1194
multidip_special::next_permutation
logical function next_permutation(p)
Construct or advance permutation.
Definition multidip_special.F90:847
multidip_tests
MULTIDIP unit tests.
Definition multidip_tests.f90:28
multidip_tests::test_2p_gint
subroutine test_2p_gint
Two-photon Green integral unit tests.
Definition multidip_tests.f90:538
multidip_tests::test_levin
subroutine test_levin
Tests of the Levin quadrature routine.
Definition multidip_tests.f90:231
multidip_tests::test_coul
subroutine test_coul
Negative-energy Coulomb function tests.
Definition multidip_tests.f90:193
multidip_tests::test_2p_cint
subroutine test_2p_cint
Two-photon Coulomb integral unit tests.
Definition multidip_tests.f90:464
multidip_tests::test_permutations
subroutine test_permutations
Permutations unit tests.
Definition multidip_tests.f90:167
multidip_tests::test_1p_eint
subroutine test_1p_eint
One-photon exponential integral unit tests.
Definition multidip_tests.f90:299
multidip_tests::run_tests
subroutine run_tests
Numerical unit tests.
Definition multidip_tests.f90:43
multidip_tests::test_1p_gint
subroutine test_1p_gint
One-photon Green integral unit tests.
Definition multidip_tests.f90:410
multidip_tests::test_sph
subroutine test_sph
Test spherical harmonics.
Definition multidip_tests.f90:67
multidip_tests::test_1p_cint
subroutine test_1p_cint
One-photon Coulomb integral unit tests.
Definition multidip_tests.f90:347
multidip_tests::test_ang_dist
subroutine test_ang_dist
Angular algebra unit tests.
Definition multidip_tests.f90:620