apply_symmetries_for_output Subroutine

    subroutine apply_symmetries_for_output(rdm, rdm_recv, spawn, open_shell)

        ! This routine will take in rdm, and output a new rdm which will have
        ! all appropriate symmetries applied so that the latter RDM can be
        ! passed to the routine to write RDMs.

        ! The input RDM must already have hermiticy symmetry applied to it.

        ! WARNING: To clarify potential confusion, we point out that this
        ! routine also applies hermiticy again, but *only* for the spatial
        ! labels, not the full spin labels. It does so specifically using a
        ! particular legacy ordering. This is not a mistake - for open shell
        ! systems we do need to apply full hermiticy, but also need to apply spatial
        ! hermiticy because spatial labels are written within each output file.

        use rdm_data_utils, only: annihilate_rdm_list

        type(rdm_list_t), intent(in) :: rdm
        type(rdm_list_t), intent(inout) :: rdm_recv
        type(rdm_spawn_t), intent(inout) :: spawn
        logical, intent(in) :: open_shell

        integer(int_rdm) :: ijkl
        integer :: ielem
        integer :: ij, kl, i, j, k, l ! spin orbitals
        integer :: p, q, r, s ! spatial orbitals
        integer :: pq_legacy, rs_legacy ! spatial orbitals
        real(dp) :: rdm_sign(rdm%sign_length)
        logical :: nearly_full, finished, all_finished
#ifdef DEBUG_
        character(*), parameter :: this_routine = "apply_symmetries_for_output"
#endif
        ASSERT(.not. tGUGA)

        ! If we're about to fill up the spawn list, perform a communication.
        nearly_full = .false.
        ! Have we finished adding RDM elements to the spawned list?
        finished = .false.
        rdm_recv%nelements = 0

        do ielem = 1, rdm%nelements
            ! If the spawned list is nearly full, perform a communication.
            if (nearly_full) then
                call communicate_rdm_spawn_t_wrapper(spawn, rdm_recv, finished, all_finished)
                nearly_full = .false.
            end if

            ijkl = rdm%elements(0, ielem)
            ! Obtain spin orbital labels and the RDM element.
            call calc_separate_rdm_labels(ijkl, ij, kl, i, j, k, l)
            call extract_sign_rdm(rdm%elements(:, ielem), rdm_sign)

            ! When there are two elements which are guaranteed to be exactly the
            ! same, we usually only want to print one of them (and to average
            ! over the equal terms). This function returns the labels we want.
            call apply_legacy_output_ordering(i, j, k, l, rdm_sign, pq_legacy, rs_legacy)

            ! For closed shell systems, want bbbb -> aaaa, baba -> abab,
            ! baab -> abba, by flipping all spins, which is a symmetry for such
            ! systems. If the first label has beta spin then we definitely want
            ! to move this RDM element to the equivalent flipped term so go ahead
            ! and flip all the spins.
            if (.not. open_shell) then
                if (is_beta(i)) then
                    ! The ab_pair macro swaps alpha and beta spins of a label
                    ! while keeping the spatial orbital unchanged.
                    i = ab_pair(i)
                    j = ab_pair(j)
                    k = ab_pair(k)
                    l = ab_pair(l)
                end if
                ! If the spatial parts of i and j are the same, and the spatial
                ! parts of k and l are *also* the same, then the RDM element won't
                ! have been added into both equivalent spin-flipped arrays
                ! because i<j and k<l is enforced), so we don't count twice.
                if (.not. (is_in_pair(i, j) .and. is_in_pair(k, l))) then
                    ! Also, if (i,j) and (k,l) have the same spatial parts, but
                    ! different spin parts ((alpha,beta) and (beta,alpha), or
                    ! vice versa) then they only occur once, again because we
                    ! enforce i<j and k<l for all stored RDM elements.
                    if (.not. (pq_legacy == rs_legacy .and. (.not. ij == kl))) then
                        rdm_sign = rdm_sign * 0.5_dp
                    end if
                end if
            end if

            call add_to_rdm_spawn_t(spawn, i, j, k, l, rdm_sign, .false., nearly_full)

            if (open_shell) then
                ! For open shell systems, if i and j have the same spatial parts,
                ! and k and l do too, then we only have baba spin signature,
                ! (because we enforce i<j, k<l) but we'd like to print out abab too.
                if (is_in_pair(i, j) .and. is_in_pair(k, l)) then
                    call add_to_rdm_spawn_t(spawn, j, i, l, k, rdm_sign, .false., nearly_full)
                end if

                ! Because we enforce hermiticy symmetry in the output, we would
                ! only print the following terms with baab. We want to print it
                ! with abba too here, so do that.
                if (pq_legacy == rs_legacy .and. (.not. ij == kl)) then
                    call add_to_rdm_spawn_t(spawn, k, l, i, j, rdm_sign, .false., nearly_full)
                end if
            end if
        end do

        finished = .true.
        ! Keep performing communications until all RDM spawnings on every
        ! processor have been communicated.
        do
            call communicate_rdm_spawn_t_wrapper(spawn, rdm_recv, finished, all_finished)
            if (all_finished) exit
        end do
        call annihilate_rdm_list(rdm_recv)

    end subroutine apply_symmetries_for_output