core: branchless minimum-image fold and flat Verlet-list cutoff table

Wave 1 of the short-range hot-path optimization:

- Branchless cuboid minimum-image fold. Encode periodicity into a masked
  inverse box length (0 for non-periodic directions) so the per-component
  fold reduces to `dx - rint(dx * inv_masked) * L` with no branch. `rint`
  maps to a single rounding instruction. Results match the previous
  round-based fold across the pair-loop input domain (separations below
  1.5 box lengths); confirmed bitwise-identical on the canonical lj and
  p3m trajectories.

- Flat per-type-pair squared-cutoff table in VerletCriterion. The
  per-candidate cutoff query in the Verlet-list build becomes a dense
  table load instead of walking the InteractionsNonBonded pointer table;
  inactive pairs store a negative sentinel so the distance comparison
  rejects them without a separate activity check.

- Hoist cuboid box parameters by value into the Verlet-build kernels via
  CuboidMinimumImage, instead of chasing the BoxGeometry reference for the
  box lengths on every candidate pair.

- Force-inline Utils::Vector operator+= / operator-= (previously outlined
  as .isra clones called from inside the pair kernel).

Canonical identity bitwise-preserved (lj, p3m); unit tests pass.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>

Author: RudolfWeeber

src/core/BoxGeometry.hpp

