multidip_special.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_special

    use, intrinsic :: iso_c_binding,   only: c_double, c_int, c_f_pointer
    use, intrinsic :: iso_fortran_env, only: input_unit, int32, int64, real64, output_unit

    use blas_lapack_gbl,  only: blasint
    use phys_const_gbl,   only: pi, imu
    use precisn_gbl,      only: wp

    implicit none

    real(wp), parameter :: rzero = 0
    real(wp), parameter :: rone = 1

#ifdef WITH_GSL

    type, bind(c) :: gsl_sf_result
        real(c_double) :: val
        real(c_double) :: err
    end type gsl_sf_result

    interface
        ! logarithm of the complex gamma function from GSL
        function gsl_sf_lngamma_complex_e (zr, zi, lnr, arg) bind(C, name='gsl_sf_lngamma_complex_e')
            use, intrinsic :: iso_c_binding, only: c_int, c_double
            import gsl_sf_result
            real(c_double), value :: zi, zr
            type(gsl_sf_result)   :: lnr, arg
            integer(c_int)        :: gsl_sf_lngamma_complex_e
        end function gsl_sf_lngamma_complex_e
    end interface

    interface
        ! Coulomb wave function from GSL
        function gsl_sf_coulomb_wave_fg_e (eta, rho, l, k, F, Fp, G, Gp, expF, expG) bind(C, name='gsl_sf_coulomb_wave_FG_e')
            use, intrinsic :: iso_c_binding, only: c_int, c_double
            import gsl_sf_result
            real(c_double), value :: eta, rho, l
            integer(c_int), value :: k
            type(gsl_sf_result)   :: f, fp, g, gp
            real(c_double)        :: expf, expg
            integer(c_int)        :: gsl_sf_coulomb_wave_fg_e
        end function gsl_sf_coulomb_wave_fg_e
    end interface

    interface
        ! irregular confluent hypergeometric function U
        pure function gsl_sf_hyperg_u (a, b, x) bind(C, name='gsl_sf_hyperg_U')
            use, intrinsic :: iso_c_binding, only: c_double
            real(c_double), value :: a, b, x
            real(c_double)        :: gsl_sf_hyperg_u
        end function gsl_sf_hyperg_u
    end interface
#endif

    interface
        ! Coulomb phase from UKRmol-out
        function cphaz (l, eta, iwrite)
            import wp
            integer  :: l, iwrite
            real(wp) :: eta, cphaz
        end function cphaz
    end interface

    interface
        ! Coulomb wave function from UKRmol-out
        subroutine coulfg (xx, eta1, xlmin, xlmax, fc, gc, fcp, gcp, mode1, kfn, ifail)
            import wp
            integer  :: mode1, kfn, ifail
            real(wp) :: xx, eta1, xlmin, xlmax, fc(*), gc(*), fcp(*), gcp(*)
        end subroutine coulfg
    end interface

    interface
        ! Decaying Whittaker function
        subroutine couln (l, Z, ERy, r, U, Up, acc, efx, ierr, nlast, fnorm)
            import wp
            integer  :: l, ierr, nlast
            real(wp) :: Z, Ery, r, U, Up, acc, efx, fnorm
        end subroutine couln
    end interface

    interface
        !Exponentially decaying solution with given (negative) energy and angular momentum
        subroutine decay (k, l, r, U, Up, iwrite)
            import wp
            integer  :: l, iwrite
            real(wp) :: k, r, U, Up
        end subroutine decay
    end interface

    interface
        ! real matrix-matrix multiplication from BLAS
        subroutine dgemm (transa, transb, m, n, k, alpha, A, lda, B, ldb, beta, C, ldc)
            import real64, blasint
            character(len=1) :: transa, transb
            integer(blasint) :: m, n, k, lda, ldb, ldc
            real(real64)     :: A(lda, *), B(ldb, *), C(ldc, *), alpha, beta
        end subroutine dgemm
    end interface

    interface
        ! real matrix-vector multiplication from BLAS
        subroutine dgemv (trans, m, n, alpha, A, lda, X, incx, beta, Y, incy)
            import real64, blasint
            character(len=1) :: trans
            integer(blasint) :: m, n, incx, incy, lda
            real(real64)     :: A(lda, *), X(*), Y(*), alpha, beta
        end subroutine dgemv
    end interface

    interface
        ! complex matrix-vector multiplication from BLAS
        subroutine zgemv (trans, m, n, alpha, A, lda, X, incx, beta, Y, incy)
            import real64, blasint
            character(len=1) :: trans
            integer(blasint) :: m, n, incx, incy, lda
            complex(real64)  :: A(lda, *), X(*), Y(*), alpha, beta
        end subroutine zgemv
    end interface

    interface
        ! complex LU decomposition from LAPACK
        subroutine zgetrf (m, n, A, lda, ipiv, info)
            import real64, blasint
            integer(blasint) :: m, n, lda, ipiv(*), info
            complex(real64)  :: A(lda, *)
        end subroutine zgetrf
    end interface

    interface
        ! complex matrix solution from LAPACK
        subroutine zgetrs (trans, n, nrhs, A, lda, ipiv, B, ldb, info)
            import real64, blasint
            character(len=1) :: trans
            integer(blasint) :: n, nrhs, lda, ldb, info, ipiv(*)
            complex(real64)  :: A(lda, *), B(ldb, *)
        end subroutine zgetrs
    end interface

    interface
        ! complex matrix inversion from LAPACK
        subroutine zgetri (n, A, lda, ipiv, work, lwork, info)
            import real64, blasint
            integer(blasint) :: n, lda, ipiv(*), lwork, info
            complex(real64)  :: A(lda, *), work(*)
        end subroutine zgetri
    end interface

