Postprocessing_module.F90

Go to the documentation of this file.

! Copyright 2019
!
! For a comprehensive list of the developers that contributed to these codes
! see the UK-AMOR website.
!
! This file is part of UKRmol-in (UKRmol+ suite).
!
!     UKRmol-in 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-in 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-in (in source/COPYING). Alternatively, you can also visit
!     <https://www.gnu.org/licenses/>.
 
module postprocessing_module
 
    use blas_lapack_gbl, only: blasint
    use cdenprop_defs,   only: civect
    use precisn,         only: shortint, wp
 
    implicit none
 
    private
 
    public postprocess
 
    ! data types used by RMT and expected in the 'molecular_data' file
    integer, parameter :: rmt_int  = shortint
    integer, parameter :: rmt_real = wp
 
    type :: outerinterface
        integer :: maxchan = 10000
        integer :: maxnmo
        integer :: ntarg
        integer :: nchan
        integer :: ismax
 
        integer  :: nfdm
        real(wp) :: rmatr
 
        real(wp),     allocatable :: r_points(:)
        type(civect), allocatable :: wamp(:)
 
        integer,  allocatable :: idtarg(:)
        integer,  allocatable :: irrchl(:)
        integer,  allocatable :: ichl(:)
        integer,  allocatable :: lchl(:)
        integer,  allocatable :: mchl(:)
        integer,  allocatable :: qchl(:)
        real(wp), allocatable :: echl(:)
        real(wp), allocatable :: a(:,:)
    contains
        procedure, public  :: init   => init_outer_interface
        procedure, public  :: deinit => deinit_outer_interface
        procedure, public  :: setup_amplitudes
 
        procedure, public  :: extract_data
        procedure, private :: get_channel_info
        procedure, private :: get_channel_couplings
        procedure, private :: get_boundary_data
 
        procedure, public  :: write_data
        procedure, private :: write_channel_info
        procedure, private :: write_boundary_data
    end type outerinterface
 
