calc_separate_rdm_labels Subroutine

    pure subroutine calc_separate_rdm_labels(ijkl, ij, kl, i, j, k, l)

        ! Decode the four orbital labels stored in the input ijkl, i.e. do the
        ! opposite of calc_combined_rdm_label.

        use SystemData, only: nbasis

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

        kl = int(mod(ijkl - 1, int(nbasis, int_rdm)**2)) + 1
        ij = int((ijkl - int(kl, int_rdm)) / int(nbasis, int_rdm)**2) + 1

        j = mod(ij - 1, nbasis) + 1
        i = (ij - j) / nbasis + 1

        l = mod(kl - 1, nbasis) + 1
        k = (kl - l) / nbasis + 1

    end subroutine calc_separate_rdm_labels