calc_pgen_tJ_model Function

    function calc_pgen_tJ_model(ilutI, ex, ic) result(pgen)
        integer(n_int), intent(in) :: ilutI(0:NIfTot)
        integer, intent(in) :: ex(2, 2), ic
        real(dp) :: pgen
#ifdef DEBUG_
        character(*), parameter :: this_routine = "calc_pgen_tJ_model"
#endif

        integer :: src(2), tgt(2)
        real(dp) :: p_elec, p_orb, cum_sum, cpt_1, cpt_2
        real(dp), allocatable :: cum_arr(:)
        integer, allocatable :: tmp_list(:)

        ASSERT(ic >= 0)

        if (ic == 0 .or. ic > 2) then
            pgen = 0.0_dp
            return
        end if

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

        ASSERT(associated(lat))

        p_elec = 1.0_dp / real(nel, dp)

        if (ic == 1) then
            ! here it is easy..
            ASSERT(is_beta(src(1)) .eqv. is_beta(tgt(1)))
            ASSERT(any(tgt(1) == lat%get_spinorb_neighbors(src(1))))
            ASSERT(IsOcc(ilutI, src(1)))
            ASSERT(IsNotOcc(ilutI, tgt(1)))

            call create_cum_list_tJ_model(ilutI, src(1), lat%get_neighbors(gtid(src(1))), &
                                          cum_arr, cum_sum, tmp_list, tgt(1), p_orb)

        else if (ic == 2) then
            ! here we have to find the correct orbital..
            ! and i just realised that we maybe have to take into account
            ! of having picked the orbitals in a different order..
            ! for HPHF reasons i have to return here if the orbitals are
            ! not "correct"
            if (same_spin(src(1), src(2)) .or. same_spin(tgt(1), tgt(2))) then
                pgen = 0.0_dp
                return
            end if
            if (.not. (is_in_pair(src(1), tgt(1)) .or. is_in_pair(src(1), tgt(2)))) then
                pgen = 0.0_dp
                return
            end if
            if (.not. (is_in_pair(src(2), tgt(1)) .or. is_in_pair(src(2), tgt(2)))) then
                pgen = 0.0_dp
                return
            end if

            call create_cum_list_tJ_model(ilutI, src(1), lat%get_neighbors(gtid(src(1))), &
                                          cum_arr, cum_sum, tmp_list, src(2), cpt_1)
            call create_cum_list_tJ_model(ilutI, src(2), lat%get_neighbors(gtid(src(2))), &
                                          cum_arr, cum_sum, tmp_list, src(1), cpt_2)

            p_orb = cpt_1 + cpt_2

#ifdef DEBUG_
            if (is_beta(src(1)) .eqv. is_beta(tgt(1))) then
                ! then those to orbitls were chosen
                ASSERT(is_beta(src(2)) .eqv. is_beta(tgt(2)))
            else if (is_beta(src(1)) .eqv. is_beta(tgt(2))) then
                ASSERT(is_beta(src(2)) .eqv. is_beta(tgt(1)))
            end if

        else
            ! something went wrong
            call stop_all(this_routine, "something went wrong!")
#endif

        end if

        pgen = p_elec * p_orb

    end function calc_pgen_tJ_model