Add support for NumPy 2.1

requirements.txt

@@ -3,7 +3,7 @@ cython>=0.29.21,<3.0.10
 setuptools>=68.1.2
 packaging>=24.0
 # required scientific packages
-numpy>=1.26.4,<2.0
+numpy>=1.26.4,<2.2
 h5py>=3.10.0
 # optional scientific packages
 scipy>=1.11.4

samples/lb_planar_couette.py

@@ -92,10 +92,9 @@ def analytical(x, t, nu, v, h, k_max):
     velocity_ref = analytical(position_ref, time, nu, v, h, k_max)
     velocity_lbf = np.copy(lbf[5, :, 0].velocity[:, 0].reshape([-1]))
     ax = plt.gca()
-    color = next(ax._get_lines.prop_cycler)['color']
-    plt.plot(velocity_ref, position_ref, '-', color=color,
-             label=f"Analytical solution at t={time:.0f}")
-    plt.plot(velocity_lbf, position_lbf, 'o', color=color,
+    line, = plt.plot(velocity_ref, position_ref, '-',
+                     label=f"Analytical solution at t={time:.0f}")
+    plt.plot(velocity_lbf, position_lbf, 'o', color=line.get_color(),
              label=f"Simulated profile at t={time:.0f}")
 
 plt.xlabel('shear velocity')

src/python/espressomd/electrokinetics.py

@@ -447,6 +447,7 @@ def shape_squeeze(shape):
             return tuple(x for x in shape if x != 1)
 
         if shape_squeeze(values.shape) != shape_squeeze(target_shape):
+            target_shape = tuple([int(x) for x in target_shape])
             raise ValueError(
                 f"Input-dimensions of '{attr}' array {values.shape} does not match slice dimensions {target_shape}")
 

src/python/espressomd/lb.py

@@ -513,6 +513,7 @@ def shape_squeeze(shape):
             return tuple(x for x in shape if x != 1)
 
         if shape_squeeze(values.shape) != shape_squeeze(target_shape):
+            target_shape = tuple([int(x) for x in target_shape])
             raise ValueError(
                 f"Input-dimensions of '{attr}' array {values.shape} does not match slice dimensions {target_shape}")
 

testsuite/python/ek_bulk_reactions.py

@@ -112,15 +112,15 @@ def detail_test_reaction(self, single_precision: bool):
 
         self.system.integrator.run(self.TIME)
 
-        domain_volume = np.prod(ek_species_product.shape)
+        domain_volume = np.prod(np.copy(ek_species_product.shape))
         analytic_time = (self.TIME + 0.5) * self.system.time_step
 
         measured_educt_densities = np.zeros(len(stoech_coeffs))
         for i, educt in enumerate(educt_species):
-            measured_educt_densities[i] = np.sum(
-                educt[:, :, :].density) / domain_volume
-        measured_product_density = np.sum(
-            ek_species_product[:, :, :].density) / domain_volume
+            measured_educt_densities[i] = np.copy(np.sum(
+                educt[:, :, :].density)) / domain_volume
+        measured_product_density = np.copy(np.sum(
+            ek_species_product[:, :, :].density)) / domain_volume
 
         analytic_educt_densities = np.zeros(len(stoech_coeffs))
         for i, coeff in enumerate(stoech_coeffs):

testsuite/python/ek_eof.py

@@ -152,7 +152,7 @@ def calc_analytic_velocity(x, integration_constant, valency,
             x_sim < -self.system.box_l[0] / 2 + offset,
             x_sim > self.system.box_l[0] / 2 - offset)] = 0.
         np.testing.assert_allclose(
-            simulated_density, analytic_density, rtol=2e-2)
+            np.copy(simulated_density), analytic_density, rtol=2e-2)
 
         analytic_velocity = calc_analytic_velocity(
             x=x_sim,
@@ -168,7 +168,7 @@ def calc_analytic_velocity(x, integration_constant, valency,
             x_sim < -self.system.box_l[0] / 2 + offset,
             x_sim > self.system.box_l[0] / 2 - offset)] = 0.
         np.testing.assert_allclose(
-            simulated_velocity,
+            np.copy(simulated_velocity),
             analytic_velocity,
             rtol=2e-2)
 

