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

    ! 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

    integer :: myproc = 1
    integer :: nprocs = 1

    real(wp), parameter :: alpha = 1/137.03599907_wp    

    character(len=1), parameter :: compt(3) = ['x', 'y', 'z']

    type kmatrix
        integer  :: mgvn, stot, nescat, nchan, nopen, ndopen, nchsq
        real(wp), allocatable :: escat(:), kmat(:,:,:)
    end type kmatrix


    type scatakcoeffs
        real(wp), allocatable :: rea(:, :, :), ima(:, :, :), evchl(:), escat(:)
        integer,  allocatable :: ichl(:), lvchl(:), mvchl(:)
        integer               :: mgvn, nesc, nchan, nstat
    end type scatakcoeffs


    type moleculardata

        ! inner transition dipoles storage
#ifdef WITH_MMAP
        real(real64),   pointer     :: dipx(:,:,:), dipy(:,:,:), dipz(:,:,:)
#else
        real(real64),   allocatable :: dipx(:,:,:), dipy(:,:,:), dipz(:,:,:)
#endif
        integer(int32), allocatable :: iidipx(:), iidipy(:), iidipz(:)
        integer(int32), allocatable :: ifdipx(:), ifdipy(:), ifdipz(:)

#ifdef WITH_MMAP
        ! helper structures for file memory mapping
        type(file_mapping) :: dipx_mmap, dipy_mmap, dipz_mmap
#endif

        ! other data
        real(real64)                :: rmatr
        real(real64),   allocatable :: crlv(:,:,:), eig(:,:), wamp(:,:,:), rg(:), etarg(:), gaunt(:, :, :)
        integer(int32), allocatable :: mgvns(:), mnp1(:), nchan(:), l2p(:,:), m2p(:,:), ichl(:,:), lm_rg(:,:)

    contains

        ! class destructor
        final :: destruct_molecular_data

    end type moleculardata

