calc_combined_rdm_label Subroutine

    pure subroutine calc_combined_rdm_label(i, j, k, l, ijkl, ij_out, kl_out)

        ! Combine the four 2-RDM orbital labels into unique integers. i and j
        ! are combined into one number, ij. k and l are combined into one
        ! number, kl. Both of these are then combined into one single
        ! number, ijkl. Assuming (i,j,k,l) are *spin* orbitals labels (which
        ! they usually will be but not necessarily), the largest value for ijkl
        ! is M^4, where M is the number of spin orbitals.

        ! The compression defined in this routine will not give a fully
        ! compressed RDM index labelling, because it allows a separate ij
        ! integer if i and j are equal, even though this RDM element is never
        ! accessed, and the same for k and l. It also doesn't take spatial
        ! symmetry into account. But this is fine if one just seeks a unique
        ! combined label for each combination of individual orbital labels.

        ! In: i, j, k, l - orbital labels for the RDM contribution.
        ! Out: ijkl - Label combining i, j, k, l.

        use SystemData, only: nbasis

        integer, intent(in) :: i, j, k, l ! spin or spatial orbitals
        integer(int_rdm), intent(out) :: ijkl
        integer, intent(out), optional :: ij_out, kl_out
        integer :: ij, kl

        ! maybe I need a change for the GUGA implementation, since
        ! we only need spatial orbitals here..
        ! todo
        ij = (i - 1) * nbasis + j
        kl = (k - 1) * nbasis + l
        ijkl = (ij - 1) * (int(nbasis, int_rdm)**2) + kl

        if (present(ij_out) .and. present(kl_out)) then
            ij_out = ij
            kl_out = kl
        end if

    end subroutine calc_combined_rdm_label