write_2rdm_hdf5 Subroutine

    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