diff --git a/maintainer/benchmarks/ghost_ab.sh b/maintainer/benchmarks/ghost_ab.sh
new file mode 100755
index 0000000000..f16dcefc56
--- /dev/null
+++ b/maintainer/benchmarks/ghost_ab.sh
@@ -0,0 +1,26 @@
+#!/usr/bin/env bash
+# Usage: ghost_ab.sh <label>   (run from the build dir)
+# Runs lj.py and p3m.py on 1/2/4/8 ranks after checking the machine is idle.
+set -euo pipefail
+label="${1:?need a label, e.g. baseline or async}"
+build="$(pwd)"
+
+# --- machine-idle gate (per project constraint) ---
+load1=$(awk '{print $1}' /proc/loadavg)
+ncpu="$(nproc)"
+others=$(who | awk -v me="$(id -un)" '$1 != me {print $1}' | sort -u | wc -l)
+if (( $(echo "${load1} > 2.0" | bc -l) )) || (( others > 0 )); then
+  echo "REFUSING to benchmark: load1=${load1}, other users=${others}, ncpu=${ncpu}." >&2
+  echo "Re-run when the machine is idle." >&2
+  exit 3
+fi
+
+out="${build}/ghost_ab_${label}.csv"
+: > "${out}"
+for ranks in 1 2 4 8; do
+  mpiexec --bind-to core -n "${ranks}" ./pypresso ../maintainer/benchmarks/lj.py \
+    --particles_per_core=1000 --volume_fraction=0.50 --output="${out}"
+  mpiexec --bind-to core -n "${ranks}" ./pypresso ../maintainer/benchmarks/p3m.py \
+    --particles_per_core=1000 --volume_fraction=0.25 --output="${out}"
+done
+echo "wrote ${out}"
diff --git a/src/core/CMakeLists.txt b/src/core/CMakeLists.txt
index ef976cb147..384bd37679 100644
--- a/src/core/CMakeLists.txt
+++ b/src/core/CMakeLists.txt
@@ -30,7 +30,9 @@ add_library(
   errorhandling.cpp
   forces.cpp
   forces_init.cpp
-  ghosts.cpp
+  ghosts/HaloExchange.cpp
+  ghosts/HaloPlanValidator.cpp
+  ghosts/particle_packing.cpp
   integrate.cpp
   npt.cpp
   particle_node.cpp
diff --git a/src/core/caliper_utils.hpp b/src/core/caliper_utils.hpp
new file mode 100644
index 0000000000..1aa1314895
--- /dev/null
+++ b/src/core/caliper_utils.hpp
@@ -0,0 +1,224 @@
+/*
+ * Copyright (C) 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 caliper_utils.hpp
+ * @brief Zero-overhead Caliper guards for the inactive (no @c CALI_CONFIG)
+ * case.
+ *
+ * When Caliper is compiled in but @c CALI_CONFIG is not set (the common
+ * production case), the standard @c CALI_MARK_BEGIN/END and
+ * @c CALI_CXX_MARK_FUNCTION macros still pay the full Caliper entry cost
+ * (thread-blackboard update, siglock acquisition).  On a hot path with
+ * many markers per step this amounts to measurable wall-clock overhead.
+ *
+ * This header provides guarded replacements:
+ * - @c espresso_cali_active() — checked once per process; in steady state this
+ *   is a single byte load + branch.
+ * - @c ESPRESSO_CALI_MARK_FUNCTION — RAII guard equivalent to
+ *   @c CALI_CXX_MARK_FUNCTION but zero-cost when inactive.
+ * - @c ESPRESSO_CALI_MARK_BEGIN(name) / @c ESPRESSO_CALI_MARK_END(name)
+ * - @c EspressoCaliLoop — RAII loop wrapper replacing the
+ *   @c ESPRESSO_CALI_MARK_LOOP_BEGIN / @c ESPRESSO_CALI_MARK_LOOP_ITERATION /
+ *   @c ESPRESSO_CALI_MARK_LOOP_END macro triplet.
+ *
+ * All macros and types are no-ops / empty when @c ESPRESSO_CALIPER is not
+ * defined.
+ *
+ * @par Activation path limitation
+ * @c espresso_cali_active() detects Caliper activation **solely** via the
+ * @c CALI_CONFIG environment variable.  Activation through
+ * @c CALI_SERVICES_ENABLE, Caliper config files, or programmatic
+ * @c cali::ConfigManager::start() is **not** detected — hot-path regions will
+ * silently produce no output for those activation paths.  @c CALI_CONFIG is
+ * the supported activation path for ESPResSo profiling (used by
+ * @c testsuite/python/caliper.py).  This restriction is intentional: checking
+ * @c CALI_CONFIG once and caching the result eliminates the Caliper entry cost
+ * (siglock + thread-blackboard update) on every hot-path marker when profiling
+ * is inactive.
+ */
+
+#ifdef ESPRESSO_CALIPER
+
+#include <caliper/cali.h>
+#include <cstdlib>
+
+/**
+ * @brief Return true if Caliper is configured for this process.
+ *
+ * Reads @c CALI_CONFIG from the environment exactly once (on the first call)
+ * and caches the result.  Subsequent calls pay only one byte load + branch.
+ *
+ * @c inline ensures a single shared static across all translation units
+ * (C++ ODR for inline functions with static locals).
+ *
+ * @note Only @c CALI_CONFIG activation is detected; see the file-level
+ *       documentation for the rationale and limitation.
+ */
+inline bool espresso_cali_active() noexcept {
+  static const bool active = (std::getenv("CALI_CONFIG") != nullptr);
+  return active;
+}
+
+/**
+ * @brief RAII region guard: begin on construction, end on destruction.
+ *
+ * Conditionally calls @c cali_begin_region / @c cali_end_region only when
+ * @c espresso_cali_active() is true.  Correctly handles all exit paths.
+ */
+struct EspressoCaliRegion {
+  const char *name_;
+  bool active_;
+
+  EspressoCaliRegion(const char *name, bool active)
+      : name_(name), active_(active) {
+    if (active_)
+      cali_begin_region(name_);
+  }
+  ~EspressoCaliRegion() {
+    if (active_)
+      cali_end_region(name_);
+  }
+  EspressoCaliRegion(const EspressoCaliRegion &) = delete;
+  EspressoCaliRegion &operator=(const EspressoCaliRegion &) = delete;
+};
+
+/**
+ * @brief RAII iteration annotation for one pass through a loop body.
+ *
+ * Returned by @c EspressoCaliLoop::iteration().  When active, calls
+ * @c cali_begin_int on the loop's iteration attribute on construction and
+ * @c cali_end_byid on destruction — identical to what @c cali::Loop::Iteration
+ * does internally.  When inactive, both operations are no-ops.
+ *
+ * The object must be stored at the top of the loop body so that its lifetime
+ * spans the loop body:
+ * @code
+ *   for (int step = 0; ...) {
+ *     auto cali_iter = cali_loop.iteration(step);
+ *     // ...
+ *   }  // cali_iter destroyed here, closing the iteration region
+ * @endcode
+ */
+struct EspressoCaliIteration {
+  cali_id_t iter_attr_;
+  bool active_;
+
+  EspressoCaliIteration(cali_id_t iter_attr, int iter, bool active)
+      : iter_attr_(iter_attr), active_(active) {
+    if (active_)
+      cali_begin_int(iter_attr_, iter);
+  }
+  ~EspressoCaliIteration() {
+    if (active_)
+      cali_end(iter_attr_);
+  }
+  EspressoCaliIteration(const EspressoCaliIteration &) = delete;
+  EspressoCaliIteration &operator=(const EspressoCaliIteration &) = delete;
+  // Moveable so it can be returned from iteration() with NRVO/move.
+  EspressoCaliIteration(EspressoCaliIteration &&other) noexcept
+      : iter_attr_(other.iter_attr_), active_(other.active_) {
+    other.active_ = false;
+  }
+  EspressoCaliIteration &operator=(EspressoCaliIteration &&) = delete;
+};
+
+/**
+ * @brief RAII loop wrapper replacing the @c CALI_CXX_MARK_LOOP_BEGIN /
+ * @c CALI_CXX_MARK_LOOP_ITERATION / @c CALI_CXX_MARK_LOOP_END macro triplet.
+ *
+ * When @c espresso_cali_active() is true at construction time, creates the
+ * same @c cali.loop region and @c iteration#name attribute as the raw Caliper
+ * macros would.  When inactive, construction, iteration(), and destruction are
+ * all no-ops with no Caliper calls.
+ *
+ * Usage:
+ * @code
+ *   EspressoCaliLoop cali_loop("Integration loop");
+ *   for (int step = 0; step < n_steps; ++step) {
+ *     auto cali_iter = cali_loop.iteration(step);  // RAII: closed at }
+ *     // ... loop body ...
+ *   }
+ *   // cali_loop destructor ends the loop region
+ * @endcode
+ */
+struct EspressoCaliLoop {
+  cali::Loop *loop_;
+  cali_id_t iter_attr_;
+
+  explicit EspressoCaliLoop(const char *name)
+      : loop_(nullptr), iter_attr_(CALI_INV_ID) {
+    if (espresso_cali_active()) {
+      loop_ = new cali::Loop(name);
+      iter_attr_ = cali_make_loop_iteration_attribute(name);
+    }
+  }
+  ~EspressoCaliLoop() {
+    if (loop_)
+      loop_->end();
+    delete loop_;
+  }
+  EspressoCaliLoop(const EspressoCaliLoop &) = delete;
+  EspressoCaliLoop &operator=(const EspressoCaliLoop &) = delete;
+
+  /**
+   * @brief Return an RAII iteration annotation for the current step.
+   *
+   * The returned @c EspressoCaliIteration must be kept alive for the duration
+   * of the loop body; store it in a named local variable.
+   */
+  EspressoCaliIteration iteration(int iter) const {
+    return EspressoCaliIteration(iter_attr_, iter, loop_ != nullptr);
+  }
+};
+
+// NOLINTBEGIN(cppcoreguidelines-macro-usage)
+
+/**
+ * @brief Guarded drop-in replacement for @c CALI_CXX_MARK_FUNCTION.
+ *
+ * Captures @c __func__ at the call site so that the region name is the
+ * enclosing function name, not the destructor name.
+ */
+#define ESPRESSO_CALI_MARK_FUNCTION                                            \
+  EspressoCaliRegion CALI_CREATE_VAR_NAME(__espresso_cali_fn, __LINE__)(       \
+      __func__, ::espresso_cali_active())
+
+/**
+ * @brief Guarded @c CALI_MARK_BEGIN — no-op when inactive.
+ */
+#define ESPRESSO_CALI_MARK_BEGIN(name)                                         \
+  do {                                                                         \
+    if (::espresso_cali_active())                                              \
+      CALI_MARK_BEGIN(name);                                                   \
+  } while (false)
+
+/**
+ * @brief Guarded @c CALI_MARK_END — no-op when inactive.
+ */
+#define ESPRESSO_CALI_MARK_END(name)                                           \
+  do {                                                                         \
+    if (::espresso_cali_active())                                              \
+      CALI_MARK_END(name);                                                     \
+  } while (false)
+
+// NOLINTEND(cppcoreguidelines-macro-usage)
+
+#endif // ESPRESSO_CALIPER
diff --git a/src/core/cell_system/AtomDecomposition.cpp b/src/core/cell_system/AtomDecomposition.cpp
index 3db794b22c..f56900cc0f 100644
--- a/src/core/cell_system/AtomDecomposition.cpp
+++ b/src/core/cell_system/AtomDecomposition.cpp
@@ -23,10 +23,15 @@
 
 #include "cell_system/Cell.hpp"
 
+#include "ghosts/HaloPlan.hpp"
+#include "ghosts/HaloPlanValidator.hpp"
+#include "ghosts/mark_boundary_cells.hpp"
+
 #include <utils/Vector.hpp>
 
 #include <boost/mpi/collectives/all_to_all.hpp>
 
+#include <cassert>
 #include <cstddef>
 #include <limits>
 #include <utility>
@@ -52,45 +57,40 @@ void AtomDecomposition::configure_neighbors() {
   local().m_neighbors = Neighbors<Cell *>(red_neighbors, black_neighbors);
 }
 
-GhostCommunicator AtomDecomposition::prepare_comm() {
-  /* no need for comm for only 1 node */
+GhostComm::HaloPlan AtomDecomposition::make_halo_plan() {
+  using GhostComm::CollectivePattern;
+  using GhostComm::CollectiveSection;
+  using GhostComm::HaloPlan;
+
+  HaloPlan plan;
+  plan.comm = m_comm;
+
+  // Single rank: no communication needed; collective section is None.
   if (m_comm.size() == 1) {
-    return GhostCommunicator{m_comm, 0};
+    plan.collective = CollectiveSection{CollectivePattern::None, {}};
+    return plan;
   }
 
-  auto ghost_comm =
-      GhostCommunicator{m_comm, static_cast<std::size_t>(m_comm.size())};
-  /* every node has its dedicated comm step */
-  for (int n = 0; n < m_comm.size(); n++) {
-    ghost_comm.communications[n].part_lists.resize(1);
-    ghost_comm.communications[n].part_lists[0] = &(cells.at(n).particles());
-    ghost_comm.communications[n].node = n;
+  // One cell pointer per rank: cells[root] is the ParticleList for that root.
+  // The engine uses op.direction to pick Broadcast (Push) or ReduceSum
+  // (Reduce) at run time, so we store Broadcast as the canonical marker that
+  // this section is active.  run_collective reads op.direction to decide which
+  // MPI collective to invoke.
+  std::vector<ParticleList *> cell_ptrs;
+  cell_ptrs.reserve(static_cast<std::size_t>(m_comm.size()));
+  for (int n = 0; n < m_comm.size(); ++n) {
+    cell_ptrs.push_back(&cells.at(static_cast<std::size_t>(n)).particles());
   }
-
-  return ghost_comm;
+  plan.collective =
+      CollectiveSection{CollectivePattern::Broadcast, std::move(cell_ptrs)};
+  return plan;
 }
 
 void AtomDecomposition::configure_comms() {
-  m_exchange_ghosts_comm = prepare_comm();
-  m_collect_ghost_force_comm = prepare_comm();
-
-  if (m_comm.size() > 1) {
-    for (int n = 0; n < m_comm.size(); n++) {
-      /* use the prefetched send buffers. Node 0 transmits first and never
-       * prefetches. */
-      if (m_comm.rank() == 0 || m_comm.rank() != n) {
-        m_exchange_ghosts_comm.communications[n].type = GHOST_BCST;
-      } else {
-        m_exchange_ghosts_comm.communications[n].type =
-            GHOST_BCST | GHOST_PREFETCH;
-      }
-      m_collect_ghost_force_comm.communications[n].type = GHOST_RDCE;
-    }
-    /* first round: all nodes except the first one prefetch their send data */
-    if (m_comm.rank() != 0) {
-      m_exchange_ghosts_comm.communications[0].type |= GHOST_PREFETCH;
-    }
-  }
+  m_halo_plan = make_halo_plan();
+  // NOTE: validation is deferred to the constructor, AFTER mark_cells() has
+  // populated local_cells()/ghost_cells(). Validating here would check empty
+  // spans (vacuously) since mark_cells() runs later.
 }
 
 void AtomDecomposition::mark_cells() {
@@ -156,6 +156,40 @@ AtomDecomposition::AtomDecomposition(boost::mpi::communicator comm,
   configure_neighbors();
   /* fill local and ghost cell lists */
   mark_cells();
+  /* classify local cells as interior or boundary.
+   *
+   * AtomDecomposition has no spatial locality: the single local cell
+   * interacts with every other rank's cell and there is no subset of
+   * particles whose force contributions are guaranteed to arrive before
+   * the velocity update.  Interior is therefore always empty and all local
+   * cells are boundary.  mark_boundary_cells() handles the ghost-neighbour
+   * case (multi-rank); the explicit loop below catches the single-rank case
+   * where there are no ghost cells and no neighbours at all.
+   */
+  GhostComm::mark_boundary_cells(AtomDecomposition::local_cells(),
+                                 AtomDecomposition::ghost_cells());
+  for (Cell *c : AtomDecomposition::local_cells()) {
+    c->m_is_boundary = true;
+  }
+#ifdef ESPRESSO_ADDITIONAL_CHECKS
+  // Validate now that local_cells()/ghost_cells() are populated by
+  // mark_cells().
+  assert(GhostComm::report_violations(
+      GhostComm::validate_halo_plan(m_halo_plan,
+                                    AtomDecomposition::local_cells(),
+                                    AtomDecomposition::ghost_cells()),
+      "AtomDecomposition"));
+  // NOTE: validate_halo_plan_symmetry is NOT called here.
+  // During checkpoint loading, decompositions are transiently rebuilt while
+  // maximal_cutoff is rank-divergent (ranks may have different cell grids for a
+  // brief window before the next consistent rebuild).  The transient plan is
+  // never used — it is immediately replaced — so the asymmetry is harmless.
+  // A construction-time collective all_to_all inside a ctor is also dangerous:
+  // if one rank aborts the others block forever in the collective.
+  // Symmetry is instead validated at FIRST USE of the plan in
+  // halo_exchange_start (see GhostComm::halo_exchange_start in
+  // HaloExchange.cpp).
+#endif
 }
 
 Utils::Vector3d AtomDecomposition::max_cutoff() const {
diff --git a/src/core/cell_system/AtomDecomposition.hpp b/src/core/cell_system/AtomDecomposition.hpp
index 03994722d5..83fb302e2a 100644
--- a/src/core/cell_system/AtomDecomposition.hpp
+++ b/src/core/cell_system/AtomDecomposition.hpp
@@ -27,7 +27,7 @@
 
 #include "BoxGeometry.hpp"
 #include "Particle.hpp"
-#include "ghosts.hpp"
+#include "ghosts/HaloPlan.hpp"
 
 #include <utils/Vector.hpp>
 
@@ -57,8 +57,10 @@ class AtomDecomposition : public ParticleDecomposition {
   std::vector<Cell *> m_local_cells;
   std::vector<Cell *> m_ghost_cells;
 
-  GhostCommunicator m_exchange_ghosts_comm;
-  GhostCommunicator m_collect_ghost_force_comm;
+  /// Holds ParticleList* into this decomposition's cells — value-copying this
+  /// object leaves these pointers dangling. TODO(Task 1.7): make non-copyable
+  /// or rebuild-on-copy.
+  GhostComm::HaloPlan m_halo_plan;
 
   BoxGeometry const &m_box;
 
@@ -68,12 +70,7 @@ class AtomDecomposition : public ParticleDecomposition {
 
   void resort(bool global_flag, std::vector<ParticleChange> &diff) override;
 
-  GhostCommunicator const &exchange_ghosts_comm() const override {
-    return m_exchange_ghosts_comm;
-  }
-  GhostCommunicator const &collect_ghost_force_comm() const override {
-    return m_collect_ghost_force_comm;
-  }
+  GhostComm::HaloPlan const *halo_plan() const override { return &m_halo_plan; }
 
   std::span<Cell *const> local_cells() const override { return m_local_cells; }
   std::span<Cell *const> ghost_cells() const override { return m_ghost_cells; }
@@ -131,7 +128,11 @@ class AtomDecomposition : public ParticleDecomposition {
   }
 
   void configure_neighbors();
-  GhostCommunicator prepare_comm();
+
+  /**
+   * @brief Build the plan-based halo plan (collective broadcast/reduce path).
+   */
+  GhostComm::HaloPlan make_halo_plan();
 
   /**
    * @brief Setup ghost communicators.
diff --git a/src/core/cell_system/Cell.hpp b/src/core/cell_system/Cell.hpp
index da36c8a005..e48ee7a5f3 100644
--- a/src/core/cell_system/Cell.hpp
+++ b/src/core/cell_system/Cell.hpp
@@ -112,4 +112,15 @@ class Cell {
    * @brief All neighbors of the cell.
    */
   neighbors_type &neighbors() { return m_neighbors; }
+
+  /**
+   * @brief Interior/boundary classification.
+   *
+   * A local cell is *boundary* iff at least one of its neighbors is a ghost
+   * cell.  Set by GhostComm::mark_boundary_cells() after neighbor setup.
+   * Default false (interior).
+   */
+  bool m_is_boundary = false;
+
+  bool is_boundary() const { return m_is_boundary; }
 };
diff --git a/src/core/cell_system/CellStructure.cpp b/src/core/cell_system/CellStructure.cpp
index 786933fad5..0b2f4237d4 100644
--- a/src/core/cell_system/CellStructure.cpp
+++ b/src/core/cell_system/CellStructure.cpp
@@ -34,6 +34,7 @@
 #include "cell_system/CellStructureType.hpp"
 #include "communication.hpp"
 #include "ghosts.hpp"
+#include "ghosts/HaloExchange.hpp"
 #include "integrators/Propagation.hpp"
 #include "lees_edwards/lees_edwards.hpp"
 #include "particle_enumeration.hpp"
@@ -46,7 +47,7 @@
 #include <utils/math/sqr.hpp>
 
 #ifdef ESPRESSO_CALIPER
-#include <caliper/cali.h>
+#include "caliper_utils.hpp"
 #endif
 
 #include <boost/mpi/collectives/all_reduce.hpp>
@@ -72,6 +73,9 @@
 #include <vector>
 
 CellStructure::~CellStructure() {
+  assert(not m_pending_ghost_reduce.has_value() &&
+         "~CellStructure: ghost force reduction still in flight at destruction "
+         "— ghosts_reduce_forces_finish() was not called");
   clear_local_properties();
   // Kokkos handle can only be freed after all Cabana containers have been freed
   m_kokkos_handle.reset();
@@ -133,7 +137,7 @@ static auto estimate_max_counts(double pair_cutoff,
 
 void CellStructure::rebuild_local_properties(double const pair_cutoff) {
 #ifdef ESPRESSO_CALIPER
-  CALI_CXX_MARK_FUNCTION;
+  ESPRESSO_CALI_MARK_FUNCTION;
 #endif
   assert(m_kokkos_handle);
   auto const num_part = get_unique_particles().size();
@@ -148,17 +152,30 @@ void CellStructure::rebuild_local_properties(double const pair_cutoff) {
   }
 #endif
   if (m_local_force) { // local properties are reallocated
-    Kokkos::realloc(get_local_force(), num_part);
-    // underlying View extent changed -> scratch buffers must be rebuilt
-    m_scatter_force.emplace(
-        Kokkos::Experimental::create_scatter_view(get_local_force()));
+    if (get_local_force().extent(0) == num_part) {
+      // Extents unchanged (always the case with a single MPI rank): zero the
+      // existing buffers in place instead of freeing and reallocating the
+      // O(n_threads * N) ScatterView scratch on every Verlet rebuild.
+      reset_local_force_buffers();
+      reset_torque_replicas_if_dirty();
+    } else {
+      Kokkos::realloc(get_local_force(), num_part);
+      // underlying View extent changed -> scratch buffers must be rebuilt
+      m_scatter_force.emplace(
+          Kokkos::Experimental::create_scatter_view(get_local_force()));
 #ifdef ESPRESSO_ROTATION
-    Kokkos::realloc(get_local_torque(), num_part);
-    // underlying View extent changed -> scratch buffers must be rebuilt
-    m_scatter_torque.emplace(
-        Kokkos::Experimental::create_scatter_view(get_local_torque()));
+      Kokkos::realloc(get_local_torque(), num_part);
+      // underlying View extent changed -> scratch buffers must be rebuilt
+      m_scatter_torque.emplace(
+          Kokkos::Experimental::create_scatter_view(get_local_torque()));
+      m_torque_replicas_dirty = false;
 #endif
-    Kokkos::realloc(get_id_to_index(), get_cached_max_local_particle_id() + 1);
+    }
+    if (get_id_to_index().extent(0) !=
+        static_cast<std::size_t>(get_cached_max_local_particle_id() + 1)) {
+      Kokkos::realloc(Kokkos::WithoutInitializing, get_id_to_index(),
+                      get_cached_max_local_particle_id() + 1);
+    }
     Kokkos::deep_copy(get_id_to_index(), -1);
     // Resize particle views using AoSoA_pack's resize method
     m_aosoa->resize(num_part);
@@ -187,35 +204,48 @@ void CellStructure::rebuild_local_properties(double const pair_cutoff) {
         std::make_unique<ListType>(0ul, num_part, max_counts);
   }
 #ifdef ESPRESSO_NPT
-  m_local_virial = std::make_unique<VirialType>("local_virial");
-  m_scatter_virial.emplace(
-      Kokkos::Experimental::create_scatter_view(*m_local_virial));
+  if (not m_local_virial) {
+    m_local_virial = std::make_unique<VirialType>("local_virial");
+    m_scatter_virial.emplace(
+        Kokkos::Experimental::create_scatter_view(*m_local_virial));
+  } else {
+    reset_virial_replicas_if_dirty();
+  }
 #endif
 }
 
-void CellStructure::reset_local_force_and_torque() {
-#ifdef ESPRESSO_CALIPER
-  CALI_CXX_MARK_FUNCTION;
-#endif
+void CellStructure::reset_local_force_buffers() {
   Kokkos::deep_copy(get_local_force(), 0.);
   m_scatter_force->reset();
-#ifdef ESPRESSO_ROTATION
-  Kokkos::deep_copy(get_local_torque(), 0.);
-  m_scatter_torque->reset();
-#endif
 }
 
-void CellStructure::reset_local_properties() {
-  Kokkos::deep_copy(get_local_force(), 0.);
-  m_scatter_force->reset();
+void CellStructure::reset_torque_replicas_if_dirty() {
 #ifdef ESPRESSO_ROTATION
-  Kokkos::deep_copy(get_local_torque(), 0.);
-  m_scatter_torque->reset();
+  if (m_torque_replicas_dirty) {
+    Kokkos::deep_copy(get_local_torque(), 0.);
+    m_scatter_torque->reset();
+    m_torque_replicas_dirty = false;
+  }
 #endif
+}
+
+void CellStructure::reset_virial_replicas_if_dirty() {
 #ifdef ESPRESSO_NPT
-  Kokkos::deep_copy(get_local_virial(), 0.);
-  m_scatter_virial->reset();
+  if (m_virial_replicas_dirty) {
+    Kokkos::deep_copy(get_local_virial(), 0.);
+    m_scatter_virial->reset();
+    m_virial_replicas_dirty = false;
+  }
+#endif
+}
+
+void CellStructure::reset_local_properties() {
+#ifdef ESPRESSO_CALIPER
+  ESPRESSO_CALI_MARK_FUNCTION;
 #endif
+  reset_local_force_buffers();
+  reset_torque_replicas_if_dirty();
+  reset_virial_replicas_if_dirty();
   Kokkos::deep_copy(get_aosoa().flags, uint8_t{0});
 }
 
@@ -259,7 +289,7 @@ void CellStructure::update_bond_storage(int &pair_count, int &angle_count,
 
 void CellStructure::set_index_map() {
 #ifdef ESPRESSO_CALIPER
-  CALI_CXX_MARK_FUNCTION;
+  ESPRESSO_CALI_MARK_FUNCTION;
 #endif
   auto &unique_particles = m_unique_particles;
   unique_particles.clear();
@@ -267,42 +297,50 @@ void CellStructure::set_index_map() {
   std::unordered_set<int> registered_index{};
   using execution_space = Kokkos::DefaultHostExecutionSpace;
   int n_threads = execution_space().concurrency();
-  std::vector<int> max_ids(n_threads);
 
   m_bond_state->reset_counts();
-  std::vector<int> pair_counts(n_threads, 0);
-  std::vector<int> angle_counts(n_threads, 0);
-  std::vector<int> dihedral_counts(n_threads, 0);
-
-  enumerate_local_particles(
-      *this, [&unique_particles, &max_ids, &pair_counts, &angle_counts,
-              &dihedral_counts](std::size_t index, Particle &p) {
-        unique_particles[index] = &p;
-        auto const thread_num = omp_get_thread_num();
-        max_ids[thread_num] = std::max(p.id(), max_ids[thread_num]);
-        for (auto const bond : p.bonds()) {
-          if (not bond.partner_ids().empty()) {
-            auto const partner_ids = bond.partner_ids();
-            if (partner_ids.size() == 1u) {
-              pair_counts[thread_num] += 1;
-            } else if (partner_ids.size() == 2u) {
-              angle_counts[thread_num] += 1;
-            } else if (partner_ids.size() == 3u) {
-              dihedral_counts[thread_num] += 1;
-            }
-          }
+  // one cache line per thread: these counters are written on every particle,
+  // so packing them into shared cache lines makes the sweep bounce lines
+  // between L3 domains
+  struct alignas(64) PerThreadCounts {
+    int max_id = 0;
+    int pair = 0;
+    int angle = 0;
+    int dihedral = 0;
+  };
+  std::vector<PerThreadCounts> thread_counts(n_threads);
+
+  enumerate_local_particles(*this, [&unique_particles, &thread_counts](
+                                       std::size_t index, Particle &p) {
+    unique_particles[index] = &p;
+    auto &counts = thread_counts[omp_get_thread_num()];
+    counts.max_id = std::max(p.id(), counts.max_id);
+    for (auto const bond : p.bonds()) {
+      if (not bond.partner_ids().empty()) {
+        auto const partner_ids = bond.partner_ids();
+        if (partner_ids.size() == 1u) {
+          counts.pair += 1;
+        } else if (partner_ids.size() == 2u) {
+          counts.angle += 1;
+        } else if (partner_ids.size() == 3u) {
+          counts.dihedral += 1;
         }
-      });
+      }
+    }
+  });
   Kokkos::fence();
-  int pair_count = std::reduce(std::begin(pair_counts), std::end(pair_counts));
-  int angle_count =
-      std::reduce(std::begin(angle_counts), std::end(angle_counts));
-  int dihedral_count =
-      std::reduce(std::begin(dihedral_counts), std::end(dihedral_counts));
+  int pair_count = 0;
+  int angle_count = 0;
+  int dihedral_count = 0;
+  int max_id = 0;
+  for (auto const &counts : thread_counts) {
+    pair_count += counts.pair;
+    angle_count += counts.angle;
+    dihedral_count += counts.dihedral;
+    max_id = std::max(counts.max_id, max_id);
+  }
   set_local_bond_numbers(pair_count, angle_count, dihedral_count);
   m_bond_state->allocate();
-
-  int max_id = *(std::max_element(max_ids.begin(), max_ids.end()));
   for (auto &p : ghost_particles()) {
     auto const *local_particle = get_local_particle(p.id());
     if (not local_particle) {
@@ -472,26 +510,89 @@ unsigned map_data_parts(unsigned data_parts) {
 #ifdef ESPRESSO_BOND_CONSTRAINT
          | ((data_parts & DATA_PART_RATTLE) ? GHOSTTRANS_RATTLE : 0u)
 #endif
-         | ((data_parts & DATA_PART_BONDS) ? GHOSTTRANS_BONDS : 0u);
+         | ((data_parts & DATA_PART_BONDS) ? GHOSTTRANS_BONDS : 0u)
+#ifdef ESPRESSO_ROTATION
+         | ((data_parts & DATA_PART_QUAT) ? GHOSTTRANS_QUAT : 0u)
+         | ((data_parts & DATA_PART_TORQUE) ? GHOSTTRANS_TORQUE : 0u)
+#endif
+         ;
   /* clang-format on */
 }
 
 void CellStructure::ghosts_count() {
-  ghost_communicator(decomposition().exchange_ghosts_comm(),
-                     *get_system().box_geo, GHOSTTRANS_PARTNUM);
+#ifdef ESPRESSO_CALIPER
+  ESPRESSO_CALI_MARK_FUNCTION;
+#endif
+  GhostComm::halo_exchange(
+      *decomposition().halo_plan(), *get_system().box_geo, GHOSTTRANS_PARTNUM,
+      {GhostComm::Direction::Push, GhostComm::Combine::Overwrite},
+      m_ghost_buffers);
 }
 void CellStructure::ghosts_update(unsigned data_parts) {
-  ghost_communicator(decomposition().exchange_ghosts_comm(),
-                     *get_system().box_geo, map_data_parts(data_parts));
+#ifdef ESPRESSO_CALIPER
+  ESPRESSO_CALI_MARK_FUNCTION;
+#endif
+  auto const parts = map_data_parts(data_parts);
+  GhostComm::halo_exchange(
+      *decomposition().halo_plan(), *get_system().box_geo, parts,
+      {GhostComm::Direction::Push, GhostComm::Combine::Overwrite},
+      m_ghost_buffers);
 }
 void CellStructure::ghosts_reduce_forces() {
-  ghost_communicator(decomposition().collect_ghost_force_comm(),
-                     *get_system().box_geo, GHOSTTRANS_FORCE);
+#ifdef ESPRESSO_CALIPER
+  ESPRESSO_CALI_MARK_FUNCTION;
+#endif
+  GhostComm::halo_exchange(
+      *decomposition().halo_plan(), *get_system().box_geo,
+      get_system().get_force_reduce_ghost_flags(),
+      {GhostComm::Direction::Reduce, GhostComm::Combine::Add}, m_ghost_buffers);
+}
+
+void CellStructure::ghosts_reduce_forces_start() {
+  assert(not m_pending_ghost_reduce.has_value() &&
+         "ghosts_reduce_forces_start: a reduction is already in flight");
+  // BEGIN fires only after emplace succeeds so a throwing start cannot
+  // leave the Caliper region open without a matching END.
+  m_pending_ghost_reduce.emplace(GhostComm::halo_exchange_start(
+      *decomposition().halo_plan(), *get_system().box_geo,
+      get_system().get_force_reduce_ghost_flags(),
+      {GhostComm::Direction::Reduce, GhostComm::Combine::Add},
+      m_ghost_buffers));
+#ifdef ESPRESSO_CALIPER
+  ESPRESSO_CALI_MARK_BEGIN("ghosts_reduce_forces");
+#endif
+}
+
+void CellStructure::ghosts_reduce_forces_finish() {
+  assert(m_pending_ghost_reduce.has_value() &&
+         "ghosts_reduce_forces_finish: no reduction is in flight");
+  try {
+    GhostComm::halo_exchange_finish(*m_pending_ghost_reduce);
+  } catch (...) {
+    // A failed finish cannot be retried: the exchange state is half-consumed
+    // (some requests waited, some buffers unpacked). Drop the pending state so
+    // a later finish attempt (e.g. the ReduceGuard in integrate.cpp) does not
+    // re-run MPI waits on completed requests.
+    m_pending_ghost_reduce.reset();
+#ifdef ESPRESSO_CALIPER
+    ESPRESSO_CALI_MARK_END("ghosts_reduce_forces");
+#endif
+    throw;
+  }
+#ifdef ESPRESSO_CALIPER
+  // END before reset so the region is closed before the optional is cleared.
+  ESPRESSO_CALI_MARK_END("ghosts_reduce_forces");
+#endif
+  m_pending_ghost_reduce.reset();
 }
 #ifdef ESPRESSO_BOND_CONSTRAINT
 void CellStructure::ghosts_reduce_rattle_correction() {
-  ghost_communicator(decomposition().collect_ghost_force_comm(),
-                     *get_system().box_geo, GHOSTTRANS_RATTLE);
+#ifdef ESPRESSO_CALIPER
+  ESPRESSO_CALI_MARK_FUNCTION;
+#endif
+  GhostComm::halo_exchange(
+      *decomposition().halo_plan(), *get_system().box_geo, GHOSTTRANS_RATTLE,
+      {GhostComm::Direction::Reduce, GhostComm::Combine::Add}, m_ghost_buffers);
 }
 #endif
 
@@ -510,6 +611,12 @@ struct UpdateParticleIndexVisitor {
 } // namespace
 
 void CellStructure::resort_particles(bool global_flag) {
+#ifdef ESPRESSO_CALIPER
+  ESPRESSO_CALI_MARK_FUNCTION;
+#endif
+  assert(not m_pending_ghost_reduce.has_value() &&
+         "resort_particles: ghost force reduction is still in flight — "
+         "call ghosts_reduce_forces_finish() first");
   invalidate_ghosts();
 
   std::vector<ParticleChange> diff;
@@ -591,6 +698,9 @@ void CellStructure::set_verlet_skin_heuristic() {
 }
 
 void CellStructure::update_ghosts_and_resort_particle(unsigned data_parts) {
+#ifdef ESPRESSO_CALIPER
+  ESPRESSO_CALI_MARK_FUNCTION;
+#endif
   /* data parts that are only updated on resort */
   auto constexpr resort_only_parts =
       Cells::DATA_PART_PROPERTIES | Cells::DATA_PART_BONDS;
diff --git a/src/core/cell_system/CellStructure.hpp b/src/core/cell_system/CellStructure.hpp
index 83b3ed9706..416b383009 100644
--- a/src/core/cell_system/CellStructure.hpp
+++ b/src/core/cell_system/CellStructure.hpp
@@ -35,6 +35,7 @@
 #include "config/config.hpp"
 #include "custom_verlet_list.hpp"
 #include "ghosts.hpp"
+#include "ghosts/HaloExchange.hpp"
 #include "system/Leaf.hpp"
 
 #include <utils/Vector.hpp>
@@ -94,7 +95,11 @@ enum DataPart : unsigned {
 #ifdef ESPRESSO_BOND_CONSTRAINT
   DATA_PART_RATTLE = 32u, /**< Particle::rattle */
 #endif
-  DATA_PART_BONDS = 64u /**< Particle::bonds */
+  DATA_PART_BONDS = 64u, /**< Particle::bonds */
+#ifdef ESPRESSO_ROTATION
+  DATA_PART_QUAT = 128u,   /**< orientation quaternion (pushed with position) */
+  DATA_PART_TORQUE = 256u, /**< torque (reduced with force) */
+#endif
 };
 } // namespace Cells
 
@@ -200,13 +205,53 @@ class CellStructure : public System::Leaf<CellStructure> {
 #ifdef ESPRESSO_ROTATION
   std::unique_ptr<ForceType> m_local_torque;
   std::optional<ScatterForce> m_scatter_torque;
+  /** @brief True if a kernel may have written to @ref m_scatter_torque since
+   *  the last reset. Cleared by the resets; when false, the torque buffers
+   *  are known to be all-zero and their O(n_threads * N) zeroing and
+   *  reduction can be skipped.
+   */
+  bool m_torque_replicas_dirty = false;
 #endif
 #ifdef ESPRESSO_NPT
   std::unique_ptr<VirialType> m_local_virial;
   std::optional<ScatterVirial> m_scatter_virial;
+  /** @brief Same contract as @ref m_torque_replicas_dirty, for the virial. */
+  bool m_virial_replicas_dirty = false;
 #endif
   std::unique_ptr<LocalBondState> m_bond_state;
   std::unique_ptr<ListType> m_verlet_list_cabana;
+  /**
+   * @brief Persistent per-neighbor buffer pool for ghost exchanges.
+   *
+   * Reused across calls to @ref ghosts_count / @ref ghosts_update /
+   * @ref ghosts_reduce_forces so that, after the first (warm-up) exchange,
+   * the underlying @c std::vector storage is retained and no per-step heap
+   * allocation occurs on the hot ghost-communication path.
+   *
+   * Mutable because the ghost methods are logically const with respect to
+   * particle data but mutate this scratch storage.
+   */
+  mutable GhostComm::ExchangeBuffers m_ghost_buffers;
+  /**
+   * @brief Scratch cell-pointer list for filtered particle iteration.
+   *
+   * Used by @c for_each_interior_particle and @c for_each_boundary_particle
+   * to hold the filtered subset of local cells.  Declared mutable so that the
+   * const-qualified iteration helpers can write to it; not thread-safe — the
+   * two filtered passes must not run concurrently (they never do: interior pass
+   * completes before the ghost-reduce finish, which precedes the boundary
+   * pass).
+   */
+  mutable std::vector<Cell *> m_filtered_cells_scratch;
+  /**
+   * @brief In-flight ghost force reduction started by
+   *        @ref ghosts_reduce_forces_start.
+   *
+   * When set, @ref ghosts_reduce_forces_finish must be called before any
+   * resort or decomposition change.  Stored as optional so that the
+   * unfinished state is detectable at runtime.
+   */
+  mutable std::optional<GhostComm::GhostExchange> m_pending_ghost_reduce;
   /** particle properties using individual Kokkos Views */
   std::unique_ptr<AoSoA_pack> m_aosoa;
   /** The local id-to-index for aosoa data */
@@ -365,6 +410,68 @@ class CellStructure : public System::Leaf<CellStructure> {
     }
   }
 
+  /**
+   * @brief Run a kernel on interior (non-boundary) local particles only.
+   *
+   * A cell is interior iff none of its neighbors is a ghost cell
+   * (@ref Cell::is_boundary returns false).  This filtered variant uses the
+   * same Kokkos parallelization as @ref for_each_local_particle — it builds
+   * a filtered cell list and delegates to @c parallel_for_each_particle_impl
+   * so that thread-level behavior is identical to the unfiltered path.
+   *
+   * The kernel is assumed to be thread-safe.
+   */
+  template <typename Callable>
+  void for_each_interior_particle(Callable &&f) const {
+    auto const all_cells = decomposition().local_cells();
+    // Build a filtered list of interior cells (those that are not boundary).
+    m_filtered_cells_scratch.clear();
+    for (auto *c : all_cells) {
+      if (not c->is_boundary())
+        m_filtered_cells_scratch.push_back(c);
+    }
+    if (m_filtered_cells_scratch.empty())
+      return;
+    std::span<Cell *const> span{m_filtered_cells_scratch};
+    if (use_parallel_for_each_local_particle()) {
+      parallel_for_each_particle_impl(span, f);
+    } else {
+      for (auto *c : span)
+        for (auto &p : c->particles())
+          f(p);
+    }
+  }
+
+  /**
+   * @brief Run a kernel on boundary local particles only.
+   *
+   * Complement of @c for_each_interior_particle: visits particles in cells
+   * where @ref Cell::is_boundary returns true.  Uses the same Kokkos
+   * parallelization structure as @ref for_each_local_particle.
+   *
+   * The kernel is assumed to be thread-safe.
+   */
+  template <typename Callable>
+  void for_each_boundary_particle(Callable &&f) const {
+    auto const all_cells = decomposition().local_cells();
+    // Build a filtered list of boundary cells.
+    m_filtered_cells_scratch.clear();
+    for (auto *c : all_cells) {
+      if (c->is_boundary())
+        m_filtered_cells_scratch.push_back(c);
+    }
+    if (m_filtered_cells_scratch.empty())
+      return;
+    std::span<Cell *const> span{m_filtered_cells_scratch};
+    if (use_parallel_for_each_local_particle()) {
+      parallel_for_each_particle_impl(span, f);
+    } else {
+      for (auto *c : span)
+        for (auto &p : c->particles())
+          f(p);
+    }
+  }
+
   /**
    * @brief Run a kernel on all ghost particles.
    * The kernel is assumed to be thread-safe.
@@ -548,6 +655,34 @@ class CellStructure : public System::Leaf<CellStructure> {
    */
   void ghosts_reduce_forces();
 
+  /**
+   * @brief Begin the split-phase ghost force reduction (non-blocking).
+   *
+   * Posts all MPI sends/receives for the force reduction and returns
+   * immediately.  The caller must later call @ref ghosts_reduce_forces_finish.
+   * Asserts that no reduction is already in flight.
+   *
+   * Only meaningful when @c comm_cart.size() > 1.  At one rank the blocking
+   * @ref ghosts_reduce_forces must be used instead (the reduction is then a
+   * pure local copy with nothing to hide).
+   */
+  void ghosts_reduce_forces_start();
+
+  /**
+   * @brief Complete the split-phase ghost force reduction.
+   *
+   * Waits for all outstanding MPI requests and unpacks/reduces force data
+   * into real particles.  Clears the pending state afterwards.
+   *
+   * Must be called exactly once after @ref ghosts_reduce_forces_start.
+   */
+  void ghosts_reduce_forces_finish();
+
+  /** @brief True when a split-phase force reduction is in flight. */
+  bool has_pending_ghost_reduce() const {
+    return m_pending_ghost_reduce.has_value();
+  }
+
   /** Set forces and torques on all ghosts to zero. */
   void ghosts_reset_forces() {
     for_each_ghost_particle([](Particle &p) { p.force_and_torque() = {}; });
@@ -654,6 +789,9 @@ class CellStructure : public System::Leaf<CellStructure> {
   /** @brief Set the particle decomposition, keeping the particles. */
   void set_particle_decomposition(
       std::unique_ptr<ParticleDecomposition> &&decomposition) {
+    assert(not m_pending_ghost_reduce.has_value() &&
+           "set_particle_decomposition: ghost force reduction is still in "
+           "flight — call ghosts_reduce_forces_finish() first");
     clear_particle_index();
 
     /* Swap in new cell system */
@@ -727,7 +865,8 @@ class CellStructure : public System::Leaf<CellStructure> {
   void set_kokkos_handle(std::shared_ptr<KokkosHandle> handle);
   void rebuild_local_properties(double pair_cutoff);
   void reset_local_properties();
-  void reset_local_force_and_torque();
+  /** @brief Zero the local force view and its scatter replicas. */
+  void reset_local_force_buffers();
 
   auto &get_id_to_index() { return *m_id_to_index; }
   auto &get_local_force() { return *m_local_force; }
@@ -735,10 +874,17 @@ class CellStructure : public System::Leaf<CellStructure> {
 #ifdef ESPRESSO_ROTATION
   auto &get_local_torque() { return *m_local_torque; }
   auto get_scatter_torque() { return *m_scatter_torque; }
+  /** @brief Declare that a kernel scattering into the torque view is about
+   *  to run. Must be called before @ref reduce via the force loop dispatch.
+   */
+  void mark_torque_replicas_dirty() { m_torque_replicas_dirty = true; }
+  auto torque_replicas_dirty() const { return m_torque_replicas_dirty; }
 #endif
 #ifdef ESPRESSO_NPT
   auto &get_local_virial() { return *m_local_virial; }
   auto get_scatter_virial() { return *m_scatter_virial; }
+  void mark_virial_replicas_dirty() { m_virial_replicas_dirty = true; }
+  auto virial_replicas_dirty() const { return m_virial_replicas_dirty; }
 #endif
 
   auto &get_aosoa() { return *m_aosoa; }
@@ -799,6 +945,13 @@ class CellStructure : public System::Leaf<CellStructure> {
   }
 
 private:
+  /** @brief Zero the torque buffers iff a kernel scattered into them since
+   *  the last reset (no-op otherwise, and without the ROTATION feature).
+   */
+  void reset_torque_replicas_if_dirty();
+  /** @brief Same contract for the virial buffers (NPT feature). */
+  void reset_virial_replicas_if_dirty();
+
   /** Non-bonded pair loop with verlet lists.
    *
    * @param pair_kernel Kernel to apply
diff --git a/src/core/cell_system/HybridDecomposition.cpp b/src/core/cell_system/HybridDecomposition.cpp
index 5eb42fab9b..0b41c670c1 100644
--- a/src/core/cell_system/HybridDecomposition.cpp
+++ b/src/core/cell_system/HybridDecomposition.cpp
@@ -27,6 +27,10 @@
 #include "BoxGeometry.hpp"
 #include "LocalBox.hpp"
 #include "ParticleList.hpp"
+#include "ghosts.hpp"
+#include "ghosts/HaloExchange.hpp"
+#include "ghosts/HaloPlanValidator.hpp"
+#include "ghosts/mark_boundary_cells.hpp"
 
 #include <utils/Vector.hpp>
 #include <utils/mpi/sendrecv.hpp>
@@ -35,6 +39,7 @@
 #include <boost/mpi/communicator.hpp>
 
 #include <algorithm>
+#include <cassert>
 #include <cstddef>
 #include <functional>
 #include <iterator>
@@ -65,20 +70,6 @@ HybridDecomposition::HybridDecomposition(boost::mpi::communicator comm,
   auto ghost_cells_n_square = m_n_square.get_ghost_cells();
   std::ranges::copy(ghost_cells_n_square, std::back_inserter(m_ghost_cells));
 
-  /* Communicators that contain communications of both child decompositions */
-  m_exchange_ghosts_comm = m_regular_decomposition.exchange_ghosts_comm();
-  auto exchange_ghosts_comm_n_square = m_n_square.exchange_ghosts_comm();
-  std::ranges::copy(exchange_ghosts_comm_n_square.communications,
-                    std::back_inserter(m_exchange_ghosts_comm.communications));
-
-  m_collect_ghost_force_comm =
-      m_regular_decomposition.collect_ghost_force_comm();
-  auto collect_ghost_force_comm_n_square =
-      m_n_square.collect_ghost_force_comm();
-  std::ranges::copy(
-      collect_ghost_force_comm_n_square.communications,
-      std::back_inserter(m_collect_ghost_force_comm.communications));
-
   /* coupling between the child decompositions via neighborship relation */
   std::vector<Cell *> additional_reds = m_n_square.get_local_cells();
   std::ranges::copy(ghost_cells_n_square, std::back_inserter(additional_reds));
@@ -90,6 +81,65 @@ HybridDecomposition::HybridDecomposition(boost::mpi::communicator comm,
     std::ranges::copy(additional_reds, std::back_inserter(red_neighbors));
     local_cell->m_neighbors = Neighbors<Cell *>(red_neighbors, black_neighbors);
   }
+
+  m_halo_plan = make_halo_plan();
+
+  /* classify local cells as interior or boundary.
+   *
+   * HybridDecomposition combines a RegularDecomposition (which already has
+   * its own wrap-aware classification) with an AtomDecomposition n-square
+   * child.  The combined cell set extends every regular local cell's
+   * neighbor list with the n-square cells, so regular cells that were
+   * interior under RegularDecomposition alone may now interact with
+   * n-square ghost cells.  Determining which cells remain genuinely
+   * interior after the coupling is non-trivial; we conservatively mark all
+   * local cells as boundary.  The compute/comm overlap degenerates
+   * gracefully to a no-op interior pass — identical to today's blocking
+   * path.
+   */
+  GhostComm::mark_boundary_cells(HybridDecomposition::local_cells(),
+                                 HybridDecomposition::ghost_cells());
+  for (Cell *c : HybridDecomposition::local_cells()) {
+    c->m_is_boundary = true;
+  }
+#ifdef ESPRESSO_ADDITIONAL_CHECKS
+  assert(GhostComm::report_violations(
+      GhostComm::validate_halo_plan(m_halo_plan,
+                                    HybridDecomposition::local_cells(),
+                                    HybridDecomposition::ghost_cells()),
+      "HybridDecomposition"));
+  // NOTE: validate_halo_plan_symmetry is NOT called here.
+  // During checkpoint loading, decompositions are transiently rebuilt while
+  // maximal_cutoff is rank-divergent (ranks may have different cell grids for a
+  // brief window before the next consistent rebuild).  The transient plan is
+  // never used — it is immediately replaced — so the asymmetry is harmless.
+  // A construction-time collective all_to_all inside a ctor is also dangerous:
+  // if one rank aborts the others block forever in the collective.
+  // Symmetry is instead validated at FIRST USE of the plan in
+  // halo_exchange_start (see GhostComm::halo_exchange_start in
+  // HaloExchange.cpp).
+#endif
+}
+
+GhostComm::HaloPlan HybridDecomposition::make_halo_plan() {
+  // Build the combined plan: p2p neighbors and local copies come from the
+  // regular child; the collective section comes from the n-square child.
+  GhostComm::HaloPlan plan;
+  plan.comm = m_comm; // use the world comm for HybridDecomposition's plan
+
+  auto const *regular_plan = m_regular_decomposition.halo_plan();
+  if (regular_plan) {
+    plan.neighbors = regular_plan->neighbors;
+    plan.local = regular_plan->local;
+  }
+
+  // Overlay the collective section from the n-square child.
+  auto const *nsq_plan = m_n_square.halo_plan();
+  if (nsq_plan && nsq_plan->collective) {
+    plan.collective = nsq_plan->collective;
+  }
+
+  return plan;
 }
 
 void HybridDecomposition::resort(bool global,
@@ -155,9 +205,12 @@ void HybridDecomposition::resort(bool global,
   m_regular_decomposition.resort(global, diff);
   m_n_square.resort(global, diff);
 
-  ghost_communicator(exchange_ghosts_comm(), m_box, GHOSTTRANS_PARTNUM);
-  ghost_communicator(exchange_ghosts_comm(), m_box,
-                     map_data_parts(m_get_global_ghost_flags()));
+  GhostComm::halo_exchange(
+      m_halo_plan, m_box, GHOSTTRANS_PARTNUM,
+      {GhostComm::Direction::Push, GhostComm::Combine::Overwrite});
+  GhostComm::halo_exchange(
+      m_halo_plan, m_box, map_data_parts(m_get_global_ghost_flags()),
+      {GhostComm::Direction::Push, GhostComm::Combine::Overwrite});
 }
 
 std::size_t HybridDecomposition::count_particles(
diff --git a/src/core/cell_system/HybridDecomposition.hpp b/src/core/cell_system/HybridDecomposition.hpp
index 51788080d4..e962c3c9c0 100644
--- a/src/core/cell_system/HybridDecomposition.hpp
+++ b/src/core/cell_system/HybridDecomposition.hpp
@@ -30,7 +30,6 @@
 #include "BoxGeometry.hpp"
 #include "LocalBox.hpp"
 #include "Particle.hpp"
-#include "ghosts.hpp"
 
 #include <utils/Vector.hpp>
 
@@ -60,8 +59,10 @@ class HybridDecomposition : public ParticleDecomposition {
   std::vector<Cell *> m_local_cells;
   std::vector<Cell *> m_ghost_cells;
 
-  GhostCommunicator m_exchange_ghosts_comm;
-  GhostCommunicator m_collect_ghost_force_comm;
+  /// Holds ParticleList* into this decomposition's cells — value-copying this
+  /// object leaves these pointers dangling. TODO(Task 1.7): make non-copyable
+  /// or rebuild-on-copy.
+  GhostComm::HaloPlan m_halo_plan;
 
   /** RegularDecomposition to hold the small particles */
   RegularDecomposition m_regular_decomposition;
@@ -92,13 +93,7 @@ class HybridDecomposition : public ParticleDecomposition {
 
   auto get_cutoff_regular() const { return m_cutoff_regular; }
 
-  GhostCommunicator const &exchange_ghosts_comm() const override {
-    return m_exchange_ghosts_comm;
-  }
-
-  GhostCommunicator const &collect_ghost_force_comm() const override {
-    return m_collect_ghost_force_comm;
-  }
+  GhostComm::HaloPlan const *halo_plan() const override { return &m_halo_plan; }
 
   std::span<Cell *const> local_cells() const override { return m_local_cells; }
   std::span<Cell *const> ghost_cells() const override { return m_ghost_cells; }
@@ -140,5 +135,11 @@ class HybridDecomposition : public ParticleDecomposition {
   }
 
 private:
+  /**
+   * @brief Build the plan-based halo plan combining the regular child's
+   *        neighbors/local with the n-square child's collective section.
+   */
+  GhostComm::HaloPlan make_halo_plan();
+
   std::size_t count_particles(std::vector<Cell *> const &local_cells) const;
 };
diff --git a/src/core/cell_system/ParticleDecomposition.hpp b/src/core/cell_system/ParticleDecomposition.hpp
index a28ea78bb0..d8eed202b9 100644
--- a/src/core/cell_system/ParticleDecomposition.hpp
+++ b/src/core/cell_system/ParticleDecomposition.hpp
@@ -22,7 +22,7 @@
 #include "cell_system/Cell.hpp"
 
 #include "BoxGeometry.hpp"
-#include "ghosts.hpp"
+#include "ghosts/HaloPlan.hpp"
 
 #include <utils/Vector.hpp>
 
@@ -74,13 +74,12 @@ class ParticleDecomposition {
   virtual void resort(bool global_flag, std::vector<ParticleChange> &diff) = 0;
 
   /**
-   * @brief Communicator for updating ghosts from the real particles.
-   */
-  virtual GhostCommunicator const &exchange_ghosts_comm() const = 0;
-  /**
-   * @brief Communicator for force reduction.
+   * @brief Topology-agnostic halo-exchange plan.
+   *
+   * Returns a @ref GhostComm::HaloPlan describing the ghost communication for
+   * this decomposition. All decompositions implement this.
    */
-  virtual GhostCommunicator const &collect_ghost_force_comm() const = 0;
+  virtual GhostComm::HaloPlan const *halo_plan() const = 0;
 
   /**
    * @brief Get pointer to local cells.
diff --git a/src/core/cell_system/RegularDecomposition.cpp b/src/core/cell_system/RegularDecomposition.cpp
index e51735b06a..18af124ad1 100644
--- a/src/core/cell_system/RegularDecomposition.cpp
+++ b/src/core/cell_system/RegularDecomposition.cpp
@@ -22,6 +22,8 @@
 #include "cell_system/RegularDecomposition.hpp"
 
 #include "cell_system/Cell.hpp"
+#include "ghosts/HaloPlanValidator.hpp"
+#include "ghosts/mark_boundary_cells.hpp"
 
 #include "communication.hpp"
 #include "error_handling/RuntimeErrorStream.hpp"
@@ -37,16 +39,25 @@
 #include <boost/mpi/collectives/all_reduce.hpp>
 #include <boost/mpi/communicator.hpp>
 #include <boost/mpi/request.hpp>
-#include <boost/range/algorithm/reverse.hpp>
 #include <boost/range/numeric.hpp>
 
+#include <Kokkos_Core.hpp>
+
 #include <algorithm>
 #include <array>
+#include <cassert>
 #include <cmath>
 #include <cstddef>
+#include <exception>
 #include <functional>
 #include <initializer_list>
 #include <iterator>
+#include <map>
+#include <mutex>
+#include <numeric>
+#include <set>
+#include <stdexcept>
+#include <string>
 #include <utility>
 #include <vector>
 
@@ -177,31 +188,124 @@ void RegularDecomposition::resort(bool global,
                                   std::vector<ParticleChange> &diff) {
   ParticleList displaced_parts;
 
-  for (auto &c : local_cells()) {
-    for (auto it = c->particles().begin(); it != c->particles().end();) {
-      fold_and_reset(*it, m_box);
+  auto const cells_span = local_cells();
+  auto const n_cells = cells_span.size();
+
+  /* Remove a misplaced particle from its cell and hand it to its target
+   * cell (or the displaced list when it left the local domain), recording
+   * the changes. Shared by the serial and the two-phase parallel sweep. */
+  auto const apply_move = [&](ParticleList &parts, Particle &&p,
+                              Cell *target_cell) {
+    diff.emplace_back(ModifiedList{parts});
+    /* Particle is not local */
+    if (target_cell == nullptr) {
+      diff.emplace_back(RemovedParticle{p.id()});
+      displaced_parts.insert(std::move(p));
+    }
+    /* Particle belongs on this node but is in the wrong cell. */
+    else {
+      target_cell->particles().insert(std::move(p));
+      diff.emplace_back(ModifiedList{target_cell->particles()});
+    }
+  };
 
-      auto target_cell = particle_to_cell(*it);
+  if (Kokkos::DefaultHostExecutionSpace().concurrency() <= 1) {
+    /* Single-threaded rank: one-pass sweep without the bookkeeping overhead
+     * of the two-phase version below. */
+    for (auto *const c : cells_span) {
+      for (auto it = c->particles().begin(); it != c->particles().end();) {
+        fold_and_reset(*it, m_box);
 
-      /* Particle is in place */
-      if (target_cell == c) {
-        std::advance(it, 1);
-        continue;
-      }
+        auto *const target_cell = particle_to_cell(*it);
 
-      auto p = std::move(*it);
-      it = c->particles().erase(it);
-      diff.emplace_back(ModifiedList{c->particles()});
+        /* Particle is in place */
+        if (target_cell == c) {
+          std::advance(it, 1);
+          continue;
+        }
 
-      /* Particle is not local */
-      if (target_cell == nullptr) {
-        diff.emplace_back(RemovedParticle{p.id()});
-        displaced_parts.insert(std::move(p));
+        auto p = std::move(*it);
+        it = c->particles().erase(it);
+        apply_move(c->particles(), std::move(p), target_cell);
       }
-      /* Particle belongs on this node but is in the wrong cell. */
-      else if (target_cell != c) {
-        target_cell->particles().insert(std::move(p));
-        diff.emplace_back(ModifiedList{target_cell->particles()});
+    }
+  } else {
+    /* Phase 1 (parallel): fold every particle position and classify it
+     * against its target cell. Only particle-local state and this cell's own
+     * move list are written, so cells can be swept concurrently.
+     * fold_and_reset() throws on image-box overflow; exceptions must not
+     * escape the parallel region, so the first error is captured and
+     * rethrown afterwards. */
+    struct Move {
+      int index;
+      Cell *target;
+    };
+    std::vector<std::vector<Move>> moves(n_cells);
+    std::mutex fold_error_mutex;
+    std::string fold_error_msg;
+    Kokkos::parallel_for(
+        "RegularDecomposition::resort::classify",
+        Kokkos::RangePolicy<Kokkos::DefaultHostExecutionSpace>(std::size_t{0},
+                                                               n_cells),
+        [&](std::size_t const ci) {
+          auto *cell = cells_span[ci];
+          int index = 0;
+          for (auto &p : cell->particles()) {
+            try {
+              fold_and_reset(p, m_box);
+            } catch (std::exception const &err) {
+              std::lock_guard<std::mutex> guard{fold_error_mutex};
+              if (fold_error_msg.empty()) {
+                fold_error_msg = err.what();
+              }
+              return;
+            }
+            if (auto *const target = particle_to_cell(p); target != cell) {
+              moves[ci].push_back({index, target});
+            }
+            ++index;
+          }
+        });
+    Kokkos::fence();
+    if (not fold_error_msg.empty()) {
+      throw std::runtime_error(fold_error_msg);
+    }
+
+    /* Phase 2 (serial): apply the moves. This replays the serial sweep
+     * exactly: ParticleList::erase() swaps the last element into the erased
+     * slot, so a slot-to-original-index map is maintained to look up the
+     * phase-1 classification of swapped-in elements. Cells without moves
+     * (the vast majority) are skipped entirely. */
+    std::vector<int> slot;
+    std::vector<Cell *> target_of;
+    for (std::size_t ci = 0; ci < n_cells; ++ci) {
+      if (moves[ci].empty()) {
+        continue;
+      }
+      auto *const c = cells_span[ci];
+      auto &parts = c->particles();
+      auto const n = static_cast<int>(parts.size());
+      target_of.assign(n, c); // target == own cell: particle is in place
+      for (auto const &move : moves[ci]) {
+        target_of[move.index] = move.target;
+      }
+      slot.resize(n);
+      std::iota(slot.begin(), slot.end(), 0);
+      int i = 0;
+      int end = n;
+      while (i < end) {
+        auto *const target_cell = target_of[slot[i]];
+
+        /* Particle is in place */
+        if (target_cell == c) {
+          ++i;
+          continue;
+        }
+
+        auto p = std::move(*(parts.begin() + i));
+        parts.erase(parts.begin() + i); // swaps the last element into slot i
+        slot[i] = slot[--end];
+        apply_move(parts, std::move(p), target_cell);
       }
     }
   }
@@ -565,146 +669,162 @@ void RegularDecomposition::init_cell_interactions() {
       }
 }
 
-namespace {
-/** Revert the order of a communicator: After calling this the
- *  communicator is working in reverted order with exchanged
- *  communication types GHOST_SEND <-> GHOST_RECV.
- */
-void revert_comm_order(GhostCommunicator &comm) {
-  /* revert order */
-  boost::reverse(comm.communications);
-
-  /* exchange SEND/RECV */
-  for (auto &c : comm.communications) {
-    if (c.type == GHOST_SEND)
-      c.type = GHOST_RECV;
-    else if (c.type == GHOST_RECV)
-      c.type = GHOST_SEND;
-    else if (c.type == GHOST_LOCL) {
-      boost::reverse(c.part_lists);
-    }
-  }
-}
-
-/** Of every two communication rounds, set the first receivers to prefetch and
- *  poststore
- */
-void assign_prefetches(GhostCommunicator &comm) {
-  for (auto it = comm.communications.begin(); it != comm.communications.end();
-       it += 2) {
-    auto next = std::next(it);
-    if (it->type == GHOST_RECV && next->type == GHOST_SEND) {
-      it->type |= GHOST_PREFETCH | GHOST_PSTSTORE;
-      next->type |= GHOST_PREFETCH | GHOST_PSTSTORE;
-    }
-  }
-}
+GhostComm::HaloPlan RegularDecomposition::make_halo_plan() {
+  using GhostComm::HaloPlan;
+  using GhostComm::LocalComm;
+  using GhostComm::NeighborComm;
+  using GhostComm::SendRegion;
 
-} // namespace
+  HaloPlan plan;
+  plan.comm = m_comm;
 
-GhostCommunicator RegularDecomposition::prepare_comm() {
-  // When running on a single MPI rank, the ghost cells are not needed,
-  // as the corresponding physical cell is available on the same MPI rank.
-  // The cell neighborships are set up such that cells across the periodic
-  // boundaries are connected as neighbors directly.
+  // Match the legacy communicator: on a single MPI rank there are no ghost
+  // cells to fill. The cell neighbourships are set up (see
+  // init_cell_interactions) so that cells across periodic boundaries are
+  // connected directly, so the plan stays empty.
   if (m_comm.size() == 1)
-    return {};
+    return plan;
 
-  int dir, lr, i, cnt, n_comm_cells[3];
-  Utils::Vector3i lc{}, hc{}, done{};
+  auto const cart_info = Utils::Mpi::cart_get<3>(m_comm);
+  auto const &node_pos = cart_info.coords;
+  auto const &node_grid = ::communicator.node_grid;
+  // Total number of MD cells along each axis across all ranks.
+  auto const global_size = hadamard_product(node_grid, cell_grid);
+  // Global (0-based) cell index of this rank's first *local* cell, i.e. of
+  // ghost-grid coordinate (1,1,1).
+  auto const global_origin = hadamard_product(node_pos, cell_grid);
+
+  // Cartesian rank owning the cell at (folded) global cell coordinate.
+  // The Cartesian communicator is always fully periodic (see
+  // Communicator::init_comm_cart), so this is well-defined for every offset.
+  auto const owner_of = [&](Utils::Vector3i const &global_cell) {
+    auto const owner_coords =
+        hadamard_division((global_cell + global_size) % global_size, cell_grid);
+    return Utils::Mpi::cart_rank<3>(m_comm, owner_coords);
+  };
 
-  auto const comm_info = Utils::Mpi::cart_get<3>(m_comm);
-  auto const node_neighbors = Utils::Mpi::cart_neighbors<3>(m_comm);
+  // Deterministic ordering key shared by both ranks of a peer pair: the linear
+  // index of the *real* cell (in the global cell grid) that a ghost mirrors.
+  auto const global_key = [&](Utils::Vector3i const &global_cell) {
+    auto const folded = (global_cell + global_size) % global_size;
+    return Utils::get_linear_index(folded, global_size);
+  };
 
-  /* calculate number of communications */
-  std::size_t num = 0;
-  for (dir = 0; dir < 3; dir++) {
-    for (lr = 0; lr < 2; lr++) {
-      /* No communication for border of non periodic direction */
-      if (comm_info.dims[dir] == 1)
-        num++;
-      else
-        num += 2;
-    }
-  }
+  // Pointer to the particle list of the cell at ghost-grid coordinate c.
+  auto const list_at = [this](Utils::Vector3i const &c) -> ParticleList * {
+    return &cells
+                .at(static_cast<std::size_t>(
+                    Utils::get_linear_index(c, ghost_cell_grid)))
+                .particles();
+  };
 
-  /* prepare communicator */
-  auto ghost_comm = GhostCommunicator{m_comm, num};
-
-  /* number of cells to communicate in a direction */
-  n_comm_cells[0] = cell_grid[1] * cell_grid[2];
-  n_comm_cells[1] = cell_grid[2] * ghost_cell_grid[0];
-  n_comm_cells[2] = ghost_cell_grid[0] * ghost_cell_grid[1];
-
-  cnt = 0;
-  /* direction loop: x, y, z */
-  for (dir = 0; dir < 3; dir++) {
-    lc[(dir + 1) % 3] = 1 - done[(dir + 1) % 3];
-    lc[(dir + 2) % 3] = 1 - done[(dir + 2) % 3];
-    hc[(dir + 1) % 3] = cell_grid[(dir + 1) % 3] + done[(dir + 1) % 3];
-    hc[(dir + 2) % 3] = cell_grid[(dir + 2) % 3] + done[(dir + 2) % 3];
-    /* lr loop: left right */
-    /* here we could in principle build in a one sided ghost
-       communication, simply by taking the lr loop only over one
-       value */
-    for (lr = 0; lr < 2; lr++) {
-      if (comm_info.dims[dir] == 1) {
-        /* just copy cells on a single node */
-        ghost_comm.communications[cnt].type = GHOST_LOCL;
-        ghost_comm.communications[cnt].node = m_comm.rank();
-
-        /* Buffer has to contain Send and Recv cells -> factor 2 */
-        ghost_comm.communications[cnt].part_lists.resize(2 * n_comm_cells[dir]);
-
-        /* fill send ghost_comm cells */
-        lc[dir] = hc[dir] = 1 + lr * (cell_grid[dir] - 1);
-
-        fill_comm_cell_lists(ghost_comm.communications[cnt].part_lists.data(),
-                             lc, hc);
-
-        /* fill recv ghost_comm cells */
-        lc[dir] = hc[dir] = 0 + (1 - lr) * (cell_grid[dir] + 1);
-
-        /* place receive cells after send cells */
-        fill_comm_cell_lists(
-            &ghost_comm.communications[cnt].part_lists[n_comm_cells[dir]], lc,
-            hc);
-
-        cnt++;
-      } else {
-        /* i: send/recv loop */
-        for (i = 0; i < 2; i++) {
-          if ((comm_info.coords[dir] + i) % 2 == 0) {
-            ghost_comm.communications[cnt].type = GHOST_SEND;
-            ghost_comm.communications[cnt].node = node_neighbors[2 * dir + lr];
-            ghost_comm.communications[cnt].part_lists.resize(n_comm_cells[dir]);
-            /* prepare folding of ghost positions */
-
-            lc[dir] = hc[dir] = 1 + lr * (cell_grid[dir] - 1);
-
-            fill_comm_cell_lists(
-                ghost_comm.communications[cnt].part_lists.data(), lc, hc);
-            cnt++;
+  // Per-peer accumulators. A "recv" pair maps one of our ghost cells to the
+  // global index of the real cell (on the peer) it mirrors. A "send" pair maps
+  // one of our real cells to its own global index (the peer will receive it
+  // into a ghost). Sorting both lists by their key makes recv[k] line up with
+  // peer.send[k] without exchanging any index arrays.
+  struct PeerBucket {
+    std::vector<std::pair<int, ParticleList *>> recv; // (key, our ghost)
+    std::vector<std::pair<int, ParticleList *>> send; // (key, our real cell)
+  };
+  std::map<int, PeerBucket> peers;
+  std::vector<std::pair<int, LocalComm>> local; // (key, self-ghost copy)
+
+  auto const this_rank = m_comm.rank();
+  auto const one = Utils::Vector3i{1, 1, 1};
+
+  // Enumerate every ghost cell exactly once (any ghost-grid coordinate with a
+  // component in the halo, i.e. == 0 or == cell_grid+1). This guarantees each
+  // ghost is a recv/dst target exactly once, regardless of how small the node
+  // grid is (dims of 1 or 2 collapse several stencil directions onto the same
+  // peer, so a local-cell x offset enumeration would double-count).
+  //
+  // For each ghost we record a *matched pair*:
+  //   * recv: this ghost, keyed by the global index of the real cell (on the
+  //     peer) it mirrors -- so it lines up with the peer's send of that cell.
+  //   * send: our own boundary real cell that the peer mirrors as its ghost in
+  //     the opposite direction, keyed by that real cell's own global index --
+  //     so it lines up with the peer's recv. The recv<->send pairing is a
+  //     bijection, which keeps send.size() == recv.size() per peer.
+  for (int gz = 0; gz < ghost_cell_grid[2]; ++gz) {
+    for (int gy = 0; gy < ghost_cell_grid[1]; ++gy) {
+      for (int gx = 0; gx < ghost_cell_grid[0]; ++gx) {
+        Utils::Vector3i const nc{gx, gy, gz};
+        // Direction sign of the halo crossing (0 if interior in a dim).
+        Utils::Vector3i side{};
+        bool is_ghost = false;
+        for (auto d = 0u; d < 3u; ++d) {
+          if (nc[d] == 0) {
+            side[d] = -1;
+            is_ghost = true;
+          } else if (nc[d] == cell_grid[d] + 1) {
+            side[d] = +1;
+            is_ghost = true;
           }
-          if ((comm_info.coords[dir] + (1 - i)) % 2 == 0) {
-            ghost_comm.communications[cnt].type = GHOST_RECV;
-            ghost_comm.communications[cnt].node =
-                node_neighbors[2 * dir + (1 - lr)];
-            ghost_comm.communications[cnt].part_lists.resize(n_comm_cells[dir]);
-
-            lc[dir] = hc[dir] = (1 - lr) * (cell_grid[dir] + 1);
+        }
+        if (not is_ghost)
+          continue; // interior (local) cell
+
+        // Global cell mirrored by this ghost and its owning peer.
+        auto const ghost_global = global_origin + (nc - one);
+        auto const peer = owner_of(ghost_global);
+        auto const recv_key = global_key(ghost_global);
+
+        if (peer == this_rank) {
+          // Periodic self-ghost (a node-grid dim equals 1): copy the matching
+          // local cell straight into the ghost (replaces GHOST_LOCL).
+          auto const src_coord = ((ghost_global + global_size) % global_size) -
+                                 global_origin + one;
+          local.emplace_back(recv_key,
+                             LocalComm{list_at(src_coord), list_at(nc), {}});
+          continue;
+        }
 
-            fill_comm_cell_lists(
-                ghost_comm.communications[cnt].part_lists.data(), lc, hc);
-            cnt++;
-          }
+        peers[peer].recv.emplace_back(recv_key, list_at(nc));
+
+        // Dual send cell: our boundary real cell mirrored by the peer's ghost
+        // in the opposite direction. Snap crossing dims to the near boundary,
+        // keep tangential dims aligned with the ghost.
+        auto mc = nc;
+        for (auto d = 0u; d < 3u; ++d) {
+          if (side[d] == -1)
+            mc[d] = 1;
+          else if (side[d] == +1)
+            mc[d] = cell_grid[d];
         }
+        auto const send_global = global_origin + (mc - one);
+        peers[peer].send.emplace_back(global_key(send_global), list_at(mc));
       }
-      done[dir] = 1;
     }
   }
 
-  return ghost_comm;
+  // Emit one NeighborComm per peer, with send/recv sorted by their shared key.
+  auto const by_key = [](auto const &a, auto const &b) {
+    return a.first < b.first;
+  };
+  for (auto &[peer, bucket] : peers) {
+    std::sort(bucket.recv.begin(), bucket.recv.end(), by_key);
+    std::sort(bucket.send.begin(), bucket.send.end(), by_key);
+    NeighborComm nc;
+    nc.peer = peer;
+    nc.recv.reserve(bucket.recv.size());
+    for (auto const &[key, cell] : bucket.recv)
+      nc.recv.push_back(cell);
+    nc.send.reserve(bucket.send.size());
+    for (auto const &[key, cell] : bucket.send)
+      nc.send.push_back(SendRegion{cell, {}});
+    plan.neighbors.push_back(std::move(nc));
+  }
+
+  // Sort the self-copies deterministically too (not required, but keeps the
+  // plan reproducible run to run).
+  std::sort(local.begin(), local.end(),
+            [](auto const &a, auto const &b) { return a.first < b.first; });
+  plan.local.reserve(local.size());
+  for (auto &[key, lc] : local)
+    plan.local.push_back(lc);
+
+  return plan;
 }
 
 RegularDecomposition::RegularDecomposition(
@@ -723,13 +843,112 @@ RegularDecomposition::RegularDecomposition(
   /* mark local and ghost cells */
   mark_cells();
 
-  /* create communicators */
-  m_exchange_ghosts_comm = prepare_comm();
-  m_collect_ghost_force_comm = prepare_comm();
-
-  /* collect forces has to be done in reverted order! */
-  revert_comm_order(m_collect_ghost_force_comm);
+  /* build the topology-agnostic direct-neighbor halo plan */
+  m_halo_plan = make_halo_plan();
+
+  /* Classify local cells as interior or boundary.
+   *
+   * Rule (a): any neighbor that is a ghost cell → boundary [base rule].
+   * Rule (b): any neighbor relation that crosses a periodic box boundary
+   *   must also make the cell boundary, even when both cells are local.
+   *   This happens on a single MPI rank (node_grid[i]==1, periodic[i]):
+   *   init_cell_interactions() wires the first and last local layer along
+   *   axis i directly without going through ghost cells, so the base rule
+   *   misses those periodic wrap-around neighbours.
+   *
+   * Implementation: precise pair-predicate that fires iff the two cells
+   * sit in the "first layer ↔ last layer" pair along a wrap axis.
+   * Ghost-grid coordinate of a cell: unpack the column-major linear index
+   *   idx = a + G[0]*(b + G[1]*c)  (G = ghost_cell_grid).
+   * First local layer along axis i: ghost coord == 1.
+   * Last local layer along axis i: ghost coord == cell_grid[i].
+   */
+  auto const &node_grid = ::communicator.node_grid;
+  auto const idx_of = [this](Cell const *c) {
+    return static_cast<int>(c - cells.data());
+  };
+  auto const ghost_coord_of = [this](int idx) -> Utils::Vector3i {
+    int const a = idx % ghost_cell_grid[0];
+    int const bc = idx / ghost_cell_grid[0];
+    int const b = bc % ghost_cell_grid[1];
+    int const c = bc / ghost_cell_grid[1];
+    return {a, b, c};
+  };
+  // Which axes wrap locally (the local domain spans the whole box along a
+  // periodic axis)?  Precomputed by value so the predicate lambda does not
+  // capture node_grid (AppleClang rejects that non-odr-use capture with
+  // -Werror,-Wunused-lambda-capture).
+  std::array<bool, 3> wrap_axis;
+  for (int i = 0; i < 3; ++i)
+    wrap_axis[i] = (node_grid[i] == 1) && m_box.periodic(i);
+  // Build the predicate only when there is at least one wrap axis; otherwise
+  // pass nullptr (no overhead in the inner loop of mark_boundary_cells).
+  bool const has_wrap_axis = wrap_axis[0] || wrap_axis[1] || wrap_axis[2];
+  std::function<bool(Cell const *, Cell const *)> wrap_pred;
+  if (has_wrap_axis) {
+    wrap_pred = [this, wrap_axis, idx_of, ghost_coord_of](
+                    Cell const *a_cell, Cell const *b_cell) -> bool {
+      auto const a_coord = ghost_coord_of(idx_of(a_cell));
+      auto const b_coord = ghost_coord_of(idx_of(b_cell));
+      for (int i = 0; i < 3; ++i) {
+        if (wrap_axis[i]) {
+          bool const a_first = (a_coord[i] == 1);
+          bool const a_last = (a_coord[i] == cell_grid[i]);
+          bool const b_first = (b_coord[i] == 1);
+          bool const b_last = (b_coord[i] == cell_grid[i]);
+          if ((a_first && b_last) || (a_last && b_first))
+            return true;
+        }
+      }
+      return false;
+    };
+  }
+  GhostComm::mark_boundary_cells(local_cells(), ghost_cells(), wrap_pred);
+
+  /* Degenerate case: node_grid[i]==1, periodic[i], cell_grid[i]==1.
+   *
+   * When cell_grid[i]==1 on a wrap axis, the single local cell layer along
+   * that axis has itself as the only periodic neighbour (both ends fold to
+   * the same global index).  init_cell_interactions() therefore excludes
+   * the self-pair (ind1==ind2 guard at line ~552), leaving neighbors().all()
+   * empty along that axis.  The wrap_predicate above is never called for
+   * that cell, so mark_boundary_cells() leaves it interior — wrong.
+   *
+   * Correct interpretation: the cell interacts with itself across the
+   * periodic boundary, so it is by definition wrap-adjacent and must be
+   * boundary.  The simplest safe fix: if any wrap axis has cell_grid[i]==1,
+   * every local cell is boundary (they all span that axis, so all are
+   * wrap-adjacent).
+   */
+  for (int i = 0; i < 3; ++i) {
+    if (wrap_axis[i] && cell_grid[i] == 1) {
+      for (Cell *c : local_cells())
+        c->m_is_boundary = true;
+      break; // one degenerate axis is enough to force all-boundary
+    }
+  }
 
-  assign_prefetches(m_exchange_ghosts_comm);
-  assign_prefetches(m_collect_ghost_force_comm);
+  /* Plan-membership pass (source 2): mark every local cell that the plan
+   * exports as a send source.  The geometric rules above (source 1) miss
+   * plan shapes such as Lees-Edwards fully-connected boundaries and ELC
+   * periodicity-change paths, where boundary cells are determined by the
+   * plan topology rather than ghost-cell adjacency.  Both sources are
+   * complementary and must both be applied. */
+  GhostComm::mark_plan_cells_boundary(m_halo_plan, local_cells());
+
+#ifdef ESPRESSO_ADDITIONAL_CHECKS
+  assert(GhostComm::report_violations(
+      GhostComm::validate_halo_plan(m_halo_plan, local_cells(), ghost_cells()),
+      "RegularDecomposition"));
+  // NOTE: validate_halo_plan_symmetry is NOT called here.
+  // During checkpoint loading, decompositions are transiently rebuilt while
+  // maximal_cutoff is rank-divergent (ranks may have different cell grids for a
+  // brief window before the next consistent rebuild).  The transient plan is
+  // never used — it is immediately replaced — so the asymmetry is harmless.
+  // A construction-time collective all_to_all inside a ctor is also dangerous:
+  // if one rank aborts the others block forever in the collective.
+  // Symmetry is instead validated at FIRST USE of the plan in
+  // halo_exchange_start (see GhostComm::halo_exchange_start in
+  // HaloExchange.cpp).
+#endif
 }
diff --git a/src/core/cell_system/RegularDecomposition.hpp b/src/core/cell_system/RegularDecomposition.hpp
index db2dd01fa5..0ecd312896 100644
--- a/src/core/cell_system/RegularDecomposition.hpp
+++ b/src/core/cell_system/RegularDecomposition.hpp
@@ -29,7 +29,7 @@
 #include "LocalBox.hpp"
 #include "Particle.hpp"
 #include "ParticleList.hpp"
-#include "ghosts.hpp"
+#include "ghosts/HaloPlan.hpp"
 
 #include <utils/Vector.hpp>
 
@@ -85,20 +85,19 @@ struct RegularDecomposition : public ParticleDecomposition {
   std::vector<Cell> cells;
   std::vector<Cell *> m_local_cells;
   std::vector<Cell *> m_ghost_cells;
-  GhostCommunicator m_exchange_ghosts_comm;
-  GhostCommunicator m_collect_ghost_force_comm;
+  /** Topology-agnostic direct-neighbor halo plan (see @c make_halo_plan).
+   *  Holds ParticleList* into this decomposition's cells — value-copying this
+   *  object leaves these pointers dangling. TODO(Task 1.7): make non-copyable
+   * or rebuild-on-copy.
+   */
+  GhostComm::HaloPlan m_halo_plan;
 
 public:
   RegularDecomposition(boost::mpi::communicator comm, double range,
                        BoxGeometry const &box_geo, LocalBox const &local_geo,
                        std::optional<std::pair<int, int>> fully_connected);
 
-  GhostCommunicator const &exchange_ghosts_comm() const override {
-    return m_exchange_ghosts_comm;
-  }
-  GhostCommunicator const &collect_ghost_force_comm() const override {
-    return m_collect_ghost_force_comm;
-  }
+  GhostComm::HaloPlan const *halo_plan() const override { return &m_halo_plan; }
 
   std::span<Cell *const> local_cells() const override { return m_local_cells; }
   std::span<Cell *const> ghost_cells() const override { return m_ghost_cells; }
@@ -226,10 +225,12 @@ struct RegularDecomposition : public ParticleDecomposition {
    */
   void init_cell_interactions();
 
-  /** Create communicators for cell structure regular decomposition (see \ref
-   *  GhostCommunicator).
+  /** @brief Build a direct-neighbor @ref GhostComm::HaloPlan.
+   *
+   *  Every one of the up-to-26 stencil neighbors is addressed as an explicit
+   *  peer (no axis-by-axis relay). The result is cached in @ref m_halo_plan.
    */
-  GhostCommunicator prepare_comm();
+  GhostComm::HaloPlan make_halo_plan();
 
   /** Maximal number of cells per node. In order to avoid memory
    *  problems due to the cell grid, one has to specify the maximal
diff --git a/src/core/electrostatics/icc.cpp b/src/core/electrostatics/icc.cpp
index 63e8f7bd97..31084816e9 100644
--- a/src/core/electrostatics/icc.cpp
+++ b/src/core/electrostatics/icc.cpp
@@ -76,7 +76,7 @@ static void force_calc_icc(
   auto const reset_kernel = [](Particle &p) { p.force_and_torque() = {}; };
   cell_structure.for_each_local_particle(reset_kernel);
   cell_structure.for_each_ghost_particle(reset_kernel);
-  cell_structure.reset_local_force_and_torque();
+  cell_structure.reset_local_force_buffers();
 
   // calc ICC forces
   cell_structure.non_bonded_loop(
@@ -221,6 +221,11 @@ void ICCStar::iteration() {
       break;
   }
 
+  // The iterations above scattered coulomb forces into the local force
+  // buffers; leave them zeroed so the force calculation that follows this
+  // charge update does not accumulate on top of stale contributions.
+  cell_structure.reset_local_force_buffers();
+
   if (global_max_rel_diff > icc_cfg.convergence) {
     runtimeErrorMsg()
         << "ICC failed to converge in the given number of maximal steps.";
diff --git a/src/core/energy.cpp b/src/core/energy.cpp
index 945f47dc4c..7dae6004e9 100644
--- a/src/core/energy.cpp
+++ b/src/core/energy.cpp
@@ -41,6 +41,10 @@
 #include "electrostatics/coulomb.hpp"
 #include "magnetostatics/dipoles.hpp"
 
+#ifdef ESPRESSO_CALIPER
+#include "caliper_utils.hpp"
+#endif
+
 #include <cmath>
 #include <cstddef>
 #include <memory>
@@ -79,7 +83,7 @@ std::shared_ptr<Observable_stat> System::calculate_energy() {
   auto const dipoles_kernel = dipoles.pair_energy_kernel();
 
 #ifdef ESPRESSO_CALIPER
-  CALI_MARK_BEGIN("cabana_short_range");
+  ESPRESSO_CALI_MARK_BEGIN("cabana_short_range");
 #endif
   // Factory instead of an eager criterion: construction fills an O(n_types^2)
   // cutoff table, so it only runs on the link-cell fallback path.
@@ -141,7 +145,7 @@ std::shared_ptr<Observable_stat> System::calculate_energy() {
   reduce_cabana_energy(local_energy, layout, obs_energy, *bonded_ias,
                        nonbonded_ias->get_max_seen_particle_type() + 1);
 #ifdef ESPRESSO_CALIPER
-  CALI_MARK_END("cabana_short_range");
+  ESPRESSO_CALI_MARK_END("cabana_short_range");
 #endif
 
 #ifdef ESPRESSO_ELECTROSTATICS
diff --git a/src/core/forces.cpp b/src/core/forces.cpp
index 941f062f1c..b1f0c59e07 100644
--- a/src/core/forces.cpp
+++ b/src/core/forces.cpp
@@ -35,6 +35,7 @@
 #include "forces_init.hpp"
 #include "forces_inline.hpp"
 #include "galilei/ComFixed.hpp"
+#include "ghosts.hpp"
 #include "immersed_boundary/ImmersedBoundaries.hpp"
 #include "integrators/Propagation.hpp"
 #include "lb/particle_coupling.hpp"
@@ -56,7 +57,7 @@
 #include <utils/math/sqr.hpp>
 
 #ifdef ESPRESSO_CALIPER
-#include <caliper/cali.h>
+#include "caliper_utils.hpp"
 #endif
 
 #include <Cabana_Core.hpp>
@@ -68,6 +69,57 @@
 #include <span>
 #include <variant>
 
+/**
+ * @brief Eligibility check for the split-phase ghost force reduction.
+ *
+ * The overlap is safe only when ALL of the following hold:
+ *
+ *  1. More than one MPI rank: at 1 rank the reduce is a local-only copy,
+ *     there is nothing to hide, so the split adds overhead for no gain.
+ *  2. ComFixed is inactive: ComFixed applies a per-type force correction that
+ *     must see the final (post-reduce) forces, which arrive on the blocking
+ * path.
+ *  3. force_cap == 0: force capping requires the final forces (post-reduce).
+ *  4. Integrator is VV or symplectic Euler: only these inertial methods use
+ *     step_2; steepest-descent, BD, SD, and NPT are not eligible (NPT virial
+ *     is not accumulated during step_2 overlap; BD/SD have no step_2
+ * half-kick).
+ *  5. NPT propagation not in use: the NPT step_2 updates the box pressure and
+ *     box length using the virial gathered from forces; it must run after the
+ *     full reduce.  NPT also implies integ_switch != NVT/SE, so condition 4
+ *     already excludes it, but we check explicitly for clarity.
+ *  6. LB-tracer arm not active: TRANS_LB_TRACER particles have forces
+ * transferred to the LB fluid in the integrate loop *after* calculate_forces
+ * and before step_2; with the split-phase path the fluid coupling would run
+ * while the reduction is in flight, leading to incorrect force accumulation.
+ */
+static bool ghost_reduce_overlap_eligible(System::System const &system) {
+  // 1. Must have more than one MPI rank.
+  if (::comm_cart.size() <= 1)
+    return false;
+  // 2. ComFixed must be inactive.
+  if (not system.comfixed->get_fixed_types().empty())
+    return false;
+  // 3. Force capping must be off.
+  if (system.get_force_cap() != 0.)
+    return false;
+  // 4. NPT propagation not in use (belt-and-suspenders after check 6).
+#ifdef ESPRESSO_NPT
+  if (system.propagation->used_propagations &
+      PropagationMode::TRANS_LANGEVIN_NPT)
+    return false;
+#endif
+  // 5. LB-tracer arm not active.
+#ifdef ESPRESSO_VIRTUAL_SITES_INERTIALESS_TRACERS
+  if (system.propagation->used_propagations & PropagationMode::TRANS_LB_TRACER)
+    return false;
+#endif
+  // 6. Integrator must be VV or symplectic Euler (inertial with step_2
+  // half-kick).
+  auto const integ = system.propagation->integ_switch;
+  return integ == INTEG_METHOD_NVT or integ == INTEG_METHOD_SYMPLECTIC_EULER;
+}
+
 static void force_capping(CellStructure &cell_structure, double force_cap) {
   if (force_cap > 0.) {
     auto const force_cap_sq = Utils::sqr(force_cap);
@@ -263,20 +315,31 @@ static ShortRangeVerletPairLoop create_specialized_verlet_pair_loop(
 
 static void reduce_cabana_forces_and_torques(System::System const &system,
                                              Utils::Vector3d *virial) {
+#ifdef ESPRESSO_CALIPER
+  ESPRESSO_CALI_MARK_FUNCTION;
+#endif
 
   auto const &unique_particles = system.cell_structure->get_unique_particles();
   auto &local_force = system.cell_structure->get_local_force();
   auto scatter_force = system.cell_structure->get_scatter_force();
   Kokkos::Experimental::contribute(local_force, scatter_force);
 #ifdef ESPRESSO_ROTATION
+  // when no kernel scattered into the torque buffers this pass, they are
+  // all-zero and both the reduction over replicas and the per-particle
+  // accumulation would only add zeros
+  auto const reduce_torque = system.cell_structure->torque_replicas_dirty();
   auto &local_torque = system.cell_structure->get_local_torque();
-  auto scatter_torque = system.cell_structure->get_scatter_torque();
-  Kokkos::Experimental::contribute(local_torque, scatter_torque);
+  if (reduce_torque) {
+    auto scatter_torque = system.cell_structure->get_scatter_torque();
+    Kokkos::Experimental::contribute(local_torque, scatter_torque);
+  }
 #endif
 #ifdef ESPRESSO_NPT
   auto &local_virial = system.cell_structure->get_local_virial();
-  auto scatter_virial = system.cell_structure->get_scatter_virial();
-  Kokkos::Experimental::contribute(local_virial, scatter_virial);
+  if (system.cell_structure->virial_replicas_dirty()) {
+    auto scatter_virial = system.cell_structure->get_scatter_virial();
+    Kokkos::Experimental::contribute(local_virial, scatter_virial);
+  }
 #endif
 
   using execution_space = Kokkos::DefaultHostExecutionSpace;
@@ -285,24 +348,22 @@ static void reduce_cabana_forces_and_torques(System::System const &system,
   Kokkos::parallel_for("reduction", policy,
                        [&local_force,
 #ifdef ESPRESSO_ROTATION
-                        &local_torque,
+                        &local_torque, reduce_torque,
 #endif
                         &unique_particles](std::size_t const i) {
                          Utils::Vector3d force{};
-#ifdef ESPRESSO_ROTATION
-                         Utils::Vector3d torque{};
-#endif
                          force[0] += local_force(i, 0);
                          force[1] += local_force(i, 1);
                          force[2] += local_force(i, 2);
-#ifdef ESPRESSO_ROTATION
-                         torque[0] += local_torque(i, 0);
-                         torque[1] += local_torque(i, 1);
-                         torque[2] += local_torque(i, 2);
-#endif
                          unique_particles.at(i)->force() += force;
 #ifdef ESPRESSO_ROTATION
-                         unique_particles.at(i)->torque() += torque;
+                         if (reduce_torque) {
+                           Utils::Vector3d torque{};
+                           torque[0] += local_torque(i, 0);
+                           torque[1] += local_torque(i, 1);
+                           torque[2] += local_torque(i, 2);
+                           unique_particles.at(i)->torque() += torque;
+                         }
 #endif
                        });
   Kokkos::fence();
@@ -318,16 +379,16 @@ static void reduce_cabana_forces_and_torques(System::System const &system,
 
 void System::System::calculate_forces() {
 #ifdef ESPRESSO_CALIPER
-  CALI_CXX_MARK_FUNCTION;
+  ESPRESSO_CALI_MARK_FUNCTION;
 #endif
 #ifdef ESPRESSO_CUDA
   {
 #ifdef ESPRESSO_CALIPER
-    CALI_MARK_BEGIN("copy_particles_to_GPU");
+    ESPRESSO_CALI_MARK_BEGIN("copy_particles_to_GPU");
 #endif
     gpu->update();
 #ifdef ESPRESSO_CALIPER
-    CALI_MARK_END("copy_particles_to_GPU");
+    ESPRESSO_CALI_MARK_END("copy_particles_to_GPU");
 #endif
   }
 #endif // ESPRESSO_CUDA
@@ -379,7 +440,7 @@ void System::System::calculate_forces() {
 #endif // ESPRESSO_ELECTROSTATICS
   init_forces_and_thermostat(*this);
 #ifdef ESPRESSO_CALIPER
-  CALI_MARK_BEGIN("calc_long_range_forces");
+  ESPRESSO_CALI_MARK_BEGIN("calc_long_range_forces");
 #endif
 #ifdef ESPRESSO_ELECTROSTATICS
   coulomb.calc_long_range_force();
@@ -388,11 +449,11 @@ void System::System::calculate_forces() {
   dipoles.calc_long_range_force();
 #endif
 #ifdef ESPRESSO_CALIPER
-  CALI_MARK_END("calc_long_range_forces");
+  ESPRESSO_CALI_MARK_END("calc_long_range_forces");
 #endif
 
 #ifdef ESPRESSO_CALIPER
-  CALI_MARK_BEGIN("cabana_short_range");
+  ESPRESSO_CALI_MARK_BEGIN("cabana_short_range");
 #endif
   auto &bs = cell_structure->bond_state();
   auto bonds_kernel_data = create_kokkos_bonds_kernel_data(*this);
@@ -410,6 +471,48 @@ void System::System::calculate_forces() {
   auto const specialized_pair_loop = create_specialized_verlet_pair_loop(
       *this, virial, elc_kernel, coulomb_kernel, dipoles_kernel);
 
+  // The specialized pair kernel scatters only into the force view. Record
+  // before the launch whether this pass can write the torque/virial scatter
+  // buffers at all, so their O(n_threads * N) zeroing and reduction can be
+  // skipped otherwise (see CellStructure::reset_local_properties). All other
+  // torque sources (Langevin rotation, virtual sites back-transfer, the
+  // dipolar solvers except dp3m, constraints) write Particle::torque()
+  // directly and never touch the scatter buffers; dp3m marks the flag at its
+  // own scatter site.
+  [[maybe_unused]] auto const generic_pair_path =
+      get_interaction_range() > 0. and
+      (not specialized_pair_loop or
+       propagation->integ_switch == INTEG_METHOD_STEEPEST_DESCENT or
+       not cell_structure->use_verlet_list);
+#ifdef ESPRESSO_ROTATION
+  // Within the generic kernel, only orientation-dependent pair potentials
+  // (Gay-Berne) and the dipolar pair kernel scatter into the torque view.
+  auto const gay_berne_active =
+      nonbonded_ias->pair_potential_active(PairPotential::GayBerne);
+  auto const dipolar_pair_kernel_active =
+#ifdef ESPRESSO_DIPOLES
+      get_ptr(dipoles_kernel) != nullptr;
+#else
+      false;
+#endif
+  if (generic_pair_path and (gay_berne_active or dipolar_pair_kernel_active)) {
+    cell_structure->mark_torque_replicas_dirty();
+    // Pair-kernel torques also land on ghost particles and only get home via
+    // the TORQUE ghost reduce. Each torque-scattering kernel therefore needs
+    // a matching arm in orientation_ghosts_needed() (System.cpp).
+    assert(get_force_reduce_ghost_flags() & GHOSTTRANS_TORQUE);
+  }
+#endif
+#ifdef ESPRESSO_NPT
+  auto const have_bonds = cell_structure->get_local_pair_bond_numbers() > 0 or
+                          cell_structure->get_local_angle_bond_numbers() > 0 or
+                          cell_structure->get_local_dihedral_bond_numbers() > 0;
+  if (generic_pair_path or
+      (bonded_ias->maximal_cutoff() >= 0. and have_bonds)) {
+    cell_structure->mark_virial_replicas_dirty();
+  }
+#endif
+
   cabana_short_range(pair_bonds_kernel, angle_bonds_kernel,
                      dihedral_bonds_kernel, first_neighbor_kernel,
                      *cell_structure, get_interaction_range(),
@@ -432,7 +535,7 @@ void System::System::calculate_forces() {
 #endif // ESPRESSO_COLLISION_DETECTION
 
 #ifdef ESPRESSO_CALIPER
-  CALI_MARK_END("cabana_short_range");
+  ESPRESSO_CALI_MARK_END("cabana_short_range");
 #endif
 
   constraints->add_forces(particles, get_sim_time());
@@ -444,18 +547,18 @@ void System::System::calculate_forces() {
   if (thermostat->lb and (propagation->used_propagations &
                           PropagationMode::TRANS_LB_MOMENTUM_EXCHANGE)) {
 #ifdef ESPRESSO_CALIPER
-    CALI_MARK_BEGIN("lb_particle_coupling");
+    ESPRESSO_CALI_MARK_BEGIN("lb_particle_coupling");
 #endif
     lb_couple_particles();
 #ifdef ESPRESSO_CALIPER
-    CALI_MARK_END("lb_particle_coupling");
+    ESPRESSO_CALI_MARK_END("lb_particle_coupling");
 #endif
   }
 
 #ifdef ESPRESSO_CUDA
   {
 #ifdef ESPRESSO_CALIPER
-    CALI_MARK_BEGIN("copy_forces_from_GPU");
+    ESPRESSO_CALI_MARK_BEGIN("copy_forces_from_GPU");
 #endif
     gpu->copy_forces_to_host(particles, this_node);
 
@@ -464,7 +567,7 @@ void System::System::calculate_forces() {
 #endif
 
 #ifdef ESPRESSO_CALIPER
-    CALI_MARK_END("copy_forces_from_GPU");
+    ESPRESSO_CALI_MARK_END("copy_forces_from_GPU");
 #endif
   }
 #endif // ESPRESSO_CUDA
@@ -483,14 +586,28 @@ void System::System::calculate_forces() {
   }
 #endif
 
-  // Communication step: ghost forces
-  cell_structure->ghosts_reduce_forces();
-
-  // should be pretty late, since it needs to zero out the total force
-  comfixed->apply(particles);
-
-  // Needs to be the last one to be effective
-  force_capping(*cell_structure, force_cap);
+  if (ghost_reduce_overlap_eligible(*this)) {
+    // Split-phase path: start the ghost force reduction now (non-blocking);
+    // integrate.cpp will run interior-cell step_2, finish the reduce, then
+    // run boundary-cell step_2.  comfixed and force_capping are inactive on
+    // this path (asserted by eligibility; the assert below is a belt-and-
+    // suspenders cross-check).
+    assert(comfixed->get_fixed_types().empty() &&
+           "ghost_reduce_overlap: comfixed must be inactive on the eligible "
+           "path");
+    assert(force_cap == 0. &&
+           "ghost_reduce_overlap: force_cap must be 0 on the eligible path");
+    cell_structure->ghosts_reduce_forces_start();
+  } else {
+    // Blocking path: finish the reduce here, then apply comfixed/capping.
+    cell_structure->ghosts_reduce_forces();
+
+    // should be pretty late, since it needs to zero out the total force
+    comfixed->apply(particles);
+
+    // Needs to be the last one to be effective
+    force_capping(*cell_structure, force_cap);
+  }
 
   // mark that forces are now up-to-date
   propagation->recalc_forces = false;
diff --git a/src/core/forces_init.cpp b/src/core/forces_init.cpp
index fc35d0c069..899b2b4cef 100644
--- a/src/core/forces_init.cpp
+++ b/src/core/forces_init.cpp
@@ -39,6 +39,8 @@
 #include "thermostats/langevin_inline.hpp"
 
 #ifdef ESPRESSO_CALIPER
+#include "caliper_utils.hpp"
+
 #include <caliper/cali.h>
 #endif
 
@@ -67,7 +69,7 @@ static ParticleForce external_force(Particle const &p) {
 /** Combined force initialization and Langevin noise application */
 void init_forces_and_thermostat(System::System const &system) {
 #ifdef ESPRESSO_CALIPER
-  CALI_CXX_MARK_FUNCTION;
+  ESPRESSO_CALI_MARK_FUNCTION;
 #endif
 
   auto &cell_structure = *system.cell_structure;
@@ -99,7 +101,8 @@ void init_forces_and_thermostat(System::System const &system) {
 #endif
     }
   });
-  cell_structure.reset_local_force_and_torque();
+  // The local force/torque scatter buffers were already zeroed this force
+  // call by prepare_verlet_list_cabana() (via update_verlet_state).
 
   // Initialize ghost forces (unchanged)
   cell_structure.ghosts_reset_forces();
diff --git a/src/core/ghosts.cpp b/src/core/ghosts.cpp
deleted file mode 100644
index ded254f9e5..0000000000
--- a/src/core/ghosts.cpp
+++ /dev/null
@@ -1,574 +0,0 @@
-/*
- * Copyright (C) 2010-2026 The ESPResSo project
- * Copyright (C) 2002,2003,2004,2005,2006,2007,2008,2009,2010
- *   Max-Planck-Institute for Polymer Research, Theory Group
- *
- * 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/>.
- */
-/** \file
- *  Ghost particles and particle exchange.
- *
- *  For more information on ghosts,
- *  see \ref ghosts.hpp "ghosts.hpp"
- *
- * Note on variable naming:
- * - a "GhostCommunicator" is always named "gcr",
- * - a "GhostCommunication" is always named "ghost_comm".
- */
-#include "ghosts.hpp"
-
-#include "BoxGeometry.hpp"
-#include "Particle.hpp"
-#include "system/System.hpp"
-
-#include <utils/serialization/memcpy_archive.hpp>
-
-#include <boost/archive/binary_iarchive.hpp>
-#include <boost/archive/binary_oarchive.hpp>
-#include <boost/iostreams/device/array.hpp>
-#include <boost/iostreams/device/back_inserter.hpp>
-#include <boost/iostreams/stream.hpp>
-#include <boost/mpi/collectives.hpp>
-#include <boost/range/numeric.hpp>
-#include <boost/serialization/vector.hpp>
-
-#include <algorithm>
-#include <cassert>
-#include <cstddef>
-#include <functional>
-#include <iterator>
-#include <limits>
-#include <span>
-#include <vector>
-
-/** Tag for ghosts communications. */
-#define REQ_GHOST_SEND 100
-
-/**
- * Class that stores marshalled data for ghost communications.
- * To store and retrieve data, use the adapter classes below.
- */
-class CommBuf {
-public:
-  /** Returns a pointer to the non-bond storage.
-   */
-  char *data() { return buf.data(); }
-  const char *data() const { return buf.data(); }
-
-  /** Returns the number of elements in the non-bond storage.
-   */
-  std::size_t size() const { return buf.size(); }
-
-  /** Resizes the underlying storage s.t. the object is capable
-   * of holding "new_size" chars.
-   * @param new_size new size
-   */
-  void resize(std::size_t new_size) { buf.resize(new_size); }
-
-  /** Returns a reference to the bond storage.
-   */
-  auto &bonds() { return bondbuf; }
-  const auto &bonds() const { return bondbuf; }
-
-  auto make_span() { return std::span(buf.data(), buf.size()); }
-
-private:
-  std::vector<char> buf;     ///< Buffer for everything but bonds
-  std::vector<char> bondbuf; ///< Buffer for bond lists
-};
-
-/** @brief Pseudo-archive to calculate the size of the serialization buffer. */
-class SerializationSizeCalculator {
-  std::size_t m_size = 0;
-
-public:
-  auto size() const { return m_size; }
-
-  template <class T> auto &operator<<(T &) {
-    m_size += sizeof(T);
-    return *this;
-  }
-
-  template <class T> auto &operator&(T &t) { return *this << t; }
-};
-
-/** @brief Type of reduction to carry out during serialization. */
-enum class ReductionPolicy {
-  /** @brief Reduction for domain-to-domain particle communication. */
-  MOVE,
-  /** @brief Reduction for cell-to-cell particle update. */
-  UPDATE,
-};
-
-/** @brief Whether to save the state to or load the state from the archive. */
-enum class SerializationDirection { SAVE, LOAD };
-
-/**
- * @brief Serialize particle data, possibly with reduction.
- * The reduction can take place during the save stage, e.g. to apply
- * a ghost shift to the particle position, or during the load stage,
- * e.g. to transfer momentum between particles in two local cells.
- */
-template <class Archive>
-static void
-serialize_and_reduce(Archive &ar, Particle &p, unsigned int data_parts,
-                     ReductionPolicy policy, SerializationDirection direction,
-                     BoxGeometry const &box_geo,
-                     Utils::Vector3d const *ghost_shift) {
-  if (data_parts & GHOSTTRANS_PROPRTS) {
-    ar & p.id() & p.mol_id() & p.type() & p.propagation();
-#ifdef ESPRESSO_ROTATION
-    ar & p.rotation();
-#ifdef ESPRESSO_ROTATIONAL_INERTIA
-    ar & p.rinertia();
-#endif
-#endif
-#ifdef ESPRESSO_MASS
-    ar & p.mass();
-#endif
-#ifdef ESPRESSO_ELECTROSTATICS
-    ar & p.q();
-#endif
-#ifdef ESPRESSO_DIPOLES
-    ar & p.dipm();
-#endif
-#ifdef ESPRESSO_LB_ELECTROHYDRODYNAMICS
-    ar & p.mu_E();
-#endif
-#ifdef ESPRESSO_VIRTUAL_SITES_RELATIVE
-    ar & p.vs_relative();
-#endif
-#ifdef ESPRESSO_THERMOSTAT_PER_PARTICLE
-    ar & p.gamma();
-#ifdef ESPRESSO_ROTATION
-    ar & p.gamma_rot();
-#endif
-#endif
-#ifdef ESPRESSO_EXTERNAL_FORCES
-    ar & p.fixed();
-    ar & p.ext_force();
-#ifdef ESPRESSO_ROTATION
-    ar & p.ext_torque();
-#endif
-#endif
-#ifdef ESPRESSO_ENGINE
-    ar & p.swimming();
-#endif
-  }
-  if (data_parts & GHOSTTRANS_POSITION) {
-    if (direction == SerializationDirection::SAVE and ghost_shift != nullptr) {
-      /* ok, this is not nice, but perhaps fast */
-      auto pos = p.pos() + *ghost_shift;
-      auto img = p.image_box();
-      box_geo.fold_position(pos, img);
-      ar & pos;
-      ar & img;
-    } else {
-      ar & p.pos();
-      ar & p.image_box();
-    }
-#ifdef ESPRESSO_ROTATION
-    ar & p.quat();
-#endif
-#ifdef ESPRESSO_BOND_CONSTRAINT
-    ar & p.pos_last_time_step();
-#endif
-  }
-  if (data_parts & GHOSTTRANS_MOMENTUM) {
-    ar & p.v();
-#ifdef ESPRESSO_ROTATION
-    ar & p.omega();
-#endif
-  }
-  if (data_parts & GHOSTTRANS_FORCE) {
-    if (policy == ReductionPolicy::UPDATE and
-        direction == SerializationDirection::LOAD) {
-      Utils::Vector3d force;
-      ar & force;
-      p.force() += force;
-    } else {
-      ar & p.force();
-    }
-#ifdef ESPRESSO_ROTATION
-    if (policy == ReductionPolicy::UPDATE and
-        direction == SerializationDirection::LOAD) {
-      Utils::Vector3d torque;
-      ar & torque;
-      p.torque() += torque;
-    } else {
-      ar & p.torque();
-    }
-#endif
-  }
-#ifdef ESPRESSO_BOND_CONSTRAINT
-  if (data_parts & GHOSTTRANS_RATTLE) {
-    if (policy == ReductionPolicy::UPDATE and
-        direction == SerializationDirection::LOAD) {
-      Utils::Vector3d correction;
-      ar & correction;
-      p.rattle_correction() += correction;
-    } else {
-      ar & p.rattle_correction();
-    }
-  }
-#endif
-}
-
-static auto calc_transmit_size(BoxGeometry const &box_geo,
-                               unsigned data_parts) {
-  SerializationSizeCalculator sizeof_archive;
-  Particle p{};
-  serialize_and_reduce(sizeof_archive, p, data_parts, ReductionPolicy::MOVE,
-                       SerializationDirection::SAVE, box_geo, nullptr);
-  return sizeof_archive.size();
-}
-
-static auto calc_transmit_size(GhostCommunication const &ghost_comm,
-                               BoxGeometry const &box_geo,
-                               unsigned int data_parts) {
-  if (data_parts & GHOSTTRANS_PARTNUM)
-    return sizeof(unsigned int) * ghost_comm.part_lists.size();
-
-  auto const n_part = boost::accumulate(
-      ghost_comm.part_lists, std::size_t{0},
-      [](std::size_t sum, auto part_list) { return sum + part_list->size(); });
-
-  return n_part * calc_transmit_size(box_geo, data_parts);
-}
-
-static void prepare_send_buffer(CommBuf &send_buffer,
-                                GhostCommunication const &ghost_comm,
-                                BoxGeometry const &box_geo,
-                                unsigned int data_parts) {
-
-  /* reallocate send buffer */
-  send_buffer.resize(calc_transmit_size(ghost_comm, box_geo, data_parts));
-  send_buffer.bonds().clear();
-
-  auto archiver = Utils::MemcpyOArchive{send_buffer.make_span()};
-
-  /* Construct archive that pushes back to the bond buffer */
-  namespace io = boost::iostreams;
-  io::stream<io::back_insert_device<std::vector<char>>> os{
-      io::back_inserter(send_buffer.bonds())};
-  boost::archive::binary_oarchive bond_archiver{os};
-
-  /* put in data */
-  for (auto part_list : ghost_comm.part_lists) {
-    if (data_parts & GHOSTTRANS_PARTNUM) {
-      assert(part_list->size() <= std::numeric_limits<unsigned int>::max());
-      auto np = static_cast<unsigned int>(part_list->size());
-      archiver << np;
-    } else {
-      for (auto &p : *part_list) {
-        serialize_and_reduce(archiver, p, data_parts, ReductionPolicy::MOVE,
-                             SerializationDirection::SAVE, box_geo,
-                             &ghost_comm.shift);
-        if (data_parts & GHOSTTRANS_BONDS) {
-          bond_archiver << p.bonds();
-        }
-      }
-    }
-  }
-
-  assert(archiver.bytes_written() == send_buffer.size());
-}
-
-static void prepare_ghost_cell(ParticleList *cell, std::size_t size) {
-  /* Adapt size */
-  cell->resize(size);
-
-  /* Mark particles as ghosts */
-  for (auto &p : *cell) {
-    p.set_ghost(true);
-  }
-}
-
-static void prepare_recv_buffer(CommBuf &recv_buffer,
-                                GhostCommunication const &ghost_comm,
-                                BoxGeometry const &box_geo,
-                                unsigned int data_parts) {
-  /* reallocate recv buffer */
-  recv_buffer.resize(calc_transmit_size(ghost_comm, box_geo, data_parts));
-  /* clear bond buffer */
-  recv_buffer.bonds().clear();
-}
-
-static void put_recv_buffer(CommBuf &recv_buffer,
-                            GhostCommunication const &ghost_comm,
-                            BoxGeometry const &box_geo,
-                            unsigned int data_parts) {
-  /* put back data */
-  auto archiver = Utils::MemcpyIArchive{recv_buffer.make_span()};
-
-  if (data_parts & GHOSTTRANS_PARTNUM) {
-    for (auto part_list : ghost_comm.part_lists) {
-      unsigned int np;
-      archiver >> np;
-      prepare_ghost_cell(part_list, np);
-    }
-  } else {
-    for (auto part_list : ghost_comm.part_lists) {
-      for (auto &p : *part_list) {
-        serialize_and_reduce(archiver, p, data_parts, ReductionPolicy::MOVE,
-                             SerializationDirection::LOAD, box_geo, nullptr);
-      }
-    }
-    if (data_parts & GHOSTTRANS_BONDS) {
-      namespace io = boost::iostreams;
-      io::stream<io::array_source> bond_stream(io::array_source{
-          recv_buffer.bonds().data(), recv_buffer.bonds().size()});
-      boost::archive::binary_iarchive bond_archiver(bond_stream);
-
-      for (auto part_list : ghost_comm.part_lists) {
-        for (auto &p : *part_list) {
-          bond_archiver >> p.bonds();
-        }
-      }
-    }
-  }
-
-  assert(archiver.bytes_read() == recv_buffer.size());
-
-  recv_buffer.bonds().clear();
-}
-
-#ifdef ESPRESSO_BOND_CONSTRAINT
-static void
-add_rattle_correction_from_recv_buffer(CommBuf &recv_buffer,
-                                       const GhostCommunication &ghost_comm) {
-  /* put back data */
-  auto archiver = Utils::MemcpyIArchive{recv_buffer.make_span()};
-  for (auto &part_list : ghost_comm.part_lists) {
-    for (Particle &part : *part_list) {
-      ParticleRattle pr;
-      archiver >> pr;
-      part.rattle_params() += pr;
-    }
-  }
-}
-#endif
-
-static void add_forces_from_recv_buffer(CommBuf &recv_buffer,
-                                        const GhostCommunication &ghost_comm) {
-  /* put back data */
-  auto archiver = Utils::MemcpyIArchive{recv_buffer.make_span()};
-  for (auto &part_list : ghost_comm.part_lists) {
-    for (Particle &part : *part_list) {
-      ParticleForce pf;
-      archiver >> pf;
-      part.force_and_torque() += pf;
-    }
-  }
-}
-
-static void cell_cell_transfer(GhostCommunication const &ghost_comm,
-                               BoxGeometry const &box_geo,
-                               unsigned int data_parts) {
-  CommBuf buffer;
-  if (!(data_parts & GHOSTTRANS_PARTNUM)) {
-    buffer.resize(calc_transmit_size(box_geo, data_parts));
-  }
-  /* transfer data */
-  auto const offset = ghost_comm.part_lists.size() / 2;
-  for (std::size_t pl = 0; pl < offset; pl++) {
-    auto *src_list = ghost_comm.part_lists[pl];
-    auto *dst_list = ghost_comm.part_lists[pl + offset];
-
-    if (data_parts & GHOSTTRANS_PARTNUM) {
-      prepare_ghost_cell(dst_list, src_list->size());
-    } else {
-      auto &src_part = *src_list;
-      auto &dst_part = *dst_list;
-      assert(src_part.size() == dst_part.size());
-
-      for (std::size_t i = 0; i < src_part.size(); i++) {
-        auto ar_out = Utils::MemcpyOArchive{buffer.make_span()};
-        auto ar_in = Utils::MemcpyIArchive{buffer.make_span()};
-        auto &p1 = src_part.begin()[i];
-        auto &p2 = dst_part.begin()[i];
-        serialize_and_reduce(ar_out, p1, data_parts, ReductionPolicy::UPDATE,
-                             SerializationDirection::SAVE, box_geo,
-                             &ghost_comm.shift);
-        serialize_and_reduce(ar_in, p2, data_parts, ReductionPolicy::UPDATE,
-                             SerializationDirection::LOAD, box_geo, nullptr);
-        if (data_parts & GHOSTTRANS_BONDS) {
-          p2.bonds() = p1.bonds();
-        }
-      }
-    }
-  }
-}
-
-static bool is_send_op(int comm_type, int node, int this_node) {
-  return ((comm_type == GHOST_SEND) || (comm_type == GHOST_RDCE) ||
-          (comm_type == GHOST_BCST && node == this_node));
-}
-
-static bool is_recv_op(int comm_type, int node, int this_node) {
-  return ((comm_type == GHOST_RECV) ||
-          (comm_type == GHOST_BCST && node != this_node) ||
-          (comm_type == GHOST_RDCE && node == this_node));
-}
-
-static bool is_prefetchable(GhostCommunication const &ghost_comm,
-                            int this_node) {
-  int const comm_type = ghost_comm.type & GHOST_JOBMASK;
-  int const prefetch = ghost_comm.type & GHOST_PREFETCH;
-  int const node = ghost_comm.node;
-  return is_send_op(comm_type, node, this_node) && prefetch;
-}
-
-static bool is_poststorable(GhostCommunication const &ghost_comm,
-                            int this_node) {
-  int const comm_type = ghost_comm.type & GHOST_JOBMASK;
-  int const poststore = ghost_comm.type & GHOST_PSTSTORE;
-  int const node = ghost_comm.node;
-  return is_recv_op(comm_type, node, this_node) && poststore;
-}
-
-void ghost_communicator(GhostCommunicator const &gcr,
-                        BoxGeometry const &box_geo, unsigned int data_parts) {
-  if (GHOSTTRANS_NONE == data_parts)
-    return;
-
-  static CommBuf send_buffer, recv_buffer;
-
-  auto const &comm = gcr.mpi_comm;
-
-  for (auto cit = gcr.communications.cbegin(); cit != gcr.communications.cend();
-       ++cit) {
-    auto const &ghost_comm = *cit;
-    int const comm_type = ghost_comm.type & GHOST_JOBMASK;
-
-    if (comm_type == GHOST_LOCL) {
-      cell_cell_transfer(ghost_comm, box_geo, data_parts);
-      continue;
-    }
-
-    int const prefetch = ghost_comm.type & GHOST_PREFETCH;
-    int const poststore = ghost_comm.type & GHOST_PSTSTORE;
-    int const node = ghost_comm.node;
-
-    /* prepare send buffer if necessary */
-    if (is_send_op(comm_type, node, comm.rank())) {
-      /* ok, we send this step, prepare send buffer if not yet done */
-      if (!prefetch) {
-        prepare_send_buffer(send_buffer, ghost_comm, box_geo, data_parts);
-      }
-      // Check prefetched send buffers (must also hold for buffers allocated
-      // in the previous lines.)
-      assert(send_buffer.size() ==
-             calc_transmit_size(ghost_comm, box_geo, data_parts));
-    } else if (prefetch) {
-      /* we do not send this time, let's look for a prefetch */
-      auto prefetch_ghost_comm =
-          std::find_if(std::next(cit), gcr.communications.cend(),
-                       [this_node = comm.rank()](auto const &other_ghost_comm) {
-                         return is_prefetchable(other_ghost_comm, this_node);
-                       });
-
-      if (prefetch_ghost_comm != gcr.communications.end())
-        prepare_send_buffer(send_buffer, *prefetch_ghost_comm, box_geo,
-                            data_parts);
-    }
-
-    /* recv buffer for recv and multinode operations to this node */
-    if (is_recv_op(comm_type, node, comm.rank()))
-      prepare_recv_buffer(recv_buffer, ghost_comm, box_geo, data_parts);
-
-    /* transfer data */
-    // Use two send/recvs in order to avoid having to serialize CommBuf
-    // (which consists of already serialized data).
-    switch (comm_type) {
-    case GHOST_RECV:
-      comm.recv(node, REQ_GHOST_SEND, recv_buffer.data(),
-                static_cast<int>(recv_buffer.size()));
-      comm.recv(node, REQ_GHOST_SEND, recv_buffer.bonds());
-      break;
-    case GHOST_SEND:
-      comm.send(node, REQ_GHOST_SEND, send_buffer.data(),
-                static_cast<int>(send_buffer.size()));
-      comm.send(node, REQ_GHOST_SEND, send_buffer.bonds());
-      break;
-    case GHOST_BCST:
-      if (node == comm.rank()) {
-        boost::mpi::broadcast(comm, send_buffer.data(),
-                              static_cast<int>(send_buffer.size()), node);
-        boost::mpi::broadcast(comm, send_buffer.bonds(), node);
-      } else {
-        boost::mpi::broadcast(comm, recv_buffer.data(),
-                              static_cast<int>(recv_buffer.size()), node);
-        boost::mpi::broadcast(comm, recv_buffer.bonds(), node);
-      }
-      break;
-    case GHOST_RDCE:
-      if (node == comm.rank())
-        boost::mpi::reduce(
-            comm, reinterpret_cast<double *>(send_buffer.data()),
-            static_cast<int>(send_buffer.size() / sizeof(double)),
-            reinterpret_cast<double *>(recv_buffer.data()), std::plus<double>{},
-            node);
-      else
-        boost::mpi::reduce(
-            comm, reinterpret_cast<double *>(send_buffer.data()),
-            static_cast<int>(send_buffer.size() / sizeof(double)),
-            std::plus<double>{}, node);
-      break;
-    }
-
-    // recv op; write back data directly, if no PSTSTORE delay is requested.
-    if (is_recv_op(comm_type, node, comm.rank())) {
-      if (!poststore) {
-        /* forces have to be added, the rest overwritten. Exception is RDCE,
-         * where the addition is integrated into the communication. */
-        if (data_parts == GHOSTTRANS_FORCE && comm_type != GHOST_RDCE)
-          add_forces_from_recv_buffer(recv_buffer, ghost_comm);
-#ifdef ESPRESSO_BOND_CONSTRAINT
-        else if (data_parts == GHOSTTRANS_RATTLE && comm_type != GHOST_RDCE)
-          add_rattle_correction_from_recv_buffer(recv_buffer, ghost_comm);
-#endif
-        else
-          put_recv_buffer(recv_buffer, ghost_comm, box_geo, data_parts);
-      }
-    } else if (poststore) {
-      /* send op; write back delayed data from last recv, when this was a
-       * prefetch send. */
-      /* find previous action where we recv and which has PSTSTORE set */
-      auto poststore_ghost_comm = std::find_if(
-          std::make_reverse_iterator(cit), gcr.communications.crend(),
-          [this_node = comm.rank()](auto const &other_ghost_comm) {
-            return is_poststorable(other_ghost_comm, this_node);
-          });
-
-      if (poststore_ghost_comm != gcr.communications.rend()) {
-        assert(recv_buffer.size() ==
-               calc_transmit_size(*poststore_ghost_comm, box_geo, data_parts));
-        /* as above */
-        if (data_parts == GHOSTTRANS_FORCE && comm_type != GHOST_RDCE)
-          add_forces_from_recv_buffer(recv_buffer, *poststore_ghost_comm);
-#ifdef ESPRESSO_BOND_CONSTRAINT
-        else if (data_parts == GHOSTTRANS_RATTLE && comm_type != GHOST_RDCE)
-          add_rattle_correction_from_recv_buffer(recv_buffer,
-                                                 *poststore_ghost_comm);
-#endif
-        else
-          put_recv_buffer(recv_buffer, *poststore_ghost_comm, box_geo,
-                          data_parts);
-      }
-    }
-  }
-}
diff --git a/src/core/ghosts.hpp b/src/core/ghosts.hpp
index 857595e3bc..9e9a155b6e 100644
--- a/src/core/ghosts.hpp
+++ b/src/core/ghosts.hpp
@@ -24,104 +24,13 @@
 /** \file
  *  Ghost particles and particle exchange.
  *
- *  In this file you find everything concerning the exchange of
- *  particle data (particles, ghosts, positions and forces) for short
- *  range interactions between the spatial domains of neighboring
- *  nodes.
- *
- *  <h2> How does this work </h2>
- *  The ghost communication transfers data from cells on one node to cells on
- *  another node during the integration process. Note that data can only be
- *  transferred from one cell to one other, and the contents of the other cell
- *  will be overwritten.
- *
- *  The particle data that has to be transferred, and especially from where to
- *  where, heavily depends on the cell system. In ESPResSo, this is abstracted
- *  in form of ghost communicators (\ref GhostCommunicator) and ghost
- *  communications (\ref GhostCommunication). The ghost communicators represent
- *  the four communications above and consist of the data to transfer (which is
- *  determined by their type) and a list of ghost communications. The data
- *  types are described by the particle data classes:
- *  - @ref GHOSTTRANS_PROPRTS transfers the @ref ParticleProperties
- *  - @ref GHOSTTRANS_POSITION transfers the @ref ParticlePosition
- *  - @ref GHOSTTRANS_MOMENTUM transfers the @ref ParticleMomentum
- *  - @ref GHOSTTRANS_FORCE transfers the @ref ParticleForce
- *  - @ref GHOSTTRANS_RATTLE transfers the @ref ParticleRattle
- *  - @ref GHOSTTRANS_PARTNUM transfers the cell sizes
- *
- *  Each ghost communication describes a single communication of the local with
- *  another node (or all other nodes). The data transferred can be any number
- *  of cells, there are five communication types:
- *  - @ref GHOST_SEND sends data to one other node
- *  - @ref GHOST_RECV recvs data from one other node. In the case of
- *    forces, they are added up, everything else is overwritten
- *  - @ref GHOST_BCST sends data to all nodes
- *  - @ref GHOST_RDCE recvs data from all nodes. In the case of
- *    forces, they are added up, everything else is overwritten
- *  - @ref GHOST_LOCL transfer data from a local cell to another local cell.
- *    In this case, the first half of the cells are the sending cells, the
- *    other half are the corresponding receivers.
- *
- *  Note that for the first four communications you have to make sure that
- *  when one node sends to another, that the sender has @ref GHOST_SEND
- *  and the receiver @ref GHOST_RECV. In the case of @ref GHOST_BCST resp.
- *  @ref GHOST_RDCE, all nodes have to have the same communication type
- *  and the same master sender/receiver (just like the MPI commands).
- *
- *  A special topic are @ref GHOST_PREFETCH and @ref GHOST_PSTSTORE. For
- *  example, if all nodes broadcast to the other, the naive implementation will
- *  be that @c n_nodes times a @ref GHOST_BCST is done with different master
- *  nodes. But this means that each time <tt>n_nodes - 1</tt> nodes wait for
- *  the master to construct its send buffer. Therefore there is the prefetch
- *  flag which can be set on a pair of recv/send operations. If the ghost
- *  communication reaches a recv operation with prefetch, the next send
- *  operation (which must have the prefetch set!!) is searched and the send
- *  buffer already created. When sending, this precreated send buffer is used.
- *  In the scenario above, all nodes create the send buffers simultaneously in
- *  the first communication step, thereby reducing the latency a little bit.
- *  The pststore is similar and postpones the write back of received data
- *  until a send operation (with a precreated send buffer) is finished.
- *
- *  The ghost communicators are created by the cell systems.
+ *  In this file you find the payload-selector flags used by the
+ *  ghost-communication engine (@ref GhostComm::halo_exchange).
  */
 
 #include <config/config.hpp>
 
-#include "BoxGeometry.hpp"
-#include "ParticleList.hpp"
-
-#include <utils/Vector.hpp>
-
-#include <boost/mpi/communicator.hpp>
-
-#include <cstddef>
-#include <utility>
-#include <vector>
-
-/** \name Transfer types, for \ref GhostCommunicator::type */
-/************************************************************/
-/**@{*/
-
-/// send to a single node
-#define GHOST_SEND 0
-/// recv from a single node
-#define GHOST_RECV 1
-/// broadcast, the node entry gives the sender
-#define GHOST_BCST 2
-/// reduce, the node entry gives the receiver
-#define GHOST_RDCE 3
-/// transfer data from cell to cell on this node
-#define GHOST_LOCL 4
-
-/// mask to the job area of the transfer type
-#define GHOST_JOBMASK 15
-/// additional flag for prefetching
-#define GHOST_PREFETCH 16
-/// additional flag for poststoring
-#define GHOST_PSTSTORE 32
-/**@}*/
-
-/** Transfer data classes, for \ref ghost_communicator */
+/** Transfer data classes, for the ghost-communication engine */
 enum : unsigned {
   GHOSTTRANS_NONE = 0u,
   /// transfer \ref ParticleProperties
@@ -138,38 +47,12 @@ enum : unsigned {
 #endif
   /// resize the receiver particle arrays to the size of the senders
   GHOSTTRANS_PARTNUM = 64u,
-  GHOSTTRANS_BONDS = 128u
-};
-
-struct GhostCommunication {
-  /** Communication type. */
-  int type;
-  /** Node to communicate with (to use with all MPI operations). */
-  int node;
-
-  /** Pointer array to particle lists to communicate. */
-  std::vector<ParticleList *> part_lists = {};
-
-  /** Position shift for ghost particles. The shift is done on the sender side.
-   */
-  Utils::Vector3d shift = {};
-};
-
-/** Properties for a ghost communication. */
-struct GhostCommunicator {
-  GhostCommunicator() = default;
-  GhostCommunicator(boost::mpi::communicator comm, std::size_t size)
-      : mpi_comm(std::move(comm)), communications(size) {}
-
-  /** Attached mpi communicator */
-  boost::mpi::communicator mpi_comm;
-
-  /** List of ghost communications. */
-  std::vector<GhostCommunication> communications;
+  GHOSTTRANS_BONDS = 128u,
+#ifdef ESPRESSO_ROTATION
+  /// transfer orientation quaternion (pushed with position;
+  /// runtime-conditional)
+  GHOSTTRANS_QUAT = 256u,
+  /// transfer torque (reduced with force; runtime-conditional)
+  GHOSTTRANS_TORQUE = 512u,
+#endif
 };
-
-/**
- * @brief Do a ghost communication with the specified data parts.
- */
-void ghost_communicator(GhostCommunicator const &gcr,
-                        BoxGeometry const &box_geo, unsigned int data_parts);
diff --git a/src/core/ghosts/HaloExchange.cpp b/src/core/ghosts/HaloExchange.cpp
new file mode 100644
index 0000000000..68b4783fed
--- /dev/null
+++ b/src/core/ghosts/HaloExchange.cpp
@@ -0,0 +1,642 @@
+/*
+ * Copyright (C) 2010-2026 The ESPResSo project
+ * Copyright (C) 2002,2003,2004,2005,2006,2007,2008,2009,2010
+ *   Max-Planck-Institute for Polymer Research, Theory Group
+ *
+ * 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/>.
+ */
+/** \file
+ *  Asynchronous, split-phase ghost-communication engine. See HaloExchange.hpp.
+ */
+
+#include "ghosts/HaloExchange.hpp"
+
+#ifdef ESPRESSO_CALIPER
+#include "caliper_utils.hpp"
+#endif
+
+#include "BoxGeometry.hpp"
+#include "ghosts.hpp"
+#include "ghosts/HaloPlan.hpp"
+#include "ghosts/HaloPlanValidator.hpp"
+#include "ghosts/particle_packing.hpp"
+
+#include <boost/mpi/collectives.hpp>
+#include <boost/mpi/nonblocking.hpp>
+
+#include <algorithm>
+#include <array>
+#include <cassert>
+#include <cstddef>
+#include <functional>
+#include <memory>
+#include <span>
+#include <vector>
+
+namespace GhostComm {
+
+namespace {
+
+/**
+ * @brief MPI tags for ghost communications.
+ *
+ * One tag per exchange-kind keeps concurrent exchanges (a future goal) from
+ * aliasing: within a single exchange each message is already uniquely matched
+ * by @c (peer, tag). The BONDS transfer needs its own tag because it is sent
+ * as a *second* message to the same peer, mirroring the legacy two-transfer.
+ */
+enum GhostTag : int {
+  TAG_POSITION = 100,
+  TAG_FORCE = 101,
+  TAG_PARTNUM = 102,
+  TAG_DEFAULT = 103,
+  TAG_BONDS = 104,
+};
+
+int main_tag(unsigned data_parts) {
+  if (data_parts & GHOSTTRANS_PARTNUM)
+    return TAG_PARTNUM;
+  if (data_parts & GHOSTTRANS_FORCE)
+    return TAG_FORCE;
+  if (data_parts & GHOSTTRANS_POSITION)
+    return TAG_POSITION;
+  return TAG_DEFAULT;
+}
+// NOTE: the (peer, tag) pair is what makes each MPI message unambiguous.
+// This scheme is safe only when each peer appears at most once in the plan
+// (the "one NeighborComm per peer" invariant). See halo_exchange_start @pre.
+
+/** @brief View a list of cell pointers as a span for the packing routines. */
+std::span<ParticleList *const> as_span(std::vector<ParticleList *> const &v) {
+  return {v.data(), v.size()};
+}
+
+/** @brief Extract the plain cell pointers from a list of send regions. */
+std::vector<ParticleList *> region_cells(std::vector<SendRegion> const &send) {
+  std::vector<ParticleList *> cells;
+  cells.reserve(send.size());
+  for (auto const &r : send)
+    cells.push_back(r.cell);
+  return cells;
+}
+
+/**
+ * @brief Pack all send regions of one NeighborComm into a single buffer.
+ *
+ * All regions are packed with a **single** @c pack_cells call so that bonds
+ * produce exactly one boost binary archive per neighbor message -- matching
+ * the single @c unpack_cells call on the receive side.  Concatenating
+ * per-region archives would corrupt the bond stream (the second archive
+ * header is misread as bond data by the receiver).
+ *
+ * Per-region shifts are honored via the common shift. All @c SendRegion.shift
+ * values within a @c NeighborComm are equal — @c RegularDecomposition::
+ * make_halo_plan() sets them all to @c {}. This lets us pack every region's
+ * cells in one @c pack_cells call (a single bond archive, matching the
+ * receiver's single @c unpack_cells).
+ *
+ * @note If a future decomposition (e.g. a Lees-Edwards refactor) ever stores
+ * DISTINCT per-region shifts, @c pack_regions must be generalized to pack
+ * each region's non-bond (flat) data with its own shift while still writing
+ * all bonds into ONE shared archive — do not simply concatenate per-region
+ * archives (that was the bond-corruption bug fixed in commit 314e7c366b).
+ */
+void pack_regions(CommBuf &buf, std::vector<SendRegion> const &regions,
+                  BoxGeometry const &box, unsigned data_parts) {
+#ifdef ESPRESSO_ADDITIONAL_CHECKS
+  // All regions in a NeighborComm must share the same shift.  If per-region
+  // shifts ever diverge the packer must be generalized (separate archives).
+  if (!regions.empty()) {
+    auto const &ref = regions.front().shift;
+    for (auto const &r : regions) {
+      assert(r.shift == ref &&
+             "pack_regions: per-region shifts differ within one NeighborComm "
+             "-- packer must be generalized for per-region shift support");
+    }
+  }
+#endif
+
+  Utils::Vector3d const common_shift =
+      regions.empty() ? Utils::Vector3d{} : regions.front().shift;
+
+  auto const cells = region_cells(regions);
+  pack_cells(buf, as_span(cells), common_shift, box, data_parts);
+}
+
+/**
+ * @brief Collective (broadcast/reduce-sum) section.
+ *
+ * Mirrors the legacy GHOST_BCST / GHOST_RDCE loop from ghost_communicator()
+ * in ghosts.cpp.  For each root rank in [0, comm.size()):
+ *
+ *   Broadcast (Push): root packs its owned cell (cells[root]) and broadcasts
+ *   the buffer to all ranks; every non-root rank unpacks into cells[root]
+ *   (their ghost copy of root's data).
+ *
+ *   ReduceSum (Reduce): every rank packs cells[root] (ghost-force contribution)
+ *   and reduces to root with std::plus on the raw double buffer; root unpacks
+ *   into cells[root] (the owned cell, so forces accumulate there).
+ *
+ * The legacy code uses one GHOST_BCST (or GHOST_RDCE) step per rank, each
+ * step addressing a single cell pointer (ghost_comm.part_lists[0] ==
+ * &cells.at(n)).  We replicate that loop exactly.
+ */
+void run_collective(HaloPlan const &plan, BoxGeometry const &box,
+                    unsigned data_parts, ExchangeOp op) {
+  // Precondition: for non-PARTNUM data parts, ghost cell sizes must already be
+  // synced by a prior GHOSTTRANS_PARTNUM exchange (same invariant as the legacy
+  // GHOST_BCST/GHOST_RDCE path). The per-root broadcast/reduce byte count is
+  // derived from the (already-sized) cells; a size mismatch across ranks would
+  // be undefined behavior.
+
+  if (!plan.collective || plan.collective->pattern == CollectivePattern::None) {
+    return;
+  }
+
+  auto const &cs = *plan.collective;
+  auto const &comm = plan.comm;
+  int const comm_size = comm.size();
+  int const my_rank = comm.rank();
+
+  // Use op.direction to select the actual MPI collective: Push -> broadcast
+  // (particles flow from owner to all ghosts); Reduce -> reduce-sum (ghost
+  // forces flow back to owner).  cs.pattern == Broadcast marks the section as
+  // active; the engine caller sets op correctly for each exchange.
+  bool const is_broadcast = (op.direction == Direction::Push);
+
+  // cs.cells has one entry per rank: cs.cells[root] is the ParticleList for
+  // that root rank (owned on root, ghost on all others).
+  assert(static_cast<int>(cs.cells.size()) == comm_size);
+
+  CommBuf buf;
+  // boost::mpi's pointer-based collectives bind a reference to the buffer even
+  // for count == 0; an empty CommBuf yields data() == nullptr, which is
+  // undefined behaviour (flagged by UBSan).  Substitute a dummy address for
+  // empty buffers -- the count of 0 guarantees it is never dereferenced.
+  alignas(double) static char empty_sentinel[sizeof(double)];
+  auto const safe_data = [](CommBuf &b) {
+    return b.size() != 0u ? b.data() : &empty_sentinel[0];
+  };
+
+  for (int root = 0; root < comm_size; ++root) {
+    ParticleList *cell = cs.cells[static_cast<std::size_t>(root)];
+    auto const cell_span = std::span<ParticleList *const>{&cell, 1};
+
+    if (is_broadcast) {
+      // Push: root broadcasts its owned particles to all other ranks.
+      if (my_rank == root) {
+        pack_cells(buf, cell_span, {}, box, data_parts);
+        boost::mpi::broadcast(comm, safe_data(buf),
+                              static_cast<int>(buf.size()), root);
+        boost::mpi::broadcast(comm, buf.bonds(), root);
+      } else {
+        buf.resize(calc_transmit_size(cell_span, box, data_parts));
+        buf.bonds().clear();
+        boost::mpi::broadcast(comm, safe_data(buf),
+                              static_cast<int>(buf.size()), root);
+        boost::mpi::broadcast(comm, buf.bonds(), root);
+        unpack_cells(buf, cell_span, box, data_parts);
+      }
+    } else {
+      // ReduceSum: every rank sends its ghost-force contribution; root
+      // accumulates into its owned cell.  Mirrors the legacy GHOST_RDCE
+      // which reduces the raw double buffer with std::plus<double>.
+      pack_cells(buf, cell_span, {}, box, data_parts);
+      auto *raw = reinterpret_cast<double *>(safe_data(buf));
+      int const count = static_cast<int>(buf.size() / sizeof(double));
+      if (my_rank == root) {
+        CommBuf recv_buf;
+        recv_buf.resize(buf.size());
+        auto *recv_raw = reinterpret_cast<double *>(safe_data(recv_buf));
+        boost::mpi::reduce(comm, raw, count, recv_raw, std::plus<double>{},
+                           root);
+        unpack_cells(recv_buf, cell_span, box, data_parts);
+      } else {
+        boost::mpi::reduce(comm, raw, count, std::plus<double>{}, root);
+      }
+    }
+  }
+}
+
+} // namespace
+
+GhostExchange halo_exchange_start(HaloPlan const &plan, BoxGeometry const &box,
+                                  unsigned data_parts, ExchangeOp op,
+                                  ExchangeBuffers &bufs) {
+  GhostExchange st;
+  st.op = op;
+  st.data_parts = data_parts;
+  st.box = &box;
+  st.plan = &plan;
+  st.bufs = &bufs;
+  if (data_parts == GHOSTTRANS_NONE)
+    return st;
+
+#ifdef ESPRESSO_ADDITIONAL_CHECKS
+  // Cross-rank symmetry check: run once per plan at first use.
+  //
+  // This check is deferred from the decomposition constructors to here because:
+  //   (a) During checkpoint loading, decompositions are transiently rebuilt
+  //       while maximal_cutoff is rank-divergent (different cell grids per rank
+  //       for a brief window before the next consistent rebuild).  The
+  //       transient plan is never used, so its asymmetry is harmless — but a
+  //       ctor-time collective would fire and abort.
+  //   (b) A collective inside a ctor risks partial-abort deadlock: if one rank
+  //       throws before entering the all_to_all, the others spin forever.
+  //
+  // All ranks enter halo_exchange_start together (ghost exchange is
+  // collective), so the all_to_all inside validate_halo_plan_symmetry is safe
+  // here. plan.symmetry_validated is mutable and rearmed to false whenever the
+  // plan object is rebuilt (every decomposition change constructs a fresh
+  // HaloPlan).
+  if (!plan.symmetry_validated) {
+    assert(GhostComm::report_violations(
+        GhostComm::validate_halo_plan_symmetry(plan),
+        "halo_exchange first use"));
+    plan.symmetry_validated = true;
+  }
+  // Peer-uniqueness is validated at plan-build time by validate_halo_plan().
+  // Op-sanity: Combine::Add is only valid for reducible parts
+  // (FORCE, FORCE|TORQUE, RATTLE).  TORQUE may only appear together with
+  // FORCE (it is never sent alone).
+  // Compute the predicate outside assert() so the #ifdef is not embedded in
+  // macro arguments (non-portable under -Werror).
+  [[maybe_unused]] bool const force_only = (data_parts == GHOSTTRANS_FORCE);
+#ifdef ESPRESSO_ROTATION
+  [[maybe_unused]] bool const force_with_torque =
+      (data_parts == (GHOSTTRANS_FORCE | GHOSTTRANS_TORQUE));
+#else
+  [[maybe_unused]] bool const force_with_torque = false;
+#endif
+  [[maybe_unused]] bool reducible = (force_only || force_with_torque);
+#ifdef ESPRESSO_BOND_CONSTRAINT
+  reducible = reducible || (data_parts == GHOSTTRANS_RATTLE);
+#endif
+  assert((op.combine != Combine::Add || reducible) &&
+         "Combine::Add only valid for reducible parts (FORCE, FORCE|TORQUE, "
+         "RATTLE)");
+#endif
+
+  auto const &comm = plan.comm;
+  auto const n = plan.neighbors.size();
+
+  // Resize the pool to the current neighbor count.  Existing entries retain
+  // their vector capacity so that, after the first call (warm-up), CommBuf
+  // ::resize on a same-sized or smaller payload is a no-op on the allocator.
+  bufs.send.resize(n);
+  bufs.recv.resize(n);
+  bufs.send_cells.resize(n);
+  bufs.recv_cells.resize(n);
+  bufs.slot_to_neighbor.resize(n);
+  // Clear and reuse the request vector (capacity retained).
+  bufs.requests.clear();
+
+  bool const push = op.direction == Direction::Push;
+  bool const bonds = (data_parts & GHOSTTRANS_BONDS) != 0u;
+  int const tag = main_tag(data_parts);
+
+  // Resolve, per neighbor, which cells we send from and which we receive into.
+  // Push:   send from the region cells, receive into the ghost (recv) cells.
+  // Reduce: roles swap -- send from ghost cells, receive into region cells,
+  //         adding on arrival.
+  for (std::size_t i = 0; i < n; ++i) {
+    auto const &nc = plan.neighbors[i];
+    if (push) {
+      bufs.send_cells[i] = region_cells(nc.send);
+      bufs.recv_cells[i] = nc.recv;
+    } else {
+      bufs.send_cells[i] = nc.recv;
+      bufs.recv_cells[i] = region_cells(nc.send);
+    }
+  }
+
+  // 1) Pack all send buffers.  Done before posting receives so that sizes are
+  //    known; also isolates the serialisation work in the ghost/pack region.
+#ifdef ESPRESSO_CALIPER
+  if (espresso_cali_active())
+    CALI_MARK_BEGIN("ghost/pack");
+#endif
+  for (std::size_t i = 0; i < n; ++i) {
+    auto const &nc = plan.neighbors[i];
+    if (push) {
+      pack_regions(bufs.send[i], nc.send, box, data_parts);
+    } else {
+      // Reduce: pack plain ghost cells, no shift.
+      pack_cells(bufs.send[i], as_span(bufs.send_cells[i]), {}, box,
+                 data_parts);
+    }
+  }
+#ifdef ESPRESSO_CALIPER
+  if (espresso_cali_active())
+    CALI_MARK_END("ghost/pack");
+#endif
+
+  // 2) Post ALL receives first (deadlock-free), then ALL sends.  Sizes are
+  //    known a-priori from the packed size of the cells we will receive into,
+  //    except PARTNUM whose ghost cells are not sized yet -- there we post a
+  //    fixed-size recv and let unpack_cells resize the cells (legacy bootstrap
+  //    path).
+#ifdef ESPRESSO_CALIPER
+  if (espresso_cali_active())
+    CALI_MARK_BEGIN("ghost/post");
+#endif
+  for (std::size_t i = 0; i < n; ++i) {
+    auto const &nc = plan.neighbors[i];
+    bufs.recv[i].resize(
+        calc_transmit_size(as_span(bufs.recv_cells[i]), box, data_parts));
+    bufs.requests.push_back(comm.irecv(nc.peer, tag, bufs.recv[i].data(),
+                                       static_cast<int>(bufs.recv[i].size())));
+  }
+  for (std::size_t i = 0; i < n; ++i) {
+    auto const &nc = plan.neighbors[i];
+    bufs.requests.push_back(comm.isend(nc.peer, tag, bufs.send[i].data(),
+                                       static_cast<int>(bufs.send[i].size())));
+  }
+
+  // BONDS (cold, resort-only path): a second per-neighbor message for the bond
+  // buffers, mirroring the legacy two-transfer. Kept entirely off the hot
+  // POSITION/FORCE path via the flag guard. boost::mpi serializes the
+  // std::vector<char> length + payload, so recv sizing is automatic here.
+  if (bonds) {
+    for (std::size_t i = 0; i < n; ++i) {
+      auto const &nc = plan.neighbors[i];
+      bufs.requests.push_back(
+          comm.irecv(nc.peer, TAG_BONDS, bufs.recv[i].bonds()));
+    }
+    for (std::size_t i = 0; i < n; ++i) {
+      auto const &nc = plan.neighbors[i];
+      bufs.requests.push_back(
+          comm.isend(nc.peer, TAG_BONDS, bufs.send[i].bonds()));
+    }
+  }
+#ifdef ESPRESSO_CALIPER
+  if (espresso_cali_active())
+    CALI_MARK_END("ghost/post");
+#endif
+
+  // Local (same-rank) copies run here, after messages are posted, so they
+  // overlap with in-flight network transfers for split-phase callers.
+  //
+  // Direction:
+  //   Push:   copy the real (src) cell into its self-ghost (dst), applying the
+  //           periodic-image shift so the ghost holds the correct replica.
+  //   Reduce: roles swap -- the self-ghost (dst) accumulated a short-range
+  //           force contribution during the pair loop that must be added back
+  //           into the real (src) cell.  Copy dst -> src (no position shift;
+  //           GHOSTTRANS_FORCE with ReductionPolicy::UPDATE accumulates via
+  //           +=).  This mirrors the legacy GHOST_LOCL entry in the reverted
+  //           collect-forces communicator; without the swap the self-wrap
+  //           force contributions (node_grid[i]==1 periodic axes) are dropped,
+  //           which silently biases forces on particles that interact across
+  //           the local periodic boundary.
+  //
+  // Correctness (no aliasing with in-flight buffers):
+  //   Push:   local_cell_copy writes ghost cells; isends carry real-cell data.
+  //   Reduce: local_cell_copy writes real (owned) cells; isends carry ghost-
+  //           cell data.  In both cases the written and sent cell sets are
+  //           disjoint.  Local copies complete before halo_exchange_finish
+  //           unpacks any receive buffer.
+  for (auto const &lc : plan.local) {
+    if (push)
+      local_cell_copy(*lc.src, *lc.dst, lc.shift, box, data_parts);
+    else
+      local_cell_copy(*lc.dst, *lc.src, Utils::Vector3d{}, box, data_parts);
+  }
+
+  return st;
+}
+
+GhostExchange halo_exchange_start(HaloPlan const &plan, BoxGeometry const &box,
+                                  unsigned data_parts, ExchangeOp op) {
+  // Convenience overload for callers that do not hold a persistent
+  // ExchangeBuffers (unit tests, cold resort paths).  Allocate a pool via
+  // unique_ptr so it is freed automatically when the GhostExchange is
+  // destroyed — on every exit path, including the GHOSTTRANS_NONE early return
+  // in halo_exchange_finish.  Use the pool overload on hot paths.
+  auto pool = std::make_unique<ExchangeBuffers>();
+  auto st = halo_exchange_start(plan, box, data_parts, op, *pool);
+  // Transfer ownership: owned keeps the pool alive; bufs already points to it.
+  st.owned = std::move(pool);
+  return st;
+}
+
+void halo_exchange_finish(GhostExchange &st) {
+  if (st.data_parts == GHOSTTRANS_NONE)
+    return;
+
+  auto &bufs = *st.bufs;
+
+  // Local (same-rank) copies were already performed in halo_exchange_start
+  // to overlap with in-flight messages; nothing to do here.
+
+  if (st.plan->collective)
+    run_collective(*st.plan, *st.box, st.data_parts, st.op);
+
+  auto const n = st.plan->neighbors.size();
+  bool const bonds = (st.data_parts & GHOSTTRANS_BONDS) != 0u;
+  bool const add = st.op.combine == Combine::Add;
+
+  // Request layout (set by halo_exchange_start):
+  //   [0..n)    = per-neighbor irecvs
+  //   [n..2n)   = per-neighbor isends
+  //   [2n..3n)  = bond irecvs   (only when bonds)
+  //   [3n..4n)  = bond isends   (only when bonds)
+
+  if (bonds) {
+    // Cold path (resort-only): bonds need *both* the flat and bond message from
+    // each neighbor before unpack_cells can run.  Keep existing wait_all-then-
+    // unpack-in-order semantics; pipelining is not safe here.
+#ifdef ESPRESSO_CALIPER
+    if (espresso_cali_active())
+      CALI_MARK_BEGIN("ghost/wait");
+#endif
+    boost::mpi::wait_all(bufs.requests.begin(), bufs.requests.end());
+#ifdef ESPRESSO_CALIPER
+    if (espresso_cali_active())
+      CALI_MARK_END("ghost/wait");
+    if (espresso_cali_active())
+      CALI_MARK_BEGIN("ghost/unpack");
+#endif
+    for (std::size_t i = 0; i < n; ++i) {
+      auto dst = as_span(bufs.recv_cells[i]);
+      unpack_cells(bufs.recv[i], dst, *st.box, st.data_parts);
+    }
+#ifdef ESPRESSO_CALIPER
+    if (espresso_cali_active())
+      CALI_MARK_END("ghost/unpack");
+#endif
+  } else if (add) {
+    // Add path (FORCE / RATTLE reduce): wait receives in *fixed neighbor order*
+    // so that floating-point addition into shared local cells is bitwise
+    // reproducible across runs.  Unpacking neighbor i overlaps delivery of
+    // i+1..n-1 (pipelining while preserving associativity).
+    //
+    // Caliper: outer ghost/wait and ghost/unpack markers always fire (even
+    // when n==0) so that label discovery works on single-rank runs.
+    // Per-neighbor inner BEGIN/END pairs interleave with the outer ones;
+    // repeated sequential begin/end accumulates — valid Caliper usage.
+#ifdef ESPRESSO_CALIPER
+    if (espresso_cali_active())
+      CALI_MARK_BEGIN("ghost/wait");
+#endif
+    for (std::size_t i = 0; i < n; ++i) {
+      bufs.requests[i].wait();
+#ifdef ESPRESSO_CALIPER
+      if (espresso_cali_active())
+        CALI_MARK_END("ghost/wait");
+      if (espresso_cali_active())
+        CALI_MARK_BEGIN("ghost/unpack");
+#endif
+      auto dst = as_span(bufs.recv_cells[i]);
+      if (st.data_parts & GHOSTTRANS_FORCE) {
+        add_forces(bufs.recv[i], dst, st.data_parts);
+      }
+#ifdef ESPRESSO_BOND_CONSTRAINT
+      else if (st.data_parts == GHOSTTRANS_RATTLE) {
+        add_rattle(bufs.recv[i], dst);
+      }
+#endif
+#ifdef ESPRESSO_CALIPER
+      if (espresso_cali_active())
+        CALI_MARK_END("ghost/unpack");
+      if (espresso_cali_active())
+        CALI_MARK_BEGIN("ghost/wait");
+#endif
+    }
+    // Wait for sends [n..2n) so buffers are safe to reuse next step.
+    // The outer ghost/wait region (opened before the loop) is still open here
+    // on the n==0 path; on n>0 the loop re-opened it after the last unpack.
+    boost::mpi::wait_all(bufs.requests.begin() + static_cast<std::ptrdiff_t>(n),
+                         bufs.requests.begin() +
+                             static_cast<std::ptrdiff_t>(2 * n));
+#ifdef ESPRESSO_CALIPER
+    if (espresso_cali_active())
+      CALI_MARK_END("ghost/wait");
+    if (espresso_cali_active())
+      CALI_MARK_BEGIN("ghost/unpack");
+    if (espresso_cali_active())
+      CALI_MARK_END("ghost/unpack");
+#endif
+  } else {
+    // Overwrite path (position/property push): arrival order is safe because
+    // the validator guarantees each ghost cell is filled by exactly one
+    // message. Use wait_any to unpack each neighbor's buffer as soon as it
+    // arrives.
+    //
+    // Bookkeeping: maintain a slot→neighbor index map so we always know which
+    // recv buffer to unpack for the completed request slot.  Completed requests
+    // are swapped to the end of the active range and the map is updated in
+    // sync.
+    //
+    // Caliper: outer ghost/wait always fires (even n==0); the wait_any loop
+    // emits nested end/begin pairs; an outer ghost/unpack pair fires after the
+    // loop so the label appears on single-rank runs.
+#ifdef ESPRESSO_CALIPER
+    if (espresso_cali_active())
+      CALI_MARK_BEGIN("ghost/wait");
+#endif
+    // Initialize the scratch map to the identity: slot i → neighbor i.
+    // The vector was pre-sized in halo_exchange_start (no allocation here).
+    for (std::size_t i = 0; i < n; ++i)
+      bufs.slot_to_neighbor[i] = i;
+
+    // active_end tracks how many recv requests are still pending.
+    std::size_t active_end = n;
+
+    while (active_end > 0) {
+      auto first = bufs.requests.begin();
+      auto last =
+          bufs.requests.begin() + static_cast<std::ptrdiff_t>(active_end);
+
+      auto [status, done_it] = boost::mpi::wait_any(first, last);
+      (void)status;
+#ifdef ESPRESSO_CALIPER
+      if (espresso_cali_active())
+        CALI_MARK_END("ghost/wait");
+      if (espresso_cali_active())
+        CALI_MARK_BEGIN("ghost/unpack");
+#endif
+
+      std::size_t const done_slot =
+          static_cast<std::size_t>(done_it - bufs.requests.begin());
+      std::size_t const neighbor_idx = bufs.slot_to_neighbor[done_slot];
+
+      {
+        auto dst = as_span(bufs.recv_cells[neighbor_idx]);
+        unpack_cells(bufs.recv[neighbor_idx], dst, *st.box, st.data_parts);
+      }
+
+      // Remove the completed slot from the active range by swapping it with
+      // the last active slot and shrinking the range.
+      --active_end;
+      if (done_slot != active_end) {
+        std::iter_swap(done_it, bufs.requests.begin() +
+                                    static_cast<std::ptrdiff_t>(active_end));
+        std::swap(bufs.slot_to_neighbor[done_slot],
+                  bufs.slot_to_neighbor[active_end]);
+      }
+#ifdef ESPRESSO_CALIPER
+      if (espresso_cali_active())
+        CALI_MARK_END("ghost/unpack");
+      if (espresso_cali_active() && active_end > 0)
+        CALI_MARK_BEGIN("ghost/wait");
+#endif
+    }
+
+    // Wait for sends [n..2n) so buffers are safe to reuse next step.
+    // Caliper balance: when n==0 the outer CALI_MARK_BEGIN("ghost/wait") posted
+    // before the loop is still open (the loop body never ran, so nothing closed
+    // it).  When n>0 the last loop iteration emitted END("ghost/unpack") and
+    // then — because active_end reached 0 — did NOT re-open "ghost/wait", so
+    // the region is closed.  The guard `if (n > 0)` therefore re-opens for
+    // n>0 and is deliberately absent for n==0 (the outer region serves).
+    // Invariant: exactly one "ghost/wait" region is open entering wait_all.
+#ifdef ESPRESSO_CALIPER
+    if (espresso_cali_active() && n > 0)
+      CALI_MARK_BEGIN("ghost/wait");
+#endif
+    boost::mpi::wait_all(bufs.requests.begin() + static_cast<std::ptrdiff_t>(n),
+                         bufs.requests.begin() +
+                             static_cast<std::ptrdiff_t>(2 * n));
+#ifdef ESPRESSO_CALIPER
+    if (espresso_cali_active())
+      CALI_MARK_END("ghost/wait");
+    if (espresso_cali_active())
+      CALI_MARK_BEGIN("ghost/unpack");
+    if (espresso_cali_active())
+      CALI_MARK_END("ghost/unpack");
+#endif
+  }
+
+  // st.owned (if non-null) will free the heap-allocated pool automatically
+  // when the GhostExchange goes out of scope after this call returns.
+}
+
+void halo_exchange(HaloPlan const &plan, BoxGeometry const &box,
+                   unsigned data_parts, ExchangeOp op, ExchangeBuffers &bufs) {
+  auto st = halo_exchange_start(plan, box, data_parts, op, bufs);
+  halo_exchange_finish(st);
+}
+
+void halo_exchange(HaloPlan const &plan, BoxGeometry const &box,
+                   unsigned data_parts, ExchangeOp op) {
+  // Convenience overload: uses the no-pool start, which allocates a temporary
+  // ExchangeBuffers on the heap and frees it in finish.
+  auto st = halo_exchange_start(plan, box, data_parts, op);
+  halo_exchange_finish(st);
+}
+
+} // namespace GhostComm
diff --git a/src/core/ghosts/HaloExchange.hpp b/src/core/ghosts/HaloExchange.hpp
new file mode 100644
index 0000000000..3d0258b627
--- /dev/null
+++ b/src/core/ghosts/HaloExchange.hpp
@@ -0,0 +1,179 @@
+/*
+ * Copyright (C) 2010-2026 The ESPResSo project
+ * Copyright (C) 2002,2003,2004,2005,2006,2007,2008,2009,2010
+ *   Max-Planck-Institute for Polymer Research, Theory Group
+ *
+ * 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
+ *  Asynchronous, split-phase ghost-communication engine.
+ *
+ *  Drives a @ref GhostComm::HaloPlan with non-blocking point-to-point MPI
+ *  (plus same-rank copies), reusing the byte-identical serialization from
+ *  particle_packing.hpp. The engine is split into a @c halo_exchange_start
+ *  phase (post all receives, pack and post all sends) and a
+ *  @c halo_exchange_finish phase (same-rank copies, wait, unpack/reduce) so
+ *  that callers can overlap other work with in-flight messages.
+ *
+ *  Deadlock-freedom is guaranteed by posting every @c irecv before any
+ *  @c isend and matching each message by @c (peer, tag); message sizes are
+ *  known a-priori from the receiving cells' packed size (see
+ *  @ref GhostComm::calc_transmit_size), except for the PARTNUM bootstrap.
+ */
+
+#include "BoxGeometry.hpp"
+#include "ghosts/HaloPlan.hpp"
+#include "ghosts/particle_packing.hpp"
+
+#include <boost/mpi/request.hpp>
+
+#include <memory>
+#include <vector>
+
+namespace GhostComm {
+
+/**
+ * @brief Persistent per-neighbor buffer pool for halo exchanges.
+ *
+ * Owns the heap-allocated send/recv buffers, MPI-request slots, and
+ * cell-pointer scratch arrays for one set of neighbor communications.
+ * Holding this object across multiple calls to @c halo_exchange allows
+ * the underlying @c std::vector storage to be reused via @c resize —
+ * after the first exchange (warm-up) the vectors retain capacity and
+ * subsequent calls incur no heap allocations on the hot path.
+ *
+ * Typical ownership: held as a member of the object that drives ghost
+ * communication (e.g. @c CellStructure::m_ghost_buffers) and passed to
+ * @c halo_exchange / @c halo_exchange_start as a mutable ref.
+ *
+ * Thread-safety: a single instance must not be used from concurrent threads.
+ * The sequential ghost-exchange loop (each exchange completes via
+ * @c halo_exchange_finish before the next starts) is therefore safe.
+ */
+struct ExchangeBuffers {
+  /** Per-neighbor packed send buffers (index-aligned with plan->neighbors). */
+  std::vector<CommBuf> send;
+  /** Per-neighbor recv buffers (index-aligned with plan->neighbors). */
+  std::vector<CommBuf> recv;
+  /** Outstanding non-blocking send/recv requests (cleared before each use). */
+  std::vector<boost::mpi::request> requests;
+  /**
+   * Scratch cell-pointer arrays for the packing routines.  For the Reduce
+   * direction send/recv roles swap; we need the plain @c ParticleList* from
+   * the @c SendRegion list.  These must outlive the pack/unpack calls.
+   */
+  std::vector<std::vector<ParticleList *>> send_cells;
+  std::vector<std::vector<ParticleList *>> recv_cells;
+  /**
+   * Scratch index map for the Overwrite (wait_any) path in
+   * @c halo_exchange_finish: maps active request slot → original neighbor
+   * index.  Sized to the current neighbor count in @c halo_exchange_start
+   * alongside the other per-neighbor vectors; capacity is retained across
+   * calls so that no heap allocation occurs on the hot position-push path
+   * after the first (warm-up) exchange.
+   */
+  std::vector<std::size_t> slot_to_neighbor;
+};
+
+/**
+ * @brief Opaque handle for one in-flight halo exchange.
+ *
+ * Records the per-call metadata (op, data parts, geometry, plan) and holds a
+ * pointer to the @ref ExchangeBuffers that back the in-flight messages.  When
+ * a persistent @ref ExchangeBuffers is supplied by the caller (via the pool
+ * overload), @c owned is null and @c bufs points to the caller's pool (no
+ * ownership).  When the no-pool convenience overload is used, @c owned holds
+ * the heap-allocated pool and @c bufs == @c owned.get(); the pool is freed
+ * automatically when the handle is destroyed — on every exit path, including
+ * the @c GHOSTTRANS_NONE early return in @c halo_exchange_finish.
+ *
+ * Created by @c halo_exchange_start and consumed exactly once by
+ * @c halo_exchange_finish.
+ */
+struct GhostExchange {
+  ExchangeOp op{};
+  unsigned data_parts = 0u;
+  BoxGeometry const *box = nullptr;
+  HaloPlan const *plan = nullptr;
+  /**
+   * Non-null only when the no-pool overload allocated the buffer pool.
+   * Destroyed (and the pool freed) when the @c GhostExchange goes out of
+   * scope, guaranteeing cleanup on every exit path.
+   */
+  std::unique_ptr<ExchangeBuffers> owned;
+  /** Non-owning pointer to the active buffer pool (caller's or @c owned). */
+  ExchangeBuffers *bufs = nullptr;
+};
+
+/**
+ * @brief Begin a halo exchange using a caller-owned buffer pool.
+ *
+ * The buffers in @p bufs are resized to the current neighbor count (retaining
+ * capacity) so that, after the first call (warm-up), no heap allocation
+ * occurs on the POSITION / FORCE hot path.
+ *
+ * @param plan        Communication plan (peers, regions, local copies).
+ * @param box         Box geometry for position folding.
+ * @param data_parts  Bitmask of GHOSTTRANS_* flags to transfer.
+ * @param op          Direction (Push/Reduce) and combine mode (Overwrite/Add).
+ * @param bufs        Persistent buffer pool (must outlive the returned handle).
+ *
+ * @pre Each peer rank appears at most once in @p plan.neighbors. Multiple
+ *      send/receive regions to the same peer **must** be folded into that
+ *      peer's single @ref NeighborComm (as the plan builders do). This is
+ *      what makes the @c (peer, data-part tag) message matching unambiguous:
+ *      without this invariant two @c NeighborComm entries sharing the same
+ *      peer and tag would cross-match, silently corrupting data or deadlocking.
+ *      The plan builders (e.g. RegularDecomposition) always produce unique
+ *      peers; the invariant is asserted at runtime when
+ *      @c ESPRESSO_ADDITIONAL_CHECKS is defined.
+ */
+GhostExchange halo_exchange_start(HaloPlan const &plan, BoxGeometry const &box,
+                                  unsigned data_parts, ExchangeOp op,
+                                  ExchangeBuffers &bufs);
+
+/**
+ * @brief Convenience overload that allocates a temporary @ref ExchangeBuffers.
+ *
+ * Use this form when buffer reuse across calls is not needed (e.g. unit tests,
+ * cold resort paths).  Equivalent to the original single-call allocation.
+ */
+GhostExchange halo_exchange_start(HaloPlan const &plan, BoxGeometry const &box,
+                                  unsigned data_parts, ExchangeOp op);
+
+/**
+ * @brief Complete a halo exchange: run same-rank copies (overlapping the
+ *        in-flight messages), wait for all requests, then unpack or reduce.
+ */
+void halo_exchange_finish(GhostExchange &state);
+
+/**
+ * @brief Blocking wrapper using a caller-owned buffer pool (no per-call alloc
+ *        after warm-up).
+ */
+void halo_exchange(HaloPlan const &plan, BoxGeometry const &box,
+                   unsigned data_parts, ExchangeOp op, ExchangeBuffers &bufs);
+
+/**
+ * @brief Blocking convenience wrapper: start + finish (allocates a temporary
+ *        @ref ExchangeBuffers; use the pool overload on hot paths).
+ */
+void halo_exchange(HaloPlan const &plan, BoxGeometry const &box,
+                   unsigned data_parts, ExchangeOp op);
+
+} // namespace GhostComm
diff --git a/src/core/ghosts/HaloPlan.hpp b/src/core/ghosts/HaloPlan.hpp
new file mode 100644
index 0000000000..e2f780ea59
--- /dev/null
+++ b/src/core/ghosts/HaloPlan.hpp
@@ -0,0 +1,75 @@
+/*
+ * Copyright (C) 2010-2026 The ESPResSo project
+ * Copyright (C) 2002,2003,2004,2005,2006,2007,2008,2009,2010
+ *   Max-Planck-Institute for Polymer Research, Theory Group
+ *
+ * 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
+#include "ParticleList.hpp"
+#include <boost/mpi/communicator.hpp>
+#include <optional>
+#include <utils/Vector.hpp>
+#include <vector>
+
+namespace GhostComm {
+enum class Direction { Push, Reduce }; // real->ghost ; ghost->real
+enum class Combine { Overwrite, Add };
+struct ExchangeOp {
+  Direction direction;
+  Combine combine;
+};
+
+struct SendRegion {
+  ParticleList *cell;
+  Utils::Vector3d shift;
+};
+struct NeighborComm {
+  int peer;
+  std::vector<SendRegion> send;
+  std::vector<ParticleList *> recv; // recv[k] <-> peer.send[k]
+};
+struct LocalComm {
+  ParticleList *src;
+  ParticleList *dst;
+  Utils::Vector3d shift;
+};
+
+enum class CollectivePattern { None, Broadcast, ReduceSum };
+struct CollectiveSection {
+  CollectivePattern pattern;
+  std::vector<ParticleList *> cells;
+};
+
+struct HaloPlan {
+  boost::mpi::communicator comm;
+  std::vector<NeighborComm> neighbors;
+  std::vector<LocalComm> local;
+  std::optional<CollectiveSection> collective;
+  /**
+   * @brief Lazily-set flag for the once-per-plan symmetry check.
+   *
+   * Set to true on the first call to halo_exchange_start for this plan
+   * (see HaloExchange.cpp, ESPRESSO_ADDITIONAL_CHECKS guard).  The plan
+   * object is rebuilt on every decomposition change, so the flag naturally
+   * rearms: a new plan always starts with symmetry_validated = false.
+   *
+   * mutable so halo_exchange_start can set it through a const HaloPlan&
+   * without requiring a non-const overload.
+   */
+  mutable bool symmetry_validated = false;
+};
+} // namespace GhostComm
diff --git a/src/core/ghosts/HaloPlanValidator.cpp b/src/core/ghosts/HaloPlanValidator.cpp
new file mode 100644
index 0000000000..05b09780db
--- /dev/null
+++ b/src/core/ghosts/HaloPlanValidator.cpp
@@ -0,0 +1,268 @@
+/*
+ * 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/>.
+ */
+/** \file
+ *  Local (single-rank) validator for GhostComm::HaloPlan.
+ *  See HaloPlanValidator.hpp.
+ */
+
+#include "ghosts/HaloPlanValidator.hpp"
+
+#include <boost/mpi/collectives/all_to_all.hpp>
+
+#include <iostream>
+#include <sstream>
+#include <unordered_map>
+#include <unordered_set>
+
+namespace GhostComm {
+
+std::vector<std::string>
+validate_halo_plan(HaloPlan const &plan, std::span<Cell *const> local_cells,
+                   std::span<Cell *const> ghost_cells) {
+  std::vector<std::string> violations;
+
+  // Build the set of expected ghost ParticleList pointers.
+  std::unordered_set<ParticleList const *> ghost_set;
+  for (Cell *c : ghost_cells) {
+    ghost_set.insert(&c->particles());
+  }
+
+  // Collective-covered cells: AtomDecomposition (collective-only) and
+  // HybridDecomposition (neighbors + collective) fill some/all ghosts via the
+  // n-square broadcast/reduce section, NOT via neighbors/local. The section's
+  // `cells` vector identifies exactly which ParticleLists it covers (one per
+  // rank; the entry for the local rank is this rank's own cell, the rest are
+  // ghost copies of every other rank's owned cell). Ghost cells found here are
+  // covered even though they are not point-to-point recv/dst targets.
+  bool const has_collective =
+      plan.collective.has_value() &&
+      plan.collective->pattern != CollectivePattern::None;
+  std::unordered_set<ParticleList const *> collective_set;
+  if (has_collective) {
+    for (ParticleList const *pl : plan.collective->cells) {
+      collective_set.insert(pl);
+    }
+  }
+
+  // Set of ghost cells that some local cell actually interacts with, i.e. that
+  // appear in a local cell's neighbor stencil. Only these ghosts must be
+  // filled by the plan; a ghost that no local cell ever references carries no
+  // physics and needs no communication. This happens on a single MPI rank:
+  // make_halo_plan() returns an empty plan there because init_cell_interactions
+  // wires periodic neighbourships directly between local cells (see
+  // RegularDecomposition), so the halo-layer cells that mark_cells() still
+  // classifies as ghosts are never referenced and must not be flagged as
+  // uncovered. Double-fill and out-of-ghost-set targets stay strict below, and
+  // the neighborship-match check still requires every *referenced* ghost to be
+  // covered, so this cannot mask a real missing-communication defect.
+  std::unordered_set<ParticleList const *> referenced_ghosts;
+  for (Cell *c : local_cells) {
+    for (Cell *n : c->neighbors().all()) {
+      ParticleList const *pl = &n->particles();
+      if (ghost_set.find(pl) != ghost_set.end()) {
+        referenced_ghosts.insert(pl);
+      }
+    }
+  }
+
+  // Check peer-uniqueness and shape; accumulate recv/dst fill counts.
+  std::unordered_map<ParticleList const *, int> fill_count;
+  std::unordered_set<int> seen_peers;
+
+  for (auto const &nc : plan.neighbors) {
+    // Peer-uniqueness
+    if (!seen_peers.insert(nc.peer).second) {
+      std::ostringstream oss;
+      oss << "peer " << nc.peer << " appears in more than one NeighborComm";
+      violations.emplace_back(oss.str());
+    }
+
+    // Shape: send.size() == recv.size()
+    if (nc.send.size() != nc.recv.size()) {
+      std::ostringstream oss;
+      oss << "NeighborComm peer=" << nc.peer
+          << " has send.size()=" << nc.send.size()
+          << " != recv.size()=" << nc.recv.size();
+      violations.emplace_back(oss.str());
+    }
+
+    // Accumulate recv fill counts; check targets are in ghost set.
+    for (ParticleList const *pl : nc.recv) {
+      if (ghost_set.find(pl) == ghost_set.end()) {
+        std::ostringstream oss;
+        oss << "NeighborComm peer=" << nc.peer
+            << " recv target is not a ghost cell";
+        violations.emplace_back(oss.str());
+      }
+      ++fill_count[pl];
+    }
+  }
+
+  // Accumulate local.dst fill counts; check targets are in ghost set.
+  for (auto const &lc : plan.local) {
+    ParticleList const *pl = lc.dst;
+    if (ghost_set.find(pl) == ghost_set.end()) {
+      violations.emplace_back("LocalComm dst target is not a ghost cell");
+    }
+    ++fill_count[pl];
+  }
+
+  // Coverage: every ghost that a local cell references must be filled.
+  // Point-to-point ghosts must appear exactly once as a recv/dst target; a
+  // double-fill is always a defect. Ghosts covered by the collective section
+  // are filled by the broadcast/reduce instead, so a zero point-to-point
+  // fill-count is fine for them (and they do not contribute to fill_count, so
+  // they cannot trigger the double-fill path). A ghost that no local cell
+  // references (e.g. the halo-layer cells on a single MPI rank, where the plan
+  // is intentionally empty) carries no physics and is not required to be
+  // filled, so a zero fill-count is fine for it too.
+  for (Cell *c : ghost_cells) {
+    ParticleList const *pl = &c->particles();
+    auto it = fill_count.find(pl);
+    int count = (it != fill_count.end()) ? it->second : 0;
+    if (count == 0) {
+      if (collective_set.find(pl) == collective_set.end() &&
+          referenced_ghosts.find(pl) != referenced_ghosts.end()) {
+        violations.emplace_back(
+            "ghost cell is never filled (missing recv/dst)");
+      }
+    } else if (count > 1) {
+      std::ostringstream oss;
+      oss << "ghost cell is filled " << count << " times (expected 1)";
+      violations.emplace_back(oss.str());
+    }
+  }
+
+  // Neighborship-match: every local cell's ghost neighbor must be covered,
+  // either as a point-to-point recv/dst target or by the collective section.
+  for (Cell *c : local_cells) {
+    for (Cell *n : c->neighbors().all()) {
+      ParticleList const *pl = &n->particles();
+      if (ghost_set.find(pl) == ghost_set.end()) {
+        continue; // not a ghost neighbor
+      }
+      if (collective_set.find(pl) != collective_set.end()) {
+        continue; // covered by the collective broadcast/reduce section
+      }
+      auto it = fill_count.find(pl);
+      int count = (it != fill_count.end()) ? it->second : 0;
+      if (count == 0) {
+        violations.emplace_back(
+            "local cell has ghost neighbor that is not a covered recv/dst "
+            "target (referenced-but-uncommunicated ghost)");
+      }
+    }
+  }
+
+  // Interior/boundary consistency: a local cell marked interior
+  // (is_boundary()==false) must have no ghost neighbor.
+  for (Cell *c : local_cells) {
+    if (c->is_boundary()) {
+      continue; // boundary cells are expected to have ghost neighbors
+    }
+    for (Cell *n : c->neighbors().all()) {
+      if (ghost_set.find(&n->particles()) != ghost_set.end()) {
+        violations.emplace_back("interior cell has a ghost neighbor");
+        break; // one violation per cell is sufficient
+      }
+    }
+  }
+
+  // Overlap-safety invariant (Task 5.3 precondition): an interior cell must
+  // not appear as a send source in any NeighborComm, and must not appear as
+  // a LocalComm src.  Send sources are real cells whose data gets copied to a
+  // peer's ghost; LocalComm srcs are the matching self-copy sources.  If an
+  // interior cell were in either, the force-reduce step would overwrite its
+  // ghost copy after the interior velocity update — violating the overlap
+  // correctness assumption.
+  //
+  // Build a set of ParticleList* that belong to interior local cells.
+  std::unordered_set<ParticleList const *> interior_set;
+  for (Cell *c : local_cells) {
+    if (!c->is_boundary()) {
+      interior_set.insert(&c->particles());
+    }
+  }
+  if (!interior_set.empty()) {
+    for (auto const &nc : plan.neighbors) {
+      for (auto const &sr : nc.send) {
+        if (interior_set.count(sr.cell)) {
+          std::ostringstream oss;
+          oss << "interior cell appears as NeighborComm send source (peer="
+              << nc.peer << ") — overlap-safety invariant violated";
+          violations.emplace_back(oss.str());
+        }
+      }
+    }
+    for (auto const &lc : plan.local) {
+      if (interior_set.count(lc.src)) {
+        violations.emplace_back(
+            "interior cell appears as LocalComm src — overlap-safety invariant "
+            "violated");
+      }
+    }
+  }
+
+  return violations;
+}
+
+std::vector<std::string> validate_halo_plan_symmetry(HaloPlan const &plan) {
+  std::vector<std::string> violations;
+
+  int const n = plan.comm.size();
+  int const me = plan.comm.rank();
+
+  // Build per-rank send and recv count arrays (size n, default 0).
+  std::vector<int> my_send_to(n, 0);
+  std::vector<int> my_recv_from(n, 0);
+  for (auto const &nc : plan.neighbors) {
+    my_send_to[nc.peer] = static_cast<int>(nc.send.size());
+    my_recv_from[nc.peer] = static_cast<int>(nc.recv.size());
+  }
+
+  // Collective all-to-all: element j of my_send_to goes to rank j.
+  // After this call peers_send_to_me[j] == rank j's send count toward me.
+  std::vector<int> peers_send_to_me(n, 0);
+  boost::mpi::all_to_all(plan.comm, my_send_to, peers_send_to_me);
+
+  // Check invariant: what I expect to receive from j == what j sends me.
+  for (int j = 0; j < n; ++j) {
+    if (j == me)
+      continue;
+    if (my_recv_from[j] != peers_send_to_me[j]) {
+      std::ostringstream oss;
+      oss << "symmetry mismatch with peer " << j << ": I expect to recv "
+          << my_recv_from[j] << " items but peer sends " << peers_send_to_me[j]
+          << " items";
+      violations.emplace_back(oss.str());
+    }
+  }
+
+  return violations;
+}
+
+bool report_violations(std::vector<std::string> const &violations,
+                       char const *context) {
+  for (auto const &v : violations) {
+    std::cerr << "halo plan validation (" << context << "): " << v << "\n";
+  }
+  return violations.empty();
+}
+
+} // namespace GhostComm
diff --git a/src/core/ghosts/HaloPlanValidator.hpp b/src/core/ghosts/HaloPlanValidator.hpp
new file mode 100644
index 0000000000..c8f1fbca01
--- /dev/null
+++ b/src/core/ghosts/HaloPlanValidator.hpp
@@ -0,0 +1,89 @@
+/*
+ * 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
+
+#include "cell_system/Cell.hpp"
+#include "ghosts/HaloPlan.hpp"
+
+#include <span>
+#include <string>
+#include <vector>
+
+namespace GhostComm {
+
+/**
+ * @brief Validate a HaloPlan for correctness.
+ *
+ * Checks:
+ * 1. Coverage: every ghost cell that a local cell references is filled.
+ *    Point-to-point ghosts appear as a recv/local.dst target exactly once; no
+ *    p2p target lies outside the ghost set; none appears twice. Ghosts filled
+ *    by the collective broadcast/reduce section (AtomDecomposition,
+ *    HybridDecomposition) are covered by that section instead of a recv/dst
+ *    target. A ghost that no local cell references (e.g. the halo-layer cells
+ *    on a single MPI rank, where the plan is intentionally empty) carries no
+ *    physics and is not required to be filled.
+ * 2. Neighborship-match: every local cell's ghost neighbor is covered, either
+ *    as a recv/dst target or by the collective section.
+ * 3. Peer-uniqueness: each peer value appears in at most one NeighborComm.
+ * 4. Shape: each NeighborComm has send.size() == recv.size().
+ * 5. Interior/boundary consistency: a cell marked interior has no ghost
+ *    neighbor (enforced by mark_boundary_cells()).
+ * 6. Overlap-safety invariant: an interior cell must not appear as a
+ *    NeighborComm send source or a LocalComm src.  This is the exact
+ *    precondition that the integrator step-2 / force-reduce overlap (Task
+ *    5.3) relies on: interior cells receive no reduce contributions and
+ *    therefore do not require the reduce to complete before their velocity
+ *    is updated.  LocalComm.dst and NeighborComm.recv are ghost cells by
+ *    construction (guaranteed by the dst-in-ghost-set coverage check #1),
+ *    so they are not checked here.
+ *
+ * @returns a human-readable violation string per problem; empty = valid.
+ */
+std::vector<std::string> validate_halo_plan(HaloPlan const &plan,
+                                            std::span<Cell *const> local_cells,
+                                            std::span<Cell *const> ghost_cells);
+
+/**
+ * @brief Cross-rank symmetry check for a HaloPlan.
+ *
+ * Uses a collective all-to-all exchange of per-rank send-counts so that every
+ * rank participates regardless of whether the neighbor sets are symmetric.
+ * This avoids the deadlock that a naive per-neighbor isend/irecv pattern
+ * would cause when a rank sends to a peer that has no matching recv posted.
+ *
+ * Invariant: for every peer P, my recv-count from P must equal P's send-count
+ * to me (i.e. my_recv_from[P] == peers_send_to_me[P]).
+ *
+ * @returns a human-readable violation string per mismatch; empty = valid.
+ */
+std::vector<std::string> validate_halo_plan_symmetry(HaloPlan const &plan);
+
+/**
+ * @brief Print violations to stderr and return whether the list was empty.
+ *
+ * Intended for use inside assert(): a bare
+ * assert(validate_halo_plan(...).empty()) aborts without showing WHICH
+ * invariant failed, which makes CI failures undebuggable.  Wrap the call:
+ * assert(GhostComm::report_violations(validate_halo_plan(...), "context")).
+ */
+bool report_violations(std::vector<std::string> const &violations,
+                       char const *context);
+
+} // namespace GhostComm
diff --git a/src/core/ghosts/mark_boundary_cells.hpp b/src/core/ghosts/mark_boundary_cells.hpp
new file mode 100644
index 0000000000..b32037c76b
--- /dev/null
+++ b/src/core/ghosts/mark_boundary_cells.hpp
@@ -0,0 +1,141 @@
+/*
+ * 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
+
+#include "cell_system/Cell.hpp"
+#include "ghosts/HaloPlan.hpp"
+
+#include <functional>
+#include <span>
+#include <unordered_map>
+#include <unordered_set>
+
+namespace GhostComm {
+
+/**
+ * @brief Classify each local cell as interior or boundary.
+ *
+ * A local cell is *boundary* for two independent reasons:
+ *
+ *   Source 1 — geometric wrap rules: the cell has a ghost neighbor (rule a)
+ *     or participates in a periodic wrap-around neighborship on a single MPI
+ *     rank (rule b, captured by the optional @p wrap_predicate).  These rules
+ *     serve future Euclidean-distance uses: wrap-neighborship must mark a cell
+ *     boundary even when it is not exported by the current plan.
+ *
+ *   Source 2 — plan membership: the cell appears as a send source in any
+ *     NeighborComm or as a LocalComm src.  Plans for geometries such as
+ *     Lees-Edwards fully-connected boundaries and ELC periodicity-change paths
+ *     can export cells that the geometric rules never see.  Those cells must
+ *     still be boundary so that the integrator step-2 / force-reduce overlap
+ *     does not update their velocities before remote force contributions
+ *     arrive.
+ *
+ * Both sources are complementary: keep both to satisfy present and future
+ * callers.  @ref mark_plan_cells_boundary applies source 2.
+ *
+ * All other local cells are *interior*.
+ *
+ * The function first resets every local cell to interior (idempotent on
+ * repeated calls / plan rebuild), then marks boundary cells.
+ *
+ * Call this right after `m_halo_plan = make_halo_plan()` in each
+ * decomposition, where `local_cells()` and `ghost_cells()` are already
+ * populated.
+ *
+ * @param local_cells  Local cell pointer span.
+ * @param ghost_cells  Ghost cell pointer span.
+ * @param wrap_predicate  Optional predicate `(Cell const *a, Cell const *b)
+ *   -> bool` that returns `true` when the neighbor relation a→b crosses a
+ *   periodic box boundary (both cells are local).  When omitted, no
+ *   additional wrapping boundary is detected.
+ */
+inline void mark_boundary_cells(
+    std::span<Cell *const> local_cells, std::span<Cell *const> ghost_cells,
+    std::function<bool(Cell const *, Cell const *)> wrap_predicate = nullptr) {
+  // Build a set of ghost ParticleList pointers for O(1) lookup.
+  std::unordered_set<ParticleList const *> ghost_set;
+  ghost_set.reserve(ghost_cells.size());
+  for (Cell *c : ghost_cells) {
+    ghost_set.insert(&c->particles());
+  }
+
+  // Reset all local cells to interior first (idempotent on rebuild).
+  for (Cell *c : local_cells) {
+    c->m_is_boundary = false;
+  }
+
+  // Mark a cell as boundary iff any of its neighbors is a ghost (rule a) or
+  // the wrap predicate fires for that neighbor pair (rule b).
+  for (Cell *c : local_cells) {
+    for (Cell *n : c->neighbors().all()) {
+      if (ghost_set.count(&n->particles()) ||
+          (wrap_predicate && wrap_predicate(c, n))) {
+        c->m_is_boundary = true;
+        break;
+      }
+    }
+  }
+}
+
+/**
+ * @brief Mark plan-exported local cells as boundary (source 2, see
+ * @ref mark_boundary_cells).
+ *
+ * Iterates the halo plan and marks every local cell that appears as a send
+ * source.  Specifically: every `NeighborComm::send[k].cell` and every
+ * `LocalComm::src`.  Recv/dst targets are ghost cells by construction and
+ * are not touched.
+ *
+ * Call this immediately after @ref mark_boundary_cells so that the combined
+ * invariant "interior ⇒ not exported by the plan" holds for any plan shape.
+ *
+ * @param plan        The halo plan produced by `make_halo_plan()`.
+ * @param local_cells Local cell pointer span (used to build the ParticleList →
+ *                    Cell reverse map in O(n)).
+ */
+inline void mark_plan_cells_boundary(HaloPlan const &plan,
+                                     std::span<Cell *const> local_cells) {
+  // Build a reverse map: ParticleList* → Cell* for O(1) lookup.
+  std::unordered_map<ParticleList const *, Cell *> pl_to_cell;
+  pl_to_cell.reserve(local_cells.size());
+  for (Cell *c : local_cells) {
+    pl_to_cell[&c->particles()] = c;
+  }
+
+  // Mark every NeighborComm send source boundary.
+  for (auto const &nc : plan.neighbors) {
+    for (auto const &sr : nc.send) {
+      auto it = pl_to_cell.find(sr.cell);
+      if (it != pl_to_cell.end()) {
+        it->second->m_is_boundary = true;
+      }
+    }
+  }
+
+  // Mark every LocalComm src boundary.
+  for (auto const &lc : plan.local) {
+    auto it = pl_to_cell.find(lc.src);
+    if (it != pl_to_cell.end()) {
+      it->second->m_is_boundary = true;
+    }
+  }
+}
+
+} // namespace GhostComm
diff --git a/src/core/ghosts/particle_packing.cpp b/src/core/ghosts/particle_packing.cpp
new file mode 100644
index 0000000000..01194abc3e
--- /dev/null
+++ b/src/core/ghosts/particle_packing.cpp
@@ -0,0 +1,362 @@
+/*
+ * Copyright (C) 2010-2026 The ESPResSo project
+ * Copyright (C) 2002,2003,2004,2005,2006,2007,2008,2009,2010
+ *   Max-Planck-Institute for Polymer Research, Theory Group
+ *
+ * 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/>.
+ */
+/** \file
+ *  Reusable particle packing/unpacking for ghost communications.
+ */
+
+#include "particle_packing.hpp"
+
+#include "BoxGeometry.hpp"
+#include "Particle.hpp"
+#include "ghosts.hpp"
+
+#include <utils/serialization/memcpy_archive.hpp>
+
+#include <boost/archive/binary_iarchive.hpp>
+#include <boost/archive/binary_oarchive.hpp>
+#include <boost/iostreams/device/array.hpp>
+#include <boost/iostreams/device/back_inserter.hpp>
+#include <boost/iostreams/stream.hpp>
+#include <boost/serialization/vector.hpp>
+
+#include <algorithm>
+#include <cassert>
+#include <cstddef>
+#include <iterator>
+#include <limits>
+#include <numeric>
+#include <span>
+#include <vector>
+
+namespace GhostComm {
+
+/** @brief Pseudo-archive to calculate the size of the serialization buffer. */
+class SerializationSizeCalculator {
+  std::size_t m_size = 0;
+
+public:
+  auto size() const { return m_size; }
+
+  template <class T> auto &operator<<(T &) {
+    m_size += sizeof(T);
+    return *this;
+  }
+
+  template <class T> auto &operator&(T &t) { return *this << t; }
+};
+
+/** @brief Type of reduction to carry out during serialization. */
+enum class ReductionPolicy {
+  /** @brief Reduction for domain-to-domain particle communication. */
+  MOVE,
+  /** @brief Reduction for cell-to-cell particle update. */
+  UPDATE,
+};
+
+/** @brief Whether to save the state to or load the state from the archive. */
+enum class SerializationDirection { SAVE, LOAD };
+
+/**
+ * @brief Serialize particle data, possibly with reduction.
+ * The reduction can take place during the save stage, e.g. to apply
+ * a ghost shift to the particle position, or during the load stage,
+ * e.g. to transfer momentum between particles in two local cells.
+ */
+template <class Archive>
+static void
+serialize_and_reduce(Archive &ar, Particle &p, unsigned int data_parts,
+                     ReductionPolicy policy, SerializationDirection direction,
+                     BoxGeometry const &box_geo,
+                     Utils::Vector3d const *ghost_shift) {
+  if (data_parts & GHOSTTRANS_PROPRTS) {
+    ar & p.id() & p.mol_id() & p.type() & p.propagation();
+#ifdef ESPRESSO_ROTATION
+    ar & p.rotation();
+#ifdef ESPRESSO_ROTATIONAL_INERTIA
+    ar & p.rinertia();
+#endif
+#endif
+#ifdef ESPRESSO_MASS
+    ar & p.mass();
+#endif
+#ifdef ESPRESSO_ELECTROSTATICS
+    ar & p.q();
+#endif
+#ifdef ESPRESSO_DIPOLES
+    ar & p.dipm();
+#endif
+#ifdef ESPRESSO_LB_ELECTROHYDRODYNAMICS
+    ar & p.mu_E();
+#endif
+#ifdef ESPRESSO_VIRTUAL_SITES_RELATIVE
+    ar & p.vs_relative();
+#endif
+#ifdef ESPRESSO_THERMOSTAT_PER_PARTICLE
+    ar & p.gamma();
+#ifdef ESPRESSO_ROTATION
+    ar & p.gamma_rot();
+#endif
+#endif
+#ifdef ESPRESSO_EXTERNAL_FORCES
+    ar & p.fixed();
+    ar & p.ext_force();
+#ifdef ESPRESSO_ROTATION
+    ar & p.ext_torque();
+#endif
+#endif
+#ifdef ESPRESSO_ENGINE
+    ar & p.swimming();
+#endif
+  }
+  if (data_parts & GHOSTTRANS_POSITION) {
+    if (direction == SerializationDirection::SAVE and ghost_shift != nullptr) {
+      /* ok, this is not nice, but perhaps fast */
+      auto pos = p.pos() + *ghost_shift;
+      auto img = p.image_box();
+      box_geo.fold_position(pos, img);
+      ar & pos;
+      ar & img;
+    } else {
+      ar & p.pos();
+      ar & p.image_box();
+    }
+#ifdef ESPRESSO_BOND_CONSTRAINT
+    ar & p.pos_last_time_step();
+#endif
+  }
+  // Wire-symmetry: pack and unpack use the same data_parts per exchange,
+  // so the layout is always symmetric; QUAT follows POSITION in the stream
+  // when both are set (position push), and TORQUE follows FORCE when both
+  // are set (force reduce).
+#ifdef ESPRESSO_ROTATION
+  if (data_parts & GHOSTTRANS_QUAT) {
+    ar & p.quat();
+  }
+#endif
+  if (data_parts & GHOSTTRANS_MOMENTUM) {
+    ar & p.v();
+#ifdef ESPRESSO_ROTATION
+    ar & p.omega();
+#endif
+  }
+  if (data_parts & GHOSTTRANS_FORCE) {
+    if (policy == ReductionPolicy::UPDATE and
+        direction == SerializationDirection::LOAD) {
+      Utils::Vector3d force;
+      ar & force;
+      p.force() += force;
+    } else {
+      ar & p.force();
+    }
+  }
+#ifdef ESPRESSO_ROTATION
+  if (data_parts & GHOSTTRANS_TORQUE) {
+    if (policy == ReductionPolicy::UPDATE and
+        direction == SerializationDirection::LOAD) {
+      Utils::Vector3d torque;
+      ar & torque;
+      p.torque() += torque;
+    } else {
+      ar & p.torque();
+    }
+  }
+#endif
+#ifdef ESPRESSO_BOND_CONSTRAINT
+  if (data_parts & GHOSTTRANS_RATTLE) {
+    if (policy == ReductionPolicy::UPDATE and
+        direction == SerializationDirection::LOAD) {
+      Utils::Vector3d correction;
+      ar & correction;
+      p.rattle_correction() += correction;
+    } else {
+      ar & p.rattle_correction();
+    }
+  }
+#endif
+}
+
+static void prepare_ghost_cell(ParticleList *cell, std::size_t size) {
+  /* Adapt size */
+  cell->resize(size);
+
+  /* Mark particles as ghosts */
+  for (auto &p : *cell) {
+    p.set_ghost(true);
+  }
+}
+
+std::size_t calc_transmit_size(BoxGeometry const &box_geo,
+                               unsigned data_parts) {
+  SerializationSizeCalculator sizeof_archive;
+  Particle p{};
+  serialize_and_reduce(sizeof_archive, p, data_parts, ReductionPolicy::MOVE,
+                       SerializationDirection::SAVE, box_geo, nullptr);
+  return sizeof_archive.size();
+}
+
+std::size_t calc_transmit_size(std::span<ParticleList *const> cells,
+                               BoxGeometry const &box_geo,
+                               unsigned data_parts) {
+  if (data_parts & GHOSTTRANS_PARTNUM)
+    return sizeof(unsigned int) * cells.size();
+
+  auto const n_part = std::accumulate(
+      cells.begin(), cells.end(), std::size_t{0},
+      [](std::size_t sum, auto part_list) { return sum + part_list->size(); });
+
+  return n_part * calc_transmit_size(box_geo, data_parts);
+}
+
+void pack_cells(CommBuf &buf, std::span<ParticleList *const> cells,
+                Utils::Vector3d const &shift, BoxGeometry const &box_geo,
+                unsigned data_parts) {
+  /* reallocate send buffer */
+  buf.resize(calc_transmit_size(cells, box_geo, data_parts));
+  buf.bonds().clear();
+
+  auto archiver = Utils::MemcpyOArchive{buf.make_span()};
+
+  /* Construct archive that pushes back to the bond buffer */
+  namespace io = boost::iostreams;
+  io::stream<io::back_insert_device<std::vector<char>>> os{
+      io::back_inserter(buf.bonds())};
+  boost::archive::binary_oarchive bond_archiver{os};
+
+  /* put in data */
+  for (auto part_list : cells) {
+    if (data_parts & GHOSTTRANS_PARTNUM) {
+      assert(part_list->size() <= std::numeric_limits<unsigned int>::max());
+      auto np = static_cast<unsigned int>(part_list->size());
+      archiver << np;
+    } else {
+      for (auto &p : *part_list) {
+        serialize_and_reduce(archiver, p, data_parts, ReductionPolicy::MOVE,
+                             SerializationDirection::SAVE, box_geo, &shift);
+        if (data_parts & GHOSTTRANS_BONDS) {
+          bond_archiver << p.bonds();
+        }
+      }
+    }
+  }
+
+  assert(archiver.bytes_written() == buf.size());
+}
+
+void unpack_cells(CommBuf &buf, std::span<ParticleList *const> cells,
+                  BoxGeometry const &box_geo, unsigned data_parts) {
+  /* put back data */
+  auto archiver = Utils::MemcpyIArchive{buf.make_span()};
+
+  if (data_parts & GHOSTTRANS_PARTNUM) {
+    for (auto part_list : cells) {
+      unsigned int np;
+      archiver >> np;
+      prepare_ghost_cell(part_list, np);
+    }
+  } else {
+    for (auto part_list : cells) {
+      for (auto &p : *part_list) {
+        serialize_and_reduce(archiver, p, data_parts, ReductionPolicy::MOVE,
+                             SerializationDirection::LOAD, box_geo, nullptr);
+      }
+    }
+    if (data_parts & GHOSTTRANS_BONDS) {
+      namespace io = boost::iostreams;
+      io::stream<io::array_source> bond_stream(
+          io::array_source{buf.bonds().data(), buf.bonds().size()});
+      boost::archive::binary_iarchive bond_archiver(bond_stream);
+
+      for (auto part_list : cells) {
+        for (auto &p : *part_list) {
+          bond_archiver >> p.bonds();
+        }
+      }
+    }
+  }
+
+  assert(archiver.bytes_read() == buf.size());
+
+  buf.bonds().clear();
+}
+
+void add_forces(CommBuf &buf, std::span<ParticleList *const> cells,
+                unsigned data_parts) {
+  /* put back data */
+  auto archiver = Utils::MemcpyIArchive{buf.make_span()};
+  for (auto &part_list : cells) {
+    for (Particle &part : *part_list) {
+      if (data_parts & GHOSTTRANS_FORCE) {
+        Utils::Vector3d force;
+        archiver >> force;
+        part.force() += force;
+      }
+#ifdef ESPRESSO_ROTATION
+      if (data_parts & GHOSTTRANS_TORQUE) {
+        Utils::Vector3d torque;
+        archiver >> torque;
+        part.torque() += torque;
+      }
+#endif
+    }
+  }
+}
+
+#ifdef ESPRESSO_BOND_CONSTRAINT
+void add_rattle(CommBuf &buf, std::span<ParticleList *const> cells) {
+  /* put back data */
+  auto archiver = Utils::MemcpyIArchive{buf.make_span()};
+  for (auto &part_list : cells) {
+    for (Particle &part : *part_list) {
+      ParticleRattle pr;
+      archiver >> pr;
+      part.rattle_params() += pr;
+    }
+  }
+}
+#endif
+
+void local_cell_copy(ParticleList &src, ParticleList &dst,
+                     Utils::Vector3d const &shift, BoxGeometry const &box_geo,
+                     unsigned data_parts) {
+  if (data_parts & GHOSTTRANS_PARTNUM) {
+    prepare_ghost_cell(&dst, src.size());
+  } else {
+    assert(src.size() == dst.size());
+    CommBuf buffer;
+    buffer.resize(calc_transmit_size(box_geo, data_parts));
+
+    for (std::size_t i = 0; i < src.size(); i++) {
+      auto ar_out = Utils::MemcpyOArchive{buffer.make_span()};
+      auto ar_in = Utils::MemcpyIArchive{buffer.make_span()};
+      auto &p1 = src.begin()[i];
+      auto &p2 = dst.begin()[i];
+      serialize_and_reduce(ar_out, p1, data_parts, ReductionPolicy::UPDATE,
+                           SerializationDirection::SAVE, box_geo, &shift);
+      serialize_and_reduce(ar_in, p2, data_parts, ReductionPolicy::UPDATE,
+                           SerializationDirection::LOAD, box_geo, nullptr);
+      if (data_parts & GHOSTTRANS_BONDS) {
+        p2.bonds() = p1.bonds();
+      }
+    }
+  }
+}
+
+} // namespace GhostComm
diff --git a/src/core/ghosts/particle_packing.hpp b/src/core/ghosts/particle_packing.hpp
new file mode 100644
index 0000000000..19c0fd5945
--- /dev/null
+++ b/src/core/ghosts/particle_packing.hpp
@@ -0,0 +1,206 @@
+/*
+ * Copyright (C) 2010-2026 The ESPResSo project
+ * Copyright (C) 2002,2003,2004,2005,2006,2007,2008,2009,2010
+ *   Max-Planck-Institute for Polymer Research, Theory Group
+ *
+ * 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
+ *  Reusable particle packing/unpacking for ghost communications.
+ *
+ *  Extracted from ghosts.cpp so that both the legacy ghost_communicator()
+ *  and future async engines can share the same byte-identical serialization.
+ */
+
+#include <config/config.hpp>
+
+#include "BoxGeometry.hpp"
+#include "ParticleList.hpp"
+
+#include <utils/Vector.hpp>
+
+#include <cstddef>
+#include <memory>
+#include <span>
+#include <vector>
+
+namespace GhostComm {
+
+namespace detail {
+/**
+ * @brief Allocator wrapper that skips value-initialization on grow.
+ *
+ * `std::vector<T, DefaultInitAllocator<T>>::resize(n)` extends storage
+ * without zero-filling when n > size(), leaving new bytes in an
+ * indeterminate state.  Use only when all grown bytes are overwritten
+ * before being read (e.g. by pack_cells / MPI irecv).
+ *
+ * All other allocator operations (allocate, deallocate, copy/move) are
+ * forwarded to `std::allocator<T>`.
+ */
+template <class T> struct DefaultInitAllocator : std::allocator<T> {
+  using value_type = T;
+  // Inherit all constructors from std::allocator<T>.
+  using std::allocator<T>::allocator;
+
+  template <class U> struct rebind {
+    using other = DefaultInitAllocator<U>;
+  };
+
+  /** Called by vector::resize for each newly appended element.
+   *  Plain `new (p) T` uses default-initialization: trivial types like
+   *  `char` are left uninitialized rather than zero-filled. */
+  void construct(T *p) noexcept { ::new (static_cast<void *>(p)) T; }
+
+  /** Forward non-default construction unchanged. */
+  template <class... Args> void construct(T *p, Args &&...args) {
+    ::new (static_cast<void *>(p)) T(std::forward<Args>(args)...);
+  }
+};
+} // namespace detail
+
+/**
+ * Class that stores marshalled data for ghost communications.
+ * To store and retrieve data, use the adapter classes below.
+ *
+ * Invariant: every byte in the non-bond buffer is written by the caller
+ * (pack_cells / MPI irecv) before it is read (unpack_cells / add_forces).
+ * The buffer therefore uses a default-init allocator so that resize()-on-grow
+ * does not zero-fill bytes that will be immediately overwritten, eliminating
+ * the `__memset_avx2` overhead visible in perf profiles at large particle
+ * counts.  The bond buffer is kept as a plain std::vector<char> because it is
+ * serialized via boost.mpi (which requires exact byte semantics) and is on the
+ * cold resort-only path.
+ */
+class CommBuf {
+  using ByteVec = std::vector<char, detail::DefaultInitAllocator<char>>;
+
+public:
+  /** Returns a pointer to the non-bond storage. */
+  char *data() { return buf.data(); }
+  const char *data() const { return buf.data(); }
+
+  /** Returns the number of elements in the non-bond storage. */
+  std::size_t size() const { return buf.size(); }
+
+  /** Resizes the underlying storage s.t. the object is capable of holding
+   *  @p new_size chars.  Bytes added by growth are left uninitialized —
+   *  the caller must write them (pack_cells / MPI irecv) before reading.
+   */
+  void resize(std::size_t new_size) { buf.resize(new_size); }
+
+  /** Returns a reference to the bond storage. */
+  auto &bonds() { return bondbuf; }
+  const auto &bonds() const { return bondbuf; }
+
+  auto make_span() { return std::span(buf.data(), buf.size()); }
+
+private:
+  ByteVec buf; ///< Buffer for everything but bonds (default-init grow)
+  std::vector<char>
+      bondbuf; ///< Buffer for bond lists (boost.mpi-serialized, cold path)
+};
+
+/**
+ * @brief Calculate the per-particle transmit size for the given data parts.
+ */
+std::size_t calc_transmit_size(BoxGeometry const &box_geo, unsigned data_parts);
+
+/**
+ * @brief Calculate the total transmit size for a list of cells and data parts.
+ *
+ * When GHOSTTRANS_PARTNUM is set, returns sizeof(unsigned int) per cell.
+ * Otherwise returns the total number of particles times the per-particle size.
+ */
+std::size_t calc_transmit_size(std::span<ParticleList *const> cells,
+                               BoxGeometry const &box_geo, unsigned data_parts);
+
+/**
+ * @brief Pack particle data from cells into a communication buffer.
+ *
+ * Equivalent to prepare_send_buffer() in the original ghosts.cpp.
+ * Handles GHOSTTRANS_PARTNUM (writes cell sizes) and GHOSTTRANS_BONDS
+ * (separate bond buffer) special cases.
+ *
+ * @param buf        Buffer to pack into (resized as needed).
+ * @param cells      Source particle lists.
+ * @param shift      Ghost shift applied to particle positions on save.
+ * @param box_geo    Box geometry for fold_position.
+ * @param data_parts Bitmask of GHOSTTRANS_* flags.
+ */
+void pack_cells(CommBuf &buf, std::span<ParticleList *const> cells,
+                Utils::Vector3d const &shift, BoxGeometry const &box_geo,
+                unsigned data_parts);
+
+/**
+ * @brief Unpack particle data from a communication buffer into cells.
+ *
+ * Equivalent to put_recv_buffer() in the original ghosts.cpp.
+ * Handles GHOSTTRANS_PARTNUM (resizes ghost cells) and GHOSTTRANS_BONDS
+ * (separate bond buffer) special cases.
+ *
+ * @param buf        Buffer to unpack from.
+ * @param cells      Destination particle lists.
+ * @param box_geo    Box geometry (unused here, kept for API symmetry).
+ * @param data_parts Bitmask of GHOSTTRANS_* flags.
+ */
+void unpack_cells(CommBuf &buf, std::span<ParticleList *const> cells,
+                  BoxGeometry const &box_geo, unsigned data_parts);
+
+/**
+ * @brief Add forces (and optionally torques) from a communication buffer
+ *        to particles in cells.
+ *
+ * @param buf        Buffer produced by a force-reduce exchange.
+ * @param cells      Destination (owned) particle lists.
+ * @param data_parts Bitmask of GHOSTTRANS_* flags; must include
+ *                   GHOSTTRANS_FORCE, may include GHOSTTRANS_TORQUE.
+ *
+ * Equivalent to add_forces_from_recv_buffer() in the original ghosts.cpp.
+ */
+void add_forces(CommBuf &buf, std::span<ParticleList *const> cells,
+                unsigned data_parts);
+
+#ifdef ESPRESSO_BOND_CONSTRAINT
+/**
+ * @brief Add rattle corrections from a communication buffer to particles.
+ *
+ * Equivalent to add_rattle_correction_from_recv_buffer() in ghosts.cpp.
+ */
+void add_rattle(CommBuf &buf, std::span<ParticleList *const> cells);
+#endif
+
+/**
+ * @brief Copy particle data from src to dst applying a ghost shift.
+ *
+ * Corresponds to one src/dst pair in cell_cell_transfer() from ghosts.cpp.
+ * For GHOSTTRANS_PARTNUM, resizes dst to match src.
+ * Otherwise serializes each particle from src and deserializes into dst,
+ * applying the shift and folding the position. Bond data is copied directly.
+ *
+ * @param src        Source particle list.
+ * @param dst        Destination particle list.
+ * @param shift      Ghost shift to apply.
+ * @param box_geo    Box geometry for fold_position.
+ * @param data_parts Bitmask of GHOSTTRANS_* flags.
+ */
+void local_cell_copy(ParticleList &src, ParticleList &dst,
+                     Utils::Vector3d const &shift, BoxGeometry const &box_geo,
+                     unsigned data_parts);
+
+} // namespace GhostComm
diff --git a/src/core/integrate.cpp b/src/core/integrate.cpp
index b6fd5609e6..067ad2ac8c 100644
--- a/src/core/integrate.cpp
+++ b/src/core/integrate.cpp
@@ -68,7 +68,7 @@
 #include <boost/mpi/collectives/all_reduce.hpp>
 
 #ifdef ESPRESSO_CALIPER
-#include <caliper/cali.h>
+#include "caliper_utils.hpp"
 #endif
 
 #ifdef ESPRESSO_VALGRIND
@@ -372,7 +372,7 @@ void walberla_agrid_sanity_checks(std::string const &method,
 
 static void resort_particles_if_needed(System::System &system) {
 #ifdef ESPRESSO_CALIPER
-  CALI_CXX_MARK_FUNCTION;
+  ESPRESSO_CALI_MARK_FUNCTION;
 #endif
   auto &cell_structure = *system.cell_structure;
   auto const offset = LeesEdwards::verlet_list_offset(
@@ -389,7 +389,7 @@ static bool integrator_step_1(CellStructure &cell_structure,
                               Propagation const &propagation,
                               System::System &system, double time_step) {
 #ifdef ESPRESSO_CALIPER
-  CALI_CXX_MARK_FUNCTION;
+  ESPRESSO_CALI_MARK_FUNCTION;
 #endif
   // steepest decent
   if (propagation.integ_switch == INTEG_METHOD_STEEPEST_DESCENT)
@@ -473,17 +473,25 @@ static bool integrator_step_1(CellStructure &cell_structure,
   return false;
 }
 
-static void integrator_step_2(CellStructure &cell_structure,
-                              Propagation const &propagation,
-                              [[maybe_unused]] System::System &system,
-                              double time_step) {
-#ifdef ESPRESSO_CALIPER
-  CALI_CXX_MARK_FUNCTION;
-#endif
-  if (propagation.integ_switch == INTEG_METHOD_STEEPEST_DESCENT)
-    return;
-
-  cell_structure.for_each_local_particle([&](Particle &p) {
+/**
+ * @brief Build the per-particle half-kick callable for step_2.
+ *
+ * Returns a lambda that captures @p propagation and @p time_step by reference
+ * and applies the velocity (and torque, if ROTATION is enabled) update for a
+ * single particle.  Virtual sites are skipped.
+ *
+ * Shared verbatim by @ref integrator_step_2 (full pass) and
+ * @ref integrator_step_2_filtered (interior / boundary passes): the lambda is
+ * constructed once, then handed to @c for_each_local_particle,
+ * @c for_each_interior_particle, or @c for_each_boundary_particle.
+ *
+ * NPT particles are intentionally absent: the NPT arm must run only on the
+ * ineligible (full-reduce-then-step_2) path and is handled separately inside
+ * @ref integrator_step_2.
+ */
+static auto make_step2_particle_kernel(Propagation const &propagation,
+                                       double time_step) {
+  return [&propagation, time_step](Particle &p) {
 #ifdef ESPRESSO_VIRTUAL_SITES
     // virtual sites are updated later in the integration loop
     if (p.is_virtual())
@@ -522,7 +530,21 @@ static void integrator_step_2(CellStructure &cell_structure,
         velocity_verlet_rotator_2(p, time_step);
 #endif
     }
-  });
+  };
+}
+
+static void integrator_step_2(CellStructure &cell_structure,
+                              Propagation const &propagation,
+                              [[maybe_unused]] System::System &system,
+                              double time_step) {
+#ifdef ESPRESSO_CALIPER
+  ESPRESSO_CALI_MARK_FUNCTION;
+#endif
+  if (propagation.integ_switch == INTEG_METHOD_STEEPEST_DESCENT)
+    return;
+
+  cell_structure.for_each_local_particle(
+      make_step2_particle_kernel(propagation, time_step));
 
 #ifdef ESPRESSO_NPT
   if ((propagation.used_propagations & PropagationMode::TRANS_LANGEVIN_NPT) and
@@ -539,9 +561,45 @@ static void integrator_step_2(CellStructure &cell_structure,
 #endif
 }
 
+/**
+ * @brief Restricted step_2 for the split-phase ghost-force-reduce overlap.
+ *
+ * Called twice per step: once with @p interior_pass = true (while the ghost
+ * reduce is in flight) and once with @p interior_pass = false (after
+ * @ref CellStructure::ghosts_reduce_forces_finish completes).  The per-particle
+ * kernel is identical to the one used by @ref integrator_step_2 — shared via
+ * @ref make_step2_particle_kernel.
+ *
+ * NPT is absent here: the eligibility check in calculate_forces() guarantees
+ * TRANS_LANGEVIN_NPT is never active when this path runs.
+ */
+static void integrator_step_2_filtered(CellStructure &cell_structure,
+                                       Propagation const &propagation,
+                                       double time_step, bool interior_pass) {
+#ifdef ESPRESSO_CALIPER
+  ESPRESSO_CALI_MARK_FUNCTION;
+#endif
+  // steepest_descent and NPT are excluded by the eligibility check; both
+  // asserts are belt-and-suspenders cross-checks.
+  assert(propagation.integ_switch != INTEG_METHOD_STEEPEST_DESCENT &&
+         "integrator_step_2_filtered: steepest-descent is ineligible");
+#ifdef ESPRESSO_NPT
+  assert((propagation.used_propagations &
+          PropagationMode::TRANS_LANGEVIN_NPT) == 0 &&
+         "integrator_step_2_filtered: NPT propagation is ineligible");
+#endif
+
+  auto const kernel = make_step2_particle_kernel(propagation, time_step);
+  if (interior_pass) {
+    cell_structure.for_each_interior_particle(kernel);
+  } else {
+    cell_structure.for_each_boundary_particle(kernel);
+  }
+}
+
 int System::System::integrate(int n_steps, int reuse_forces) {
 #ifdef ESPRESSO_CALIPER
-  CALI_CXX_MARK_FUNCTION;
+  ESPRESSO_CALI_MARK_FUNCTION;
 #endif
   auto &propagation = *this->propagation;
 #ifdef ESPRESSO_VIRTUAL_SITES_RELATIVE
@@ -576,7 +634,7 @@ int System::System::integrate(int n_steps, int reuse_forces) {
       ((reuse_forces != INTEG_REUSE_FORCES_ALWAYS) and
        propagation.recalc_forces)) {
 #ifdef ESPRESSO_CALIPER
-    CALI_MARK_BEGIN("Initial Force Calculation");
+    ESPRESSO_CALI_MARK_BEGIN("Initial Force Calculation");
 #endif
     thermostat->lb_coupling_deactivate();
 
@@ -596,6 +654,14 @@ int System::System::integrate(int n_steps, int reuse_forces) {
 
     calculate_forces();
 
+    // If calculate_forces started a split-phase ghost reduce, finish it now:
+    // the initial-force path has no step_2 to overlap with (n_steps may be 0,
+    // or the forces are for the previous step's record), so we just complete
+    // the reduce immediately.
+    if (cell_structure->has_pending_ghost_reduce()) {
+      cell_structure->ghosts_reduce_forces_finish();
+    }
+
     if (propagation.integ_switch != INTEG_METHOD_STEEPEST_DESCENT) {
 #ifdef ESPRESSO_ROTATION
       convert_initial_torques(cell_structure->local_particles());
@@ -603,7 +669,7 @@ int System::System::integrate(int n_steps, int reuse_forces) {
     }
 
 #ifdef ESPRESSO_CALIPER
-    CALI_MARK_END("Initial Force Calculation");
+    ESPRESSO_CALI_MARK_END("Initial Force Calculation");
 #endif
   }
 
@@ -636,12 +702,13 @@ int System::System::integrate(int n_steps, int reuse_forces) {
 #endif
   // Integration loop
 #ifdef ESPRESSO_CALIPER
-  CALI_CXX_MARK_LOOP_BEGIN(integration_loop, "Integration loop");
+  EspressoCaliLoop espresso_cali_integration_loop("Integration loop");
 #endif
   int integrated_steps = 0;
   for (int step = 0; step < n_steps; step++) {
 #ifdef ESPRESSO_CALIPER
-    CALI_CXX_MARK_LOOP_ITERATION(integration_loop, step);
+    auto espresso_cali_integration_iter =
+        espresso_cali_integration_loop.iteration(step);
 #endif
 
 #ifdef ESPRESSO_BOND_CONSTRAINT
@@ -717,11 +784,51 @@ int System::System::integrate(int n_steps, int reuse_forces) {
 #ifdef ESPRESSO_VIRTUAL_SITES_INERTIALESS_TRACERS
     if (thermostat->lb and
         (propagation.used_propagations & PropagationMode::TRANS_LB_TRACER)) {
+      // LB-tracer arm is ineligible for the split path; this block only runs
+      // on the blocking path where has_pending_ghost_reduce() is false.
+      assert(not cell_structure->has_pending_ghost_reduce() &&
+             "LB-tracer arm must be inactive on the split-phase path");
       lb_tracers_add_particle_force_to_fluid(*cell_structure, *box_geo,
                                              *local_geo, lb);
     }
 #endif
-    integrator_step_2(*cell_structure, propagation, *this, time_step);
+    if (cell_structure->has_pending_ghost_reduce()) {
+      // Split-phase path: interior half-kick runs while the ghost reduce is
+      // in flight; boundary half-kick runs after the reduce finishes.
+      // NPT, BD, steepest-descent, and LB-tracer are ineligible and never
+      // reach this branch (asserted inside integrator_step_2_filtered).
+      //
+      // RAII guard: if the interior pass throws (bad_alloc, Kokkos error, user
+      // callback), finish the pending reduce on unwind so MPI requests are not
+      // left dangling and the next start-assert does not fire.
+      // ESPResSo normally propagates errors via runtimeErrorMsg rather than
+      // exceptions, so this guard fires only in exceptional circumstances; the
+      // normal path calls finish() explicitly and the guard becomes a no-op.
+      struct ReduceGuard {
+        CellStructure *cs;
+        bool active;
+        ~ReduceGuard() {
+          if (active and cs->has_pending_ghost_reduce()) {
+            try {
+              cs->ghosts_reduce_forces_finish();
+            } catch (...) {
+              // The guard only runs during unwind from another exception;
+              // letting a second one escape this (implicitly noexcept)
+              // destructor would call std::terminate. Keep the original.
+            }
+          }
+        }
+      } guard{cell_structure.get(), true};
+
+      integrator_step_2_filtered(*cell_structure, propagation, time_step,
+                                 /*interior_pass=*/true);
+      cell_structure->ghosts_reduce_forces_finish();
+      guard.active = false; // normal path: finish already called above
+      integrator_step_2_filtered(*cell_structure, propagation, time_step,
+                                 /*interior_pass=*/false);
+    } else {
+      integrator_step_2(*cell_structure, propagation, *this, time_step);
+    }
     if (propagation.integ_switch == INTEG_METHOD_BD) {
       resort_particles_if_needed(*this);
     }
@@ -758,19 +865,19 @@ int System::System::integrate(int n_steps, int reuse_forces) {
           propagation.lb_skipped_md_steps = 0;
           propagation.ek_skipped_md_steps = 0;
 #ifdef ESPRESSO_CALIPER
-          CALI_MARK_BEGIN("lb_propagation");
+          ESPRESSO_CALI_MARK_BEGIN("lb_propagation");
 #endif
           lb.propagate();
           lb.ghost_communication_vel();
 #ifdef ESPRESSO_CALIPER
-          CALI_MARK_END("lb_propagation");
+          ESPRESSO_CALI_MARK_END("lb_propagation");
 #endif
 #ifdef ESPRESSO_CALIPER
-          CALI_MARK_BEGIN("ek_propagation");
+          ESPRESSO_CALI_MARK_BEGIN("ek_propagation");
 #endif
           ek.propagate();
 #ifdef ESPRESSO_CALIPER
-          CALI_MARK_END("ek_propagation");
+          ESPRESSO_CALI_MARK_END("ek_propagation");
 #endif
         }
       } else if (lb_active) {
@@ -779,11 +886,11 @@ int System::System::integrate(int n_steps, int reuse_forces) {
         if (propagation.lb_skipped_md_steps >= md_steps_per_lb_step) {
           propagation.lb_skipped_md_steps = 0;
 #ifdef ESPRESSO_CALIPER
-          CALI_MARK_BEGIN("lb_propagation");
+          ESPRESSO_CALI_MARK_BEGIN("lb_propagation");
 #endif
           lb.propagate();
 #ifdef ESPRESSO_CALIPER
-          CALI_MARK_END("lb_propagation");
+          ESPRESSO_CALI_MARK_END("lb_propagation");
 #endif
         }
       } else if (ek_active) {
@@ -792,11 +899,11 @@ int System::System::integrate(int n_steps, int reuse_forces) {
         if (propagation.ek_skipped_md_steps >= md_steps_per_ek_step) {
           propagation.ek_skipped_md_steps = 0;
 #ifdef ESPRESSO_CALIPER
-          CALI_MARK_BEGIN("ek_propagation");
+          ESPRESSO_CALI_MARK_BEGIN("ek_propagation");
 #endif
           ek.propagate();
 #ifdef ESPRESSO_CALIPER
-          CALI_MARK_END("ek_propagation");
+          ESPRESSO_CALI_MARK_END("ek_propagation");
 #endif
         }
       }
@@ -809,14 +916,14 @@ int System::System::integrate(int n_steps, int reuse_forces) {
       if (thermostat->lb and
           (propagation.used_propagations & PropagationMode::TRANS_LB_TRACER)) {
 #ifdef ESPRESSO_CALIPER
-        CALI_MARK_BEGIN("lb_tracers_propagation");
+        ESPRESSO_CALI_MARK_BEGIN("lb_tracers_propagation");
 #endif
         if (lb_active) {
           lb.ghost_communication_vel();
         }
         lb_tracers_propagate(*cell_structure, lb, time_step);
 #ifdef ESPRESSO_CALIPER
-        CALI_MARK_END("lb_tracers_propagation");
+        ESPRESSO_CALI_MARK_END("lb_tracers_propagation");
 #endif
       }
 #endif
@@ -847,9 +954,7 @@ int System::System::integrate(int n_steps, int reuse_forces) {
     lb.ghost_communication();
   }
   lees_edwards->update_box_params(*box_geo, sim_time);
-#ifdef ESPRESSO_CALIPER
-  CALI_CXX_MARK_LOOP_END(integration_loop);
-#endif
+  // espresso_cali_integration_loop destructor ends the Caliper loop region.
 
 #ifdef ESPRESSO_VALGRIND
   CALLGRIND_STOP_INSTRUMENTATION;
diff --git a/src/core/lb/particle_coupling.cpp b/src/core/lb/particle_coupling.cpp
index f47a840503..260fde92fa 100644
--- a/src/core/lb/particle_coupling.cpp
+++ b/src/core/lb/particle_coupling.cpp
@@ -36,7 +36,7 @@
 #include <boost/mpi.hpp>
 
 #ifdef ESPRESSO_CALIPER
-#include <caliper/cali.h>
+#include "caliper_utils.hpp"
 #endif
 
 #include <cassert>
@@ -347,7 +347,7 @@ static void lb_coupling_sanity_checks(Particle const &p) {
 
 void System::System::lb_couple_particles() {
 #ifdef ESPRESSO_CALIPER
-  CALI_CXX_MARK_FUNCTION;
+  ESPRESSO_CALI_MARK_FUNCTION;
 #endif
   assert(thermostat->lb != nullptr);
   if (thermostat->lb->couple_to_md) {
diff --git a/src/core/magnetostatics/dp3m_heffte.impl.hpp b/src/core/magnetostatics/dp3m_heffte.impl.hpp
index 9d9f6afecc..5a5f767b07 100644
--- a/src/core/magnetostatics/dp3m_heffte.impl.hpp
+++ b/src/core/magnetostatics/dp3m_heffte.impl.hpp
@@ -261,6 +261,7 @@ template <int cao> struct AssignTorques {
 
     auto const n_part = dp3m.inter_weights.size();
     auto const &unique_particles = cell_structure.get_unique_particles();
+    cell_structure.mark_torque_replicas_dirty();
     auto scatter_torque = cell_structure.get_scatter_torque();
     kokkos_parallel_range_for(
         "AssignTorques", std::size_t{0u}, n_part, [&](std::size_t p_index) {
diff --git a/src/core/nonbonded_interactions/nonbonded_interaction_data.hpp b/src/core/nonbonded_interactions/nonbonded_interaction_data.hpp
index 8c79583b98..3bd43c8e65 100644
--- a/src/core/nonbonded_interactions/nonbonded_interaction_data.hpp
+++ b/src/core/nonbonded_interactions/nonbonded_interaction_data.hpp
@@ -486,6 +486,11 @@ class InteractionsNonBonded : public System::Leaf<InteractionsNonBonded> {
   /** @brief Bitfield of currently active pair potentials. */
   auto combined_active_pair_mask() const { return m_combined_active_pair_mask; }
 
+  /** @brief Whether @p potential is active for any type pair. */
+  bool pair_potential_active(PairPotential potential) const {
+    return (m_combined_active_pair_mask & pair_potential_bit(potential)) != 0u;
+  }
+
 #ifdef ESPRESSO_THOLE
   /** @brief Whether any type pair has Thole damping configured. */
   auto any_thole_configured() const { return m_any_thole_configured; }
diff --git a/src/core/short_range_cabana.hpp b/src/core/short_range_cabana.hpp
index b2148a910a..6f26edee56 100644
--- a/src/core/short_range_cabana.hpp
+++ b/src/core/short_range_cabana.hpp
@@ -32,7 +32,7 @@
 #include <Cabana_NeighborList.hpp>
 
 #ifdef ESPRESSO_CALIPER
-#include <caliper/cali.h>
+#include "caliper_utils.hpp"
 #endif
 
 #include <functional>
@@ -193,7 +193,7 @@ update_cabana_state(CellStructure &cell_structure,
                     auto const &make_verlet_criterion, double const pair_cutoff,
                     auto const integ_switch) {
 #ifdef ESPRESSO_CALIPER
-  CALI_CXX_MARK_FUNCTION;
+  ESPRESSO_CALI_MARK_FUNCTION;
 #endif
   using policy_type = Kokkos::RangePolicy<execution_space>;
   auto const rebuild = cell_structure.prepare_verlet_list_cabana(pair_cutoff);
@@ -209,7 +209,7 @@ update_cabana_state(CellStructure &cell_structure,
     // Fill particle storage (full commit)
     // ===================================================
 #ifdef ESPRESSO_CALIPER
-    CALI_MARK_BEGIN("AoSoA commit full");
+    ESPRESSO_CALI_MARK_BEGIN("AoSoA commit full");
 #endif
     int pair_count = 0;
     int angle_count = 0;
@@ -265,7 +265,7 @@ update_cabana_state(CellStructure &cell_structure,
                          });
     Kokkos::fence();
 #ifdef ESPRESSO_CALIPER
-    CALI_MARK_END("AoSoA commit full");
+    ESPRESSO_CALI_MARK_END("AoSoA commit full");
 #endif
 
     // ===================================================
@@ -274,7 +274,7 @@ update_cabana_state(CellStructure &cell_structure,
     bool rebuild_vl = (integ_switch != INTEG_METHOD_STEEPEST_DESCENT and
                        cell_structure.use_verlet_list);
 #ifdef ESPRESSO_CALIPER
-    CALI_MARK_BEGIN("Verlet list creation");
+    ESPRESSO_CALI_MARK_BEGIN("Verlet list creation");
 #endif
     cell_structure.rebuild_verlet_list_cabana(
         [&](std::span<Cell *const> cells, BoxGeometry const &box,
@@ -303,14 +303,14 @@ update_cabana_state(CellStructure &cell_structure,
         },
         rebuild_vl);
 #ifdef ESPRESSO_CALIPER
-    CALI_MARK_END("Verlet list creation");
+    ESPRESSO_CALI_MARK_END("Verlet list creation");
 #endif
   } else {
     // ===================================================
     // Fill particle storage (partial update)
     // ===================================================
 #ifdef ESPRESSO_CALIPER
-    CALI_MARK_BEGIN("AoSoA commit partial");
+    ESPRESSO_CALI_MARK_BEGIN("AoSoA commit partial");
 #endif
 #ifdef ESPRESSO_EXCLUSIONS
     // commit_particle accumulates the any-exclusion aggregate below; clear it
@@ -325,7 +325,7 @@ update_cabana_state(CellStructure &cell_structure,
         });
     Kokkos::fence();
 #ifdef ESPRESSO_CALIPER
-    CALI_MARK_END("AoSoA commit partial");
+    ESPRESSO_CALI_MARK_END("AoSoA commit partial");
 #endif
   }
 }
@@ -372,7 +372,7 @@ void cabana_short_range(auto const &pair_bonds_kernel,
 
   if (bond_cutoff >= 0.) {
 #ifdef ESPRESSO_CALIPER
-    CALI_MARK_BEGIN("cabana_bond_loop");
+    ESPRESSO_CALI_MARK_BEGIN("cabana_bond_loop");
 #endif
     auto const n_pair_bonds = cell_structure.get_local_pair_bond_numbers();
     auto const n_angle_bonds = cell_structure.get_local_angle_bond_numbers();
@@ -403,14 +403,14 @@ void cabana_short_range(auto const &pair_bonds_kernel,
       Kokkos::fence();
     }
 #ifdef ESPRESSO_CALIPER
-    CALI_MARK_END("cabana_bond_loop");
+    ESPRESSO_CALI_MARK_END("cabana_bond_loop");
 #endif
   }
 
   // Cabana short range loop
   if (pair_cutoff > 0.) {
 #ifdef ESPRESSO_CALIPER
-    CALI_MARK_BEGIN("cabana_pair_loop");
+    ESPRESSO_CALI_MARK_BEGIN("cabana_pair_loop");
 #endif
     if (integ_switch != INTEG_METHOD_STEEPEST_DESCENT and
         cell_structure.use_verlet_list) {
@@ -445,7 +445,7 @@ void cabana_short_range(auto const &pair_bonds_kernel,
     }
     Kokkos::fence();
 #ifdef ESPRESSO_CALIPER
-    CALI_MARK_END("cabana_pair_loop");
+    ESPRESSO_CALI_MARK_END("cabana_pair_loop");
 #endif
   }
 }
diff --git a/src/core/short_range_loop.hpp b/src/core/short_range_loop.hpp
index 2c2e9734ce..57ddbdfa4d 100644
--- a/src/core/short_range_loop.hpp
+++ b/src/core/short_range_loop.hpp
@@ -24,7 +24,7 @@
 #include "cell_system/CellStructure.hpp"
 
 #ifdef ESPRESSO_CALIPER
-#include <caliper/cali.h>
+#include "caliper_utils.hpp"
 #endif
 
 #include <cassert>
@@ -46,7 +46,7 @@ void short_range_loop(BondKernel bond_kernel, PairKernel pair_kernel,
                       double bond_cutoff,
                       VerletCriterion const &verlet_criterion = {}) {
 #ifdef ESPRESSO_CALIPER
-  CALI_CXX_MARK_FUNCTION;
+  ESPRESSO_CALI_MARK_FUNCTION;
 #endif
 
   assert(cell_structure.get_resort_particles() == Cells::RESORT_NONE);
diff --git a/src/core/system/System.cpp b/src/core/system/System.cpp
index a22ecc3263..0a6175fe32 100644
--- a/src/core/system/System.cpp
+++ b/src/core/system/System.cpp
@@ -36,6 +36,7 @@
 #include "communication.hpp"
 #include "electrostatics/icc.hpp"
 #include "errorhandling.hpp"
+#include "ghosts.hpp"
 #include "integrators/steepest_descent.hpp"
 #include "nonbonded_interactions/VerletCriterion.hpp"
 #include "npt.hpp"
@@ -535,6 +536,95 @@ void System::on_integration_start() {
   on_observable_calc();
 }
 
+/**
+ * @brief Returns true when any active physics requires orientation of ghost
+ *        particles.  Used by both get_global_ghost_flags (QUAT push) and
+ *        get_force_reduce_ghost_flags (TORQUE reduce).
+ *
+ * Conservative: include a reader when in doubt.  Non-rotating plain LJ
+ * must yield false so that both bits stay OFF for that common case.
+ *
+ * Readers of ghost particle quat / director / calc_dip:
+ *  - short_range_cabana commit_particle: quat → director for Gay-Berne /
+ *    Dipoles pair kernels (ghost particles in AoSoA)
+ *  - vs_relative_update_particles: p_ref.quat() where p_ref may be a ghost
+ *  - LB particle coupling (swimmer): p.calc_director() on ghost particles
+ *    when ESPRESSO_ENGINE and LB are both active
+ *  - HomogeneousMagneticField constraint: p.calc_dip() (local particles only
+ *    — does NOT read ghosts; included conservatively via dipole check)
+ *  - dipolar solvers (dp3m, dds, dlc, scafacos): p.calc_dip() on
+ *    unique_particles which includes ghost particles that have no local
+ *    counterpart (see CellStructure::set_index_map, lines 310-322); covered by
+ *    the dipolar-solver-set condition
+ *  - ShapeBasedConstraint / calc_non_central_force: reads p.quat() on LOCAL
+ *    particles only (Constraints iterate local_particles) — safe without the
+ *    bit
+ *  - Stoner-Wohlfarth integrate_magnetodynamics: calls p_ref->calc_director()
+ *    where p_ref comes from get_reference_particle / get_local_particle, which
+ *    CAN return a ghost; covered by the explicit #ifdef below
+ *  - rotational integrators: operate on local particles only (OK)
+ *  - ICC blocking reduce: resets force_and_torque on local particles and may
+ *    carry a zero-valued TORQUE payload on the wire when orientation physics is
+ *    concurrently active — harmless (bytes only, value is zero)
+ */
+#ifdef ESPRESSO_ROTATION
+static bool orientation_ghosts_needed(System const &sys) {
+  // Rotational propagation active: ghost quat needed for vs_relative position
+  // update and for any kernel that uses the director of a ghost particle.
+  if (sys.propagation->used_propagations &
+      (PropagationMode::ROT_EULER | PropagationMode::ROT_LANGEVIN |
+       PropagationMode::ROT_BROWNIAN | PropagationMode::ROT_STOKESIAN |
+       PropagationMode::ROT_VS_RELATIVE |
+       PropagationMode::ROT_VS_INDEPENDENT)) {
+    return true;
+  }
+
+  // Virtual sites relative uses p_ref.quat() where p_ref may be a ghost.
+  if (sys.propagation->used_propagations &
+      (PropagationMode::TRANS_VS_RELATIVE | PropagationMode::ROT_VS_RELATIVE |
+       PropagationMode::ROT_VS_INDEPENDENT)) {
+    return true;
+  }
+
+#ifdef ESPRESSO_DIPOLES
+  // Dipolar solver set: particles have a dipole moment derived from quat;
+  // short_range_cabana reads quat → director for dipole pair kernels on ghosts.
+  if (sys.dipoles.impl && sys.dipoles.impl->solver) {
+    return true;
+  }
+#endif
+
+#ifdef ESPRESSO_GAY_BERNE
+  // Gay-Berne anisotropic nonbonded interaction configured: short_range_cabana
+  // commits ghost quat → director for the Cabana pair kernel.
+  if (sys.nonbonded_ias->pair_potential_active(PairPotential::GayBerne)) {
+    return true;
+  }
+#endif
+
+#ifdef ESPRESSO_ENGINE
+  // LB active with ENGINE compiled in: coupling loop includes ghost particles
+  // and calls p.calc_director() for swimmer_force_on_fluid particles.
+  if (sys.lb.is_solver_set()) {
+    return true;
+  }
+#endif
+
+#ifdef ESPRESSO_THERMAL_STONER_WOHLFARTH
+  // Stoner-Wohlfarth integrate_magnetodynamics calls
+  // get_reference_particle / get_local_particle which can return a ghost, then
+  // reads p_ref->calc_director() on it.  SW requires the Langevin thermostat
+  // (enforced in integrate.cpp integrator_sanity_checks); use that as the
+  // cheapest reliable system-level runtime signal.
+  if (sys.thermostat->thermo_switch & THERMO_LANGEVIN) {
+    return true;
+  }
+#endif
+
+  return false;
+}
+#endif // ESPRESSO_ROTATION
+
 /**
  * @brief Returns the ghost flags required for running pair
  *        kernels for the global state, e.g. the force calculation.
@@ -561,9 +651,25 @@ unsigned System::get_global_ghost_flags() const {
   }
 #endif
 
+#ifdef ESPRESSO_ROTATION
+  if (orientation_ghosts_needed(*this)) {
+    data_parts |= Cells::DATA_PART_QUAT;
+  }
+#endif
+
   return data_parts;
 }
 
+unsigned System::get_force_reduce_ghost_flags() const {
+  unsigned flags = GHOSTTRANS_FORCE;
+#ifdef ESPRESSO_ROTATION
+  if (orientation_ghosts_needed(*this)) {
+    flags |= GHOSTTRANS_TORQUE;
+  }
+#endif
+  return flags;
+}
+
 #ifdef ESPRESSO_NPT
 bool System::System::has_npt_enabled() const {
   return (propagation->integ_switch == INTEG_METHOD_NPT_ISO_AND) or
diff --git a/src/core/system/System.hpp b/src/core/system/System.hpp
index 1354795cb7..4e56d301ba 100644
--- a/src/core/system/System.hpp
+++ b/src/core/system/System.hpp
@@ -147,6 +147,14 @@ class System : public std::enable_shared_from_this<System> {
 
   unsigned get_global_ghost_flags() const;
 
+  /**
+   * @brief Returns the ghost flags for the force-reduce step.
+   *
+   * Always includes GHOSTTRANS_FORCE; adds GHOSTTRANS_TORQUE when orientation
+   * physics is active (same whitelist as get_global_ghost_flags for QUAT).
+   */
+  unsigned get_force_reduce_ghost_flags() const;
+
   /** Check electrostatic and magnetostatic methods are properly initialized.
    *  @return true if sanity checks failed.
    */
diff --git a/src/core/unit_tests/CMakeLists.txt b/src/core/unit_tests/CMakeLists.txt
index 1719a76045..98d9eea45a 100644
--- a/src/core/unit_tests/CMakeLists.txt
+++ b/src/core/unit_tests/CMakeLists.txt
@@ -90,6 +90,25 @@ espresso_unit_test(SRC bonded_interactions_map_test.cpp DEPENDS espresso::core
 espresso_unit_test(SRC bond_breakage_test.cpp DEPENDS espresso::core
                    Kokkos::kokkos)
 espresso_unit_test(SRC fene_test.cpp DEPENDS espresso::core Kokkos::kokkos)
+espresso_unit_test(SRC Cell_boundary_test.cpp DEPENDS espresso::core
+                   Kokkos::kokkos)
+espresso_unit_test(SRC HaloPlan_test.cpp DEPENDS espresso::core Kokkos::kokkos
+                   Boost::mpi MPI::MPI_CXX)
+espresso_unit_test(SRC RegularHaloPlan_test.cpp DEPENDS espresso::core
+                   Kokkos::kokkos Boost::mpi MPI::MPI_CXX NUM_PROC 4)
+espresso_unit_test(SRC HaloExchange_test.cpp DEPENDS espresso::core
+                   Kokkos::kokkos Boost::mpi MPI::MPI_CXX NUM_PROC 2)
+# Register the validator test on both 4 ranks (exercises face/edge/corner peers
+# on a non-trivial node grid) and 1 rank (where make_halo_plan is intentionally
+# empty and the halo-layer cells are unreferenced ghosts).
+espresso_unit_test_executable(
+  NAME HaloPlanValidator_test SRC HaloPlanValidator_test.cpp DEPENDS
+  espresso::core Kokkos::kokkos Boost::mpi MPI::MPI_CXX)
+foreach(TEST_NUM_PROC 1 4)
+  espresso_unit_test_register(
+    NAME HaloPlanValidator_test_${TEST_NUM_PROC}_mpi_ranks TARGET
+    HaloPlanValidator_test NUM_PROC ${TEST_NUM_PROC})
+endforeach()
 espresso_unit_test(SRC TabulatedPotential_test.cpp DEPENDS espresso::core
                    Kokkos::kokkos)
 if(ESPRESSO_BUILD_WITH_WALBERLA)
diff --git a/src/core/unit_tests/Cell_boundary_test.cpp b/src/core/unit_tests/Cell_boundary_test.cpp
new file mode 100644
index 0000000000..ccc28eb04c
--- /dev/null
+++ b/src/core/unit_tests/Cell_boundary_test.cpp
@@ -0,0 +1,325 @@
+/*
+ * 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/>.
+ */
+
+/**
+ * @file
+ * TDD for Task 3.1: interior/boundary cell classification.
+ *
+ * Tests:
+ *  1. mark_boundary_cells: a local cell with only local neighbors → interior.
+ *  2. mark_boundary_cells: a local cell with a ghost neighbor → boundary.
+ *  3. mark_boundary_cells is idempotent (second call gives same result).
+ *  4. validate_halo_plan consistency check: fires when a cell is mis-marked
+ *     interior but has a ghost neighbor.
+ *  5. validate_halo_plan consistency check: silent for a correctly-marked
+ *     boundary cell.
+ */
+
+#define BOOST_TEST_MODULE Cell_boundary test
+#define BOOST_TEST_DYN_LINK
+#include <boost/test/unit_test.hpp>
+
+#include "cell_system/Cell.hpp"
+#include "ghosts/HaloPlan.hpp"
+#include "ghosts/HaloPlanValidator.hpp"
+#include "ghosts/mark_boundary_cells.hpp"
+
+#include <span>
+#include <vector>
+
+// ── mark_boundary_cells ──────────────────────────────────────────────────────
+
+// A local cell whose only neighbors are other local cells must be interior.
+BOOST_AUTO_TEST_CASE(local_only_neighbors_is_interior) {
+  Cell local_a, local_b, local_c;
+  // local_a has two neighbors, both local (not ghost)
+  local_a.m_neighbors =
+      Neighbors<Cell *>(std::vector<Cell *>{&local_b, &local_c}, {});
+
+  std::vector<Cell *> locals{&local_a, &local_b, &local_c};
+  std::vector<Cell *> ghosts{};
+
+  GhostComm::mark_boundary_cells(locals, ghosts);
+
+  BOOST_CHECK_MESSAGE(!local_a.is_boundary(),
+                      "cell with only local neighbors must be interior");
+  BOOST_CHECK_MESSAGE(!local_b.is_boundary(),
+                      "local_b (no neighbors listed) must be interior");
+  BOOST_CHECK_MESSAGE(!local_c.is_boundary(),
+                      "local_c (no neighbors listed) must be interior");
+}
+
+// A local cell that has at least one ghost neighbor must be boundary.
+BOOST_AUTO_TEST_CASE(ghost_neighbor_makes_boundary) {
+  Cell local_a, local_b, ghost;
+  // local_a has one local neighbor and one ghost neighbor
+  local_a.m_neighbors =
+      Neighbors<Cell *>(std::vector<Cell *>{&local_b, &ghost}, {});
+  // local_b has only a local neighbor
+  local_b.m_neighbors = Neighbors<Cell *>(std::vector<Cell *>{&local_a}, {});
+
+  std::vector<Cell *> locals{&local_a, &local_b};
+  std::vector<Cell *> ghosts{&ghost};
+
+  GhostComm::mark_boundary_cells(locals, ghosts);
+
+  BOOST_CHECK_MESSAGE(local_a.is_boundary(),
+                      "cell with a ghost neighbor must be boundary");
+  BOOST_CHECK_MESSAGE(!local_b.is_boundary(),
+                      "cell without ghost neighbors must be interior");
+}
+
+// Calling mark_boundary_cells a second time (rebuild) must yield the same
+// result (reset + re-mark = idempotent).
+BOOST_AUTO_TEST_CASE(mark_boundary_cells_is_idempotent) {
+  Cell local, ghost;
+  local.m_neighbors = Neighbors<Cell *>(std::vector<Cell *>{&ghost}, {});
+  std::vector<Cell *> locals{&local};
+  std::vector<Cell *> ghosts{&ghost};
+
+  GhostComm::mark_boundary_cells(locals, ghosts);
+  bool const first_result = local.is_boundary();
+
+  GhostComm::mark_boundary_cells(locals, ghosts);
+  bool const second_result = local.is_boundary();
+
+  BOOST_CHECK_EQUAL(first_result, second_result);
+  BOOST_CHECK(first_result); // should be boundary
+}
+
+// ── wrap-predicate (Task 5.2)
+// ─────────────────────────────────────────────────
+
+// A pair of local cells connected by a wrap-predicate that always fires must
+// both be marked boundary, even without any ghost neighbors.
+BOOST_AUTO_TEST_CASE(wrap_predicate_makes_boundary) {
+  Cell local_a, local_b;
+  // local_a and local_b are mutual local neighbors (no ghosts).
+  local_a.m_neighbors = Neighbors<Cell *>(std::vector<Cell *>{&local_b}, {});
+  local_b.m_neighbors = Neighbors<Cell *>(std::vector<Cell *>{&local_a}, {});
+
+  std::vector<Cell *> locals{&local_a, &local_b};
+  std::vector<Cell *> ghosts{};
+
+  // Predicate that fires for every pair: simulates the wrap axis case.
+  auto always_wrap = [](Cell const *, Cell const *) { return true; };
+
+  GhostComm::mark_boundary_cells(locals, ghosts, always_wrap);
+
+  BOOST_CHECK_MESSAGE(
+      local_a.is_boundary(),
+      "cell with a wrap-predicate neighbor must be boundary (no ghost needed)");
+  BOOST_CHECK_MESSAGE(
+      local_b.is_boundary(),
+      "cell with a wrap-predicate neighbor must be boundary (no ghost needed)");
+}
+
+// A pair of fully-local cells with NO ghost neighbors and a wrap predicate
+// that never fires must remain interior.
+BOOST_AUTO_TEST_CASE(no_wrap_no_ghost_stays_interior) {
+  Cell local_a, local_b;
+  local_a.m_neighbors = Neighbors<Cell *>(std::vector<Cell *>{&local_b}, {});
+  local_b.m_neighbors = Neighbors<Cell *>(std::vector<Cell *>{&local_a}, {});
+
+  std::vector<Cell *> locals{&local_a, &local_b};
+  std::vector<Cell *> ghosts{};
+
+  // Predicate that never fires: no wrap axis.
+  auto never_wrap = [](Cell const *, Cell const *) { return false; };
+
+  GhostComm::mark_boundary_cells(locals, ghosts, never_wrap);
+
+  BOOST_CHECK_MESSAGE(
+      !local_a.is_boundary(),
+      "cell with only local non-wrap neighbors must be interior");
+  BOOST_CHECK_MESSAGE(
+      !local_b.is_boundary(),
+      "cell with only local non-wrap neighbors must be interior");
+}
+
+// ── validator overlap-safety invariant (Task 5.2) ────────────────────────────
+
+// An interior cell that appears as a NeighborComm send source must trigger
+// the overlap-safety invariant violation.
+BOOST_AUTO_TEST_CASE(validator_fires_for_interior_cell_in_send_region) {
+  using namespace GhostComm;
+  // Two local cells: local_a is interior (no ghost neighbors),
+  // local_b is boundary (ghost neighbor).
+  Cell local_a, local_b, ghost;
+  local_a.m_neighbors = Neighbors<Cell *>(std::vector<Cell *>{&local_b}, {});
+  local_b.m_neighbors =
+      Neighbors<Cell *>(std::vector<Cell *>{&local_a, &ghost}, {});
+
+  std::vector<Cell *> locals{&local_a, &local_b};
+  std::vector<Cell *> ghosts{&ghost};
+
+  GhostComm::mark_boundary_cells(locals, ghosts);
+  BOOST_REQUIRE(!local_a.is_boundary()); // pre-condition: local_a is interior
+  BOOST_REQUIRE(local_b.is_boundary());  // pre-condition: local_b is boundary
+
+  // Build a plan that covers the ghost (so coverage checks pass), but
+  // incorrectly lists local_a (interior) as a send source.
+  HaloPlan plan;
+  plan.neighbors.push_back(
+      NeighborComm{1,
+                   {SendRegion{&local_a.particles(), {}}, // ← interior cell!
+                    SendRegion{&local_b.particles(), {}}},
+                   {&ghost.particles(), &ghost.particles()}});
+  // Fix double-fill by using a separate ghost for the second recv slot.
+  Cell ghost2;
+  ghosts.push_back(&ghost2);
+  plan.neighbors[0].recv = {&ghost.particles(), &ghost2.particles()};
+
+  auto violations = validate_halo_plan(plan, locals, ghosts);
+  bool found_overlap =
+      std::ranges::any_of(violations, [](std::string const &v) {
+        return v.find("overlap-safety invariant") != std::string::npos;
+      });
+  BOOST_CHECK_MESSAGE(found_overlap,
+                      "expected overlap-safety invariant violation; got: " +
+                          (violations.empty() ? "(none)" : violations.front()));
+}
+
+// An interior cell that appears as a LocalComm src must trigger the
+// overlap-safety invariant violation.
+BOOST_AUTO_TEST_CASE(validator_fires_for_interior_cell_in_local_comm_src) {
+  using namespace GhostComm;
+  Cell local_interior, local_boundary, ghost;
+  local_interior.m_neighbors =
+      Neighbors<Cell *>(std::vector<Cell *>{&local_boundary}, {});
+  local_boundary.m_neighbors =
+      Neighbors<Cell *>(std::vector<Cell *>{&local_interior, &ghost}, {});
+
+  std::vector<Cell *> locals{&local_interior, &local_boundary};
+  std::vector<Cell *> ghosts{&ghost};
+
+  GhostComm::mark_boundary_cells(locals, ghosts);
+  BOOST_REQUIRE(!local_interior.is_boundary());
+  BOOST_REQUIRE(local_boundary.is_boundary());
+
+  // Plan with a LocalComm that uses local_interior as src → violation.
+  HaloPlan plan;
+  plan.local.push_back(
+      LocalComm{&local_interior.particles(), &ghost.particles(), {}});
+
+  auto violations = validate_halo_plan(plan, locals, ghosts);
+  bool found_overlap =
+      std::ranges::any_of(violations, [](std::string const &v) {
+        return v.find("overlap-safety invariant") != std::string::npos;
+      });
+  BOOST_CHECK_MESSAGE(found_overlap,
+                      "expected overlap-safety invariant violation; got: " +
+                          (violations.empty() ? "(none)" : violations.front()));
+}
+
+// A correctly-marked boundary cell as send source must NOT trigger the
+// overlap-safety invariant (sanity check that we don't over-fire).
+BOOST_AUTO_TEST_CASE(validator_silent_for_boundary_cell_in_send_region) {
+  using namespace GhostComm;
+  Cell local_boundary, ghost;
+  local_boundary.m_neighbors =
+      Neighbors<Cell *>(std::vector<Cell *>{&ghost}, {});
+
+  std::vector<Cell *> locals{&local_boundary};
+  std::vector<Cell *> ghosts{&ghost};
+
+  GhostComm::mark_boundary_cells(locals, ghosts);
+  BOOST_REQUIRE(local_boundary.is_boundary());
+
+  HaloPlan plan;
+  plan.neighbors.push_back(NeighborComm{
+      1, {SendRegion{&local_boundary.particles(), {}}}, {&ghost.particles()}});
+
+  auto violations = validate_halo_plan(plan, locals, ghosts);
+  bool found_overlap =
+      std::ranges::any_of(violations, [](std::string const &v) {
+        return v.find("overlap-safety invariant") != std::string::npos;
+      });
+  BOOST_CHECK_MESSAGE(!found_overlap,
+                      "boundary cell as send source must not trigger overlap "
+                      "invariant; violations: " +
+                          (violations.empty() ? "(none)" : violations.front()));
+  BOOST_CHECK_MESSAGE(violations.empty(),
+                      "expected no violations for a correct plan");
+}
+
+// ── validator consistency check ──────────────────────────────────────────────
+
+// A local cell whose is_boundary()==false but that has a ghost neighbor must
+// trigger the "interior cell has a ghost neighbor" consistency violation.
+BOOST_AUTO_TEST_CASE(validator_fires_for_mismarked_interior_cell) {
+  using namespace GhostComm;
+  Cell local, ghost;
+  // Wire ghost as a neighbor of local, but do NOT call mark_boundary_cells,
+  // so local.is_boundary() stays false (the default).
+  local.m_neighbors = Neighbors<Cell *>(std::vector<Cell *>{&ghost}, {});
+
+  std::vector<Cell *> locals{&local};
+  std::vector<Cell *> ghosts{&ghost};
+
+  // Build a minimal valid plan that covers the ghost so that the coverage
+  // checks pass — we only want to test the consistency check in isolation.
+  HaloPlan plan;
+  plan.neighbors.push_back(NeighborComm{
+      1, {SendRegion{&local.particles(), {}}}, {&ghost.particles()}});
+
+  auto violations = validate_halo_plan(plan, locals, ghosts);
+  bool found_consistency =
+      std::ranges::any_of(violations, [](std::string const &v) {
+        return v.find("interior cell has a ghost neighbor") !=
+               std::string::npos;
+      });
+  BOOST_CHECK_MESSAGE(
+      found_consistency,
+      "expected 'interior cell has a ghost neighbor' violation; "
+      "got: " +
+          (violations.empty() ? "(none)" : violations.front()));
+}
+
+// A correctly-marked boundary cell must NOT trigger the consistency violation.
+BOOST_AUTO_TEST_CASE(validator_silent_for_correctly_marked_boundary_cell) {
+  using namespace GhostComm;
+  Cell local, ghost;
+  local.m_neighbors = Neighbors<Cell *>(std::vector<Cell *>{&ghost}, {});
+
+  std::vector<Cell *> locals{&local};
+  std::vector<Cell *> ghosts{&ghost};
+
+  // Properly mark the cell as boundary first.
+  GhostComm::mark_boundary_cells(locals, ghosts);
+  BOOST_REQUIRE(local.is_boundary()); // pre-condition for this test
+
+  HaloPlan plan;
+  plan.neighbors.push_back(NeighborComm{
+      1, {SendRegion{&local.particles(), {}}}, {&ghost.particles()}});
+
+  auto violations = validate_halo_plan(plan, locals, ghosts);
+  bool found_consistency =
+      std::ranges::any_of(violations, [](std::string const &v) {
+        return v.find("interior cell has a ghost neighbor") !=
+               std::string::npos;
+      });
+  BOOST_CHECK_MESSAGE(
+      !found_consistency,
+      "expected no consistency violation for a correctly-marked "
+      "boundary cell; got: " +
+          (violations.empty() ? "(none)" : violations.front()));
+  BOOST_CHECK_MESSAGE(violations.empty(),
+                      "expected no violations at all for a correct plan");
+}
diff --git a/src/core/unit_tests/HaloExchange_test.cpp b/src/core/unit_tests/HaloExchange_test.cpp
new file mode 100644
index 0000000000..a4f352e4df
--- /dev/null
+++ b/src/core/unit_tests/HaloExchange_test.cpp
@@ -0,0 +1,351 @@
+/*
+ * Copyright (C) 2010-2026 The ESPResSo project
+ * Copyright (C) 2002,2003,2004,2005,2006,2007,2008,2009,2010
+ *   Max-Planck-Institute for Polymer Research, Theory Group
+ *
+ * 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/>.
+ */
+#define BOOST_TEST_NO_MAIN
+#define BOOST_TEST_MODULE HaloExchange test
+#define BOOST_TEST_ALTERNATIVE_INIT_API
+#define BOOST_TEST_DYN_LINK
+#include <boost/test/unit_test.hpp>
+
+#include <boost/mpi.hpp>
+
+#include "BondList.hpp"
+#include "BoxGeometry.hpp"
+#include "Particle.hpp"
+#include "ParticleList.hpp"
+#include "ghosts.hpp"
+#include "ghosts/HaloExchange.hpp"
+#include "ghosts/HaloPlan.hpp"
+
+#include <array>
+
+namespace utf = boost::unit_test;
+
+/*
+ * On two ranks, each rank owns a single particle tagged by its rank. A Push
+ * of GHOSTTRANS_POSITION must overwrite the local ghost cell with the peer's
+ * position. This exercises the deadlock-free irecv/isend protocol and the
+ * a-priori recv sizing.
+ */
+BOOST_AUTO_TEST_CASE(push_positions_between_two_ranks,
+                     *utf::precondition([](utf::test_unit_id) {
+                       return boost::mpi::communicator{}.size() == 2;
+                     })) {
+  using namespace GhostComm;
+  boost::mpi::communicator world;
+  int const me = world.rank();
+  int const other = 1 - me;
+
+  BoxGeometry box;
+  box.set_length({10., 10., 10.});
+
+  ParticleList local, ghost;
+  local.resize(1);
+  local.begin()->id() = me; // tag by owner rank
+  local.begin()->pos() = {double(me), 2. * double(me), 3. * double(me)};
+  ghost.resize(1); // sized as if ghosts_count() had already run
+
+  HaloPlan plan;
+  plan.comm = world;
+  plan.neighbors.push_back(
+      NeighborComm{other, {SendRegion{&local, {}}}, {&ghost}});
+
+  halo_exchange(plan, box, GHOSTTRANS_POSITION,
+                {Direction::Push, Combine::Overwrite});
+
+  // ghost cell now holds the peer's (folded) position.
+  BOOST_CHECK_CLOSE(ghost.begin()->pos()[0], double(other), 1e-12);
+  BOOST_CHECK_CLOSE(ghost.begin()->pos()[1], 2. * double(other), 1e-12);
+  BOOST_CHECK_CLOSE(ghost.begin()->pos()[2], 3. * double(other), 1e-12);
+  // the local (owned) particle must be untouched by a Push.
+  BOOST_CHECK_CLOSE(local.begin()->pos()[0], double(me), 1e-12);
+}
+
+/*
+ * A Reduce with Combine::Add of GHOSTTRANS_FORCE folds ghost forces back onto
+ * the owner. The send/recv roles swap: we send from the *ghost* cell and add
+ * into the *owned* cell. Each rank seeds a distinct force on its ghost, so the
+ * owner's force must end up as (its own pre-existing force + the peer's ghost
+ * force).
+ */
+BOOST_AUTO_TEST_CASE(reduce_forces_between_two_ranks,
+                     *utf::precondition([](utf::test_unit_id) {
+                       return boost::mpi::communicator{}.size() == 2;
+                     })) {
+  using namespace GhostComm;
+  boost::mpi::communicator world;
+  int const me = world.rank();
+  int const other = 1 - me;
+
+  BoxGeometry box;
+  box.set_length({10., 10., 10.});
+
+  ParticleList local, ghost;
+  local.resize(1);
+  ghost.resize(1);
+  // Pre-existing force on the owned particle; Reduce must add to it.
+  local.begin()->force() = {1., 0., 0.};
+  // Force accumulated on this rank's ghost (from short-range interactions).
+  // Encode the ghost's owner (== other) so the sum is checkable.
+  ghost.begin()->force() = {0., 10. * double(me + 1), 100. * double(me + 1)};
+
+  HaloPlan plan;
+  plan.comm = world;
+  plan.neighbors.push_back(
+      NeighborComm{other, {SendRegion{&local, {}}}, {&ghost}});
+
+  halo_exchange(plan, box, GHOSTTRANS_FORCE, {Direction::Reduce, Combine::Add});
+
+  // Owner force == own pre-existing force + peer's ghost force.
+  BOOST_CHECK_CLOSE(local.begin()->force()[0], 1., 1e-12);
+  BOOST_CHECK_CLOSE(local.begin()->force()[1], 10. * double(other + 1), 1e-12);
+  BOOST_CHECK_CLOSE(local.begin()->force()[2], 100. * double(other + 1), 1e-12);
+}
+
+/*
+ * The cold, resort-only BONDS path is exchanged as a *second* per-neighbor
+ * message (the CommBuf bond buffer) alongside the flat particle buffer. This
+ * checks that a bond placed on the owner round-trips byte-correct into the
+ * peer's ghost cell over the async engine.
+ */
+BOOST_AUTO_TEST_CASE(push_bonds_between_two_ranks,
+                     *utf::precondition([](utf::test_unit_id) {
+                       return boost::mpi::communicator{}.size() == 2;
+                     })) {
+  using namespace GhostComm;
+  boost::mpi::communicator world;
+  int const me = world.rank();
+  int const other = 1 - me;
+
+  BoxGeometry box;
+  box.set_length({10., 10., 10.});
+
+  ParticleList local, ghost;
+  local.resize(1);
+  local.begin()->id() = me;
+  std::array<int, 2> partners{me, me + 5};
+  local.begin()->bonds().insert(BondView{me + 1, partners});
+  ghost.resize(1);
+
+  HaloPlan plan;
+  plan.comm = world;
+  plan.neighbors.push_back(
+      NeighborComm{other, {SendRegion{&local, {}}}, {&ghost}});
+
+  halo_exchange(plan, box, GHOSTTRANS_PROPRTS | GHOSTTRANS_BONDS,
+                {Direction::Push, Combine::Overwrite});
+
+  BOOST_REQUIRE_EQUAL(ghost.begin()->bonds().size(), 1);
+  auto it = ghost.begin()->bonds().begin();
+  BOOST_CHECK_EQUAL((*it).bond_id(), other + 1);
+  BOOST_REQUIRE_EQUAL((*it).partner_ids().size(), 2u);
+  BOOST_CHECK_EQUAL((*it).partner_ids()[0], other);
+  BOOST_CHECK_EQUAL((*it).partner_ids()[1], other + 5);
+}
+
+/*
+ * Collective section: on two ranks each rank owns one cell.  A Broadcast (Push)
+ * must copy each rank's owned cell data into the other rank's ghost cell (the
+ * entry for that rank in the collective cells vector).  The Broadcast loop runs
+ * root=0 then root=1; after root=0's broadcast, rank 1's cells[0] holds rank
+ * 0's particle; after root=1's broadcast, rank 0's cells[1] holds rank 1's
+ * particle.
+ *
+ * A subsequent ReduceSum (Reduce) must sum ghost-force contributions back to
+ * each owner: rank 0's cells[0] and rank 1's cells[1] accumulate the sums.
+ */
+BOOST_AUTO_TEST_CASE(collective_broadcast_and_reduce,
+                     *utf::precondition([](utf::test_unit_id) {
+                       return boost::mpi::communicator{}.size() == 2;
+                     })) {
+  using namespace GhostComm;
+  boost::mpi::communicator world;
+  int const me = world.rank();
+
+  BoxGeometry box;
+  box.set_length({10., 10., 10.});
+
+  // cells[0] = local cell for rank 0, cells[1] = local cell for rank 1.
+  // Each rank owns cells[me] and uses cells[1 - me] as ghost storage.
+  ParticleList cell0, cell1;
+  cell0.resize(1);
+  cell1.resize(1);
+
+  // Seed the owned cell with a distinguishable position.
+  if (me == 0) {
+    cell0.begin()->pos() = {1.0, 2.0, 3.0};
+    cell0.begin()->id() = 0;
+  } else {
+    cell1.begin()->pos() = {4.0, 5.0, 6.0};
+    cell1.begin()->id() = 1;
+  }
+
+  HaloPlan plan;
+  plan.comm = world;
+  // cells[0] -> rank 0's cell, cells[1] -> rank 1's cell.
+  plan.collective =
+      CollectiveSection{CollectivePattern::Broadcast, {&cell0, &cell1}};
+
+  // --- Broadcast (Push): every rank's owned cell is broadcast to all. ---
+  halo_exchange(plan, box, GHOSTTRANS_POSITION,
+                {Direction::Push, Combine::Overwrite});
+
+  // After broadcast, each rank must have both particles.
+  BOOST_CHECK_CLOSE(cell0.begin()->pos()[0], 1.0, 1e-12);
+  BOOST_CHECK_CLOSE(cell0.begin()->pos()[1], 2.0, 1e-12);
+  BOOST_CHECK_CLOSE(cell0.begin()->pos()[2], 3.0, 1e-12);
+  BOOST_CHECK_CLOSE(cell1.begin()->pos()[0], 4.0, 1e-12);
+  BOOST_CHECK_CLOSE(cell1.begin()->pos()[1], 5.0, 1e-12);
+  BOOST_CHECK_CLOSE(cell1.begin()->pos()[2], 6.0, 1e-12);
+
+  // --- ReduceSum (Reduce): ghost forces are summed back to owner. ---
+  // Seed ghost forces: each rank seeds a force on the other's ghost cell.
+  if (me == 0) {
+    cell0.begin()->force() = {10.0, 0.0, 0.0}; // owned: pre-existing
+    cell1.begin()->force() = {0.0, 20.0, 0.0}; // ghost of rank1
+  } else {
+    cell0.begin()->force() = {0.0, 30.0, 0.0}; // ghost of rank0
+    cell1.begin()->force() = {40.0, 0.0, 0.0}; // owned: pre-existing
+  }
+
+  halo_exchange(plan, box, GHOSTTRANS_FORCE, {Direction::Reduce, Combine::Add});
+
+  // root=0: reduce sums rank0's cell0.force = {10,0,0} + rank1's cell0.force
+  // = {0,30,0} -> {10,30,0} on rank 0.
+  // root=1: reduce sums rank0's cell1.force = {0,20,0} + rank1's cell1.force
+  // = {40,0,0} -> {40,20,0} on rank 1.
+  // Only the root's copy of each cell is meaningful after the reduce.
+  if (me == 0) {
+    BOOST_CHECK_CLOSE(cell0.begin()->force()[0], 10.0, 1e-12);
+    BOOST_CHECK_CLOSE(cell0.begin()->force()[1], 30.0, 1e-12);
+  } else {
+    BOOST_CHECK_CLOSE(cell1.begin()->force()[0], 40.0, 1e-12);
+    BOOST_CHECK_CLOSE(cell1.begin()->force()[1], 20.0, 1e-12);
+  }
+}
+
+/*
+ * Same-rank (plan.local) path: Push direction.
+ *
+ * A LocalComm{src=real, dst=ghost, shift} represents a periodic-wrap copy on a
+ * node_grid==1 axis.  For Direction::Push the engine must call
+ * local_cell_copy(src, dst, shift) so the ghost receives the real particle's
+ * position, folded by the shift.  This test pins that semantics and would fail
+ * if the arguments were swapped (ghost would receive a zero-shifted copy of
+ * its own position instead of the real cell's shifted position).
+ *
+ * Runs on 1 rank (plan.local requires no MPI).
+ */
+BOOST_AUTO_TEST_CASE(local_push_position_shift,
+                     *utf::precondition([](utf::test_unit_id) {
+                       return boost::mpi::communicator{}.size() == 1;
+                     })) {
+  using namespace GhostComm;
+
+  BoxGeometry box;
+  box.set_length({10., 10., 10.});
+
+  ParticleList real_cell, ghost_cell;
+  real_cell.resize(1);
+  ghost_cell.resize(1);
+  real_cell.begin()->pos() = {3.0, 4.0, 5.0};
+  // Ghost starts at a distinct position so we know overwrite occurred.
+  ghost_cell.begin()->pos() = {0.0, 0.0, 0.0};
+
+  // shift of -10 on x: ghost should see folded position = 3 + (-10) = -7 ->
+  // folded inside [0,10) gives 3.0 (no folding needed for this range), but the
+  // raw shift is applied by the SAVE path so we check with the raw expected
+  // value.  Use shift == {0,0,0} to keep the position arithmetic trivial and
+  // focus on src/dst identity.
+  HaloPlan plan;
+  plan.comm = boost::mpi::communicator{};
+  plan.local.push_back(LocalComm{&real_cell, &ghost_cell, {0., 0., 0.}});
+
+  halo_exchange(plan, box, GHOSTTRANS_POSITION,
+                {Direction::Push, Combine::Overwrite});
+
+  // Ghost must now hold the real cell's position.
+  BOOST_CHECK_CLOSE(ghost_cell.begin()->pos()[0], 3.0, 1e-12);
+  BOOST_CHECK_CLOSE(ghost_cell.begin()->pos()[1], 4.0, 1e-12);
+  BOOST_CHECK_CLOSE(ghost_cell.begin()->pos()[2], 5.0, 1e-12);
+  // Real cell must be untouched.
+  BOOST_CHECK_CLOSE(real_cell.begin()->pos()[0], 3.0, 1e-12);
+  BOOST_CHECK_CLOSE(real_cell.begin()->pos()[1], 4.0, 1e-12);
+  BOOST_CHECK_CLOSE(real_cell.begin()->pos()[2], 5.0, 1e-12);
+}
+
+/*
+ * Same-rank (plan.local) path: Reduce direction — regression for commit
+ * 3c62b1353e ("core/ghosts: reduce self-copy ghost forces back to owner").
+ *
+ * On a node_grid==1 periodic axis the plan builder emits a LocalComm with
+ * src=real_cell, dst=ghost_cell.  During the force-reduce step the engine must
+ * ADD the ghost force into the real cell, i.e. call
+ *   local_cell_copy(*lc.dst, *lc.src, {}, box, GHOSTTRANS_FORCE)
+ * (dst=ghost is the new src, src=real is the new dst).
+ *
+ * The pre-fix code called local_cell_copy(*lc.src, *lc.dst, ...) for both
+ * directions, which silently dropped the ghost's force (it added the real
+ * cell's force back into the ghost instead of the other way around).
+ *
+ * Failure signature when the fix is reverted:
+ *   real_cell.force == {1,0,0} (unchanged — ghost force never reached it)
+ *   ghost_cell.force == {1,20,300} (real force accumulated into ghost instead)
+ *
+ * This test runs on 1 rank (plan.local requires no MPI) so it is always
+ * exercised, regardless of the MPI geometry used in the test suite.
+ */
+BOOST_AUTO_TEST_CASE(local_reduce_force_role_swap,
+                     *utf::precondition([](utf::test_unit_id) {
+                       return boost::mpi::communicator{}.size() == 1;
+                     })) {
+  using namespace GhostComm;
+
+  BoxGeometry box;
+  box.set_length({10., 10., 10.});
+
+  ParticleList real_cell, ghost_cell;
+  real_cell.resize(1);
+  ghost_cell.resize(1);
+  // Pre-existing force on the owned real particle.
+  real_cell.begin()->force() = {1., 0., 0.};
+  // Force accumulated on the ghost during the pair-interaction loop.
+  ghost_cell.begin()->force() = {0., 20., 300.};
+
+  HaloPlan plan;
+  plan.comm = boost::mpi::communicator{};
+  // LocalComm is always stored as {src=real, dst=ghost} by the plan builder.
+  plan.local.push_back(LocalComm{&real_cell, &ghost_cell, {0., 0., 0.}});
+
+  halo_exchange(plan, box, GHOSTTRANS_FORCE, {Direction::Reduce, Combine::Add});
+
+  // Real cell must hold pre-existing + ghost force.
+  BOOST_CHECK_CLOSE(real_cell.begin()->force()[0], 1.0, 1e-12);
+  BOOST_CHECK_CLOSE(real_cell.begin()->force()[1], 20.0, 1e-12);
+  BOOST_CHECK_CLOSE(real_cell.begin()->force()[2], 300.0, 1e-12);
+  // Ghost cell is NOT the target; its value is not checked here (the engine
+  // may or may not modify it, but the real cell result is what matters).
+}
+
+int main(int argc, char **argv) {
+  boost::mpi::environment mpi_env(argc, argv);
+
+  return boost::unit_test::unit_test_main(init_unit_test, argc, argv);
+}
diff --git a/src/core/unit_tests/HaloPlanValidator_test.cpp b/src/core/unit_tests/HaloPlanValidator_test.cpp
new file mode 100644
index 0000000000..bf86290618
--- /dev/null
+++ b/src/core/unit_tests/HaloPlanValidator_test.cpp
@@ -0,0 +1,393 @@
+/*
+ * 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/>.
+ */
+
+#define BOOST_TEST_MODULE HaloPlanValidator test
+#define BOOST_TEST_DYN_LINK
+#include "EspressoCoreGlobalConfig.hpp"
+
+#include "BoxGeometry.hpp"
+#include "LocalBox.hpp"
+#include "cell_system/AtomDecomposition.hpp"
+#include "cell_system/Cell.hpp"
+#include "cell_system/HybridDecomposition.hpp"
+#include "cell_system/RegularDecomposition.hpp"
+#include "communication.hpp"
+#include "ghosts/HaloPlan.hpp"
+#include "ghosts/HaloPlanValidator.hpp"
+#include "ghosts/mark_boundary_cells.hpp"
+
+#include <boost/mpi.hpp>
+#include <boost/test/unit_test.hpp>
+
+#include <functional>
+#include <optional>
+#include <set>
+#include <span>
+#include <vector>
+
+BOOST_TEST_GLOBAL_FIXTURE(EspressoCoreGlobalConfig);
+
+namespace {
+BoxGeometry make_periodic_box(double box_l) {
+  BoxGeometry box;
+  box.set_length(Utils::Vector3d::broadcast(box_l));
+  box.set_periodic(0u, true);
+  box.set_periodic(1u, true);
+  box.set_periodic(2u, true);
+  return box;
+}
+
+RegularDecomposition make_dd(Utils::Vector3i const &node_grid, double box_l,
+                             double range) {
+  ::communicator.set_node_grid(node_grid);
+  auto const box = make_periodic_box(box_l);
+  auto const local_box = LocalBox::make_regular_decomposition(
+      box.length(), ::communicator.calc_node_index(), ::communicator.node_grid);
+  return {::communicator.comm, range, box, local_box, std::nullopt};
+}
+
+// The AtomDecomposition (n-square) covers every ghost via its collective
+// broadcast/reduce section, so the validator must recognise collective-covered
+// ghosts as covered. The box lifetime must outlive the returned decomposition;
+// AtomDecomposition holds a BoxGeometry const& (see AtomDecomposition.hpp).
+AtomDecomposition make_atom_dd(BoxGeometry const &box) {
+  return {::communicator.comm, box};
+}
+
+// The HybridDecomposition combines the regular child's point-to-point section
+// with the n-square child's collective section, so its plan exercises BOTH
+// coverage paths at once.
+HybridDecomposition make_hybrid_dd(Utils::Vector3i const &node_grid,
+                                   BoxGeometry const &box,
+                                   double cutoff_regular, double skin,
+                                   std::set<int> const &n_square_types) {
+  ::communicator.set_node_grid(node_grid);
+  auto const local_box = LocalBox::make_regular_decomposition(
+      box.length(), ::communicator.calc_node_index(), ::communicator.node_grid);
+  return {
+      ::communicator.comm, cutoff_regular, skin, []() { return false; }, box,
+      local_box,           n_square_types};
+}
+} // namespace
+
+namespace utf = boost::unit_test;
+static bool has_4_mpi_ranks() { return boost::mpi::communicator{}.size() == 4; }
+static bool has_1_mpi_rank() { return boost::mpi::communicator{}.size() == 1; }
+
+// ── Task 2.1: single-rank local checks ──────────────────────────────────────
+// These use hand-made cells and are rank-independent, so they pass at NUM_PROC
+// 4 without modification.
+
+BOOST_AUTO_TEST_CASE(detects_defects) {
+  using namespace GhostComm;
+  // 1 local cell whose only ghost neighbor is `ghost`.
+  Cell local, ghost;
+  local.m_neighbors = Neighbors<Cell *>(std::vector<Cell *>{&ghost}, {});
+  std::vector<Cell *> locals{&local}, ghosts{&ghost};
+
+  // Mark interior/boundary so the consistency check in validate_halo_plan
+  // can correctly classify `local` as boundary (has a ghost neighbor).
+  mark_boundary_cells(locals, ghosts);
+
+  HaloPlan good;
+  good.neighbors.push_back(NeighborComm{
+      1, {SendRegion{&local.particles(), {}}}, {&ghost.particles()}});
+  BOOST_CHECK(validate_halo_plan(good, locals, ghosts).empty());
+
+  // uncovered ghost (neighborship-match + coverage failure)
+  HaloPlan uncovered; // empty plan, ghost never filled
+  BOOST_CHECK(!validate_halo_plan(uncovered, locals, ghosts).empty());
+
+  // duplicate peer
+  HaloPlan duppeer = good;
+  duppeer.neighbors.push_back(NeighborComm{
+      1, {SendRegion{&local.particles(), {}}}, {&ghost.particles()}});
+  BOOST_CHECK(!validate_halo_plan(duppeer, locals, ghosts).empty());
+
+  // double-filled ghost
+  HaloPlan dbl = good;
+  dbl.neighbors.push_back(NeighborComm{
+      2, {SendRegion{&local.particles(), {}}}, {&ghost.particles()}});
+  BOOST_CHECK(!validate_halo_plan(dbl, locals, ghosts).empty());
+
+  // send/recv size mismatch
+  HaloPlan badshape = good;
+  badshape.neighbors[0].recv.clear();
+  BOOST_CHECK(!validate_halo_plan(badshape, locals, ghosts).empty());
+
+  // recv target outside ghost set
+  Cell alien;
+  HaloPlan alien_recv = good;
+  alien_recv.neighbors[0].recv[0] = &alien.particles();
+  auto violations = validate_halo_plan(alien_recv, locals, ghosts);
+  BOOST_CHECK_MESSAGE(!violations.empty(),
+                      "Expected recv-outside-ghost violation");
+
+  // isolated shape mismatch (send.size() != recv.size())
+  HaloPlan shape_only = good;
+  shape_only.neighbors[0].send.push_back(SendRegion{&local.particles(), {}});
+  // send.size()==2, recv.size()==1 -> shape violation; ghost still covered once
+  auto shape_violations = validate_halo_plan(shape_only, locals, ghosts);
+  BOOST_CHECK_MESSAGE(!shape_violations.empty(),
+                      "Expected shape-mismatch violation");
+}
+
+// ── Task 2.3: real-plan local checks ────────────────────────────────────────
+// Build a real RegularDecomposition on 4 ranks and assert that the local
+// validator (coverage + neighborship-match + peer-uniqueness + shape) passes.
+// This is the guard that the wiring in RegularDecomposition.cpp won't
+// false-trip in CI (which builds with ESPRESSO_ADDITIONAL_CHECKS on).
+BOOST_AUTO_TEST_CASE(local_checks_pass_real_decomposition,
+                     *utf::precondition([](utf::test_unit_id) {
+                       return has_4_mpi_ranks();
+                     })) {
+  using namespace GhostComm;
+  auto const dd = make_dd({2, 2, 1}, 8., 1.2);
+  auto const *plan = dd.halo_plan();
+  BOOST_REQUIRE(plan != nullptr);
+  auto violations =
+      validate_halo_plan(*plan, dd.local_cells(), dd.ghost_cells());
+  BOOST_CHECK_MESSAGE(
+      violations.empty(),
+      "Expected no local-check violations for real decomposition");
+}
+
+// Single-rank real RegularDecomposition. make_halo_plan() returns an empty
+// plan on one rank (periodic neighbourships are wired directly between local
+// cells in init_cell_interactions), yet mark_cells() still classifies the
+// halo-layer cells as ghosts. Those ghosts are never referenced by a local
+// cell's neighbor stencil, so the validator must NOT flag them as uncovered.
+// This is the regression that aborted particle_reduction_test,
+// p3m_test_1_mpi_ranks, EspressoSystemStandAlone_test_1_mpi_ranks etc. under a
+// checks-on (RelWithAssert) build. Runs only when launched on a single rank.
+BOOST_AUTO_TEST_CASE(local_checks_pass_single_rank_decomposition,
+                     *utf::precondition([](utf::test_unit_id) {
+                       return has_1_mpi_rank();
+                     })) {
+  using namespace GhostComm;
+  auto const dd = make_dd({1, 1, 1}, 8., 1.2);
+  auto const *plan = dd.halo_plan();
+  BOOST_REQUIRE(plan != nullptr);
+  // Sanity: the single-rank plan is intentionally empty (no p2p, no local),
+  // but there are ghost cells present.
+  BOOST_CHECK(plan->neighbors.empty());
+  BOOST_CHECK(plan->local.empty());
+  BOOST_CHECK(not dd.ghost_cells().empty());
+  auto violations =
+      validate_halo_plan(*plan, dd.local_cells(), dd.ghost_cells());
+  BOOST_CHECK_MESSAGE(violations.empty(),
+                      "Expected no local-check violations for a single-rank "
+                      "RegularDecomposition (unreferenced halo-layer ghosts)");
+}
+
+// ── Fix 2 guard: collective-covered ghosts ──────────────────────────────────
+// A real AtomDecomposition (n-square) covers ALL of its ghosts via the
+// collective broadcast/reduce section, not via point-to-point recv/dst
+// targets. Before the fix the validator false-positived on every such ghost
+// ("never filled"). These cases run the validator directly in the checks-OFF
+// build, so they catch the false-positives without needing a checks-on build.
+BOOST_AUTO_TEST_CASE(atom_decomposition_collective_coverage,
+                     *utf::precondition([](utf::test_unit_id) {
+                       return has_4_mpi_ranks();
+                     })) {
+  using namespace GhostComm;
+  ::communicator.set_node_grid({4, 1, 1});
+  auto const box = make_periodic_box(8.);
+  auto const dd = make_atom_dd(box);
+  auto const *plan = dd.halo_plan();
+  BOOST_REQUIRE(plan != nullptr);
+  auto violations =
+      validate_halo_plan(*plan, dd.local_cells(), dd.ghost_cells());
+  BOOST_CHECK_MESSAGE(violations.empty(),
+                      "Expected no local-check violations for real "
+                      "AtomDecomposition (collective-covered ghosts)");
+  BOOST_CHECK(validate_halo_plan_symmetry(*plan).empty());
+}
+
+// A real HybridDecomposition covers the regular ghosts via point-to-point and
+// the n-square ghosts via the collective section — both paths at once.
+BOOST_AUTO_TEST_CASE(hybrid_decomposition_mixed_coverage,
+                     *utf::precondition([](utf::test_unit_id) {
+                       return has_4_mpi_ranks();
+                     })) {
+  using namespace GhostComm;
+  auto const box = make_periodic_box(8.);
+  auto const dd = make_hybrid_dd({2, 2, 1}, box, 1.2, 0.1, std::set<int>{0});
+  auto const *plan = dd.halo_plan();
+  BOOST_REQUIRE(plan != nullptr);
+  auto violations =
+      validate_halo_plan(*plan, dd.local_cells(), dd.ghost_cells());
+  BOOST_CHECK_MESSAGE(violations.empty(),
+                      "Expected no local-check violations for real "
+                      "HybridDecomposition (mixed p2p + collective coverage)");
+  BOOST_CHECK(validate_halo_plan_symmetry(*plan).empty());
+}
+
+// ── Task 5.2: wrap-aware boundary classification on 1 rank ──────────────────
+// On a single MPI rank all periodic axes are "wrap axes": no ghost cells are
+// generated for them (init_cell_interactions wires local cells directly across
+// the periodic boundary), so the base ghost-neighbour rule would leave those
+// cells interior.  The Task 5.2 wrap predicate must catch this and mark the
+// first and last local layers along each periodic axis as boundary.
+//
+// Correctness argument: every cell in the first or last local layer (ghost-grid
+// coord 1 or cell_grid[i] along axis i) has at least one neighbour in the
+// opposite boundary layer.  After mark_boundary_cells with the wrap predicate
+// those cells must be boundary.  Cells strictly in the interior (all ghost-grid
+// coords in [2, cell_grid[i]-1] along every axis) have no wrap neighbour and
+// must remain interior — PROVIDED the cell grid is large enough to have such
+// interior cells (cell_grid[i] >= 3 on every axis).
+BOOST_AUTO_TEST_CASE(single_rank_wrap_axis_cells_are_boundary,
+                     *utf::precondition([](utf::test_unit_id) {
+                       return has_1_mpi_rank();
+                     })) {
+  using namespace GhostComm;
+  // Use a box large enough for cell_grid >= 3 along every axis so that
+  // genuinely interior cells exist (they would be missed by the wrap rule).
+  // range=1.0 with box_l=6.0 gives cell_grid={6,6,6}.
+  auto const dd = make_dd({1, 1, 1}, 6., 1.0);
+
+  // Every cell in the first or last local layer along any axis must be
+  // boundary due to the periodic wrap.
+  auto const &gcg = dd.ghost_cell_grid;
+  auto const &cg = dd.cell_grid;
+
+  bool any_wrap_layer_interior = false;
+
+  for (Cell const *c : dd.local_cells()) {
+    // Compute ghost-grid coordinate.
+    int const idx = static_cast<int>(c - dd.cells.data());
+    int const gx = idx % gcg[0];
+    int const gy = (idx / gcg[0]) % gcg[1];
+    int const gz = idx / (gcg[0] * gcg[1]);
+
+    bool const in_wrap_layer = (gx == 1 || gx == cg[0] || gy == 1 ||
+                                gy == cg[1] || gz == 1 || gz == cg[2]);
+
+    if (in_wrap_layer && !c->is_boundary()) {
+      any_wrap_layer_interior = true;
+    }
+  }
+
+  BOOST_CHECK_MESSAGE(
+      !any_wrap_layer_interior,
+      "single-rank: all cells in first/last local layer must be boundary "
+      "(periodic-wrap rule)");
+
+  // Sanity: cells strictly in the interior must remain interior.
+  if (cg[0] >= 3 && cg[1] >= 3 && cg[2] >= 3) {
+    bool found_interior = false;
+    for (Cell const *c : dd.local_cells()) {
+      int const idx = static_cast<int>(c - dd.cells.data());
+      int const gx = idx % gcg[0];
+      int const gy = (idx / gcg[0]) % gcg[1];
+      int const gz = idx / (gcg[0] * gcg[1]);
+      bool const strictly_interior =
+          (gx >= 2 && gx <= cg[0] - 1 && gy >= 2 && gy <= cg[1] - 1 &&
+           gz >= 2 && gz <= cg[2] - 1);
+      if (strictly_interior) {
+        BOOST_CHECK_MESSAGE(!c->is_boundary(),
+                            "strictly interior cell must not be boundary");
+        found_interior = true;
+      }
+    }
+    BOOST_CHECK_MESSAGE(found_interior,
+                        "expected at least one strictly-interior cell "
+                        "(cell_grid >= 3 on all axes)");
+  }
+}
+
+// Degenerate wrap case: cell_grid==1 along at least one periodic axis.
+//
+// When node_grid[i]==1 && periodic[i] && cell_grid[i]==1 the single cell
+// layer on that axis is its own periodic neighbour.  init_cell_interactions()
+// drops the self-pair (ind1==ind2), so neighbors().all() has no entry along
+// that axis and the wrap predicate in mark_boundary_cells() is never invoked.
+// Without the post-predicate fixup every local cell would be left INTERIOR —
+// wrong, because each cell interacts with itself across the periodic boundary.
+//
+// box_l=1.0, range=2.0 → cell_grid[i]=floor(1/2)=0 → clamped to 1 on all
+// three axes.  Every local cell must be boundary.
+BOOST_AUTO_TEST_CASE(single_rank_cell_grid_one_all_cells_are_boundary,
+                     *utf::precondition([](utf::test_unit_id) {
+                       return has_1_mpi_rank();
+                     })) {
+  using namespace GhostComm;
+  // range > box_l → cell_grid clamps to {1,1,1}
+  auto const dd = make_dd({1, 1, 1}, 1.0, 2.0);
+  // Verify the precondition: cell_grid should be {1,1,1}.
+  BOOST_REQUIRE_EQUAL(dd.cell_grid[0], 1);
+  BOOST_REQUIRE_EQUAL(dd.cell_grid[1], 1);
+  BOOST_REQUIRE_EQUAL(dd.cell_grid[2], 1);
+  // Every local cell must be boundary (self-interaction across periodic
+  // boundary).
+  for (Cell const *c : dd.local_cells()) {
+    BOOST_CHECK_MESSAGE(c->is_boundary(),
+                        "cell_grid==1 wrap axis: every local cell must be "
+                        "boundary (self-interaction across periodic boundary)");
+  }
+}
+
+// ── Task 2.2: cross-rank symmetry checks ────────────────────────────────────
+
+// Good case: real RegularDecomposition on 4 ranks must be symmetric.
+BOOST_AUTO_TEST_CASE(symmetry_good_real_decomposition,
+                     *utf::precondition([](utf::test_unit_id) {
+                       return has_4_mpi_ranks();
+                     })) {
+  using namespace GhostComm;
+  auto const dd = make_dd({2, 2, 1}, 8., 1.2);
+  auto const *plan = dd.halo_plan();
+  BOOST_REQUIRE(plan != nullptr);
+  auto violations = validate_halo_plan_symmetry(*plan);
+  BOOST_CHECK_MESSAGE(violations.empty(),
+                      "Expected no symmetry violations for real decomposition");
+}
+
+// Asymmetric case: every rank posts a NeighborComm to peer=(me+1)%n with
+// send=[1 cell] and recv=[] (expects 0 from peer), but peer sends 1 to me.
+// Each rank should detect a mismatch for peer (me-1+n)%n.
+BOOST_AUTO_TEST_CASE(symmetry_detects_mismatch,
+                     *utf::precondition([](utf::test_unit_id) {
+                       return has_4_mpi_ranks();
+                     })) {
+  using namespace GhostComm;
+  boost::mpi::communicator comm;
+  int const n = comm.size();
+  int const me = comm.rank();
+
+  // Dummy cell for the send region.
+  Cell dummy;
+
+  HaloPlan asym;
+  asym.comm = comm;
+  // Every rank sends 1 item to (me+1)%n but declares recv=0 from that peer.
+  // So rank (me+1)%n will receive my send count=1 via all_to_all but see that
+  // my declared recv count from (me+1)%n is 0, which doesn't match what
+  // (me+1)%n actually sends.  Each rank therefore detects a mismatch.
+  int const send_peer = (me + 1) % n;
+  asym.neighbors.push_back(
+      NeighborComm{send_peer,
+                   {SendRegion{&dummy.particles(), {}}}, // send.size() == 1
+                   {}});                                 // recv.size() == 0
+
+  auto violations = validate_halo_plan_symmetry(asym);
+  BOOST_CHECK_MESSAGE(!violations.empty(),
+                      "Expected symmetry violation: recv=0 but peer sends 1");
+}
diff --git a/src/core/unit_tests/HaloPlan_test.cpp b/src/core/unit_tests/HaloPlan_test.cpp
new file mode 100644
index 0000000000..b0b9c94190
--- /dev/null
+++ b/src/core/unit_tests/HaloPlan_test.cpp
@@ -0,0 +1,36 @@
+/*
+ * Copyright (C) 2010-2026 The ESPResSo project
+ * Copyright (C) 2002,2003,2004,2005,2006,2007,2008,2009,2010
+ *   Max-Planck-Institute for Polymer Research, Theory Group
+ *
+ * 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/>.
+ */
+#define BOOST_TEST_MODULE HaloPlan test
+#define BOOST_TEST_DYN_LINK
+#include <boost/test/unit_test.hpp>
+
+#include "ParticleList.hpp"
+#include "ghosts/HaloPlan.hpp"
+
+BOOST_AUTO_TEST_CASE(neighbor_comm_shapes) {
+  using namespace GhostComm;
+  ParticleList a, b;
+  NeighborComm nc{/*peer*/ 1, {SendRegion{&a, {}}}, {&b}};
+  BOOST_CHECK_EQUAL(nc.peer, 1);
+  BOOST_CHECK_EQUAL(nc.send.size(), nc.recv.size()); // recv[k] <-> peer.send[k]
+  ExchangeOp op{Direction::Reduce, Combine::Add};
+  BOOST_CHECK(op.direction == Direction::Reduce);
+}
diff --git a/src/core/unit_tests/RegularHaloPlan_test.cpp b/src/core/unit_tests/RegularHaloPlan_test.cpp
new file mode 100644
index 0000000000..474ddabeeb
--- /dev/null
+++ b/src/core/unit_tests/RegularHaloPlan_test.cpp
@@ -0,0 +1,154 @@
+/*
+ * Copyright (C) 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/>.
+ */
+
+#define BOOST_TEST_MODULE RegularHaloPlan test
+#define BOOST_TEST_DYN_LINK
+#include <boost/test/unit_test.hpp>
+
+#include "EspressoCoreGlobalConfig.hpp"
+
+#include "BoxGeometry.hpp"
+#include "LocalBox.hpp"
+#include "ParticleList.hpp"
+#include "cell_system/RegularDecomposition.hpp"
+#include "communication.hpp"
+#include "ghosts/HaloPlan.hpp"
+
+#include <utils/Vector.hpp>
+
+#include <boost/mpi.hpp>
+
+#include <cstddef>
+#include <map>
+#include <optional>
+#include <set>
+#include <vector>
+
+namespace utf = boost::unit_test;
+
+BOOST_TEST_GLOBAL_FIXTURE(EspressoCoreGlobalConfig);
+
+// Only meaningful on 4 ranks: exercises face/edge/corner peers on a
+// non-trivial node grid.
+static bool has_4_mpi_ranks() { return boost::mpi::communicator{}.size() == 4; }
+
+namespace {
+/**
+ * @brief Build a RegularDecomposition directly on the global Cartesian
+ * communicator (mirrors the setup done by the cell-structure façade).
+ */
+RegularDecomposition make_dd(Utils::Vector3i const &node_grid, double box_l,
+                             double range) {
+  ::communicator.set_node_grid(node_grid);
+  BoxGeometry box;
+  box.set_length(Utils::Vector3d::broadcast(box_l));
+  box.set_periodic(0u, true);
+  box.set_periodic(1u, true);
+  box.set_periodic(2u, true);
+  auto const local_box = LocalBox::make_regular_decomposition(
+      box.length(), ::communicator.calc_node_index(), ::communicator.node_grid);
+  return {::communicator.comm, range, box, local_box, std::nullopt};
+}
+
+/**
+ * @brief Assert that every ghost cell is the target of exactly one recv/dst
+ * entry, and that the union of recv/dst targets equals the ghost-cell set.
+ */
+void check_coverage(RegularDecomposition const &dd) {
+  auto const *plan = dd.halo_plan();
+  BOOST_REQUIRE(plan != nullptr);
+
+  // Count how often each ParticleList is filled as a recv (peer) or dst
+  // (local) target.
+  std::map<ParticleList const *, unsigned> filled;
+  for (auto const &nc : plan->neighbors) {
+    // recv[k] must line up with peer.send[k].
+    BOOST_CHECK_EQUAL(nc.send.size(), nc.recv.size());
+    for (auto const *c : nc.recv) {
+      filled[c] += 1u;
+    }
+  }
+  for (auto const &lc : plan->local) {
+    filled[lc.dst] += 1u;
+  }
+
+  // Peers must be unique.
+  std::set<int> peers;
+  for (auto const &nc : plan->neighbors) {
+    BOOST_CHECK(peers.insert(nc.peer).second);
+    BOOST_CHECK(nc.peer != ::communicator.comm.rank());
+  }
+
+  // (a) no ghost cell is filled twice; (b) coverage: every ghost cell is
+  // filled exactly once, and nothing outside the ghost set is filled.
+  std::set<ParticleList const *> ghost_set;
+  for (auto const *g : dd.ghost_cells()) {
+    ghost_set.insert(&g->particles());
+    BOOST_CHECK_EQUAL(filled[&g->particles()], 1u);
+  }
+  // The union of recv/dst equals the ghost-cell set (no extras).
+  BOOST_CHECK_EQUAL(filled.size(), ghost_set.size());
+  for (auto const &[cell, count] : filled) {
+    BOOST_CHECK_EQUAL(count, 1u);
+    BOOST_CHECK(ghost_set.count(cell) == 1u);
+  }
+
+  // Cross-rank consistency: for the recv[k] <-> peer.send[k] contract to be
+  // satisfiable, my recv-from-P count must equal P's send-to-me count (and
+  // vice-versa). This is exactly what the async engine relies on for a-priori
+  // recv sizing, so verify it here rather than deferring entirely to 1.5.
+  auto comm = plan->comm;
+  std::map<int, std::size_t> my_recv, my_send;
+  for (auto const &nc : plan->neighbors) {
+    my_recv[nc.peer] = nc.recv.size();
+    my_send[nc.peer] = nc.send.size();
+  }
+  std::vector<boost::mpi::request> reqs;
+  std::map<int, std::size_t> peer_send, peer_recv;
+  for (auto const &nc : plan->neighbors) {
+    reqs.push_back(comm.isend(nc.peer, 0, my_send[nc.peer]));
+    reqs.push_back(comm.isend(nc.peer, 1, my_recv[nc.peer]));
+    reqs.push_back(comm.irecv(nc.peer, 0, peer_send[nc.peer]));
+    reqs.push_back(comm.irecv(nc.peer, 1, peer_recv[nc.peer]));
+  }
+  boost::mpi::wait_all(reqs.begin(), reqs.end());
+  for (auto const &nc : plan->neighbors) {
+    // My recv from P must equal P's send to me, and my send must equal P's
+    // recv.
+    BOOST_CHECK_EQUAL(my_recv[nc.peer], peer_send[nc.peer]);
+    BOOST_CHECK_EQUAL(my_send[nc.peer], peer_recv[nc.peer]);
+  }
+}
+} // namespace
+
+BOOST_AUTO_TEST_CASE(every_ghost_cell_covered_once_2x2x1,
+                     *utf::precondition([](utf::test_unit_id) {
+                       return has_4_mpi_ranks();
+                     })) {
+  auto const dd = make_dd({2, 2, 1}, 8., 1.2);
+  check_coverage(dd);
+}
+
+BOOST_AUTO_TEST_CASE(every_ghost_cell_covered_once_4x1x1,
+                     *utf::precondition([](utf::test_unit_id) {
+                       return has_4_mpi_ranks();
+                     })) {
+  auto const dd = make_dd({4, 1, 1}, 8., 1.2);
+  check_coverage(dd);
+}
diff --git a/src/core/virtual_sites/relative.cpp b/src/core/virtual_sites/relative.cpp
index 3ca277b34e..08493f77be 100644
--- a/src/core/virtual_sites/relative.cpp
+++ b/src/core/virtual_sites/relative.cpp
@@ -120,8 +120,15 @@ static bool is_vs(Particle const &p) {
 
 void vs_relative_update_particles(CellStructure &cell_structure,
                                   BoxGeometry const &box_geo) {
-  cell_structure.ghosts_update(Cells::DATA_PART_POSITION |
-                               Cells::DATA_PART_MOMENTUM);
+  cell_structure.ghosts_update(
+      Cells::DATA_PART_POSITION |
+      Cells::DATA_PART_MOMENTUM
+#ifdef ESPRESSO_ROTATION
+      // Orientation of the reference particle is needed
+      // for vs_relative position/orientation update
+      | Cells::DATA_PART_QUAT
+#endif
+  );
 
   cell_structure.for_each_local_particle([&](Particle &p) {
     if (!is_vs(p)) {
diff --git a/testsuite/python/caliper.py b/testsuite/python/caliper.py
index 09068d7830..f80f02d205 100644
--- a/testsuite/python/caliper.py
+++ b/testsuite/python/caliper.py
@@ -28,8 +28,26 @@
 
 EXPECTED_LABELS = """
 integrate
+  update_ghosts_and_resort_particle
+    resort_particles
+    ghosts_count
+      ghost/pack
+      ghost/post
+      ghost/wait
+      ghost/unpack
+    ghosts_update
+      ghost/pack
+      ghost/post
+      ghost/wait
+      ghost/unpack
   update_cabana_state
   Initial Force Calculation
+    update_ghosts_and_resort_particle
+      ghosts_update
+        ghost/pack
+        ghost/post
+        ghost/wait
+        ghost/unpack
     calculate_forces
       copy_particles_to_GPU
       update_cabana_state
@@ -37,9 +55,26 @@
       calc_long_range_forces
       cabana_short_range
       copy_forces_from_GPU
+      ghosts_reduce_forces
+        ghost/pack
+        ghost/post
+        ghost/wait
+        ghost/unpack
   Integration loop
     integrator_step_1
     resort_particles_if_needed
+    update_ghosts_and_resort_particle
+      resort_particles
+      ghosts_count
+        ghost/pack
+        ghost/post
+        ghost/wait
+        ghost/unpack
+      ghosts_update
+        ghost/pack
+        ghost/post
+        ghost/wait
+        ghost/unpack
     calculate_forces
       copy_particles_to_GPU
       update_cabana_state
@@ -47,8 +82,19 @@
       calc_long_range_forces
       cabana_short_range
       copy_forces_from_GPU
+      ghosts_reduce_forces
+        ghost/pack
+        ghost/post
+        ghost/wait
+        ghost/unpack
     integrator_step_2
 calc_energies
+  update_ghosts_and_resort_particle
+    ghosts_update
+      ghost/pack
+      ghost/post
+      ghost/wait
+      ghost/unpack
   update_cabana_state
   cabana_short_range
 """
@@ -76,7 +122,7 @@ def test_runtime_report(self):
         header = "Path\tMin time/rank\tAvg time/rank\tMax time/rank\tTime %"
         self.assertEqual(lines[0].split(), header.split(),
                          msg=f"Caliper summary should start with '{header}'")
-        labels = [line[:36].rstrip() for line in lines[1:]]
+        labels = [line[:38].rstrip() for line in lines[1:]]
         labels_ref = [x.rstrip() for x in EXPECTED_LABELS.strip().split("\n")
                       if x.rstrip() and ("GPU" not in x.upper() or has_cuda)]
         for l in labels_ref: assert l in labels, f"Label {l} not in profile"