calc_rdm_prop Subroutine

    subroutine calc_rdm_prop(rdm, rdm_prop)

        use rdm_data, only: rdm_list_t
        use rdm_integral_fns, only: GetPropInts
        use SystemData, only: nel
        use LoggingData, only: iNumPropToEst

        type(rdm_list_t), intent(in) :: rdm
        real(dp), intent(out) :: rdm_prop(iNumPropToEst, rdm%sign_length)

        integer(int_rdm) :: ijkl
        integer :: ielem, iprop, ij, kl, i, j, k, l
        real(dp) :: rdm_sign(rdm%sign_length)

        rdm_prop = 0.0_dp

        ! Loop over all elements in the 2-RDM.
        do ielem = 1, rdm%nelements
            ijkl = rdm%elements(0, ielem)
            call extract_sign_rdm(rdm%elements(:, ielem), rdm_sign)

            ! Decode pqrs label into p, q, r and s labels.
            call calc_separate_rdm_labels(ijkl, ij, kl, i, j, k, l)

            do iprop = 1, iNumPropToEst
                ! We get dot product of the 1-RDM and one-electron property integrals:
                if (i == k) rdm_prop(iprop, :) = rdm_prop(iprop, :) + rdm_sign * GetPropInts(j, l, iprop) / (nel - 1)
                if (j == l) rdm_prop(iprop, :) = rdm_prop(iprop, :) + rdm_sign * GetPropInts(i, k, iprop) / (nel - 1)
                if (i == l) rdm_prop(iprop, :) = rdm_prop(iprop, :) - rdm_sign * GetPropInts(j, k, iprop) / (nel - 1)
                if (j == k) rdm_prop(iprop, :) = rdm_prop(iprop, :) - rdm_sign * GetPropInts(i, l, iprop) / (nel - 1)
            end do
        end do
    end subroutine calc_rdm_prop