pickOrbs_real_hubbard_single Subroutine

    subroutine pickOrbs_real_hubbard_single(ilut, nI, csf_i, excitInfo, pgen)
        ! write a specialized orbital picker for the real-space hubbard
        ! implementation, since we do not need all the symmetry stuff and
        ! we have to take the correct TMAT values for the ml-spin values
        integer(n_int), intent(in) :: ilut(0:nifguga)
        integer, intent(in) :: nI(nel)
        type(CSF_Info_t), intent(in) :: csf_i
        type(ExcitationInformation_t), intent(out) :: excitInfo
        real(dp), intent(out) :: pgen
        character(*), parameter :: this_routine = "pickOrbs_real_hubbard_single"

        integer :: elec, orb_i, orb_a
        real(dp) :: cum_arr(nSpatOrbs), elec_factor, cum_sum, r

        unused_var(ilut)

        ! first pick electron randomly:
        elec = 1 + floor(genrand_real2_dSFMT() * nel)

        orb_i = gtID(nI(elec))

        ! still use cum arrays to enable applying guga restrictions!
        ! but all orbitals are possible now of course

        select case (csf_i%stepvector(orb_i))

        case (1)
            ! i need a switch possibility for other d = 1 values
            call gen_cum_list_real_hub_1(csf_i, orb_i, cum_arr)

            elec_factor = 1.0_dp

        case (2)
            ! i need a switch possibility for other d = 2 values
            call gen_cum_list_real_hub_2(csf_i, orb_i, cum_arr)

            elec_factor = 1.0_dp

        case (3)
            ! no restrictions actually
            call gen_cum_list_real_hub_3(csf_i, orb_i, cum_arr)

            ! but twice the chance to have picked this spatial orbital:
            elec_factor = 2.0_dp

        case default
            call stop_all(this_routine, "should not have picked empty orb!")

        end select

        cum_sum = cum_arr(nSpatOrbs)

        if (near_zero(cum_sum)) then
            orb_a = 0
            excitInfo%valid = .false.
            return
        else
            r = genrand_real2_dSFMT() * cum_sum
            orb_a = binary_search_first_ge(cum_arr, r)

            if (orb_a == 1) then
                pgen = cum_arr(1) / cum_sum
            else
                pgen = (cum_arr(orb_a) - cum_arr(orb_a - 1)) / cum_sum

            end if
        end if

        ASSERT(orb_a /= orb_i)

        pgen = pgen * elec_factor / real(nel, dp)

        if (orb_a < orb_i) then
            ! raising generator
            excitInfo = assign_excitInfo_values_single(gen_type%R, orb_a, orb_i, orb_a, orb_i)

        else
            ! lowering generator
            excitInfo = assign_excitInfo_values_single(gen_type%L, orb_a, orb_i, orb_i, orb_a)

        end if

    end subroutine pickOrbs_real_hubbard_single