contains

    subroutine read_kmatrices (km, lukmt, nkset)

        type(KMatrix), allocatable, intent(inout) :: km(:)
        integer,                    intent(in)    :: lukmt(:), nkset(:)

        integer :: i, k, n, ifail, u, key, nset, nrec, ninfo, ndata, ntarg, nvib, ndis, maxne, ie, gutot, ion
        integer :: nopen, ndopen, nchsq, a, b, c, mye, nene
        logical :: iwarn
        real(wp) :: dE1, dE2, r, rmass, en
        real(wp), allocatable :: kmat(:)
        character(len=80) :: title

        do i = 1, size(km)

            u = lukmt(i)
            nset = nkset(i)

            print '(/,A,I0,A,I0,/)', 'Reading K-matrices from unit ', u, ', set ', nset

            call getset(u, nset, 11, 'UNFORMATTED', ifail)

            read (u) key, nset, nrec, ninfo, ndata
            read (u) title
            read (u) km(i) % mgvn, km(i) % stot, gutot, ion, r, rmass
            read (u) ntarg, nvib, ndis, km(i) % nchan, maxne

            print '(2x,A,I0)', 'mgvn  = ', km(i) % mgvn
            print '(2x,A,I0)', 'stot  = ', km(i) % 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(i) % nchan
            print '(2x,A,I0)', 'maxne = ', maxne

            print '(/,2x,A,/)', 'Energy ranges:'

            km(i) % nescat = 0
            do ie = 1, maxne
                read (u) k, n, de1, de2
                print '(4x,I4,I8,F8.3,F8.3)', k, n, de1, de2
                km(i) % nescat = km(i) % nescat + n
            end do

            nene = (km(i) % nescat + nprocs - 1) / nprocs
            allocate (km(i) % kmat(km(i) % nchan, km(i) % nchan, nene), &
                      km(i) % escat(km(i) % nescat), &
                      kmat(km(i) % nchan * (km(i) % nchan + 1)))

            iwarn = .true.
            km(i) % kmat = 0
            mye = 0
            do ie = 1, km(i) % nescat
                read (u) nopen, ndopen, nchsq, en, kmat(1:nchsq)
                km(i) % escat(ie) = en/2  ! Ry -> a.u.
                c = 0
                if (mod(ie - 1, nprocs) + 1 /= myproc) cycle
                mye = mye + 1
                do a = 1, nopen
                    do b = 1, a
                        c = c + 1
                        km(i) % kmat(a, b, mye) = kmat(c)
                        km(i) % kmat(b, a, mye) = kmat(c)
                    end do
                end do
                if (iwarn .and. nopen < km(i) % nchan) then
                    print '(/,2x,A,I0,A)', 'WARNING: Unit ', u, ' 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 do

    end subroutine read_kmatrices


    subroutine read_wfncoeffs (ak, lusct)

        type(ScatAkCoeffs), allocatable, intent(inout) :: ak(:)
        integer, intent(in) :: lusct(:)

        real(wp), allocatable :: ReA(:), ImA(:)
        integer :: i, j, u, keysc, nset, nrec, ninfo, ndata, stot, gutot, nscat, ich, ie, mye, nene
        real(wp) :: rr

        character(len=80) :: header

        do i = 1, size(ak)

            u = lusct(i)

            print '(/,A,I0,/)', 'Reading wave function coefficients from unit ', u

            open (u, status = 'old', action = 'read', form = 'unformatted')

            read (u) keysc, nset, nrec, ninfo, ndata
            read (u) header
            read (u) nscat, ak(i) % mgvn, stot, gutot, ak(i) % nstat, ak(i) % nchan, ak(i) % nesc

            print '(2x,A,I0)', 'mgvn  = ', ak(i) % mgvn
            print '(2x,A,I0)', 'stot  = ', stot
            print '(2x,A,I0)', 'nstat = ', ak(i) % nstat
            print '(2x,A,I0)', 'nchan = ', ak(i) % nchan
            print '(2x,A,I0)', 'nesc  = ', ak(i) % nesc

            read (u) rr

            allocate (ak(i) % ichl(ak(i) % nchan), &
                      ak(i) % lvchl(ak(i) % nchan), &
                      ak(i) % mvchl(ak(i) % nchan), &
                      ak(i) % evchl(ak(i) % nchan))

            read (u) ak(i) % ichl, ak(i) % lvchl, ak(i) % mvchl, ak(i) % evchl

            print '(/,2x,A)', 'channel table'
            do j = 1, ak(i) % nchan
                print '(2x,4I6,F15.7)', j, ak(i) % ichl(j), ak(i) % lvchl(j), ak(i) % mvchl(j), ak(i) % evchl(j)
            end do

            nene = (ak(i) % nesc + nprocs - 1) / nprocs
            allocate (ak(i) % ReA(ak(i) % nstat, ak(i) % nchan, nene), &
                      ak(i) % ImA(ak(i) % nstat, ak(i) % nchan, nene), &
                      ak(i) % escat(ak(i) % nesc), &
                      rea(ak(i) % nstat), &
                      ima(ak(i) % nstat))

            mye = 0
            do ie = 1, ak(i) % nesc
                if (mod(ie - 1, nprocs) + 1 == myproc) mye = mye + 1
                do ich = 1, ak(i) % nchan
                    read (u) ak(i) % escat(ie), j, rea
                    read (u) ak(i) % escat(ie), j, ima
                    if (mod(ie - 1, nprocs) + 1 == myproc) then
                        ak(i) % ReA(:, ich, mye) = rea
                        ak(i) % ImA(:, ich, mye) = ima
                    end if
                end do
            end do

            close (u)

        end do

    end subroutine read_wfncoeffs


    subroutine read_molecular_data (moldat)

        type(MolecularData), intent(inout) :: moldat

        integer(int32), allocatable :: ltarg(:), starg(:), nconat(:)
        real(real64),   allocatable :: r_points(:), cf(:, :, :), wmat2(:, :, :)

        integer(int32) :: i, j, u, s, s1, s2, ntarg, nrg, nelc, nz, lrang2, lamax, inast, nchmx, nstmx, lmaxp1, nfdm
        integer(int32) :: lrgl, nspn, npty
        integer(int64) :: offset, length
        real(real64)   :: bbloch, delta_r

        character(*), parameter :: filename = 'molecular_data'

        open (newunit = u, file = filename, access = 'stream', status = 'old', action = 'read')

        print '(/,A,/)', 'Reading file "molecular_data"'

        read (u) s, s1, s2
        length = s * int(s1, int64) * int(s2, int64) * wp_bytes
        allocate (moldat % iidipy(s), moldat % ifdipy(s))
        read (u) moldat % iidipy, moldat % ifdipy
        inquire (u, pos = offset)
