gen_excit_rs_hubbard Subroutine

    subroutine gen_excit_rs_hubbard(nI, ilutI, nJ, ilutJ, exFlag, ic, &
                                    ex, tParity, pGen, hel, store, run)
        !! An API interfacing function for generate_excitation to the rest of NECI:
        !!
        !! Requires guga_bitRepOps::current_csf_i to be set according to the ilutI.
        integer, intent(in) :: nI(nel), exFlag
        integer(n_int), intent(in) :: ilutI(0:NIfTot)
        integer, intent(out) :: nJ(nel), ic, ex(2, maxExcit)
        integer(n_int), intent(out) :: ilutJ(0:NifTot)
        real(dp), intent(out) :: pGen
        logical, intent(out) :: tParity
        HElement_t(dp), intent(out) :: hel
        type(excit_gen_store_type), intent(inout), target :: store
        integer, intent(in), optional :: run

        character(*), parameter :: this_routine = "gen_excit_rs_hubbard"

        integer :: ind, elec, src, orb
        integer, allocatable :: neighbors(:)
        real(dp), allocatable :: cum_arr(:)
        real(dp) :: cum_sum, p_elec, p_orb

        type(ExcitationInformation_t) :: excitInfo
        integer(n_int) :: ilutGj(0:nifguga)

        unused_var(exFlag)
        hel = h_cast(0.0_dp)
#ifdef WARNING_WORKAROUND_
        if (present(run)) then
            unused_var(run)
        end if
#endif
        unused_var(store)

        ASSERT(associated(lat))

        ic = 1
        ! i only have single excitations in the hubbard model
        ! the first plan is to choose an electron at random
        elec = 1 + int(genrand_real2_dsfmt() * nel)

        p_elec = 1.0_dp / real(nel, dp)
        ! and then from the neighbors of this electron we pick an empty
        ! spinorbital randomly, since all have the same matrix element
        src = nI(elec)

        ! now get neighbors
        neighbors = lat%get_spinorb_neighbors(src)

        call create_cum_list_rs_hubbard(ilutI, src, neighbors, cum_arr, cum_sum)

        if (cum_sum < EPS) then
            nJ(1) = 0
            pgen = 0.0_dp
            return
        end if

        call pick_from_cum_list(cum_arr, cum_sum, ind, p_orb)

        orb = neighbors(ind)

        pgen = p_elec * p_orb

        call make_single(nI, nJ, elec, orb, ex, tParity)

        ilutJ = make_ilutJ(ilutI, ex, 1)

        ! change for the mixed guga implementation
        if (tgen_guga_crude) then

            if (nJ(1) == 0) then
                pgen = 0.0_dp
                return
            end if

            call convert_ilut_toGUGA(ilutJ, ilutGj)

            if (.not. isProperCSF_ilut(ilutGJ, .true.)) then
                nJ(1) = 0
                pgen = 0.0_dp
            end if

            ASSERT(is_compatible(ilutI, current_csf_i))
            call calc_guga_matrix_element(ilutI, current_csf_i, ilutJ, CSF_Info_t(ilutJ), excitInfo, hel, .true.)

            if (abs(hel) < EPS) then
                nJ(1) = 0
                pgen = 0.0_dp
            end if

            global_excitinfo = excitInfo

            return
        end if

    end subroutine gen_excit_rs_hubbard