calc_pgen_heisenberg_model Function

    function calc_pgen_heisenberg_model(ilutI, ex, ic) result(pgen)
        integer(n_int), intent(in) :: ilutI(0:NIfTot)
        integer, intent(in) :: ex(2, ic), ic
        real(dp) :: pgen
#ifdef DEBUG_
        character(*), parameter :: this_routine = "calc_pgen_heisenberg_model"
#endif
        integer :: src(2), tgt(2)
        real(dp) :: p_elec, cum_sum, cpt_1, cpt_2
        real(dp), allocatable :: cum_arr(:)

        ASSERT(ic >= 0)
        ASSERT(associated(lat))
        if (ic /= 2) then
            pgen = 0.0_dp
            return
        end if

        src = get_src(ex)
        tgt = get_tgt(ex)

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

        p_elec = 1.0_dp / real(nel, dp)

        if (is_beta(src(1)) .eqv. is_beta(tgt(1))) then
            ASSERT(is_in_pair(src(1), tgt(2)))
            ASSERT(is_in_pair(src(2), tgt(1)))
            ASSERT(same_spin(src(2), tgt(2)))

            call create_cum_list_heisenberg(ilutI, src(1), lat%get_spinorb_neighbors(src(1)), &
                                            cum_arr, cum_sum, tgt(1), cpt_1)
            call create_cum_list_heisenberg(ilutI, src(2), lat%get_spinorb_neighbors(src(2)), &
                                            cum_arr, cum_sum, tgt(2), cpt_2)

        else if (is_beta(src(1)) .eqv. is_beta(tgt(2))) then
            ASSERT(is_in_pair(src(1), tgt(1)))
            ASSERT(is_in_pair(src(2), tgt(2)))
            ASSERT(same_spin(src(2), tgt(1)))
            call create_cum_list_heisenberg(ilutI, src(1), lat%get_spinorb_neighbors(src(1)), &
                                            cum_arr, cum_sum, tgt(2), cpt_1)
            call create_cum_list_heisenberg(ilutI, src(2), lat%get_spinorb_neighbors(src(2)), &
                                            cum_arr, cum_sum, tgt(1), cpt_2)
#ifdef DEBUG_
        else
            call stop_all(this_routine, "something went wrong!")
#endif
        end if

        pgen = p_elec * (cpt_1 + cpt_2)

    end function calc_pgen_heisenberg_model