#ifdef WITH_MMAP
        call moldat % dipy_mmap % init(filename, offset - 1, length)
        call c_f_pointer(moldat % dipy_mmap % ptr, moldat % dipy, [s1, s2, s])
#else
        allocate (moldat % dipy(s1, s2, s))
        read (u) moldat % dipy
#endif
        offset = offset + length
        print '(2x,A,*(1x,I0))', '[m = -1] iidip =', moldat % iidipy
        print '(2x,A,*(1x,I0))', '         ifdip =', moldat % ifdipy

        read (u, pos = offset) s, s1, s2
        length = s * int(s1, int64) * int(s2, int64) * wp_bytes
        allocate (moldat % iidipz(s), moldat % ifdipz(s))
        read (u) moldat % iidipz, moldat % ifdipz
        inquire (u, pos = offset)
#ifdef WITH_MMAP
        call moldat % dipz_mmap % init(filename, offset - 1, length)
        call c_f_pointer(moldat % dipz_mmap % ptr, moldat % dipz, [s1, s2, s])
#else
        allocate (moldat % dipz(s1, s2, s))
        read (u) moldat % dipz
#endif
        offset = offset + length
        print '(2x,A,*(1x,I0))', '[m =  0] iidip =', moldat % iidipz
        print '(2x,A,*(1x,I0))', '         ifdip =', moldat % ifdipz

        read (u, pos = offset) s, s1, s2
        length = s * int(s1, int64) * int(s2, int64) * wp_bytes
        allocate (moldat % iidipx(s), moldat % ifdipx(s))
        read (u) moldat % iidipx, moldat % ifdipx
        inquire (u, pos = offset)
#ifdef WITH_MMAP
        call moldat % dipx_mmap % init(filename, offset - 1, length)
        call c_f_pointer(moldat % dipx_mmap % ptr, moldat % dipx, [s1, s2, s])
#else
        allocate (moldat % dipx(s1, s2, s))
        read (u) moldat % dipx
#endif
        offset = offset + length
        print '(2x,A,*(1x,I0))', '[m = +1] iidip =', moldat % iidipx
        print '(2x,A,*(1x,I0))', '         ifdip =', moldat % ifdipx

        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) nelc, 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 (r_points(nfdm + 1))
        read (u) r_points
        print '(2x,A,*(1x,F8.3))', 'sampling radii:', 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(nchmx, nstmx, inast), cf(nchmx, nchmx, lamax), moldat % ichl(nchmx, inast), &
                  moldat % mgvns(inast), moldat % nchan(inast), moldat % mnp1(inast))

        do i = 1, inast

            read (u) lrgl, nspn, npty, moldat % nchan(i), moldat % mnp1(i)
            moldat % mgvns(i) = lrgl - 1
            print '(4x,I3,I8)', moldat % mgvns(i), moldat % mnp1(i)
            read (u) nconat, moldat % l2p(:, i), moldat % m2p(:, i)
            read (u) moldat % eig(:, i), moldat % wamp(:, :, i), cf, moldat % ichl(:, i)
            read (u) s1, s2
            allocate (wmat2(s1, s2, nfdm))
            read (u) wmat2
            deallocate (wmat2)

        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,I3,A,I6,I3,SP,I6)', j, ':', moldat % ichl(j, i), moldat % l2p(j, i), moldat % m2p(j, i)
            end do
        end do

        close (u)

        call setup_angular_integrals(moldat)

    end subroutine read_molecular_data


    subroutine destruct_molecular_data (moldat)

        type(MolecularData), intent(inout) :: moldat

