calc_pgen_k_space_hubbard_triples Function

    function calc_pgen_k_space_hubbard_triples(nI, ilutI, ex, ic) result(pgen)
        integer, intent(in) :: nI(nel), ex(:, :), ic
        integer(n_int), intent(in) :: ilutI(0:niftot)
        real(dp) :: pgen
#ifdef DEBUG_
        real(dp) :: test
#endif
        real(dp) :: p_elec, p_orb(2), cum_arr(nbasis / 2), cum_sum
        integer :: orb_list(nbasis / 2, 2), sum_ms, orb_a, orbs(2)

        if (ic /= 3) then
            pgen = 0.0_dp
            return
        end if

        sum_ms = sum(get_spin_pn(ex(1, :)))

        ! check spins
        if (.not. (sum_ms == 1 .or. sum_ms == -1) .or. sum_ms /= sum(get_spin_pn(ex(2, :)))) then
            pgen = 0.0_dp
            return
        end if

        ! get the probabilites for the electrons and orbital (a)
        if (sum_ms == 1) then
            p_elec = 1.0_dp / real(nel * (nel - 1), dp) * &
                     (1.0_dp / real(nOccBeta, dp) + 2.0_dp / real(nel - 2, dp))

            p_orb(1) = 1.0_dp / real(nbasis / 2 - nOccBeta, dp)

        else
            p_elec = 1.0_dp / real(nel * (nel - 1), dp) * &
                     (1.0_dp / real(nOccAlpha, dp) + 2.0_dp / real(nel - 2, dp))

            p_orb(1) = 1.0_dp / real(nbasis / 2 - nOccAlpha, dp)

        end if

        ! for this i need the minority spin orbital (a)
        orb_a = find_minority_spin(ex(2, :))

        orbs = pack(ex(2, :), ex(2, :) /= orb_a)

        call create_bc_list_hubbard(nI, ilutI, ex(1, :), orb_a, orb_list, cum_arr, &
                                    cum_sum, orbs(1), p_orb(2))

        pgen = p_elec * product(p_orb) * 2.0_dp

#ifdef DEBUG_
        call create_bc_list_hubbard(nI, ilutI, ex(1, :), orb_a, orb_list, cum_arr, &
                                    cum_sum, orbs(2), test)

        if (abs(test - p_orb(2)) > 1.e-8) then
            print *, "pgen assumption wrong:!"
            print *, "p_orb: ", p_orb(2)
            print *, "test: ", test
            print *, "ex(2,:): ", ex(2, :)
        end if
#endif

    end function calc_pgen_k_space_hubbard_triples