calc_1rdms_from_2rdms Subroutine

    subroutine calc_1rdms_from_2rdms(rdm_defs, one_rdms, two_rdms, rdm_trace, open_shell)

        ! For each 2-RDM in two_rdms, if open_shell is true then calculate the
        ! full spinned 1-RDM, otherwise calculate the spinfree 1-RDM. The
        ! former case is defined by:
        !
        ! \gamma_{i,j} = \frac{1}{N-1} \sum_k \Gamma_{ik,jk}
        !
        ! Here, i, j and k are spatial labels. N is the number of electrons.
        !
        ! The spinfree case is then a contraction over the spin labels of the
        ! spinned 1-RDM:
        !
        ! \gamma^{spinfree}_{p,q} = \gamma_{p\alpha,q\alpha) + \gamma_{p\beta,q\beta)
        !
        ! where p and q are spatial labels.

        ! 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 rdm_data, only: rdm_definitions_t
        use RotateOrbsData, only: SymLabelListInv_rot
        use SystemData, only: nel
        use UMatCache, only: spatial

        type(rdm_definitions_t), intent(in) :: rdm_defs
        type(one_rdm_t), intent(inout) :: one_rdms(:)
        type(rdm_list_t), intent(in) :: two_rdms
        real(dp), intent(in) :: rdm_trace(:)
        logical, intent(in) :: open_shell

        integer(int_rdm) :: ijkl
        integer :: ielem, irdm, ierr
        integer :: ij, kl, i, j, k, l ! spin orbitals
        integer :: p, q, r, s ! spatial orbitals
        real(dp) :: rdm_sign(two_rdms%sign_length)
        real(dp), allocatable :: temp_rdm(:, :)
        logical :: is_transition_rdm
#ifdef DEBUG_
        character(*), parameter :: this_routine = "calc_1rdms_from_2rdms"
#endif
        ASSERT(.not. tGUGA)

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

        ! Loop over all elements of the 2-RDM, \Gamma_{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
            ijkl = two_rdms%elements(0, ielem)
            ! Obtain spin orbital labels and the RDM element.
            call calc_separate_rdm_labels(ijkl, ij, kl, i, j, k, l)

            ! For closed shell systems we work with spatial orbitals, to
            ! calculate spin-free 1RDMs.
            if (open_shell) then
                p = i; q = j;
                r = k; s = l;
            else
                p = spatial(i); q = spatial(j);
                r = spatial(k); s = spatial(l);
            end if

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

            ! If abba or baab term - swap last two indices and sign.
            if (.not. same_spin(i, k)) then
                if (open_shell) then
                    r = l; s = k;
                else
                    r = spatial(l); s = spatial(k);
                end if
                rdm_sign = -rdm_sign
            end if

            associate(ind => SymLabelListInv_rot)

                ! An element of the form \Gamma_{aq,as}.
                if (p == r) then
                    do irdm = 1, size(one_rdms)
                        one_rdms(irdm)%matrix(ind(q), ind(s)) = one_rdms(irdm)%matrix(ind(q), ind(s)) + rdm_sign(irdm)
                    end do
                end if
                ! 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

                ! The below cases give contributions by swapping one pair of
                ! indices. Only include these contributions if we have aaaa or
                ! bbbb terms. This because if we had a term with spin signature
                ! abab (for example), then swapping as below would give abba
                ! or baab terms, which don't contribute to the 1-RDM.
                if (same_spin(k, l)) then
                    ! An element of the form \Gamma_{pa,as}.
                    if (p == s) then
                        do irdm = 1, size(one_rdms)
                            one_rdms(irdm)%matrix(ind(q), ind(r)) = one_rdms(irdm)%matrix(ind(q), ind(r)) - rdm_sign(irdm)
                        end do
                    end if
                    ! An element of the form \Gamma_{aq,ra}.
                    if (q == r) then
                        do irdm = 1, size(one_rdms)
                            one_rdms(irdm)%matrix(ind(p), ind(s)) = one_rdms(irdm)%matrix(ind(p), ind(s)) - rdm_sign(irdm)
                        end do
                    end if
                end if

            end associate
        end do

        ! Allocate a temporary RDM array.
        allocate(temp_rdm(size(one_rdms(1)%matrix, 1), size(one_rdms(1)%matrix, 2)), stat=ierr)

        ! Make non-transition RDMs symmetric.
        do irdm = 1, size(one_rdms)
            is_transition_rdm = rdm_defs%state_labels(1, irdm) /= rdm_defs%state_labels(2, irdm)
            if (.not. is_transition_rdm) then
                ! Use temp_rdm as temporary space for the transpose, to (hopefully)
                ! prevent a temporary array being created in the sum below.
                temp_rdm = transpose(one_rdms(irdm)%matrix)
                one_rdms(irdm)%matrix = (one_rdms(irdm)%matrix + temp_rdm) / 2.0_dp
            end if
        end do

        ! 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_2rdms