Support numpy >= 2.4 and relax the numpy requirement to <3

Two changes are needed to work with numpy 2.4, both backward compatible down to numpy 1.x: 1. PyACECalculator/PyACEEnsembleCalculator filled the ASE results dict via np.float64(energy.reshape(-1,)). numpy 2.4 expired the 'conversion of ndim > 0 arrays to scalars' deprecation (present since numpy 1.25), so the scalar constructor now returns a shape-(1,) array instead of collapsing it to a scalar, and atoms.get_potential_energy() leaked a 1-element array of the total energy. Scalar results now go through float(), which behaves identically on every numpy version, and per-atom arrays through astype(np.float64), which also keeps energies_dev/forces_dev arrays for single-atom structures on older numpy. 2. pyace.radial used np.trapz, which numpy 2.0 kept only as a deprecated alias of np.trapezoid and numpy 2.4 removed. The module now picks whichever of the two exists. A regression assertion in tests/test_PyACECalculator.py guards the scalar energy output. Verified by building the package and running tests/test_PyACECalculator.py against numpy 1.26.4, 2.3.5 and 2.4.6. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_015Ui6Drmop9gZXo1fgn8yVg
srmnitc · Sep 4, 2026 · Jul 16, 2026 · Sep 4, 2026 · Jul 16, 2026 · 35b2620e7864feb4ffc2c9a725bc8b2a39bade7a
commit 35b2620e7864feb4ffc2c9a725bc8b2a39bade7a
diff --git a/setup.py b/setup.py
@@ -192,7 +192,7 @@ def build_extension(self, ext: CMakeExtension) -> None:
                                           build_ext=CMakeBuild)),
     zip_safe=False,
     url='https://github.com/ICAMS/python-ace',
-    install_requires=['numpy<=1.26.4',
+    install_requires=['numpy<3',
                       'ase',
                       'pandas<=2.0',
                       'ruamel.yaml',

diff --git a/src/pyace/asecalc.py b/src/pyace/asecalc.py
@@ -164,8 +164,8 @@ def calculate(self, atoms=None, properties=('energy', 'forces', 'stress', 'energ
         self.energies = np.array(self.ace.energies)
 
         self.results = {
-            'energy': np.float64(self.energy.reshape(-1, )),
-            'free_energy': np.float64(self.energy.reshape(-1, )),
+            'energy': float(self.energy),
+            'free_energy': float(self.energy),
             'forces': self.forces.astype(np.float64),
             'energies': self.energies.astype(np.float64),
             'gamma': np.array(self.ace.gamma_grade, dtype=np.float64)
@@ -305,21 +305,21 @@ def calculate(self, atoms=None, properties=('energy', 'forces', 'stress', 'energ
 
         self.results = {
             # mean
-            'energy': np.float64(self.energy.reshape(-1, )),
-            'free_energy': np.float64(self.energy.reshape(-1, )),
+            'energy': float(self.energy),
+            'free_energy': float(self.energy),
             'forces': self.forces.astype(np.float64),
             'energies': self.energies.astype(np.float64),
 
             # std
-            'energy_std': np.float64(self.energy_std.reshape(-1, )),
-            'free_energy_std': np.float64(self.energy_std.reshape(-1, )),
+            'energy_std': float(self.energy_std),
+            'free_energy_std': float(self.energy_std),
             'forces_std': self.forces_std.astype(np.float64),
             'energies_std': self.energies_std.astype(np.float64),
 
             # dev
-            'energy_dev': np.float64(self.energy_dev),
-            'energies_dev': np.float64(self.energies_dev),
-            'forces_dev': np.float64(self.forces_dev)
+            'energy_dev': float(self.energy_dev),
+            'energies_dev': self.energies_dev.astype(np.float64),
+            'forces_dev': self.forces_dev.astype(np.float64)
         }
 
         if self.atoms.number_of_lattice_vectors == 3:

diff --git a/src/pyace/radial.py b/src/pyace/radial.py
@@ -6,6 +6,9 @@
 import numpy as np
 from typing import Union
 
+# np.trapz was renamed to np.trapezoid in numpy 2.0 and removed in numpy 2.4
+_trapezoid = np.trapezoid if hasattr(np, "trapezoid") else np.trapz
+
 from pyace import ACECTildeBasisSet, ACEBBasisSet, BBasisConfiguration
 
 
@@ -17,7 +20,7 @@ def integrate(xs, table):
         frs = np.abs(table)
         sum_frs = np.sum(frs, axis=(1, 2))
         integrand = sum_frs * xs ** 2
-        integral = np.trapz(integrand, x=xs)
+        integral = _trapezoid(integrand, x=xs)
         return integral
     else:
         return 0

diff --git a/tests/test_PyACECalculator.py b/tests/test_PyACECalculator.py
@@ -48,6 +48,10 @@ def test_setup():
     f1 = a.get_forces()
     print(e1)
     print(f1)
+    # regression check for numpy >= 2.4, where the energy leaked out as a
+    # 1-element array instead of a scalar
+    assert isinstance(e1, float)
+    assert f1.shape == (2, 3)
 
 
 def test_load_YAML():