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

    ! 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

contains

    subroutine multidip_main

        use multidip_io,     only: moleculardata, kmatrix, scatakcoeffs, read_wfncoeffs, read_kmatrices, read_molecular_data, &
                                   myproc, nprocs
        use multidip_params, only: nmaxphotons
        use multidip_tests,  only: run_tests

        type(MolecularData)             :: moldat
        type(KMatrix),      allocatable :: km(:)
        type(ScatAkCoeffs), allocatable :: ak(:)

        character(len=64) :: arg
        integer           :: order, nirr, lusct(8), lukmt(8), nkset(8), iarg, narg, input, i
        logical           :: verbose
        real(wp)          :: omega(nMaxPhotons)
        complex(wp)       :: polar(3, nMaxPhotons)

        namelist /mdip/ order, lusct, lukmt, nkset, polar, omega, verbose

        order = 1
        lusct = 0
        lukmt = 0
        nkset = 1
        polar = 0
        omega = 0
        verbose = .false.

        call print_ukrmol_header(output_unit)

        print '(/,A)', 'Program MULTIDIP: Calculation of multi-photon ionization'

        iarg = 1
        narg = command_argument_count()
        input = input_unit

        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
                print '(3A)', 'Unknown command line argument "', arg, '"'
                stop 1
            end if
            iarg = iarg + 1
        end do

        read (input, nml = mdip)

        if (order < 1) then
            print '(A)', 'Order must be positive'
            stop 1
        end if

        if (order > nmaxphotons) then
            print '(A,I0)', 'Order must be <= ', nmaxphotons
        end if

        polar(:, order + 1:nmaxphotons) = 0
        omega(order + 1:nmaxphotons) = 0

        if (input /= input_unit) close (input)

#ifdef WITH_COARRAYS
        myproc = this_image()
        nprocs = num_images()
        print '(/,A,I0,A,I0)', 'Parallel mode; image ', myproc, ' of ', nprocs
