multidip_levin Module Reference

Levin quadrature for numerical integration. More...

Functions/Subroutines
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) More...
complex(wp) function	levin_adapt (Z, ra, rb, c, m, s1, s2, l1, l2, k1, k2)
	One-dimensional adaptive Levin integration. More...
subroutine	levin_prepare (Z, ra, rb, s1, s2, l1, l2, k1, k2, depth, converged, Hl_buffer)
	Precalculate Coulomb-Hankel functions for subsequent Levin integration. More...
subroutine	levin_eval (Z, r, s1, s2, l1, l2, k1, k2, Hl)
	Evaluate a pair of Coulomb-Hankel functions in given point. More...
recursive subroutine	levin_improve (Z, ra, rb, c, m, l1, l2, k1, k2, depth, interval_idx, tgt_depth, converged, Hl, estimates)
	Perform one subdivision iteration of adaptive Levin integration and update estimates. More...
complex(wp) function	levin_integrate (Z, ra, rb, c, m, l1, l2, k1, k2, Hla, Hlb)
	Fixed-order Levin integration. More...
complex(wp) function	levin_integrate_2x2 (Z, ra, rb, c, m, lc, lo, kc, ko, wa, wb)
	Fixed-order Levin integration (dim 2) More...
complex(wp) function	levin_integrate_4x4 (Z, ra, rb, c, m, l1, l2, k1, k2, Hla, Hlb)
	Fixed-order Levin integration (dim 4) More...

Detailed Description

Levin quadrature for numerical integration.

Author

J Benda

Date

2021

This module implements the efficient Levin quadrature for numerical integration of oscillatory functions. It is used by multidip_integ to integrate products of Coulomb-Hankel functions. Its advantage lies particularly in the fact that it only requires evaluation of the Coulomb functions at the endpoints of the integration range; within the integration range it works with recurrence formulas only. Given how expensive the evaluation of the Coulomb-Hankel function is, using Levin integration in place of the straightforward Romberg integration speeds up the calculations by several orders of magnitude.

See the following papers for more details on the theory:

• D. Levin, Fast integration of rapidly oscillatory functions, J. Comput. Appl. Math. 67 (1996) 95-101.
• J. L. Powell, Recurrence formulas for Coulomb functions, Phys. Rev. 72 (1947) 626.

Function/Subroutine Documentation

levin_adapt()

complex(wp) function multidip_levin::levin_adapt	(	real(wp), intent(in)	Z,
		real(wp), intent(in)	ra,
		real(wp), intent(in)	rb,
		real(wp), intent(in)	c,
		integer, intent(in)	m,
		integer, intent(in)	s1,
		integer, intent(in)	s2,
		integer, intent(in)	l1,
		integer, intent(in)	l2,
		complex(wp), intent(in)	k1,
		complex(wp), intent(in)	k2
	)

One-dimensional adaptive Levin integration.

Author

J Benda

Date

2021 - 2023

As needed, recursively subdivide the complete integration range and apply a fixed-order Levin quadrature to the subintervals until further subdivision changes the resulting integral only a little (within tolerance).

A scheme illustrating the indexing of nodes at different subdivision levels:

  (ra)                                    (rb)
   1                                       0
   |---------------------------------------|              (depth 0)

   2                   3                   0
   |-------------------|-------------------|              (depth 1)

   4         5         6         7         0
   |---------|---------|---------|---------|              (depth 2)

   8    9   10   11   12   13   14   15    0
   |----|----|----|----|----|----|----|----|              (depth 3)

   ... etc ...

The advantage of this scheme is that Nodes whose indices differ only by an even multiplicative factor are equivalent (collocated) and can be uniquely indexed by the further irreducible odd number obtained by repeated division by two. The repeated division by two can be efficiently implemented as a right bit shift (shiftr) by offset given by the number of trailing zeros (trailz) in the binary representation of the index, often together amounting to just two machine instructions.

Intervals have the same index as their left point. That is, the top-level interval is the interval 1, consisting of two subintervals 2 and 3, each of which is further composed of further subintervals.

The positive-energy Coulomb wave functions are evaluated only in nodes adjacent to intervals where the integral has not yet converged. The negative-energy Coulomb functions as well as the rest of the integrand is repeatedly evaluated in Chebyshev nodes within each subinterval, see levin_integrate_2x2 and levin_integrate_4x4.

Parameters

