multidip_routines.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_routines
 
    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
 
#ifdef HAVE_NO_FINDLOC
    ! supply missing intrinsic functions
    use algorithms,  only: findloc
#endif
 
    ! GBTOlib
    use blas_lapack_gbl,  only: blasint
    use phys_const_gbl,   only: pi, imu, to_ev
    use precisn_gbl,      only: wp
 
    implicit none
 
    type intermediatestate
 
        integer :: order
        integer :: mgvn
        integer :: dcomp
 
        real(wp), allocatable :: ck(:, :, :)
 
        real(wp), allocatable :: ap(:, :, :)
 
        real(wp), allocatable :: dip(:, :, :)
 
        type(intermediatestate), pointer :: parent => null()
 
        type(intermediatestate), pointer :: prev => null()
        type(intermediatestate), pointer :: next => null()
 
    end type intermediatestate
 
    ! for debugging purposes
    logical, parameter :: test_expansion = .false.
 
contains
 
    subroutine multidip_main
 
        use, intrinsic :: iso_c_binding, only: c_ptr
 
        use mpi_gbl,          only: comm_rank => myrank, comm_size => nprocs, mpi_mod_finalize, mpi_mod_start, mpi_xermsg
        use multidip_io,      only: moleculardata, kmatrix, scatakcoeffs, read_wfncoeffs, read_kmatrices, read_molecular_data, &
                                    myproc, nprocs
        use multidip_params,  only: nmaxphotons, read_input_namelist, num_integ_algo
#ifdef WITH_GSL
        use multidip_special, only: gsl_set_error_handler_off
#endif
        use multidip_tests,   only: run_tests
 
        type(moleculardata)             :: moldat
        type(kmatrix),      allocatable :: km(:)
        type(scatakcoeffs), allocatable :: ak(:)
 
        character(len=256) :: arg, rmt_data, raw, msg
        type(c_ptr)        :: dummy
 
        integer     :: order, nirr, lusct(8), lukmt(8), lubnd, nkset(8), iarg, narg, input, i, erange(2)
        logical     :: verbose, mpiio
        real(wp)    :: omega(nMaxPhotons), first_IP, rasym
        complex(wp) :: polar(3, nMaxPhotons)
 
        call print_ukrmol_header(output_unit)
 
        print '(/,A)', 'Program MULTIDIP: Calculation of multi-photon ionization'
 
        iarg = 1
        narg = command_argument_count()
        input = input_unit
 
#ifdef WITH_GSL
        dummy = gsl_set_error_handler_off()
