rdm_hdf5.F90 Source File


Contents

Source Code


Source Code

#include "macros.h"

module rdm_hdf5

    use MPI_wrapper, only: mpiInfoNull, CommGlobal
    use Parallel_neci, only: iProcIndex, stop_all, MPIBCast
    use constants, only: dp, stdout, hdf_err
#ifdef USE_HDF_
    use parallel_hdf5_utils, only: read_data_phdf5, write_data_phdf5
    use hdf5, only: h5open_f, h5pset_fapl_mpio_f, h5gcreate_f, h5pclose_f, &
                    h5close_f, h5garbage_collect_f, h5pcreate_f, H5P_FILE_ACCESS_F, &
                    hid_t, hsize_t, h5fclose_f, h5gclose_f, H5F_ACC_TRUNC_F, &
                    h5fcreate_f
#endif
    use fortran_strings, only: str

    implicit none
    private
    public :: write_rdms_hdf5

contains

    !> Write all RDMs specified in the input to an HDF5 archive.
    subroutine write_rdms_hdf5(rdm_defs, rdm, rdm_trace, one_rdms)
        use rdm_data, only: rdm_definitions_t, rdm_list_t, one_rdm_t
        use LoggingData, only: tWriteSpinFreeRDM, tPrint1RDM
        !> Type carrying information such as number and type of RDM (regular/transition)
        type(rdm_definitions_t), intent(in) :: rdm_defs
        !> 2RDM data distributed over all MPI ranks
        type(rdm_list_t), intent(in) :: rdm
        !> Normalisation of the 2RDM
        real(dp), intent(in) :: rdm_trace(rdm%sign_length)
        !> Array carrying the 1RDM belonging to each root
        type(one_rdm_t), intent(inout), optional :: one_rdms(:)
        !> Required for the stop all to check for HDF5 library
        character(*), parameter :: this_routine = 'write_hdf5_rdms'
#ifdef USE_HDF_
        !> File and group handles
        integer(hid_t) :: file_id, root_id, rdm_id, plist_id
        !> HDF5 error code
        integer(hdf_err) :: err
        !> Name of the HDF5 archive
        character(255) :: filename
        !> The iroot'th eigenvector of Hamiltonian, e.g. 1 for ground state,
        !> 2 for first excited state, ...
        integer :: iroot

        if (rdm_defs%nrdms_transition > 0) then
            call stop_all(this_routine, "Transition RDM support pending.")
        end if

        do iroot = 1, rdm_defs%nrdms
            if (iProcIndex == 0) then
                filename = 'fciqmc.rdms.' // str(iroot) // '.h5'
                write(stdout, *) "============== Writing HDF5 RDMs =============="
                write(stdout, *) "File name: ", trim(filename)
                write(stdout, *) "Regular RDMs saved in /archive/rdms/AA00/ where &
                                     &A denotes the number of fermion operators."
            end if
            call MPIBCast(filename)

            ! All HDF5 APIs that modify structural metadata are collective
            call h5open_f(err)
            call h5pcreate_f(H5P_FILE_ACCESS_F, plist_id, err)
            ! setup file access property list for MPI-IO access.
            call h5pset_fapl_mpio_f(plist_id, CommGlobal, mpiInfoNull, err)
            ! create the file collectively, H5F_ACC_TRUNC_F = overwrite existing file
            call h5fcreate_f(filename, H5F_ACC_TRUNC_F, file_id, err, access_prp=plist_id)
            call h5pclose_f(plist_id, err)
            call h5gcreate_f(file_id, 'archive', root_id, err)
            call h5gcreate_f(root_id, 'rdms', rdm_id, err)

            if (tPrint1RDM) then
                call write_1rdm_hdf5(rdm_id, one_rdms(iroot)%matrix)
                write(stdout, *) "1RDM written to file."
            end if
            if (tWriteSpinFreeRDM) then
                call write_2rdm_hdf5(rdm_id, rdm, rdm_trace, iroot)
                write(stdout, *) "2RDM written to file."
            end if

            if (iProcIndex == 0) write(stdout, *) "closing RDM file."
            call h5gclose_f(rdm_id, err)
            call h5gclose_f(root_id, err)
            call h5fclose_f(file_id, err)
            call h5close_f(err)
            call h5garbage_collect_f(err)
            if (iProcIndex == 0) write(stdout, *) "RDM file write successful."
        end do
#else
      unused_var(rdm_defs)
      unused_var(rdm)
      unused_var(rdm_trace)
      unused_var(one_rdms)
      call stop_all(this_routine, 'HDF5 support not enabled at compile time.')
#endif
    end subroutine write_rdms_hdf5

    !> Write the 1RDM to an HDF5 archive.
    subroutine write_1rdm_hdf5(parent, one_rdm)
        use RotateOrbsData, only: ind => SymLabelListInv_rot
#ifndef USE_HDF_
        !> HDF5 file handle of the parent directory.
        integer, intent(in) :: parent
        !> 1RDM data redundant on each MPI rank
        real(dp), intent(inout) :: one_rdm(:, :)
        !> Required for the stop all to check for HDF5 library
        character(*), parameter :: this_routine = 'write_1rdm_hdf5'