@@ -35,25 +35,25 @@
 namespace detail {
 /**
  * @brief Get the minimum-image distance between two coordinates.
- * @param a               Coordinate of the terminal point.
- * @param b               Coordinate of the initial point.
- * @param box_length      Box length.
- * @param box_length_inv  Inverse box length
- * @param box_length_half Half box length
- * @param periodic        Box periodicity.
+ *
+ * Branchless fold: the periodicity is encoded in the masked inverse box
+ * length (0 for non-periodic directions, where <tt>rint</tt> then yields a
+ * zero image shift). Uses <tt>rint</tt> (round half to even) rather than
+ * <tt>round</tt>, so it maps to a single rounding instruction; the two only
+ * differ for separations of exactly half a box length, where both images
+ * are equidistant.
+ *
+ * @param a                      Coordinate of the terminal point.
+ * @param b                      Coordinate of the initial point.
+ * @param box_length             Box length.
+ * @param box_length_inv_masked  Inverse box length if periodic, 0 otherwise.
  * @return Shortest distance from @p b to @p a across periodic images,
  *         i.e. <tt>a - b</tt>. Can be negative.
  */
 template <typename T>
-T get_mi_coord(T a, T b, T box_length, T box_length_inv, T box_length_half,
-               bool periodic) {
+T get_mi_coord_masked(T a, T b, T box_length, T box_length_inv_masked) {
   auto const dx = a - b;
-
-  if (periodic && (std::abs(dx) > box_length_half)) {
-    return dx - std::round(dx * box_length_inv) * box_length;
-  }
-
-  return dx;
+  return dx - std::rint(dx * box_length_inv_masked) * box_length;
 }
 
 /**
@@ -66,8 +66,8 @@ T get_mi_coord(T a, T b, T box_length, T box_length_inv, T box_length_half,
  *         i.e. <tt>a - b</tt>. Can be negative.
  */
 template <typename T> T get_mi_coord(T a, T b, T box_length, bool periodic) {
-  return get_mi_coord(a, b, box_length, 1. / box_length, 0.5 * box_length,
-                      periodic);
+  return get_mi_coord_masked(a, b, box_length,
+                             periodic ? T{1.} / box_length : T{0.});
 }
 
 /** @brief Calculate image box shift vector.
@@ -95,6 +95,37 @@ inline auto unfolded_position(Utils::Vector3d const &pos,
 
 enum class BoxType { CUBOID = 0, LEES_EDWARDS = 1 };
 
+/**
+ * @brief Cuboid minimum-image fold parameters for hot pair loops.
+ *
+ * Capture an instance by value in a kernel to hoist the box data out of the
+ * pair loop: member loads then come from the kernel's own frame and the
+ * compiler can keep them in registers, instead of re-reading them through a
+ * @ref BoxGeometry reference for every pair. Only valid for cuboid boxes;
+ * Lees-Edwards boxes need the full @ref BoxGeometry::get_mi_vector.
+ */
+class CuboidMinimumImage {
+  Utils::Vector3d m_length;
+  Utils::Vector3d m_length_inv_masked;
+
+public:
+  CuboidMinimumImage(Utils::Vector3d const &length,
+                     Utils::Vector3d const &length_inv_masked)
+      : m_length(length), m_length_inv_masked(length_inv_masked) {}
+
+  /** @brief Squared minimum-image distance between two coordinates. */
+  ESPRESSO_ATTR_ALWAYS_INLINE inline double
+  dist2(Utils::Vector3d const &a, Utils::Vector3d const &b) const {
+    double acc = 0.;
+    for (auto c = 0u; c < 3u; ++c) {
+      auto const dx = detail::get_mi_coord_masked(a[c], b[c], m_length[c],
+                                                  m_length_inv_masked[c]);
+      acc += dx * dx;
+    }
+    return acc;
+  }
+};
+
 class BoxGeometry {
 public:
   BoxGeometry() {
@@ -121,6 +152,10 @@ class BoxGeometry {
   Utils::Vector3d m_length = {1., 1., 1.};
   /** Inverse side lengths of the box */
   Utils::Vector3d m_length_inv = {1., 1., 1.};
+  /** Inverse side lengths for periodic directions, 0 for non-periodic ones.
+   *  Folding the periodicity into the inverse length makes the cuboid
+   *  minimum-image fold branchless (see @ref detail::get_mi_coord_masked). */
+  Utils::Vector3d m_length_inv_masked = {1., 1., 1.};
   /** Half side lengths of the box */
   Utils::Vector3d m_length_half = {0.5, 0.5, 0.5};
 
@@ -134,7 +169,10 @@ class BoxGeometry {
    * @param coord The coordinate to set the periodicity for.
    * @param val True if this direction should be periodic.
    */
-  void set_periodic(unsigned coord, bool val) { m_periodic.set(coord, val); }
+  void set_periodic(unsigned coord, bool val) {
+    m_periodic.set(coord, val);
+    m_length_inv_masked[coord] = val ? m_length_inv[coord] : 0.;
+  }
 
   /**
    * @brief Check periodicity in direction.
@@ -173,6 +211,9 @@ class BoxGeometry {
     assert(box_l > Utils::Vector3d::broadcast(0.));
     m_length = box_l;
     m_length_inv = {1. / box_l[0], 1. / box_l[1], 1. / box_l[2]};
+    for (auto c = 0u; c < 3u; ++c) {
+      m_length_inv_masked[c] = m_periodic[c] ? m_length_inv[c] : 0.;
+    }
     m_length_half = 0.5 * box_l;
   }
 
@@ -193,8 +234,14 @@ class BoxGeometry {
   template <typename T> T inline get_mi_coord(T a, T b, unsigned coord) const {
     assert(coord <= 2u);
 
-    return detail::get_mi_coord(a, b, m_length[coord], m_length_inv[coord],
-                                m_length_half[coord], m_periodic[coord]);
+    return detail::get_mi_coord_masked(
+        a, b, static_cast<T>(m_length[coord]),
+        static_cast<T>(m_length_inv_masked[coord]));
+  }
+
+  /** @brief Cuboid minimum-image fold parameters for hoisting into kernels. */
+  auto cuboid_minimum_image() const {
+    return CuboidMinimumImage{m_length, m_length_inv_masked};
   }
 
   /**

src/core/nonbonded_interactions/VerletCriterion.hpp

@@ -30,6 +30,11 @@
 #include <utils/index.hpp>
 #include <utils/math/sqr.hpp>
 
+#include <algorithm>
+#include <cassert>
+#include <cstddef>
+#include <vector>
+
 struct GetNonbondedCutoff {
   GetNonbondedCutoff(System::System const &system) : m_system{system} {}
   auto operator()(int type_i, int type_j) const {
@@ -54,7 +59,14 @@ template <typename CutoffGetter = GetNonbondedCutoff> class VerletCriterion {
       return inactive_cutoff;
     return Utils::sqr(x + m_skin);
   }
-  CutoffGetter get_nonbonded_cutoff;
+  /** Dense row-major table of squared effective (cutoff + skin) values per
+   *  type pair, @ref inactive_cutoff for inactive pairs. The per-type-pair
+   *  cutoff query runs once per candidate pair in the Verlet-list build, so
+   *  it must be a plain load instead of a walk through the
+   *  @ref InteractionsNonBonded pointer table.
+   */
+  std::vector<double> m_eff_cut2_table;
+  int m_n_types;
 
 public:
   VerletCriterion(System::System const &system, double skin, double max_cut,
@@ -63,8 +75,27 @@ template <typename CutoffGetter = GetNonbondedCutoff> class VerletCriterion {
       : m_skin(skin), m_eff_max_cut2(eff_cutoff_sqr(max_cut)),
         m_eff_coulomb_cut2(eff_cutoff_sqr(coulomb_cut)),
         m_eff_dipolar_cut2(eff_cutoff_sqr(dipolar_cut)),
-        m_collision_cut2(eff_cutoff_sqr(collision_detection_cutoff)),
-        get_nonbonded_cutoff(system) {}
+        m_collision_cut2(eff_cutoff_sqr(collision_detection_cutoff)) {
+    CutoffGetter const get_nonbonded_cutoff(system);
+    auto const max_type = system.nonbonded_ias->get_max_seen_particle_type();
+    m_n_types = std::max(max_type + 1, 1);
+    m_eff_cut2_table.assign(static_cast<std::size_t>(m_n_types) *
+                                static_cast<std::size_t>(m_n_types),
+                            inactive_cutoff);
+    for (int type_i = 0; type_i <= max_type; ++type_i) {
+      for (int type_j = type_i; type_j <= max_type; ++type_j) {
+        auto const cut = get_nonbonded_cutoff(type_i, type_j);
+        auto const eff_cut2 =
+            (cut == inactive_cutoff) ? inactive_cutoff : Utils::sqr(cut + skin);
+        m_eff_cut2_table[static_cast<std::size_t>(type_i) *
+                             static_cast<std::size_t>(m_n_types) +
+                         static_cast<std::size_t>(type_j)] = eff_cut2;
+        m_eff_cut2_table[static_cast<std::size_t>(type_j) *
+                             static_cast<std::size_t>(m_n_types) +
+                         static_cast<std::size_t>(type_i)] = eff_cut2;
+      }
+    }
+  }
 
   bool operator()(const Particle &p1, const Particle &p2, double dist2) const {
     if (dist2 > m_eff_max_cut2)
@@ -88,9 +119,15 @@ template <typename CutoffGetter = GetNonbondedCutoff> class VerletCriterion {
       return true;
 #endif
 
-    // Within short-range distance (including dpd and the like)
-    auto const ia_cut = get_nonbonded_cutoff(p1.type(), p2.type());
-    return (ia_cut != inactive_cutoff) &&
-           (dist2 <= Utils::sqr(ia_cut + m_skin));
+    // Within short-range distance (including dpd and the like). Inactive
+    // pairs hold inactive_cutoff (negative) in the table, so the comparison
+    // rejects them without a separate activity check.
+    auto const type_i = p1.type();
+    auto const type_j = p2.type();
+    assert(type_i >= 0 and type_i < m_n_types);
+    assert(type_j >= 0 and type_j < m_n_types);
+    return dist2 <= m_eff_cut2_table[static_cast<std::size_t>(type_i) *
+                                         static_cast<std::size_t>(m_n_types) +
+                                     static_cast<std::size_t>(type_j)];
   }
 };

src/core/short_range_cabana.hpp

@@ -121,14 +121,27 @@ link_cell_kokkos(std::span<Cell *const> cells, BoxGeometry const &box_geo,
     }
   }
 
+  // Hoist the cuboid minimum-image parameters by value: the fold runs once
+  // per candidate pair and must not chase the BoxGeometry reference for its
+  // box lengths every time. Lees-Edwards boxes take the full BoxGeometry
+  // path (shear offset handling).
+  bool const has_lees_edwards = box_geo.type() == BoxType::LEES_EDWARDS;
+  auto const cuboid_minimum_image = box_geo.cuboid_minimum_image();
+  auto const minimum_image_dist2 =
+      [&box_geo, has_lees_edwards, cuboid_minimum_image](
+          Utils::Vector3d const &a, Utils::Vector3d const &b) {
+        return has_lees_edwards ? box_geo.get_mi_dist2(a, b)
+                                : cuboid_minimum_image.dist2(a, b);
+      };
+
   // Iterate the cells' store-ROW bags directly and REBIND two cached views
   // (p1 + partner) per work item via Particle::attach_to_store, instead of
   // driving RowParticleRange iterators (each embeds a Particle by value, so
   // std::next(it) and per-neighbour range begin()/end() would build fresh
   // Particles). One reused view per role per work item (one cell per Kokkos
   // work item) is thread-safe. Iteration ORDER is unchanged.
-  auto intra_kernel = [&cells, &box_geo, &verlet_criterion, &id_to_index,
-                       &intra_operator, &interleaved_positions,
+  auto intra_kernel = [&cells, minimum_image_dist2, &verlet_criterion,
+                       &id_to_index, &intra_operator, &interleaved_positions,
                        max_id](const int i) {
     auto &store = cells[i]->store();
     // Contiguous store-row range; clean store, so index directly.
@@ -183,7 +196,7 @@ link_cell_kokkos(std::span<Cell *const> cells, BoxGeometry const &box_geo,
                   Utils::Vector3d{p2_base[0u], p2_base[pos_comp_stride],
                                   p2_base[2u * pos_comp_stride]};
               if (verlet_criterion(p1, p2,
-                                   box_geo.get_mi_dist2(p1_pos, p2_pos))) {
+                                   minimum_image_dist2(p1_pos, p2_pos))) {
                 auto const jj = id_to_index(id_column[row_b]);
                 if (jj >= 0) {
                   intra_operator(ii, jj);
@@ -196,8 +209,8 @@ link_cell_kokkos(std::span<Cell *const> cells, BoxGeometry const &box_geo,
     }
   };
 
-  auto inter_kernel = [&cells, &box_geo, &verlet_criterion, &id_to_index,
-                       &inter_operator, &interleaved_positions,
+  auto inter_kernel = [&cells, minimum_image_dist2, &verlet_criterion,
+                       &id_to_index, &inter_operator, &interleaved_positions,
                        max_id](const int i) {
     auto &store = cells[i]->store();
     // Contiguous store-row range; clean store, so index directly.
@@ -247,7 +260,7 @@ link_cell_kokkos(std::span<Cell *const> cells, BoxGeometry const &box_geo,
                     Utils::Vector3d{p2_base[0u], p2_base[pos_comp_stride],
                                     p2_base[2u * pos_comp_stride]};
                 if (verlet_criterion(p1, p2,
-                                     box_geo.get_mi_dist2(p1_pos, p2_pos))) {
+                                     minimum_image_dist2(p1_pos, p2_pos))) {
                   auto const jj = id_to_index(id_column[row_k]);
                   if (jj >= 0) {
                     inter_operator(ii, jj);

src/utils/include/utils/Vector.hpp

@@ -30,6 +30,7 @@
 #include <boost/qvm/vec_traits.hpp>
 
 #include "utils/Array.hpp"
+#include "utils/attributes.hpp"
 
 #include <algorithm>
 #include <cassert>
@@ -270,7 +271,8 @@ auto operator+(Vector<T, N> const &a, Vector<U, N> const &b) {
 }
 
 template <std::size_t N, typename T>
-auto &operator+=(Vector<T, N> &a, Vector<T, N> const &b) {
+ESPRESSO_ATTR_ALWAYS_INLINE inline auto &operator+=(Vector<T, N> &a,
+                                                    Vector<T, N> const &b) {
   std::ranges::transform(a, b, std::begin(a), std::plus<T>());
   return a;
 }
@@ -288,7 +290,8 @@ Vector<T, N> operator-(Vector<T, N> const &a) {
 }
 
 template <std::size_t N, typename T>
-Vector<T, N> &operator-=(Vector<T, N> &a, Vector<T, N> const &b) {
+ESPRESSO_ATTR_ALWAYS_INLINE inline Vector<T, N> &
+operator-=(Vector<T, N> &a, Vector<T, N> const &b) {
   std::ranges::transform(a, b, std::begin(a), std::minus<T>());
   return a;
 }

src/utils/include/utils/attributes.hpp

@@ -0,0 +1,35 @@
+/*
+ * Copyright (C) 2010-2026 The ESPResSo project
+ *
+ * This file is part of ESPResSo.
+ *
+ * ESPResSo is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * ESPResSo is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#pragma once
+
+/** \file
+ *  Compiler-attribute macros shared across utils headers.
+ *
+ *  Uses the same macro name and guard as the core attributes header, so
+ *  either header can be included first without redefinition.
+ */
+
+#ifndef ESPRESSO_ATTR_ALWAYS_INLINE
+#if defined(__GNUG__) or defined(__clang__)
+#define ESPRESSO_ATTR_ALWAYS_INLINE [[gnu::always_inline]]
+#else
+#define ESPRESSO_ATTR_ALWAYS_INLINE
+#endif
+#endif