fix factorial2

docs/example/single_point/h2.py

@@ -15,7 +15,7 @@
 
 # define the wave function
 wf = SlaterJastrow(mol, kinetic='jacobi',
-                   configs='ground_state', jastrow=jastrow).gto2sto()
+                   configs='ground_state', jastrow=jastrow) #.gto2sto()
 
 # sampler
 sampler = Metropolis(nwalkers=1000, nstep=1000, step_size=0.25,

qmctorch/wavefunction/orbitals/norm_orbital.py

@@ -1,6 +1,6 @@
 import torch
 import numpy as np
-
+from scipy.special import factorial2
 
 def atomic_orbital_norm(basis):
     """Computes the norm of the atomic orbitals
@@ -141,16 +141,24 @@ def norm_gaussian_cartesian(a, b, c, exp):
         torch.tensor: normalization factor
     """
 
-    from scipy.special import factorial2 as f2
-
     pref = torch.as_tensor((2 * exp / np.pi) ** (0.75))
     am1 = (2 * a - 1).astype("int")
     x = (4 * exp) ** (a / 2) / torch.sqrt(torch.as_tensor(f2(am1)))
 
     bm1 = (2 * b - 1).astype("int")
     y = (4 * exp) ** (b / 2) / torch.sqrt(torch.as_tensor(f2(bm1)))
-
+    
     cm1 = (2 * c - 1).astype("int")
-    z = (4 * exp) ** (c / 2) / torch.sqrt(torch.as_tensor(f2(cm1)))
+    z = (4 * exp) ** (c / 2) / torch.sqrt(torch.as_tensor(f2(cm1))) 
 
     return (pref * x * y * z).type(torch.get_default_dtype())
+
+def f2(x):
+    """Returns the f2 of x with f2(x<1) = 1 as implemented in scipy 1.10.
+    """
+    # compute the x!!
+    out = factorial2(x)
+
+    # set all the elements lower than 1 to 1 
+    out[out<1] = 1
+    return out