find_open_shell_indices Function

    function find_open_shell_indices(nI) result(open_shells)
        ! find the indices of the open-shell orbitals of a nI rep.
        integer, intent(in) :: nI(:)
        integer, allocatable :: open_shells(:)

        integer(n_int) :: ilutI(0:niftot)
        integer :: n_open, i, cnt

        call EncodeBitDet(nI, ilutI)

        n_open = count_open_orbs(ilutI)

        allocate(open_shells(n_open))
        open_shells = 0

        cnt = 0
        do i = 1, size(nI)
            if (.not. IsDoub(ilutI, nI(i))) then
                cnt = cnt + 1
                open_shells(cnt) = i
            end if
        end do

    end function find_open_shell_indices