multidip_io.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_io
 
    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, iostat_end
 
    ! GBTOlib
    use blas_lapack_gbl,  only: blasint
#ifdef WITH_MMAP
    use file_mapping_gbl, only: file_mapping
#endif
    use phys_const_gbl,   only: pi, imu, to_ev
    use precisn_gbl,      only: wp, wp_bytes
 
    implicit none
 
#ifdef WITH_SCALAPACK
    integer(blasint), external :: numroc
#endif
 
    integer :: myproc = 1
    integer :: nprocs = 1
 
    type kmatrix
        integer  :: mgvn, stot, nescat, nchan, nopen, ndopen, nchsq, lukmt, nkset, maxne
        real(wp), allocatable :: escat(:)
    contains
        procedure :: get_kmatrix
        procedure :: read_kmatrix_file
        procedure :: reset_kmatrix_position
    end type kmatrix
 
 
    type scatakcoeffs
        real(wp), allocatable :: rea(:, :, :), ima(:, :, :), evchl(:)
        integer,  allocatable :: ichl(:), lvchl(:), mvchl(:)
        integer               :: mgvn, nesc, nchan, nstat, lusct
    contains
        procedure :: get_wfncoeffs
        procedure :: read_wfncoeffs_file
        procedure :: reset_wfncoeffs_position
    end type scatakcoeffs
 
 
    type mappedmatrix
        real(real64), pointer :: mat(:, :) => null()  
#ifdef WITH_MMAP
        type(file_mapping)    :: mapping
#endif
        logical          :: distributed  = .false.    
        integer(blasint) :: desc(9)      = 0
        integer(blasint) :: row_context  = 0
        integer(blasint) :: blk_context  = 0
        integer(blasint) :: block_size   = 0
    contains
        procedure :: load => load_mapped_matrix
        final     :: destruct_mapped_matrix
    end type mappedmatrix
 
 
    type moleculardata
 
        ! inner transition dipoles storage
        type(mappedmatrix), allocatable :: dipx(:), dipy(:), dipz(:)
        integer(int32),     allocatable :: iidipx(:), iidipy(:), iidipz(:)
        integer(int32),     allocatable :: ifdipx(:), ifdipy(:), ifdipz(:)
 
        ! boundary amplitudes
        type(mappedmatrix), allocatable :: wamp(:)
        type(mappedmatrix), allocatable :: wmat2(:, :)
 
        ! number of electrons, number of protons
        integer(int32) :: nelc, nz
 
        ! other data
        real(real64)                :: rmatr
        real(real64),   allocatable :: crlv(:,:,:), eig(:,:), rg(:), etarg(:), gaunt(:, :, :), r_points(:)
        integer(int32), allocatable :: mgvns(:), stot(:), mnp1(:), nchan(:), l2p(:,:), m2p(:,:), ichl(:,:), lm_rg(:,:)
 
    end type moleculardata
 
contains
 
    subroutine load_mapped_matrix (this, filename, u, offset, rows, cols)
 
#ifdef WITH_SCALAPACK
        use mpi
        use mpi_gbl, only: mpi_xermsg
 
        integer, parameter :: mpiint = kind(mpi_comm_world)
        integer, parameter :: mpiofs = mpi_offset_kind
 
        character(len=200) :: msg
 
        integer(mpiint) :: ndims, gsizes(2), distrb(2), dargs(2), psizes(2), locsiz, fh, ierr, datype
        integer(mpiint) :: comm_size, comm_rank, stat(MPI_STATUS_SIZE), mtwo = 2, dims(2)
        integer(mpiofs) :: offs
#endif
 
        class(mappedmatrix), intent(inout) :: this
        character(len=*),    intent(in)    :: filename
        integer(int32),      intent(in)    :: u, rows, cols
        integer(int64),      intent(in)    :: offset
 
        integer(int64)   :: length
        integer(blasint) :: zero = 0, one = 1, nprow, npcol, lrows, lcols, min_block_size, max_block_size, ld, info, myprow, mypcol
        integer(blasint) :: brows, bcols
 
