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/>.
!
module multidip_outer
 
    use blas_lapack_gbl, only: blasint
    use phys_const_gbl,  only: pi, imu
    use precisn_gbl,     only: wp
 
    implicit none
 
contains
 
    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
 
 
    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
 
 
    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 and transform it to the generalized angular momentum basis
 
        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