contains

    subroutine calculate_s_matrix (Km, Sm, nopen)

        complex(wp), allocatable, intent(inout) :: Sm(:, :)
        real(wp),                 intent(in)    :: Km(:, :)
        integer,                  intent(in)    :: nopen

        complex(wp), allocatable :: IpiK(:, :), ImiK(:, :)
        integer :: i, nchan

        nchan = size(sm, 1)

        sm = 0

        allocate (ipik(nchan, nchan), imik(nchan, nchan))

        ipik = 0; forall (i = 1:nchan) ipik(i, i) = 1
        imik = 0; forall (i = 1:nchan) imik(i, i) = 1
        ipik(1:nopen, 1:nopen) = ipik(1:nopen, 1:nopen) + imu*km(1:nopen, 1:nopen)
        imik(1:nopen, 1:nopen) = imik(1:nopen, 1:nopen) - imu*km(1:nopen, 1:nopen)

        call invert_matrix(imik)

        sm = matmul(ipik, imik)

    end subroutine calculate_s_matrix


    subroutine scatak (rmatr, eig, etot, w, Ek, l, Km, Ak)

        complex(wp), allocatable, intent(inout) :: Ak(:,:)
        real(wp), intent(in) :: rmatr, eig(:), etot, Ek(:), w(:,:), Km(:,:)
        integer,  intent(in) :: l(:)

        complex(wp), allocatable :: T(:,:)
        real(wp),    allocatable :: V(:,:), TT(:,:,:), Sp(:,:), Cp(:,:), P(:), XX(:,:,:), wFp(:,:,:), k(:)

        integer(blasint) :: i, j, nchan, nstat, nchan2, nopen, nopen2
        real(wp)         :: F, Fp, G, Gp

        nchan = size(ek)
        nchan2 = 2*nchan
        nopen = count(ek > 0)
        nopen2 = 2*nopen
        nstat = size(eig)
        k = sqrt(2*abs(ek))

        allocate (v(nchan, nchan), t(nopen, nopen), sp(nchan, nchan), cp(nchan, nchan), tt(nchan, nopen, 2), &
                  xx(nchan, nopen, 2), wfp(nstat, nopen, 2))

        t = 0; v = 0; sp = 0; cp = 0

        ! construct, factorize and invert T
        forall (i = 1:nopen) t(i,i) = 1
        t = t + imu*km(1:nopen, 1:nopen)
        call invert_matrix(t)
        tt(1:nopen, 1:nopen, 1) = real(T, wp)
        tt(1:nopen, 1:nopen, 2) = aimag(t)

        ! evaluate the Coulomb functions and the normalization factor
        do i = 1, nchan
            call coul(l(i), ek(i), rmatr, f, fp, g, gp)
            v(i,i) = sqrt(2/(pi*k(i)))
            sp(i,i) = fp
            cp(i,i) = gp
        end do

        ! photoionization coefficient
        call dgemm('N', 'N', nchan, nopen,  nchan, rone, cp, nchan, km, nchan, rone,  sp,  nchan)  ! S' + C' K
        call dgemm('N', 'N', nchan, nopen2, nopen, rone, sp, nchan, tt, nchan, rzero, xx,  nchan)  ! (S' + C' K) T
        call dgemm('N', 'N', nchan, nopen2, nchan, rone, v,  nchan, xx, nchan, rzero, tt,  nchan)  ! V (S' + C' K) T
        call dgemm('T', 'N', nstat, nopen2, nchan, rone, w,  nchan, tt, nchan, rzero, wfp, nstat)  ! wT V (S' + C' K) T
        p = 0.5 / (eig - etot)
        do i = 1, nopen
            ak(:, i) = p * cmplx(wfp(:, i, 1), wfp(:, i, 2), wp)
        end do
        do i = nopen + 1, nchan
            ak(:, i) = 0
        end do

    end subroutine scatak


    real(wp) function cphase (l, k) result (sigma)

        integer,  intent(in) :: l
        real(wp), intent(in) :: k