#ifdef WITH_SCALAPACK
        if (this % distributed) then
            call mpi_comm_size(mpi_comm_world, comm_size, ierr)
            call assert_status('Failed to obtain world communicator size: ', int(ierr))
            call mpi_comm_rank(mpi_comm_world, comm_rank, ierr)
            call assert_status('Failed to obtain world communicator rank: ', int(ierr))
 
            dims = 0
            call mpi_dims_create(comm_size, mtwo, dims, ierr)
            call assert_status('Failed to reshape world communicator to rectangular grid: ', int(ierr))
            nprow = maxval(dims)
            npcol = minval(dims)
 
            if (npcol * nprow /= nprocs) then
                write (msg, '(a,3(i0,a))') 'BLACS grid ', nprow, 'x', npcol, ' does not use all ', nprocs, ' processes. Fix it.'
                call mpi_xermsg('multidip_io', 'load_mapped_matrix', trim(msg), 1, 1);
            end if
 
            call blacs_get(zero, zero, this % row_context)
            call blacs_get(zero, zero, this % blk_context)
            call blacs_gridinit(this % row_context, 'R', one, nprow * npcol)
            call blacs_gridinit(this % blk_context, 'R', nprow, npcol)
            call blacs_gridinfo(this % blk_context, nprow, npcol, myprow, mypcol)
 
            min_block_size = 1
            max_block_size = 512  ! rather arbitrary, but this particular choice coincides with a common page size (4 kiB)
 
            this % block_size = min(rows / nprow, cols / npcol)
            this % block_size = min(this % block_size, max_block_size)
            this % block_size = max(this % block_size, min_block_size)
 
            brows = rows
            bcols = cols
            lrows = numroc(brows, this % block_size, myprow, zero, nprow)
            lcols = numroc(bcols, this % block_size, mypcol, zero, npcol)
            ld = max(one, lrows)
 
            if (int(lrows, int64)*int(lcols, int64) > huge(1_int32)) then
                print '(a,i0,a,i0,a)', 'WARNING: Local portion of the dipole matrix has size ', lrows, ' x ', lcols, '.'
                print '(9x,a)', 'This means that it cannot be indexed with 4-byte integers.'
                print '(9x,a,/)', 'If reading fails, try distributing the calculation more.'
            end if
 
            call descinit(this % desc, brows, bcols, this % block_size, this % block_size, zero, zero, this % blk_context, ld, info)
 
            allocate (this % mat(ld, lcols))
 
            offs   = offset - 1
            ndims  = 2
            gsizes = [ rows, cols ]
            distrb = [ mpi_distribute_cyclic, mpi_distribute_cyclic ]
            dargs  = [ this % block_size, this % block_size ]
            psizes = [ nprow, npcol ]
            locsiz = lrows * lcols
 
            call mpi_file_open(mpi_comm_world, filename, mpi_mode_rdonly, mpi_info_null, fh, ierr)
            call assert_status('Failed to open file for MPI-IO: ', int(ierr))
            call mpi_type_create_darray(comm_size, comm_rank, ndims, gsizes, distrb, dargs, &
                                  psizes, mpi_order_fortran, mpi_real8, datype, ierr)
            call assert_status('Failed to define distributed array data type: ', int(ierr))
            call mpi_type_commit(datype, ierr)
            call assert_status('Failed to commit distributed array data type: ', int(ierr))
            call mpi_file_set_view(fh, offs, mpi_real8, datype, 'native', mpi_info_null, ierr)
            call assert_status('Failed to define file view: ', int(ierr))
            call mpi_file_read_all(fh, this % mat, locsiz, mpi_real8, stat, ierr)
            call assert_status('Failed to read from file: ', int(ierr))
            call mpi_type_free(datype, ierr)
            call assert_status('Failed to release distributed array data type: ', int(ierr))
            call mpi_file_close(fh, ierr)
            call assert_status('Failed to close file: ', int(ierr))
        else
#endif
#ifdef WITH_MMAP
            length = int(rows, int64) * int(cols, int64) * wp_bytes
            call this % mapping % init(filename, offset - 1, length)
            call c_f_pointer(this % mapping % ptr, this % mat, [rows, cols])
#else
            allocate (this % mat(rows, cols))
            read (u) this % mat
#endif
#ifdef WITH_SCALAPACK
        end if
#endif
 
    end subroutine load_mapped_matrix
 
 
    subroutine destruct_mapped_matrix (this)
 
        type(mappedmatrix), intent(inout) :: this
 
#ifdef WITH_SCALAPACK
        if (this % distributed) then
            deallocate (this % mat)
        else
#endif
#ifndef WITH_MMAP
            if (associated(this % mat)) then
                deallocate (this % mat)
            end if
#endif
#ifdef WITH_SCALAPACK
        end if
