Refactor short range loop and cells (#3874)

Rewrite of #3787

Description of changes:
- simplify short range loop
   - remove `for_each_pair.hpp ` and `verlet_ia.hpp` helper functions
   - delegate distance calculation and cell iteration to the `CellStructure` class
- hide implementation details of `CellStructure`
   - make `local_cells()` private
   - encapsulate debug functions
   - remove global variables `n_verlet_updates` and `rebuild_verletlist`

src/core/AtomDecomposition.hpp

@@ -90,7 +90,9 @@ class AtomDecomposition : public ParticleDecomposition {
   Utils::Vector3d max_range() const override;
   /* Return true if minimum image convention is
    * needed for distance calculation. */
-  bool minimum_image_distance() const override { return true; }
+  boost::optional<BoxGeometry> minimum_image_distance() const override {
+    return m_box;
+  }
 
 private:
   /**

src/core/BoxGeometry.hpp

@@ -25,6 +25,7 @@
 
 #include <bitset>
 #include <cassert>
+#include <cmath>
 
 class BoxGeometry {
 public:
@@ -92,6 +93,9 @@ template <typename T> T get_mi_coord(T a, T b, T box_length, bool periodic) {
 
 /**
  * @brief Get the minimum-image vector between two coordinates.
+ *
+ * @tparam T Floating point type.
+ *
  * @param a     Coordinate of the terminal point.
  * @param b     Coordinate of the initial point.
  * @param box   Box parameters (side lengths, periodicity).

src/core/CMakeLists.txt

@@ -6,7 +6,6 @@ set(EspressoCore_SRC
     comfixed_global.cpp
     communication.cpp
     constraints.cpp
-    debug.cpp
     dpd.cpp
     energy.cpp
     errorhandling.cpp

src/core/CellStructure.cpp

@@ -26,6 +26,52 @@
 
 #include <utils/contains.hpp>
 
+#include <stdexcept>
+
+void CellStructure::check_particle_index() {
+  auto const max_id = get_max_local_particle_id();
+
+  for (auto const &p : local_particles()) {
+    auto const id = p.identity();
+
+    if (id < 0 || id > max_id) {
+      throw std::runtime_error("Particle id out of bounds.");
+    }
+
+    if (get_local_particle(id) != &p) {
+      throw std::runtime_error("Invalid local particle index entry.");
+    }
+  }
+
+  /* checks: local particle id */
+  int local_part_cnt = 0;
+  for (int n = 0; n < get_max_local_particle_id() + 1; n++) {
+    if (get_local_particle(n) != nullptr) {
+      local_part_cnt++;
+      if (get_local_particle(n)->p.identity != n) {
+        throw std::runtime_error("local_particles part has corrupted id.");
+      }
+    }
+  }
+
+  if (local_part_cnt != local_particles().size()) {
+    throw std::runtime_error(
+        std::to_string(local_particles().size()) + " parts in cells but " +
+        std::to_string(local_part_cnt) + " parts in local_particles");
+  }
+}
+
+void CellStructure::check_particle_sorting() {
+  for (auto cell : local_cells()) {
+    for (auto const &p : cell->particles()) {
+      if (particle_to_cell(p) != cell) {
+        throw std::runtime_error("misplaced particle with id " +
+                                 std::to_string(p.identity()));
+      }
+    }
+  }
+}
+
 Cell *CellStructure::particle_to_cell(const Particle &p) {
   return decomposition().particle_to_cell(p);
 }
@@ -166,6 +212,13 @@ void CellStructure::resort_particles(int global_flag) {
   for (auto d : diff) {
     boost::apply_visitor(UpdateParticleIndexVisitor{this}, d);
   }
+
+  m_rebuild_verlet_list = true;
+
+#ifdef ADDITIONAL_CHECKS
+  check_particle_index();
+  check_particle_sorting();
+#endif
 }
 
 void CellStructure::set_atom_decomposition(boost::mpi::communicator const &comm,
@@ -180,4 +233,4 @@ void CellStructure::set_domain_decomposition(
   set_particle_decomposition(
       std::make_unique<DomainDecomposition>(comm, range, box, local_geo));
   m_type = CELL_STRUCTURE_DOMDEC;
-}
+}

src/core/CellStructure.hpp

@@ -30,11 +30,13 @@
 #include "ParticleDecomposition.hpp"
 #include "ParticleList.hpp"
 #include "ParticleRange.hpp"
+#include "algorithm/link_cell.hpp"
 #include "bond_error.hpp"
 #include "ghosts.hpp"
 
 #include <boost/algorithm/cxx11/any_of.hpp>
 #include <boost/container/static_vector.hpp>
+#include <boost/iterator/indirect_iterator.hpp>
 #include <boost/range/algorithm/find_if.hpp>
 #include <boost/range/algorithm/transform.hpp>
 
@@ -80,6 +82,32 @@ inline ParticleRange particles(Utils::Span<Cell *> cells) {
 }
 } // namespace Cells
 
