multidip_integ.f90

Go to the documentation of this file.

! Copyright 2020
!
! For a comprehensive list of the developers that contributed to these codes
! see the UK-AMOR website.
!
! This file is part of UKRmol-out (UKRmol+ suite).
!
!     UKRmol-out is free software: you can redistribute it and/or modify
!     it under the terms of the GNU General Public License as published by
!     the Free Software Foundation, either version 3 of the License, or
!     (at your option) any later version.
!
!     UKRmol-out is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!     GNU General Public License for more details.
!
!     You should have received a copy of the GNU General Public License
!     along with  UKRmol-out (in source/COPYING). Alternatively, you can also visit
!     <https://www.gnu.org/licenses/>.
!
!> \brief   Special integrals needed by MULTIDIP
!> \author  J Benda
!> \date    2020 - 2023
!>
!> Module containing routines for calculation of the multi-dimensional dipole integrals used in outer correction
!> of transition dipole elements in MULTIDIP.
!>
module multidip_integ
 
    use iso_fortran_env, only: error_unit
    use multidip_util,   only: find_column, append_column
    use phys_const_gbl,  only: pi, imu
    use precisn_gbl,     only: wp
 
    implicit none
 
contains
 
    !> \brief   Multi-dimensional triangular integral of exponentials and powers
    !> \author  J Benda
    !> \date    2020 - 2023
    !>
    !> Evaluate the nested many-dimensional triangular integral
    !> \f[
    !>    \int\limits_a^{+\infty} \dots \int\limits_{r_3}^{+\infty} \mathrm{e}^{\mathrm{i}s_4 \theta_4(r_2)}
    !>    r_2^{m_2} \mathrm{e}^{\mathrm{i}s_3 \theta_3(r_2)} \int\limits_{r_2}^{+\infty}
    !>    \mathrm{e}^{\mathrm{i}s_2 \theta_2(r_1)} r_1^{m_1} \mathrm{e}^{\mathrm{i}s_1 \theta_1(r_1)}
    !>    \mathrm{d}r_1 \mathrm{d}r_2 \mathrm{d}r_3 \dots
    !> \f]
    !> where
    !> \f[
    !>    \theta_1(r_1) = k_1 r_1 + \log(2 k_1 r_1)/k_1 - \pi l_1/2 + \sigma_{l_1}(k_1)
    !> \f]
    !> is the asymptotic phase of a Coulomb function. The arrays passed to this function
    !> need to be ordered from the inner-most (right-most) integral outward. The function
    !> result does not contain the overall phase and damp factor, which needs to be added
    !> manually.
    !>
    !> \param Z  Residual ion charge
    !> \param a  Lower bound
    !> \param c  Damping factor (additional exp(-c*r) added to all r^m functions)
    !> \param N  Dimension (number of integration variables)
    !> \param m  Array of integer powers of length N
    !> \param s  Array of integer signs (+1 or -1) in exponents of length 2*N
    !> \param k  Array of linear momenta of length 2*N
    !>
    complex(wp) function nested_exp_integ (Z, a, c, N, m, s, k) result (val)
 
        use multidip_params, only: ntermsppt
 
        integer,     intent(in) :: n, m(0:n-1), s(0:2*n-1)
        real(wp),    intent(in) :: z, a, c
        complex(wp), intent(in) :: k(0:2*n-1)
 
        integer, parameter :: npp = ntermsppt
 
        complex(wp) :: i(0:n*(npp+1)), f(0:n*(npp+1)-1), fg(0:n*(npp+1)-1), g(0:(n+1)*(npp+1)-1), j(0:n*(npp+1)-1)
        complex(wp) :: t(0:npp, 0:n*(npp + 1)-1)
 
        real(wp) :: b(0:n*(npp+1)-1, 0:n*(npp+1))
 
        complex(wp) :: u, v
        integer     :: o, p, q, lvl
 
        if (n < 1) then
            val = 0
            return
        end if
 
        ! For a given expansion depth of the integral (Npp) we need that many
        ! derivatives of the base function (here of the coordinate power r^m).
        ! With N levels of nesting, this combines to N * Npp derivatives in total
        ! for the first base function. Further base functions are not needed
        ! to such high orders.
 
        i  = 0
        f  = 0
        fg = 0
        g  = 0
        j  = 0
        t  = 0
 
        ! Construct Pascal triangle of binomial coefficients.
 
        b = 0
        do p = lbound(b, 1), ubound(b, 1)
            b(p, 0) = 1
            do q = 1, p - 1
                b(p, q) = b(p - 1, q - 1) + b(p - 1, q)
            end do
            b(p, p) = 1
        end do
 
        ! Constants of the phase function. This accumulates all phases along the nested dimensions.
 
        u = 0.
        v = 0.
 
        ! Initialize the "deepest-nested" integral, so that f*I = f.
 
        i(0) = 1.
 
        ! Recursively evaluate levels.
 
        do lvl = 0, n - 1
 
            ! 1. Evaluate the base function for this level, including all its derivatives.
            !    The base function is simply a coordinate power, r^m[lvl], so its
            !    derivatives can be calculated directly.
 
            f(0) = a**m(lvl)
 
            do p = 1, (n - lvl)*(npp + 1) - 1
                f(p) = (m(lvl) + 1 - p) * f(p - 1) / a
            end do
 
            ! 2. Combine the base function with the integral from the previous level,
            !    resulting in the effective base function fg for this level.
 
            do p = 0, (n - lvl)*(npp + 1) - 1
                fg(p) = 0
                do q = 0, p
                    fg(p) = fg(p) + b(p,q) * f(p - q) * i(q)
                end do
            end do
 
            ! 3. Evaluate the phase function at this level. This is required for construction
            !    of the individual terms of the asymptotic expansion, as well as for finalization
            !    of the integral after the terms are summer. The highest derivative degree
            !    needed is Npp plus the number of derivatives of the resulting integral that
            !    we are calculating at this level.
 
            u = u + s(2*lvl+1)*k(2*lvl+1) + s(2*lvl+0)*k(2*lvl+0) + imu*c
            v = v + z*s(2*lvl+1)/abs(k(2*lvl+1)) + z*s(2*lvl+0)/abs(k(2*lvl+0))
 
            g(0) = imu*a / (u*a + v)
            g(1) = imu*v / ((u*a + v)*(u*a + v))
 
            do p = 2, (n - lvl + 1)*(npp + 1) - 1
                g(p) = g(p - 1) * (-p) * u / (u*a + v)
            end do
 
            ! 4. Evaluate all terms of the asymptotic expansion of the integral, as well as their
            !    derivatives up to the order required for the integral itself, that is (N - lvl)*(Npp + 1).
 
            do p = 0, (n - lvl)*(npp + 1) - 1
                t(0, p) = fg(p)
            end do
 
            do o = 1, npp
                do p = 0, (n - lvl)*(npp + 1) - o - 1
                    t(o, p) = 0
                    do q = 0, p + 1
                        t(o, p) = t(o, p) + b(p + 1, q) * g(p + 1 - q) * t(o - 1, q)
                    end do
                end do
            end do
 
            ! 5. Sum the terms of the asymptotic expansion to get the reduced integral J (and its derivatives).
 
            do p = 0, (n - lvl - 1)*(npp + 1)
                j(p) = 0
                do o = 0, npp
                    j(p) = j(p) + t(o, p)
                end do
            end do
 
            ! 5. Evaluate the resulting integral (and derivatives) at this level.
 
            do p = 0, (n - lvl - 1)*(npp + 1)
                i(p) = 0
                do q = 0, p
                    i(p) = i(p) + b(p, q) * g(p - q) * j(q)
                end do
            end do
 
        end do
 
        val = i(0)
 
    end function nested_exp_integ
 
 
    !> \brief   Multi-dimensional triangular integral of Coulomb-Hankel functions and powers
    !> \author  J Benda
    !> \date    2020 - 2023
    !>
    !> Evaluate the nested many-dimensional triangular integral
    !> \f[
    !>    \int\limits_a^{+\infty} \dots \int\limits_{r_3}^{+\infty} H_4(r_2) r_2^{m_2} H_3(r_2)
    !>      \int\limits_{r_2}^{+\infty} H_2(r_1) r_1^{m_1} H_1(r_1) \mathrm{d}r_1 \mathrm{d}r_2 \mathrm{d}r_3 \dots
    !> \f]
    !> using the asymptotic form of Coulomb-Hankel functions. The arrays passed to this function
    !> need to be ordered from the inner-most (right-most) integral outward. For each term, use
    !> \ref nested_exp_integ.
    !>
    !> \param Z  Residual ion charge
    !> \param a  Lower bound
    !> \param c  Damping factor (additional exp(-c*r) added to all r^m functions)
    !> \param N  Dimension (number of integration variables)
    !> \param m  Array of integer powers of length N
    !> \param s  Array of integer signs (+1 or -1) in exponents of length 2*N
    !> \param l  Array of angular momenta of length 2*N
    !> \param k  Array of linear momenta of length 2*N
    !>
    complex(wp) function nested_coul_integ (Z, a, c, N, m, s, l, k) result (val)
 
        use multidip_params,  only: ntermsasy
        use multidip_special, only: cphase, kahan_add
 
        integer,     intent(in) :: n, m(1:n), s(1:2*n), l(1:2*n)
        real(wp),    intent(in) :: z, a, c
        complex(wp), intent(in) :: k(1:2*n)
 
        integer     :: nn(2*n), ms(n), p
        complex(wp) :: an(2*n), bn(2*n), factor(2*n)
 
        logical     :: at_end
        real(wp)    :: total_phase, total_exponent, kr, ki
        complex(wp) :: integ, err
 
        at_end = .false.
        val = 0
        err = 0
        factor = 1.
        nn = 0
 
        ! initial values of the Pochhammer symbols in the asymptotic series for the Coulomb-Hankel (or real Whittaker) functions
        do p = 1, 2*n
            an(p) = 1 + l(p) - z*s(p)*imu/k(p)
            bn(p) = -l(p) - z*s(p)*imu/k(p)
        end do
 
        ! sum all terms of the Cartesian product of all Coulomb-Hankel (or real Whittaker) functions expansions
        do while (.not. at_end)
 
            ! total power in r^m
            do p = 1, n
                ms(p) = m(p) - nn(2*p-1) - nn(2*p-0)
            end do
 
            ! radial integral
            integ = nested_exp_integ(z, a, c, n, ms, s, k)
 
            ! combine all coefficients from the expansions of all Coulomb-Hankel (or real Whittaker) functions
            do p = 1, 2*n
                integ = integ * factor(p)
            end do
 
            ! add the integral value using the floating-point-error compensating routine
            call kahan_add(val, integ, err)
 
            ! pick up the next single term from the product of expansions of all Coulomb-Hankel (or real Whittaker) functions
            do p = 1, 2*n
                nn(p) = nn(p) + 1
                if (nn(p) == ntermsasy) then
                    nn(p) = 0
                    factor(p) = 1.
                    at_end = .true.
                else
                    factor(p) = factor(p) * (an(p) + nn(p) - 1) * (bn(p) + nn(p) - 1) / (2 * imu * s(p) * nn(p) * k(p))
                    at_end = .false.
                    exit
                end if
            end do
 
        end do
 
        ! calculate the overall exponential and phase factor
        total_phase = 0
        total_exponent = 0
        do p = 1, 2*n
            kr = real(k(p), wp)
            ki = aimag(k(p))
            if (ki == 0) then
                ! positive energy solution -> Coulomb-Hankel function
                total_phase = total_phase + s(p)*(kr*a + z*log(2*kr*a)/kr - pi*l(p)/2 + cphase(z, l(p), kr))
            else
                ! negative energy solution -> Whittaker function
                total_exponent = total_exponent + s(p)*(-ki*a + z*log(2*ki*a)/ki)
            end if
        end do
 
        ! multiply the sum of expansion terms with the combined exponential factor
        val = val * cmplx(cos(total_phase), sin(total_phase), wp) * exp(total_exponent - n*c*a)
 
    end function nested_coul_integ
 
 
    !> \brief   Multi-dimensional integral of Coulomb-Hankel and Coulomb-Green functions and powers
    !> \author  J Benda
    !> \date    2020 - 2024
    !>
    !> Evaluate the many-dimensional integral
    !> \f[
    !>     \int\limits_a^{+\infty} \dots \int\limits_a^{+\infty} H_N(r_N) \dots g_3(r_3, r_2) r_2^{m_2}
    !>        g_2(r_2, r_1) r_1^{m_1} H_1(r_1) \mathrm{d}r_1 \dots
    !> \f]
    !> using the asymptotic form of Coulomb-Hankel functions. The arrays passed to this function
    !> need to be ordered from the right-most integral to the right. This function iterates over
    !> all possible orderings of \f$ (r_1, r_2, \dots, r_N) \f$ and for each of these it splits the
    !> integral into (hyper-)triangular integrals and integrates those using \ref nested_coul_integ.
    !>
    !> For small values of "a" it may increase the accuracy if the leading interval (a,r0) (or more specifically the set
    !> Q = (a,+∞)^N \ (b,+∞)^N) is integrated numerically. For such use, provide r0 > a. Otherwise set r0 to zero.
    !>
    !> When the global paramter `coulomb_check` is on, the integral is not attempted if any of the Coulomb functions
    !> in the integrand is not sufficiently well approximated by the asymptotic form at the R-matrix radius (or at the
    !> asymptotic radius if given). In such a case, the function returns the NaN constant.
    !>
    !> The one-dimensional case with c = 0, m = 0, l1 = l2 and real momenta is treated in a special way using a closed-form
    !> formula from the article: *H. F. Arnoldus, T. F. George: Analytical evaluation of elastic Coulomb integrals,
    !> J. Math. Phys. 33 (1992) 578–583*.
    !>
    !> The (adapted) formula reads:
    !> \f[
    !>     \lim_{c \rightarrow 0+} \int\limits_a^{+\infty} X_l(k_2; r) Y_l(k_1; r) \exp(-cr) dr
    !>     = \frac{X_l'(k_2; a) Y_l(k_1; a) - X_l(k_2; a) Y_l'(k_1; a)}{k_2^2 - k_1^2} \,,
    !> \f]
    !> where \f$ X_l \f$ and \f$ Y_l \f$ stand for arbitrary Coulomb functions and primes denote derivatives with respect
    !> to the radius.
    !>
    !> \param Z  Residual ion charge
    !> \param a  Lower bound of the integration
    !> \param r0 Optional radius from which to apply asymptotic integrals
    !> \param c  Damping factor (additional exp(-c*r) added to all r^m functions)
    !> \param N  Dimension (number of integration variables)
    !> \param sa Sign of the right-most Coulomb-Hankel function
    !> \param sb Sign of the left-most Coulomb-Hankel function
    !> \param m  Array of integer powers of length N
    !> \param l  Array of angular momenta of length N + 1
    !> \param k  Array of linear momenta of length N + 1
    !>
    complex(wp) function nested_cgreen_integ (Z, a, r0, c, N, sa, sb, m, l, k) result (val)
 
        use ieee_arithmetic,  only: ieee_value, ieee_quiet_nan
        use mpi_gbl,          only: mpi_xermsg
        use multidip_levin,   only: nested_cgreen_correct_levin
        use multidip_params,  only: cache_integrals, coultol, coulomb_check, num_integ_algo, print_warnings, ion_ion_analytic
        use multidip_romberg, only: nested_cgreen_correct_romberg, nested_cgreen_eval
        use multidip_special, only: next_permutation, kahan_add, coul, coulh, coulh_asy
 
        integer,     intent(in) :: n, sa, sb, m(:), l(:)
        real(wp),    intent(in) :: z, a, c, r0
        complex(wp), intent(in) :: k(:)
 
        logical     :: acc(n)
        integer     :: order(n), mp(n), lp(2*n), sp(2*n), p, signs
        real(wp)    :: b, rs(n), f(2, 2), g(2, 2)
        complex(wp) :: kp(2*n), err, integ, ix, ia, h(2, 2)
 
        character(len=100) :: msg
 
        val = 0
        err = 0
        b = a
 
        ! integrate the ion-ion transition analytically
        if (ion_ion_analytic .and. n == 1 .and. m(1) == 0 .and. l(1) == l(2) .and. aimag(k(1)) == 0 .and. aimag(k(2)) == 0) then
            call coul(z, l(1), real(k(1))**2/2, a, f(1,1), f(2,1), g(1,1), g(2,1))
            call coul(z, l(2), real(k(2))**2/2, a, f(1,2), f(2,2), g(1,2), g(2,2))
            h(:, 1) = cmplx(g(:, 1), sa*f(:, 1), wp)
            h(:, 2) = cmplx(g(:, 2), sb*f(:, 2), wp)
            val = (h(2,2)*h(1,1) - h(1,2)*h(2,1))/(real(k(2))**2 - real(k(1))**2)
            return
        end if
 
        ! test the accuracy of the integrand at the R-matrix radius "a" (NOTE: numerical correction implemented only for N = 1)
        if (a < r0 .and. n == 1) then
            rs = a
            ix = nested_cgreen_eval(z, c, n, sa, sb, m, l, k, rs, .false.)
            ia = nested_cgreen_eval(z, c, n, sa, sb, m, l, k, rs, .true.)
            if (abs(ix) > 1e-10 .and. abs(ia - ix) > coultol * abs(ix)) then
                ! R-matrix radius too small, prefer the user-given asymptotic radius "r0"
                b = r0
            end if
        end if
 
        ! test the accuracy of the integrand at the asymptotic radius "b" (NOTE: nested_cgreen_eval implemented only for N = 1)
        if (coulomb_check .and. n == 1) then
            rs = b
            ix = nested_cgreen_eval(z, c, n, sa, sb, m, l, k, rs, .false.)
            ia = nested_cgreen_eval(z, c, n, sa, sb, m, l, k, rs, .true.)
            if (abs(ix) > 1e-10 .and. abs(ia - ix) > abs(ix)) then
                if (print_warnings) then
                    ! asymptotic radius too small, do not bother to evaluate the integrals at all
                    !$omp critical (print_warning)
                    write (*, '(5x,A,3(I0,1x))', advance = 'no') 'Warning: Giving up Coul-Green integ '
                    write (*, '(3(A,I0,SP))', advance = 'no') 'm = ', m, ', sa = ', sa, ', sb = ', sb
                    write (*, '(A,*(I0,","))', advance = 'no') ', l = ', l
                    write (*, '(A,SP,*(F0.4,","))', advance = 'no') ' E = ', real(k**2/2, wp)
                    write (*, '(SP,2(A,2E10.4),A)') ' Ix = ', ix, 'i, Ia = ', ia, 'i'
                    !$omp end critical (print_warning)
                end if
                return
            end if
        end if
 
        ! initial ordering of coordinates, rN < ... < r2 < r1 (order(1) = index of the largest coordinate)
        order = -1
 
        ! for all orderings of the coordinates { rN, ..., r2, r1 }
        do while (next_permutation(order))
 
            ! reorder quantum numbers according to coordinate order
            do p = 1, n
                mp(p) = m(order(p))
            end do
 
            ! for all regular Coulomb function splittings, F = (H⁺ - H⁻)/2i
            do signs = 0, 2**(n - 1) - 1
 
                ! Each bit of 'signs' (starting from the least significant one) corresponds to the sign of the Coulomb-Hankel
                ! function H arising from the F function in the corresponding Coulomb-Green's function. If 'j' is the position
                ! of the bit (starting from 0 for the least significant one), then the sign belongs to the Coulomb-Green's
                ! function g(r(j+1), r(j)). There are N-1 Coulomb-Green's functions in total in the integrand, so we are iterating
                ! only over all combinations of the first N-1 bits.
 
                ! these (original) coordinate indices have been processed already
                acc = .false.
 
                ! assemble quantum numbers for coordinates sorted from largest to smallest
                do p = 1, n
 
                    ! For each coordinate, there are two Coulomb-Hankel functions. For r0, one of the Coulomb-Hankel functions
                    ! is H[sa]. For rN one of them is H[sb]. Others come from separation of g(r(j+1), r(j)). Here we are processing
                    ! coordinate r(j), where j = order[n]. This coordinate is connected by the Coulomb-Green's functions to the
                    ! previous coordinate r(j-1) and to the next coordinate r(j+1). The former connection is realized by (j-1)-th
                    ! Coulomb-Green's function, while the latter by the j-th one. We compare the current coordinate r(j) to its
                    ! neighbours r(j-1) and r(j+1) simply by checking whether we have already processed them, because we are
                    ! processing the coordinates in order of the descending order. Based on this comparison, either the appropriate
                    ! F sign (if this coordinate is smaller than neighbour) or the +1 sign for H+ (if the this coordinate is greater
                    ! then the neighbour) is used.
 
                    sp(2*p - 1) = sa
                    sp(2*p - 0) = sb
 
                    if (order(p) > 1) sp(2*p - 1) = merge(merge(-1, +1, btest(signs, order(p) - 2)), +1, acc(order(p) - 1))
                    if (order(p) < n) sp(2*p - 0) = merge(merge(-1, +1, btest(signs, order(p) - 1)), +1, acc(order(p) + 1))
 
                    acc(order(p)) = .true.
 
                    kp(2*p - 1) = k(order(p) - 0)
                    kp(2*p - 0) = k(order(p) + 1)
 
                    lp(2*p - 1) = l(order(p) - 0)
                    lp(2*p - 0) = l(order(p) + 1)
 
                end do
 
                ! evaluate integral for this hyper-triangle
                integ = (-1)**popcnt(signs) * nested_coul_integ(z, b, c, n, mp, sp, lp, kp)
 
                ! add the integral value using the floating-point-error compensating routine
                call kahan_add(val, integ, err)
 
            end do
 
        end do
 
        ! add the Green's function factors
        do p = 2, n
            val = val * imu / k(p)
        end do
 
        ! integrate over Q = (a,+∞)^N \ (b,+∞)^N numerically
        if (a < b) then
            select case (num_integ_algo)
                case (1)
                    call nested_cgreen_correct_romberg(z, a, b, c, n, sa, sb, m, l, k, val)
                case (2)
                    call nested_cgreen_correct_levin(z, a, b, c, n, sa, sb, m, l, k, val)
                case default
                    write (msg, '(a,i0)') 'Unknown numerical integration algorithm ', num_integ_algo
                    call mpi_xermsg('multidip_integ', 'nested_cgreen_integ', trim(msg), 1, 1)
            end select
        end if
 
    end function nested_cgreen_integ
 
end module multidip_integ