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/>.
!
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
 
    type :: nested_cgreen_integ_cache_internal_t
 
        integer,     allocatable :: l(:, :)         ! known sets of angular momentum
        integer,     allocatable :: m(:, :)         ! known sets of angular momentum projection
        complex(wp), allocatable :: k(:, :)         ! per stored integral: sequence of momenta
        integer,     allocatable :: idx(:, :)       ! per stored integral: 1. 'sa', 2. 'sb', 3. 'm' idx, 4. 'l' idx, 5. 'k' idx
        complex(wp), allocatable :: values(:, :)    ! per stored integral: its value (first dimension is always 1)
 
    end type
 
    type :: nested_cgreen_integ_cache_t
 
        real(wp) :: a       ! R-matrix radius
        real(wp) :: r0      ! asymptotic radius
        real(wp) :: c       ! dipole damping coefficient
 
        integer :: ie       ! which energy to use
 
        type(nested_cgreen_integ_cache_internal_t), allocatable :: data(:, :)  ! one data structure per photon order and energy
 
    contains
 
        procedure :: init => init_nested_cgreen_integ_cache_t
        procedure :: get  => get_nested_cgreen_integ_cache_t
        procedure :: set  => set_nested_cgreen_integ_cache_t
 
    end type
 
    type(nested_cgreen_integ_cache_t) :: nested_cgreen_integ_cache
 
contains
 
    subroutine init_nested_cgreen_integ_cache_t (this, N, ne, a, r0, c)
 
        class(nested_cgreen_integ_cache_t), intent(inout) :: this
        integer,                            intent(in)    :: N, ne
        real(wp),                           intent(in)    :: a, r0, c
 
        if (allocated(this % data)) deallocate(this % data)
 
        allocate (this % data(n, ne))
 
        this % a  = a
        this % r0 = r0
        this % c  = c
 
    end subroutine init_nested_cgreen_integ_cache_t
 
 
    logical function get_nested_cgreen_integ_cache_t (this, a, r0, c, N, sa, sb, m, l, k, val) result (status)
 
        class(nested_cgreen_integ_cache_t), intent(in) :: this
 
        real(wp),    intent(in)    :: a, r0, c
        integer,     intent(in)    :: n, sa, sb, m(:), l(:)
        complex(wp), intent(in)    :: k(:)
        complex(wp), intent(inout) :: val
 
        integer :: i, im, il, ik
 
        status = .false.
        if (allocated(this % data)) then
            if (this % a == a .and. this % r0 == r0 .and. this % c == c .and. size(this % data) >= n) then
                if (allocated(this % data(n, this % ie) % values)) then
                    im = find_column(this % data(n, this % ie) % m, m);  if (im == 0) return
                    il = find_column(this % data(n, this % ie) % l, l);  if (il == 0) return
                    ik = find_column(this % data(n, this % ie) % k, k);  if (ik == 0) return
                    do i = 1, size(this % data(n, this % ie) % idx, 2)
                        if (all(this % data(n, this % ie) % idx(:, i) == [sa, sb, im, il, ik])) then
                            val = this % data(n, this % ie) % values(1, i)
                            status = .true.
                            return
                        end if
                    end do
                end if
            end if
        end if
 
    end function get_nested_cgreen_integ_cache_t
 
 
    subroutine set_nested_cgreen_integ_cache_t (this, a, r0, c, N, sa, sb, m, l, k, val)
 
        class(nested_cgreen_integ_cache_t), intent(inout) :: this
 
        real(wp),    intent(in) :: a, r0, c
        integer,     intent(in) :: N, sa, sb, m(:), l(:)
        complex(wp), intent(in) :: k(:)
        complex(wp), intent(in) :: val
 
        integer     :: im, il, ik, ii
        complex(wp) :: orig_val
 
        if (this % get(a, r0, c, n, sa, sb, m, l, k, orig_val)) then
            write (error_unit, '(a)') 'Error: Overwriting existing integral cache entry'
            stop 1
        end if
 
        if (this % a /= a .or. this % r0 /= r0 .or. this % c /= c) then
            write(error_unit, '(a)') 'Error: Inserting value into cache for incompatible setup'
            stop 1
        end if
 
        im = find_column(this % data(n, this % ie) % m, m);  if (im == 0) im = append_column(this % data(n, this % ie) % m, m)
        il = find_column(this % data(n, this % ie) % l, l);  if (il == 0) il = append_column(this % data(n, this % ie) % l, l)
        ik = find_column(this % data(n, this % ie) % k, k);  if (ik == 0) ik = append_column(this % data(n, this % ie) % k, k)
 
        ii = append_column(this % data(n, this % ie) % idx,    [sa, sb, im, il, ik])
        ii = append_column(this % data(n, this % ie) % values, [val])
 
    end subroutine set_nested_cgreen_integ_cache_t
 
 
    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
 
 
    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
 
 
    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
        use multidip_romberg, only: nested_cgreen_correct_romberg, nested_cgreen_eval
        use multidip_special, only: next_permutation, kahan_add, 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)
        complex(wp) :: kp(2*n), err, integ, ix, ia
 
        character(len=100) :: msg
 
        val = 0
        err = 0
        b = a
 
        ! check if we computed this integral before
        if (nested_cgreen_integ_cache % get(a, r0, c, n, sa, sb, m, l, k, val)) return
 
        ! 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
 
        ! remember this integral
        if (cache_integrals) call nested_cgreen_integ_cache % set(a, r0, c, n, sa, sb, m, l, k, val)
 
    end function nested_cgreen_integ
 
end module multidip_integ