pure subroutine calc_orbital_pgen_contr_pchb(this, csf_i, occ_orbs, cpt_a, cpt_b)
class(GugaAliasSampler_t), intent(in) :: this
type(CSF_Info_t), intent(in) :: csf_i
integer, intent(in) :: occ_orbs(2)
real(dp), intent(out) :: cpt_a, cpt_b
integer :: ij
unused_var(csf_i)
! this function can in theory be called with both i < j and i > j..
! to take the correct values here!
! although for the fuseIndex function the order is irrelevant
! i think i have to consider both the above and below contribution..
! but i am not so sure how.. in this PCHB case..
! maybe in PCHB those 2 are just the same.. i am confused
ij = fuseIndex(gtID(occ_orbs(1)), gtID(occ_orbs(2)))
! and both I and J are electron and hole indices here
cpt_a = this%alias_sampler%get_prob(ij, ij) / 2.0_dp
! but since they will get added in later routines i actually have
! do divide by 2 here..
! and in the PCHB there is no difference between the 2!
cpt_b = cpt_a
end subroutine calc_orbital_pgen_contr_pchb