multidip_outer.f90

Go to the documentation of this file.

! Copyright 2023
!
! 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/>.
!
!> \brief   Routines related to outer region asymptotic channels
!> \authors J Benda
!> \date    2023
!>
module multidip_outer
 
    use blas_lapack_gbl, only: blasint
    use phys_const_gbl,  only: pi, imu
    use precisn_gbl,     only: wp
 
    implicit none
 
contains
 
    !> \brief   Evaluate wfn in channels for inside and outside of the R-matrix sphere
    !> \authors J Benda
    !> \date    2023
    !>
    !> Obtain the single-electron radial wave function amplitudes in channels by projecting the inner wave function at all sampling
    !> radii contained in molecular_data. Continue these to twice the R-matrix radius by evaluating the outer expansion in the
    !> outer region. The outer expansion assumes only outgoing (or exponentially decaying real) radial functions.
    !>
    !> The results are saved into a text file, which has the radial distance in the first column and sampled radial wave-functions
    !> in the remaining columns, one per partial wave channel. Every first column gives real part, every second the imaginary part.
    !>
    !> \param[in]  filename  Name of a text file to write.
    !> \param[in]  moldat    Molecular data class.
    !> \param[in]  irr       Inde xof irreducible representation in molecular_data.
    !> \param[in]  ck        Inner region expansion coefficients (in terms of inner region eigenstates).
    !> \param[in]  ap        Outer region expansion coefficients (in terms of partial wave channel eigenfunctions).
    !> \param[in]  Etot      Total energy of the system.
    !>
    subroutine test_outer_expansion (filename, moldat, irr, ck, ap, Etot)
 
        use multidip_io, only: moleculardata
 
        character(len=*),       intent(in) :: filename
        type(moleculardata),    intent(in) :: moldat
        integer,                intent(in) :: irr
        real(wp),               intent(in) :: ck(:, :), ap(:, :), Etot
 
        real(wp),    allocatable :: f(:, :), Ek(:), S(:), C(:), dSdr(:), dCdr(:)
 
        integer  :: u, j, nfdm, nchan
        real(wp) :: r, dr
 
        nfdm = size(moldat % r_points) - 1
        nchan = moldat % nchan(irr)
        dr = 0.1
 
        allocate (f(nchan, 2), ek(nchan), s(nchan), c(nchan), dsdr(nchan), dcdr(nchan))
 
        open (newunit = u, file = filename, action = 'write', form = 'formatted')
 
        do j = 1, nfdm
            if (.not. associated(moldat % wmat2(j, irr) % mat)) then
                print '(a)', 'Error: wmat2 has not been read from molecular_data'
                stop 1
            else if (moldat % wmat2(j, irr) % distributed) then
                print '(a)', 'Error: test_outer_expansion not implemented in MPI-IO mode'
                stop 1
            end if
            f = matmul(moldat % wmat2(j, irr) % mat, ck)
            write (u, '(*(e25.15))') moldat % r_points(j), transpose(f)
        end do
 
        ek = etot - moldat % etarg(moldat % ichl(:, irr))
 
        do j = 0, nint(moldat % rmatr / dr)
            r = moldat % rmatr + j*dr
            call evaluate_fundamental_solutions(moldat, r, irr, nchan, ek, s, c, dsdr, dcdr, sqrtknorm = .false.)
            write (u, '(*(e25.15))') r, cmplx(ap(:, 1), ap(:, 2), wp) * (c + imu*s)
        end do
 
        close (u)
 
    end subroutine test_outer_expansion
 
 
    !> \brief   Evaluate asymptotic solutions at boundary
    !> \authors J Benda
    !> \date    2023
    !>
    !> Obtain amplitudes and derivatives of the fundamental solutions at region boundary.
    !>
    !> The evaluated functions use the same normalization as `ASYWFN` in *cfasym.f*. That is, the amplitudes contain an extra
    !> factor of 1/sqrt(k) in addition to the standard normalization of Coulomb functions. Similarly, the derivatives contain
    !> the same factor of 1/sqrt(k). Because the functions are differentiated with respect to 'r', this 1/sqrt(k)
    !> factor combines with the 'k' coming from derivative of the argument, giving sqrt(k) compared to standard normalization
    !> of the derivative of the Coulomb functions.
    !>
    !> \param[in] moldat  Molecular data class.
    !> \param[in] r       Evaluation radius.
    !> \param[in] irr     Irreducible representation index.
    !> \param[in] nopen   Restrict evaluation of solutions to this number of (open or closed) partial wave channels.
    !> \param[in] Ek      Photoelectron kinetic energy in each channel.
    !> \param[inout] S    Regular solution amplitude per partial wave channel.
    !> \param[inout] C    Irregular solution amplitude per partial wave channel.
    !> \param[inout] Sp   Regular solution derivative per partial wave channel.
    !> \param[inout] Cp   Irregular solution derivative per partial wave channel.
    !> \param[in] sqrtknorm   Set to .false. to avoid adding the 1/sqrt(k) factor discussed above.
    !>
    subroutine evaluate_fundamental_solutions (moldat, r, irr, nopen, Ek, S, C, Sp, Cp, sqrtknorm)
 
        use multidip_io,      only: moleculardata
        use multidip_special, only: coul
 
        type(moleculardata),    intent(in)    :: moldat
        real(wp),               intent(in)    :: r, Ek(:)
        integer,                intent(in)    :: irr, nopen
        real(wp),               intent(inout) :: S(:), C(:), Sp(:), Cp(:)
        logical, optional,      intent(in)    :: sqrtknorm
 
        real(wp) :: k, kfactor, Z, F, Fp, G, Gp
        integer  :: ichan, nchan, l
 
        z = moldat % nz - moldat % nelc
        nchan = moldat % nchan(irr)
 
        !$omp parallel do default(none) private(ichan, k, l, kfactor, F, Fp, G, Gp) &
        !$omp& shared(moldat, irr, Z, r, Ek, S, C, Sp, Cp, sqrtknorm, nchan, nopen)
        do ichan = 1, min(nchan, nopen)
 
            k = sqrt(2*abs(ek(ichan)))
            l = moldat % l2p(ichan, irr)
 
            ! k-dependent prefactor optionally added to the evaluated solutions
            kfactor = 1/sqrt(k)
            if (present(sqrtknorm)) then
                if (.not. sqrtknorm) then
                    kfactor = 1
                end if
            end if
 
            call coul(z, l, ek(ichan), r, f, fp, g, gp)
 
            s(ichan) = f * kfactor
            c(ichan) = g * kfactor
            sp(ichan) = fp * kfactor
            cp(ichan) = gp * kfactor
 
        end do
 
    end subroutine evaluate_fundamental_solutions
 
 
    !> \brief   Calculate (generalized) K-matrix
    !> \authors J Benda
    !> \date    2023
    !>
    !> Alternative to RSOLVE. Calculates the K-matrix without propagation, simply by projecting the inner wave function
    !> on the asymptotic channels.
    !>
    !> First, the R-matrix is calculated from boundary amplitudes and R-matrix poles.
    !> Finally, the generalized K-matrix is obtained by solution of the standard matrix equation
    !> \f[
    !>      (C - R C') K = R S' - S \,.
    !> \f]
    !> In this equation the diagonal matrices S and C consist of the regular and irregular solutions of the dipole-coupled
    !> outer region asymptotic equations.
    !>
    !> \param[in] moldat  Molecular data class.
    !> \param[in] nopen   Number of open channels.
    !> \param[in] irr     Irreducible representation index.
    !> \param[in] Etot    Total energy of the whole system.
    !> \param[in] S       Value of the regular asymptotic solution per partial wave channel.
    !> \param[in] C       Value of the irregular asymptotic solution per partial wave channel.
    !> \param[in] Sp      Derivative of the regular asymptotic solution per partial wave channel.
    !> \param[in] Cp      Derivative of the irregular asymptotic solution per partial wave channel.
    !> \param[inout] Kmat Kmatrix to calculate.
    !>
    subroutine calculate_k_matrix (moldat, nopen, irr, Etot, S, C, Sp, Cp, Kmat)
 
        use multidip_io,      only: moleculardata
        use multidip_params,  only: rzero
        use multidip_special, only: calculate_r_matrix, blas_dgetrf => dgetrf, blas_dgetrs => dgetrs
 
        type(moleculardata), intent(in) :: moldat
 
        integer,  intent(in)    :: nopen, irr
        real(wp), intent(in)    :: Etot
        real(wp), intent(in)    :: S(:), C(:), Sp(:), Cp(:)
        real(wp), intent(inout) :: Kmat(:, :)
 
        real(wp),         allocatable :: Amat(:, :)     ! matrix of equations to be composed and solved
        real(wp),         allocatable :: Rmat(:, :)     ! R-matrix in the standard partial wave channels
        integer(blasint), allocatable :: ipiv(:)        ! permutation array for use in zgetrf/s
 
        integer(blasint) :: ldk, info, n, nrhs
        integer          :: nstat, nchan, ichan, jchan
 
        nchan = moldat % nchan(irr)
        nstat = moldat % mnp1(irr)
 
        ! 1. Calculate the standard R-matrix
 
        allocate (rmat(nchan, nchan))
 
        call calculate_r_matrix(nchan, nstat, moldat % wamp(irr) % mat, moldat % eig(1:nstat, irr), etot, rmat)
 
        ! 2. Solve for the K-matrix
 
        allocate (amat(nchan, nchan), ipiv(nchan))
 
        !$omp parallel do collapse(2) private(ichan, jchan)
        do jchan = 1, nchan
            do ichan = 1, nchan
                amat(ichan, jchan) = merge(c(jchan), rzero, ichan == jchan) - rmat(ichan, jchan)*cp(jchan)
            end do
        end do
        !$omp parallel do collapse(2) private(ichan, jchan)
        do jchan = 1, nopen
            do ichan = 1, nchan
                kmat(ichan, jchan) = rmat(ichan, jchan)*sp(jchan) - merge(s(jchan), rzero, ichan == jchan)
            end do
        end do
        !$omp parallel do collapse(2) private(ichan, jchan)
        do jchan = nopen + 1, nchan
            do ichan = 1, nchan
                kmat(ichan, jchan) = rzero
            end do
        end do
 
        n    = int(nchan, blasint)
        nrhs = int(nopen, blasint)
        ldk  = int(size(kmat, 1), blasint)
 
        call blas_dgetrf(n, n, amat, n, ipiv, info)
        call blas_dgetrs('N', n, nrhs, amat, n, ipiv, kmat, ldk, info)
 
    end subroutine calculate_k_matrix
 
end module multidip_outer