extract_2_rdm_ind Subroutine

    pure subroutine extract_2_rdm_ind(ijkl, i, j, k, l, &
                                      ij_out, kl_out, excit_lvl, excit_typ)
        ! the inverse routine of the function above.
        ! it is actually practical to have ij and kl also available at
        ! times, since it can identify diagonal entries of the two-RDM
        integer(int_rdm), intent(in) :: ijkl
        integer, intent(out) :: i, j, k, l
        integer(int_rdm), intent(out), optional :: ij_out, kl_out
        integer, intent(out), optional :: excit_lvl, excit_typ

        integer(int_rdm) :: ij, kl

        integer(int_rdm) :: ijkl_

        ijkl_ = iand(ijkl, rdm_ind_bitmask)

        kl = mod(ijkl_ - 1, int(nSpatOrbs, int_rdm)**2) + 1
        ij = (ijkl_ - kl) / (nSpatOrbs**2) + 1

        call extract_1_rdm_ind(ij, i, j)
        call extract_1_rdm_ind(kl, k, l)

        if (present(ij_out)) ij_out = ij
        if (present(kl_out)) kl_out = kl

        if (present(excit_lvl)) then
            excit_lvl = extract_excit_lvl_rdm(ijkl)
        end if
        if (present(excit_typ)) then
            excit_typ = extract_excit_type_rdm(ijkl)
        end if

    end subroutine extract_2_rdm_ind