contains
 
    subroutine postprocess (SCATCI_input, solutions)
 
        use class_molecular_properties_data, only: molecular_properties_data
        use options_module,                  only: optionsset
        use solutionhandler_module,          only: solutionhandler
 
        type(optionsset),                   intent(in)    :: scatci_input
        type(solutionhandler), allocatable, intent(inout) :: solutions(:)
 
        type(molecular_properties_data) :: properties_all
 
        call redistribute_solutions(scatci_input, solutions)
        call cdenprop_properties(scatci_input, solutions, properties_all)
        call outer_interface(scatci_input, solutions, properties_all)
 
    end subroutine postprocess
 
 
    subroutine redistribute_solutions (SCATCI_input, solutions)
 
        use consts_mpi_ci,          only: pass_to_cdenprop, no_diagonalization
        use mpi_gbl,                only: mpiint, myrank, mpi_mod_bcast
        use options_module,         only: optionsset
        use solutionhandler_module, only: solutionhandler
        use parallelization_module, only: grid => process_grid
 
        type(optionsset),                   intent(in)    :: SCATCI_input
        type(solutionhandler), allocatable, intent(inout) :: solutions(:)
 
        type(solutionhandler) :: oldsolution
        integer(mpiint)       :: srcrank
        integer               :: i, nnuc, nstat
 
        ! nothing to do when the solutions are already uniformly distributed everywhere
        if (grid % sequential) return
 
        ! set up the uniformly distributed solutions
        do i = 1, size(scatci_input % opts)
            if (iand(scatci_input % opts(i) % vector_storage_method, pass_to_cdenprop) /= 0 .and. &
                scatci_input % opts(i) % diagonalization_flag /= no_diagonalization) then
 
                ! 1. Share metadata for the solutions from the master of group that found that solution to everyone.
 
                srcrank = grid % group_master_world_rank(grid % which_group_is_work(i))
 
                call mpi_mod_bcast(solutions(i) % core_energy,  srcrank)
                call mpi_mod_bcast(solutions(i) % vec_dimen,    srcrank)
 
                call mpi_mod_bcast(solutions(i) % ci_vec % nstat,           srcrank)
                call mpi_mod_bcast(solutions(i) % ci_vec % mgvn,            srcrank)
                call mpi_mod_bcast(solutions(i) % ci_vec % s,               srcrank)
                call mpi_mod_bcast(solutions(i) % ci_vec % sz,              srcrank)
                call mpi_mod_bcast(solutions(i) % ci_vec % nnuc,            srcrank)
                call mpi_mod_bcast(solutions(i) % ci_vec % e0,              srcrank)
                call mpi_mod_bcast(solutions(i) % ci_vec % CV_is_scalapack, srcrank)
                call mpi_mod_bcast(solutions(i) % ci_vec % mat_dimen_r,     srcrank)
                call mpi_mod_bcast(solutions(i) % ci_vec % mat_dimen_c,     srcrank)
 
                nnuc  = solutions(i) % ci_vec % nnuc
                nstat = solutions(i) % ci_vec % nstat
 
                if (.not. allocated(solutions(i) % ci_vec % ei)) then
                    allocate (solutions(i) % ci_vec % ei(nstat))
                end if
 
                call mpi_mod_bcast(solutions(i) % ci_vec % cname(1:nnuc),   srcrank)
                call mpi_mod_bcast(solutions(i) % ci_vec % charge(1:nnuc),  srcrank)
                call mpi_mod_bcast(solutions(i) % ci_vec % xnuc(1:nnuc),    srcrank)
                call mpi_mod_bcast(solutions(i) % ci_vec % ynuc(1:nnuc),    srcrank)
                call mpi_mod_bcast(solutions(i) % ci_vec % znuc(1:nnuc),    srcrank)
                call mpi_mod_bcast(solutions(i) % ci_vec % ei(1:nstat),     srcrank)
 
                ! 2. Backup the original solutions, which are distributed over the MPI group.
 
                oldsolution = solutions(i)
 
                if (.not. grid % is_my_group_work(i)) then
                    oldsolution % ci_vec % blacs_context = -1
                    oldsolution % ci_vec % descr_CV_mat = -1
                end if
 
                ! 3. Redistribute the eigenvectors over whole MPI world.
 
                solutions(i) % ci_vec % blacs_context   = grid % wcntxt
                solutions(i) % ci_vec % nprow           = grid % wprows
                solutions(i) % ci_vec % npcol           = grid % wpcols
 
                call solutions(i) % ci_vec % init_CV(int(oldsolution % ci_vec % mat_dimen_r), &
                                                     int(oldsolution % ci_vec % mat_dimen_c))
                call solutions(i) % ci_vec % redistribute(oldsolution % ci_vec)
 
            end if
        end do
 
    end subroutine redistribute_solutions
 
 
    subroutine cdenprop_properties (SCATCI_input, solutions, properties_all)
 
        use class_molecular_properties_data, only: molecular_properties_data
        use class_namelists,                 only: cdenprop_namelists
        use consts_mpi_ci,                   only: pass_to_cdenprop, no_diagonalization, no_ci_target
        use cdenprop_procs,                  only: cdenprop_drv
        use cdenprop_defs,                   only: ir_max
        use mpi_gbl,                         only: myrank, master, mpi_reduceall_sum_cfp
        use options_module,                  only: optionsset
        use solutionhandler_module,          only: solutionhandler
 
        type(optionsset),                   intent(in)    :: SCATCI_input
        type(solutionhandler), allocatable, intent(inout) :: solutions(:)
        type(molecular_properties_data),    intent(inout) :: properties_all
 
        type(cdenprop_namelists) :: namelist_ij
        type(civect)             :: oldblock
 
        integer :: i, j, ntgsym, ninputs, maxpole, nblock, mxstat
 
        ninputs = size(scatci_input % opts)
        maxpole = 2
 
        ! for now, skip evaluation of properties if RMT data not required
        if (.not. scatci_input % write_dip .and. .not. scatci_input % write_rmt) return
 
        ! prepare free slots for the distributed dense property matrices
        call properties_all % clean
        call properties_all % preallocate_property_blocks(ninputs**2 * (maxpole + 1)**2)
 
        ! calculate properties for all pairs of symmetries
        sym_i_loop: do i = 1, ninputs
            if (iand(scatci_input % opts(i) % vector_storage_method, pass_to_cdenprop) /= 0 .and. &
                scatci_input % opts(i) % diagonalization_flag /= no_diagonalization) then
 
                sym_j_loop: do j = i, ninputs
                    if (iand(scatci_input % opts(j) % vector_storage_method, pass_to_cdenprop) /= 0 .and. &
                        scatci_input % opts(j) % diagonalization_flag /= no_diagonalization) then
 
                        ! only calculate self-overlaps when required
                        if (i == j .and. .not. scatci_input % all_props) cycle
 
                        ! virtual CDENPROP namelist
                        call namelist_ij % init(2)
 
                        ! assign values obtained from SCATCI namelists
                        namelist_ij % lucsf(1:2)   = (/ scatci_input % opts(i) % megul, &
                                                        scatci_input % opts(j) % megul /)
                        namelist_ij % nciset(1:2)  = (/ scatci_input % opts(i) % output_ci_set_number, &
                                                        scatci_input % opts(j) % output_ci_set_number /)
                        namelist_ij % nstat        = 0      ! use all states
                        namelist_ij % lupropw      = 0      ! no output needed now, will be done later
                        namelist_ij % lucivec      = 0      ! no CI units needed, will be given as arguments
                        namelist_ij % ukrmolp_ints = scatci_input % opts(i) % use_UKRMOL_integrals
                        namelist_ij % luprop       = merge(scatci_input % opts(i) % integral_unit, &
                                                           scatci_input % luprop, &
                                                           scatci_input % luprop < 0)
 
                        ntgsym = scatci_input % opts(i) % num_target_sym
 
                        ! following depend on whether this is target or scattering calculation
                        if (scatci_input % opts(i) % ci_target_flag == no_ci_target) then
                            namelist_ij % numtgt(1:ir_max) = 0
                            namelist_ij % max_multipole    = 2  ! dipole + quadrupole
                        else
                            namelist_ij % max_multipole    = 1  ! dipole only
                            namelist_ij % numtgt(1:ir_max) = 0
                            namelist_ij % numtgt(1:ntgsym) = scatci_input % opts(i) % num_target_state_sym(1:ntgsym)
                        end if
 
                        ! use BLACS context of the first vector for communication (they are all of the same shape anyway)
                        solutions(j) % ci_vec % blacs_context   = solutions(i) % ci_vec % blacs_context
                        solutions(j) % ci_vec % descr_CV_mat(2) = solutions(i) % ci_vec % descr_CV_mat(2)
 
                        ! get the properties in CDENPROP
                        call cdenprop_drv(solutions(i) % ci_vec, solutions(j) % ci_vec, namelist_ij, properties_all)
 
                        ! cleanup
                        call namelist_ij % dealloc
 
                    end if
                end do sym_j_loop
 
            end if
        end do sym_i_loop
 
        ! prefer SWINTERF-compatible format of the property file (so that it can be, e.g., read back to CDENPROP_ALL)
        properties_all % swintf_format = .true.
 
        ! write calculated properties
        if (scatci_input % write_dip .and. properties_all % no_symmetries > 0) then
            if (scatci_input % all_props) call properties_all % sort_by_energy
            call properties_all % write_properties(624, scatci_input % opts(1) % use_UKRMOL_integrals)
        end if
 
        ! inflate all property matrices to uniform size for easy output for RMT, which wants it that way
        if (scatci_input % write_rmt) then
            nblock = properties_all % no_blocks
            mxstat = max(maxval(properties_all % dense_blocks(1:nblock) % mat_dimen_r), &
                         maxval(properties_all % dense_blocks(1:nblock) % mat_dimen_c))
            do i = 1, nblock
                associate(d => properties_all % dense_blocks(i))
                    if (d % mat_dimen_r /= mxstat .or. d % mat_dimen_c /= mxstat) then
                        oldblock = d
                        call d % init_CV(mxstat, mxstat)
                        d % CV = 0
                        call d % redistribute(oldblock)
                        call oldblock % final_CV
                    end if
                end associate
            end do
        end if
 
    end subroutine cdenprop_properties
 
 
    subroutine outer_interface (SCATCI_input, solutions, inner_properties)
 
        use class_molecular_properties_data, only: molecular_properties_data
        use consts_mpi_ci,                   only: pass_to_cdenprop, no_diagonalization
        use options_module,                  only: optionsset
        use solutionhandler_module,          only: solutionhandler
 
        type(optionsset),                   intent(in) :: SCATCI_input
        type(solutionhandler), allocatable, intent(in) :: solutions(:)
        type(molecular_properties_data),    intent(in) :: inner_properties
 
        type(molecular_properties_data)   :: target_properties
        type(outerinterface), allocatable :: intf(:)
 
        integer :: i, nsym, nset
 
        nset = 0                          ! number of sets in the output files fort.10 and fort.21
        nsym = size(scatci_input % opts)  ! number of inputs (symmetries) that we dealt with
 
        if (.not. scatci_input % write_amp .and. .not. scatci_input % write_rmt) return
 
        allocate (intf(nsym))
 
        call target_properties % read_properties(24, 1)
 
        ! run SWITERF for all symmetries
        do i = 1, nsym
            call intf(i) % init(scatci_input)
            if (iand(scatci_input % opts(i) % vector_storage_method, pass_to_cdenprop) /= 0 .and. &
                scatci_input % opts(i) % diagonalization_flag /= no_diagonalization) then
                nset = nset + 1
                call intf(i) % setup_amplitudes(scatci_input % opts(i))
                call intf(i) % extract_data(scatci_input % opts(i), target_properties, solutions(i))
                if (scatci_input % write_amp) then
                    call intf(i) % write_data(i, scatci_input, target_properties, solutions(i))
                end if
            end if
        end do
 
        if (scatci_input % write_rmt) then
            call write_rmt_data(scatci_input, inner_properties, target_properties, solutions, intf)
        end if
 
    end subroutine outer_interface
 
 
    subroutine init_outer_interface (this, input)
 
        use mpi_gbl,        only: mpi_xermsg
        use options_module, only: optionsset
 
        class(outerinterface), intent(inout) :: this
        type(optionsset),      intent(in)    :: input
 
        integer :: i, ierr
 
        this % ntarg  = 0
        this % nchan  = 0
        this % ismax  = 0
        this % maxnmo = 0
 
        this % nfdm    = input % nfdm
        this % rmatr   = input % rmatr
 
        allocate (this % r_points(this % nfdm + 1))
 
        do i = 1, this % nfdm + 1
            this % r_points(i) = this % rmatr - (this % nfdm + 1 - i) * input % delta_r
        end do
 
        allocate (this % irrchl(this % maxchan), &
                  this % ichl(this % maxchan),   &
                  this % lchl(this % maxchan),   &
                  this % mchl(this % maxchan),   &
                  this % qchl(this % maxchan),   &
                  this % echl(this % maxchan), stat = ierr)
 
        if (ierr /= 0) then
            call mpi_xermsg('OuterInterface', 'init_outer_interface', 'Failed to allocate memory for channel data.', 1, 1)
        end if
 
        if (allocated(input % idtarg)) then
            this % ntarg  = size(input % idtarg)
            this % idtarg = input % idtarg
        end if
 
    end subroutine init_outer_interface
 
 
    subroutine deinit_outer_interface (this)
 
        class(outerinterface), intent(inout) :: this
 
        integer :: i
 
        ! release the no-longer-needed BLACS context
        if (allocated(this % wamp)) then
            do i = 1, size(this % wamp)
                call this % wamp(i) % final_CV
#if defined(usempi) && defined(scalapack)
                if (i == 1 .and. this % wamp(this % nfdm + 1) % blacs_context >= 0) then
                    call blacs_gridexit(this % wamp(this % nfdm + 1) % blacs_context)
                end if
