write_1rdm_hdf5 Subroutine

private subroutine write_1rdm_hdf5(parent, one_rdm)

Write the 1RDM to an HDF5 archive.

Arguments

Type IntentOptional Attributes Name
integer(kind=hid_t), intent(in) :: parent

HDF5 file handle of the parent directory.

real(kind=dp), intent(inout) :: one_rdm(:,:)

1RDM data redundant on each MPI rank


Contents

Source Code


Source Code

    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