precompute_pgen Subroutine

    subroutine precompute_pgen()

        ! set the number of unoccupied alpha/beta
        nUnoccAlpha = nBasis / 2 - nOccAlpha
        nUnoccBeta = nBasis / 2 - nOccBeta

        ! the number of valid triple excitations is just given by the binomial coefficients
        pgen3B = nOccBeta * (nOccBeta - 1) * (nOccBeta - 2)
        pgen3B = scaleInvert(same_spin_perm, pgen3B)

        pgen2B = nOccBeta * (nOccBeta - 1) * nOccAlpha
        pgen2B = scaleInvert(opp_spin_perm, pgen2B)

        pgen1B = nOccBeta * nOccAlpha * (nOccAlpha - 1)
        pgen1B = scaleInvert(opp_spin_perm, pgen1B)

        pgen0B = nOccAlpha * (nOccAlpha - 1) * (nOccAlpha - 2)
        pgen0B = scaleInvert(same_spin_perm, pgen0B)

    contains
        pure function scaleInvert(scl, p) result(sp)
            real(dp), intent(in) :: scl, p
            real(dp) :: sp

            if (p > eps) then
                sp = scl / p
            else
                sp = 0.0_dp
            end if
        end function scaleInvert
    end subroutine precompute_pgen