#endif
 
        do while (iarg <= narg)
            call get_command_argument(iarg, arg)
            if (arg == '--test') then
                call run_tests
                return
            else if (arg == '--input') then
                iarg = iarg + 1
                call get_command_argument(iarg, arg)
                open (input, file = trim(arg), action = 'read', form = 'formatted')
            else if (arg == '--help') then
                print '(/,A,/)', 'Available options:'
                print '(A)', '  --help          display this help and exit'
                print '(A)', '  --input FILE    use FILE as the input file (instead of standard input)'
                print '(A)', '  --test          run internal sanity checks and exit'
                print *
                return
            else
                write (msg, '(3a)') 'Unknown command line argument "', trim(arg), '"'
                call mpi_xermsg('multidip_routines', 'multidip_main', trim(msg), 1, 1)
            end if
            iarg = iarg + 1
        end do
 
        call read_input_namelist(input, order, lusct, lukmt, lubnd, nkset, polar, omega, verbose, rmt_data, first_ip, rasym, raw, &
                                 erange, mpiio)
 
        if (input /= input_unit) close (input)
 
        call mpi_mod_start
 
        myproc = comm_rank + 1
        nprocs = comm_size
        print '(/,A,I0,A,I0)', 'Parallel mode; image ', myproc, ' of ', nprocs
 
        print '(/,A,/)', 'Absorbed photons'
        do i = 1, order
            print '(2x,A,I0,3x,A,SP,2F6.3,A,2F6.3,A,2F6.3,A,F0.3,A)', &
                '#', i, 'eX: ', polar(1, i), 'i, eY: ', polar(2, i), 'i, eZ: ', polar(3, i), 'i, energy: ', omega(i) * to_ev, ' eV'
        end do
 
        if (first_ip > 0) then
            print '(/,A,F0.5,A)', 'Override first IP to ', first_ip * to_ev, ' eV'
        else
            print '(/,A)', 'Use calculated IP'
        end if
 
        nirr = count(lukmt > 0)
 
        allocate (km(nirr))
 
        if (any(lusct > 0)) then
            allocate (ak(nirr))
            call read_wfncoeffs(ak, lusct)
        end if
 
        call read_kmatrices(km, lukmt, nkset)
        call read_molecular_data(moldat, trim(rmt_data), mpiio, test_expansion)
 
        if (any(erange /= 0)) then
            if (erange(1) < 1 .or. erange(1) > erange(2) .or. erange(2) > minval(km % nescat)) then
                print '(/,A,I0)', 'Error: Energy range must satisfy 1 <= erange(1) <= erange(2) <= ', minval(km % nescat)
            end if
        else
            erange(1) = 1
            erange(2) = minval(km % nescat)
        end if
 
        if (rasym > moldat % rmatr) then
            if (order <= 2) then
                write (*, '(/,A,F0.1,A,F0.1,A)', advance = 'no') 'Interval ', moldat % rmatr, ' to ', rasym, &
                    ' will be integrated numerically using '
                select case (num_integ_algo)
                    case (1);  print '(a)', 'Romberg quadrature'
                    case (2);  print '(a)', 'Levin quadrature'
                end select
            else
                call mpi_xermsg('multidip_routines', 'multidip_main', &
                                'Error: Numerical integration (rasym) is currently only implemented for the second order.', 1, 1)
            end if
        end if
 
        call multidip_driver(order, moldat, km, ak, lubnd, omega(1:order), polar(:, 1:order), verbose, first_ip, rasym, raw, erange)
 
        call mpi_mod_finalize
 
    end subroutine multidip_main
 
 
    subroutine multidip_driver (order, moldat, km, ak, lubnd, omega, polar, verbose, first_IP, r0, raw, erange)
 
        use multidip_integ,  only: nested_cgreen_integ_cache
        use multidip_io,     only: moleculardata, kmatrix, scatakcoeffs, nprocs, get_diptrans, write_energy_grid
        use multidip_params, only: cache_integrals, extend_istate
 
        integer,                         intent(in) :: order, erange(2), lubnd
        logical,                         intent(in) :: verbose
        type(moleculardata),             intent(in) :: moldat
        type(kmatrix),      allocatable, intent(in) :: km(:)
        type(scatakcoeffs), allocatable, intent(in) :: ak(:)
        real(wp),                        intent(in) :: omega(:), first_IP, r0
        complex(wp),                     intent(in) :: polar(:, :)
        character(len=*),                intent(in) :: raw
 
        type(intermediatestate), pointer :: states, state
 
        integer,     allocatable :: iidip(:), ifdip(:)
        real(wp),    allocatable :: escat(:)
 
        integer  :: j, s, mgvni, mgvnn, mgvn1, mgvn2, nesc, irri, iki, icomp, nstati, nchani, max_dim
        real(wp) :: Ei
 
        iki    = 1                                       ! which K-matrix file (and Ak file) to use for the initial state
        mgvni  = km(iki) % mgvn                          ! IRR of the initial state
        nesc   = erange(2) - erange(1) + 1               ! number of scattering energies
        irri   = findloc(moldat % mgvns, mgvni, 1)       ! internal index of initial IRR in molecular_data
        nstati = moldat % mnp1(irri)                     ! number of inner region states in initial IRR
        nchani = moldat % nchan(irri)                    ! number of outer region channels in initial IRR
 
        allocate (states)
        allocate (states % ck(nstati, 2, (nesc + nprocs - 1) / nprocs))
        allocate (states % ap(nchani, 2, (nesc + nprocs - 1) / nprocs))
 
        ! initialize Coulomb-Green integral caching
        if (cache_integrals) then
            max_dim = merge(order, order - 1, extend_istate)
            call nested_cgreen_integ_cache % init(max_dim, (nesc + nprocs - 1) / nprocs, moldat % rmatr, r0, 0._wp)
        end if
 
        ! prepate the initial state
        call setup_initial_state(states, moldat, irri, lubnd, ei)
 
        ! for all absorbed photons
        do j = 1, order
 
            print '(/,2(A,I0))', 'Photon ', j, ' of ', order
 
            ! for all components of the dipole operator
            do icomp = 1, 3
 
                call get_diptrans(moldat, icomp, iidip, ifdip)
 
                ! for all transitions mediated by this component that are stored in molecular_data
                do s = 1, size(iidip)
 
                    mgvn1 = iidip(s) - 1
                    mgvn2 = ifdip(s) - 1
 
                    ! apply this transition to all relevant previous intermediate states
                    state => states
                    do while (associated(state))
 
                        ! process states of previous order
                        if (state % order + 1 == j) then
 
                            ! IRR of the previous intermediate state (i.e. last element in MGVN history chain)
                            mgvnn = state % mgvn
 
                            ! skip this transition if it is not applicable to the current IRR or there are no K-matrices for it
                            if ((mgvnn == mgvn1 .or. mgvnn == mgvn2) .and. &
                                any(km(:) % mgvn == mgvn1) .and. &
                                any(km(:) % mgvn == mgvn2)) then
 
                                ! either calculate the next intermediate state or evaluate the dipole elements
                                if (j < order) then
                                    call solve_intermediate_state(moldat, order, omega, icomp, s, mgvnn, &
                                                                  mgvn1, mgvn2, km, state, verbose, ei, first_ip, r0, erange)
                                else
                                    call extract_dipole_elements(moldat, order, omega, icomp, s, mgvnn, &
                                                                 mgvn1, mgvn2, km, ak, state, verbose, ei, first_ip, r0, erange)
                                end if
 
                            end if
 
                        end if
 
                        ! move on to the next state
                        state => state % next
 
                    end do
 
                end do
 
            end do
 
        end do
 
        ! calculate observables
        escat = km(iki) % escat(erange(1):erange(2))
        call write_energy_grid(escat)
        call calculate_pw_transition_elements(moldat, order, states, escat, ei, first_ip, omega, polar, erange)
        call calculate_asymmetry_parameters(moldat, order, states, escat, ei, first_ip, omega, raw, erange)
 
        ! clean up the linked list
        state => states
        do while (associated(state % next))
            state => state % next
            deallocate (state % prev)
        end do
        deallocate (state)
 
    end subroutine multidip_driver
 
 
    subroutine setup_initial_state (states, moldat, irr, lubnd, Ei)
 
        use multidip_io,     only: moleculardata, apply_boundary_amplitudes, read_bndcoeffs
        use multidip_outer,  only: test_outer_expansion, evaluate_fundamental_solutions
        use multidip_params, only: extend_istate
 
        type(intermediatestate), pointer, intent(inout) :: states
        type(moleculardata),              intent(in)    :: moldat
        integer,                          intent(in)    :: irr, lubnd
        real(wp),                         intent(out)   :: Ei
 
        complex(wp), allocatable :: ap(:)
        real(wp),    allocatable :: Ek(:), bnd(:, :), wbnd(:, :), fc(:), gc(:), fcp(:), gcp(:)
        real(wp)                 :: r
        integer                  :: mye, nmye, nchan, nstat, ichan, stot, mgvn, nbnd
 
        nmye = size(states % ck, 3)
        mgvn = moldat % mgvns(irr)
        stot = moldat % stot(irr)
        nchan = moldat % nchan(irr)
        nstat = moldat % mnp1(irr)
        nbnd = 1
        r = moldat % rmatr
        ei = moldat % eig(1, irr)
 
        ! construct representation of the initial state in the inner region
        states % order = 0
        states % mgvn = mgvn
        states % dcomp = 0
 
        allocate (bnd(nstat, nbnd))
 
        if (lubnd <= 0) then
            bnd(:, 1) = 0    ! all inner region eigenstates are unpopulated...
            bnd(1, 1) = 1    ! ... except for the first one, which has the coefficient 1.0
        else
            call read_bndcoeffs(bnd(:, 1), ei, lubnd, mgvn, stot)   ! get calculated coefs from BOUND
        end if
 
        do mye = 1, nmye
            states % ck(:, 1, mye) = bnd(:, 1)
            states % ck(:, 2, mye) = 0
        end do
 
        print '(/,a,/)', 'Bound state information'
        print '(4x,a,e25.15)', 'First R-matrix pole       = ', moldat % eig(1, irr)
        print '(4x,a,e25.15)', 'Calculated initial energy = ', ei
 
        ! construct representation of the initial state in the outer region
        if (extend_istate) then
 
            allocate (wbnd(nchan, nbnd), fc(nchan), gc(nchan), fcp(nchan), gcp(nchan), ap(nchan))
 
            call apply_boundary_amplitudes(moldat, irr, 'N', bnd, wbnd)
 
            ! "photoelectron" "kinetic" energies (actually negative closed channel thresholds)
            ek = ei - moldat % etarg(moldat % ichl(: , irr))
 
            call evaluate_fundamental_solutions(moldat, r, irr, nchan, ek, fc, gc, fcp, gcp, sqrtknorm = .false.)
 
            ! obtain outer coeffs by projecting the inner state on channels at boundary and dividing by decaying Whittaker function
            ap = wbnd(:, 1) / gc
 
            do mye = 1, nmye
                states % ap(:, 1, mye) = real(ap)
                states % ap(:, 2, mye) = aimag(ap)
            end do
 
            print '(4x,a,e25.15)', 'Total inner norm          = ', norm2(bnd(:, 1))
            print '(/,a,/)', 'Bound state outer region channel coefficients and amplitudes'
            print '(4x,a,4x,a,23x,a,20x,a,20x,a)', 'chan', 'Ek', 'Re ap', 'Im ap', 'f_p'
            do ichan = 1, nchan
                print '(i6,a,4e25.15)', ichan, ': ', ek(ichan), ap(ichan), wbnd(ichan, 1)
            end do
 
        else
 
            states % ap(:, :, :) = 0    ! outer region is empty
 
        end if
 
        if (test_expansion) then
            call test_outer_expansion('initial-state.txt', moldat, irr, states % ck(:, :, 1), states % ap(:, :, 1), ei)
        end if
 
    end subroutine setup_initial_state
 
 
    subroutine solve_intermediate_state (moldat, order, Ephoton, icomp, s, mgvnn, mgvn1, mgvn2, &
                                         km, state, verbose, calc_Ei, first_IP, r0, erange)
 
        use multidip_io,      only: moleculardata, kmatrix, scatakcoeffs, myproc, nprocs, apply_boundary_amplitudes, &
                                    apply_dipole_matrix, scale_boundary_amplitudes
        use multidip_outer,   only: evaluate_fundamental_solutions
        use multidip_params,  only: closed_interm, closed_range, compt
        use multidip_special, only: coul, solve_complex_system
 
        type(moleculardata),    intent(in) :: moldat
        type(kmatrix),          intent(in) :: km(:)
 
        type(intermediatestate), pointer, intent(inout) :: state
        type(intermediatestate), pointer :: last
 
        integer,  intent(in) :: order, icomp, s, mgvnn, mgvn1, mgvn2, erange(2)
        logical,  intent(in) :: verbose
        real(wp), intent(in) :: Ephoton(:), calc_Ei, first_IP, r0
 
        integer  :: nstatn, nstatf, nchann, nchanf, nopen, mgvnf, irrn, irrf, ikn, ikf, ipw, nesc, ie, mye, t, dt, i
        real(wp) :: Ei, Etotf
        character(len=1) :: transp
        character(len=1024) :: filename
 
        complex(wp), allocatable :: beta(:), app(:), H(:, :), lambda(:), hc(:), hcp(:)
        real(wp),    allocatable :: fc(:), fcp(:), gc(:), gcp(:)
 
        integer,  allocatable :: lf(:)
        real(wp), allocatable :: Ef(:), P(:), rhs(:, :), omega(:), echl(:)
        real(wp), allocatable :: Pw(:, :), wPw(:, :), wPD(:, :), IwPD(:, :), wckp(:, :), Dpsi(:, :, :), corr(:, :), ckp(:, :)
 
        if (mgvnn == mgvn1) then
            mgvnf = mgvn2
            transp = 'T'
        else
            mgvnf = mgvn1
            transp = 'N'
        end if
 
        irrn = findloc(moldat % mgvns, mgvnn, 1)
        irrf = findloc(moldat % mgvns, mgvnf, 1)
 
        ikn = findloc(km(:) % mgvn, mgvnn, 1)
        ikf = findloc(km(:) % mgvn, mgvnf, 1)
 
        nstatn = moldat % mnp1(irrn)
        nstatf = moldat % mnp1(irrf)
 
        nchann = moldat % nchan(irrn)
        nchanf = moldat % nchan(irrf)
 
        nesc = erange(2) - erange(1) + 1
        mye = 0
        t = 0
 
        allocate (pw(nstatf, nchanf), wpd(nchanf, 2), iwpd(nchanf, 2), wpw(nchanf, nchanf), wckp(nchanf, 2), omega(order), &
                  dpsi(nstatf, 2, (nesc + nprocs - 1) / nprocs), fc(nchanf), gc(nchanf), fcp(nchanf), gcp(nchanf), &
                  hc(nchanf), hcp(nchanf), corr(nchanf, 2), echl(nchanf), &
                  ef(nchanf), lf(nchanf), ckp(nstatf, 2), h(nchanf, nchanf), lambda(nchanf), beta(nchanf), &
                  app(nchanf))
 
        echl(1:nchanf) = moldat % etarg(moldat % ichl(1:nchanf, irrf)) - moldat % etarg(1)
 
        ! set up a new intermediate state at the end of the storage chain
        last => state
        do while (associated(last % next))
            last => last % next
        end do
        allocate (last % next)
        last % next % prev => last
        last => last % next
        last % order = state % order + 1
        last % mgvn = mgvnf
        last % dcomp = icomp
        last % parent => state
        allocate (last % ck(nstatf, 2, (nesc + nprocs - 1) / nprocs))
        allocate (last % ap(nchanf, 2, (nesc + nprocs - 1) / nprocs))
 
        call print_transition_header(last)
 
        ! apply the inner dipoles matrix on the previous inner region solution expansions
        call reset_timer(t, dt)
        call apply_dipole_matrix(moldat, icomp, s, transp, nstatf, nstatn, state % ck, dpsi)
        call reset_timer(t, dt)
        print '(4x,A,I0,A)', 'inner dipoles applied in ', dt, ' s'
 
        ! for all scattering energies
        do ie = erange(1) + myproc - 1, erange(2), nprocs
 
            mye = mye + 1
 
            ! calculate photon energies and the initial state energy (take into account user-specified first ionization potential)
            ei = calc_ei
            call calculate_photon_energies(first_ip, km(ikf) % escat(ie), moldat % etarg(1), ei, ephoton, omega)
 
            ! quantum numbers of the final state
            etotf = ei + sum(omega(1 : last % order))
            ef = km(ikf) % escat(ie) - echl - sum(omega(last % order + 1 :))
            lf = moldat % l2p(1:nchanf, irrf)
            nopen = merge(count(ef + closed_range > 0), count(ef > 0), closed_interm)
 
            print '(4x,A,*(I0,A))', 'energy ', ie, ' of ', km(ikf) % nescat, &
                                    ' (', count(ef > 0), ' open / ', nopen, ' open + weakly closed channels of ', nchanf, ')'
 
            ! evaluate Coulomb functions at this energy for all partial waves in the final IRR
            fc = 0; fcp = 0; gc = 0; gcp = 0; hc = 0; hcp = 0; lambda = 0
            call evaluate_fundamental_solutions(moldat, moldat % rmatr, irrf, nopen, ef, fc, gc, fcp, gcp, sqrtknorm = .false.)
            do ipw = 1, nopen
                hc(ipw)  = gc(ipw) + imu*fc(ipw)
                hcp(ipw) = gcp(ipw) + imu*fcp(ipw)
                lambda(ipw) = hcp(ipw) / hc(ipw)
            end do
 
            ! inner region part of the right-hand side of the master equation
            rhs = dpsi(:, :, mye)
 
            ! outer region correction to the right-hand side
            call multiint(moldat, r0, ei, km(ikf) % escat(ie) - sum(omega(last % order + 1 :)), omega, mye, last, +1, beta)
            beta = -2 * beta / sqrt(2*abs(ef))
            where (ef < 0) beta = beta / imu
            corr(:, 1) = real(beta * (fcp - lambda * fc), wp)
            corr(:, 2) = aimag(beta * (fcp - lambda * fc))
            call apply_boundary_amplitudes(moldat, irrf, 'T', corr, dpsi(:, :, mye))
            rhs = rhs - 0.5 * dpsi(:, :, mye)
 
            ! solve the master equation using the Sherman-Morrison-Woodbury formula
            p = 1 / (etotf - moldat % eig(1:nstatf, irrf))
            ckp(:, 1) = p * rhs(:, 1)
            ckp(:, 2) = p * rhs(:, 2)
 
            if (nopen > 0) then
 
                ! compose the reduced inverse Hamiltonian
                h = 0
                call scale_boundary_amplitudes(moldat, irrf, p, pw)
                call apply_boundary_amplitudes(moldat, irrf, 'N', pw, wpw)
                h(1:nopen, 1:nopen) = wpw(1:nopen, 1:nopen)
                do i = 1, nopen
                    h(i, i) = h(i, i) + 2/lambda(i)
                end do
                do i = nopen + 1, nchanf
                    h(i, i) = 1
                end do
                call apply_boundary_amplitudes(moldat, irrf, 'N', ckp, wpd)
 
                ! solve the reduced system
                iwpd = 0
                call solve_complex_system(nopen, h, wpd, iwpd)
 
                ! expand the reduced solution to the full solution
                call apply_boundary_amplitudes(moldat, irrf, 'T', iwpd, ckp)
                ckp(:, 1) = p * (rhs(:, 1) - ckp(:, 1))
                ckp(:, 2) = p * (rhs(:, 2) - ckp(:, 2))
 
            end if
 
            ! extract the outer region channel amplitudes
            call apply_boundary_amplitudes(moldat, irrf, 'N', ckp, wckp)
            app = cmplx(wckp(:, 1), wckp(:, 2), wp)
            where (hc /= 0)
                app = (app - beta*fc) / hc
            elsewhere
                app = 0
            end where
 
            ! for debugging purposes
            if (test_expansion) then
                wckp(:, 1) = real(app)
                wckp(:, 2) = aimag(app)
                write (filename, '(a,i0,a,i0,a)') 'intermediate-state-', irrf, '-', ie, '.txt'
                call test_outer_expansion(trim(filename), moldat, irrf, ckp, wckp, etotf)
            end if
 
            if (verbose) print '(4x,A,*(E25.15))', ' - beta: ', beta(1:nopen)
            if (verbose) print '(4x,A,*(E25.15))', ' - ap:   ', app(1:nopen)
 
            last % ck(:, 1, mye) = ckp(:, 1)
            last % ck(:, 2, mye) = ckp(:, 2)
 
            last % ap(:, 1, mye) = real(app, wp)
            last % ap(:, 2, mye) = aimag(app)
 
        end do
 
        call reset_timer(t, dt)
        print '(4x,A,I0,A)', 'intermediate states solved in ', dt, ' s'
 
    end subroutine solve_intermediate_state
 
 
    subroutine extract_dipole_elements (moldat, order, Ephoton, icomp, s, mgvnn, mgvn1, mgvn2, &
                                        km, ak, state, verbose, calc_Ei, first_IP, r0, erange)
 
        use multidip_io,      only: moleculardata, kmatrix, scatakcoeffs, myproc, nprocs, apply_dipole_matrix
        use multidip_outer,   only: evaluate_fundamental_solutions, calculate_k_matrix
        use multidip_params,  only: cone, czero, compt, maxtarg, custom_kmatrices
        use multidip_special, only: calculate_s_matrix, calculate_t_matrix, blas_zgemv => zgemv
 
        type(scatakcoeffs), allocatable, intent(in) :: ak(:)
 
        type(moleculardata),    intent(in) :: moldat
        type(kmatrix),          intent(in) :: km(:)
 
        type(intermediatestate), pointer, intent(inout) :: state
        type(intermediatestate), pointer                :: last
 
        integer,  intent(in) :: order, icomp, s, mgvnn, mgvn1, mgvn2, erange(2)
        logical,  intent(in) :: verbose
        real(wp), intent(in) :: Ephoton(:), calc_Ei, first_IP, r0
 
        character(len=1) :: transp
        character(len=1024) :: filename
 
        complex(wp), allocatable :: Af(:, :), dm(:), dp(:), Sm(:, :), tmat(:, :)
        real(wp), allocatable :: Dpsi(:, :, :), Ef(:), echlf(:), d_inner(:, :), d_outer(:, :), Re_Af(:, :), Im_Af(:, :)
        real(wp), allocatable :: omega(:), d_total(:, :), kmat(:, :), Sf(:), Cf(:), Sfp(:), Cfp(:)
        integer,  allocatable :: lf(:)
 
        real(wp) :: Etotf, Ei, Ek
        integer  :: mgvnf, irrn, irrf, ikn, ikf, nstatn, nstatf, nchann, nchanf, nesc, ie, nopen, nochf, mye, nene, maxchan, i
        integer  :: t, dt
 
        integer(blasint) :: m, n, lds, inc = 1
 
        if (mgvnn == mgvn1) then
            mgvnf = mgvn2
            transp = 'T'
        else
            mgvnf = mgvn1
            transp = 'N'
        end if
 
        irrn = findloc(moldat % mgvns, mgvnn, 1)
        irrf = findloc(moldat % mgvns, mgvnf, 1)
 
        ikn = findloc(km(:) % mgvn, mgvnn, 1)
        ikf = findloc(km(:) % mgvn, mgvnf, 1)
 
        nstatn = moldat % mnp1(irrn)
        nstatf = moldat % mnp1(irrf)
 
        nchann = km(ikn) % nchan
        nchanf = km(ikf) % nchan
 
        nesc = erange(2) - erange(1) + 1
        nene = (nesc + nprocs - 1) / nprocs
        mye = 0
        t = 0
 
        maxchan = nchanf
        if (maxtarg > 0) then
            maxchan = count(moldat % ichl(1:nchanf, irrf) <= maxtarg)
        end if
 
        allocate (dpsi(nstatf, 2, nene), omega(order), ef(nchanf), lf(nchanf), echlf(nchanf), &
                  re_af(nstatf, maxchan), im_af(nstatf, maxchan), af(nstatf, maxchan), sm(nchanf, nchanf), &
                  d_inner(maxchan, 2), d_outer(maxchan, 2), d_total(maxchan, 2), dm(maxchan), dp(nchanf), &
                  sf(nchanf), cf(nchanf), sfp(nchanf), cfp(nchanf), kmat(nchanf, nchanf), tmat(nchanf, nchanf))
 
        echlf(1:nchanf) = moldat % etarg(moldat % ichl(1:nchanf, irrf)) - moldat % etarg(1)
 
        ! set up a final state at the end of the storage chain
        last => state
        do while (associated(last % next))
            last => last % next
        end do
        allocate (last % next)
        last % next % prev => last
        last => last % next
        last % order = order
        last % mgvn = mgvnf
        last % dcomp = icomp
        last % parent => state
        allocate (last % dip(maxchan, 2, nene))
 
        call print_transition_header(last)
 
        ! apply the inner dipoles matrix on the previous inner region solution expansions
        call reset_timer(t, dt)
        call apply_dipole_matrix(moldat, icomp, s, transp, nstatf, nstatn, state % ck, dpsi)
        call reset_timer(t, dt)
        print '(4x,A,I0,A)', 'inner dipoles applied in ', dt, ' s'
 
        ! position K-matrix file pointer at the beginning of (this process') K-matrix data
        call km(ikf) % reset_kmatrix_position(skip = myproc - 1)
        if (allocated(ak)) call ak(ikf) % reset_wfncoeffs_position(skip = myproc - 1)
 
        ! for all scattering energies
        do ie = myproc, erange(2), nprocs
 
            ! read the K-matrix and Ak coeffs from the file and then skip the next 'nprocs - 1' records used by other images
            call km(ikf) % get_kmatrix(kmat, skip = nprocs - 1)
            if (allocated(ak)) call ak(ikf) % get_wfncoeffs(ek, re_af, im_af, skip = nprocs - 1)
            if (ie < erange(1)) cycle
 
            mye = mye + 1
            last % dip(:, :, mye) = 0
 
            ! calculate photon energies and the initial state energy (take into account user-specified first ionization potential)
            ei = calc_ei
            call calculate_photon_energies(first_ip, km(ikf) % escat(ie), moldat % etarg(1), ei, ephoton, omega)
 
            ! quantum numbers of the final state
            etotf = ei + sum(omega)
            ef = km(ikf) % escat(ie) - echlf        ! photoelectron energies in all (open and closed) channels
            lf = moldat % l2p(1:nchanf, irrf)       ! photoelectron angular momentum in all (open and closed) channels
            nopen = count(ef > 0)                   ! number of all open channels
            nochf = min(nopen, maxchan)             ! number of open channels for which we want to obtain dipoles
 
            print '(4x,A,*(I0,A))', 'energy ', ie, ' of ', km(ikf) % nescat, ' (', nopen, ' open channels of ', nchanf, ')'
 
            ! get values and derivatives of the fundamental solutions of uncoupled (or dipole-coupled) outer region equations
            call evaluate_fundamental_solutions(moldat, moldat % rmatr, irrf, nchanf, ef, sf, cf, sfp, cfp, sqrtknorm = .true.)
 
            ! if desired, recalculate the K-matrix here (possibly to include long-range dipole)
            if (custom_kmatrices) then
                call calculate_k_matrix(moldat, nopen, irrf, etotf, sf, cf, sfp, cfp, kmat)
            end if
            call calculate_t_matrix(kmat, tmat, nopen)
 
            ! contract D*psi with complex-conjugated wave function coefficients Ak to get inner region contribution
            if (.not. allocated(ak)) then
                call apply_ak_coefficients_multidip(dpsi(:, :, mye), d_inner, moldat, nopen, irrf, etotf, &
                                                    sfp, cfp, kmat, tmat, .true.)
            else
                call apply_ak_coefficiens_compak(dpsi(:, :, mye), d_inner, re_af(:, 1:nochf), im_af(:, 1:nochf))
            end if
 
            ! outer region correction
            dm = 0
            dp = 0
            call multiint(moldat, r0, ei, km(ikf) % escat(ie), omega, mye, last, -1, dm(1:nochf))
            call multiint(moldat, r0, ei, km(ikf) % escat(ie), omega, mye, last, +1, dp(1:nopen))
            dm(1:nochf) = dm(1:nochf) * imu / sqrt(2*pi*sqrt(2*ef(1:nochf)))
            dp(1:nopen) = dp(1:nopen) * imu / sqrt(2*pi*sqrt(2*ef(1:nopen)))
            if (verbose) print '(5x,A,1x,*(1x,E25.15))', 'dm = ', dm(1:nochf)
            if (verbose) print '(5x,A,1x,*(1x,E25.15))', 'dp = ', dp(1:nopen)
            if (any(dp /= 0) .or. test_expansion) then
                m   = int(nopen, blasint)
                n   = int(nochf, blasint)
                lds = int(nchanf, blasint)
                call calculate_s_matrix(tmat, sm, nopen)
                call blas_zgemv('T', m, n, -cone, sm, lds, dp, inc, cone, dm, inc)  ! dm = dm - S^T dp
            end if
            d_outer(:, 1) = real(dm)
            d_outer(:, 2) = aimag(dm)
 
            ! for debugging purposes
            if (test_expansion) then
                write (filename, '(a,i0,a,i0,a)') 'final-state-', irrf, '-', ie, '.txt'
                call test_final_expansion(trim(filename), moldat, irrf, nopen, etotf, ef, sfp, cfp, kmat, tmat)
            end if
 
            ! combine and store the partial wave dipoles
            d_total = d_inner + d_outer
            last % dip(:, :, mye) = d_total
 
            if (verbose) print '(5x,A,1x,*(1x,E25.15))', 'd_inner = ', (d_inner(i, 1), d_inner(i, 2), i = 1, nochf)
            if (verbose) print '(5x,A,1x,*(1x,E25.15))', 'd_outer = ', (d_outer(i, 1), d_outer(i, 2), i = 1, nochf)
            if (verbose) print '(5x,A,1x,*(1x,E25.15))', 'd_total = ', (d_total(i, 1), d_total(i, 2), i = 1, nochf)
 
        end do
 
        call reset_timer(t, dt)
        print '(4x,A,I0,A)', 'matrix elements calculated in ', dt, ' s'
 
    end subroutine extract_dipole_elements
 
 
    subroutine print_transition_header (state)
 
        use multidip_params, only: compt
 
        type(intermediatestate), pointer, intent(in) :: state
        type(intermediatestate), pointer :: ptr
 
        integer, allocatable :: irreds(:), compts(:)
        integer :: order, i
 
        order = state % order
 
        allocate (irreds(order + 1), compts(order))
 
        ptr => state
        irreds(order + 1) = ptr % mgvn
        do while (associated(ptr % parent))
            compts(order) = ptr % dcomp
            ptr => ptr % parent
            irreds(order) = ptr % mgvn
            order = order - 1
        end do
 
        write (*, '(2x,a)', advance = 'no') 'Transition'
        do i = 1, state % order
            write (*, '(1x,i0,1x,3a)', advance = 'no') irreds(i), '-[', compt(compts(i)), ']->'
        end do
        write (*, '(1x,i0)') irreds(state % order + 1)
 
    end subroutine print_transition_header
 
 
    subroutine apply_ak_coefficiens_compak (psi, Apsi, ReAk, ImAk)
 
        use multidip_params,  only: rone, rzero
        use multidip_special, only: blas_dgemv => dgemv
 
        real(wp), intent(in)    :: psi(:, :), ReAk(:, :), ImAk(:, :)
        real(wp), intent(inout) :: Apsi(:, :)
 
        integer(blasint) :: m, n, inc = 1
 
        m = int(size(reak, 1), blasint)
        n = int(size(reak, 2), blasint)
 
        call blas_dgemv('T', m, n,  rone, reak, m, psi(:, 1), inc, rzero, apsi(:, 1), inc)
        call blas_dgemv('T', m, n, +rone, imak, m, psi(:, 2), inc, +rone, apsi(:, 1), inc)
        call blas_dgemv('T', m, n,  rone, reak, m, psi(:, 2), inc, rzero, apsi(:, 2), inc)
        call blas_dgemv('T', m, n, -rone, imak, m, psi(:, 1), inc, +rone, apsi(:, 2), inc)
 
    end subroutine apply_ak_coefficiens_compak
 
 
    subroutine apply_ak_coefficients_multidip (psi, Apsi, moldat, nopen, irr, Etot, Sp, Cp, kmat, tmat, conj)
 
        use multidip_io,      only: moleculardata, apply_boundary_amplitudes
        use multidip_params,  only: czero, rzero
 
        type(moleculardata), intent(in)    :: moldat
        real(wp),            intent(in)    :: psi(:, :), Etot
        real(wp),            intent(inout) :: Apsi(:, :)
        integer,             intent(in)    :: nopen, irr
        real(wp),            intent(in)    :: Sp(:), Cp(:), kmat(:, :)
        complex(wp),         intent(in)    :: tmat(:, :)
        logical,             intent(in)    :: conj
 
        real(wp),    allocatable :: tmps(:, :), tmpc(:, :)
        complex(wp), allocatable :: tmpo(:)
 
        integer(blasint) :: n, ldk, one = 1
 
        integer     :: nchan, nstat, i, j
        real(wp)    :: Fij
        complex(wp) :: z, T
 
        nstat = moldat % mnp1(irr)
        nchan = moldat % nchan(irr)
        n = nopen
        ldk = size(kmat, 1)
 
        allocate (tmps(nstat, 2), tmpc(nchan, 2), tmpo(nopen))
 
        ! multiply solution by R-matrix poles
        tmps(:, 1) = psi(:, 1) / (moldat % eig(1:nstat, irr) - etot)
        tmps(:, 2) = psi(:, 2) / (moldat % eig(1:nstat, irr) - etot)
 
        ! contract with the boundary amplitudes
        call apply_boundary_amplitudes(moldat, irr, 'N', tmps, tmpc)
 
        ! multiply by complex-conjugated matrix F' = (TS' + TC'K) / √2π
        do j = 1, nopen
            tmpo(j) = 0
            do i = 1, nchan
                fij = (merge(sp(i), rzero, i == j) + cp(i)*kmat(i, j)) / sqrt(2*pi)
                tmpo(j) = tmpo(j) + cmplx(tmpc(i, 1), tmpc(i, 2), wp)*fij
            end do
        end do
 
        ! multiply by complex-conjugated T (nopen-by-nopen)
        do j = 1, min(nopen, size(apsi, 1))
            apsi(j, 1) = 0
            apsi(j, 2) = 0
            do i = 1, nopen
                t = merge(1, 0, i == j) + conjg(tmat(i, j))/2  ! = (1 + iK)⁻¹
                if (conj) t = conjg(t)
                z = tmpo(i) * t
                apsi(j, 1) = apsi(j, 1) + real(z)
                apsi(j, 2) = apsi(j, 2) + aimag(z)
            end do
        end do
 
        ! erase the closed-channel amplitudes
        do j = nopen + 1, min(nchan, size(apsi, 1))
            apsi(j, 1) = 0
            apsi(j, 2) = 0
        end do
 
    end subroutine apply_ak_coefficients_multidip
 
 
    subroutine test_final_expansion (filename, moldat, irr, nopen, Etot, Ek, Sp, Cp, kmat, tmat)
 
        use multidip_io,    only: moleculardata
        use multidip_outer, only: evaluate_fundamental_solutions
 
        character(len=*),       intent(in) :: filename
        type(moleculardata),    intent(in) :: moldat
        integer,                intent(in) :: irr, nopen
        real(wp),               intent(in) :: Etot, Ek(:), Sp(:), Cp(:), kmat(:, :)
        complex(wp),            intent(in) :: tmat(:, :)
 
        real(wp),    allocatable :: psi(:, :), Apsi(:, :), S(:), C(:), dSdr(:), dCdr(:)
        complex(wp), allocatable :: Fqp(:, :), Hp(:), Hm(:)
 
        real(wp)    :: dr, r
        integer     :: i, nfdm, u, nchan, p, q, nstat
 
        nfdm = size(moldat % r_points) - 1
        nchan = moldat % nchan(irr)
        nstat = moldat % mnp1(irr)
        dr = 0.1
 
        allocate (fqp(nchan, nopen), psi(nstat, 2), apsi(nchan, 2), s(nchan), c(nchan), hp(nchan), hm(nchan), dsdr(nchan), &
                  dcdr(nchan))
 
        open (newunit = u, file = filename, action = 'write', form = 'formatted')
 
        ! inner region
        do i = 1, nfdm
            write (u, '(e25.15)', advance = 'no') moldat % r_points(i)
            do q = 1, nchan
                if (.not. associated(moldat % wmat2(i, irr) % mat)) then
                    print '(a)', 'Error: wmat2 has not been read from molecular_data'
                    stop 1
                else if (moldat % wmat2(i, irr) % distributed) then
                    print '(a)', 'Error: test_outer_expansion not implemented in MPI-IO mode'
                    stop 1
                end if
                ! the "wave-function" is a row of a boundary amplitudes matrix to reduce the Ak coefficient with
                psi(:, 1) = moldat % wmat2(i, irr) % mat(q, :)
                psi(:, 2) = 0
                ! perform the reduction w*A
                call apply_ak_coefficients_multidip(psi, apsi, moldat, nopen, irr, etot, sp, cp, kmat, tmat, .false.)
                ! combine components to a single complex number
                fqp(q, 1:nopen) = cmplx(apsi(1:nopen, 1), apsi(1:nopen, 2), wp)
            end do
            write (u, '(*(e25.15))') fqp(1:nchan, 1:nopen)
        end do
 
        ! outer region
        do i = 0, nint(moldat % rmatr/dr)
            r = moldat % rmatr + i*dr
            write (u, '(e25.15)', advance = 'no') r
            call evaluate_fundamental_solutions(moldat, r, irr, nchan, ek, s, c, dsdr, dcdr)
            hp = (c + imu*s) * (-imu) / sqrt(2*pi)
            hm = (c - imu*s) * (-imu) / sqrt(2*pi)
            do p = 1, nopen
                do q = 1, nchan
                    if (q == p) fqp(q, p) = hp(q)
                    if (q /= p) fqp(q, p) = 0
                    fqp(q, p) = fqp(q, p) - hm(q) * conjg(1 + tmat(q, p))
                    write (u, '(*(e25.15))', advance = 'no') fqp(q, p)
                end do
            end do
            write (u, '()')
        end do
 
        close (u)
 
    end subroutine test_final_expansion
 
 
    subroutine reset_timer (t, dt)
 
        use mpi_gbl, only: mpi_mod_barrier
 
        integer, intent(inout) :: t
        integer, intent(out)   :: dt
        integer                :: clk_now, clk_rate, clk_max, ierr
 
        call mpi_mod_barrier(ierr)
 
        call system_clock(clk_now, clk_rate, clk_max)
 
        if (t < clk_now) then
            dt = (clk_now - t) / clk_rate
        else
            ! Cater for possible system clock wrap-around. Also pay attention to integer overflow; clk_max
            ! is typically the upper limit of the given integer type.
            dt = ((clk_max - t) + clk_now) / clk_rate
        end if
 
        t = clk_now
 
    end subroutine reset_timer
 
 
    subroutine calculate_photon_energies (first_IP, escat, etarg, Ei, Ephoton, omega)
 
        real(wp),              intent(in)    :: first_IP, escat, etarg, Ephoton(:)
        real(wp),              intent(inout) :: Ei
        real(wp), allocatable, intent(inout) :: omega(:)
        real(wp)                             :: IP, Ef
 
        ! calculate total final energy
        ef = etarg + escat
 
        ! obtain the calculated ionization potential (including non-zero photoelectron energy)
        ip = ef - ei
 
        ! set the requested ionization potential by shifting the initial state energy
        if (first_ip > 0) then
            ip = first_ip + escat
            ei = ef - ip
        end if
 
        ! calculate photon energies
        where (ephoton /= 0)
            omega = ephoton
        elsewhere
            omega = (ip - sum(ephoton)) / (size(omega) - count(ephoton /= 0))
        end where
 
    end subroutine calculate_photon_energies
 
 
    subroutine multiint (moldat, r0, Ei, esc, omega, ie, state, sb, dip)
 
        use mpi_gbl,         only: mpi_xermsg
        use multidip_integ,  only: nested_cgreen_integ_cache
        use multidip_io,     only: moleculardata
        use multidip_params, only: cache_integrals, closed_interm, closed_range
 
        type(moleculardata),              intent(in) :: moldat
        type(intermediatestate), pointer, intent(in) :: state
 
        complex(wp),              intent(inout) :: dip(:)
        real(wp),    allocatable, intent(inout) :: omega(:)
        real(wp),                 intent(in)    :: Ei, esc, r0
        integer,                  intent(in)    :: ie, sb
 
        complex(wp), allocatable :: k(:)
        integer,     allocatable :: l(:), m(:)
 
        real(wp)    :: c, a, Ekf, ebase, echl, etarg
        integer     :: nphot, mgvn, ichan, irr, ichl
 
        dip = 0
 
        if (cache_integrals) nested_cgreen_integ_cache % ie = ie    ! use integral cache for this energy
 
        if (.not. associated(state % parent)) then      ! state with no parent should not be a final state in the first place
            call mpi_xermsg('multidip_routines', 'multiint', 'Runtime error: unexpected parent state', 1, 1)
        end if
 
        a = moldat % rmatr                              ! R-matrix radius
        c = 0                                           ! dipole damping factor
        nphot = state % order                           ! number of photons absorbed by the final state
        mgvn = state % mgvn                             ! MGVN of the final state
        irr = findloc(moldat % mgvns, mgvn, 1)          ! index of this spin-symmetry in molecular_data
        ebase = moldat % etarg(1)                       ! energy of the first ion state
 
        !$omp parallel if(.not. cache_integrals) default(none) shared(moldat, state, omega, dip, closed_range, closed_interm) &
        !$omp& private(k, l, m, ichan, ichl, etarg, echl, Ekf) &
        !$omp& firstprivate(r0, c, esc, sb, ie, Ei, ebase, irr, nphot)
 
        allocate (k(nphot + 1), l(nphot + 1), m(nphot))
 
        !$omp do schedule(dynamic)
 
        do ichan = 1, size(dip)
 
            ! set up quantum numbers for this channel
            ichl = moldat % ichl(ichan, irr)                        ! index of ion state this partial wave is coupled to
            etarg = moldat % etarg(ichl)                            ! energy of that ion state
            echl = etarg - ebase                                    ! excitation threshold of this ion state w.r.t. ion ground
            ekf = esc - echl                                        ! always positive for final states (but not for intermediate)
            if (ekf + closed_range <= 0) cycle                      ! ignore strongly closed channels
            if (ekf < 0 .and. .not. closed_interm) cycle            ! optionally ignore also weakly closed channels
            k(1) = sqrt(2*abs(ekf)); if (ekf < 0) k(1) = k(1) * imu
            l(1) = moldat % l2p(ichan, irr)
 
            ! evaluate the correction dipole integral for this final channel
            dip(ichan) = multiint_chain(moldat, r0, ei, esc - omega(nphot), omega, ie, c, 1, state, ichan, sb, k, l, m)
 
        end do
 
        !$omp end do
        !$omp end parallel
 
    end subroutine multiint
 
 
    recursive complex(wp) function multiint_chain (moldat, r0, Ei, esc, omega, ie, c, N, state, ichanf, sb, k, l, m) result (integ)
 
        use multidip_integ,  only: nested_cgreen_integ
        use multidip_io,     only: moleculardata
        use multidip_params, only: closed_interm, closed_range
 
        type(moleculardata),              intent(in) :: moldat
        type(intermediatestate), pointer, intent(in) :: state
        type(intermediatestate), pointer             :: parent
 
        complex(wp), allocatable, intent(inout) :: k(:)
        real(wp),    allocatable, intent(inout) :: omega(:)
        integer,     allocatable, intent(inout) :: l(:), m(:)
        real(wp),                 intent(in)    :: r0, ei, esc
        integer,                  intent(in)    :: ie, sb, n, ichanf
 
        complex(wp) :: ap, kr(n + 1)
        real(wp)    :: z, c, a, ek, qion, qpws, ebase, etarg, echl
        integer     :: mgvni, mgvnf, nchani, irri, irrf, ichani, ichl, nphot, dcomp, lr(n + 1), mr(n)
 
        integ = 0
 
        ! no need to recurse further if the current state is the initial state
        if (.not. associated(state % parent)) return
 
        parent => state % parent                        ! previous intermediate state
        nphot = parent % order                          ! number of photons absorbed by the previous intermediate state
        ebase = moldat % etarg(1)                       ! energy of the ground ion state
        a = moldat % rmatr                              ! R-matrix radius
        z = moldat % nz - moldat % nelc                 ! Residual ion charge
        dcomp = state % dcomp                           ! last polarisation component absorbed by the final state
        mgvnf = state % mgvn                            ! irreducible representation of the current state
        mgvni = parent % mgvn                           ! irreducible representation of the previos state
        irrf = findloc(moldat % mgvns, mgvnf, 1)        ! index of the current spin-symmetry in molecular_data
        irri = findloc(moldat % mgvns, mgvni, 1)        ! index of the previous spin-symmetry in molecular_data
        nchani = moldat % nchan(irri)                   ! number of channels in the previous state irreducible representation
 
        ! loop over all channels of the intermediate state that are dipole-coupled to ichanf
        do ichani = 1, nchani
 
            ! set up quantum numbers for this channel
            ichl = moldat % ichl(ichani, irri)                          ! index of ion state this pw is coupled to
            etarg = moldat % etarg(ichl)                                ! energy of that state
            echl = etarg - ebase                                        ! excitation threshold of this ion state w.r.t. ground ion
            ek = esc - echl                                             ! can become negative starting from some channel
            if (ek + closed_range <= 0) exit                            ! ignore strongly closed channels
            if (ek < 0 .and. .not. closed_interm) exit                  ! optionally ignore also weakly closed channels
            k(n + 1) = sqrt(2*abs(ek)); if (ek < 0) k(n + 1) = k(n + 1) * imu
            l(n + 1) = moldat % l2p(ichani, irri)
            kr(1 : n + 1) = k(n + 1 : 1 : -1)
            lr(1 : n + 1) = l(n + 1 : 1 : -1)
 
            ! dipole angular integrals
            qion = channel_coupling_ion(moldat, dcomp, irrf, irri, ichanf, ichani)
            qpws = channel_coupling_pws(moldat, dcomp, irrf, irri, ichanf, ichani)
 
            ! outer region expansion coefficient for this partial wave and energy
            ap = cmplx(parent % ap(ichani, 1, ie), parent % ap(ichani, 2, ie), wp)
 
            ! evaluate Coulomb-Green's integrals for dipole-coupled ions
            if (qion /= 0) then
                m(n) = 0
                mr(1:n) = m(n:1:-1)
                if (ap /= 0) integ = integ + qion * ap * nested_cgreen_integ(z, a, r0, c, n, +1, sb, mr, lr, kr)
                if (nphot > 0) integ = integ + qion * multiint_chain(moldat, r0, ei, esc - omega(nphot), &
                                                      omega, ie, c, n + 1, parent, ichani, sb, k, l, m)
            end if
 
            ! evaluate Coulomb-Green's integrals for dipole-coupled partial waves
            if (qpws /= 0) then
                m(n) = 1
                mr(1:n) = m(n:1:-1)
                if (ap /= 0) integ = integ + qpws * ap * nested_cgreen_integ(z, a, r0, c, n, +1, sb, mr, lr, kr)
                if (nphot > 0) integ = integ + qpws * multiint_chain(moldat, r0, ei, esc - omega(nphot), &
                                                      omega, ie, c, n + 1, parent, ichani, sb, k, l, m)
            end if
 
        end do
 
    end function multiint_chain
 
 
    real(wp) function channel_coupling_ion (moldat, dcomp, irrf, irri, ichanf, ichani) result (qion)
 
        use multidip_io,     only: moleculardata
        use multidip_params, only: carti
 
        type(moleculardata), intent(in) :: moldat
        integer,             intent(in) :: dcomp, irrf, irri, ichanf, ichani
 
        integer :: itargi, itargf, li, lf, mi, mf
 
        lf = moldat % l2p(ichanf, irrf); mf = moldat % m2p(ichanf, irrf); itargf = moldat % ichl(ichanf, irrf)
        li = moldat % l2p(ichani, irri); mi = moldat % m2p(ichani, irri); itargi = moldat % ichl(ichani, irri)
 
        if (lf /= li .or. mf /= mi) then
            qion = 0
        else
            qion = moldat % crlv(itargf, itargi, carti(dcomp))
        end if
 
    end function channel_coupling_ion
 
 
    real(wp) function channel_coupling_pws (moldat, dcomp, irrf, irri, ichanf, ichani) result (qpws)
 
        use multidip_io,     only: moleculardata
        use multidip_params, only: cartm
 
        type(moleculardata), intent(in) :: moldat
        integer,             intent(in) :: dcomp, irrf, irri, ichanf, ichani
 
        integer :: itargi, itargf, li, lf, mi, mf, ipwi, ipwf
 
        itargf = moldat % ichl(ichanf, irrf)
        itargi = moldat % ichl(ichani, irri)
 
        if (itargf /= itargi) then
            qpws = 0
        else
            lf = moldat % l2p(ichanf, irrf); mf = moldat % m2p(ichanf, irrf); ipwf = lf*(lf + 1) + mf + 1
            li = moldat % l2p(ichani, irri); mi = moldat % m2p(ichani, irri); ipwi = li*(li + 1) + mi + 1
            qpws = sqrt(4*pi/3) * moldat % gaunt(ipwf, ipwi, cartm(dcomp))
        end if
 
    end function channel_coupling_pws
 
 
    subroutine calculate_pw_transition_elements (moldat, order, state, escat, calc_Ei, first_IP, Ephoton, polar, erange)
 
        use multidip_io,      only: moleculardata, myproc, nprocs, write_partial_wave_moments, write_cross_section
        use multidip_params,  only: alpha, maxtarg
        use multidip_special, only: cphase
 
        type(moleculardata),              intent(in) :: moldat
        type(intermediatestate), pointer, intent(in) :: state
        type(intermediatestate), pointer             :: ptr, parent
 
        integer,     intent(in) :: order, erange(2)
        real(wp),    intent(in) :: escat(:), calc_Ei, first_IP, Ephoton(:)
        complex(wp), intent(in) :: polar(:, :)
        integer                 :: nesc, ntarg, nchan, mxchan, lf, ichan, itarg, ie, irr, mgvn, mye, nirr, nene
        real(wp)                :: Z, Ef, Ei, kf, sigma, Q
        complex(wp)             :: d, ej
        complex(wp), allocatable :: M(:, :, :)
        real(wp),    allocatable :: cs(:, :), omega(:)
 
        z = moldat % nz - moldat % nelc         ! residual charge
        nirr = size(moldat % mgvns)             ! number of irreducible representations in molecular_data
        nesc = size(escat)                      ! total number of scattering energies
        mxchan = maxval(moldat % nchan)         ! maximal number of channels per irreducible representation
        ntarg = size(moldat % etarg)            ! number of targes (residual ions)
        nene = (nesc + nprocs - 1) / nprocs     ! number of scattering energies managed by this image
 
        ! if not all final targets are required, reduce some of the above limits
        if (maxtarg > 0) then
            mxchan = 0
            do irr = 1, nirr
                nchan = moldat % nchan(irr)
                mxchan = max(mxchan, count(moldat % ichl(1:nchan, irr) <= maxtarg))
            end do
            ntarg = min(ntarg, maxtarg)
        end if
 
        ! set up final storage
        allocate (omega(order), m(mxchan, nene, 0:7), cs(2 + ntarg, nene))
        m = 0
        cs = 0
 
        ! the dipole elements stored in molecular_data do not contain the charge, so for each transition we need a factor of -1
        q = (-1)**order
 
        ! calculate fixed-in-space transition matrix elements
        ptr => state
        do while (associated(ptr))
            if (ptr % order == order) then
 
                ! get polarisation factor for this absorption chain
                ej = polar(ptr % dcomp, ptr % order)
                parent => ptr % parent
                do while (associated(parent % parent))
                    ej = ej * polar(parent % dcomp, parent % order)
                    parent => parent % parent
                end do
 
                ! add transition element contribution from this absorption chain
                mye = 0
                do ie = myproc, nesc, nprocs
                    mye = mye + 1
                    do ichan = 1, size(ptr % dip, 1)
                        d = ej * cmplx(ptr % dip(ichan, 1, mye), ptr % dip(ichan, 2, mye), wp)
                        m(ichan, mye, ptr % mgvn) = m(ichan, mye, ptr % mgvn) + q * d
                    end do
                end do
 
            end if
            ptr => ptr % next
        end do
 
        ! write partial wave transition moments without Coulomb phase factor
        call write_partial_wave_moments(moldat, m, nesc, '')
 
        ! process transition matrix elements
        mye = 0
        do ie = myproc, nesc, nprocs
 
            mye = mye + 1
            ei = calc_ei
            call calculate_photon_energies(first_ip, escat(ie), moldat % etarg(1), ei, ephoton, omega)
 
            ! add phase factors to transition matrix elements
            do irr = 1, nirr
                mgvn = moldat % mgvns(irr)
                nchan = min(int(moldat % nchan(irr)), mxchan)
                do ichan = 1, nchan
                    lf = moldat % l2p(ichan, irr)
                    ef = escat(ie) - moldat % etarg(moldat % ichl(ichan, irr)) + moldat % etarg(1)
                    if (ef > 0) then
                        kf = sqrt(2*ef)
                        sigma = cphase(z, lf, kf)
                        m(ichan, mye, mgvn) = m(ichan, mye, mgvn) * imu**(-lf) * (cos(sigma) + imu*sin(sigma))
                    end if
                end do
            end do
 
            ! calculate oriented cross section
            cs(1, mye) = omega(1) * to_ev
            do itarg = 1, ntarg
                cs(2 + itarg, mye) = 0
                do irr = 1, nirr
                    mgvn = moldat % mgvns(irr)
                    nchan = min(int(moldat % nchan(irr)), mxchan)
                    do ichan = 1, nchan
                        if (itarg == moldat % ichl(ichan, irr)) then
                            d = m(ichan, mye, mgvn)
                            cs(2 + itarg, mye) = cs(2 + itarg, mye) + real(d * conjg(d))
                        end if
                    end do
                end do
            end do
            cs(2, mye) = sum(cs(3:, mye))
            cs(2:, mye) = cs(2:, mye) * 2*pi * (2*pi*alpha)**order * product(omega)
 
        end do
 
        ! write partial wave transition moments with Coulomb phase factor
        call write_partial_wave_moments(moldat, m, nesc, '+cphase')
 
        ! write fixed-in-space cross section
        call write_cross_section(cs, nesc, erange, 'gen_photo_xsec.txt')
 
    end subroutine calculate_pw_transition_elements
 
 
    subroutine calculate_asymmetry_parameters (moldat, order, state, escat, calc_Ei, first_IP, Ephoton, raw, erange)
 
        use multidip_io,       only: moleculardata, myproc, nprocs, write_raw_dipoles, write_cross_section
        use multidip_params,   only: alpha, maxtarg
        use multidip_special,  only: cphase
 
        character(len=*),                 intent(in) :: raw
        type(moleculardata),              intent(in) :: moldat
        type(intermediatestate), pointer, intent(in) :: state
        type(intermediatestate), pointer             :: ptr, parent
 
        integer,     intent(in) :: order, erange(2)
        real(wp),    intent(in) :: escat(:), calc_Ei, first_IP, Ephoton(:)
 
        integer                  :: nesc, ii, li, mi, nchain, ichain
        integer                  :: ntarg, ichan, ie, mye, irr, L, nene, maxl, pw
        integer,     allocatable :: chains(:, :)
        real(wp),    allocatable :: omega(:), cs(:, :)
        complex(wp), allocatable :: M_xyz(:, :, :, :), M_sph(:, :, :, :), beta(:, :)
        complex(wp)              :: d
        real(wp)                 :: Z, Ek, k, sigma, Ei, Q
        character(len=50)        :: filename
 
        z = moldat % nz - moldat % nelc
        nesc = size(escat)                      ! total number of scattering energies
        nene = (nesc + nprocs - 1) / nprocs     ! number of scattering energies managed by this image
        maxl = maxval(moldat % l2p)             ! highest partial wave angular momentum
        ntarg = size(moldat % etarg)            ! number of targes (residual ions)
 
        ! if not all final targets are required, reduce some of the above limits
        if (maxtarg > 0) then
            ntarg = min(ntarg, maxtarg)
            maxl = maxval(moldat % l2p, moldat % ichl <= maxtarg)
        end if
 
        ! find out how many absorption chains we have
        nchain = 0
        ptr => state
        do while (associated(ptr))
            if (ptr % order == order) then
                nchain = nchain + 1
            end if
            ptr => ptr % next
        end do
 
        allocate (omega(order), chains(order, nchain), cs(ntarg, nene), beta(2 + ntarg, nesc), &
                  m_xyz(nene, (maxl + 1)**2, nchain, ntarg), &
                  m_sph(nene, (maxl + 1)**2, nchain, ntarg))
 
 
        m_xyz = 0   ! multiphoton transition matrix elements in Cartesian basis
        m_sph = 0   ! multiphoton transition matrix elements in spherical basis
 
        ! the dipole elements stored in molecular_data do not contain the charge, so for each transition we need a factor of -1
        q = (-1)**order
 
        ! assemble multi-photon molecular transition elements in Cartesian basis
        ichain = 0
        ptr => state
        do while (associated(ptr))
            if (ptr % order == order) then
 
                ! find this irreducible representation in molecular_data
                irr = findloc(moldat % mgvns, ptr % mgvn, 1)
 
                ! assemble polarisation components for this absorption chain
                ichain = ichain + 1
                parent => ptr
                do while (associated(parent % parent))
                    chains(parent % order, ichain) = mod(1 + parent % dcomp, 3) - 1  ! y ~ -1, z ~ 0, x ~ +1
                    parent => parent % parent
                end do
 
                ! copy all Cartesian matrix elements for this irreducible representation and absorption history to M_xyz
                do ichan = 1, size(ptr % dip, 1)
                    ii = moldat % ichl(ichan, irr)
                    li = moldat % l2p(ichan, irr)
                    mi = moldat % m2p(ichan, irr)
                    pw = li*li + li + mi + 1
                    mye = 0
                    do ie = myproc, nesc, nprocs
                        mye = mye + 1
                        ek = escat(ie) - moldat % etarg(ii) + moldat % etarg(1)
                        if (ek > 0) then
                            k = sqrt(2*ek)
                            sigma = cphase(z, li, k)
                            d = cmplx(ptr % dip(ichan, 1, mye), ptr % dip(ichan, 2, mye), wp)
                            if (.not. d == d) d = 0  ! clear NaNs
                            m_xyz(mye, pw, ichain, ii) = imu**(-li) * (cos(sigma) + imu*sin(sigma)) * q * d
                        end if
                    end do
                end do
 
            end if
            ptr => ptr % next
        end do
 
        if (raw == 'xyz' .or. raw == 'both') then
            call write_raw_dipoles(m_xyz, chains, nesc, 'xyz')
        end if
 
        ! transform the multiphoton transition matrix elements from Cartesian basis (M_xyz) to spherical basis (M_sph)
        call convert_xyz_to_sph (m_xyz, m_sph, maxl, chains)
 
        if (raw == 'sph' .or. raw == 'both') then
            call write_raw_dipoles(m_sph, chains, nesc, 'sph')
        end if
 
        ! evaluate and write the asymmetry parameters for all possible orders
        do l = 0, 2*order
 
            print '(/,A,I0)', 'Evaluating asymmetry parameter for L = ', l
 
            ! calculate the absolute asymmetry parameter
            call calculate_quadratic_dipole_sph(beta, l, maxl, chains, chains, ntarg, nesc, m_sph, m_sph)
 
            ! sum partial cross section, add prefactors
            mye = 0
            do ie = myproc, nesc, nprocs
                mye = mye + 1
                ei = calc_ei
                call calculate_photon_energies(first_ip, escat(ie), moldat % etarg(1), ei, ephoton, omega)
                beta(1, mye) = omega(1) * to_ev
                beta(2:, mye) = beta(2:, mye) * 2*pi * (2*pi*alpha)**order * product(abs(omega))
                if (l == 0) then
                    cs(:, mye) = real(beta(3:, mye)) ! partial cross sections per target
                    beta(2, mye) = sum(cs(:, mye))   ! total cross section
                else
                    where (cs(:, mye) > 0)
                        beta(3:, mye) = beta(3:, mye) / cs(:, mye)  ! scale partial asymmetry parameter by the partial cross section
                    end where
                    beta(2, :) = 0                                  ! there is no "total beta" for L > 0
                end if
            end do
 
            ! write the asymmetry parameter
            write (filename, '(A,I0,A,I0,A)') 'gen_', order, 'photo_beta_', l, '.txt'
            call write_cross_section(real(beta, wp), nesc, erange, filename)
 
        end do
 
    end subroutine calculate_asymmetry_parameters
 
 
    subroutine convert_xyz_to_sph (M_xyz, M_sph, maxl, chains)
 
        use dipelm_special_functions, only: sph_basis_transform_elm
 
        complex(wp), allocatable, intent(in)    :: M_xyz(:, :, :, :)
        complex(wp), allocatable, intent(inout) :: M_sph(:, :, :, :)
        integer,                  intent(in)    :: maxl, chains(:, :)
 
        integer     :: order, ntarg, nchain, i, ii, l, mi, mj, pwi, pwj, ichain, jchain
        complex(wp) :: cpl
 
        order = size(chains, 1)
        nchain = size(chains, 2)
        ntarg = size(m_xyz, 4)
 
        m_sph = 0
 
        do ii = 1, ntarg
            do ichain = 1, nchain
                do jchain = 1, nchain
                    do l = 0, maxl
                        do mi = -l, l
                            do mj = -l, l
                                cpl = conjg(sph_basis_transform_elm(l, mj, mi, 'Slm'))
                                do i = 1, order
                                    cpl = cpl * sph_basis_transform_elm(1, chains(i, jchain), chains(i, ichain), 'Slm')
                                end do
                                pwi = l*l + l + mi + 1
                                pwj = l*l + l + mj + 1
                                m_sph(:, pwj, jchain, ii) = m_sph(:, pwj, jchain, ii) + cpl * m_xyz(:, pwi, ichain, ii)
                            end do
                        end do
                    end do
                end do
            end do
        end do
 
    end subroutine convert_xyz_to_sph
 
 
    subroutine calculate_quadratic_dipole_sph (beta, L, maxl, chains1, chains2, ntarg, nesc, M1, M2)
 
        use mpi_gbl,          only: mpi_reduceall_inplace_sum_wp
        use multidip_io,      only: myproc, nprocs
        use multidip_special, only: beta_contraction_tensor
 
        complex(wp), allocatable, intent(inout) :: beta(:, :)
        complex(wp), allocatable, intent(in)    :: M1(:, :, :, :), M2(:, :, :, :)
        integer,                  intent(in)    :: L, maxl, chains1(:, :), chains2(:, :), ntarg, nesc
 
        complex(wp)           :: MM
        real(wp), allocatable :: T(:, :, :, :), buffer(:)
        integer               :: ie, mye, itarg, idx, li, mi, lj, mj, pwi, pwj, ichain, jchain, nchain1, nchain2, order1, order2
        integer               :: qi(size(chains1, 1)), qj(size(chains2, 1)), p(size(chains1, 1))
 
        p = 0  ! hard-code laboratory-frame polarisations to zero
        beta = 0
        order1 = size(chains1, 1); nchain1 = size(chains1, 2)
        order2 = size(chains2, 1); nchain2 = size(chains2, 2)
 
        if (order1 /= order2) then
            print '(A)', 'WARNING: calculate_quadratic_dipole_sph is implemented only equal absorption orders'
            return
        end if
 
        allocate (t((maxl + 1)**2, (maxl + 1)**2, nchain1, nchain2))
 
        ! erase the angular integrals storage
        !$omp parallel do default(none) private(pwi, pwj, ichain, jchain) shared(maxl, nchain1, nchain2, T) collapse(4)
        do jchain = 1, nchain2
            do ichain = 1, nchain1
                do pwj = 1, (maxl + 1)**2
                    do pwi = 1, (maxl + 1)**2
                        t(pwi, pwj, ichain, jchain) = 0
                    end do
                end do
            end do
        end do
 
        ! calculate angular integrals (distribute over images and threads)
        !$omp parallel do default(none) schedule(dynamic) private(idx, pwi, pwj, li, lj, mi, mj, ichain, jchain, qi, qj) &
        !$omp& shared(nchain1, nchain2, maxl, chains1, chains2, nprocs, myproc, L, p, T, order1)
        do idx = myproc, (maxl + 1)**4 * nchain1 * nchain2, nprocs
 
            ! unpack idx = ((pwi*(maxl + 1)**2 + pwj)*nchain1 + ichain - 1)*nchain2 + jchain
            pwi    = 1 +     (idx - 1) / (nchain2 * nchain1 * (maxl + 1)**2)    ! = 1, ..., (maxl + 1)^2
            pwj    = 1 + mod((idx - 1) / (nchain2 * nchain1), (maxl + 1)**2)    ! = 1, ..., (maxl + 1)^2
            ichain = 1 + mod((idx - 1) /  nchain2,  nchain1)                    ! = 1, ..., nchain1
            jchain = 1 + mod( idx - 1,    nchain2)                              ! = 1, ..., nchain2
 
            ! unpack pwi = li*li + li + mi + 1
            li = floor(sqrt(pwi - 1._wp))
            mi = pwi - li*li - li - 1
            qi = chains1(:, ichain)
 
            ! unpack pwj = lj*lj + lj + mj + 1
            lj = floor(sqrt(pwj - 1._wp))
            mj = pwj - lj*lj - lj - 1
            qj = chains2(:, jchain)
 
            t(pwi, pwj, ichain, jchain) = beta_contraction_tensor(l, order1, p, li, mi, qi, lj, mj, qj)
 
        end do
 
        ! allreduce the contraction tensor
        buffer = reshape(t, [size(t)])
        call mpi_reduceall_inplace_sum_wp(buffer, size(buffer))
        t = reshape(buffer, shape(t))
 
        ! calculate the asymmetry parameter for this image's energies (distribute over threads)
        !$omp parallel do default(none) private(idx, itarg, pwi, pwj, ichain, jchain, li, mi, lj, mj, mye, MM) reduction(+:beta) &
        !$omp& shared(ntarg, maxl, nchain1, nchain2, myproc, nesc, nprocs, M1, M2, T)
        do idx = 1, ntarg * (maxl + 1)**4 * nchain1 * nchain2
 
            ! unpack idx = (((itarg - 1)*(maxl + 1)**2 + pwi)*(maxl + 1)**2 + pwj)*nchain1 + ichain - 1)*nchain2 + jchain
            itarg  = 1 +     (idx - 1) / (nchain2 * nchain1 * (maxl + 1)**4)                    ! = 1, ..., ntarg
            pwi    = 1 + mod((idx - 1) / (nchain2 * nchain1 * (maxl + 1)**2), (maxl + 1)**2)    ! = 1, ..., (maxl + 1)^2
            pwj    = 1 + mod((idx - 1) / (nchain2 * nchain1), (maxl + 1)**2)                    ! = 1, ..., (maxl + 1)^2
            ichain = 1 + mod((idx - 1) /  nchain2,  nchain1)                                    ! = 1, ..., nchain1
            jchain = 1 + mod( idx - 1,    nchain2)                                              ! = 1, ..., nchain2
 
            ! unpack pwi = li*li + li + mi
            li = floor(sqrt(pwi - 1._wp))
            mi = pwi - li*li - li - 1
 
            ! unpack pwj = lj*lj + lj + mj
            lj = floor(sqrt(pwj - 1._wp))
            mj = pwj - lj*lj - lj - 1
 
            mye = 0
            do ie = myproc, nesc, nprocs
                mye = mye + 1
                mm = conjg(m1(mye, pwi, ichain, itarg)) * m2(mye, pwj, jchain, itarg)
                beta(2 + itarg, mye) = beta(2 + itarg, mye) + t(pwi, pwj, ichain, jchain) * mm
            end do
 
        end do
 
    end subroutine calculate_quadratic_dipole_sph
 
 
    subroutine calculate_quadratic_dipole_xyz (beta, L, maxl, chains1, chains2, ntarg, nesc, M1, M2)
 
        use multidip_io,      only: myproc, nprocs
        use multidip_special, only: cartesian_vector_component_product_average
 
        complex(wp), allocatable, intent(inout) :: beta(:, :)
        complex(wp), allocatable, intent(in)    :: M1(:, :, :, :), M2(:, :, :, :)
        integer,                  intent(in)    :: L, maxl, chains1(:, :), chains2(:, :), ntarg, nesc
 
        complex(wp) :: MM
        real(wp)    :: T
        integer     :: ie, mye, itarg, li, mi, lj, mj, pwi, pwj, qi(size(chains1, 1)), qj(size(chains2, 1)), ichain, jchain
 
        beta = 0
 
        if (l /= 0) then
            print '(A)', 'WARNING: calculate_quadratic_dipole_xyz is implemented only for L = 0'
            return
        end if
 
        ! calculate the asymmetry parameter
        !$omp parallel do collapse(2) default(none) private(ichain, jchain, qi, qj, T, pwi, mye, ie, MM) reduction(+:beta) &
        !$omp& shared(chains1, chains2, maxl, myproc, nesc, nprocs, ntarg, M1, M2)
        do ichain = 1, size(chains1, 2)
            do jchain = 1, size(chains2, 2)
                qi = chains1(:, ichain)
                qj = chains2(:, jchain)
                t = cartesian_vector_component_product_average([qi, qj])
                do pwi = 1, (maxl + 1)**2
                    mye = 0
                    do ie = myproc, nesc, nprocs
                        mye = mye + 1
                        do itarg = 1, ntarg
                            mm = conjg(m1(mye, pwi, ichain, itarg)) * m2(mye, pwi, jchain, itarg)
                            beta(2 + itarg, mye) = beta(2 + itarg, mye) + t * mm
                        end do
                    end do
                end do
            end do
        end do
 
    end subroutine calculate_quadratic_dipole_xyz
 
end module multidip_routines