#ifdef WITH_GSL
        type(gsl_sf_result) :: lnr, arg
        integer(c_int)      :: err
        real(c_double)      :: zr, zi

        zr = l + 1
        zi = -1/k
        err = gsl_sf_lngamma_complex_e(zr, zi, lnr, arg)
        sigma = arg % val
#else
        sigma = cphaz(l, -1/k, 0)
#endif

    end function cphase


    subroutine coul (l, Ek, r, F, Fp, G, Gp)

        integer,  intent(in)    :: l
        real(wp), intent(in)    :: Ek, r
        real(wp), intent(inout) :: F, Fp, G, Gp
        real(wp) :: Z = 1
#ifdef WITH_GSL
        call coul_gsl(l, ek, r, f, fp, g, gp)
#else
        call coul_ukrmol(l, ek, r, f, fp, g, gp)
#endif
    end subroutine coul


    subroutine coul_gsl (l, Ek, r, F, Fp, G, Gp)

        integer,  intent(in)    :: l
        real(wp), intent(in)    :: Ek, r
        real(wp), intent(inout) :: F, Fp, G, Gp
#ifdef WITH_GSL
        real(wp) :: Z = 1
        type(gsl_sf_result) :: resF, resFp, resG, resGp
        real(c_double)      :: expF, expG, eta, rho, lc, k, a, b, x, U, Up
        integer(c_int)      :: err, zero = 0

        if (ek > 0) then
            ! positive-energy solution
            k = sqrt(2*ek)
            rho = k*r
            eta = -z/k
            lc = l
            err = gsl_sf_coulomb_wave_fg_e(eta, rho, lc, zero, resf, resfp, resg, resgp, expf, expg)
            f  = resf  % val
            fp = resfp % val * k
            g  = resg  % val
            gp = resgp % val * k
        else
            ! negative-energy solution, Whittaker function W (WARNING: charged only)
            k = sqrt(-2*ek)
            a = l + 1 - z/k
            b = 2*l + 2
            x = 2*k*r
            u = gsl_sf_hyperg_u(a, b, x)
            up = -2*k*a * gsl_sf_hyperg_u(a + 1, b + 1, x)  ! DLMF §13.3.22
            f = 0
            fp = 0
            g = exp(-x/2) * x**(l + 1) * u  ! DLMF §13.14.3
            gp = exp(-x/2) * x**(l + 1) * (-k * u + (l + 1)/r * u + up)
        end if
#else
        f = 0; fp = 0
        g = 0; gp = 0
