multidip_romberg.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_romberg
 
    use precisn_gbl, only: wp
 
    implicit none
 
contains
 
    subroutine nested_cgreen_correct_romberg (Z, a, b, c, N, sa, sb, m, l, k, v)
 
        use multidip_params,  only: max_romb_level, epsrel
        use multidip_special, only: kahan_add
 
        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
 
        integer  :: level, i
        real(wp) :: h, r
 
        complex(wp) :: romb(max_romb_level)  ! sequence of Romberg estimations, with the highest order placed in the first element
        complex(wp) :: integ, y, ya, yb, error, prev
 
        h = b - a
        ya = nested_cgreen_eval(z, c, n, sa, sb, m, l, k, [ a ], .false.)
        yb = nested_cgreen_eval(z, c, n, sa, sb, m, l, k, [ b ], .false.)
        romb(1) = h * (ya + yb) / 2
 
        do level = 2, max_romb_level
 
            ! remember the previous best Romberg estimate
            prev = romb(1)
 
            ! calculate next trapezoidal estimate from the previous one
            integ = 0
            error = 0
            h = h / 2
            do i = 1, 2**(level - 2)
                r = a + (2*i - 1)*h
                y = nested_cgreen_eval(z, c, n, sa, sb, m, l, k, [ r ], .false.)
                call kahan_add(integ, y, error)
            end do
            romb(level) = romb(level - 1)/2 + h*integ
 
            ! perform the Romberg extrapolation step
            do i = level - 1, 1, -1
                romb(i) = (4**(level - 1) * romb(i + 1) - romb(i)) / (4**(level - 1) - 1)
            end do
 
            ! abort on NaN
            if (.not. romb(1) == romb(1)) then
                print '(a)', 'WARNING: NaN encountered during Romberg quadrature'
                exit
            end if
 
            ! compare the current best estimate to the previous one
            if (abs(romb(1) - prev) <= epsrel * (abs(romb(1)) + abs(prev))) then
                exit
            end if
 
        end do
 
        if (level > max_romb_level) then
            print '(A,I0,A)', 'WARNING: Romberg quadrature did not converge in ', max_romb_level, ' subdivisions'
        end if
 
        v = v + romb(1)
 
    end subroutine nested_cgreen_correct_romberg
 
 
    complex(wp) function nested_cgreen_eval (Z, c, N, sa, sb, m, l, k, rs, asy) result (val)
 
        use mpi_gbl,          only: mpi_xermsg
        use multidip_special, only: coulh, coulh_asy
 
        integer,     intent(in) :: n, sa, sb, m(:), l(:)
        real(wp),    intent(in) :: z, c, rs(:)
        complex(wp), intent(in) :: k(:)
        logical,     intent(in) :: asy
 
        real(wp)    :: ek(n+1)
        complex(wp) :: h1, h2
 
        ek = real(k(1:n+1)**2, wp)/2
 
        ! WARNING: Only one-dimensional at the moment (no Green's functions)
        if (n > 1) then
            call mpi_xermsg('multidip_romberg', 'nested_cgreen_eval', 'nested_cgreen_eval not implemented for higher orders', 1, 1)
        end if
 
        if (asy) then
            h1 = coulh_asy(z, sa, l(1), ek(1), rs(1))
            h2 = coulh_asy(z, sb, l(2), ek(2), rs(1))
        else
            h1 = coulh(z, sa, l(1), ek(1), rs(1))
            h2 = coulh(z, sb, l(2), ek(2), rs(1))
        end if
 
        val = h2 * rs(1)**m(1) * exp(-c*rs(1)) * h1
 
    end function nested_cgreen_eval
 
end module multidip_romberg