diff --git a/setup.py b/setup.py
index c35dadb..d02b181 100644
--- 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
index 359d6f9..5228759 100644
--- 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
index c4cadb3..0999064 100644
--- 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
index 2f71236..c5f5a52 100644
--- 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():