#endif
        nullify (this % mat)
 
    end subroutine destruct_mapped_matrix
 
 
    subroutine read_kmatrices (km, lukmt, nkset)
 
        use mpi_gbl, only: mpi_xermsg
 
        type(kmatrix), allocatable, intent(inout) :: km(:)
        integer,                    intent(in)    :: lukmt(:), nkset(:)
 
        character(len=200) :: msg
        integer :: i
 
        do i = 1, size(km)
 
            km(i) % lukmt = lukmt(i)
            km(i) % nkset = nkset(i)
 
            call km(i) % read_kmatrix_file
 
            if (km(i) % nescat /= km(1) % nescat) then
                write (msg, '(a,2(i0,a))') 'ERROR: K-matrix units ', lukmt(1), ' and ', lukmt(i), &
                                           ' hold different number of scattering energies!'
                call mpi_xermsg('multidip_io', 'read_kmatrices', trim(msg), 1, 1)
            end if
 
            if (any(km(i) % escat(1 : km(i) % nescat) /= km(1) % escat(1 : km(1) % nescat))) then
                write (msg, '(a,2(i0,a))') 'ERROR: K-matrix units ', lukmt(1), ' and ', lukmt(i), &
                                           ' hold different scattering energies!'
                call mpi_xermsg('multidip_io', 'read_kmatrices', trim(msg), 1, 1)
            end if
 
        end do
 
    end subroutine read_kmatrices
 
 
    subroutine read_kmatrix_file (km)
 
        use mpi_gbl, only: mpi_xermsg
 
        class(kmatrix), intent(inout) :: km
 
        integer :: k, n, ifail, key, nset, nrec, ninfo, ndata, ntarg, nvib, ndis, ie, gutot, ion
        integer :: nopen, ndopen, nchsq, a, b, c, mye, nene, ierr
        logical :: iwarn
        real(wp) :: dE1, dE2, r, rmass, en
        real(wp), allocatable :: kmat(:)
        character(len=80) :: title, msg
 
        print '(/,A,I0,A,I0,/)', 'Reading K-matrices from unit ', km % lukmt, ', set ', km % nkset
 
        call getset(km % lukmt, km % nkset, 11, 'UNFORMATTED', ifail)
 
        read (km % lukmt) key, nset, nrec, ninfo, ndata
        read (km % lukmt) title
        read (km % lukmt) km % mgvn, km % stot, gutot, ion, r, rmass
        read (km % lukmt) ntarg, nvib, ndis, km % nchan, km % maxne
 
        print '(2x,A,I0)', 'mgvn  = ', km % mgvn
        print '(2x,A,I0)', 'stot  = ', km % stot
        print '(2x,A,I0)', 'ntarg = ', ntarg
        print '(2x,A,I0)', 'nvib  = ', nvib
        print '(2x,A,I0)', 'ndis  = ', ndis
        print '(2x,A,I0)', 'nchan = ', km % nchan
        print '(2x,A,I0)', 'maxne = ', km % maxne
 
        print '(/,2x,A,/)', 'Energy ranges:'
 
        km % nescat = 0
        do ie = 1, km % maxne
            read (km % lukmt) k, n, de1, de2
            print '(4x,I4,I8,F8.3,F8.3)', k, n, de1, de2
            km % nescat = km % nescat + n
        end do
 
        nene = (km % nescat + nprocs - 1) / nprocs
        allocate (km % escat(km % nescat), kmat(km % nchan * (km % nchan + 1)))
 
        iwarn = .true.
        mye = 0
        do ie = 1, km % nescat
            read (km % lukmt, iostat = ierr) nopen, ndopen, nchsq, en!, kmat(1:nchsq)
            if (ierr == iostat_end) then
                print '(/,2x,A,I0,A,I0,A)', 'WARNING: Unit ', km % lukmt, ' contains only ', ie - 1, ' K-matrices.'
                print '(11x,A)', 'This means that CFASYM failed during RSOLVE for some energies.'
                print '(11x,A)', 'The scattering energies MUST be consistent across symmetries!'
                km % nescat = ie - 1
                exit
            else if (ierr /= 0) then
                write (msg, '(2x,a,i0,a)') 'ERROR ', ierr, ' while reading K-matrix file.'
                call mpi_xermsg('multidip_io', 'read_kmatrix_file', trim(msg), 1, 1)
            end if
            km % escat(ie) = en/2  ! Ry -> a.u.
            if (iwarn .and. nopen < km % nchan) then
                print '(/,2x,A,I0,A)', 'WARNING: Unit ', km % lukmt, ' contains only open-open subset of the full K-matrix.'
                print '(11x,A)', 'Use IKTYPE = 1 in RSOLVE to save full K-matrices.'
                print '(11x,A)', 'You may get spurious below-threshold noise.'
                iwarn = .false.
            end if
        end do
 
        deallocate (kmat)
 
    end subroutine read_kmatrix_file
 
 
    subroutine reset_kmatrix_position (km, skip)
 
        class(kmatrix),    intent(in) :: km
        integer, optional, intent(in) :: skip
 
        integer :: i
 
        ! rewind to the beginning of the K-matrix unit
        call getset(km % lukmt, km % nkset, 11, 'UNFORMATTED', i)
 
        read (km % lukmt)
        read (km % lukmt)
        read (km % lukmt)
        read (km % lukmt)
 
        do i = 1, km % maxne
            read (km % lukmt)
        end do
 
        if (present(skip)) then
            do i = 1, skip
                read (km % lukmt)
            end do
        end if
 
    end subroutine reset_kmatrix_position
 
 
    subroutine get_kmatrix (km, kmat, skip)
 
        use mpi_gbl, only: mpi_xermsg
 
        class(kmatrix),        intent(in)    :: km
        real(wp), allocatable, intent(inout) :: kmat(:, :)
        integer,  optional,    intent(in)    :: skip
 
        real(wp), allocatable :: buffer(:)
        character(len=200)    :: msg
 
        real(wp) :: en
        integer  :: a, b, c, i, ierr, nopen, ndopen, nchsq
 
        allocate (buffer(km % nchan * (km % nchan + 1)))
 
        read (km % lukmt, iostat = ierr) nopen, ndopen, nchsq, en, buffer(1:nchsq)
 
        if (ierr == iostat_end) then
            write (msg, '(a,i0,a)') 'ERROR: Unit ', km % lukmt, ' ended while reading K-matrix!'
            call mpi_xermsg('multidip_io', 'get_kmatrix', trim(msg), 1, 1)
        else if (ierr /= 0) then
            write (msg, '(2x,a,2(i0,a))') 'ERROR ', ierr, ' while reading K-matrix file unit ', km % lukmt, '!'
            call mpi_xermsg('multidip_io', 'get_kmatrix', trim(msg), 1, 1)
        end if
 
        if (allocated(kmat)) then
            if (size(kmat, 1) /= km % nchan .or. size(kmat, 2) /= km % nchan) then
                deallocate (kmat)
            end if
        end if
 
        if (.not. allocated(kmat)) then
            allocate (kmat(km % nchan, km % nchan))
        end if
 
        c = 0
        do a = 1, nopen
            do b = 1, a
                c = c + 1
                kmat(a, b) = buffer(c)
                kmat(b, a) = buffer(c)
            end do
        end do
 
        if (present(skip)) then
            do i = 1, skip
                read (km % lukmt, iostat = ierr)  ! ignore EOF
            end do
        end if
 
    end subroutine get_kmatrix
 
 
    subroutine read_wfncoeffs (ak, lusct)
 
        type(scatakcoeffs), allocatable, intent(inout) :: ak(:)
        integer, intent(in) :: lusct(:)
        integer :: i
 
        do i = 1, size(ak)
            ak(i) % lusct = lusct(i)
            call ak(i) % read_wfncoeffs_file
        end do
 
    end subroutine read_wfncoeffs
 
 
    subroutine read_bndcoeffs (bnd, Ei, lubnd, mgvn0, stot0)
 
        real(wp), intent(out) :: bnd(:), Ei
        integer,  intent(in)  :: lubnd, mgvn0, stot0
 
        real(wp), allocatable :: xvec(:), etot(:), vtemp(:)
        real(wp)              :: rr
        integer               :: nbset, nchan, mgvn, stot, gutot, nstat, nbound, iprint, iwrite, ifail
        character(len=11)     :: bform
 
        bform = 'unformatted'
        nbset = 1
        iprint = 0
        iwrite = output_unit
 
        open (lubnd, action = 'read', form = bform)
 
        call readbh(lubnd, nbset, nchan, mgvn, stot, gutot, nstat, nbound, rr, bform, iprint, iwrite, ifail)
 
        if (mgvn /= mgvn0 .or. stot /= stot0) then
            print '(*(a,i0))', 'ERROR: Bound coefficients at unit ', lubnd, ' are for MGVN = ', mgvn, ', STOT = ', stot, &
                               ', but I am looking for MGVN = ', mgvn0, ', STOT = ', stot0
            stop 1
        end if
 
        nbound = 1  ! read just the first bound state, ignore the rest
        allocate (xvec(nbound*nchan), vtemp(nbound), etot(nbound))
 
        call readbc(nstat, etot, vtemp, bnd, nbound, nchan, xvec)
 
        ei = etot(1)
 
        close (lubnd)
 
    end subroutine read_bndcoeffs
 
 
    subroutine read_wfncoeffs_file (ak)
 
        class(scatakcoeffs), intent(inout) :: ak
 
        real(wp), allocatable :: ReA(:), ImA(:)
        real(wp)              :: rr
        integer               :: i, j, keysc, nset, nrec, ninfo, ndata, stot, gutot, nscat, ich, ie, mye, nene
        character(len=80)     :: header
 
        print '(/,A,I0,/)', 'Reading wave function coefficients from unit ', ak % lusct
 
        open (ak % lusct, status = 'old', action = 'read', form = 'unformatted')
 
        read (ak % lusct) keysc, nset, nrec, ninfo, ndata
        read (ak % lusct) header
        read (ak % lusct) nscat, ak % mgvn, stot, gutot, ak % nstat, ak % nchan, ak % nesc
 
        print '(2x,A,I0)', 'mgvn  = ', ak % mgvn
        print '(2x,A,I0)', 'stot  = ', stot
        print '(2x,A,I0)', 'nstat = ', ak % nstat
        print '(2x,A,I0)', 'nchan = ', ak % nchan
        print '(2x,A,I0)', 'nesc  = ', ak % nesc
 
        read (ak % lusct) rr
 
        if (allocated(ak % ichl))  deallocate (ak % ichl)
        if (allocated(ak % lvchl)) deallocate (ak % lvchl)
        if (allocated(ak % mvchl)) deallocate (ak % mvchl)
        if (allocated(ak % evchl)) deallocate (ak % evchl)
 
        allocate (ak % ichl(ak % nchan))
        allocate (ak % lvchl(ak % nchan))
        allocate (ak % mvchl(ak % nchan))
        allocate (ak % evchl(ak % nchan))
 
        read (ak % lusct) ak % ichl, ak % lvchl, ak % mvchl, ak % evchl
 
        print '(/,2x,A)', 'channel table'
        do j = 1, ak % nchan
            print '(2x,4I6,F15.7)', j, ak % ichl(j), ak % lvchl(j), ak % mvchl(j), ak % evchl(j)
        end do
 
        close (ak % lusct)
 
    end subroutine read_wfncoeffs_file
 
 
    subroutine reset_wfncoeffs_position (ak, skip)
 
        class(scatakcoeffs), intent(in) :: ak
        integer, optional,   intent(in) :: skip
 
        integer :: i, ich
 
        ! rewind to the beginning of the K-matrix unit
        call getset(ak % lusct, 1, 88, 'UNFORMATTED', i)
 
        read (ak % lusct)
        read (ak % lusct)
        read (ak % lusct)
        read (ak % lusct)
        read (ak % lusct)
 
        if (present(skip)) then
            do i = 1, skip
                do ich = 1, ak % nchan
                    read (ak % lusct)
                    read (ak % lusct)
                end do
            end do
        end if
 
    end subroutine reset_wfncoeffs_position
 
 
    subroutine get_wfncoeffs (ak, E, Re_Ak, Im_Ak, skip)
 
        use mpi_gbl, only: mpi_xermsg
 
        class(scatakcoeffs),   intent(in)    :: ak
        real(wp), allocatable, intent(inout) :: Re_Ak(:, :), Im_Ak(:, :)
        integer,  optional,    intent(in)    :: skip
        real(wp),              intent(inout) :: E
 
        character(len=250) :: msg, iomsg
 
        integer  :: i, ich, ierr
 
        do ich = 1, ak % nchan
            read (ak % lusct, iostat = ierr, iomsg = iomsg) e, i, re_ak(:, ich)
            read (ak % lusct, iostat = ierr, iomsg = iomsg) e, i, im_ak(:, ich)
            if (ierr == iostat_end) then
                write (msg, '(2x,a,i0,a)') 'ERROR: Unit ', ak % lusct, ' ended while reading Ak-coeffs!'
                call mpi_xermsg('multidip_io', 'get_wfncoeffs', trim(msg), 1, 1)
            else if (ierr /= 0) then
                write (msg, '(2x,a,2(i0,a),a)') 'ERROR ', ierr, ' while reading Ak-coeffs file unit ', ak % lusct, ': ', iomsg
                call mpi_xermsg('multidip_io', 'get_wfncoeffs', trim(msg), 1, 1)
            end if
        end do
 
        if (present(skip)) then
            do i = 1, skip
                do ich = 1, ak % nchan
                    read (ak % lusct, iostat = ierr)  ! ignore EOF
                    read (ak % lusct, iostat = ierr)  ! ignore EOF
                end do
            end do
        end if
 
    end subroutine get_wfncoeffs
 
 
    subroutine read_molecular_data (moldat, filename, mpiio, read_wmat2)
 
        use multidip_params, only: check_dipoles
 
        type(moleculardata), intent(inout) :: moldat
        character(*),        intent(in)    :: filename
        logical,             intent(in)    :: mpiio, read_wmat2
 
        integer(int32), allocatable :: ltarg(:), starg(:), nconat(:)
        real(real64),   allocatable :: r_points(:), cf(:, :, :)
 
        integer(int32) :: i, j, u, s, s1, s2, ntarg, nrg, lrang2, lamax, inast, nchmx, nstmx, lmaxp1, nfdm
        integer(int32) :: lrgl, nspn, npty
        integer(int64) :: offset, length
        real(real64)   :: bbloch, delta_r
 
        open (newunit = u, file = filename, access = 'stream', status = 'old', action = 'read')
 
        print '(/,3A,/)', 'Reading file "', filename, '"'
 
        read (u) s, s1, s2
        length = int(s1, int64) * int(s2, int64) * wp_bytes
        allocate (moldat % iidipy(s), moldat % ifdipy(s), moldat % dipy(s))
        read (u) moldat % iidipy, moldat % ifdipy
        inquire (u, pos = offset)
        do i = 1, s
            moldat % dipy(i) % distributed = mpiio
            call moldat % dipy(i) % load(filename, u, offset, s1, s2)
            offset = offset + length
        end do
        print '(2x,A,*(1x,I0))', '[m = -1] iidip =', moldat % iidipy
        print '(2x,A,*(1x,I0))', '         ifdip =', moldat % ifdipy
 
        if (check_dipoles) then
            do i = 1, s
                if (sum(abs(moldat % dipy(i) % mat)) == 0) then
                    print '(/,2x,A,/)', 'WARNING: Norm of some Y dipoles is zero. Something went wrong in CDENPROP?'
                end if
            end do
        end if
 
        read (u, pos = offset) s, s1, s2
        length = int(s1, int64) * int(s2, int64) * wp_bytes
        allocate (moldat % iidipz(s), moldat % ifdipz(s), moldat % dipz(s))
        read (u) moldat % iidipz, moldat % ifdipz
        inquire (u, pos = offset)
        do i = 1, s
            moldat % dipz(i) % distributed = mpiio
            call moldat % dipz(i) % load(filename, u, offset, s1, s2)
            offset = offset + length
        end do
        print '(2x,A,*(1x,I0))', '[m =  0] iidip =', moldat % iidipz
        print '(2x,A,*(1x,I0))', '         ifdip =', moldat % ifdipz
 
        if (check_dipoles) then
            do i = 1, s
                if (sum(abs(moldat % dipz(i) % mat)) == 0) then
                    print '(/,2x,A,/)', 'WARNING: Norm of some Z dipoles is zero. Something went wrong in CDENPROP?'
                end if
            end do
        end if
 
        read (u, pos = offset) s, s1, s2
        length = int(s1, int64) * int(s2, int64) * wp_bytes
        allocate (moldat % iidipx(s), moldat % ifdipx(s), moldat % dipx(s))
        read (u) moldat % iidipx, moldat % ifdipx
        inquire (u, pos = offset)
        do i = 1, s
            moldat % dipx(i) % distributed = mpiio
            call moldat % dipx(i) % load(filename, u, offset, s1, s2)
            offset = offset + length
        end do
        print '(2x,A,*(1x,I0))', '[m = +1] iidip =', moldat % iidipx
        print '(2x,A,*(1x,I0))', '         ifdip =', moldat % ifdipx
 
        if (check_dipoles) then
            do i = 1, s
                if (sum(abs(moldat % dipx(i) % mat)) == 0) then
                    print '(/,2x,A,/)', 'WARNING: Norm of some X dipoles is zero. Something went wrong in CDENPROP?'
                end if
            end do
        end if
 
        read (u, pos = offset) ntarg
        print '(/,2x,A,I0)', 'ntarg = ', ntarg
        allocate (moldat % crlv(ntarg, ntarg, 3))
        read (u) moldat % crlv
 
        read (u) nrg
        allocate (moldat % rg(nrg), moldat % lm_rg(6, nrg))
        read (u) moldat % rg, moldat % lm_rg
 
        read (u) moldat % nelc, moldat % nz, lrang2, lamax, ntarg, inast, nchmx, nstmx, lmaxp1
        read (u) moldat % rmatr, bbloch
        allocate (moldat % etarg(ntarg), ltarg(ntarg), starg(ntarg))
        read (u) moldat % etarg, ltarg, starg
        read (u) nfdm, delta_r
        allocate (moldat % r_points(nfdm + 1))
        read (u) moldat % r_points
        print '(2x,A,*(1x,F8.3))', 'sampling radii:', moldat % r_points
 
        print '(/,2x,A)', 'number of inner eigenstates per symmetry'
 
        allocate (nconat(ntarg), moldat % l2p(nchmx, inast), moldat % m2p(nchmx, inast), moldat % eig(nstmx, inast), &
                  moldat % wamp(inast), cf(nchmx, nchmx, lamax), moldat % ichl(nchmx, inast), &
                  moldat % mgvns(inast), moldat % stot(inast), moldat % nchan(inast), moldat % mnp1(inast), &
                  moldat % wmat2(nfdm, inast))
 
        inquire (u, pos = offset)
 
        length = int(nchmx, int64) * int(nstmx, int64) * wp_bytes
 
        do i = 1, inast
 
            read (u, pos = offset) lrgl, nspn, npty, moldat % nchan(i), moldat % mnp1(i)
            moldat % mgvns(i) = lrgl - 1
            moldat % stot(i) = nspn
            print '(4x,I3,I8)', moldat % mgvns(i), moldat % mnp1(i)
            read (u) nconat, moldat % l2p(:, i), moldat % m2p(:, i)
            read (u) moldat % eig(:, i)
            inquire (u, pos = offset)
            call moldat % wamp(i) % load(filename, u, offset, nchmx, nstmx)
            offset = offset + length
            read (u, pos = offset) cf, moldat % ichl(:, i)
            read (u) s1, s2
            inquire (u, pos = offset)
 
            ! read or map wmat2 only when explicitly required; it is large and used only for debugging
            if (read_wmat2) then
                do j = 1, nfdm
                    call moldat % wmat2(j, i) % load(filename, u, offset, s1, s2)
                    offset = offset + int(s1, int64) * int(s2, int64) * wp_bytes
                end do
            else
                offset = offset + int(s1, int64) * int(s2, int64) * nfdm * wp_bytes
            end if
 
        end do
 
        do i = 1, inast
            print '(/,2x,A,I0,A,I0,A)', 'channels in symmetry #', i, ' (IRR ', moldat % mgvns(i), '):'
            do j = 1, moldat % nchan(i)
                print '(2x,I4,A,I6,I3,SP,I6)', j, ':', moldat % ichl(j, i), moldat % l2p(j, i), moldat % m2p(j, i)
                if (moldat % ichl(j, i) > ntarg) then
                    print '(4(A,I0),A)', 'Warning: Channel ', j, ' in symmetry ', moldat % mgvns(i), ' is coupled to state ', &
                                 moldat % ichl(j, i), ', which is beyond the number of targets ', ntarg, ' (bug in mpi-scatci!)'
                    moldat % ichl(j, i) = ntarg
                end if
            end do
        end do
 
        close (u)
 
        call setup_angular_integrals(moldat)
 
    end subroutine read_molecular_data
 
 
    subroutine setup_angular_integrals (moldat)
 
        use mpi_gbl, only: mpi_xermsg
 
        type(moleculardata), intent(inout) :: moldat
 
        integer :: maxl, i, l1, m1, l2, m2, l3, m3
 
        maxl = max(0_int32, maxval(moldat % lm_rg))
 
        allocate (moldat % gaunt((maxl + 1)**2, (maxl + 1)**2, -1:1))
 
        moldat % gaunt = 0
 
        do i = 1, size(moldat % rg)
            l1 = moldat % lm_rg(1, i); m1 = moldat % lm_rg(2, i)
            l2 = moldat % lm_rg(3, i); m2 = moldat % lm_rg(4, i)
            l3 = moldat % lm_rg(5, i); m3 = moldat % lm_rg(6, i)
            if (l3 /= 1 .or. abs(m3) > 1) then
                call mpi_xermsg('multidip_io', 'setup_angular_integrals', &
                                'Unexpected layout of the angular integral storage in molecular_data', 1, 1)
            end if
            moldat % gaunt(l1*(l1+1) + m1 + 1, l2*(l2+1) + m2 + 1, m3) = moldat % rg(i)
        end do
 
    end subroutine setup_angular_integrals
 
 
    subroutine get_diptrans (moldat, I, iidip, ifdip)
 
        type(moleculardata),  intent(in)    :: moldat
        integer,              intent(in)    :: I
        integer, allocatable, intent(inout) :: iidip(:), ifdip(:)
 
        if (allocated(iidip)) deallocate (iidip)
        if (allocated(ifdip)) deallocate (ifdip)
 
        select case (i)
            case (1)
                iidip = moldat % iidipx
                ifdip = moldat % ifdipx
            case (2)
                iidip = moldat % iidipy
                ifdip = moldat % ifdipy
            case (3)
                iidip = moldat % iidipz
                ifdip = moldat % ifdipz
        end select
 
    end subroutine get_diptrans
 
 
    subroutine apply_dipole_matrix (moldat, component, irrpair, transp, nf, nn, X, Y)
 
        use multidip_params,  only: rone, rzero
        use multidip_special, only: blas_dgemm => dgemm
 
        integer(blasint), parameter :: one = 1, zero = 0
 
        type(moleculardata), intent(in)    :: moldat
        integer,             intent(in)    :: component, irrpair, nf, nn
        character(len=1),    intent(in)    :: transp
        real(wp),            intent(inout) :: X(:, :, :), Y(:, :, :)
 
        integer(blasint) :: lda, ldb, ldc, m, n, k, info, Xdesc(9), Xdesc2D(9), Ydesc(9), Ydesc2D(9), desc(9), gn
        integer(blasint) :: Xlocr, Xlocc, Ylocr, Ylocc, ldX, ldY, rctx, bctx, rbksz, cbksz, myprow, mypcol, nprow, npcol
        logical          :: dist
 
        real(wp), allocatable :: Xloc(:, :), Yloc(:, :)
        real(wp), pointer     :: dips(:, :)
 
        ldb = size(x, 1)
        ldc = size(y, 1)
 
        m = nf
        n = size(x, 2) * size(x, 3)
        k = nn
 
        ! pick the correct inner dipole block (or its local portion if MPI-IO was used)
        select case (component)
            case (1)
                dist = moldat % dipx(irrpair) % distributed
                desc = moldat % dipx(irrpair) % desc
                rctx = moldat % dipx(irrpair) % row_context
                bctx = moldat % dipx(irrpair) % blk_context
                lda  = size(moldat % dipx(irrpair) % mat, 1)
                dips => moldat % dipx(irrpair) % mat(:, :)
            case (2)
                dist = moldat % dipy(irrpair) % distributed
                desc = moldat % dipy(irrpair) % desc
                rctx = moldat % dipy(irrpair) % row_context
                bctx = moldat % dipy(irrpair) % blk_context
                lda  = size(moldat % dipy(irrpair) % mat, 1)
                dips => moldat % dipy(irrpair) % mat(:, :)
            case (3)
                dist = moldat % dipz(irrpair) % distributed
                desc = moldat % dipz(irrpair) % desc
                rctx = moldat % dipz(irrpair) % row_context
                bctx = moldat % dipz(irrpair) % blk_context
                lda  = size(moldat % dipz(irrpair) % mat, 1)
                dips => moldat % dipz(irrpair) % mat(:, :)
        end select
 