testsuite/python/ek_interface.py

@@ -288,17 +288,18 @@ def test_grid_index(self):
             reactants=[ek_reactant], coefficient=1.5, lattice=self.lattice, tau=self.params["tau"])
         # check ranges and out-of-bounds access
         shape = np.around(self.system.box_l / self.params["agrid"]).astype(int)
+        int_shape = [int(x) for x in shape]  # cast away numpy integer types
         for i in range(3):
             n = [0, 0, 0]
             n[i] -= shape[i]
             ek_reaction[n[0], n[1], n[2]] = True
             self.assertTrue(ek_reaction[0, 0, 0])
             self.assertEqual(ek_reaction[tuple(n)], ek_reaction[0, 0, 0])
             self.assertEqual(ek_species[tuple(n)], ek_species[0, 0, 0])
-            for offset in (shape[i] + 1, -(shape[i] + 1)):
+            for offset in (int_shape[i] + 1, -(int_shape[i] + 1)):
                 n = [0, 0, 0]
                 n[i] += offset
-                err_msg = rf"provided index \[{str(n)[1:-1]}\] is out of range for shape \[{str(list(shape))[1:-1]}\]"  # nopep8
+                err_msg = rf"provided index \[{str(n)[1:-1]}\] is out of range for shape \[{str(int_shape)[1:-1]}\]"  # nopep8
                 with self.assertRaisesRegex(IndexError, err_msg):
                     ek_reaction[tuple(n)]
                 with self.assertRaisesRegex(IndexError, err_msg):

testsuite/python/lb.py

@@ -481,15 +481,16 @@ def test_grid_index(self):
         self.system.lb = lbf
         # check ranges and out-of-bounds access
         shape = lbf.shape
+        int_shape = [int(x) for x in shape]  # cast away numpy integer types
         for i in range(3):
             n = [0, 0, 0]
             n[i] -= shape[i]
             lbf[n[0], n[1], n[2]].velocity
             self.assertEqual(lbf[tuple(n)], lbf[0, 0, 0])
-            for offset in (shape[i] + 1, -(shape[i] + 1)):
+            for offset in (int_shape[i] + 1, -(int_shape[i] + 1)):
                 n = [0, 0, 0]
                 n[i] += offset
-                err_msg = rf"provided index \[{str(n)[1:-1]}\] is out of range for shape \[{str(list(shape))[1:-1]}\]"  # nopep8
+                err_msg = rf"provided index \[{str(n)[1:-1]}\] is out of range for shape \[{str(int_shape)[1:-1]}\]"  # nopep8
                 with self.assertRaisesRegex(IndexError, err_msg):
                     lbf[tuple(n)].velocity
         # node index
@@ -760,7 +761,7 @@ def test_ext_force_density(self):
         # Check node velocities
         for node_velocity in lbf[:, :, :].velocity.reshape((-1, 3)):
             np.testing.assert_allclose(
-                node_velocity, fluid_velocity, atol=1E-6)
+                np.copy(node_velocity), fluid_velocity, atol=1E-6)
 
     @utx.skipIfMissingFeatures("EXTERNAL_FORCES")
     def test_unequal_time_step(self):

testsuite/python/lb_lees_edwards_particle_coupling.py

@@ -292,10 +292,11 @@ def test_viscous_coupling_with_shear_vel(self):
         system.thermostat.set_lb(LB_fluid=lbf, seed=123, gamma=1)
         system.integrator.run(5000)
         for n in lbf[:, :, :]:
-            np.testing.assert_allclose(n.velocity[1:], [0, 0], atol=1E-8)
+            np.testing.assert_allclose(
+                np.copy(n.velocity[1:]), [0, 0], atol=1E-8)
         pos = np.random.random(3) * box_l
         p = system.part.add(pos=pos, v=lbf.get_interpolated_velocity(pos=pos))
-        np.testing.assert_allclose(p.v[1:], [0, 0], atol=1E-8)
+        np.testing.assert_allclose(np.copy(p.v)[1:], [0, 0], atol=1E-8)
         for _ in range(1000):
             system.integrator.run(1, reuse_forces=True)
             np.testing.assert_allclose(np.copy(p.f), np.zeros(3), atol=2E-6)

testsuite/python/lb_momentum_conservation.py

@@ -66,7 +66,7 @@ def test(self):
         self.system.lb = self.lbf
         self.system.thermostat.set_lb(LB_fluid=self.lbf, gamma=GAMMA, seed=1)
         np.testing.assert_allclose(
-            self.lbf.ext_force_density,
+            np.copy(self.lbf.ext_force_density),
             LB_PARAMS["ext_force_density"])
 
         # Initial momentum before integration = 0

testsuite/python/lees_edwards.py

@@ -416,15 +416,15 @@ def test_interactions(self):
             system.non_bonded_inter[0, 0].soft_sphere.set_params(
                 a=k_non_bonded / 2, n=-2, cutoff=r_cut)
             system.integrator.run(0)
-            r_12 = system.distance_vec(p1, p2)
+            r_12 = np.copy(system.distance_vec(p1, p2))
 
             np.testing.assert_allclose(
                 k_non_bonded * r_12, np.copy(p1.f))
             np.testing.assert_allclose(np.copy(p1.f), -np.copy(p2.f))
 
             np.testing.assert_allclose(
                 np.copy(system.analysis.pressure_tensor()["non_bonded"]),
-                np.outer(r_12, p2.f) / system.volume())
+                np.outer(r_12, np.copy(p2.f)) / system.volume())
 
             np.testing.assert_almost_equal(
                 system.analysis.energy()["non_bonded"],
@@ -456,13 +456,13 @@ def test_virt_sites(self):
         system.lees_edwards.set_boundary_conditions(
             shear_direction="x", shear_plane_normal="y", protocol=lin_protocol)
         # Test position and velocity of VS with Le shift
-        old_p3_pos = p3.pos
+        old_p3_pos = np.copy(p3.pos)
         expected_p3_pos = old_p3_pos - \
             np.array((get_lin_pos_offset(system.time, **params_lin), 0, 0))
         system.integrator.run(0, recalc_forces=True)
         np.testing.assert_allclose(np.copy(p3.pos_folded), expected_p3_pos)
         np.testing.assert_allclose(
-            p3.v, p1.v + np.array((params_lin["shear_velocity"], 0, 0)))
+            np.copy(p3.v), np.copy(p1.v) + np.array((params_lin["shear_velocity"], 0, 0)))
 
         # Check distances
         np.testing.assert_allclose(

testsuite/python/propagation_brownian.py

@@ -163,14 +163,14 @@ def test_fix_rotation(self):
         system.integrator.set_brownian_dynamics()
         system.integrator.run(3)
         np.testing.assert_allclose(
-            part.omega_lab, [0, 0, 1.3 / 1.5], atol=1e-14)
+            np.copy(part.omega_lab), [0, 0, 1.3 / 1.5], atol=1e-14)
 
         # noise only
         part.ext_torque = 3 * [0.]
         system.thermostat.set_brownian(
             kT=1, gamma=1, gamma_rotation=1.5, seed=41)
         system.integrator.run(3)
