calc_1rdms_from_spinfree_2rdms Subroutine

    subroutine calc_1rdms_from_spinfree_2rdms(one_rdms, two_rdms, rdm_trace)

        ! For each 2-RDM in two_rdms, calculate the corresponding spin-free
        ! 1-RDM:
        !
        ! \gamma^{spinfree}_{p,q} = \frac{1}{N-1} \sum_a \Gamma^{spinfree}_{pa,qa}
        !
        ! Here, p, q and a are spatial labels. N is the number of electrons.

        ! The spinfree 1-RDM is defined in terms of the spinned 1-RDM by:
        !
        ! \gamma^{spinfree}_{p,q} = \gamma_{p\alpha,q\alpha) + \gamma_{p\beta,q\beta)

        ! IMPORTANT: This routine should *only* be used by taking *spin-free*
        ! 2-RDMs as the input. Specifically, it takes spin-free RDMs as returned
        ! by the create_spinfree_2rdm routine, which does *not* restrict the
        ! labels allowed. Inputting 2-RDMs in other will give incorrect results!

        ! The output 1-RDM elements are sorted in the standard form: elements
        ! are indexed using the SymLabelListInv_rot array, so that the 1-RDMs
        ! will be in block-diagonal form, with elements within each symmetry
        ! block stored together.

        use Parallel_neci, only: MPISumAll
        use RotateOrbsData, only: SymLabelListInv_rot
        use SystemData, only: nel
        use UMatCache, only: spatial

        type(one_rdm_t), intent(inout) :: one_rdms(:)
        type(rdm_list_t), intent(in) :: two_rdms
        real(dp), intent(in) :: rdm_trace(:)

        integer(int_rdm) :: pqrs
        integer :: pq, rs, p, q, r, s ! spatial orbitals
        integer :: ielem, irdm, ierr
        real(dp) :: rdm_sign(two_rdms%sign_length)
        real(dp), allocatable :: temp_rdm(:, :)

        do irdm = 1, size(one_rdms)
            one_rdms(irdm)%matrix = 0.0_dp
        end do

        ! Loop over all elements of the 2-RDM, \Gamma^{spinfree}_{pq,rs}, where
        ! p, q, r and s are spatial labels. If at least two spatial indices are
        ! the same then we have a contribution to the 1-RDM.
        do ielem = 1, two_rdms%nelements
            pqrs = two_rdms%elements(0, ielem)
            ! Obtain spin orbital labels and the RDM element.
            if (tGUGA) then
                call extract_2_rdm_ind(pqrs, p, q, r, s)
            else
                call calc_separate_rdm_labels(pqrs, pq, rs, r, s, q, p)
            end if

            call extract_sign_rdm(two_rdms%elements(:, ielem), rdm_sign)

            associate(ind => SymLabelListInv_rot)

                if (tGUGA) then
                    if (r == s) then
                        do irdm = 1, size(one_rdms)
                            one_rdms(irdm)%matrix(ind(p), ind(q)) = &
                                one_rdms(irdm)%matrix(ind(p), ind(q)) + rdm_sign(irdm) / 2.0_dp
                        end do
                    end if
                    ! if I count both, I do not need the factor 2..
                    if (p == q) then
                        do irdm = 1, size(one_rdms)
                            one_rdms(irdm)%matrix(ind(r), ind(s)) = &
                                one_rdms(irdm)%matrix(ind(r), ind(s)) + rdm_sign(irdm) / 2.0_dp
                        end do
                    end if
                else
                    ! An element of the form \Gamma_{pa,ra}.
                    if (q == s) then
                        do irdm = 1, size(one_rdms)
                            one_rdms(irdm)%matrix(ind(p), ind(r)) = &
                                one_rdms(irdm)%matrix(ind(p), ind(r)) + rdm_sign(irdm)
                        end do
                    end if
                end if
            end associate
        end do

        ! Allocate a temporary array in which to receive the MPI communication.
        allocate(temp_rdm(size(one_rdms(1)%matrix, 1), size(one_rdms(1)%matrix, 2)), stat=ierr)

        ! Perform a sum over all processes, for each 1-RDM being sampled.
        do irdm = 1, size(one_rdms)
            call MPISumAll(one_rdms(irdm)%matrix, temp_rdm)
            ! Copy summed RDM back to the main array, and normalise.
            one_rdms(irdm)%matrix = temp_rdm / (rdm_trace(irdm) * real(nel - 1, dp))
        end do

        deallocate(temp_rdm, stat=ierr)

    end subroutine calc_1rdms_from_spinfree_2rdms