gen_excit_rs_hubbard_spin_dependent_transcorr Subroutine

    subroutine gen_excit_rs_hubbard_spin_dependent_transcorr(nI, ilutI, nJ, ilutJ, exFlag, ic, &
                                                             ex, tParity, pGen, hel, store, run)
        ! new excitation generator for the real-space hubbard model with
        ! the hopping transcorrelation, which leads to double excitations
        ! and long-range single excitations in the real-space hubbard..
        ! this complicates things alot!
        ! this is the specific implementation for spin-dependent
        ! trans-correlation

        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_transcorr"

        integer :: ind, elec, src, orb
        real(dp) :: cum_arr_t(nBasis / 2)
        ! i have to resolve this conflict:
        real(dp), allocatable :: cum_arr_o(:)
        integer, allocatable :: neighbors(:)
        real(dp) :: cum_sum, p_elec, p_orb

        unused_var(exFlag)
        unused_var(run)
        unused_var(store)

        ilutJ = 0_n_int
        ic = 0
        nJ(1) = 0
        hel = h_cast(0.0_dp)

#ifdef WARNING_WORKAROUND_
        hel = 0.0_dp
        if (present(run)) then
            unused_var(run)
        end if
#endif
        unused_var(store)
        ASSERT(associated(lat))

        ic = 1

        ! pick the electron randomly

        ! still 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)

        ! and for alpha-electrons we have trans-correlation
        if (is_alpha(src)) then
            ! now we can have more than only nearest neighbor hopping!
            ! so implement a new cum-list creator
            call create_cum_list_rs_hubbard_transcorr_single(nI, ilutI, src, cum_arr_t, cum_sum)
        else
            ! only hopping to neighbors allowed

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

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

        end if

        ! the rest stays the same i guess..
        if (cum_sum < EPS) then
            nJ(1) = 0
            pgen = 0.0_dp
            return
        end if

        if (is_alpha(src)) then
            call pick_from_cum_list(cum_arr_t, cum_sum, ind, p_orb)
            ! we know it is alpha
            orb = 2 * ind
        else
            call pick_from_cum_list(cum_arr_o, cum_sum, ind, p_orb)
            orb = neighbors(ind)
        end if

        pgen = p_elec * p_orb

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

        ilutJ = make_ilutJ(ilutI, ex, 1)

    end subroutine gen_excit_rs_hubbard_spin_dependent_transcorr