Z	Residual ion charge
ra	Lower bound
rb	Upper bound
c	Damping factor (additional exp(-c*r) added to all r^m functions)
m	Radial coordinate power
s1	Sign of the first Coulomb-Hankel function
s2	Sign of the second Coulomb-Hankel function
l1	Angular momentum of the first Coulomb-Hankel function
l2	Angular momentum of the second Coulomb-Hankel function
k1	Complex momentum of the first Coulomb-Hankel function
k2	Complex momentum of the second Coulomb-Hankel function

Definition at line 133 of file multidip_levin.f90.

Here is the call graph for this function:

Here is the caller graph for this function:

levin_eval()

subroutine multidip_levin::levin_eval	(	real(wp), intent(in)	Z,
		real(wp), intent(in)	r,
		integer, intent(in)	s1,
		integer, intent(in)	s2,
		integer, intent(in)	l1,
		integer, intent(in)	l2,
		complex(wp), intent(in)	k1,
		complex(wp), intent(in)	k2,
		complex(wp), dimension(2, 2), intent(inout)	Hl
	)

Evaluate a pair of Coulomb-Hankel functions in given point.

Author

J Benda

Date

2021 - 2023

Evaluate elements of the Coulomb-Hankel functions matrix:

\[ \mathsf{H} = \left(\begin{matrix} H_{l_1 }^{s_1}(-1/k_1, k_1 r) & H_{l_2 }^{s_2}(-1/k_2, k_2 r) \\ H_{l_1 + 1}^{s_1}(-1/k_1, k_1 r) & H_{l_2 + 1}^{s_2}(-1/k_2, k_2 r) \end{matrix}\right) \]

This subroutine will not distinguish between positive- and negative-energy Coulomb functions (the input momenta are complex), but only the positive-energy values are actually later used in the rest of the code.

Parameters

Z	Residual ion charge
r	Evaluation radius
s1	Sign of the first Coulomb-Hankel function
s2	Sign of the second Coulomb-Hankel function
l1	Angular momentum of the first Coulomb-Hankel function
l2	Angular momentum of the second Coulomb-Hankel function
k1	Complex momentum of the first Coulomb-Hankel function
k2	Complex momentum of the second Coulomb-Hankel function
Hl	Coulomb-Hankel functions evaluated at the evaluation radius (for l1, l1+1, l2 and l2+1)

Definition at line 285 of file multidip_levin.f90.

Here is the call graph for this function:

Here is the caller graph for this function:

levin_improve()

recursive subroutine multidip_levin::levin_improve	(	real(wp), intent(in)	Z,
		real(wp), intent(in)	ra,
		real(wp), intent(in)	rb,
		real(wp), intent(in)	c,
		integer, intent(in)	m,
		integer, intent(in)	l1,
		integer, intent(in)	l2,
		complex(wp), intent(in)	k1,
		complex(wp), intent(in)	k2,
		integer, intent(in)	depth,
		integer, intent(in)	interval_idx,
		integer, intent(in)	tgt_depth,
		logical, dimension(:), intent(inout)	converged,
		complex(wp), dimension(2, 2, *), intent(in)	Hl,
		complex(wp), dimension(:), intent(inout)	estimates
	)

Perform one subdivision iteration of adaptive Levin integration and update estimates.

Author

J Benda

Date

2021 - 2023

Recurse to the deepest target depth and evaluate integral estimates in all sub-intervals marked as "not converged". Use these new values to improve estimates of all parent intervals (on all parent levels). If the updated estimate of any interval changes within the builtin tolerance, mark it (and its subintervals) as converged.

This subroutine checks the convergence using a relative and absolute tolerance. For convergence, either the relative change between the base and improved estimate has to meet the relative tolerance, or the absolute difference between these two estimates has to meet the absolute tolerance. The absolute tolerance is calculated as a product of the relative tolerance and the estimate of the top-level integral.

Parameters

Z	Residual ion charge
ra	Lower bound
rb	Upper bound
c	Damping factor (additional exp(-c*r) added to all r^m functions)
m	Radial coordinate power
l1	Angular momentum of the first Coulomb-Hankel function
l2	Angular momentum of the second Coulomb-Hankel function
k1	Complex momentum of the first Coulomb-Hankel function
k2	Complex momentum of the second Coulomb-Hankel function
depth	Subdivision depth that gave rise to the current integration interval
interval_idx	Index of the current integration interval
tgt_depth	Current subdivision depth resulting in 2**depth subintervals
converged	Logical flags per interval at all subdivision depths indicating whether its value needs improving
Hl	Coulomb-Hankel functions evaluated at unique subdivision nodes (for l1, l1+1, l2 and l2+1)
estimates	Integral estimates for all intervals in all subdivision levels

