gen_all_singles_rs_hub_default Subroutine

    subroutine gen_all_singles_rs_hub_default(nI, n_excits, det_list, sign_list)
        integer, intent(in) :: nI(nel)
        integer, intent(out) :: n_excits
        integer(n_int), intent(out), allocatable :: det_list(:, :)
        real(dp), intent(out), allocatable, optional :: sign_list(:)
#if defined(DEBUG_) || defined(CMPLX_)
        character(*), parameter :: this_routine = "gen_all_singles_rs_hub_default"
#endif
        integer(n_int) :: ilut(0:NIfTot), ilutJ(0:NIfTot)
        integer :: n_bound, i, src, j, neigh, ex(2, 2), pos, nJ(nel)
        integer(n_int), allocatable :: temp_list(:, :)
        integer, allocatable :: neighbors(:)
        HElement_t(dp) :: elem
        real(dp), allocatable :: temp_sign(:)
        logical :: t_sign, tpar

        ASSERT(associated(lat))

        call EncodeBitDet(nI, ilut)

        n_excits = 1

        n_bound = nel * (nbasis - nel)
        allocate(temp_list(0:NIfTot, n_bound))
        temp_list = 0_n_int

        if (present(sign_list)) then
            t_sign = .true.
            allocate(temp_sign(n_bound), source=0.0_dp)
        else
            t_sign = .false.
        end if

        ex = 0
        ! loop over electrons
        do i = 1, nel
            src = nI(i)
            ex(1, 1) = src
            neighbors = lat%get_spinorb_neighbors(src)

            ! and loop over the neighboring sites of this electron
            do j = 1, size(neighbors)

                neigh = neighbors(j)

                ex(2, 1) = neigh

                ASSERT(same_spin(src, neigh))

                ! if it is not occupied it should be a possible excitation
                if (IsNotOcc(ilut, neighbors(j))) then
                    ! but to be sure, check the matrix element:
                    ! but use the lattice get_helement_lattice to
                    ! avoid circular dependencies
                    call make_single(nI, nJ, i, neighbors(j), ex, tpar)
                    elem = get_helement_lattice(nI, nJ)

                    if (abs(elem) > EPS) then

                        ilutJ = make_ilutJ(ilut, ex, 1)

                        ! actually a search is not really necessary.. since
                        ! all the single excitations are unique.. but
                        ! just to be sure
                        pos = binary_search_ilut(temp_list(:, 1:n_excits), ilutJ, nifd + 1)

                        if (pos < 0) then

                            temp_list(:, n_excits) = ilutJ

                            if (t_sign) then
#ifdef CMPLX_
                                call stop_all(this_routine, "not implemented for complex")
#else
                                temp_sign(n_excits) = sign(1.0_dp, elem)
                                call sort(temp_list(:, 1:n_excits), temp_sign(1:n_excits))
#endif
                            else
                                call sort(temp_list(:, 1:n_excits), ilut_lt, ilut_gt)
                            end if
                            n_excits = n_excits + 1
                            ! damn.. i have to sort everytime i guess..
                        end if
                    end if
                end if
            end do
        end do

        n_excits = n_excits - 1
        allocate(det_list(0:NIfTot, n_excits), source=temp_list(:, 1:n_excits))

        if (t_sign) then
            allocate(sign_list(n_excits), source=temp_sign(1:n_excits))

            call sort(det_list, sign_list)
        else
            call sort(det_list, ilut_lt, ilut_gt)
        end if

    end subroutine gen_all_singles_rs_hub_default