multidip_levin.f90

Go to the documentation of this file.

! Copyright 2021
!
! 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_levin
 
    use precisn_gbl,    only: wp
 
    implicit none
 
contains
 
    subroutine nested_cgreen_correct_levin (Z, a, b, c, N, sa, sb, m, l, k, v)
 
        integer,     intent(in)    :: N, sa, sb, m(:), l(:)
        real(wp),    intent(in)    :: Z, a, b, c
        complex(wp), intent(in)    :: k(:)
        complex(wp), intent(inout) :: v
 
        ! WARNING: Implemented for 1D integrals only
        v = v + levin_adapt(z, a, b, c, m(1), sa, sb, l(1), l(2), k(1), k(2))
 
    end subroutine nested_cgreen_correct_levin
 
 
    complex(wp) function levin_adapt (Z, ra, rb, c, m, s1, s2, l1, l2, k1, k2) result (integ)
 
        use multidip_params, only: max_levin_level
 
        real(wp),    intent(in) :: z, ra, rb, c
        complex(wp), intent(in) :: k1, k2
        integer,     intent(in) :: m, s1, s2, l1, l2
 
        integer, parameter :: max_intervals = 2**max_levin_level
 
        integer     :: depth
        logical     :: converged(max_intervals)
        complex(wp) :: estimates(max_intervals)
        complex(wp) :: hl_buffer(2, 2, max_intervals + 1)
 
        ! at the beginning the integral over the whole integration range is definitely not converged
        converged(1) = .false.
 
        ! for all binary subdivision orders
        do depth = 0, max_levin_level - 1
 
            ! precalculate Coulomb function in subinterval endpoints for this subdivision of the integration interval
            call levin_prepare(z, ra, rb, s1, s2, l1, l2, k1, k2, depth, converged, hl_buffer)
 
            ! improve accuracy of not-yet-converged subintervals from previous subdivision
            call levin_improve(z, ra, rb, c, m, l1, l2, k1, k2, 0, 1, depth, converged, hl_buffer, estimates)
 
            ! abort on NaN
            if (.not. estimates(1) == estimates(1)) then
                print '(a)', 'WARNING: NaN encountered during Levin quadrature'
                exit
            end if
 
            ! if the top-level integral converged, terminate
            if (converged(1)) exit
 
        end do
 
        ! if full subdivision allowance was exhausted without convergence, complain
        if (depth == max_levin_level) then
            print '(A,I0,A)', 'WARNING: Adaptive Levin integration did not converge in ', max_levin_level, ' subdivisions'
        end if
 
        ! return the best stimate of the top-level integration interval
        integ = estimates(1)
 
    end function levin_adapt
 
 
    subroutine levin_prepare (Z, ra, rb, s1, s2, l1, l2, k1, k2, depth, converged, Hl_buffer)
 
        integer,     intent(in)    :: s1, s2, l1, l2, depth
        real(wp),    intent(in)    :: Z, ra, rb
        complex(wp), intent(in)    :: k1, k2
        complex(wp), intent(inout) :: Hl_buffer(2, 2, *)
        logical,     intent(inout) :: converged(:)
 
        real(wp) :: r
        integer  :: num_intervals, first_interval_idx, interval_idx, parent_interval_idx, evaluation_point_idx, storage_point_idx
 
        num_intervals = 2**depth
        first_interval_idx = 2**depth
 
        ! when called for the first time, evaluate Coulomb functions at the end points of the total integration range
        if (depth == 0) then
            call levin_eval(z, ra, s1, s2, l1, l2, k1, k2, hl_buffer(:, :, 2))
            call levin_eval(z, rb, s1, s2, l1, l2, k1, k2, hl_buffer(:, :, 1))
            return
        end if
 
        ! evaluate Coulomb functions in all other needed points at this subdivision level
        do interval_idx = first_interval_idx, first_interval_idx + num_intervals - 1, 2
 
            ! check if the parent interval (if any) of this sub-interval is converged
            if (depth > 0) then
 
                parent_interval_idx = interval_idx / 2
 
                ! if yes, mark this sub-interval as converged, too, and move on
                if (converged(parent_interval_idx)) then
                    converged(interval_idx) = .true.
                    converged(interval_idx + 1) = .true.
                    cycle
                end if
 
            end if
 
            ! otherwise mark it as unconverged
            converged(interval_idx) = .false.
            converged(interval_idx + 1) = .false.
 
            ! only odd index (new in this subdivision level) needs to be precomputed, which is the right endpoint of this interval
            evaluation_point_idx = interval_idx + 1
            storage_point_idx = (evaluation_point_idx + 3)/2
 
            ! calculate position
            r = ra + (evaluation_point_idx - first_interval_idx) * (rb - ra) / num_intervals
 
            ! perform the evaluation of Coulomb functions
            call levin_eval(z, r, s1, s2, l1, l2, k1, k2, hl_buffer(:, :, storage_point_idx))
 
        end do
 
    end subroutine levin_prepare
 
 
    subroutine levin_eval (Z, r, s1, s2, l1, l2, k1, k2, Hl)
 
        use multidip_special, only: coulh
 
        real(wp),    intent(in)    :: Z, r
        integer,     intent(in)    :: s1, s2, l1, l2
        complex(wp), intent(in)    :: k1, k2
        complex(wp), intent(inout) :: Hl(2, 2)
 
        real(wp) :: Ek1, Ek2
 
        ek1 = real(k1**2, wp)/2
        ek2 = real(k2**2, wp)/2
 
        hl(1, 1) = coulh(z, s1, l1,     ek1, r)
        hl(2, 1) = coulh(z, s1, l1 + 1, ek1, r)
        hl(1, 2) = coulh(z, s2, l2,     ek2, r)
        hl(2, 2) = coulh(z, s2, l2 + 1, ek2, r)
 
    end subroutine levin_eval
 
 
    recursive subroutine levin_improve (Z, ra, rb, c, m, l1, l2, k1, k2, depth, interval_idx, tgt_depth, converged, Hl, estimates)
 
        real(wp),    intent(in)    :: z, ra, rb, c
        complex(wp), intent(in)    :: k1, k2
        integer,     intent(in)    :: m, l1, l2, depth, interval_idx, tgt_depth
        complex(wp), intent(in)    :: hl(2, 2, *)
        logical,     intent(inout) :: converged(:)
        complex(wp), intent(inout) :: estimates(:)
 
        integer  :: num_local_intervals, first_interval_idx, max_global_point_idx, left_point_global_idx, right_point_global_idx
        integer  :: left_point_local_idx, right_point_local_idx, left_point_unique_idx, right_point_unique_idx
        integer  :: left_point_storage_idx, right_point_storage_idx, left_subinterval_idx, right_subinterval_idx
        real(wp) :: a, b, epsrel, epsabs, delta, reference
 
        complex(wp) :: best_estimate,improved_estimate
 
        ! at the target depth simply evaluate the integrals over the given subinterval and store the estimates
        if (depth == tgt_depth) then
 
            num_local_intervals = 2**depth
            first_interval_idx = 2**depth
            max_global_point_idx = 2**(depth + 1)
 
            left_point_global_idx = interval_idx;
            right_point_global_idx = left_point_global_idx + 1
 
            left_point_local_idx = interval_idx - first_interval_idx
            right_point_local_idx = left_point_local_idx + 1
 
            left_point_unique_idx = shiftr(left_point_global_idx, trailz(left_point_global_idx))
            right_point_unique_idx = shiftr(mod(right_point_global_idx, max_global_point_idx), trailz(right_point_global_idx))
 
            left_point_storage_idx = (left_point_unique_idx + 3)/2
            right_point_storage_idx = merge(1, (right_point_unique_idx + 3)/2, right_point_unique_idx == 0)
 
            a = ra + left_point_local_idx * (rb - ra) / num_local_intervals
            b = ra + right_point_local_idx * (rb - ra) / num_local_intervals
 
            estimates(interval_idx) = levin_integrate(z, a, b, c, m, l1, l2, k1, k2, &
                                                      hl(:, :, left_point_storage_idx), &
                                                      hl(:, :, right_point_storage_idx))
 
            return
 
        end if
 
        ! get subinterval indices
        left_subinterval_idx = 2*interval_idx;
        right_subinterval_idx = 2*interval_idx + 1;
 
        ! update all non-converged sub-intervals
        if (.not. converged(left_subinterval_idx)) then
            call levin_improve(z, ra, rb, c, m, l1, l2, k1, k2, &
                               depth + 1, left_subinterval_idx, tgt_depth, converged, hl, estimates)
        end if
        if (.not. converged(right_subinterval_idx)) then
            call levin_improve(z, ra, rb, c, m, l1, l2, k1, k2, &
                               depth + 1, right_subinterval_idx, tgt_depth, converged, hl, estimates)
        end if
 
        ! check convergence by comparing estimate for this interval to the sum of the estimates for the two sub.intervals
        best_estimate = estimates(interval_idx)
        improved_estimate = estimates(left_subinterval_idx) + estimates(right_subinterval_idx)
        epsrel = 1e-7
        epsabs = epsrel * abs(estimates(1))  ! TODO ... take into account there are more unconverged subintervals in this depth
        delta = abs(best_estimate - improved_estimate)
        reference = 0.5 * (abs(best_estimate) + abs(improved_estimate))
        if (delta < epsabs * reference .or. delta < epsabs) then
            converged(interval_idx) = .true.
            converged(left_subinterval_idx) = .true.
            converged(right_subinterval_idx) = .true.
        end if
 
        ! update the estimate for this interval with the sum of estimates for the two subintervals
        estimates(interval_idx) = improved_estimate
 
    end subroutine levin_improve
 
 
    complex(wp) function levin_integrate (Z, ra, rb, c, m, l1, l2, k1, k2, Hla, Hlb) result (integ)
 
        use mpi_gbl, only: mpi_xermsg
 
        real(wp),    intent(in)    :: z, ra, rb, c
        complex(wp), intent(in)    :: k1, k2
        integer,     intent(in)    :: m, l1, l2
        complex(wp), intent(in)    :: hla(2, 2), hlb(2, 2)
 
        if (aimag(k1) == 0 .and. aimag(k2) == 0) then
            ! free-free transition, both Coulomb functions have positive energy and oscillate -> use 4-dimensional Levin integrator
            integ = levin_integrate_4x4(z, ra, rb, c, m, l1, l2, real(k1, wp), real(k2, wp), hla, hlb)
        else if (aimag(k1) == 0 .and. aimag(k2) /= 0) then
            ! free-closed transition, only H1 has positive energy and oscillates -> use 2-dimensional Levin integrator
            integ = levin_integrate_2x2(z, ra, rb, c, m, l2, l1, aimag(k2), real(k1, wp), hla(:, 1), hlb(:, 1))
        else if (aimag(k1) /= 0 .and. aimag(k2) == 0) then
            ! closed-free transition, only H2 has positive energy and oscillates -> use 2-dimensional Levin integrator
            integ = levin_integrate_2x2(z, ra, rb, c, m, l1, l2, aimag(k1), real(k2, wp), hla(:, 2), hlb(:, 2))
        else
            call mpi_xermsg('multidip_levin', 'levin_integrate', 'Levin integrator cannot integrate closed-closed integrals.', 1, 1)
        end if
 
    end function levin_integrate
 
 
    complex(wp) function levin_integrate_2x2 (Z, ra, rb, c, m, lc, lo, kc, ko, wa, wb) result (integ)
 
        use blas_lapack_gbl,  only: blasint
        use mpi_gbl,          only: mpi_xermsg
        use multidip_params,  only: rone, rhalf, cheb_order
        use multidip_special, only: blas_dgetrf => dgetrf, blas_dgetrs => dgetrs, coulh
        use phys_const_gbl,   only: pi
 
        real(wp),    intent(in)    :: z, ra, rb, c, kc, ko
        integer,     intent(in)    :: m, lc, lo
        complex(wp), intent(in)    :: wa(2), wb(2)
 
        integer, parameter :: levin_dim = 2
        integer, parameter :: mat_dim = levin_dim * (cheb_order + 1)
 
        integer  :: ipoly, inode, a, b, i, j, row, col, sgn, ifun, flips
 
        real(wp) :: x, r, rl, elem, cx
        real(wp) :: cheb_node(0 : cheb_order)
        real(wp) :: cheb_value(0 : cheb_order, 0 : cheb_order)
        real(wp) :: cheb_deriv(0 : cheb_order, 0 : cheb_order)
        real(wp) :: aa(levin_dim, levin_dim)
        real(wp) :: ab(levin_dim, levin_dim)
        real(wp) :: mat(mat_dim, mat_dim)
        real(wp) :: coef(mat_dim)
        real(wp) :: pa(levin_dim), pb(levin_dim)
 
        complex(wp) :: hlna, hlnb, hlpa, hlpb, wl
        complex(wp) :: contrib_a, contrib_b, delta
 
        integer(blasint) :: n, pivots(mat_dim), info
 
        ! only allow positive non-empty intervals
        if (ra >= rb) then
            call mpi_xermsg('multidip_levin', 'levin_integrate_2x2', &
                            'Levin integration not implemented for non-positive intervals', 1, 1)
        end if
 
        ! populate split matrices of the recurrence relations (total matrix is LevinA = Aa/r + Ab)
        rl = sqrt(ko*ko + z*z/((lo + 1)*(lo + 1)))
 
        aa(1, 1) = + lo + 1;  aa(1, 2) = 0
        aa(2, 1) = 0;         aa(2, 2) = - lo - 1
 
        ab(1, 1) = -z/(lo + 1);     ab(1, 2) = -rl
        ab(2, 1) = +rl;             ab(2, 2) = +z/(lo + 1)
 
        ! evaluate Chebyshev nodes as well as Chebyshev polynomials and their derivatives in these nodes
        do inode = 0, cheb_order
 
            x = cos((inode + rhalf) * pi / (cheb_order + 1))
 
            cheb_node(inode) = x
 
            cheb_value(inode, 0) = 1;   cheb_deriv(inode, 0) = 0
            cheb_value(inode, 1) = x;   cheb_deriv(inode, 1) = 1
 
            do ipoly = 1, cheb_order - 1
                cheb_value(inode, ipoly + 1) = 2*x*cheb_value(inode, ipoly) - cheb_value(inode, ipoly - 1)
                cheb_deriv(inode, ipoly + 1) = 2*(ipoly + 1)*cheb_value(inode, ipoly) &
                                                + (ipoly + 1)*cheb_deriv(inode, ipoly - 1) / max(1, ipoly - 1)
            end do
 
        end do
 
        ! construct the matrix of the equations
        do a = 1, levin_dim
            do i = 0, cheb_order
                row = (a - 1) * (cheb_order + 1) + i + 1
                do b = 1, levin_dim
                    do j = 0, cheb_order
                        col = (b - 1) * (cheb_order + 1) + j + 1
                        x = cheb_node(i)
                        r = ra + (x + 1) / 2 * (rb - ra)
                        elem = aa(b, a)/r + ab(b, a)
                        mat(row, col) = elem * cheb_value(i, j)
                        if (a == b) then
                            mat(row, col) = mat(row, col) + 2/(rb - ra) * cheb_deriv(i, j)
                        end if
                    end do
                end do
            end do
        end do
 
        ! construct the right-hand side (only populate elements pertaining to the first function)
        do a = 1, levin_dim
            do i = 0, cheb_order
                row = (a - 1) * (cheb_order + 1) + i + 1
                x = cheb_node(i)
                r = ra + (x + 1) * (rb - ra) / 2
                if (a == 1) then
                    wl = real(coulh(z, 0, lc, -kc*kc/2, r), wp)
                    coef(row) = r**m * exp(-c*r) * wl
                else
                    coef(row) = 0
                end if
            end do
        end do
 
        ! solve the set of equations
        n = mat_dim
        call blas_dgetrf(n, n, mat, n, pivots, info)
        if (info /= 0) then
            call mpi_xermsg('multidip_levin', 'levin_integrate_2x2', &
                            'LU decomposition in Levin integration failed', int(info), 1)
        end if
        call blas_dgetrs('N', n, 1_blasint, mat, n, pivots, coef, n, info)
        if (info /= 0) then
            call mpi_xermsg('multidip_levin', 'levin_integrate_2x2', &
                            'LU backsubstitution in Levin integration failed', int(info), 1)
        end if
 
        ! evaluate the boundary terms
        do a = 1, levin_dim
            pa(a) = 0
            pb(a) = 0
            do i = 0, cheb_order
                row = (a - 1) * (cheb_order + 1) + i + 1
                cx = coef(row)
                sgn = merge(+1, -1, mod(i, 2) == 0)
                pa(a) = pa(a) + sgn*cx
                pb(a) = pb(a) + cx
            end do
        end do
 
        ! evaluate the integral
        integ = 0
        do ifun = 1, levin_dim
            contrib_a = wa(ifun) * pa(ifun)
            contrib_b = wb(ifun) * pb(ifun)
            delta = contrib_b - contrib_a
            integ = integ + delta
        end do
 
    end function levin_integrate_2x2
 
 
    complex(wp) function levin_integrate_4x4 (Z, ra, rb, c, m, l1, l2, k1, k2, Hla, Hlb) result (integ)
 
        use blas_lapack_gbl,  only: blasint
        use mpi_gbl,          only: mpi_xermsg
        use multidip_params,  only: rone, rhalf, cheb_order
        use multidip_special, only: blas_dgetrf => dgetrf, blas_dgetrs => dgetrs
        use phys_const_gbl,   only: pi
 
        real(wp),    intent(in)    :: z, ra, rb, c, k1, k2
        integer,     intent(in)    :: m, l1, l2
        complex(wp), intent(in)    :: hla(2, 2), hlb(2, 2)
 
        integer, parameter :: levin_dim = 4
        integer, parameter :: mat_dim = levin_dim * (cheb_order + 1)
 
        integer  :: ipoly, inode, a, b, i, j, row, col, sgn, ifun, flips
 
        real(wp) :: x, r, r1, r2, elem, cx
        real(wp) :: cheb_node(0 : cheb_order)
        real(wp) :: cheb_value(0 : cheb_order, 0 : cheb_order)
        real(wp) :: cheb_deriv(0 : cheb_order, 0 : cheb_order)
        real(wp) :: aa(levin_dim, levin_dim)
        real(wp) :: ab(levin_dim, levin_dim)
        real(wp) :: mat(mat_dim, mat_dim)
        real(wp) :: coef(mat_dim)
        real(wp) :: pa(levin_dim), pb(levin_dim)
 
        complex(wp) :: hln1a, hln2a, hln1b, hln2b, hlp1a, hlp2a, hlp1b, hlp2b, wa(levin_dim), wb(levin_dim)
        complex(wp) :: contrib_a, contrib_b, delta
 
        integer(blasint) :: n, pivots(mat_dim), info
 
        ! only allow positive non-empty intervals
        if (ra >= rb) then
            call mpi_xermsg('multidip_levin', 'levin_integrate_4x4', &
                            'Levin integration not implemented for non-positive intervals', 1, 1)
        end if
 
        ! populate split matrices of the recurrence relations (total matrix is LevinA = Aa/r + Ab)
        aa = 0
        ab = 0
 
        r1 = sqrt(k1*k1 + z/((l1 + 1)*(l1 + 1)))
        r2 = sqrt(k2*k2 + z/((l2 + 1)*(l2 + 1)))
 
        aa(1, 1) = + l1 + l2 + 2
        aa(2, 2) = + l1 - l2
        aa(3, 3) = - l1 + l2
        aa(4, 4) = - l1 - l2 - 2
 
        ab(1, 1) = - z/(l1 + 1) - z/(l2 + 1)
        ab(2, 2) = - z/(l1 + 1) + z/(l2 + 1)
        ab(3, 3) = + z/(l1 + 1) - z/(l2 + 1)
        ab(4, 4) = + z/(l1 + 1) + z/(l2 + 1)
 
        ab(1, 2) = -r2;  ab(1, 3) = -r1
        ab(2, 1) = +r2;  ab(2, 4) = -r1
        ab(3, 1) = +r1;  ab(3, 4) = -r2
        ab(4, 2) = +r1;  ab(4, 3) = +r2
 
        ! evaluate Chebyshev nodes as well as Chebyshev polynomials and their derivatives in these nodes
        do inode = 0, cheb_order
 
            x = cos((inode + rhalf) * pi / (cheb_order + 1))
 
            cheb_node(inode) = x
 
            cheb_value(inode, 0) = 1;   cheb_deriv(inode, 0) = 0
            cheb_value(inode, 1) = x;   cheb_deriv(inode, 1) = 1
 
            do ipoly = 1, cheb_order - 1
                cheb_value(inode, ipoly + 1) = 2*x*cheb_value(inode, ipoly) - cheb_value(inode, ipoly - 1)
                cheb_deriv(inode, ipoly + 1) = 2*(ipoly + 1)*cheb_value(inode, ipoly) &
                                                + (ipoly + 1)*cheb_deriv(inode, ipoly - 1) / max(1, ipoly - 1)
            end do
 
        end do
 
        ! construct the matrix of the equations
        do a = 1, levin_dim
            do i = 0, cheb_order
                row = (a - 1) * (cheb_order + 1) + i + 1
                do b = 1, levin_dim
                    do j = 0, cheb_order
                        col = (b - 1) * (cheb_order + 1) + j + 1
                        x = cheb_node(i)
                        r = ra + (x + 1) / 2 * (rb - ra)
                        elem = aa(b, a)/r + ab(b, a)
                        mat(row, col) = elem * cheb_value(i, j)
                        if (a == b) then
                            mat(row, col) = mat(row, col) + 2/(rb - ra) * cheb_deriv(i, j)
                        end if
                    end do
                end do
                !print '(8x,SP,*(E10.3,1x))', real(mat(row, :), wp)
            end do
        end do
 
        ! construct the right-hand side (only populate elements pertaining to the first function)
        do a = 1, levin_dim
            do i = 0, cheb_order
                row = (a - 1) * (cheb_order + 1) + i + 1
                x = cheb_node(i)
                r = ra + (x + 1) * (rb - ra) / 2
                if (a == 1) then
                    coef(row) = r**m * exp(-c*r)
                else
                    coef(row) = 0
                end if
            end do
        end do
 
        ! solve the set of equations
        n = mat_dim
        call blas_dgetrf(n, n, mat, n, pivots, info)
        if (info /= 0) then
            call mpi_xermsg('multidip_levin', 'levin_integrate_4x4', &
                            'LU decomposition in Levin integration failed', int(info), 1)
        end if
        call blas_dgetrs('N', n, 1_blasint, mat, n, pivots, coef, n, info)
        if (info /= 0) then
            call mpi_xermsg('multidip_levin', 'levin_integrate_4x4', &
                            'LU backsubstitution in Levin integration failed', int(info), 1)
        end if
 
        hln1a = hla(1, 1); hlp1a = hla(2, 1)
        hln2a = hla(1, 2); hlp2a = hla(2, 2)
        hln1b = hlb(1, 1); hlp1b = hlb(2, 1)
        hln2b = hlb(1, 2); hlp2b = hlb(2, 2)
 
        wa = [ hln1a * hln2a, hln1a * hlp2a, hlp1a * hln2a, hlp1a * hlp2a ]
        wb = [ hln1b * hln2b, hln1b * hlp2b, hlp1b * hln2b, hlp1b * hlp2b ]
 
        ! evaluate the boundary terms
        do a = 1, levin_dim
            pa(a) = 0
            pb(a) = 0
            do i = 0, cheb_order
                row = (a - 1) * (cheb_order + 1) + i + 1
                cx = coef(row)
                sgn = merge(+1, -1, mod(i, 2) == 0)
                pa(a) = pa(a) + sgn*cx
                pb(a) = pb(a) + cx
            end do
        end do
 
        ! evaluate the integral
        integ = 0
        do ifun = 1, levin_dim
            contrib_a = wa(ifun) * pa(ifun)
            contrib_b = wb(ifun) * pb(ifun)
            delta = contrib_b - contrib_a
            integ = integ + delta
        end do
 
    end function levin_integrate_4x4
 
end module multidip_levin