+/**
+ * @brief Distance vector and length handed to pair kernels.
+ */
+struct Distance {
+  explicit Distance(Utils::Vector3d const &vec21)
+      : vec21(vec21), dist2(vec21.norm2()) {}
+
+  Utils::Vector3d vec21;
+  double dist2;
+};
+namespace detail {
+struct MinimalImageDistance {
+  const BoxGeometry box;
+
+  Distance operator()(Particle const &p1, Particle const &p2) const {
+    return Distance(get_mi_vector(p1.r.p, p2.r.p, box));
+  }
+};
+
+struct EuclidianDistance {
+  Distance operator()(Particle const &p1, Particle const &p2) const {
+    return Distance(p1.r.p - p2.r.p);
+  }
+};
+} // namespace detail
+
 /** Describes a cell structure / cell system. Contains information
  *  about the communication of cell contents (particles, ghosts, ...)
  *  between different nodes and the relation between particle
@@ -103,6 +131,9 @@ struct CellStructure {
 public:
   bool use_verlet_list = true;
 
+  bool m_rebuild_verlet_list = true;
+  std::vector<std::pair<Particle *, Particle *>> m_verlet_list;
+
   /**
    * @brief Update local particle index.
    *
@@ -212,11 +243,14 @@ struct CellStructure {
   /** Maximal pair range supported by current cell system. */
   Utils::Vector3d max_range() const;
 
-  /** Return the global local_cells */
+private:
   Utils::Span<Cell *> local_cells();
+
+public:
   ParticleRange local_particles();
   ParticleRange ghost_particles();
 
+private:
   /** Cell system dependent function to find the right cell for a
    *  particle.
    *  \param  p Particle.
@@ -225,6 +259,7 @@ struct CellStructure {
    */
   Cell *particle_to_cell(const Particle &p);
 
+public:
   /**
    * @brief Add a particle.
    *
@@ -293,7 +328,6 @@ struct CellStructure {
     return assert(m_decomposition), *m_decomposition;
   }
 
-private:
   ParticleDecomposition &decomposition() {
     return assert(m_decomposition), *m_decomposition;
   }
@@ -327,6 +361,7 @@ struct CellStructure {
    * Cells::DataPart
    */
   void ghosts_update(unsigned data_parts);
+
   /**
    * @brief Add forces from ghost particles to real particles.
    */
@@ -355,7 +390,6 @@ struct CellStructure {
     return partners;
   }
 
