gen_excit_heisenberg_model Subroutine

    subroutine gen_excit_heisenberg_model(nI, ilutI, nJ, ilutJ, exFlag, ic, &
                                          ex, tParity, pGen, hel, store, run)
        ! the heisenberg excitation generator is only a small modification of
        ! the t-J excitation generator without single excitation hoppings,
        ! due to half-filling

        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
#ifdef DEBUG_
        character(*), parameter :: this_routine = "gen_excit_heisenberg_model"
#endif
        integer :: elec, src, ind, tgt, src_opp, tgt_opp, elec_opp
        real(dp) :: p_elec, cum_sum, r, p_orb, cpt_opp
        integer, allocatable :: neighbors(:)
        real(dp), allocatable :: cum_arr(:)

        unused_var(exFlag)
        unused_var(store)
        if (present(run)) then
            unused_var(run)
        end if

        hel = h_cast(0.0_dp)

        ASSERT(associated(lat))
        ASSERT(nel == nbasis / 2)

        ic = 2

        ! still pick the first electron at random
        elec = 1 + int(genrand_real2_dsfmt() * nel)

        p_elec = 1.0_dp / real(nel, dp)

        src = nI(elec)

        ! in the heisenberg model i am sure that every site is occupied..
        ! so i could pick the spin-orbit neighbors and check if the spin-orbital
        ! is empty, this would indicate opposite neighboring spins..
        neighbors = lat%get_spinorb_neighbors(src)

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

        if (cum_sum < EPS) then
            nJ(1) = 0
            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

        ! and then i have to check for the opposite order generation prob
        if (is_beta(src)) then
            src_opp = neighbors(ind) + 1
            tgt_opp = src + 1
        else
            src_opp = neighbors(ind) - 1
            tgt_opp = src - 1
        end if

        ASSERT(IsOcc(ilutI, src_opp))
        ASSERT(IsNotOcc(ilutI, tgt_opp))

        call create_cum_list_heisenberg(ilutI, src_opp, lat%get_spinorb_neighbors(src_opp), &
                                        cum_arr, cum_sum, tgt_opp, cpt_opp)

        pgen = p_elec * (p_orb + cpt_opp)

        elec_opp = find_elec_in_ni(nI, src_opp)

        call make_double(nI, nJ, elec, elec_opp, tgt, tgt_opp, ex, tParity)

        ilutJ = make_ilutJ(ilutI, ex, 2)

    end subroutine gen_excit_heisenberg_model