Make qchem compatible with scipy factorial2

doc/releases/changelog-dev.md

@@ -30,6 +30,10 @@
 * QNode transforms in `qml.qinfo` now support custom wire labels.
   [#4331](https://github.com/PennyLaneAI/pennylane/pull/4331)
 
+* The `qchem` functions `primitive_norm` and `contracted_norm` are modified to be compatible with
+  higher versions of scipy. The private function `_fac2` for computing double factorials is added. 
+  [#4321](https://github.com/PennyLaneAI/pennylane/pull/4321)
+
 * The default label for a `StatePrep` operator is now `|Ψ⟩`.
   [(#4340)](https://github.com/PennyLaneAI/pennylane/pull/4340)
 
@@ -111,6 +115,7 @@
 
 This release contains contributions from (in alphabetical order):
 
+Soran Jahangiri,
 Isaac De Vlugt,
 Edward Jiang,
 Christina Lee,

pennylane/qchem/integrals.py

@@ -56,10 +56,11 @@ def primitive_norm(l, alpha):
     array([1.79444183])
     """
     lx, ly, lz = l
+
     n = (
         (2 * alpha / np.pi) ** 0.75
         * (4 * alpha) ** (sum(l) / 2)
-        / qml.math.sqrt(fac2(2 * lx - 1) * fac2(2 * ly - 1) * fac2(2 * lz - 1))
+        / qml.math.sqrt(_fac2(2 * lx - 1) * _fac2(2 * ly - 1) * _fac2(2 * lz - 1))
     )
     return n
 
@@ -102,14 +103,30 @@ def contracted_norm(l, alpha, a):
     0.39969026908800853
     """
     lx, ly, lz = l
-    c = np.pi**1.5 / 2 ** sum(l) * fac2(2 * lx - 1) * fac2(2 * ly - 1) * fac2(2 * lz - 1)
+
+    c = np.pi**1.5 / 2 ** sum(l) * _fac2(2 * lx - 1) * _fac2(2 * ly - 1) * _fac2(2 * lz - 1)
     s = (
         (a.reshape(len(a), 1) * a) / ((alpha.reshape(len(alpha), 1) + alpha) ** (sum(l) + 1.5))
     ).sum()
     n = 1 / qml.math.sqrt(c * s)
     return n
 
 
+def _fac2(n):
+    """Compute the double factorial of an integer.
+
+    The function uses the definition :math:`(-1)!! = 1`.
+
+    Args:
+        n (int): number for which the double factorial is computed
+
+    Returns:
+        int: the computed double factorial
+
+    """
+    return int(fac2(n) if n != -1 else 1)
+
+
 def _generate_params(params, args):
     """Generate basis set parameters. The default values are used for the non-differentiable
     parameters and the user-defined values are used for the differentiable ones.

tests/qchem/test_integrals.py

@@ -152,6 +152,23 @@ def test_hermite_coulomb(self, t, u, v, n, p, dr, h_ref):
         h = qchem.integrals._hermite_coulomb(t, u, v, n, p, dr)
         assert np.allclose(h, h_ref)
 
+    @pytest.mark.parametrize(
+        ("n", "result"),
+        [
+            (0, 1),
+            (-1, 1),
+            (-2, 0),
+            (-3, 0),
+            (2, 2),
+            (5, 15),
+            (8, 384),
+        ],
+    )
+    def test_fac2(self, n, result):
+        r"""Test that the _fac2 function returns a correct value."""
+        value = qchem.integrals._fac2(n)
+        assert value == result
+
 
 class TestOverlap:
     """Tests for overlap integrals"""