-public:
   /**
    * @brief Execute kernel for every bond on particle.
    * @tparam Handler Callable, which can be invoked with
@@ -387,9 +421,10 @@ struct CellStructure {
     }
   }
 
-private:
-  /** Go through ghost cells and remove the ghost entries from the
-      local particle index. */
+  /**
+   * @brief Go through ghost cells and remove the ghost entries from the
+   * local particle index.
+   */
   void invalidate_ghosts() {
     for (auto const &p : ghost_particles()) {
       if (get_local_particle(p.identity()) == &p) {
@@ -443,11 +478,130 @@ struct CellStructure {
                                 double range, BoxGeometry const &box,
                                 LocalBox<double> const &local_geo);
 
+public:
+  template <class BondKernel> void bond_loop(BondKernel const &bond_kernel) {
+    for (auto &p : local_particles()) {
+      execute_bond_handler(p, bond_kernel);
+    }
+  }
+
+private:
   /**
-   * @brief Return true if minimum image convention is
-   *        needed for distance calculation. */
-  bool minimum_image_distance() const {
-    return m_decomposition->minimum_image_distance();
+   * @brief Run link_cell algorithm for local cells.
+   *
+   * @tparam Kernel Needs to be callable with (Particle, Particle, Distance).
+   * @param kernel Pair kernel functor.
+   */
+  template <class Kernel> void link_cell(Kernel kernel) {
+    auto const maybe_box = decomposition().minimum_image_distance();
+    auto const first = boost::make_indirect_iterator(local_cells().begin());
+    auto const last = boost::make_indirect_iterator(local_cells().end());
+
+    if (maybe_box) {
+      Algorithm::link_cell(
+          first, last,
+          [&kernel, df = detail::MinimalImageDistance{*maybe_box}](
+              Particle &p1, Particle &p2) { kernel(p1, p2, df(p1, p2)); });
+    } else {
+      Algorithm::link_cell(
+          first, last,
+          [&kernel, df = detail::EuclidianDistance{}](
+              Particle &p1, Particle &p2) { kernel(p1, p2, df(p1, p2)); });
+    }
+  }
+
+public:
+  /** Non-bonded pair loop with potential use
+   * of verlet lists.
+   * @param pair_kernel Kernel to apply
+   */
+  template <class PairKernel> void non_bonded_loop(PairKernel pair_kernel) {
+    link_cell(pair_kernel);
+  }
+
+  /** Non-bonded pair loop with potential use
+   * of verlet lists.
+   * @param pair_kernel Kernel to apply
+   * @param verlet_criterion Filter for verlet lists.
+   */
+  template <class PairKernel, class VerletCriterion>
+  void non_bonded_loop(PairKernel &&pair_kernel,
+                       const VerletCriterion &verlet_criterion) {
+    /* In this case the verlet list update is attached to
+     * the pair kernel, and the verlet list is rebuilt as
+     * we go. */
+    if (use_verlet_list && m_rebuild_verlet_list) {
+      m_verlet_list.clear();
+
+      link_cell([&pair_kernel, &verlet_criterion,
+                 this](Particle &p1, Particle &p2, Distance const &d) {
+        if (verlet_criterion(p1, p2, d)) {
+          m_verlet_list.emplace_back(&p1, &p2);
+          pair_kernel(p1, p2, d);
+        }
+      });
+
+      m_rebuild_verlet_list = false;
+    } else if (use_verlet_list && not m_rebuild_verlet_list) {
+      auto const maybe_box = decomposition().minimum_image_distance();
+      /* In this case the pair kernel is just run over the verlet list. */
+      if (maybe_box) {
+        auto const distance_function = detail::MinimalImageDistance{*maybe_box};
+        for (auto &pair : m_verlet_list) {
+          pair_kernel(*pair.first, *pair.second,
+                      distance_function(*pair.first, *pair.second));
+        }
+      } else {
+        auto const distance_function = detail::EuclidianDistance{};
+        for (auto &pair : m_verlet_list) {
+          pair_kernel(*pair.first, *pair.second,
+                      distance_function(*pair.first, *pair.second));
+        }
+      }
+    } else {
+      /* No verlet lists, just run the kernel with pairs from the cells. */
+      link_cell(pair_kernel);
+    }
+  }
+
+private:
+  /**
+   * @brief Check that particle index is commensurate with particles.
+   *
+   * For each local particles is checked that has a correct entry
+   * in the particles index, and that there are no excess (non-existing)
+   * particles in the index.
+   */
+  void check_particle_index();
+
+  /**
+   * @brief Check that particles are in the correct cell.
+   *
+   * This checks for all local particles that the result
+   * of particles_to_cell is the cell the particles is
+   * actually in, e.g. that the particles are sorted according
+   * to particles_to_cell.
+   */
+  void check_particle_sorting();
+
+public:
+  /**
+   * @brief Find cell a particle is stored in.
+   *
+   * For local particles, this returns the cell they
+   * are stored in, otherwise nullptr is returned.
+   *
+   * @param p Particle to find cell for
+   * @return Cell for particle or nullptr.
+   */
+  Cell *find_current_cell(const Particle &p) {
+    assert(not get_resort_particles());
+
+    if (p.l.ghost) {
+      return nullptr;
+    }
+
+    return particle_to_cell(p);
   }
 };
 

src/core/DomainDecomposition.cpp

@@ -224,6 +224,7 @@ void DomainDecomposition::resort(bool global,
     }
   }
 }
+
 void DomainDecomposition::mark_cells() {
   int cnt_c = 0;
 

src/core/DomainDecomposition.hpp

@@ -81,7 +81,6 @@ struct DomainDecomposition : public ParticleDecomposition {
                       const BoxGeometry &box_geo,
                       const LocalBox<double> &local_geo);
 
-public:
   GhostCommunicator const &exchange_ghosts_comm() const override {
     return m_exchange_ghosts_comm;
   }
@@ -100,10 +99,13 @@ struct DomainDecomposition : public ParticleDecomposition {
     return position_to_cell(p.r.p);
   }
 
-  bool minimum_image_distance() const override { return false; }
   void resort(bool global, std::vector<ParticleChange> &diff) override;
   Utils::Vector3d max_range() const override;
 
+  boost::optional<BoxGeometry> minimum_image_distance() const override {
+    return {};
+  }
+
 private:
   /** Fill local_cells list and ghost_cells list for use with domain
    *  decomposition.