#ifdef WITH_MMAP
        call moldat % dipx_mmap % final
        call moldat % dipy_mmap % final
        call moldat % dipz_mmap % final
#endif
    end subroutine destruct_molecular_data


    subroutine setup_angular_integrals (moldat)

        type(MolecularData), intent(inout) :: moldat

        integer :: maxl, i, l1, m1, l2, m2, l3, m3

        maxl = 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
                print '(A)', 'Unexpected layout of the angular integral storage in molecular_data'
                stop
            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, I, s, transp, nf, nn, X, Y)

        use multidip_params,  only: rone, rzero
        use multidip_special, only: blas_dgemm => dgemm

        type(MolecularData), intent(in)    :: moldat
        integer,             intent(in)    :: I, s, nf, nn
        character(len=1),    intent(in)    :: transp
        real(wp),            intent(inout) :: X(:, :, :), Y(:, :, :)

        integer(blasint) :: lda, ldb, ldc, m, n, k

        ldb = size(x, 1)
        ldc = size(y, 1)

        m = nf
        n = size(x, 2) * size(x, 3)
        k = nn

        select case (i)
            case (1)
                lda = size(moldat % dipx, 1)
                call blas_dgemm(transp, 'N', m, n, k, rone, moldat % dipx(:,:,s), lda, x, ldb, rzero, y, ldc)
            case (2)
                lda = size(moldat % dipy, 1)
                call blas_dgemm(transp, 'N', m, n, k, rone, moldat % dipy(:,:,s), lda, x, ldb, rzero, y, ldc)
            case (3)
                lda = size(moldat % dipz, 1)
                call blas_dgemm(transp, 'N', m, n, k, rone, moldat % dipz(:,:,s), lda, x, ldb, rzero, y, ldc)
        end select

    end subroutine apply_dipole_matrix


    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

        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, 1)
        ldb = size(x, 1)
        ldc = size(y, 1)

        call blas_dgemm(transp, 'N', m, n, k, alpha, moldat % wamp(:, :, irr), 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)

        forall (p = 1:nchan, k = 1:nstat)
            vw(k, p) = v(k) * moldat % wamp(p, k, irr)
        end forall

    end subroutine scale_boundary_amplitudes


    subroutine write_partial_wave_moments (moldat, M, suffix)

        type(MolecularData),      intent(in)    :: moldat
#ifdef WITH_COARRAYS
        complex(wp), allocatable, intent(inout) :: M(:, :, :)[:]
#else
        complex(wp), allocatable, intent(inout) :: M(:, :, :)
#endif
        character(len=*),         intent(in)    :: suffix

        integer :: i, irr, mgvn, u, ichan, ichlf, lf, mf
        character(len=100) :: filename

#ifdef WITH_COARRAYS
        ! gather coarray from all images to master
        sync all
        if (myproc /= 1) return
        do i = 2, nprocs
            m(:, :, :) = m(:, :, :) + m(:, :, :)[i]  ! this operation may require unlimited stack
        end do
#endif

        ! write to files
        do irr = 1, size(moldat % mgvns)
            mgvn = moldat % mgvns(irr)
            do ichan = 1, moldat % nchan(irr)
                if (any(m(ichan, :, mgvn) /= 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)') m(ichan, :, mgvn)
                    close (u)
                end if
            end do
        end do

    end subroutine write_partial_wave_moments


    subroutine write_cross_section (cs)

#ifdef WITH_COARRAYS
        real(wp), allocatable, intent(inout) :: cs(:, :)[:]
#else
        real(wp), allocatable, intent(inout) :: cs(:, :)
#endif
        integer :: u, i

#ifdef WITH_COARRAYS
        ! gather coarray from all images to master
        sync all
        if (myproc /= 1) return
        do i = 2, nprocs
            cs(:, :) = cs(:, :) + cs(:, :)[i]  ! this operation may require unlimited stack
        end do
#endif

        ! write to file
        open (newunit = u, file = 'gen_photo_xsec.txt', action = 'write', form = 'formatted')
        write (u, '(2E25.15)') cs
        close (u)

    end subroutine write_cross_section

end module multidip_io