gen_excit_rs_hubbard_transcorr_uniform Subroutine

    subroutine gen_excit_rs_hubbard_transcorr_uniform(nI, ilutI, nJ, ilutJ, exFlag, ic, &
                                                      ex, tParity, pGen, hel, store, run)
        ! also create an uniform excitation generator for the hopping
        ! transcorrelated hubbard. mainly to test where the instabilities
        ! in the weighted come from

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

        integer :: elecs(2), orbs(2), src(2), spin
        real(dp) :: 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 pick_spin_opp_elecs(nI, elecs, p_elec)

            src = nI(elecs)

            ASSERT(.not. same_spin(src(1), src(2)))

            call pick_spin_opp_holes(ilutI, orbs, p_orb)

            ASSERT(.not. same_spin(orbs(1), orbs(2)))

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

            ! do i need a factor of 2? since orbitals could be switched
            ! the other way around?
            pgen = p_elec * p_orb * pDoubles

            call make_double(nI, nJ, elecs(1), elecs(2), orbs(1), orbs(2), ex, tParity)
            ilutJ = make_ilutJ(ilutI, ex, 2)

            ! to be compatible with my test-suite actually calculate the
            ! matrix element here..
            if (abs(get_double_helem_rs_hub_transcorr(ex, .false.)) < EPS) then
                nJ(1) = 0
                pgen = 0.0_dp
                return
            end if
        else

            ic = 1

            elecs(1) = 1 + int(genrand_real2_dsfmt() * nel)

            src(1) = nI(elecs(1))

            p_elec = 1.0_dp / real(nel, dp)

            spin = get_spin(src(1)) - 1
            ! and now pick a spin-parallel hole!
            call pick_random_hole(ilutI, orbs(1), p_orb, spin)

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

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

            call make_single(nI, nJ, elecs(1), orbs(1), ex, tParity)
            ilutJ = make_ilutJ(ilutI, ex, 1)

            if (abs(get_single_helem_rs_hub_transcorr(nI, ex, .false.)) < EPS) then
                nJ(1) = 0
                pgen = 0.0_dp
                return
            end if
        end if

#ifdef DEBUG_
        temp_pgen = calc_pgen_rs_hubbard_transcorr_uniform(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_uniform