#endif
            end do
            deallocate (this % wamp)
        end if
 
        if (allocated(this % irrchl))   deallocate (this % irrchl)
        if (allocated(this % ichl))     deallocate (this % ichl)
        if (allocated(this % lchl))     deallocate (this % lchl)
        if (allocated(this % mchl))     deallocate (this % mchl)
        if (allocated(this % qchl))     deallocate (this % qchl)
        if (allocated(this % echl))     deallocate (this % echl)
        if (allocated(this % a))        deallocate (this % a)
        if (allocated(this % r_points)) deallocate (this % r_points)
 
        this % nchan  = 0
        this % ismax  = 0
        this % maxnmo = 0
        this % nfdm   = 0
 
    end subroutine deinit_outer_interface
 
 
    subroutine setup_amplitudes (this, opts)
 
        use options_module,       only: options
        use ukrmol_interface_gbl, only: read_ukrmolp_basis, eval_amplitudes
 
        class(outerinterface), intent(inout) :: this
        type(options),         intent(in)    :: opts
 
        if (opts % use_UKRMOL_integrals) then
            call read_ukrmolp_basis(opts % integral_unit)
            call eval_amplitudes(this % rmatr, .true.)
        end if
 
    end subroutine setup_amplitudes
 
 
    subroutine extract_data (this, opts, prop, solution)
 
        use class_molecular_properties_data, only: molecular_properties_data
        use options_module,                  only: options
        use solutionhandler_module,          only: solutionhandler
 
        class(outerinterface),           intent(inout) :: this
        type(molecular_properties_data), intent(in)    :: prop
        type(options),                   intent(in)    :: opts
        type(solutionhandler),           intent(in)    :: solution
 
        integer,  parameter :: ismax  = 2
        real(wp), parameter :: alpha0 = 0
        real(wp), parameter :: alpha2 = 0
 
        logical :: use_pol
        integer :: lamax
 
        lamax = ismax
        use_pol = .false.  ! if set to .true., then channel_couplings will take into account the polarizability
 
        if (alpha0 /= 0 .or. alpha2 /= 0) then
           use_pol = .true.
           lamax = max(ismax, 3)  ! polarizability (1/r^4) corresponds to lambda = 3
        end if
 
        this % ismax = lamax
 
        call this % get_channel_info(opts, prop)
        call this % get_channel_couplings(lamax, prop % no_states, prop, alpha0, alpha2, use_pol)
        call this % get_boundary_data(opts, prop, solution)
 
        ! redefine magnetic quantum number by multiplying in the projectile charge
        this % mchl(1:this % nchan) = this % mchl(1:this % nchan) * this % qchl(1:this % nchan)
 
    end subroutine extract_data
 
 
    subroutine write_data (this, i, opts, prop, solution)
 
        use class_molecular_properties_data, only: molecular_properties_data
        use mpi_gbl,                         only: myrank, master
        use options_module,                  only: optionsset
        use solutionhandler_module,          only: solutionhandler
 
        class(outerinterface),           intent(in) :: this
        type(molecular_properties_data), intent(in) :: prop
        type(optionsset),                intent(in) :: opts
        type(solutionhandler),           intent(in) :: solution
        integer,                         intent(in) :: i
 
        if (myrank == master) then
            call this % write_channel_info(i, opts % opts(i), opts % cform, prop)
            call this % write_boundary_data(i, opts % opts(i), opts % rform, prop, solution)
        end if
 
    end subroutine write_data
 
 
    subroutine get_channel_info (this, opts, prop)
 
        use algorithms,                      only: findloc
        use class_molecular_properties_data, only: molecular_properties_data
        use const_gbl,                       only: stdout
        use global_utils,                    only: mprod
        use mpi_gbl,                         only: mpi_xermsg
        use options_module,                  only: options
        use precisn,                         only: wp
        use ukrmol_interface_gbl,            only: ukp_preamp
 
        class(outerinterface),           intent(inout) :: this
        type(options),                   intent(in)    :: opts
        type(molecular_properties_data), intent(in)    :: prop
 
        integer  :: i, j, k, mgvn, spin, tgt_stat, tgt_symm, tgt_mgvn, tgt_spin, nch, prj_mgvn, congen_tgt_id, denprop_tgt_id, isym
        real(wp) :: tgt_enrg, ebase
 
        ebase  = prop % energies(1)
        mgvn   = opts % lambda
        spin   = nint(2 * opts % spin + 1)
 
        this % maxnmo = sum(opts % num_orbitals_sym)
        this % nchan  = 0
 
        ! construct default IDTARG, i.e. denprop (energy-order) index for each congen (symmetry-order) target state
        if (.not. allocated(this % idtarg)) then
            this % ntarg = opts % total_num_target_states
            allocate (this % idtarg(this % ntarg))
            congen_tgt_id = 0
            congen_symm_loop: do j = 1, size(opts % num_target_state_sym)
                congen_targ_loop: do k = 1, opts % num_target_state_sym(j)
                    congen_tgt_id = congen_tgt_id + 1
                    denprop_tgt_id = 0
                    tgt_mgvn = opts % target_spatial(congen_tgt_id)
                    tgt_spin = opts % target_multiplicity(congen_tgt_id)
                    denprop_targ_loop: do i = 1, this % ntarg
                        isym = prop % energies_index(2, i)
                        if (prop % symmetries_index(1, isym) == tgt_mgvn .and. &
                            prop % symmetries_index(2, isym) == tgt_spin) then
                            denprop_tgt_id = denprop_tgt_id + 1
                        end if
                        if (denprop_tgt_id == k) then
                            this % idtarg(congen_tgt_id) = i
                            exit denprop_targ_loop
                        end if
                    end do denprop_targ_loop
                end do congen_targ_loop
            end do congen_symm_loop
        end if
 
        write (stdout, '(/,1x,"Target state indices sorted by energy (IDTARG): ")')
        write (stdout, '(  1x,20I4)') this % idtarg
 
        do i = 1, prop % no_states
 
            ! skip states not referenced in IDTARG
            if (findloc(this % idtarg, i, 1) == 0) cycle
 
            ! retrieve target data from the properties structure
            tgt_stat = prop % energies_index(1, i)
            tgt_symm = prop % energies_index(2, i)
            tgt_mgvn = prop % symmetries_index(1, tgt_symm)
            tgt_spin = prop % symmetries_index(2, tgt_symm)
            tgt_enrg = prop % energies(i)
 
            ! angular momentum composition rule -- must allow for projectile spin
            if (abs(spin - tgt_spin) /= 1) cycle
 
            ! get projectile total symmetry IRR from the multiplication table
            prj_mgvn = mprod(tgt_mgvn + 1, mgvn + 1, 0, stdout)
 
            ! retrieve channel descriptions from the integral library
            if (opts % use_UKRMOL_integrals) then
                call ukp_preamp(prj_mgvn, this % nchan + 1, this % lchl, this % mchl, this % qchl, nch, this % maxnmo)
            else
                call mpi_xermsg('Postprocessing_module', 'get_channel_info', &
                                'Interface for SWINTERF integrals not implemented.', 1, 1)
            end if
 
            ! fill in the remaining channel identifiers
            do j = 1, nch
                this % ichl(this % nchan + j) = tgt_stat
                this % echl(this % nchan + j) = 2 * (tgt_enrg - ebase)
                this % irrchl(this % nchan + j) = prj_mgvn
            end do
 
            ! update total number of channels
            this % nchan = this % nchan + nch
 
        end do
 
        write (stdout, '(/," Channel Target  l  m  q       Irr  Energy   Energy(eV)")')
 
        do i = 1, this % nchan
            write (stdout, '(I5,I8,I5,2I3,2x,A4,I3,2F10.6)') i, this % ichl(i), this % lchl(i), this % mchl(i), this % qchl(i), &
                                                             'N/A', this % irrchl(i), this % echl(i), this % echl(i) / 0.073500d0
        end do
 
    end subroutine get_channel_info
 
 
    subroutine get_boundary_data (this, opts, prop, solution)
 
        use class_molecular_properties_data, only: molecular_properties_data
        use const_gbl,                       only: stdout
        use global_utils,                    only: mprod
        use mpi_gbl,                         only: master, myrank
        use options_module,                  only: options
        use solutionhandler_module,          only: solutionhandler
        use ukrmol_interface_gbl,            only: eval_amplitudes, ukp_readamp
 
        class(outerinterface),           intent(inout)  :: this
        type(options),                   intent(in)     :: opts
        type(molecular_properties_data), intent(in)     :: prop
        type(solutionhandler),           intent(in)     :: solution
 
        type(civect) :: orb_amps, dist_orb_amps, sol_amps
 
        integer,  allocatable :: mcont(:), idchl(:)
        real(wp), allocatable :: ampls(:)
 
        integer :: i, j, k, ir, ncontmo(8), iprnt = 0, nchan, nstat, ierr, nnoncontmo, norbs, morbs, ncnt, offs, nvo(8) = 0
 
        nchan = this % nchan
        nstat = solution % ci_vec % nstat
        morbs = maxval(opts % num_orbitals_sym)  ! safe upper bound on maximal number of continuum orbitals per channel...
        norbs = dot_product(opts % num_target_state_sym, opts % num_continuum_orbitals_target)  ! ... and total
 
        allocate(ampls(morbs), idchl(nchan))
 
        ! set up energy order channel map
        k = 0
        do j = 1, this % ntarg
            do i = 1, nchan
                if (this % ichl(i) == this % idtarg(j)) then
                    k = k + 1
                    idchl(k) = i
                end if
            end do
        end do
 
        write (stdout, '(/,1x,"Channel target state energy map (IDCHL): ")')
        write (stdout, '(  1x,20I4)') idchl
 
        ! set up storage for the raw boundary amplitudes (master-only ScaLAPACK matrix)
        orb_amps % nprow = 1
        orb_amps % npcol = 1
        orb_amps % blacs_context   = -1
        orb_amps % CV_is_scalapack = solution % ci_vec % CV_is_scalapack
