Multidip  1.0
Multi-photon matrix elements
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 !
28 module multidip_tests
29 
30  use precisn_gbl, only: wp
31 
32  implicit none
33 
34 contains
35 
42  subroutine run_tests
43 
44  print '(/,A)', 'Running all tests'
45 
46  call test_permutations
47  call test_coul
48  call test_levin
49  call test_1p_eint
50  call test_1p_cint
51  call test_1p_gint
52  call test_2p_cint
53  call test_2p_gint
54  call test_ang_dist
55 
56  end subroutine run_tests
57 
58 
66  subroutine test_permutations
67 
68  use multidip_special, only: next_permutation
69 
70  integer :: i, order(4)
71 
72  order = -1
73 
74  print '(/,A,/)', 'Testing permutations of 4 elements...'
75 
76  i = 0
77  do while (next_permutation(order))
78  i = i + 1
79  print '(A,I2,A,*(1x,I0))', ' ', i, ':', order
80  end do
81 
82  end subroutine test_permutations
83 
84 
92  subroutine test_coul
93 
94  use multidip_params, only: ntermsasy
95  use multidip_special, only: coul_gsl, coul_ukrmol
96 
97  real(wp) :: r = 100, z = 1, f, fp, g1, g1p, g2, g2p, g3, g3p, k
98  real(wp) :: Ek(4) = [ -1.000, -0.100, -0.010, -0.001 ]
99  integer :: ls(3) = [ 0, 1, 2 ], ipw, ie, n, l
100 
101  print '(/,A)', 'Testing Whittaker functions'
102  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''', &
103  'asy W', 'asy W'''
104 
105  do ie = 1, size(ek)
106  k = sqrt(-2*ek(ie))
107  do ipw = 1, size(ls)
108  l = ls(ipw)
109  call coul_gsl(z, l, ek(ie), r, f, fp, g1, g1p)
110  call coul_ukrmol(z, l, ek(ie), r, f, fp, g2, g2p)
111  f = 1; g3 = 0; g3p = 0
112  do n = 0, ntermsasy
113  if (n > 0) f = f * (l + 1 - 1/k + n - 1) * (-l - 1/k + n - 1) / n / (-2*k*r)
114  g3 = g3 + f
115  g3p = g3p + (1/k - n)/r * f
116  end do
117  g3p = exp(-k*r) * (2*k*r)**(1/k) * (-k*g3 + g3p)
118  g3 = exp(-k*r) * (2*k*r)**(1/k) * g3
119  print '(I3,F8.3,6E25.15)', l, ek(ie), g1, g1p, g2, g2p, g3, g3p
120  end do
121  end do
122 
123  end subroutine test_coul
124 
125 
130  subroutine test_levin
131 
132  use multidip_romberg, only: nested_cgreen_correct_romberg
133  use multidip_levin, only: nested_cgreen_correct_levin
134 
135  real :: t0, t1, t2
136  real(wp) :: ra = 15, rb = 100, c = 0, z = 1
137  integer :: m, s1 = +1, s2, l1, l2
138  complex(wp) :: k1, k2, romb, levn
139 
140  print '(/,A)', 'Testing Levin quadrature'
141  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)', &
142  'Ra', 'Rb', 'Re k1', 'Im k1', 'Re k2', 'Im k2', 'm', 's1', 's2', 'l1', 'l2', &
143  're romberg', 'im romberg', 're levin', 'im levin', μ's romberg', μ's levin'
144 
145  k1 = (0.005, 0.000)
146  k2 = (0.500, 0.000)
147 
148  do l1 = 0, 1
149  do m = 0, 1
150  do s2 = -1, +1, 2
151 
152  l2 = l1 + 1
153 
154  call cpu_time(t0)
155  romb = 0; call nested_cgreen_correct_romberg(z, ra, rb, c, 1, s1, s2, [ m ], [ l1, l2 ], [ k1, k2 ], romb)
156  call cpu_time(t1)
157  levn = 0; call nested_cgreen_correct_levin (z, ra, rb, c, 1, s1, s2, [ m ], [ l1, l2 ], [ k1, k2 ], levn)
158  call cpu_time(t2)
159 
160  print '(2(F4.0,1x),4(F8.3,1x),I4,SP,2I4,SS,2I4,4E25.15,2I12)', &
161  ra, rb, k1, k2, m, s1, s2, l1, l2, romb, levn, int(1e6*(t1 - t0)), int(1e6*(t2 - t1))
162 
163  end do
164  end do
165  end do
166 
167  k1 = (0.000, 0.005)
168  k2 = (0.500, 0.000)
169 
170  do l1 = 0, 1
171  do m = 0, 1
172  do s2 = -1, +1, 2
173 
174  l2 = l1 + 1
175 
176  call cpu_time(t0)
177  romb = 0; call nested_cgreen_correct_romberg(z, ra, rb, c, 1, s1, s2, [ m ], [ l1, l2 ], [ k1, k2 ], romb)
178  call cpu_time(t1)
179  levn = 0; call nested_cgreen_correct_levin (z, ra, rb, c, 1, s1, s2, [ m ], [ l1, l2 ], [ k1, k2 ], levn)
180  call cpu_time(t2)
181 
182  print '(2(F4.0,1x),4(F8.3,1x),I4,SP,2I4,SS,2I4,4E25.15,2I12)', &
183  ra, rb, k1, k2, m, s1, s2, l1, l2, romb, levn, int(1e6*(t1 - t0)), int(1e6*(t2 - t1))
184 
185  end do
186  end do
187  end do
188 
189  end subroutine test_levin
190 
191 
198  subroutine test_1p_eint
199 
200  use multidip_integ, only: nested_exp_integ
201 
202  real(wp) :: Z = 1, a = 150, c = 0
203  real(wp) :: k1(2) = [ 0.011080185297_wp, 0.680747949322_wp ]
204  real(wp) :: k2(2) = [ 1.069105413537_wp, 1.438602233160_wp ]
205 
206  integer :: ie, m, s2, ms(1), ss(2)
207 
208  complex(wp) :: ks(2), integ
209  complex(wp) :: reference(2,3,2) = reshape([(-4.772929e-01_wp, +8.885647e+01_wp), (-8.424650e-04_wp, +5.923864e-01_wp), &
210  (+9.979744e-06_wp, +3.949186e-03_wp), (+1.307372e+00_wp, +3.240230e+02_wp), &
211  (+3.965361e-02_wp, +2.158463e+00_wp), (+4.697322e-04_wp, +1.437268e-02_wp), &
212  (-2.211202e-01_wp, +7.029786e+01_wp), (-9.899674e-06_wp, +4.686525e-01_wp), &
213  (+9.695004e-06_wp, +3.124290e-03_wp), (-1.752969e+00_wp, +1.992839e+02_wp), &
214  (+8.061387e-05_wp, +1.328557e+00_wp), (+7.894722e-05_wp, +8.855647e-03_wp)], &
215  [2, 3, 2])
216 
217  print '(/,A)', 'Testing one-photon exponential integrals'
218  print '(/,2x,A,1x,A,1x,A,2x,A,12x,A,13x,A,2x,A,3x,A,4x,A)', 'm', 's2', 's1', 'k2', 'k1', &
219  're calculated', 'im calculated', 're reference', 'im reference'
220 
221  do ie = 1, size(k1)
222  do m = 0, 2
223  do s2 = 1, 2
224 
225  ms = [ 1 - m ]
226  ss = [ (-1)**(s2+1), +1 ]
227  ks = [ k1(ie), k2(ie) ]
228 
229  integ = nested_exp_integ(z, a, c, 1, ms, ss, ks)
230 
231  print '(SP,3I3,SS,2F14.10,SP,4E16.7)', ms(1), ss(2), ss(1), real(ks(2),wp), real(ks(1),wp), &
232  real(integ, wp), aimag(integ), real(reference(s2,m+1,ie), wp), aimag(reference(s2,m+1,ie))
233  end do
234  end do
235  end do
236 
237  end subroutine test_1p_eint
238 
239 
246  subroutine test_1p_cint
247 
248  use multidip_integ, only: nested_coul_integ
249 
250  real(wp) :: a = 150, c = 0, z = 1
251  real(wp) :: k1(2) = [ 0.011080185297_wp, 0.680747949322_wp ]
252  real(wp) :: k2(2) = [ 1.069105413537_wp, 1.438602233160_wp ]
253 
254  integer :: ie, ipw, s2, ms(1), ss(2), ls(2)
255  integer :: l1(8) = [ 1, 1, 3, 1, 3, 3, 3, 3 ]
256  integer :: l2(8) = [ 0, 2, 2, 2, 2, 4, 2, 4 ]
257 
258  complex(wp) :: ks(2), integ
259  complex(wp) :: reference(2,8,2) = reshape([(-3.454030e+13_wp, +1.290212e+14_wp), (+4.518674e+14_wp, +1.770164e+14_wp), &
260  (-1.054015e+14_wp, -8.206518e+13_wp), (-3.390989e+14_wp, +3.473927e+14_wp), &
261  (-1.012103e+14_wp, -8.818320e+13_wp), (-3.209952e+14_wp, +3.673575e+14_wp), &
262  (-1.054015e+14_wp, -8.206518e+13_wp), (-3.390989e+14_wp, +3.473927e+14_wp), &
263  (-1.012103e+14_wp, -8.818320e+13_wp), (-3.209952e+14_wp, +3.673575e+14_wp), &
264  (+1.307948e+14_wp, +3.038993e+13_wp), (+1.111367e+14_wp, -4.753884e+14_wp), &
265  (-1.012103e+14_wp, -8.818320e+13_wp), (-3.209952e+14_wp, +3.673575e+14_wp), &
266  (+1.307948e+14_wp, +3.038993e+13_wp), (+1.111367e+14_wp, -4.753884e+14_wp), &
267  (-6.382646e+01_wp, +2.799838e+01_wp), (+7.071896e+01_wp, +1.844497e+02_wp), &
268  (+1.584632e+01_wp, -6.787748e+01_wp), (-1.900621e+02_wp, -5.395706e+01_wp), &
269  (+5.054597e+01_wp, +4.803432e+01_wp), (+4.962959e+01_wp, +1.911980e+02_wp), &
270  (+1.584632e+01_wp, -6.787748e+01_wp), (-1.900621e+02_wp, -5.395706e+01_wp), &
271  (+5.054597e+01_wp, +4.803432e+01_wp), (+4.962959e+01_wp, +1.911980e+02_wp), &
272  (-6.443654e+01_wp, -2.668020e+01_wp), (-1.150300e+02_wp, -1.606738e+02_wp), &
273  (+5.054597e+01_wp, +4.803432e+01_wp), (+4.962959e+01_wp, +1.911980e+02_wp), &
274  (-6.443654e+01_wp, -2.668020e+01_wp), (-1.150300e+02_wp, -1.606738e+02_wp)], &
275  [2, 8, 2])
276 
277  print '(/,A)', 'Testing one-photon Coulomb integrals'
278  print '(/,1x,A,1x,A,1x,A,2x,A,12x,A,13x,A,2x,A,3x,A,4x,A)', 's2', 'l2', 'l1', 'k2', 'k1', &
279  're calculated', 'im calculated', 're reference', 'im reference'
280 
281  do ie = 1, size(k1)
282  do ipw = 1, size(l1)
283  do s2 = 1, 2
284 
285  ms = [ 1 ]
286  ss = [ (-1)**(s2+1), +1 ]
287  ls = [ l1(ipw), l2(ipw) ]
288  ks = [ k1(ie), k2(ie) ]
289 
290  integ = nested_coul_integ(z, a, c, 1, ms, ss, ls, ks)
291 
292  print '(SP,I3,SS,2I3,2F14.10,SP,4E16.7)', (-1)**(s2+1), l2(ipw), l1(ipw), k2(ie), k1(ie), &
293  real(integ, wp), aimag(integ), real(reference(s2,ipw,ie), wp), aimag(reference(s2,ipw,ie))
294  end do
295  end do
296  end do
297 
298  end subroutine test_1p_cint
299 
300 
309  subroutine test_1p_gint
310 
311  use multidip_integ, only: nested_cgreen_integ
312 
313  real(wp) :: a = 50, b = 150, c = 0, z = 1
314 
315  real(wp) :: k1(2) = [ 0.011080185297_wp, 0.680747949322_wp ]
316  real(wp) :: k2(2) = [ 1.069105413537_wp, 1.438602233160_wp ]
317 
318  integer :: l1(8) = [ 1, 1, 3, 1, 3, 3, 3, 3 ]
319  integer :: l2(8) = [ 0, 2, 2, 2, 2, 4, 2, 4 ]
320 
321  integer :: ie, ipw, s1, s2, sa, sb, ms(1), ls(2)
322 
323  complex(wp) :: ks(2), integ_ai, integ_bi, integ_ab
324 
325  print '(/,A)', 'Testing one-photon Green integrals (numerical integration)'
326  print '(/,4(1x,A),2x,A,12x,A,12x,A,3x,A,3x,A,3x,A)', 's2', 's1', 'l2', 'l1', 'k2', 'k1', &
327  're numerical ', 'im numerical ', 're asymptotic', 'im asymptotic'
328 
329  do ie = 1, size(k1)
330  do ipw = 1, size(l1)
331  do s2 = 1, 2
332  do s1 = 1, 2
333 
334  ms = [ 1 ]
335  ls = [ l1(ipw), l2(ipw) ]
336  ks = [ k1(ie), k2(ie) ]
337 
338  sa = (-1)**(s1 + 1)
339  sb = (-1)**(s2 + 1)
340 
341  integ_ai = nested_cgreen_integ(z, a, a, c, 1, sa, sb, ms, ls, ks)
342  integ_bi = nested_cgreen_integ(z, b, b, c, 1, sa, sb, ms, ls, ks)
343  integ_ab = nested_cgreen_integ(z, a, b, c, 1, sa, sb, ms, ls, ks)
344 
345  print '(SP,2I3,SS,2I3,2F14.10,SP,8E16.7)', sb, sa, l2(ipw), l1(ipw), k2(ie), k1(ie), &
346  real(integ_ab - integ_bi, wp), aimag(integ_ab - integ_bi), &
347  real(integ_ai - integ_bi, wp), aimag(integ_ai - integ_bi)
348 
349  end do
350  end do
351  end do
352  end do
353 
354  end subroutine test_1p_gint
355 
356 
363  subroutine test_2p_cint
364 
365  use multidip_integ, only: nested_coul_integ
366 
367  real(wp) :: a = 150, c = 0, z = 1
368 
369  real(wp) :: k1(2) = [ 0.011080185297_wp, 0.680747949322_wp ]
370  real(wp) :: k2(2) = [ 1.069105413537_wp, 1.438602233160_wp ]
371  real(wp) :: k3(2) = [ 1.511902774652_wp, 1.917220644579_wp ]
372 
373  integer :: l1(8) = [ 1, 1, 3, 1, 3, 3, 3, 3 ]
374  integer :: l2(8) = [ 0, 2, 2, 2, 2, 4, 2, 4 ]
375  integer :: l3(8) = [ 1, 1, 1, 3, 3, 3, 3, 3 ]
376 
377  integer :: ie, ipw, s4, ms(2), ss(4), ls(4)
378 
379  complex(wp) :: ks(4), integ1, integ2, integ
380  complex(wp) :: reference(2,8,2) = reshape([(-1.465883e+16_wp, -4.950701e+15_wp), (-3.842765e+16_wp, 8.620168e+16_wp), &
381  (-1.466146e+16_wp, -4.951551e+15_wp), (-3.843793e+16_wp, 8.623263e+16_wp), &
382  (-1.443664e+16_wp, -5.780723e+15_wp), (-3.378946e+16_wp, 8.865409e+16_wp), &
383  ( 1.520251e+16_wp, -2.902864e+15_wp), (-8.964058e+15_wp, -9.400575e+16_wp), &
384  ( 1.541476e+16_wp, -2.070261e+15_wp), (-1.420349e+16_wp, -9.382680e+16_wp), &
385  ( 1.542117e+16_wp, -2.071208e+15_wp), (-1.422236e+16_wp, -9.390101e+16_wp), &
386  ( 1.541476e+16_wp, -2.070261e+15_wp), (-1.420349e+16_wp, -9.382680e+16_wp), &
387  ( 1.542117e+16_wp, -2.071208e+15_wp), (-1.422236e+16_wp, -9.390101e+16_wp), &
388  (-4.322904e+03_wp, -4.890496e+03_wp), (-2.301503e+04_wp, 1.829823e+04_wp), &
389  (-4.323345e+03_wp, -4.890988e+03_wp), (-2.301886e+04_wp, 1.830140e+04_wp), &
390  ( 6.409175e+03_wp, -1.254518e+03_wp), ( 2.742317e+04_wp, 1.058896e+04_wp), &
391  ( 5.914821e+03_wp, 2.763232e+03_wp), ( 1.382079e+04_wp, -2.596191e+04_wp), &
392  (-5.378829e+03_wp, 3.705000e+03_wp), (-2.937519e+04_wp, 1.216512e+03_wp), &
393  (-5.380092e+03_wp, 3.705885e+03_wp), (-2.938678e+04_wp, 1.217206e+03_wp), &
394  (-5.378829e+03_wp, 3.705000e+03_wp), (-2.937519e+04_wp, 1.216512e+03_wp), &
395  (-5.380092e+03_wp, 3.705885e+03_wp), (-2.938678e+04_wp, 1.217206e+03_wp)], &
396  [2, 8, 2])
397 
398  print '(/,A)', 'Testing two-photon Coulomb integrals'
399  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', &
400  're calculated', 'im calculated', 're reference', 'im reference'
401 
402  do ie = 1, size(k1)
403  do ipw = 1, size(l1)
404  do s4 = 1, 2
405 
406  ms = [ 1, 1 ]
407 
408  ss = [ (-1)**(s4+1), +1, -1, +1 ]
409  ls = [ l1(ipw), l2(ipw), l2(ipw), l3(ipw) ]
410  ks = [ k1(ie), k2(ie), k2(ie), k3(ie) ]
411 
412  integ1 = nested_coul_integ(z, a, c, 2, ms, ss, ls, ks)
413 
414  ss = [ +1, +1, -1, (-1)**(s4+1) ]
415  ls = [ l3(ipw), l2(ipw), l2(ipw), l1(ipw) ]
416  ks = [ k3(ie), k2(ie), k2(ie), k1(ie) ]
417 
418  integ2 = nested_coul_integ(z, a, c, 2, ms, ss, ls, ks)
419 
420  integ = integ1 + integ2
421 
422  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), &
423  real(integ, wp), aimag(integ), real(reference(s4,ipw,ie), wp), aimag(reference(s4,ipw,ie))
424  end do
425  end do
426  end do
427 
428  end subroutine test_2p_cint
429 
430 
437  subroutine test_2p_gint
438 
439  use multidip_integ, only: nested_cgreen_integ
440 
441  real(wp) :: a = 150, c = 0, z = 1
442 
443  real(wp) :: k1(2) = [ 0.011080185297_wp, 0.680747949322_wp ]
444  real(wp) :: k2(2) = [ 1.069105413537_wp, 1.438602233160_wp ]
445  real(wp) :: k3(2) = [ 1.511902774652_wp, 1.917220644579_wp ]
446 
447  integer :: l1(8) = [ 1, 1, 3, 1, 3, 3, 3, 3 ]
448  integer :: l2(8) = [ 0, 2, 2, 2, 2, 4, 2, 4 ]
449  integer :: l3(8) = [ 1, 1, 1, 3, 3, 3, 3, 3 ]
450 
451  integer :: ie, ipw, s4, ms(2), ls(3)
452 
453  complex(wp) :: ks(3), integ1, integ2
454  complex(wp) :: reference(2,2,8,2) = reshape([(-8.310927e+15_wp, 7.523384e+15_wp), (-5.909915e+16_wp, 6.931934e+16_wp), &
455  ( 5.679523e+16_wp, 4.672527e+16_wp), ( 8.225897e+24_wp, 1.848188e+24_wp), &
456  (-6.506509e+15_wp, 2.066669e+16_wp), (-1.618385e+16_wp, 4.224605e+15_wp), &
457  ( 1.050096e+17_wp, 4.554189e+16_wp), (-4.312479e+24_wp, -7.240598e+24_wp), &
458  (-7.673705e+15_wp, 2.037601e+16_wp), (-1.647265e+16_wp, 3.342092e+15_wp), &
459  ( 1.078882e+17_wp, 3.987695e+16_wp), (-4.728274e+24_wp, -7.026109e+24_wp), &
460  (-4.503787e+15_wp, -2.119630e+16_wp), ( 1.619053e+16_wp, 4.287343e+15_wp), &
461  (-1.138495e+17_wp, 1.201762e+16_wp), ( 2.556995e+23_wp, 8.432496e+24_wp), &
462  (-3.344463e+15_wp, -2.151751e+16_wp), ( 1.600768e+16_wp, 5.199000e+15_wp), &
463  (-1.135684e+17_wp, 1.836697e+16_wp), ( 7.238276e+23_wp, 8.446832e+24_wp), &
464  (-8.612323e+15_wp, -1.723199e+16_wp), ( 3.535435e+16_wp, 4.008848e+16_wp), &
465  (-9.807095e+16_wp, 3.755844e+16_wp), ( 7.275593e+24_wp, 4.341652e+24_wp), &
466  (-3.344463e+15_wp, -2.151751e+16_wp), ( 1.600768e+16_wp, 5.199000e+15_wp), &
467  (-1.135684e+17_wp, 1.836697e+16_wp), ( 7.238276e+23_wp, 8.446832e+24_wp), &
468  (-8.612323e+15_wp, -1.723199e+16_wp), ( 3.535435e+16_wp, 4.008848e+16_wp), &
469  (-9.807095e+16_wp, 3.755844e+16_wp), ( 7.275593e+24_wp, 4.341652e+24_wp), &
470  (-1.247516e+03_wp, 3.118370e+03_wp), (-2.763140e+04_wp, 1.784564e+04_wp), &
471  ( 1.332880e+04_wp, 9.920923e+03_wp), ( 6.345060e+03_wp, 4.578023e+04_wp), &
472  (-3.896753e+03_wp, 9.080720e+02_wp), (-6.662204e+03_wp, 3.014237e+04_wp), &
473  ( 6.717303e+03_wp, 1.712697e+04_wp), ( 4.397613e+04_wp, -1.193480e+04_wp), &
474  ( 1.189054e+03_wp, -3.822275e+03_wp), (-2.263210e+04_wp, -2.099869e+04_wp), &
475  ( 1.136951e+04_wp, -1.445428e+04_wp), (-3.253996e+04_wp, -3.189015e+04_wp), &
476  ( 3.212357e+03_wp, -2.385969e+03_wp), ( 1.813020e+04_wp, -2.499891e+04_wp), &
477  (-1.298806e+04_wp, -1.303291e+04_wp), (-4.510261e+04_wp, -6.582287e+03_wp), &
478  ( 4.327121e+02_wp, 3.979850e+03_wp), ( 1.239242e+04_wp, 2.828936e+04_wp), &
479  (-4.667920e+03_wp, 1.779027e+04_wp), ( 1.713954e+04_wp, 4.222912e+04_wp), &
480  ( 1.147174e+03_wp, 5.355515e+03_wp), ( 1.429665e+04_wp, 1.753491e+04_wp), &
481  (-5.010378e+03_wp, 2.217364e+04_wp), ( 3.781810e+04_wp, 1.923663e+04_wp), &
482  ( 4.327121e+02_wp, 3.979850e+03_wp), ( 1.239242e+04_wp, 2.828936e+04_wp), &
483  (-4.667920e+03_wp, 1.779027e+04_wp), ( 1.713954e+04_wp, 4.222912e+04_wp), &
484  ( 1.147174e+03_wp, 5.355515e+03_wp), ( 1.429665e+04_wp, 1.753491e+04_wp), &
485  (-5.010378e+03_wp, 2.217364e+04_wp), ( 3.781810e+04_wp, 1.923663e+04_wp)], &
486  [2,2,8,2])
487 
488  print '(/,A)', 'Testing two-photon Green integrals'
489  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', &
490  're calculated', 'im calculated', 're reference', 'im reference'
491 
492  do ie = 1, size(k1)
493  do ipw = 1, size(l1)
494  do s4 = 1, 2
495 
496  ms = [ 1, 1 ]
497  ls = [ l1(ipw), l2(ipw), l3(ipw) ]
498  ks = [ k1(ie), k2(ie), k3(ie) ]
499 
500  integ1 = nested_cgreen_integ(z, a, a, c, 2, (-1)**(s4+1), +1, ms, ls, ks)
501  integ2 = nested_cgreen_integ(z, a, a, c, 2, (-1)**(s4+1), -1, ms, ls, ks)
502 
503  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),&
504  real(integ1, wp), aimag(integ1), real(reference(1,s4,ipw,ie), wp), aimag(reference(1,s4,ipw,ie))
505  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),&
506  real(integ2, wp), aimag(integ2), real(reference(2,s4,ipw,ie), wp), aimag(reference(2,s4,ipw,ie))
507  end do
508  end do
509  end do
510  end subroutine test_2p_gint
511 
512 
519  subroutine test_ang_dist
520 
521  use dipelm_special_functions, only: threej
522  use multidip_params, only: rone
523  use multidip_special, only: beta_contraction_tensor, beta_2p_demekhin, beta_2p_arb_pol_sph_harm
524 
525  integer :: J = 0, li = 1, lj = 1, mi = 0, mj = 0, ki = 0, kj = 0, qi = 0, qj = 0, p = 0
526  real(wp) :: beta, ref, beta_ap
527 
528  print '(/,A)', 'Testing one-photon angular integrals'
529  print '(/,8A4,2A12)', 'J', 'p', 'li', 'lj', 'mi', 'mj', 'qi', 'qj', 'TJij', 'reference'
530 
531  do j = 0, 2
532  do p = 0, +1
533  do mi = -1, +1
534  do mj = -1, +1
535  do qi = -1, +1
536  do qj = -1, +1
537  beta = beta_contraction_tensor(j, 1, [ p ], li, mi, [ qi ], lj, mj, [ qj ])
538  ref = (-1)**(p + mi + qi) * (2*j + 1) * sqrt((2*li + rone)*(2*lj + rone)) &
539  * threej(2*li, 0, 2*lj, 0, 2*j, 0) &
540  * threej(2*1, 2*p, 2*1, -2*p, 2*j, 0) &
541  * threej(2*li, 2*mi, 2*lj, -2*mj, 2*j, -2*(mi-mj)) &
542  * threej(2*1, 2*qi, 2*1, -2*qj, 2*j, -2*(mi-mj))
543  if (abs(beta) > 1e-10 .or. abs(ref) > 1e-10) then
544  print '(I4,SP,I4,SS,2I4,SP,4I4,2F12.7,A8)', &
545  j, p, li, lj, mi, mj, qi, qj, beta, ref, &
546  merge(' ok', 'WRONG!', abs(beta - ref) / (abs(beta) + abs(ref)) < 1e-5)
547  end if
548  end do
549  end do
550  end do
551  end do
552  end do
553  end do
554 
555  print '(/,A)', 'Testing two-photon angular integrals'
556  print '(/,10A4,2A12)', 'J', 'p', 'li', 'lj', 'mi', 'mj', 'ki', 'kj', 'qi', 'qj', 'TJij', 'reference'
557 
558  do j = 0, 2
559  do p = 0, +1
560  do mi = -1, +1
561  do mj = -1, +1
562  do ki = -1, +1
563  do kj = -1, +1
564  do qi = -1, +1
565  do qj = -1, +1
566  beta = beta_contraction_tensor(j, 2, [ p, p ], li, mi, [ ki, qi ], lj, mj, [ kj, qj ])
567  ref = beta_2p_demekhin(j, p, p, li, mi, ki, qi, lj, mj, kj, qj)
568  if (abs(beta) > 1e-10 .or. abs(ref) > 1e-10) then
569  print '("1st test ",I4,SP,I4,SS,2I4,SP,6I4,2F12.7,A8)', &
570  j, p, li, lj, mi, mj, ki, kj, qi, qj, beta, ref, &
571  merge(' ok', 'WRONG!', abs(beta - ref) / (abs(beta) + abs(ref)) < 1e-5)
572  end if
573 
574  beta_ap = beta_2p_arb_pol_sph_harm&
575  &(j, 0, 2, [ p, p ], li, mi, [ ki, qi ], [ p, p ], lj, mj, [ kj, qj ])
576  if (abs(beta_ap) > 1e-10 .or. abs(ref) > 1e-10) then
577  print '("2nd test ",I4,SP,I4,SS,2I4,SP,6I4,2F12.7,A8)', &
578  j, p, li, lj, mi, mj, ki, kj, qi, qj, beta_ap, ref, &
579  merge(' ok', 'WRONG!', abs(beta_ap - ref) / (abs(beta_ap) + abs(ref)) < 1e-5)
580  end if
581  end do
582  end do
583  end do
584  end do
585  end do
586  end do
587  end do
588  end do
589 
590  end subroutine test_ang_dist
591 
592 end 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:486
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:211
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:370
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::rone
real(wp), parameter rone
Definition: multidip_params.F90:57
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_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:1123
multidip_special::coul_ukrmol
subroutine coul_ukrmol(Z, l, Ek, r, F, Fp, G, Gp)
Coulomb functions (UKRmol)
Definition: multidip_special.F90:661
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_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:991
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:1191
multidip_special::next_permutation
logical function next_permutation(p)
Construct or advance permutation.
Definition: multidip_special.F90:844
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:438
multidip_tests::test_levin
subroutine test_levin
Tests of the Levin quadrature routine.
Definition: multidip_tests.f90:131
multidip_tests::test_coul
subroutine test_coul
Negative-energy Coulomb function tests.
Definition: multidip_tests.f90:93
multidip_tests::test_2p_cint
subroutine test_2p_cint
Two-photon Coulomb integral unit tests.
Definition: multidip_tests.f90:364
multidip_tests::test_permutations
subroutine test_permutations
Permutations unit tests.
Definition: multidip_tests.f90:67
multidip_tests::test_1p_eint
subroutine test_1p_eint
One-photon exponential integral unit tests.
Definition: multidip_tests.f90:199
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:310
multidip_tests::test_1p_cint
subroutine test_1p_cint
One-photon Coulomb integral unit tests.
Definition: multidip_tests.f90:247
multidip_tests::test_ang_dist
subroutine test_ang_dist
Angular algebra unit tests.
Definition: multidip_tests.f90:520