Definition at line 336 of file multidip_levin.f90.

Here is the call graph for this function:

Here is the caller graph for this function:

levin_integrate()

complex(wp) function multidip_levin::levin_integrate	(	real(wp), intent(in)	Z,
		real(wp), intent(in)	ra,
		real(wp), intent(in)	rb,
		real(wp), intent(in)	c,
		integer, intent(in)	m,
		integer, intent(in)	l1,
		integer, intent(in)	l2,
		complex(wp), intent(in)	k1,
		complex(wp), intent(in)	k2,
		complex(wp), dimension(2, 2), intent(in)	Hla,
		complex(wp), dimension(2, 2), intent(in)	Hlb
	)

Fixed-order Levin integration.

Author

J Benda

Date

2021 - 2024

Integrate product of two Coulomb-Hankel functions, coordinate power and radial exponential,

\[ I = \int\limits_{r_a}^{r_b} H_{l_1}^{s_1}(\eta_1, k_1 r) r^{m} \mathrm{e}^{-cr} H_{l_2}^{s_2}(\eta_2, k_2 r) \mathrm{d}r \]

(accepting both positive- and negative-energy Coulomb functions) using the method from "D. Levin, Fast integration of rapidly oscillatory functions, J. Comput. Appl. Math. 67 (1996) 95-101".

Coulomb functions at end points are not evaluated in this function and need to be provided via arguments Hla and Hlb, see levin_eval.

Parameters

Z	Residual ion charge
ra	Lower bound
rb	Upper bound
c	Damping factor (additional exp(-c*r) added to all r^m functions)
m	Radial coordinate power
l1	Angular momentum of the first Coulomb-Hankel function
l2	Angular momentum of the second Coulomb-Hankel function
k1	Complex momentum of the first Coulomb-Hankel function
k2	Complex momentum of the second Coulomb-Hankel function
Hla	Coulomb-Hankel functions evaluated at ra (for l1, l1+1, l2 and l2+1)
Hlb	Coulomb-Hankel functions evaluated at rb (for l1, l1+1, l2 and l2+1)

Definition at line 442 of file multidip_levin.f90.

Here is the call graph for this function:

Here is the caller graph for this function:

levin_integrate_2x2()

complex(wp) function multidip_levin::levin_integrate_2x2	(	real(wp), intent(in)	Z,
		real(wp), intent(in)	ra,
		real(wp), intent(in)	rb,
		real(wp), intent(in)	c,
		integer, intent(in)	m,
		integer, intent(in)	lc,
		integer, intent(in)	lo,
		real(wp), intent(in)	kc,
		real(wp), intent(in)	ko,
		complex(wp), dimension(2), intent(in)	wa,
		complex(wp), dimension(2), intent(in)	wb
	)

Fixed-order Levin integration (dim 2)

Author

J Benda

Date

2021 - 2024

Integrate product of a Coulomb-Hankel function, decreasing real Whittaker function, coordinate power and radial exponential,

\[ I = \int\limits_{r_a}^{r_b} H_{l_1}^{s_1}(\eta_1, k_1 r) r^{m} \mathrm{e}^{-cr} W_{1/|k_2|,l_2+1/2}(2 |k_2| r) \mathrm{d}r \]

using the method from "D. Levin, Fast integration of rapidly oscillatory functions, J. Comput. Appl. Math. 67 (1996) 95-101".

The Levin's objective function \( \mathbf{w}(r) \) has the following two components, where the first of these it the one of interest:

\[ \mathbf{w}(r) = \left(\begin{matrix} H_{l_1 }^{s_1}(\eta_1, k_1 r) \\ H_{l_1 + 1}^{s_1}(\eta_1, k_1 r) \end{matrix}\right) . \]

The Levin matrix \( \mathsf{A}(r) = \mathsf{a}/r + \mathsf{b} \) is constructed from the known Coulomb function recursion relations,

\[ \frac{\mathrm{d}}{\mathrm{d}r} H_{l+1}(\eta, kr) = k R_{l+1}(\eta) H_{l}(\eta, kr) - k S_{l+1}(\eta, kr) H_{l+1}(\eta, kr) \,, \]

\[ \frac{\mathrm{d}}{\mathrm{d}r} H_{l}(\eta, kr) = k S_{l+1}(\eta, kr) H_{l}(\eta, kr) - k R_{l+1}(\eta) H_{l+1}(\eta, kr) \,. \]

where

\[ R_{l}(\eta) = \sqrt{1 + \eta^2/l^2} \,, \qquad S_{l}(\eta, \rho) = l/\rho + \eta/l \,. \]