#if defined(usempi) && defined(scalapack)
        if (orb_amps % CV_is_scalapack) then
            call blacs_get(-1_blasint, 0_blasint, orb_amps % blacs_context)
            call blacs_gridinit(orb_amps % blacs_context, 'R', 1_blasint, 1_blasint)
        end if
#endif
        call orb_amps % init_CV(nchan, norbs)
 
        ! set up storage for the evaluated boundary amplitudes
        if (.not. allocated(this % wamp)) then
            allocate (this % wamp(this % nfdm + 1))
        end if
 
        ! most of them are distributed...
        this % wamp(:) % nprow = solution % ci_vec % nprow
        this % wamp(:) % npcol = solution % ci_vec % npcol
        this % wamp(:) % blacs_context = solution % ci_vec % blacs_context
        this % wamp(:) % CV_is_scalapack = solution % ci_vec % CV_is_scalapack
 
        ! ...but the last is master-only (needed for write_boundary_data)
        this % wamp(this % nfdm + 1) % nprow = 1
        this % wamp(this % nfdm + 1) % npcol = 1
        this % wamp(this % nfdm + 1) % blacs_context = orb_amps % blacs_context
 
        ! allocate the matrices
        do ir = 1, this % nfdm + 1
            call this % wamp(ir) % init_CV(nchan, nstat)
        end do
 
        ! for all evaluation radii
        do ir = 1, this % nfdm + 1
 
            ! get raw boundary amplitudes of the continuum orbitals
            if (myrank == master) then
                ! retrieve orbital amplitudes from the integral library
                call eval_amplitudes(this % r_points(ir), .false.)
                call ukp_readamp(orb_amps % CV, this % nchan, this % irrchl, this % lchl, this % mchl, this % qchl, ncontmo, &
                                 opts % lambda_continuum_orbitals_target, iprnt)
 
                ! re-align the orbital amplitude data to their absolute positions in the list of all virt+cont orbitals:
                !
                !              1st target used in channels...........2nd target used in channels...........etc....
                !    orb_amps: [virtual orbitals][continuum orbitals][virtual orbitals][continuum orbitals]etc...
                offs = 0
                do i = 1, nchan
                    j = idchl(i)
                    if (i > 1) then
                        if (this % ichl(j) /= this % ichl(idchl(i - 1))) then
                            offs = offs + ncnt + nvo(this % irrchl(j))
                        end if
                    end if
                    ncnt = ncontmo(this % irrchl(j))
                    ampls(1:ncnt) = orb_amps % CV(j, 1:ncnt)
                    orb_amps % CV(j, 1:norbs) = 0
                    orb_amps % CV(j, offs + 1:offs + ncnt) = ampls(1:ncnt)
                end do
            end if
 
            ! multiply the raw boundary amplitudes with the eigenvectors
            if (.not. solution % ci_vec % CV_is_scalapack) then
                call this % wamp(ir) % A_B_matmul(orb_amps, solution % ci_vec, 'N', 'N')
            else
                ! distribute raw boundary amplitudes over solution context
                dist_orb_amps % nprow = solution % ci_vec % nprow
                dist_orb_amps % npcol = solution % ci_vec % npcol
                dist_orb_amps % blacs_context   = solution % ci_vec % blacs_context
                dist_orb_amps % CV_is_scalapack = solution % ci_vec % CV_is_scalapack
                call dist_orb_amps % init_CV(nchan, norbs)
                call dist_orb_amps % redistribute(orb_amps)
 
                ! multiply distributed matrices
                if (ir < this % nfdm + 1) then
                    call this % wamp(ir) % A_B_matmul(dist_orb_amps, solution % ci_vec, 'N', 'N')
                else
                    ! this is needed by write_boundary_data, which cannot use distributed array
                    sol_amps % nprow = solution % ci_vec % nprow
                    sol_amps % npcol = solution % ci_vec % npcol
                    sol_amps % blacs_context   = solution % ci_vec % blacs_context
                    sol_amps % CV_is_scalapack = solution % ci_vec % CV_is_scalapack
                    call sol_amps % init_CV(nchan, nstat)
                    call sol_amps % A_B_matmul(dist_orb_amps, solution % ci_vec, 'N', 'N')
                    call this % wamp(ir) % redistribute(sol_amps, solution % ci_vec % blacs_context)
                end if
            end if
 
            ! change normalization convention to the RMT one
            this % wamp(ir) % CV = sqrt(2.0_wp) * this % wamp(ir) % CV
 
        end do
 
    end subroutine get_boundary_data
 
 
    subroutine write_channel_info (this, nchset, opts, cform, prop)
 
        use algorithms,                      only: findloc
        use class_molecular_properties_data, only: molecular_properties_data
        use consts,                          only: amu
        use consts_mpi_ci,                   only: symtype_d2h
        use options_module,                  only: options
 
        class(outerinterface), intent(in) :: this
        integer,               intent(in) :: nchset
        type(options),         intent(in) :: opts
        character(len=1),      intent(in) :: cform
 
        type(molecular_properties_data), intent(in) :: prop
 
        integer, parameter :: luchan = 10, keych = 10
 
        integer :: i, k, nvd, nrec, ninfo, ndata, nvib, ndis, stot, gutot, ntarg, ivtarg(1), iv(1), ion, mtarg, starg, gutarg, nnuc
        integer :: isymm
        real(wp) :: r, rmass, xmass, etarg
 
        nnuc  = prop % no_nuclei
        nvib  = 0
        ndis  = 0
        gutot = symtype_d2h  ! to indicate polyatomic code
        ion   = nint(sum(prop % charge(1:nnuc))) - (opts % num_electrons - 1)
        r     = 0
        ntarg = this % ntarg
        stot  = nint(2 * opts % spin + 1)
        rmass = sum(1 / prop % mass(1:nnuc), prop % mass(1:nnuc) > 0)
 
        if (rmass > 0) then
            rmass = amu / rmass
        end if
 
        nvd   = nvib + ndis
        nrec  = 3 + ntarg + this % nchan + nvd
        ninfo = 1
        ndata = nrec - ninfo
 
        if (cform == 'F') then
            write (luchan, '(16I5)') keych, nchset, nrec, ninfo, ndata
            write (luchan, '(A80)') opts % diag_name
            write (luchan, '(16I5)') ntarg, nvib, ndis, this % nchan
            write (luchan, '(4I5,2D20.12)') opts % lambda, stot, gutot, ion, r, rmass
 
            k = 0
            do i = 1, prop % no_states
                if (findloc(this % idtarg, i, 1) /= 0) then
                    k = k + 1
                    isymm  = prop % energies_index(2, i)
                    mtarg  = prop % symmetries_index(1, isymm)
                    starg  = prop % symmetries_index(2, isymm)
                    gutarg = 0
                    etarg  = prop % energies(i)
                    write (luchan, '(4I5,2D20.12)') k, mtarg, starg, gutarg, etarg
                end if
            end do
            do i = 1, this % nchan
                write (luchan, '(4I5,2D20.12)') i, this % ichl(i), this % lchl(i), this % mchl(i), this % echl(i)
            end do
            do i = 1, nvd
                write (luchan, '(16I5)') i, ivtarg(i), iv(i)
            end do
        else
            write (luchan) keych, nchset, nrec, ninfo, ndata
            write (luchan) opts % diag_name
            write (luchan) ntarg, nvib, ndis, this % nchan
            write (luchan) opts % lambda, stot, gutot, ion, r, rmass
 
            k = 0
            do i = 1, prop % no_states
                if (findloc(this % idtarg, i, 1) > 0) then
                    k = k + 1
                    isymm  = prop % energies_index(2, i)
                    mtarg  = prop % symmetries_index(1, isymm)
                    starg  = prop % symmetries_index(2, isymm)
                    gutarg = 0
                    etarg  = prop % energies(i)
                    write (luchan) k, mtarg, starg, gutarg, etarg
                end if
            end do
            do i = 1, this % nchan
                write (luchan) i, this % ichl(i), this % lchl(i), this % mchl(i), this % echl(i)
            end do
            do i = 1, nvd
                write (luchan) i, ivtarg(i), iv(i)
            end do
        end if
 
    end subroutine write_channel_info
 
 
    subroutine write_boundary_data (this, nrmset, opts, rform, prop, solution)
 
        use cdenprop_defs,                   only: amass
        use class_molecular_properties_data, only: molecular_properties_data
        use consts,                          only: amu
        use consts_mpi_ci,                   only: symtype_d2h
        use options_module,                  only: options
        use solutionhandler_module,          only: solutionhandler
 
        class(outerinterface),           intent(in) :: this
        integer,                         intent(in) :: nrmset
        type(options),                   intent(in) :: opts
        character(len=1),                intent(in) :: rform
        type(molecular_properties_data), intent(in) :: prop
        type(solutionhandler),           intent(in) :: solution
 
        integer, parameter :: keyrm = 11
        integer, parameter :: lurmt = 21
 
        integer  :: i, j, k, nchsq, nrec, ninfo, ndata, ntarg, nvib, ndis, nchan, stot, gutot, ion, ibut, iex, nocsf, npole, nnuc
        integer  :: ismax, nstat
        real(wp) :: r, rmass, rmatr
 
        nnuc  = prop % no_nuclei
        iex   = 0
        ibut  = 0
        nocsf = 0
        npole = 0
        ismax = 2
        rmatr = this % rmatr
        ntarg = this % ntarg
        nstat = solution % ci_vec % nstat
        nvib  = 0
        ndis  = 0
        nchan = this % nchan
        stot  = nint(2 * opts % spin + 1)
        gutot = symtype_d2h  ! to indicate polyatomic code
        ion   = nint(sum(prop % charge(1:nnuc))) - (opts % num_electrons - 1)
        r     = 0
        nchsq = nchan * (nchan + 1) / 2
        nrec  = 4
        rmass = sum(1 / prop % mass(1:nnuc), prop % mass(1:nnuc) > 0)
 
        if (rmass > 0) then
            rmass = amu / rmass
        end if
 
        if (rform == 'F') then
            if (ismax > 0) nrec = nrec + (nchsq * ismax + 3) / 4
            nrec = nrec + (nstat + 3) / 4
            nrec = nrec + (nchan * nstat + 3) / 4
            if (npole > 0) nrec = nrec + (nocsf * npole + 3) / 4
            if (ibut > 0) nrec = nrec + (3 * nchan + 3) / 4
            if (abs(ibut) > 1) nrec = nrec + 1 + (nchan + 3) / 4 + (iex + 3) / 4 + (nchan * iex + 3) / 4
        else
            nrec = nrec + 2
            if (ismax > 0) nrec = nrec + 1
            if (npole > 0) nrec = nrec + 1
            if (ibut > 0) nrec = nrec + 1
            if (abs(ibut) > 1) nrec = nrec + 4
        end if
 
        ninfo = 1
        ndata = nrec - ninfo
 
        if (rform == 'F') then
            write (lurmt, '(10I7)') keyrm, nrmset, nrec, ninfo, ndata
            write (lurmt, '(A80)') opts % diag_name
            write (lurmt, '(16I5)') ntarg, nvib, ndis, nchan
            write (lurmt, '(4I5,2D20.13)') opts % lambda, stot, gutot, ion, r, rmass
            write (lurmt, '(4I5,2D20.13)') ismax, nstat, npole, ibut, rmatr
           !if (abs(ibut) > 1) write(lurmt, '(i10,d20.13,i5)') nocsf, ezero, iex
            if (ismax > 0) write(lurmt, '(4D20.13)') ((this % a(i,k),i=1,nchsq),k=1,ismax)
            write (lurmt, '(4D20.13)') (solution % ci_vec % ei(i) + solution % core_energy,i=1,nstat)
            write (lurmt, '(4D20.13)') ((this % wamp(this % nfdm + 1) % cv(i,j)/sqrt(2.0_wp),i=1,nchan),j=1,nstat)
            if (npole > 0) write (lurmt, '(4D20.13)') ((solution % ci_vec % cv(i,j),i=1,nocsf),j=1,npole)   ! TODO: all ranks need to write their part
           !if (ibut > 0) write (lurmt, '(4D20.13)') ((bcoef(i,j),i=1,3),j=1,nchan)
           !if (abs(ibut) > 1) write (lurmt, '(4D20.13)') (sfac(j),j=1,nchan)
           !if (abs(ibut) > 1) write (lurmt, '(4D20.13)') (ecex(j),j=1,iex)
           !if (abs(ibut) > 1) write (lurmt, '(4D20.13)') ((rcex(i,j),i=1,nchan),j=1,iex)
        else
            write(lurmt) keyrm, nrmset, nrec, ninfo, ndata
            write(lurmt) opts % diag_name
            write(lurmt) ntarg, nvib, ndis, nchan
            write(lurmt) opts % lambda, stot, gutot, ion, r, rmass
            write(lurmt) ismax, nstat, npole, ibut, rmatr
           !if (abs(ibut) > 1) write(lurmt) nocsf, ezero, iex
            if (ismax > 0) write(lurmt) ((this % a(i,k),i=1,nchsq),k=1,ismax)
            write (lurmt) (solution % ci_vec % ei(i) + solution % core_energy,i=1,nstat)
            write (lurmt) ((this % wamp(this % nfdm + 1) % cv(i,j)/sqrt(2.0_wp),i=1,nchan),j=1,nstat)
            if (npole > 0) write (lurmt) ((solution % ci_vec % cv(i,j),i=1,nocsf),j=1,npole)                ! TODO: all ranks need to write their part
           !if (ibut > 0) write (lurmt) ((bcoef(i,j),i=1,3),j=1,nchan)
           !if (abs(ibut) > 1) write (lurmt) (sfac(j),j=1,nchan)
           !if (abs(ibut) > 1) write (lurmt) (ecex(j),j=1,iex)
           !if (abs(ibut) > 1) write (lurmt) ((rcex(i,j),i=1,nchan),j=1,iex)
        end if
 
    end subroutine write_boundary_data
 
 
    subroutine get_channel_couplings (this, ismax, ntarg, target_properties, alpha0, alpha2, use_pol)
 
        use coupling_obj_gbl,                only: couplings_type
        use class_molecular_properties_data, only: molecular_properties_data
        use mpi_gbl,                         only: mpi_xermsg
        use precisn,                         only: wp
 
        class(outerinterface),           intent(inout) :: this
        type(molecular_properties_data), intent(in)    :: target_properties
 
        integer,  intent(in)    :: ismax, ntarg
        real(wp), intent(in)    :: alpha0, alpha2
        logical,  intent(in)    :: use_pol
 
        real(wp), parameter :: fourpi  = 12.5663706143591729538505735331180115_wp
        real(wp), parameter :: roneh   =  0.70710678118654752440084436210484903_wp !sqrt(0.5)
        real(wp), parameter :: rothree =  0.57735026918962576450914878050195745_wp !1/sqrt(3)
 
        ! The Gaunt coefficients for the real spherical harmonics defined below are needed to express the xx, yy, zz angular
        ! behaviour in terms of the real spherical harmonics.
 
        ! x^2 = x*x ~ X_{11}*X_{11} = x2_X00*X_{00} + x2_X20*X_{20} + x2_X22*X_{22}
        real(wp), parameter :: x2_X00 =  0.282094791773878e+00_wp
        real(wp), parameter :: x2_X20 = -0.126156626101008e+00_wp
        real(wp), parameter :: x2_X22 =  0.218509686118416e+00_wp
        ! y^2 = y*y ~ X_{1-1}*X_{1-1} = y2_X00*X_{00} + y2_X20*X_{20} + y2_X22*X_{22}
        real(wp), parameter :: y2_X00 =  0.282094791773878e+00_wp
        real(wp), parameter :: y2_X20 = -0.126156626101008e+00_wp
        real(wp), parameter :: y2_X22 = -0.218509686118416e+00_wp
        ! z^2 = z*z ~ X_{10}*X_{10} = z2_X00*X_{00} + z2_X20*X_{20}
        real(wp), parameter :: z2_X00 =  0.282094791773878e+00_wp
        real(wp), parameter :: z2_X20 =  0.252313252202016e+00_wp
 
        type(couplings_type) :: couplings
 
        real(wp), allocatable :: prop(:,:,:)
 
        real(wp) :: cpl, sph_cpl, fac
        integer  :: l1, m1, q1, l2, m2, q2, no_cpl, lqt, isq, it1, it2
        integer  :: ch_1, ch_2, lambda, mlambda, lmin, lmax
        logical  :: use_alpha2
 
        if (use_pol .and. ntarg > 1) then
            write (*, '("WARNING: adding polarization potential while more than one target state is present in the outer region.")')
        end if
 
        ! get a convenient dense matrix set with the properties
        call target_properties % decompress(1, ismax, prop)
 
        ! find size parameters based on polarizability setup
        use_alpha2 = (use_pol .and. alpha2 /= 0.0_wp)
        lmax = merge(max(ismax, 3), ismax, use_pol)     ! maximal angular momentum transfer (polarizability corresponds to lambda=3)
        no_cpl = this % nchan * (this % nchan + 1) / 2  ! total number of unique combinations of the scattering channels
 
        ! allocate memory for the coupling coefficients
        if (allocated(this % a)) deallocate (this % a)
        allocate (this % a(no_cpl, lmax))               ! indices (:, 0) are never needed, though
        this % a(:,:) = 0.0_wp
 
        do ch_1 = 1, this % nchan
 
            l1 = this % lchl(ch_1)  ! the l, |m| values correspond to the l,|m| values of the real spherical harmonics
            m1 = this % mchl(ch_1)
            q1 = this % qchl(ch_1)  ! q is sign(m) or 0 if m=0.
            m1 = m1 * q1            ! determine m
            it1 = this % ichl(ch_1) ! sequence number of the target state corresponding to this channel
 
            do ch_2 = 1, ch_1
 
                l2 = this % lchl(ch_2)  ! the l, |m| values correspond to the l,|m| values of the real spherical harmonics
                m2 = this % mchl(ch_2)
                q2 = this % qchl(ch_2)  ! q is sign(m) or 0 if m=0.
                m2 = m2 * q2            ! determine m
                it2 = this % ichl(ch_2) ! sequence number of the target state corresponding to this channel
 
                ! use the selection rules for the real spherical harmonics to determine the range of lambda values
                ! which may give non-zero real Gaunt coefficients
                call couplings % bounds_rg(l1, l2, m1, m2, lmin, lmax)
 
                !don't include potentials with lambda > ismax
                if (lmax > ismax) lmax = ismax
 
                ! lambda = 0 would correspond to the monopole contribution, i.e. taking into account the total (perhaps nonzero)
                ! molecular charge. This is taken into account separately in RSOLVE.
                if (lmin == 0) lmin = 2
 
                ! linear index corresponding to the current combination of the channels (ch_1,ch_2).
                lqt = ch_1 * (ch_1 - 1) / 2 + ch_2
 
                ! Loop over the multipole moments which may be non-zero: we loop in steps of two since the selection rules
                ! for the r.s.h. imply that only if the sum of all L values is even the coupling may be non-zero.
                do lambda = lmin, lmax, 2
 
                    ! see the formula for V_{ij}: effectively this converts the inner region property from solid harmonic
                    ! normalization to spherical harmonic normalization as needed by the Leg. expansion.
                    fac = sqrt(fourpi / (2 * lambda + 1.0_wp))
 
                    do mlambda = -lambda, lambda
 
                        ! rgaunt: the coupling coefficient for the real spherical harmonics (l1,m1,l2,m2,lambda,mlambda);
                        ! The factor 2.0_wp converts the units of the inner region
                        ! properties from Hartree to Rydberg since that's the energy
                        ! unit used in the outer region.
                        cpl = 2.0_wp * fac * couplings % rgaunt(l1, lambda, l2, m1, mlambda, m2)
 
                        ! linear index corresponding to the current (lambda,mlambda) values. lambda*lambda is the number
                        ! of all (lambda1,mlambda1) combinations for lambda1 = lambda-1.
                        isq = lambda * lambda + lambda + mlambda
 
                        ! increment the value of the coefficient for the multipole coupling potential of order lambda between
                        ! the target states it1 and it2
                        this % a(lqt, lambda) = this % a(lqt, lambda) + cpl * prop(it1, it2, isq)
 
                    end do !mlambda
 
                end do !lambda
 
                ! add polarizabilities for the ground state channels, but only if use_pol == .true.
                if (use_pol .and. it1 == 1 .and. it2 == 1) then
 
                    ! the radial dependence of the polarization potential is r^{-4} which corresponds to lambda = 3
                    lambda = 3
 
                    ! obviously, the spherical polarizability (l=0) couples only the channels with identical angular behaviour
                    if (l1 == l2 .and. m1 == m2) then
                        sph_cpl = 1.0_wp
                        write(*, '("Spherical part of the polarization potential will be added to the &
                                   &channel (target state,l,m), coefficient value: (",3i5,") ",e25.15)') &
                                   it1, l1, m1, -sph_cpl * alpha0
                        this % a(lqt, lambda) = this % a(lqt, lambda) - sph_cpl * alpha0  ! add the spherical part of the polarizability
                    else
                        sph_cpl = 0.0_wp
                    end if
 
                    ! Add the non-spherical part of the polarizability; the polarization potential here corresponds to:
                    ! alpha_{2}*C_{i}*r_{-4}, where the angular behaviour of C_{i} = xy, xz, yz, xx, yy, zz
                    ! Note that we assume below that the polarizability tensor has the same values (alpha2) for all components.
                    ! Also note that as long as l1,m1 and l2,m2 belong to the same IR, only the totally symmetric components
                    ! of the pol. tensor (xx,yy,zz) may contribute.
                    if (use_alpha2) then
 
                        write(*, '("Non-spherical part of the polarization potential &
                                   &(target state,l1,m2,l2,m2): ",5i5)') it1, l1, m1, l2, m2
 
                        ! x^2 component
                        cpl = x2_x00 * couplings % rgaunt(l1, 0, l2, m1, 0, m2) + &
                              x2_x20 * couplings % rgaunt(l1, 2, l2, m1, 0, m2) + &
                              x2_x22 * couplings % rgaunt(l1, 2, l2, m1, 2, m2)
                        write(*, '("x^2 coefficient: ",e25.15)') -cpl * alpha2
                        this % a(lqt, lambda) = this % a(lqt, lambda) - cpl * alpha2
 
                        ! y^2 component
                        cpl = y2_x00 * couplings % rgaunt(l1, 0, l2, m1, 0, m2) + &
                              y2_x20 * couplings % rgaunt(l1, 2, l2, m1, 0, m2) + &
                              y2_x22 * couplings % rgaunt(l1, 2, l2, m1, 2, m2)
                        write(*, '("y^2 coefficient: ",e25.15)') -cpl * alpha2
                        this % a(lqt, lambda) = this % a(lqt, lambda) - cpl * alpha2
 
                        ! z^2 component
                        cpl = z2_x00 * couplings % rgaunt(l1, 0, l2, m1, 0, m2) + &
                              z2_x20 * couplings % rgaunt(l1, 2, l2, m1, 0, m2)
                        write(*, '("z^2 coefficient: ",e25.15)') -cpl * alpha2
                        this % a(lqt, lambda) = this % a(lqt, lambda) - cpl * alpha2
 
                        ! xz component
                        cpl = couplings % rgaunt(l1, 2, l2, m1, 1, m2) * rothree
                        write(*, '("xz coefficient: ",e25.15)') -cpl * alpha2
                        this % a(lqt, lambda) = this % a(lqt, lambda) - cpl * alpha2
 
                        ! yz component
                        cpl = couplings % rgaunt(l1, 2, l2, m1, -1, m2) * rothree
                        write(*, '("yz coefficient: ",e25.15)') -cpl * alpha2
                        this % a(lqt, lambda) = this % a(lqt, lambda) - cpl * alpha2
 
                        ! xy component
                        cpl = couplings % rgaunt(l1, 2, l2, m1, -2, m2) * rothree
                        write(*, '("xy coefficient: ",e25.15)') -cpl * alpha2
                        this % a(lqt, lambda) = this % a(lqt, lambda) - cpl * alpha2
 
                  end if  ! use_alpha2
 
               end if  ! use_pol
 
            end do  ! ch_2
 
         end do  ! ch_1
 
    end subroutine get_channel_couplings
 
 
    subroutine write_rmt_data (input, inner_properties, target_properties, solutions, intf)
 
        use algorithms,                      only: findloc, insertion_sort
        use class_molecular_properties_data, only: molecular_properties_data
        use const_gbl,                       only: stdout
        use consts_mpi_ci,                   only: pass_to_cdenprop, no_diagonalization
        use mpi_gbl,                         only: myrank, master, mpiint, mpi_xermsg, mpi_mod_file_open_write, &
                                                   mpi_mod_file_close, mpi_mod_file_set_size, mpi_mod_file_write, mpi_mod_barrier
        use options_module,                  only: optionsset
        use solutionhandler_module,          only: solutionhandler
 
        type(optionsset),                   intent(in) :: input
        type(molecular_properties_data),    intent(in) :: inner_properties
        type(molecular_properties_data),    intent(in) :: target_properties
        type(solutionhandler), allocatable, intent(in) :: solutions(:)
        type(outerinterface),  allocatable, intent(in) :: intf(:)
 
        real(rmt_real),   allocatable :: cf(:,:,:), wmat(:,:), crlv(:,:), etarg(:), rg(:), eig(:), r_points(:)
        integer(rmt_int), allocatable :: iidip(:), ifdip(:), lm_rg(:,:), ltarg(:), starg(:), nconat(:), l2p(:), m2p(:), ichl(:), &
                                         propblocks(:)
        integer,          allocatable :: ions(:)
 
        real(rmt_real)   :: bbloch = 0, prop, rmatr, delta_r
        integer          :: i, j, k, l, nsym, m, irr_map(8,8), iblock, u, v, nsymm, ind, ip, irr1, irr2, it1, it2, lu, neig, nnuc, &
                            ierr, nblock, sym1, sym2, ref
        integer(mpiint)  :: fh
        integer(rmt_int) :: s, s1, s2, ntarg, n_rg, nelc, nz, lran2, lamax, inast, nchmx, nstmx, lmaxp1, nchan, lrgl, nspn, npty, &
                            mnp1, mxstat, lrang2, lmax, nfdm
 
        ! number of N + 1 irrs for which we have inner solution and dipoles
        inast   = count(iand(input % opts(:) % vector_storage_method, pass_to_cdenprop) /= 0 .and. &
                        input % opts(:) % diagonalization_flag /= no_diagonalization)
 
        ref     = maxloc(intf(:) % ntarg, 1)                ! irr with the first non-zero number of outer region states
        nnuc    = target_properties % no_nuclei             ! number of nuclei in the molecule
        nfdm    = input % nfdm                              ! number of boundary amplitudes evaluation radii *inside* R-matrix sphere
        delta_r = input % delta_r                           ! spacing between the evaluation distances
        rmatr   = input % rmatr                             ! radius of the inner region
        nsymm   = size(input % opts(:))                     ! number of N + 1 irrs for which we have inputs (= size(intf(:)))
        ntarg   = intf(ref) % ntarg                         ! number of target ("ionic") states
        nelc    = input % opts(1) % num_electrons - 1       ! number of electrons of the target (ion)
        nz      = sum(target_properties % charge(1:nnuc))   ! number of protons in the molecule
 
        mxstat  = maxval(inner_properties % dense_blocks(1:inner_properties % no_blocks) % mat_dimen)  ! largest dipole block size
 
        ! avoid the routine altogether when there are no data for any reason
        if (nsymm == 0 .or. inast == 0 .or. ntarg == 0 .or. mxstat <= 0) then
            write (stdout, *) 'Nothing to do in RMT interface'
            return
        end if
 
        lrang2  = 0         ! angular momentum limit
        lamax   = 0         ! angular momentum limit
        nchmx   = 0         ! maximal number of channels across irreducible representations
        nstmx   = 0         ! maximal number of inner region eigenstates across representations
 
        ! get energy-sorted list of ionic states for use in outer region (only a subset of all if IDTARG is incomplete)
        ions = intf(ref) % idtarg
        call insertion_sort(ions)
 
        ! open the output file and reset its contents
        call mpi_mod_file_open_write('molecular_data', fh, ierr)
        call mpi_mod_file_set_size(fh, 0)
 
        allocate (propblocks(inner_properties % no_blocks))
 
        write (stdout, '(/,"Writing RMT data file")')
        write (stdout, '(  "=====================")')
        write (stdout, '(/,"  Inner region dipole block max size: ",I0)') mxstat
        write (stdout, '(/,"  Ionic states used in outer region: ")')
        write (stdout, '(/,"    ",20I5)') ions
 
        ! store all (N + 1)-electron dipole blocks
        do m = -1, 1
 
            ! find out which IRR pairs yield non-zero matrix elements of Y[1,m]
            irr_map = 0
            nblock = 0
            do ip = 1, inner_properties % no_blocks
                if (inner_properties % dense_index(3, ip) == 1 .and. &
                    inner_properties % dense_index(4, ip) == m) then
                    sym1 = inner_properties % dense_index(1, ip)
                    sym2 = inner_properties % dense_index(2, ip)
                    irr1 = inner_properties % symmetries_index(1, sym1) + 1
                    irr2 = inner_properties % symmetries_index(1, sym2) + 1
                    nblock = nblock + 1
                    propblocks(nblock) = ip
                    irr_map(irr1, irr2) = nblock
                end if
            end do
 
            ! get some free memory
            if (allocated(iidip)) deallocate (iidip)
            if (allocated(ifdip)) deallocate (ifdip)
            allocate (iidip(nblock))
            allocate (ifdip(nblock))
            iidip(:) = -1
            ifdip(:) = -1
 
            ! compose sequence of bra and ket symmetry indices
            do ip = 1, nblock
                sym1 = inner_properties % dense_index(1, propblocks(ip))
                sym2 = inner_properties % dense_index(2, propblocks(ip))
                irr1 = inner_properties % symmetries_index(1, sym1) + 1
                irr2 = inner_properties % symmetries_index(1, sym2) + 1
                iidip(ip) = irr1
                ifdip(ip) = irr2
            end do
 
            ! write dipole metadata
            call mpi_mod_file_write(fh, int(nblock, rmt_int))
            call mpi_mod_file_write(fh, mxstat)
            call mpi_mod_file_write(fh, mxstat)
            call mpi_mod_file_write(fh, iidip, nblock)
            call mpi_mod_file_write(fh, ifdip, nblock)
 
            ! collective write of the ScaLAPACK matrices to the stream file
            do ip = 1, nblock
                associate(d => inner_properties % dense_blocks(propblocks(ip)))
                    if (d % CV_is_scalapack) then
                        call mpi_mod_file_write(fh, int(d % mat_dimen_r), int(d % mat_dimen_c), int(d % nprow), int(d % npcol), &
                                                    int(d % scal_block_size), int(d % scal_block_size), d % CV, &
                                                    int(d % local_row_dimen), int(d % local_col_dimen))
                        call mpi_mod_barrier(ierr)
                    else
                        call mpi_mod_file_write(fh, d % CV, int(d % mat_dimen_r), int(d % mat_dimen_c))
                    end if
                end associate
            end do
        end do
 
        ! set up angular momentum limits
        do i = 1, nsymm
            if (iand(input % opts(i) % vector_storage_method, pass_to_cdenprop) /= 0 .and. &
                input % opts(i) % diagonalization_flag /= no_diagonalization) then
                if (intf(i) % nchan > 0) then
                    lrang2 = max(int(lrang2), maxval(intf(i) % lchl(1:intf(i) % nchan)))
                    lamax  = max(int(lamax), intf(i) % ismax)
                    nchmx  = max(int(nchmx), intf(i) % nchan)
                    lmaxp1 = lrang2
                end if
                nstmx = max(int(nstmx), solutions(i) % vec_dimen)
            end if
        end do
 
        ! number of target (ion) states
        call mpi_mod_file_write(fh, ntarg)
 
        ! allocate memory for auxiliary arrays
        allocate (etarg(ntarg), ltarg(ntarg), starg(ntarg), nconat(ntarg), l2p(nchmx), m2p(nchmx), eig(nstmx), &
                  wmat(nchmx, nstmx), cf(nchmx, nchmx, lamax), ichl(nchmx), crlv(ntarg, ntarg), r_points(nfdm + 1), stat = ierr)
        if (ierr /= 0) then
            call mpi_xermsg('Postprocessing_module', 'write_rmt_data', 'Memory allocation failure.', 1, 1)
        end if
 
        ! target data
        do i = 1, ntarg
            etarg(i) = target_properties % energies(ions(i))
            ltarg(i) = target_properties % symmetries_index(1, target_properties % energies_index(2, ions(i))) + 1
            starg(i) = target_properties % symmetries_index(2, target_properties % energies_index(2, ions(i)))
        end do
 
        ! finite difference points inside (and at) the R-matrix sphere
        do i = 1, nfdm + 1
            r_points(i) = rmatr - (nfdm + 1 - i) * delta_r
        end do
 
        ! store all N-electron dipole blocks
        do m = -1, 1
            crlv(:,:) = 0
            do ip = 1, target_properties % non_zero_properties
                if (target_properties % properties_index(3, ip) == 1 .and. &
                    target_properties % properties_index(4, ip) == m) then
                    ! get state information for this property
                    it1  = target_properties % properties_index(1, ip)
                    it2  = target_properties % properties_index(2, ip)
                    prop = target_properties % properties(ip)
                    ! check that these two states are requested in the outer region
                    it1 = findloc(ions, it1, 1)
                    it2 = findloc(ions, it2, 1)
                    if (it1 > 0 .and. it2 > 0) then
                        crlv(it1, it2) = prop
                        crlv(it2, it1) = prop
                    end if
                end if
            end do
 
            call mpi_mod_file_write(fh, crlv, int(ntarg), int(ntarg))
        end do
 
        ! generate the coupling coefficients needed in RMT to construct the laser-target and laser-electron asymptotic potentials
        lmax = 0
        do i = 1, nsymm
            if (iand(input % opts(i) % vector_storage_method, pass_to_cdenprop) /= 0 .and. &
                input % opts(i) % diagonalization_flag /= no_diagonalization) then
                if (intf(i) % nchan > 0) then
                    lmax = max(int(lmax), maxval(intf(i) % lchl(1:intf(i) % nchan)))
                end if
            end if
        end do
        call generate_couplings(lmax, n_rg, rg, lm_rg)
 
        call mpi_mod_file_write(fh, n_rg)
        call mpi_mod_file_write(fh, rg, int(n_rg))
        call mpi_mod_file_write(fh, lm_rg, 6, int(n_rg))
        call mpi_mod_file_write(fh, nelc)
        call mpi_mod_file_write(fh, nz)
        call mpi_mod_file_write(fh, lrang2)
        call mpi_mod_file_write(fh, lamax)
        call mpi_mod_file_write(fh, ntarg)
        call mpi_mod_file_write(fh, inast)
        call mpi_mod_file_write(fh, nchmx)
        call mpi_mod_file_write(fh, nstmx)
        call mpi_mod_file_write(fh, lmaxp1)
        call mpi_mod_file_write(fh, rmatr)
        call mpi_mod_file_write(fh, bbloch)
        call mpi_mod_file_write(fh, etarg, int(ntarg))
        call mpi_mod_file_write(fh, ltarg, int(ntarg))
        call mpi_mod_file_write(fh, starg, int(ntarg))
        call mpi_mod_file_write(fh, nfdm)
        call mpi_mod_file_write(fh, delta_r)
        call mpi_mod_file_write(fh, r_points, int(nfdm + 1))
 
        ! write channel information (per symmetry)
        do i = 1, nsymm
            if (iand(input % opts(i) % vector_storage_method, pass_to_cdenprop) /= 0 .and. &
                input % opts(i) % diagonalization_flag /= no_diagonalization) then
 
                lrgl  = 1 + input % opts(i) % lambda
                nspn  = int(2 * input % opts(i) % spin + 1, rmt_int)
                npty  = 0
                nchan = intf(i) % nchan
                mnp1  = solutions(i) % vec_dimen
 
                do j = 1, ntarg
                    nconat(j) = count(intf(i) % ichl(1:nchan) == j)
                end do
 
                l2p(:)  = 0
                m2p(:)  = 0
                ichl(:) = 0
 
                do j = 1, nchan
                    l2p(j) = intf(i) % lchl(j)
                    m2p(j) = intf(i) % mchl(j)
                    ichl(j) = intf(i) % ichl(j)
                end do
 
                call mpi_mod_file_write(fh, lrgl)
                call mpi_mod_file_write(fh, nspn)
                call mpi_mod_file_write(fh, npty)
                call mpi_mod_file_write(fh, nchan)
                call mpi_mod_file_write(fh, mnp1)
                call mpi_mod_file_write(fh, nconat, int(ntarg))
                call mpi_mod_file_write(fh, l2p, int(nchmx))
                call mpi_mod_file_write(fh, m2p, int(nchmx))
 
                eig(:) = 0
                eig(1:mnp1) = solutions(i) % ci_vec % ei(1:mnp1) + solutions(i) % core_energy
 
                wmat(:,:) = 0
                if (myrank == master) wmat(1:nchan, 1:mnp1) = intf(i) % wamp(nfdm + 1) % CV(1:nchan, 1:mnp1)
 
                call mpi_mod_file_write(fh, eig, int(nstmx))
                call mpi_mod_file_write(fh, wmat, int(nchmx), int(nstmx))
 
                cf(:,:,:) = 0
                do l = 1, intf(i) % ismax
                    do k = 1, nchan
                        do j = 1, nchan
                            u = max(j, k)
                            v = min(j, k)
                            ind = u * (u - 1) / 2 + v
                            cf(j, k, l) = intf(i) % a(ind, l)
                        end do
                    end do
                end do
 
                call mpi_mod_file_write(fh, cf, int(nchmx), int(nchmx), int(lamax))
                call mpi_mod_file_write(fh, ichl, int(nchmx))
 
                s1 = intf(i) % wamp(1) % mat_dimen_r
                s2 = intf(i) % wamp(1) % mat_dimen_c
 
                call mpi_mod_file_write(fh, s1)
                call mpi_mod_file_write(fh, s2)
 
                ! collective write of the ScaLAPACK matrix to the stream file
                do j = 1, nfdm
                    associate(w => intf(i) % wamp(j))
                        if (w % CV_is_scalapack) then
                            call mpi_mod_file_write(fh, int(w % mat_dimen_r), int(w % mat_dimen_c), &
                                                        int(w % nprow), int(w % npcol), &
                                                        int(w % scal_block_size), int(w % scal_block_size), w % CV, &
                                                        int(w % local_row_dimen), int(w % local_col_dimen))
                        else
                            call mpi_mod_file_write(fh, w % CV, int(w % mat_dimen_r), int(w % mat_dimen_c))
                        end if
                    end associate
                end do
 
            end if
        end do
 
        call mpi_mod_file_close(fh, ierr)
 
    end subroutine write_rmt_data
 
 
    subroutine generate_couplings (maxl, n_rg, rg, lm_rg)
 
        use coupling_obj_gbl, only: couplings_type
 
        integer(rmt_int),              intent(in)    :: maxl
        integer(rmt_int),              intent(out)   :: n_rg
        real(rmt_real),   allocatable, intent(inout) :: rg(:)
        integer(rmt_int), allocatable, intent(inout) :: lm_rg(:,:)
 
        type(couplings_type) :: cpl
 
        integer  :: l1, l2, l3, m1, m2, m3, err
        real(wp) :: tmp
 
        call cpl % prec_cgaunt(int(maxl))
 
        l3 = 1
 
        ! How many non-zero dipole couplings there are
        n_rg = 0
        do l1 = 0, maxl
            do m1 = -l1, l1
                do l2 = 0, maxl
                    do m2 = -l2, l2
                        do m3 = -l3, l3
                            tmp = cpl % rgaunt(l1, l2, l3, m1, m2, m3)
                            if (tmp /= 0) n_rg = n_rg + 1
                        end do !m3
                    end do !m2
                end do !l2
            end do !m1
        end do !l1
 
        if (allocated(rg)) deallocate (rg)
        if (allocated(lm_rg)) deallocate (lm_rg)
        n_rg = max(n_rg, 1_rmt_int)  ! make sure that even if all couplings are zero this routine allocates rg, lm_rg
        allocate (rg(n_rg), lm_rg(6, n_rg), stat = err)
        rg = 0
        lm_rg = -2
 
        ! Save them
        n_rg = 0
        do l1 = 0, maxl
            do m1 = -l1, l1
                do l2 = 0, maxl
                    do m2 = -l2, l2
                        do m3 = -l3, l3
                            tmp = cpl % rgaunt(l1, l2, l3, m1, m2, m3)
                            if (tmp /= 0) then
                                n_rg = n_rg + 1
                                rg(n_rg) = tmp
                                lm_rg(1:6, n_rg) = (/ l1, m1, l2, m2, l3, m3 /)
                            end if
                        end do !m3
                    end do !m2
                end do !l2
            end do !m1
        end do !l1
 
   end subroutine generate_couplings
 
end module postprocessing_module