Codes using factorial2 function

[yingxingcheng]

Describe the bug
The factorial2 function in SciPy returns different values for negative inputs depending on the version. For versions older than 1.11.0, it returns 1, while for versions 1.11.0 and newer, it returns 0.

Code that utilizes the factorial2 function may encounter potential issues, such as those found in the contractions.py file:

    @property
    def norm_prim_cart(self):
        r"""Return the normalization constants of the Cartesian Gaussian primitives.

        For a Cartesian primitive with exponent :math:`\alpha_i`, the normalization constant is:

        .. math::
           N(\alpha_i, \vec{a}) = \sqrt {
           \left(\frac{2\alpha_i}{\pi}\right)^\frac{3}{2}
           \frac{(4\alpha_i)^{a_x + a_y + a_z}}{(2a_x - 1)!! (2a_y - 1)!! (2a_z - 1)!!}}

        Returns
        -------
        norm_prim_cart : np.ndarray(L, K)
            The normalization constants of the Cartesian Gaussian primitives.
            `L` is the number of contracted Cartesian Gaussian functions for the given angular
            momentum, i.e. :math:`(\ell + 1) * (\ell + 2) / 2`
            `K` is the number of exponents (i.e. primitives).

        """
        exponents = self.exps[np.newaxis, :]
        angmom_components_cart = self.angmom_components_cart[:, :, np.newaxis]

        return (
            (2 * exponents / np.pi) ** (3 / 4)
            * ((4 * exponents) ** (self.angmom / 2))
            / np.sqrt(np.prod(factorial2(2 * angmom_components_cart - 1), axis=1))
        )

An unwanted value will be returned when angmom_components_cart is 0, i.e., $s$ orbitals.

To Reproduce
Use SciPy with version $\geq 1.11.0$

[tovrstra]

This is related to the following:

• ENH: extensions of factorial{,2,k} for complex domains, recurrences, ... scipy/scipy#18409
• MBIS charges from fchk denspart#15

[PaulWAyers]

As a very tentative change, we can just specify the maximum scipy version as 1.10.1 or lower. That's a stopgap measure.

#131