-        self.assertGreater(np.linalg.norm(part.omega_lab), 0.)
+        self.assertGreater(np.linalg.norm(np.copy(part.omega_lab)), 0.)
 
     @utx.skipIfMissingFeatures(["MASS",
                                 "ROTATIONAL_INERTIA",

testsuite/python/test_checkpoint.py

@@ -440,17 +440,17 @@ def test_particle_properties(self):
         if espressomd.has_features(['EXTERNAL_FORCES', 'ROTATION']):
             np.testing.assert_allclose(np.copy(p3.ext_torque), [0.3, 0.5, 0.7])
         if espressomd.has_features('ROTATIONAL_INERTIA'):
-            np.testing.assert_allclose(p3.rinertia, [2., 3., 4.])
+            np.testing.assert_allclose(np.copy(p3.rinertia), [2., 3., 4.])
         if espressomd.has_features('THERMOSTAT_PER_PARTICLE'):
             gamma = 2.
             if espressomd.has_features('PARTICLE_ANISOTROPY'):
                 if 'THERM.LB' in modes:
                     gamma = np.array([2., 2., 2.])
                 else:
                     gamma = np.array([2., 3., 4.])
-            np.testing.assert_allclose(p4.gamma, gamma)
+            np.testing.assert_allclose(np.copy(p4.gamma), gamma)
             if espressomd.has_features('ROTATION'):
-                np.testing.assert_allclose(p3.gamma_rot, 2. * gamma)
+                np.testing.assert_allclose(np.copy(p3.gamma_rot), 2. * gamma)
         if espressomd.has_features('ENGINE'):
             self.assertEqual(
                 p3.swimming,
@@ -466,11 +466,13 @@ def test_particle_properties(self):
             from scipy.spatial.transform import Rotation as R
             q_ind = ([1, 2, 3, 0],)  # convert from scalar-first to scalar-last
             vs_id, vs_dist, vs_quat = p2.vs_relative
-            d = p2.pos - p1.pos
+            d = np.copy(p2.pos - p1.pos)
+            vs_quat = np.copy(vs_quat)
+            p_quat = np.copy(p1.quat)
             theta = np.arccos(d[2] / np.linalg.norm(d))
             assert abs(theta - 3. * np.pi / 4.) < 1e-8
             q = np.array([0., 0., np.sin(theta / 2.), -np.cos(theta / 2.)])
-            r = R.from_quat(p1.quat[q_ind]) * R.from_quat(vs_quat[q_ind])
+            r = R.from_quat(p_quat[q_ind]) * R.from_quat(vs_quat[q_ind])
             self.assertEqual(vs_id, p1.id)
             np.testing.assert_allclose(vs_dist, np.sqrt(2.))
             np.testing.assert_allclose(q[q_ind], r.as_quat(), atol=1e-10)
@@ -761,7 +763,7 @@ def test_virtual_sites(self):
         self.assertEqual(p_real.vs_relative[1], 0.)
         self.assertEqual(p_virt.vs_relative[1], np.sqrt(2.))
         np.testing.assert_allclose(
-            p_real.vs_relative[2], [1., 0., 0., 0.], atol=1e-10)
+            np.copy(p_real.vs_relative[2]), [1., 0., 0., 0.], atol=1e-10)
 
     def test_mean_variance_calculator(self):
         acc_mean_variance = system.auto_update_accumulators[0]
@@ -1034,8 +1036,8 @@ def test_union(self):
         p1 = system.part.add(pos=[1., 1.6, 0.], type=6)
         p2 = system.part.add(pos=[system.box_l[0] - 1., 1.6, 0.], type=6)
         system.integrator.run(0, recalc_forces=True)
-        np.testing.assert_allclose(p1.f, [0., 1e8, 0.], atol=1e-3)
-        np.testing.assert_allclose(p2.f, [0., 1e8, 0.], atol=1e-3)
+        np.testing.assert_allclose(np.copy(p1.f), [0., 1e8, 0.], atol=1e-3)
+        np.testing.assert_allclose(np.copy(p2.f), [0., 1e8, 0.], atol=1e-3)
         p1.remove()
         p2.remove()
         system.non_bonded_inter[2, 6].reset()