See for example "J. L. Powell, Recurrence formulas for Coulomb functions, Phys. Rev. 72 (1947) 626" or DLMF §33.4.

The two-component Levin auxiliary non-oscillatory function \( \mathbf{p}(r) \) is expanded in Chebyshev polynomials of order 5 and the expansion coefficients obtained from collocation equation of rank 6. Collocation points are chosen to be idential to Chebyshev nodes. This should result in the best possible interpolation of \( \mathbf{p}(r) \) for given Chebyshev order.

Coulomb functions at end points are not evaluated in this function and need to be provided via arguments wa and wb.

Parameters

Z	Residual ion charge
ra	Lower bound
rb	Upper bound
c	Damping factor (additional exp(-c*r) added to all r^m functions)
m	Radial coordinate power
lc	Angular momentum of the negative-energy Coulomb-Hankel function
lo	Angular momentum of the positive-energy Coulomb-Hankel function
kc	Magnitude of the imaginary momentum of the negative-energy Coulomb-Hankel function
ko	Magnitude of the real momentum of the positive-energy Coulomb-Hankel function
wa	Positive-energy Coulomb-Hankel function evaluated at ra (for lo and lo+1)
wb	Positive-energy Coulomb-Hankel function evaluated at rb (for lo and lo+1)

Definition at line 524 of file multidip_levin.f90.

Here is the call graph for this function:

Here is the caller graph for this function:

levin_integrate_4x4()

complex(wp) function multidip_levin::levin_integrate_4x4	(	real(wp), intent(in)	Z,
		real(wp), intent(in)	ra,
		real(wp), intent(in)	rb,
		real(wp), intent(in)	c,
		integer, intent(in)	m,
		integer, intent(in)	l1,
		integer, intent(in)	l2,
		real(wp), intent(in)	k1,
		real(wp), intent(in)	k2,
		complex(wp), dimension(2, 2), intent(in)	Hla,
		complex(wp), dimension(2, 2), intent(in)	Hlb
	)

Fixed-order Levin integration (dim 4)

Author

J Benda

Date

2021 - 2024

Integrate product of two Coulomb-Hankel functions, coordinate power and radial exponential,

\[ I = \int\limits_{r_a}^{r_b} H_{l_1}^{s_1}(\eta_1, k_1 r) r^{m} \mathrm{e}^{-cr} H_{l_2}^{s_2}(\eta_2, k_2 r) \mathrm{d}r \]

using the method from "D. Levin, Fast integration of rapidly oscillatory functions, J. Comput. Appl. Math. 67 (1996) 95-101".

The Levin's objective function \( \mathbf{w}(r) \) has the following four components, where the first of these it the one of interest:

\[ \mathbf{w}(r) = \left(\begin{matrix} H_{l_1 }^{s_1}(\eta_1, k_1 r) H_{l_2 }^{s_2}(\eta_2, k_2 r) \\ H_{l_1 }^{s_1}(\eta_1, k_1 r) H_{l_2 + 1}^{s_2}(\eta_2, k_2 r) \\ H_{l_1 + 1}^{s_1}(\eta_1, k_1 r) H_{l_2 }^{s_2}(\eta_2, k_2 r) \\ H_{l_1 + 1}^{s_1}(\eta_1, k_1 r) H_{l_2 + 1}^{s_2}(\eta_2, k_2 r) \end{matrix}\right) . \]

The Levin matrix \( \mathsf{A}(r) = \mathsf{a}/r + \mathsf{b} \) is constructed from the known Coulomb function recursion relations,

\[ \frac{\mathrm{d}}{\mathrm{d}r} H_{l+1}(\eta, kr) = k R_{l+1}(\eta) H_{l}(\eta, kr) - k S_{l+1}(\eta, kr) H_{l+1}(\eta, kr) \,, \]

\[ \frac{\mathrm{d}}{\mathrm{d}r} H_{l}(\eta, kr) = k S_{l+1}(\eta, kr) H_{l}(\eta, kr) - k R_{l+1}(\eta) H_{l+1}(\eta, kr) \,. \]

where

\[ R_{l}(\eta) = \sqrt{1 + \eta^2/l^2} \,, \qquad S_{l}(\eta, \rho) = l/\rho + \eta/l \,. \]

See for example "J. L. Powell, Recurrence formulas for Coulomb functions, Phys. Rev. 72 (1947) 626" or DLMF §33.4.