#ifdef WITH_SCALAPACK
        if (dist) then
            ! total number of columns in X, Y (summed over all parallel images)
            gn = n * nprocs
 
            ! get information about the BLACS process grid
            call blacs_gridinfo(bctx, nprow, npcol, myprow, mypcol)
 
            ! set up descriptors for the MULTIDIP distribution of vectors X and Y (round-robin in energies)
            call descinit(xdesc, k, gn, ldb, n, zero, zero, rctx, ldb, info)
            call assert_status('Failed to set up X BLACS descriptor: ', int(info))
            call descinit(ydesc, m, gn, ldc, n, zero, zero, rctx, ldc, info)
            call assert_status('Failed to set up Y BLACS descriptor: ', int(info))
 
            ! pick a reasonable block size for good load balancing (aim at one ScaLAPACK block per process at least)
            rbksz = clip(min(nf, nn) / int(nprow), 1, 64)
            cbksz = clip(int(gn / npcol), 1, 64)
 
            ! find out the number of rows and columns in the local portion of the redistributed X matrix
            xlocr = numroc(k, rbksz, myprow, zero, nprow);  ldx = max(one, xlocr)
            xlocc = numroc(gn, cbksz, mypcol, zero, npcol)
 
            ! find out the number of rows and columns in the local portion of the redistributed Y matrix
            ylocr = numroc(m, rbksz, myprow, zero, nprow);  ldy = max(one, ylocr)
            ylocc = numroc(gn, cbksz, mypcol, zero, npcol)
 
            ! allocate the local portions of the redistributed X and Y matrices
            allocate (xloc(ldx, xlocc), yloc(ldy, ylocc))
 
            ! set up descriptors for the ScaLAPACK distribution of vectors X and Y (proper two-dimensional block-cyclic)
            call descinit(xdesc2d, k, gn, rbksz, cbksz, zero, zero, bctx, ldx, info)
            call assert_status('Failed to set up X2D BLACS descriptor: ', int(info))
            call descinit(ydesc2d, m, gn, rbksz, cbksz, zero, zero, bctx, ldy, info)
            call assert_status('Failed to set up Y2D BLACS descriptor: ', int(info))
 
            ! redistribute X -> Xloc; multiply; redistribute Yloc -> Y
            call pdgemr2d(k, gn, x, one, one, xdesc, xloc, one, one, xdesc2d, bctx)
            call pdgemm(transp, 'N', m, gn, k, rone, dips, one, one, desc, xloc, one, one, xdesc2d, rzero, yloc, one, one, ydesc2d)
            call pdgemr2d(m, gn, yloc, one, one, ydesc2d, y, one, one, ydesc, bctx)
        else
