function calc_pgen_k_space_hubbard_uniform_transcorr(ex, ic) result(pgen)
! need a calc pgen functionality for the uniform transcorrelated
! excitation generator
integer, intent(in) :: ex(:, :), ic
real(dp) :: pgen
#ifdef DEBUG_
character(*), parameter :: this_routine = "calc_pgen_k_space_hubbard_uniform_transcorr"
#endif
real(dp) :: p_elec, p_orb
integer :: sum_ms
pgen = 0.0_dp
if (ic == 2) then
pgen = pDoubles
if (same_spin(ex(1, 1), ex(1, 2))) then
pgen = pgen * pParallel
if (is_beta(ex(1, 1))) then
p_elec = 1.0_dp / real(nOccBeta * (nOccBeta - 1), dp)
p_orb = 2.0_dp / real(nbasis / 2 - nOccBeta, dp)
else
p_elec = 1.0_dp / real(nOccAlpha * (nOccAlpha - 1), dp)
p_orb = 2.0_dp / real(nbasis / 2 - nOccAlpha, dp)
end if
else
pgen = pgen * (1.0_dp - pParallel)
p_elec = 1.0_dp / real(nOccBeta * nOccAlpha, dp)
p_orb = 2.0_dp / real(nbasis - nel, dp)
end if
else
pgen = 1.0_dp - pDoubles
sum_ms = sum(get_spin_pn(ex(1, :)))
ASSERT(sum_ms == 1 .or. sum_ms == -1)
if (sum_ms == 1) then
p_elec = 2.0_dp / real(nel * (nel - 1), dp) * &
(1.0_dp / real(nOccBeta, dp) + 2.0_dp / real(nel - 2, dp))
p_orb = 1.0_dp / real(nbasis / 2 - nOccBeta, dp) * &
2.0_dp / real(nbasis / 2 - nOccAlpha, dp)
else if (sum_ms == -1) then
p_elec = 2.0_dp / real(nel * (nel - 1), dp) * &
(1.0_dp / real(nOccAlpha, dp) + 2.0_dp / real(nel - 2, dp))
p_orb = 1.0_dp / real(nbasis / 2 - nOccAlpha, dp) * &
2.0_dp / real(nBasis / 2 - nOccBeta, dp)
end if
end if
pgen = pgen * p_elec * p_orb
end function calc_pgen_k_space_hubbard_uniform_transcorr