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 phys_const_gbl, only: pi, imu
    use precisn_gbl,    only: wp

    implicit none

contains

    complex(wp) function nested_exp_integ (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) :: 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 + s(2*lvl+1)/abs(k(2*lvl+1)) + 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 (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) :: 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) - s(p)*imu/k(p)
            bn(p) = -l(p) - 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(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 + log(2*kr*a)/kr - pi*l(p)/2 + cphase(l(p), kr))
            else
                ! negative energy solution -> Whittaker function
                total_exponent = total_exponent + s(p)*(-ki*a + 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 (a, c, N, sa, sb, m, l, k) result (val)

        use multidip_special, only: next_permutation, kahan_add

        integer,     intent(in) :: n, sa, sb, m(:), l(:)
        real(wp),    intent(in) :: a, c
        complex(wp), intent(in) :: k(:)

        logical     :: acc(n)
        integer     :: order(n), mp(n), lp(2*n), sp(2*n), p, signs
        complex(wp) :: kp(2*n), err, integ

        val = 0
        err = 0

        ! 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(a, 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

    end function nested_cgreen_integ

end module multidip_integ