apply_legacy_output_ordering Subroutine

    pure subroutine apply_legacy_output_ordering(i, j, k, l, rdm_sign, pq_legacy, rs_legacy)

        ! Enforce the symmetries of RDMs to only keep certain combinations of
        ! i, j, k, l spin labels, where a redundancy exists. Whenever we have
        ! an unused combination, flip/swap labels (and the sign if necessary).

        ! For example the 2-RDM is hermitian, so if ij /= kl, then we only need
        ! to print either \Gamma_{ij,kl} or \Gamma{kl,ij}, but not both. Which
        ! combinations we decide to print is decided below, which is purely a
        ! legacy decision (as far as I know!). See comments below for definitions
        ! of what we keep.

        use SystemData, only: nbasis
        use UMatCache, only: spatial

        integer, intent(inout) :: i, j, k, l ! spin orbitals
        real(dp), intent(inout) :: rdm_sign(:)
        integer, intent(out) :: pq_legacy, rs_legacy ! spatial orbitals

        integer :: i_temp, j_temp ! spin orbitals
        integer :: p, q, r, s ! spatial_orbitals

        ! RDMs are output in files labelled by their spin signatures:
        ! aaaa, abab, abba, bbbb, baba or baab.
        ! Within each file, therefore, only spatial orbital labels are printed.
        ! Thus, we need to use spatial orbitals to determine which RDM elements
        ! are  to kept, and which transformed.
        p = spatial(i); q = spatial(j);
        r = spatial(k); s = spatial(l);
        ! When we calculate the combined labels, pq and rs, we would
        ! usually have p and q swapped below, and similarly with r and s.
        ! However, the old RDM files prints only RDM elements with pq < rs,
        ! where pq and rs are defined as follows.
        pq_legacy = (q - 1) * nbasis + p
        rs_legacy = (s - 1) * nbasis + r

        ! Apply symmetry (for *real* RDMs), to only print elements from one
        ! half of the RDM, using the legacy ordering.
        if (pq_legacy > rs_legacy) then
            i_temp = i; j_temp = j;
            i = k; j = l;
            k = i_temp; l = j_temp;
        end if

        ! If either i and j have the same spatial part, of k and l have the
        ! same spatial part, and we have a spin signature with 2 alphas and
        ! 2 betas, then the convention is to output it as either abab or
        ! baba, but *not* as abba or baab. If we have abba or baab in this
        ! case then we have to swap two indices and introduce a minus sign.
        ! Because we enforce i<j and k<l in all RDM elements, there are only
        ! two possibilities to consider:
        if (is_in_pair(i, j) .and. is_beta(i) .and. is_alpha(j) .and. &
            is_alpha(k) .and. is_beta(l)) then
            i = ab_pair(i)
            j = ab_pair(j)
            rdm_sign = -rdm_sign
        else if (is_in_pair(k, l) .and. is_alpha(i) .and. is_beta(j) .and. &
                 is_beta(k) .and. is_alpha(l)) then
            k = ab_pair(k)
            l = ab_pair(l)
            rdm_sign = -rdm_sign
        end if

    end subroutine apply_legacy_output_ordering