#endif

        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

        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)

        call multidip_driver(order, moldat, km, ak, omega, polar, verbose)

    end subroutine multidip_main


    subroutine multidip_driver (order, moldat, km, ak, omega, polar, verbose)

        use multidip_io, only: moleculardata, kmatrix, scatakcoeffs, nprocs, get_diptrans

        integer,                         intent(in) :: order
        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(:)
        complex(wp),                     intent(in) :: polar(:, :)

        type(IntermediateState), pointer :: states, state

        integer, allocatable :: iidip(:), ifdip(:)

        integer :: j, s, mgvni, mgvnn, mgvn1, mgvn2, nesc, irri, iki, icomp, nstati, nchani

        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   = km(iki) % nescat                        ! number of scattering energies
        irri   = minloc(abs(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))

        ! construct representation of the initial state
        states % order = 0
        states % mgvn = mgvni
        states % dcomp = 0
        states % ck(:, :, :) = 0
        states % ck(1, 1, :) = 1  ! only the first inner state is populated (at all scattering energies)
        states % ap(:, :, :) = 0

        ! 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, irri, icomp, s, &
                                                                  mgvnn, mgvn1, mgvn2, km, state, verbose)
                                else
                                    call extract_dipole_elements(moldat, order, omega, irri, icomp, s, &
                                                                 mgvnn, mgvn1, mgvn2, km, ak, state, verbose)
                                end if

                            end if

                        end if

                        ! move on to the next state
                        state => state % next

                    end do

                end do

            end do

        end do

        ! calculate partial wave dipoles and cross sections
        call calculate_observables(moldat, order, states, km(iki) % escat, moldat % eig(1, irri), omega, polar)

        ! 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 solve_intermediate_state (moldat, order, Ephoton, irri, icomp, s, mgvnn, mgvn1, mgvn2, km, state, verbose)

        use multidip_io,      only: moleculardata, kmatrix, scatakcoeffs, myproc, nprocs, apply_boundary_amplitudes, &
                                    apply_dipole_matrix, scale_boundary_amplitudes
        use multidip_params,  only: closed_interm, compt
        use multidip_special, only: coul, solve_complex_system

        type(MolecularData),        intent(in)    :: moldat
        type(KMatrix), allocatable, intent(in)    :: km(:)

        type(IntermediateState), pointer, intent(inout) :: state
        type(IntermediateState), pointer :: last

        integer,  intent(in) :: order, irri, icomp, s, mgvnn, mgvn1, mgvn2
        logical,  intent(in) :: verbose
        real(wp), intent(in) :: Ephoton(:)

        integer  :: nstatn, nstatf, nchann, nchanf, nopen, mgvnf, irrn, irrf, ikn, ikf, ipw, nesc, ie, mye
        real(wp) :: Ei, IP, Etotf, dE = 0.1
        character(len=1) :: transp

        complex(wp), allocatable :: beta(:), app(:), H(:, :), lambda(:), hc(:), hcp(:)

        integer,  allocatable :: lf(:)
        real(wp), allocatable :: Ef(:), fc(:), fcp(:), gc(:), gcp(:), P(:), rhs(:, :), omega(:)
        real(wp), allocatable :: Pw(:, :), wPw(:, :), wPD(:, :), IwPD(:, :), wckp(:, :), Dpsi(:, :, :), corr(:, :), ckp(:, :)

        integer :: i

        if (mgvnn == mgvn1) then
            mgvnf = mgvn2
            transp = 'T'
        else
            mgvnf = mgvn1
            transp = 'N'
        end if

        irrn = minloc(abs(moldat % mgvns - mgvnn), 1)  ! = findloc(moldat % mgvns, mgvnn, 1)
        irrf = minloc(abs(moldat % mgvns - mgvnf), 1)  ! = findloc(moldat % mgvns, mgvnf, 1)

        ikn = minloc(abs(km(:) % mgvn - mgvnn), 1)  ! = findloc(km(:) % mgvn, mgvnn, 1)
        ikf = minloc(abs(km(:) % mgvn - mgvnf), 1)  ! = findloc(km(:) % mgvn, mgvnf, 1)

        nstatn = moldat % mnp1(irrn)
        nstatf = moldat % mnp1(irrf)

        nchann = moldat % nchan(irrn)
        nchanf = moldat % nchan(irrf)

        ei = moldat % eig(1, irri)
        nesc = km(ikf) % nescat
        mye = 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), &
                  ef(nchanf), lf(nchanf), ckp(nstatf, 2), h(nchanf, nchanf), lambda(nchanf), beta(nchanf), &
                  app(nchanf))

        print '(2x,A,I0,A,I0,5A)', 'Transition ', mgvnn, ' -> ', mgvnf, ' (', compt(icomp), ', ', transp, ')'

        ! 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))

        ! apply the inner dipoles matrix on the previous inner region solution expansions
        call apply_dipole_matrix(moldat, icomp, s, transp, nstatf, nstatn, state % ck, dpsi)

        ! for all scattering energies
        do ie = myproc, nesc, nprocs

            ! quantum numbers of the final state
            mye = mye + 1
            ip = km(ikf) % escat(ie) + moldat % etarg(1) - ei
            call calculate_photon_energies(ip, ephoton, omega)
            etotf = ei + sum(omega(1 : last % order))
            ef = etotf - moldat % etarg(moldat % ichl(1:nchanf, irrf))
            lf = moldat % l2p(1:nchanf, irrf)
            nopen = count(ef > 0)

            if (verbose) print '(4x,*(A,I0))', 'energy ', ie, ' of ', nesc
            if (verbose) print '(4x,A,I0)', ' - open channels: ', nopen

            ! redefine the number of open channels to include also weakly closed ones
            if (closed_interm) then
                nopen = count(ef + de > 0)
                if (verbose) print '(4x,A,I0)', ' - open + weakly closed channels: ', nopen
            end if

            ! 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
            do ipw = 1, nopen
                call coul(lf(ipw), ef(ipw), moldat % rmatr, fc(ipw), fcp(ipw), gc(ipw), gcp(ipw))
                hc(ipw)  = cmplx(gc(ipw),  fc(ipw),  wp)
                hcp(ipw) = cmplx(gcp(ipw), fcp(ipw), wp)
                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, ei, 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)
                forall (i = 1:nopen) h(i, i) = h(i, i) + 2/lambda(i)
                forall (i = nopen + 1:nchanf) h(i, i) = 1
                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

            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

    end subroutine solve_intermediate_state


    subroutine extract_dipole_elements (moldat, order, Ephoton, irri, icomp, s, mgvnn, mgvn1, mgvn2, km, ak, state, verbose)

        use multidip_io,      only: moleculardata, kmatrix, scatakcoeffs, myproc, nprocs, apply_dipole_matrix
        use multidip_params,  only: rone, rzero, cone, czero, compt
        use multidip_special, only: scatak,calculate_s_matrix, blas_dgemv => dgemv, blas_zgemv => zgemv

        type(MolecularData),             intent(in)    :: moldat
        type(KMatrix),      allocatable, intent(in)    :: km(:)
        type(ScatAkCoeffs), allocatable, intent(in)    :: ak(:)

        type(IntermediateState), pointer, intent(inout) :: state
        type(IntermediateState), pointer                :: last

        integer,  intent(in) :: irri, order, icomp, s, mgvnn, mgvn1, mgvn2
        logical,  intent(in) :: verbose
        real(wp), intent(in) :: Ephoton(:)

        character(len=1) :: transp

        complex(wp), allocatable :: Af(:, :), dm(:), dp(:), Sm(:, :)
        real(wp), allocatable :: Dpsi(:, :, :), Ef(:), echlf(:), d_inner(:, :), d_outer(:, :), Re_Af(:, :), Im_Af(:, :)
        real(wp), allocatable :: omega(:), d_total(:, :)
        integer,  allocatable :: lf(:)

        real(wp) :: Etotf, Ei, IP
        integer  :: mgvnf, irrn, irrf, ikn, ikf, nstatn, nstatf, nchann, nchanf, nesc, ie, nopen, mye

        integer(blasint) :: m, n, inc = 1

        if (mgvnn == mgvn1) then
            mgvnf = mgvn2
            transp = 'T'
        else
            mgvnf = mgvn1
            transp = 'N'
        end if

        irrn = minloc(abs(moldat % mgvns - mgvnn), 1)  ! = findloc(moldat % mgvns, mgvnn, 1)
        irrf = minloc(abs(moldat % mgvns - mgvnf), 1)  ! = findloc(moldat % mgvns, mgvnf, 1)

        ikn = minloc(abs(km(:) % mgvn - mgvnn), 1)  ! = findloc(km(:) % mgvn, mgvnn, 1)
        ikf = minloc(abs(km(:) % mgvn - mgvnf), 1)  ! = findloc(km(:) % mgvn, mgvnf, 1)

        nstatn = moldat % mnp1(irrn)
        nstatf = moldat % mnp1(irrf)

        nchann = km(ikn) % nchan
        nchanf = km(ikf) % nchan

        ei = moldat % eig(1, irri)
        nesc = km(ikf) % nescat
        mye = 0

        allocate (dpsi(nstatf, 2, (nesc + nprocs - 1) / nprocs), d_inner(nchanf, 2), d_outer(nchanf, 2), dm(nchanf), dp(nchanf), &
                  sm(nchanf, nchanf), omega(order), d_total(nchanf, 2), ef(nchanf), lf(nchanf), echlf(nchanf))

        print '(2x,A,I0,A,I0,5A)', 'Transition ', mgvnn, ' -> ', mgvnf, ' (', compt(icomp), ', ', transp, ')'

        ! 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(nchanf, 2, nesc))

        ! apply the inner dipoles matrix on the previous inner region solution expansions
        call apply_dipole_matrix(moldat, icomp, s, transp, nstatf, nstatn, state % ck, dpsi)

        ! for all scattering energies
        do ie = myproc, nesc, nprocs

            ! quantum numbers of the final state
            mye = mye + 1
            etotf = km(ikf) % escat(ie) + moldat % etarg(1)
            ip = etotf - ei
            call calculate_photon_energies(ip, ephoton, omega)
            echlf = moldat % etarg(moldat % ichl(1:nchanf, irrf)) - moldat % etarg(1)
            ef = km(ikf) % escat(ie) - echlf
            nopen = count(ef > 0)
            lf = moldat % l2p(1:nchanf, irrf)

            if (verbose) print '(4x,*(A,I0))', 'energy ', ie, ' of ', nesc

            ! obtain wave function coefficients in the final IRR from the Ak file data structure, or calculate them on the fly
            if (allocated(ak)) then
                re_af = ak(ikf) % ReA(:, :, mye)
                im_af = ak(ikf) % ImA(:, :, mye)
            else
                if (.not. allocated(af)) allocate (af(nstatf, nchanf))

                call scatak(moldat % rmatr, &
                            moldat % eig(1:nstatf, irrf), &
                            etotf, &
                            moldat % wamp(1:nchanf, 1:nstatf, irrf), &
                            ef, &
                            lf, &
                            km(ikf) % kmat(:, :, mye), &
                            af)

                re_af = real(Af, wp)
                im_af = aimag(af)
            end if

            ! contract D*psi with complex-conjugated wave function coefficients Ak
            m = nstatf
            n = nchanf
            call blas_dgemv('T', m, n,  rone, re_af, m, dpsi(:, 1, mye), inc, rzero, d_inner(:, 1), inc)
            call blas_dgemv('T', m, n, +rone, im_af, m, dpsi(:, 2, mye), inc, +rone, d_inner(:, 1), inc)
            call blas_dgemv('T', m, n,  rone, re_af, m, dpsi(:, 2, mye), inc, rzero, d_inner(:, 2), inc)
            call blas_dgemv('T', m, n, -rone, im_af, m, dpsi(:, 1, mye), inc, +rone, d_inner(:, 2), inc)

            ! calculate S-matrix for this energy
            call calculate_s_matrix(km(ikf) % kmat(:, :, mye), sm, nopen)

            ! outer region correction
            call multiint(moldat, ei, omega, mye, last, -1, dm)
            call multiint(moldat, ei, omega, mye, last, +1, dp)
            where (ef > 0)
                dm = dm * imu / sqrt(2*pi*sqrt(2*ef))
                dp = dp * imu / sqrt(2*pi*sqrt(2*ef))
            elsewhere
                dm = 0
                dp = 0
            end where
            call blas_zgemv('T', n, n, -cone, sm, n, dp, inc, cone, dm, inc)  ! dm = dm - S^T dp
            d_outer(:, 1) = real(dm, wp)
            d_outer(:, 2) = aimag(dm)

            ! combine and store the partial wave dipoles
            d_total = d_inner + d_outer
            last % dip(:, :, ie) = d_total

            if (verbose) print '(5x,A,1x,*(1x,E25.15))', 'd_inner = ', (d_inner(n, 1), d_inner(n, 2), n = 1, nopen)
            if (verbose) print '(5x,A,1x,*(1x,E25.15))', 'd_outer = ', (d_outer(n, 1), d_outer(n, 2), n = 1, nopen)
            if (verbose) print '(5x,A,1x,*(1x,E25.15))', 'd_total = ', (d_total(n, 1), d_total(n, 2), n = 1, nopen)

        end do

    end subroutine extract_dipole_elements


    subroutine calculate_photon_energies (IP, Ephoton, omega)

        real(wp),              intent(in)    :: IP, Ephoton(:)
        real(wp), allocatable, intent(inout) :: omega(:)

        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, Ei, omega, ie, state, sb, dip)

        use multidip_io,     only: moleculardata
        use multidip_params, only: closed_interm

        type(MolecularData),              intent(in) :: moldat
        type(IntermediateState), pointer, intent(in) :: state

        complex(wp), allocatable, intent(inout) :: dip(:)
        real(wp),    allocatable, intent(inout) :: omega(:)
        real(wp),                 intent(in)    :: Ei
        integer,                  intent(in)    :: ie, sb

        complex(wp), allocatable :: k(:)
        integer,     allocatable :: l(:), m(:)

        real(wp)    :: c, a, Ekf, Etot
        integer     :: nphot, mgvn, nchan, ichan, irr

        dip = 0

        if (.not. associated(state % parent)) stop 1    ! state with no parent should not be a final state in the first place

        a = moldat % rmatr                              ! R-matrix radius
        c = 0                                           ! dipole damping factor
        nphot = state % order                           ! number of photons absorbed by the final state
        etot = ei + sum(omega(1:nphot))                 ! total energy of the final state
        mgvn = state % mgvn                             ! MGVN of the final state
        irr = minloc(abs(moldat % mgvns - mgvn), 1)     ! = findloc(moldat % mgvns, mgvn, 1)
        nchan = moldat % nchan(irr)                     ! number of channels in the final irreducible representation

        allocate (k(nphot), l(nphot), m(nphot - 1))

        do ichan = 1, nchan

            ! set up quantum numbers for this channel
            ekf = etot - moldat % etarg(moldat % ichl(ichan, irr))
            if (ekf < 0 .and. .not. closed_interm) cycle  ! ignore 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, ei, omega, ie, c, 1, state, ichan, sb, k, l, m)

        end do

    end subroutine multiint


    recursive complex(wp) function multiint_chain (moldat, Ei, 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

        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)    :: ei
        integer,                  intent(in)    :: ie, sb, n, ichanf

        complex(wp) :: ap, kr(n + 1)
        real(wp)    :: c, a, ek, qion, qpws, etot
        integer     :: i, mgvni, mgvnf, nchani, irri, ichani, nphot, dcomp, lr(n + 1), mr(n)

        integ = 0

        if (state % parent % order == 0) return         ! if the parent is the initial state, there are not tails to integrate

        a = moldat % rmatr                              ! R-matrix radius
        dcomp = state % dcomp                           ! last polarisation component absorbed by the final state
        parent => state % parent                        ! previous intermediate state
        nphot = parent % order                          ! number of photons absorbed by the previous intermediate state
        mgvnf = state % mgvn                            ! irreducible representation of the current state
        mgvni = parent % mgvn                           ! irreducible representation of the previos state
        irri = minloc(abs(moldat % mgvns - mgvni), 1)   ! = findloc(moldat % mgvns, mgvni, 1)
        nchani = moldat % nchan(irri)                   ! number of channels in the previous state irreducible representation
        etot = ei + sum(omega(1:nphot))                 ! total energy of the previous intermediate state

        ! 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
            ek = etot - moldat % etarg(moldat % ichl(ichani, irri))
            if (ek < 0 .and. .not. closed_interm) cycle  ! ignore 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)
            forall (i = 1:n+1) kr(i) = k(n + 2 - i)
            forall (i = 1:n+1) lr(i) = l(n + 2 - i)

            ! dipole angular integrals
            qion = channel_coupling_ion(moldat, dcomp, mgvnf, mgvni, ichanf, ichani)
            qpws = channel_coupling_pws(moldat, dcomp, mgvnf, mgvni, 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
                forall (i = 1:n) mr(i) = m(n + 1 - i)
                integ = integ + qion * ap * nested_cgreen_integ(a, c, n, +1, sb, mr, lr, kr)
                integ = integ + qion * multiint_chain(moldat, ei, 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
                forall (i = 1:n) mr(i) = m(n + 1 - i)
                integ = integ + qpws * ap * nested_cgreen_integ(a, c, n, +1, sb, mr, lr, kr)
                integ = integ + qpws * multiint_chain(moldat, ei, 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, mgvnf, mgvni, ichanf, ichani) result (Qion)

        use multidip_io, only: moleculardata

        type(moleculardata), intent(in) :: moldat
        integer,             intent(in) :: dcomp, mgvnf, mgvni, ichanf, ichani

        integer :: tr(3) = [ 3, 1, 2 ]  ! how crlv is stored: x is last, y first, z in between
        integer :: irri, irrf, itargi, itargf, li, lf, mi, mf

        irrf = minloc(abs(moldat % mgvns - mgvnf), 1)  ! = findloc(moldat % mgvns, mgvnf, 1)
        irri = minloc(abs(moldat % mgvns - mgvni), 1)  ! = findloc(moldat % mgvns, mgvni, 1)

        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, tr(dcomp))
        end if

    end function channel_coupling_ion


    real(wp) function channel_coupling_pws (moldat, dcomp, mgvnf, mgvni, ichanf, ichani) result (Qpws)

        use multidip_io, only: moleculardata

        type(moleculardata), intent(in) :: moldat
        integer,             intent(in) :: dcomp, mgvnf, mgvni, ichanf, ichani

        integer :: m(3) = [ +1, -1, 0 ]  ! real spherical harmonics Y1m associated with directions x,y,z
        integer :: irri, irrf, itargi, itargf, li, lf, mi, mf, ipwi, ipwf

        irrf = minloc(abs(moldat % mgvns - mgvnf), 1)  ! = findloc(moldat % mgvns, mgvnf, 1)
        irri = minloc(abs(moldat % mgvns - mgvni), 1)  ! = findloc(moldat % mgvns, mgvni, 1)

        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, m(dcomp))
        end if

    end function channel_coupling_pws


    subroutine calculate_observables (moldat, order, state, escat, Ei, Ephoton, polar)

        use multidip_io,      only: moleculardata, myproc, nprocs, write_partial_wave_moments, write_cross_section
        use multidip_params,  only: alpha
        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
        real(wp),    intent(in) :: escat(:)
        complex(wp), intent(in) :: polar(:, :)
        real(wp),    intent(in) :: Ei, Ephoton(:)

        integer                  :: i, nesc, mxchan, lf, ichan, ie, irr, mgvn
