gen_excit_rs_hubbard_transcorr Subroutine

    subroutine gen_excit_rs_hubbard_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!

        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(nBasis / 2)
        real(dp) :: cum_sum, p_elec, p_orb
#ifdef DEBUG_
        real(dp) :: temp_pgen
#endif

        unused_var(exFlag)
        unused_var(store)

        ilutJ = 0_n_int
        ic = 0
        nJ(1) = 0
        hel = h_cast(0.0_dp)
#ifdef WARNING_WORKAROUND_
        if (present(run)) then
            unused_var(run)
        end if
#endif

        ASSERT(associated(lat))

        if (genrand_real2_dsfmt() < pDoubles) then
            ic = 2

            call gen_double_excit_rs_hub_transcorr(nI, ilutI, nJ, ilutJ, ex, tParity, pgen)

            pgen = pDoubles * pgen

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

        else
            ic = 1

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

            ! 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, cum_sum)

            ! the rest stays the same i guess..
            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)

            ! all orbitals are possible i guess, so make cum_arr for all
            ! orbitals as ind already. we "just" have to fix the spin
            if (is_beta(src)) then
                orb = 2 * ind - 1
            else
                orb = 2 * ind
            end if

            pgen = p_elec * p_orb * (1.0_dp - pDoubles)

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

        end if

        ilutJ = make_ilutJ(ilutI, ex, ic)

#ifdef DEBUG_
        temp_pgen = calc_pgen_rs_hubbard_transcorr(nI, ilutI, ex, ic)
        if (abs(pgen - temp_pgen) > EPS) then
            root_print  "calculated pgen differ for exitation: "
            root_print  "nI: ", nI
            root_print  "ex: ", ex
            root_print  "ic: ", ic
            root_print  "pgen: ", pgen
            root_print  "calc. pgen: ", temp_pgen
            root_print  "H_ij: ", get_helement_lattice(nI, nJ, ic)
        end if
#endif

    end subroutine gen_excit_rs_hubbard_transcorr