Refactor C++ code

Replace Utils functionality by equivalent C++20 features.
Remove unused or duplicated code. Improve code coverage.

src/core/MpiCallbacks.hpp

@@ -36,8 +36,6 @@
  */
 
 #include <utils/NumeratedContainer.hpp>
-#include <utils/tuple.hpp>
-#include <utils/type_traits.hpp>
 
 #include <boost/mpi/collectives/broadcast.hpp>
 #include <boost/mpi/communicator.hpp>
@@ -70,10 +68,6 @@ using is_allowed_argument =
                                (!std::is_const_v<std::remove_reference_t<T>> &&
                                 std::is_lvalue_reference_v<T>))>;
 
-template <class... Args>
-using are_allowed_arguments =
-    typename Utils::conjunction<is_allowed_argument<Args>...>::type;
-
 /**
  * @brief Invoke a callable with arguments from an mpi buffer.
  *
@@ -87,14 +81,14 @@ using are_allowed_arguments =
  */
 template <class F, class... Args>
 auto invoke(F f, boost::mpi::packed_iarchive &ia) {
-  static_assert(are_allowed_arguments<Args...>::value,
+  static_assert(std::conjunction_v<is_allowed_argument<Args>...>,
                 "Pointers and non-const references are not allowed as "
                 "arguments for callbacks.");
 
   /* This is the local receive buffer for the parameters. We have to strip
      away const so we can actually deserialize into it. */
   std::tuple<std::remove_const_t<std::remove_reference_t<Args>>...> params;
-  Utils::for_each([&ia](auto &e) { ia >> e; }, params);
+  std::apply([&ia](auto &&...e) { ((ia >> e), ...); }, params);
 
   /* We add const here, so that parameters can only be by value
      or const reference. Output parameters on callbacks are not
@@ -361,8 +355,8 @@ class MpiCallbacks {
     oa << id;
 
     /* Pack the arguments into a packed mpi buffer. */
-    Utils::for_each([&oa](auto &&e) { oa << e; },
-                    std::forward_as_tuple(std::forward<Args>(args)...));
+    std::apply([&oa](auto &&...e) { ((oa << e), ...); },
+               std::forward_as_tuple(std::forward<Args>(args)...));
 
     boost::mpi::broadcast(m_comm, oa, 0);
   }

src/core/TabulatedPotential.hpp

@@ -38,7 +38,7 @@ struct TabulatedPotential {
   double minval = -1.0;
   /** Position on the x-axis of the last tabulated value. */
   double maxval = -1.0;
-  /** Distance on the x-axis between tabulated values. */
+  /** Inverse distance on the x-axis between tabulated values. */
   double invstepsize = 0.0;
   /** Tabulated forces. */
   std::vector<double> force_tab;

src/core/bonded_interactions/angle_common.hpp

@@ -29,19 +29,10 @@
 
 #include <utils/Vector.hpp>
 
+#include <algorithm>
 #include <cmath>
 #include <tuple>
 
-namespace detail {
-inline double sanitize_cosine(double cosine) {
-  if (cosine > TINY_COS_VALUE)
-    cosine = TINY_COS_VALUE;
-  if (cosine < -TINY_COS_VALUE)
-    cosine = -TINY_COS_VALUE;
-  return cosine;
-}
-} // namespace detail
-
 /** Compute the cosine of the angle between three particles.
  *
  *  @param[in]  vec1             Vector from central particle to left particle.
@@ -58,7 +49,7 @@ inline double calc_cosine(Utils::Vector3d const &vec1,
   /* cosine of the angle between vec1 and vec2 */
   auto cos_phi = (vec1 * vec2) / std::sqrt(vec1.norm2() * vec2.norm2());
   if (sanitize_cosine) {
-    cos_phi = detail::sanitize_cosine(cos_phi);
+    cos_phi = std::clamp(cos_phi, -TINY_COS_VALUE, TINY_COS_VALUE);
   }
   return cos_phi;
 }
@@ -84,7 +75,7 @@ angle_generic_force(Utils::Vector3d const &vec1, Utils::Vector3d const &vec2,
   auto const d2 = vec2.norm();
   auto cos_phi = (vec1 * vec2) / (d1 * d2);
   if (sanitize_cosine) {
-    cos_phi = detail::sanitize_cosine(cos_phi);
+    cos_phi = std::clamp(cos_phi, -TINY_COS_VALUE, TINY_COS_VALUE);
   }
   /* force factor */
   auto const fac = forceFactor(cos_phi);

src/core/constraints/ShapeBasedConstraint.cpp

@@ -71,15 +71,15 @@ double ShapeBasedConstraint::min_dist(BoxGeometry const &box_geo,
   auto const local_mindist = std::accumulate(
       particles.begin(), particles.end(),
       std::numeric_limits<double>::infinity(),
-      [this, &box_geo](double min, Particle const &p) {
+      [this, &box_geo](double acc, Particle const &p) {
         auto const &ia_params = get_ia_param(p.type());
         if (is_active(ia_params)) {
           double dist;
           Utils::Vector3d vec;
           m_shape->calculate_dist(box_geo.folded_position(p.pos()), dist, vec);
-          return std::min(min, dist);
+          acc = std::min(acc, dist);
         }
-        return min;
+        return acc;
       });
   boost::mpi::reduce(comm_cart, local_mindist, global_mindist,
                      boost::mpi::minimum<double>(), 0);

src/core/cuda/init_cuda.cu

@@ -24,6 +24,7 @@
 #include <cuda_runtime.h>
 
 #include <cstring>
+#include <memory>
 #include <string>
 
 #if defined(OMPI_MPI_H) || defined(_MPI_H)
@@ -98,25 +99,14 @@ int cuda_get_device() {
 }
 
 bool cuda_test_device_access() {
-  int *d = nullptr;
+  auto const deleter = [](int *p) { cudaFree(reinterpret_cast<void *>(p)); };
+  int *ptr = nullptr;
   int h = 42;
-  cudaError_t err;
-
-  err = cudaMalloc((void **)&d, sizeof(int));
-  if (err != cudaSuccess) {
-    throw cuda_runtime_error_cuda(err);
-  }
-  err = cudaMemcpy(d, &h, sizeof(int), cudaMemcpyHostToDevice);
-  if (err != cudaSuccess) {
-    cudaFree(d);
-    throw cuda_runtime_error_cuda(err);
-  }
+  CUDA_CHECK(cudaMalloc(reinterpret_cast<void **>(&ptr), sizeof(int)));
+  std::unique_ptr<int, decltype(deleter)> d(ptr, deleter);
+  CUDA_CHECK(cudaMemcpy(d.get(), &h, sizeof(int), cudaMemcpyHostToDevice));
   h = 0;
-  err = cudaMemcpy(&h, d, sizeof(int), cudaMemcpyDeviceToHost);
-  cudaFree(d);
-  if (err != cudaSuccess) {
-    throw cuda_runtime_error_cuda(err);
-  }
+  CUDA_CHECK(cudaMemcpy(&h, d.get(), sizeof(int), cudaMemcpyDeviceToHost));
   return h != 42;
 }
 

src/core/electrostatics/elc.cpp

@@ -1147,7 +1147,7 @@ template <ChargeProtocol protocol, typename combined_range>
 void modify_p3m_sums(elc_data const &elc, CoulombP3M &solver,
                      combined_range const &p_q_pos_range) {
 
-  auto local_n = 0;
+  auto local_n = std::size_t{0u};
   auto local_q2 = 0.0;
   auto local_q = 0.0;
   for (auto zipped : p_q_pos_range) {
@@ -1180,7 +1180,7 @@ void modify_p3m_sums(elc_data const &elc, CoulombP3M &solver,
     }
   }
 
-  auto global_n = 0;
+  auto global_n = std::size_t{0u};
   auto global_q2 = 0.;
   auto global_q = 0.;
   boost::mpi::all_reduce(comm_cart, local_n, global_n, std::plus<>());

src/core/electrostatics/p3m.cpp

@@ -97,7 +97,7 @@
 
 template <typename FloatType, Arch Architecture>
 void CoulombP3MImpl<FloatType, Architecture>::count_charged_particles() {
-  auto local_n = 0;
+  auto local_n = std::size_t{0u};
   auto local_q2 = 0.0;
   auto local_q = 0.0;
 
@@ -463,6 +463,9 @@ double CoulombP3MImpl<FloatType, Architecture>::long_range_kernel(
 #else
   auto constexpr npt_flag = false;
 #endif
+  if (p3m.sum_qpart == 0u) {
+    return 0.;
+  }
 
   if (p3m.sum_q2 > 0.) {
     if (not has_actor_of_type<ElectrostaticLayerCorrection>(

src/core/electrostatics/p3m.hpp

@@ -49,6 +49,7 @@
 #include <utils/math/AS_erfc_part.hpp>
 
 #include <cmath>
+#include <cstddef>
 #include <numbers>
 
 /** @brief P3M solver. */
@@ -89,7 +90,7 @@ struct CoulombP3M : public Coulomb::Actor<CoulombP3M> {
    * the sum of the squared charges.
    */
   virtual void count_charged_particles() = 0;
-  virtual void count_charged_particles_elc(int, double, double) = 0;
+  virtual void count_charged_particles_elc(std::size_t, double, double) = 0;
   virtual void adapt_epsilon_elc() = 0;
 
   /**

src/core/electrostatics/p3m.impl.hpp

@@ -35,6 +35,7 @@
 
 #include <utils/Vector.hpp>
 
+#include <cstddef>
 #include <memory>
 #include <optional>
 #include <stdexcept>
@@ -45,11 +46,11 @@ template <typename FloatType>
 struct p3m_data_struct_coulomb : public p3m_data_struct<FloatType> {
   using p3m_data_struct<FloatType>::p3m_data_struct;
 
-  /** number of charged particles (only on head node). */
-  int sum_qpart = 0;
-  /** Sum of square of charges (only on head node). */
+  /** number of charged particles. */
+  std::size_t sum_qpart = 0;
+  /** Sum of square of charges. */
   double sum_q2 = 0.;
-  /** square of sum of charges (only on head node). */
+  /** square of sum of charges. */
   double square_sum_q = 0.;
 
   p3m_interpolation_cache inter_weights;
@@ -105,7 +106,7 @@ struct CoulombP3MImpl : public CoulombP3M {
   }
   void tune() override;
   void count_charged_particles() override;
-  void count_charged_particles_elc(int n, double sum_q2,
+  void count_charged_particles_elc(std::size_t n, double sum_q2,
                                    double square_sum_q) override {
     p3m.sum_qpart = n;
     p3m.sum_q2 = sum_q2;

src/core/electrostatics/p3m_gpu_cuda.cu

@@ -549,8 +549,7 @@ void assign_forces(P3MGpuData const &params,
 void p3m_gpu_init(std::shared_ptr<P3MGpuParams> &data, int cao,
                   Utils::Vector3i const &mesh, double alpha,
                   Utils::Vector3d const &box_l, std::size_t n_part) {
-  if (mesh == Utils::Vector3i::broadcast(-1))
-    throw std::runtime_error("P3M: invalid mesh size");
+  assert(mesh != Utils::Vector3i::broadcast(-1));
 
   if (not data) {
     data = std::make_shared<P3MGpuParams>();

src/core/electrostatics/scafacos_impl.cpp

@@ -121,10 +121,9 @@ void CoulombScafacosImpl::tune_r_cut() {
       break;
     }
 
+    r_min += dr;
     if (t_mid > t_min) {
-      r_max = r_min += dr;
-    } else {
-      r_min += dr;
+      r_max = r_min;
     }
   }
   assert(r_opt >= 0.);

src/core/forces_inline.hpp

@@ -274,15 +274,15 @@ inline void add_non_bonded_pair_force(
 
 /** Compute the bonded interaction force between particle pairs.
  *
+ *  @param[in] iaparams    Bonded parameters for the interaction.
  *  @param[in] p1          First particle.
  *  @param[in] p2          Second particle.
- *  @param[in] iaparams    Bonded parameters for the interaction.
  *  @param[in] dx          Vector between @p p1 and @p p2.
  *  @param[in] kernel      Coulomb force kernel.
  */
 inline std::optional<Utils::Vector3d> calc_bond_pair_force(
-    Particle const &p1, Particle const &p2,
-    Bonded_IA_Parameters const &iaparams, Utils::Vector3d const &dx,
+    Bonded_IA_Parameters const &iaparams, Particle const &p1,
+    Particle const &p2, Utils::Vector3d const &dx,
     Coulomb::ShortRangeForceKernel::kernel_type const *kernel) {
   if (auto const *iap = boost::get<FeneBond>(&iaparams)) {
     return iap->force(dx);
@@ -334,7 +334,7 @@ inline bool add_bonded_two_body_force(
       return false;
     }
   } else {
-    auto result = calc_bond_pair_force(p1, p2, iaparams, dx, kernel);
+    auto result = calc_bond_pair_force(iaparams, p1, p2, dx, kernel);
     if (result) {
       p1.force() += result.value();
       p2.force() -= result.value();

src/core/magnetostatics/dp3m.cpp

@@ -83,7 +83,7 @@
 
 template <typename FloatType, Arch Architecture>
 void DipolarP3MImpl<FloatType, Architecture>::count_magnetic_particles() {
-  int local_n = 0;
+  auto local_n = std::size_t{0u};
   double local_mu2 = 0.;
 
   for (auto const &p : get_system().cell_structure->local_particles()) {

src/core/magnetostatics/dp3m.impl.hpp

@@ -35,6 +35,7 @@
 
 #include <utils/Vector.hpp>
 
+#include <cstddef>
 #include <memory>
 #include <optional>
 #include <stdexcept>
@@ -46,9 +47,9 @@ template <typename FloatType>
 struct p3m_data_struct_dipoles : public p3m_data_struct<FloatType> {
   using p3m_data_struct<FloatType>::p3m_data_struct;
 
-  /** number of dipolar particles (only on head node). */
-  int sum_dip_part = 0;
-  /** Sum of square of magnetic dipoles (only on head node). */
+  /** number of dipolar particles. */
+  std::size_t sum_dip_part = 0;
+  /** Sum of square of magnetic dipoles. */
   double sum_mu2 = 0.;
 
   /** position shift for calculation of first assignment mesh point. */

src/core/p3m/TuningLogger.hpp

@@ -67,7 +67,7 @@ class TuningLogger {
   }
 
   void tuning_goals(double accuracy, double prefactor, double box_l,
-                    int n_particles, double sum_prop) const {
+                    std::size_t n_particles, double sum_prop) const {
     if (m_verbose) {
       std::string particle_trait;
       std::string particle_property;
@@ -82,7 +82,7 @@ class TuningLogger {
         break;
       }
       std::printf("%s tune parameters: Accuracy goal = %.5e prefactor = %.5e\n"
-                  "System: box_l = %.5e # %s part = %d %s = %.5e\n",
+                  "System: box_l = %.5e # %s part = %zu %s = %.5e\n",
                   m_name.c_str(), accuracy, prefactor, box_l,
                   particle_trait.c_str(), n_particles,
                   particle_property.c_str(), sum_prop);