#endif
#ifdef USE_HDF_
        !> HDF5 file handle of the parent directory.
        integer(hid_t), intent(in) :: parent
        !> 1RDM data redundant on each MPI rank
        real(dp), intent(inout) :: one_rdm(:, :)
        !> Used to accumulate the RDM indices before writing
        integer, allocatable :: index(:), indices(:,:)
        !> Accumulate the non-zero RDM values respectively
        real(dp), allocatable :: values(:)
        !> Spatial orbital indices
        integer :: p, q
        !> HDF5 file handle of 1RDM "1100" group
        integer(hid_t) :: grp_id
        !> HDF5 error code
        integer(hdf_err) :: err

        values = [real(dp) ::]
        index = [integer ::]
        ! All the 1RDM values are stored on all ranks redundantly anyway.
        ! Here the values and indices are only accumulated on the head node,
        ! to prevent writing the same data multiple times to the HDF5 archive.
        if (iProcIndex == 0) then
            do p = 1, size(one_rdm, dim=1)
                do q = 1, size(one_rdm, dim=1)
                    if (p >= q) then
                        if (abs(one_rdm(ind(p), ind(q))) > 1.e-12_dp) then
                          index = [index, p]; index = [index, q]
                          values = [values, one_rdm(ind(p), ind(q))]
                        end if
                    end if
                end do
            end do
        end if
        indices = reshape(index, [2, size(values)])

        call h5gcreate_f(parent, '1100', grp_id, err)
        call write_data_phdf5(values, 'values', grp_id)
        call write_data_phdf5(indices, 'indices', grp_id)
        call h5gclose_f(grp_id, err)
#else
      unused_var(parent)
      unused_var(one_rdm)
      call stop_all(this_routine, 'HDF5 support not enabled at compile time.')
#endif
    end subroutine write_1rdm_hdf5

    !> Write the 2RDM to an HDF5 archive.
    subroutine write_2rdm_hdf5(parent, rdm, rdm_trace, iroot)
        use rdm_data, only: rdm_list_t, int_rdm
        use rdm_data_utils, only: calc_separate_rdm_labels, extract_sign_rdm
#ifndef USE_HDF_
        !> HDF5 file handle of the parent directory.
        integer, intent(in) :: parent
        !> 2RDM data distributed over all MPI ranks
        type(rdm_list_t), intent(in) :: rdm
        !> Normalisation of the 2RDM
        real(dp), intent(in) :: rdm_trace(rdm%sign_length)
        !> The iroot'th eigenvector of Hamiltonian, e.g. 1 for ground state,
        !> 2 for first excited state, ...
        integer, intent(in) :: iroot
        character(*), parameter :: this_routine = 'write_2rdm_hdf5'
#endif
#ifdef USE_HDF_
        !> HDF5 file handle of the parent directory.
        integer(hid_t), intent(in) :: parent
        !> 2RDM data distributed over all MPI ranks
        type(rdm_list_t), intent(in) :: rdm
        !> Normalisation of the 2RDM
        real(dp), intent(in) :: rdm_trace(rdm%sign_length)
        !> The 2RDM storage with symmetry introduces a sign that has to kept
        !> track of, refer to "rdm_data.F90" for further information
        real(dp) :: rdm_sign(rdm%sign_length)
        !> The iroot'th eigenvector of Hamiltonian, e.g. 1 for ground state,
        !> 2 for first excited state, ...
        integer, intent(in) :: iroot
        !> HDF5 file handle of 1RDM "1100" group
        integer(hid_t) :: grp_id
        !> HDF5 error code
        integer(hdf_err) :: err
        !> Folded 2RDM loop index, refer to stochastic GUGA-CASSCF paper 2021 SI
        integer(int_rdm) :: pqrs
        !> Spatial orbital and loop indices
        integer :: ielem, p, q, r, s, pq, rs
        !> Used to accumulate the RDM indices before writing
        integer, allocatable :: index(:), indices(:,:)
        !> Accumulate the non-zero RDM values respectively
        real(dp), allocatable :: values(:)

        ! create dynamic "indices" and "value" arrays and append to them in loop
        values = [real(dp) ::]
        index = [integer ::]
        do ielem = 1, rdm%nelements
            pqrs = rdm%elements(0, ielem)
            call extract_sign_rdm(rdm%elements(:, ielem), rdm_sign)
            rdm_sign = rdm_sign / rdm_trace
            call calc_separate_rdm_labels(pqrs, pq, rs, p, s, q, r)
            if (abs(rdm_sign(iroot)) > 1.e-12_dp) then
                if (p >= q .and. pq >= rs .and. p >= r .and. p >= s) then
                    index = [index, p]; index = [index, q]
                    index = [index, r]; index = [index, s]
                    values = [values, rdm_sign(iroot)]
                end if
            end if
        end do
        indices = reshape(index, [4, size(values)])

        call h5gcreate_f(parent, '2200', grp_id, err)
        call write_data_phdf5(values, 'values', grp_id)
        call write_data_phdf5(indices, 'indices', grp_id)
        call h5gclose_f(grp_id, err)
#else
      unused_var(parent)
      unused_var(rdm)
      unused_var(rdm_trace)
      unused_var(iroot)
      call stop_all(this_routine, 'HDF5 support not enabled at compile time.')
#endif
    end subroutine write_2rdm_hdf5

end module rdm_hdf5