#endif
    end subroutine coul_gsl


    subroutine coul_ukrmol (l, Ek, r, F, Fp, G, Gp)

        integer,  intent(in)    :: l
        real(wp), intent(in)    :: Ek, r
        real(wp), intent(inout) :: F, Fp, G, Gp
        real(wp) :: Z = 1
        integer  :: ifail, mode = 1, kfn = 0, nterms, iwrite = 0
        real(wp) :: fc(l + 1), gc(l + 1), fcp(l + 1), gcp(l + 1), lc, efx, fnorm, acc = 1e-10_wp, k

        lc = l
        if (ek > 0) then
            ! positive-energy solution
            k = sqrt(2*ek)
            call coulfg(k*r, -z/k, lc, lc, fc, gc, fcp, gcp, mode, kfn, ifail)
            f = fc(l + 1)
            g = gc(l + 1)
            fp = fcp(l + 1) * k
            gp = gcp(l + 1) * k
        else
            ! negative energy solution (charged or neutral)
            if (z /= 0) then
                call couln(l, z, 2*ek, r, g, gp, acc, efx, ifail, nterms, fnorm)
            else
                call decay(2*ek, l, r, g, gp, iwrite)
            end if
            f = 0
            fp = 0
        end if

    end subroutine coul_ukrmol


    subroutine invert_matrix (T)

        complex(wp), allocatable, intent(inout) :: T(:, :)

        complex(wp), allocatable :: work(:)
        complex(wp)              :: nwork(1)

        integer(blasint), allocatable :: ipiv(:)
        integer(blasint)              :: m, info, lwork

        m = size(t, 1)

        ! calculate LU decomposition
        allocate (ipiv(m))
        call zgetrf(m, m, t, m, ipiv, info)

        ! query for workspace size and allocate memory
        lwork = -1
        call zgetri(m, t, m, ipiv, nwork, lwork, info)
        lwork = int(nwork(1))
        allocate (work(lwork))

        ! compute the inverse
        call zgetri(m, t, m, ipiv, work, lwork, info)

    end subroutine invert_matrix


    subroutine solve_complex_system (n, A, X, Y)

        complex(wp), allocatable, intent(inout) :: A(:, :)
        real(wp),    allocatable, intent(inout) :: X(:, :), Y(:, :)
        integer,                  intent(in)    :: n

        integer(blasint)              :: m, info, nrhs, lda
        integer(blasint), allocatable :: ipiv(:)
        complex(wp),      allocatable :: XX(:)

        allocate (ipiv(n))

        ! calculate LU decomposition
        m = n
        lda = size(a, 1)
        call zgetrf(m, m, a, lda, ipiv, info)

        ! combine real and imaginary part of the right-hand side
        xx = cmplx(x(1:n, 1), x(1:n, 2), wp)

        ! solve the equation
        nrhs = 1
        call zgetrs('N', m, nrhs, a, lda, ipiv, xx, m, info)

        y(1:n, 1) = real(XX, wp)
        y(1:n, 2) = aimag(xx)

    end subroutine solve_complex_system


    logical function next_permutation (p) result (ok)

        integer, intent(inout) :: p(:)
        integer                :: i, j, k, n

        n = size(p)

        ! no permutation in empty array
        if (n == 0) then
            ok = .false.
            return
        end if

        ! initialize a new permutation
        if (any(p < 1)) then
            forall (i = 1 : n) p(i) = i
            ok = .true.
            return
        end if

        ! one-element permutation cannot be advanced further
        if (n == 1) then
            ok = .false.
            return
        end if

        i = n

        do
            k = i
            i = i - 1
            if (p(i) < p(k)) then
                j = n
                do while (p(i) >= p(j))
                    j = j - 1
                end do
                call swap(p(i), p(j))
                call reverse(p(k:n))
                ok = .true.
                return
            else if (i == 1) then
                ok = .false.
                return
            end if
        end do

    end function next_permutation


    subroutine swap (a, b)

        integer, intent(inout) :: a, b
        integer                :: c

        c = a
        a = b
        b = c

    end subroutine swap


    subroutine reverse (a)

        integer, intent(inout) :: a(:)
        integer                :: i, n

        n = size(a)

        do i = 1, n/2
            call swap(a(i), a(n + 1 - i))
        end do

    end subroutine reverse


    subroutine kahan_add (X, dX, err)

        complex(wp), intent(inout) :: X, err
        complex(wp), intent(in)    :: dX

        complex(wp) :: Y, Z

        y = dx - err
        z = x + y
        err = (z - x) - y  ! the parenthesis must be obeyed, otherwise we get rubbish
        x = z

    end subroutine kahan_add

end module multidip_special