#ifdef WITH_COARRAYS
        real(wp),    allocatable :: cs(:, :)[:]
        complex(wp), allocatable :: M(:, :, :)[:]
#else
        real(wp),    allocatable :: cs(:, :)
        complex(wp), allocatable :: M(:, :, :)
#endif
        real(wp),    allocatable :: omega(:)
        real(wp)                 :: Ef, kf, sigma, navg, IP, Q
        complex(wp)              :: ej, d

        nesc = size(escat)
        mxchan = maxval(moldat % nchan)

        ! set up final storage
#ifdef WITH_COARRAYS
        allocate (cs(2, nesc)[*], m(mxchan, nesc, 0:7)[*])
#else
        allocate (cs(2, nesc), m(mxchan, nesc, 0:7))
#endif
        allocate (omega(order))

        cs = 0
        m = 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

        ! how many photons have unspecified (i.e. zero) polarization?
        navg = 0
        do i = 1, order
            if (all(polar(:, i) == 0)) then
                navg = navg + 1
            end if
        end do

        ! process contributions from all final states
        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
                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
                do ie = myproc, nesc, nprocs
                    do ichan = 1, size(ptr % dip, 1)
                        d = cmplx(ptr % dip(ichan, 1, ie), ptr % dip(ichan, 2, ie), wp)
                        m(ichan, ie, ptr % mgvn) = m(ichan, ie, ptr % mgvn) + q * ej * d
                        cs(2, ie) = cs(2, ie) + real(d*conjg(d), wp)
                    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, '')

        ! add prefactors
        do ie = myproc, nesc, nprocs
            ip = escat(ie) + moldat % etarg(1) - ei
            call calculate_photon_energies(ip, ephoton, omega)
            cs(1, ie) = omega(1) * to_ev
            do irr = 1, size(moldat % mgvns)
                mgvn = moldat % mgvns(irr)
                do ichan = 1, moldat % nchan(irr)
                    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(lf, kf)
                        m(ichan, ie, mgvn) = m(ichan, ie, mgvn) * imu**(-lf) * (cos(sigma) + imu*sin(sigma))
                    end if
                end do
            end do
            cs(2, ie) = cs(2, ie) * 2*pi * (2*pi*alpha)**order * product(omega) / 3**navg
        end do

        ! write partial wave transition moments with Coulomb phase factor
        call write_partial_wave_moments(moldat, m, '+cphase')

        ! write orientation-averaged cross section
        call write_cross_section(cs)

    end subroutine calculate_observables

end module multidip_routines