#endif
            call blas_dgemm(transp, 'N', m, n, k, rone, dips, lda, x, ldb, rzero, y, ldc)
#ifdef WITH_SCALAPACK
        end if
#endif
 
    end subroutine apply_dipole_matrix
 
 
    subroutine assert_status (message, errcode)
 
        use mpi_gbl, only: mpi_xermsg
 
        character(len=*), intent(in) :: message
        integer,          intent(in) :: errcode
 
        if (errcode /= 0) then
            call mpi_xermsg('multidip_io', 'assert_status', trim(message), errcode, 1)
        end if
 
    end subroutine assert_status
 
 
    integer function clip (x, a, b) result (y)
 
        integer :: x, a, b
 
        y = max(a, min(x, b))
 
    end function clip
 
 
    subroutine apply_boundary_amplitudes (moldat, irr, transp, X, Y)
 
        use multidip_special, only: blas_dgemm => dgemm
 
        type(moleculardata),   intent(in)    :: moldat
        integer,               intent(in)    :: irr
        character(len=1),      intent(in)    :: transp
        real(wp),              intent(inout) :: X(:,:), Y(:,:)
 
        integer(blasint)  :: m, n, k, lda, ldb, ldc
        real(wp)          :: alpha = 1, beta = 0
        real(wp), pointer :: wamp(:, :)
 
        if (transp == 'T') then
            m = moldat % mnp1(irr)
            n = size(x, 2)
            k = moldat % nchan(irr)
        else
            m = moldat % nchan(irr)
            n = size(x, 2)
            k = moldat % mnp1(irr)
        end if
 
        lda = size(moldat % wamp(irr) % mat, 1)
        ldb = size(x, 1)
        ldc = size(y, 1)
 
        if (lda < merge(m, k, transp == 'N') .or. ldb < k .or. ldc < m) then
            print '(a)', 'ERROR: Incompatible matrix dimensions in apply_boundary_amplitudes'
            print '(7x,a,6(", ",a,1x,i0))', transp, 'm =', m, 'n =', n, 'k =', k, 'lda =', lda, 'ldb =', ldb, 'ldc =', ldc
            stop 1
        end if
 
        wamp => moldat % wamp(irr) % mat(:, :)
 
        call blas_dgemm(transp, 'N', m, n, k, alpha, wamp, lda, x, ldb, beta, y, ldc)
 
    end subroutine apply_boundary_amplitudes
 
 
    subroutine scale_boundary_amplitudes (moldat, irr, v, vw)
 
        type(moleculardata),   intent(in)    :: moldat
        integer,               intent(in)    :: irr
        real(wp), allocatable, intent(in)    :: v(:)
        real(wp), allocatable, intent(inout) :: vw(:, :)
 
        integer :: k, nstat, p, nchan
 
        nstat = moldat % mnp1(irr)
        nchan = moldat % nchan(irr)
 
        !$omp parallel default(none) private(p, k) firstprivate(nchan, nstat, irr) shared(vw, v, moldat)
        !$omp do collapse(2)
        do p = 1, nchan
            do k = 1, nstat
                vw(k, p) = v(k) * moldat % wamp(irr) % mat(p, k)
            end do
        end do
        !$omp end do
        !$omp end parallel
 
    end subroutine scale_boundary_amplitudes
 
 
    subroutine write_energy_grid (escat)
 
        real(wp), intent(in) :: escat(:)
 
        integer :: u
 
        open (newunit = u, file = 'pe_energies.txt', action = 'write', form = 'formatted')
        write (u, '(E25.15)') escat
        close (u)
 
    end subroutine write_energy_grid
 
 
    subroutine write_partial_wave_moments (moldat, M, nesc, suffix)
 
        use mpi_gbl, only: mpi_reduce_inplace_sum_wp
 
        type(moleculardata),      intent(in) :: moldat
        character(len=*),         intent(in) :: suffix
        integer,                  intent(in) :: nesc
        complex(wp), allocatable, intent(in) :: M(:, :, :)
 
        complex(wp), allocatable :: val(:)
        real(wp),    allocatable :: reval(:), imval(:)
 
        integer            :: ie, irr, mgvn, u, ichan, ichlf, lf, mf, mye
        character(len=100) :: filename
 
        allocate (val(nesc))
 
        ! there is one file per irreducible representation and channel
        do irr = 1, size(moldat % mgvns)
            mgvn = moldat % mgvns(irr)
            do ichan = 1, size(m, 1)
                val = 0
                mye = 0
                do ie = myproc, nesc, nprocs
                    mye = mye + 1
                    val(ie) = m(ichan, mye, mgvn)
                end do
 
                ! combine data from all images to master
                reval = real(val)
                imval = aimag(val)
                call mpi_reduce_inplace_sum_wp(reval, nesc)
                call mpi_reduce_inplace_sum_wp(imval, nesc)
                val = cmplx(reval, imval, wp)
 
                ! the first image performs the writing
                if (myproc == 1 .and. any(val /= 0)) then
                    ichlf = moldat % ichl(ichan, irr)
                    lf = moldat % l2p(ichan, irr)
                    mf = moldat % m2p(ichan, irr)
                    write (filename, '(3(A,I0),A,SP,I0,3A)') 'pwdip-', mgvn, '-(', ichlf, ',', lf, ',', mf, ')', suffix, '.txt'
                    open (newunit = u, file = filename, action = 'write', form = 'formatted')
                    write (u, '(2E25.15)') val
                    close (u)
                end if
            end do
        end do
 
    end subroutine write_partial_wave_moments
 
 
    subroutine write_raw_dipoles (M, chains, nesc, stem)
 
        use mpi_gbl, only: mpi_reduce_inplace_sum_wp
 
        character(len=*),         intent(in) :: stem
        integer,     allocatable, intent(in) :: chains(:, :)
        complex(wp), allocatable, intent(in) :: M(:, :, :, :)
        integer,                  intent(in) :: nesc
 
        character(len=100)          :: filename, history
        character(len=1), parameter :: xyz(-1:1) = [ 'y', 'z', 'x' ]
        character(len=1), parameter :: sph(-1:1) = [ 'm', '0', 'p' ]
 
        complex(wp), allocatable :: val(:)
        real(wp),    allocatable :: reval(:), imval(:)
 
        integer :: u, itarg, ntarg, ichain, nchain, ipw, npw, icomp, lf, mf, ie, mye
 
        npw = size(m, 2)
        nchain = size(m, 3)
        ntarg = size(m, 4)
 
        allocate (val(nesc))
 
        ! there is one output file per target, absorption chain and partial wave
        do itarg = 1, ntarg
            do ichain = 1, nchain
                do ipw = 1, npw
 
                    val = 0
                    mye = 0
                    do ie = myproc, nesc, nprocs
                        mye = mye + 1
                        val(ie) = m(mye, ipw, ichain, itarg)
                    end do
 
                    ! combine data from all images to master
                    reval = real(val)
                    imval = aimag(val)
                    call mpi_reduce_inplace_sum_wp(reval, nesc)
                    call mpi_reduce_inplace_sum_wp(imval, nesc)
                    val = cmplx(reval, imval, wp)
 
                    ! the first image performs the writing
                    if (myproc == 1 .and. any(val /= 0)) then
                        lf = int(sqrt(real(ipw - 1, wp)))
                        mf = ipw - lf*lf - lf - 1
                        history = ''
                        do icomp = 1, size(chains, 1)
                            select case (trim(stem))
                                case ('xyz'); history = xyz(chains(icomp, ichain)) // trim(history)
                                case ('sph'); history = sph(chains(icomp, ichain)) // trim(history)
                            end select
                        end do
                        write (filename, '(3A,I0,A,I0,A,SP,I0,A)') 'rawdip-', trim(history), '-(', itarg, ',', lf, ',', mf, ').txt'
                        open (newunit = u, file = filename, action = 'write', form = 'formatted')
                        write (u, '(2E25.15)') val
                        close (u)
                    end if
 
                end do
            end do
        end do
 
    end subroutine write_raw_dipoles
 
 
    subroutine write_cross_section (cs, nesc, erange, filename)
 
        use mpi_gbl, only: mpi_reduce_inplace_sum_wp
 
        character(len=*), intent(in) :: filename
        real(wp),         intent(in) :: cs(:, :)
        integer,          intent(in) :: nesc, erange(2)
 
        real(wp), allocatable :: csg(:, :), buffer(:)
 
        integer :: u, ie, mye, nene, ntgt
 
        ntgt = size(cs, 1)
        nene = size(cs, 2)
 
        ! a local copy of 'cs' for gathering data from parallel images
        allocate (csg(ntgt, nesc))
        csg = 0
        mye = 0
        do ie = myproc, nesc, nprocs
            mye = mye + 1
            csg(:, ie) = cs(:, mye)
        end do
 
        ! combine data from all images to master
        buffer = reshape(csg, [size(csg)])
        call mpi_reduce_inplace_sum_wp(buffer, size(buffer))
        csg = reshape(buffer, shape(csg))
 
        ! the first image now writes (non-zero) cross sections to file
        if (myproc == 1 .and. any(csg(2:, :) /= 0)) then
            open (newunit = u, file = filename, action = 'write', form = 'formatted')
            do ie = max(1, erange(1)), min(nesc, erange(2))
                write (u, '(*(E25.15))') csg(:, ie)
            end do
            close (u)
        end if
 
    end subroutine write_cross_section
 
end module multidip_io