pick_orbitals_guga_tJ Subroutine

    subroutine pick_orbitals_guga_tJ(ilut, nI, csf_i, excitInfo, orb_pgen)
        ! orbital picking routine for the GUGA t-J model.
        ! the most effective way would be to pick a hole first, instead
        ! of an random electron, so the chance of a succesful hop
        ! increases.. but I might be too lazy for now to do that
        ! actually.
        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) :: orb_pgen

        integer :: elec, id, ind, tgt
        real(dp) :: p_elec, cum_sum, r, p_orb
        integer, allocatable :: neighbors(:)
        real(dp), allocatable :: cum_arr(:)

        unused_var(ilut)

        elec = 1 + int(genrand_real2_dsfmt() * nel)
        p_elec = 1.0_dp / real(nel, dp)

        id = gtID(nI(elec))

        neighbors = lat%get_neighbors(lat%get_site_index(id))

        call gen_guga_tJ_cum_list(csf_i, id, cum_arr)

        cum_sum = cum_arr(size(neighbors))

        if (cum_sum < EPS) then
            excitInfo%valid = .false.
            orb_pgen = 0.0_dp
            return
        end if

        r = genrand_real2_dsfmt() * cum_sum
        ind = binary_search_first_ge(cum_arr, r)

        tgt = neighbors(ind)

        if (ind == 1) then
            p_orb = cum_arr(1) / cum_sum
        else
            p_orb = (cum_arr(ind) - cum_arr(ind - 1)) / cum_sum
        end if

        if (tgt < id) then
            excitInfo = assign_excitInfo_values_single(gen_type%R, tgt, id, tgt, id)
        else
            excitInfo = assign_excitInfo_values_single(gen_type%L, tgt, id, id, tgt)
        end if

        orb_pgen = p_elec * p_orb

    end subroutine pick_orbitals_guga_tJ