The four-component Levin auxiliary non-oscillatory function \( \mathbf{p}(r) \) is expanded in Chebyshev polynomials of order 5 and the expansion coefficients obtained from collocation equation of rank 6. Collocation points are chosen to be idential to Chebyshev nodes. This should result in the best possible interpolation of \( \mathbf{p}(r) \) for given Chebyshev order.

Coulomb functions at end points are not evaluated in this function and need to be provided via arguments Hla and Hlb.

Parameters

Z	Residual ion charge.
ra	Lower bound
rb	Upper bound
c	Damping factor (additional exp(-c*r) added to all r^m functions)
m	Radial coordinate power
l1	Angular momentum of the first Coulomb-Hankel function
l2	Angular momentum of the second Coulomb-Hankel function
k1	Complex momentum of the first Coulomb-Hankel function
k2	Complex momentum of the second Coulomb-Hankel function
Hla	Coulomb-Hankel functions evaluated at ra (for l1, l1+1, l2 and l2+1)
Hlb	Coulomb-Hankel functions evaluated at rb (for l1, l1+1, l2 and l2+1)

Definition at line 720 of file multidip_levin.f90.

Here is the caller graph for this function:

levin_prepare()

subroutine multidip_levin::levin_prepare	(	real(wp), intent(in)	Z,
		real(wp), intent(in)	ra,
		real(wp), intent(in)	rb,
		integer, intent(in)	s1,
		integer, intent(in)	s2,
		integer, intent(in)	l1,
		integer, intent(in)	l2,
		complex(wp), intent(in)	k1,
		complex(wp), intent(in)	k2,
		integer, intent(in)	depth,
		logical, dimension(:), intent(inout)	converged,
		complex(wp), dimension(2, 2, *), intent(inout)	Hl_buffer
	)

Precalculate Coulomb-Hankel functions for subsequent Levin integration.

Author

J Benda

Date

2021 - 2023

Set up the given subdivision depth for subsequent use. This involves evaluation of the Coulomb functions at the endpoints of the subintervals of the not-yet-converged parent intervals, as well as marking these subintervals of the not-yet-converged parent intervals as "not converged", while marking subintervals of the already converged parent intervals as "converged" (and not evaluating Coulomb functions for their sake).

Parameters

Z	Residual ion charge
ra	Lower bound
rb	Upper bound
s1	Sign of the first Coulomb-Hankel function
s2	Sign of the second Coulomb-Hankel function
l1	Angular momentum of the first Coulomb-Hankel function
l2	Angular momentum of the second Coulomb-Hankel function
k1	Complex momentum of the first Coulomb-Hankel function
k2	Complex momentum of the second Coulomb-Hankel function
depth	Current subdivision depth resulting in 2**depth subintervals
converged	Logical flags per interval at all subdivision depths indicating whether its value needs improving
Hl_buffer	Coulomb-Hankel functions evaluated at unique subdivision nodes (for l1, l1+1, l2 and l2+1)

Definition at line 204 of file multidip_levin.f90.

Here is the call graph for this function:

Here is the caller graph for this function:

nested_cgreen_correct_levin()

subroutine multidip_levin::nested_cgreen_correct_levin	(	real(wp), intent(in)	Z,
		real(wp), intent(in)	a,
		real(wp), intent(in)	b,
		real(wp), intent(in)	c,
		integer, intent(in)	N,
		integer, intent(in)	sa,
		integer, intent(in)	sb,
		integer, dimension(:), intent(in)	m,
		integer, dimension(:), intent(in)	l,
		complex(wp), dimension(:), intent(in)	k,
		complex(wp), intent(inout)	v
	)

Numerically correct asymptotic approximation of the Coulomb-Green integral (Levin integration)

Author

J Benda

Date

2021 - 2023

This function computes the integral of the Coulomb-Green's integrand over Q = (a,+∞)^N \ (b,+∞)^N numerically. Currently, Levin integration based on the Chebyshev interpolation is used.

Parameters

Z	Residual ion charge
a	Lower bound
b	Upper bound
c	Damping factor (additional exp(-c*r) added to all r^m functions)
N	Dimension (number of integration variables)
sa	Sign of the right-most Coulomb-Hankel function
sb	Sign of the left-most Coulomb-Hankel function
m	Array of integer powers of length N
l	Array of angular momenta of length N + 1
k	Array of linear momenta of length N + 1
v	Value of the asymptotic integral to correct (updated on exit from this subroutine)

Warning

This is currently only implemented for one-dimensional case containing no Green's function at all.

Todo:

Generalize for arbitrary dimension!

Definition at line 69 of file multidip_levin.f90.

Here is the call graph for this function:

Here is the caller graph for this function: