fill_molcas_rdms Subroutine

    subroutine fill_molcas_rdms(rdm, rdm_trace, irdm, psmat, pamat, dmat)

        !! Populate the Molcas RDM arrays PSMAT/PAMAT/DMAT.

        type(rdm_list_t), intent(in) :: rdm
            !! Stores RDMs as 1D lists whose elements can be accessed through
            !! a has table.
        real(dp), intent(in) :: rdm_trace(rdm%sign_length)
            !! Trace of RDMs required for normalisation of sampled arrays.
        integer, intent(in) :: irdm
            !! loop index over the different states to be sampled;
            !! 2 replicas are required per state.
        real(dp), intent(out), allocatable :: psmat(:), pamat(:), dmat(:)
            !! Molcas RDM arrays to be populated.

        integer :: n_one_rdm, n_two_rdm, ielem
        integer :: pqrs_m, pq_m, rs_m, p, q, r, s, p_m, q_m, r_m, s_m
        integer(int_rdm) :: pqrs
        real(dp) :: rdm_sign(rdm%sign_length), rdm_sign_
        real(dp), allocatable :: dmat_loc(:), psmat_loc(:), pamat_loc(:)

        n_one_rdm = nSpatorbs * (nSpatorbs + 1) .div. 2
        n_two_rdm = n_one_rdm * (n_one_rdm + 1) .div. 2

        allocate(dmat_loc(n_one_rdm), source=0.0_dp)
        allocate(psmat_loc(n_two_rdm), source=0.0_dp)
        allocate(pamat_loc(n_two_rdm), source=0.0_dp)

        do ielem = 1, rdm%nelements
            pqrs = rdm%elements(0, ielem)
            call extract_2_rdm_ind(pqrs, p, q, r, s)
            pqrs_m = contract_molcas_2_rdm_index(p, q, r, s)
            call extract_molcas_2_rdm_index(pqrs_m, &
                                            p_m, q_m, r_m, s_m, pq_m, rs_m)

            call extract_sign_rdm(rdm%elements(:, ielem), rdm_sign)
            rdm_sign_ = rdm_sign(irdm) / rdm_trace(irdm)

            ! now make the fill logic
            ! the molcas rdm elements are given by
            ! if r /= s (and probably here p /= q)
            ! psmat_loc(pqrs) = (two_rdm(pqrs) + two_rdm(pqsr)) / 2
            ! pamat_loc(pqrs) = (two_rdm(pqrs) - two_rdm(pqsr)) / 2
            ! if r == s (and probably .or. p == q)
            ! psmat_loc(pqrs) = 2 * two_rdm(pqrs)
            ! es geht eigentlich nur drum wann das element
            ! negativ zur anti-symmetrischen beiträgt..
            ! und wenn es nur zur diagonalen beiträgt..

            if (pq_m == rs_m) then
                if (p_m == q_m) then
                    psmat_loc(pqrs_m) = psmat_loc(pqrs_m) + rdm_sign_ / 2.0_dp
                else
                    psmat_loc(pqrs_m) = psmat_loc(pqrs_m) + rdm_sign_ / 4.0_dp
                    pamat_loc(pqrs_m) = pamat_loc(pqrs_m) + &
                                        molcas_sign(p, q, r, s) * rdm_sign_ / 4.0_dp
                end if
            else
                if (p_m == q_m) then
                    psmat_loc(pqrs_m) = psmat_loc(pqrs_m) + rdm_sign_ / 4.0_dp
                else
                    psmat_loc(pqrs_m) = psmat_loc(pqrs_m) + rdm_sign_ / 8.0_dp
                    if (r_m /= s_m) then
                        pamat_loc(pqrs_m) = pamat_loc(pqrs_m) + &
                                            molcas_sign(p, q, r, s) * rdm_sign_ / 8.0_dp
                    end if
                end if
            end if

            pq_m = contract_molcas_1_rdm_index(p, q)
            rs_m = contract_molcas_1_rdm_index(r, s)
            ! convert to 1-RDM
            if (r == s .and. p == q) then
                dmat_loc(pq_m) = dmat_loc(pq_m) + rdm_sign_
                dmat_loc(rs_m) = dmat_loc(rs_m) + rdm_sign_
            else
                if (r == s) then
                    dmat_loc(pq_m) = dmat_loc(pq_m) + rdm_sign_ / 2.0_dp
                else if (p == q) then
                    dmat_loc(rs_m) = dmat_loc(rs_m) + rdm_sign_ / 2.0_dp
                end if
            end if
        end do

        dmat_loc(:) = dmat_loc(:) / (2.0_dp * real(nel - 1, dp))

        allocate(dmat(n_one_rdm), source=0.0_dp)
        allocate(psmat(n_two_rdm), source=0.0_dp)
        allocate(pamat(n_two_rdm), source=0.0_dp)

        call MPISumAll(dmat_loc, dmat)
        call MPISumAll(psmat_loc, psmat)
        call MPISumAll(pamat_loc, pamat)

    end subroutine fill_molcas_rdms