Support recent CMake and NumPy versions (#4992)

Description of changes:
- add experimental support for NumPy 2
- modernize CMake build system

CMakeLists.txt

@@ -20,18 +20,19 @@
 #
 
 cmake_minimum_required(VERSION 3.25.1)
-message(STATUS "CMake version: ${CMAKE_VERSION}")
-if(POLICY CMP0076)
-  # make target_sources() convert relative paths to absolute
-  cmake_policy(SET CMP0076 NEW)
-endif()
-if(POLICY CMP0025)
-  # make CXX_COMPILER_ID return "AppleClang" instead of "Clang" for Apple Clang
-  cmake_policy(SET CMP0025 NEW)
+cmake_policy(VERSION 3.25.1)
+if(POLICY CMP0167)
+  # use BoostConfig.cmake shipped with Boost 1.70+ instead of the one in CMake
+  cmake_policy(SET CMP0167 NEW)
 endif()
+message(STATUS "CMake version: ${CMAKE_VERSION}")
 # CMake modules/macros are in a subdirectory to keep this file cleaner
 set(CMAKE_MODULE_PATH ${CMAKE_CURRENT_SOURCE_DIR}/cmake)
 
+# project info
+project(ESPResSo)
+set(PROJECT_VERSION "4.3-dev")
+
 # C++ standard
 enable_language(CXX)
 set(CMAKE_CXX_STANDARD 20 CACHE STRING "C++ standard to be used")
@@ -54,16 +55,10 @@ espresso_minimal_compiler_version("Intel" 2021.9)
 espresso_minimal_compiler_version("IntelLLVM" 2023.1)
 
 include(FeatureSummary)
-project(ESPResSo)
 include(GNUInstallDirs)
+include(FetchContent)
 include(espresso_option_enum)
 include(espresso_enable_avx2_support)
-if(POLICY CMP0074)
-  # make find_package() use <PackageName>_ROOT variables
-  cmake_policy(SET CMP0074 NEW)
-endif()
-
-set(PROJECT_VERSION "4.3-dev")
 
 #
 # CMake internal vars
@@ -80,6 +75,31 @@ set(CMAKE_CXX_FLAGS_RELWITHASSERT "${CMAKE_CXX_FLAGS_RELWITHASSERT} -O3 -g")
 # On Mac OS X, first look for other packages, then frameworks
 set(CMAKE_FIND_FRAMEWORK LAST)
 
+# avoid applying patches twice with FetchContent
+set(FETCHCONTENT_UPDATES_DISCONNECTED ON)
+
+#
+# External dependencies declarations
+#
+
+# cmake-format: off
+FetchContent_Declare(
+  walberla
+  GIT_REPOSITORY https://i10git.cs.fau.de/walberla/walberla.git
+  GIT_TAG        b0842e1a493ce19ef1bbb8d2cf382fc343970a7f
+)
+FetchContent_Declare(
+  stokesian_dynamics
+  GIT_REPOSITORY https://github.com/hmenke/espresso-stokesian-dynamics.git
+  GIT_TAG        862a7537a366f0c32f0c25e46bd107bea590faea
+)
+FetchContent_Declare(
+  caliper
+  GIT_REPOSITORY https://github.com/LLNL/Caliper.git
+  GIT_TAG        v2.10.0
+)
+# cmake-format: on
+
 # ##############################################################################
 # User input options
 # ##############################################################################
@@ -316,17 +336,10 @@ endif()
 
 if(ESPRESSO_BUILD_WITH_STOKESIAN_DYNAMICS)
   set(CMAKE_INSTALL_LIBDIR "${ESPRESSO_INSTALL_LIBDIR}")
-  include(FetchContent)
-  FetchContent_Declare(
-    stokesian_dynamics
-    GIT_REPOSITORY https://github.com/hmenke/espresso-stokesian-dynamics.git
-    GIT_TAG 862a7537a366f0c32f0c25e46bd107bea590faea)
   FetchContent_GetProperties(stokesian_dynamics)
   set(STOKESIAN_DYNAMICS 1)
   if(NOT stokesian_dynamics_POPULATED)
-    FetchContent_Populate(stokesian_dynamics)
-    add_subdirectory(${stokesian_dynamics_SOURCE_DIR}
-                     ${stokesian_dynamics_BINARY_DIR})
+    FetchContent_MakeAvailable(stokesian_dynamics)
   endif()
 endif()
 
@@ -600,18 +613,7 @@ endif()
 #
 
 if(ESPRESSO_BUILD_WITH_WALBERLA)
-  # cmake-format: off
-  include(FetchContent)
-  FetchContent_Declare(
-    walberla
-    GIT_REPOSITORY https://i10git.cs.fau.de/walberla/walberla.git
-    GIT_TAG        b0842e1a493ce19ef1bbb8d2cf382fc343970a7f
-  )
-  # workaround for https://gitlab.kitware.com/cmake/cmake/-/issues/21146
-  if(NOT DEFINED walberla_SOURCE_DIR OR NOT EXISTS "${walberla_SOURCE_DIR}")
-    FetchContent_Populate(walberla)
-  endif()
-  # cmake-format: on
+  FetchContent_GetProperties(walberla)
   string(REGEX REPLACE "([/\\]walberla)-src$" "\\1-build" walberla_BINARY_DIR
                        "${walberla_SOURCE_DIR}")
   set(WALBERLA_BUILD_TESTS off CACHE BOOL "")
@@ -637,7 +639,9 @@ if(ESPRESSO_BUILD_WITH_WALBERLA)
     set(WALBERLA_BUILD_WITH_FFTW off CACHE BOOL "")
   endif()
   set(WALBERLA_BUILD_WITH_FASTMATH off CACHE BOOL "")
-  add_subdirectory("${walberla_SOURCE_DIR}" "${walberla_BINARY_DIR}")
+  if(NOT walberla_POPULATED)
+    FetchContent_MakeAvailable(walberla)
+  endif()
   set(WALBERLA_LIBS
       walberla::core walberla::domain_decomposition walberla::blockforest
       walberla::boundary walberla::field walberla::lbm walberla::timeloop
@@ -664,23 +668,15 @@ if(ESPRESSO_BUILD_WITH_WALBERLA)
 endif()
 
 if(ESPRESSO_BUILD_WITH_CALIPER)
-  # cmake-format: off
-  include(FetchContent)
-  FetchContent_Declare(
-    caliper
-    GIT_REPOSITORY https://github.com/LLNL/Caliper.git
-    GIT_TAG        v2.10.0
-  )
-  if(NOT DEFINED caliper_SOURCE_DIR OR NOT EXISTS "${caliper_SOURCE_DIR}")
-    FetchContent_Populate(caliper)
-  endif()
-  # cmake-format: on
+  FetchContent_GetProperties(caliper)
   set(CALIPER_OPTION_PREFIX on CACHE BOOL "")
   set(CALIPER_WITH_MPI on CACHE BOOL "")
   set(CALIPER_WITH_NVTX off CACHE BOOL "")
   set(CALIPER_WITH_CUPTI off CACHE BOOL "")
   set(CALIPER_BUILD_SHARED_LIBS on CACHE BOOL "")
-  add_subdirectory("${caliper_SOURCE_DIR}")
+  if(NOT caliper_POPULATED)
+    FetchContent_MakeAvailable(caliper)
+  endif()
   target_compile_options(
     caliper-services
     PRIVATE

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",