# calc_pgen_k_space_hubbard_triples Function

## public function calc_pgen_k_space_hubbard_triples(nI, ilutI, ex, ic) result(pgen)

### Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: nI(nel)
integer(kind=n_int), intent(in) :: ilutI(0:niftot)
integer, intent(in) :: ex(:,:)
integer, intent(in) :: ic

## Source Code

    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