Roll out new particle interface (#4460)

Description of changes:
- use the new `Particle` setters and getters in the core, except in electrostatics/magnetostatics
- hide `ext_flag` and `rotation` implementation details from the script interface
- make `lees_edwards` and `bond_breakage` submodules of `src/core` for consistency with `src/script_interface`
- make the bond breakage global variable static and unit test bond breakage actions
- cleanup comments (and clang-format tags) in source files modified by this rollout

src/core/AtomDecomposition.cpp

@@ -99,9 +99,9 @@ void AtomDecomposition::mark_cells() {
 void AtomDecomposition::resort(bool global_flag,
                                std::vector<ParticleChange> &diff) {
   for (auto &p : local().particles()) {
-    fold_position(p.r.p, p.l.i, m_box);
+    fold_position(p.pos(), p.image_box(), m_box);
 
-    p.l.p_old = p.r.p;
+    p.pos_at_last_verlet_update() = p.pos();
   }
 
   /* Local updates are a NoOp for this decomposition. */

src/core/BoxGeometry.hpp

@@ -19,8 +19,8 @@
 #ifndef CORE_BOX_GEOMETRY_HPP
 #define CORE_BOX_GEOMETRY_HPP
 
-#include "LeesEdwardsBC.hpp"
 #include "algorithm/periodic_fold.hpp"
+#include "lees_edwards/LeesEdwardsBC.hpp"
 
 #include <utils/Vector.hpp>
 #include <utils/math/sgn.hpp>
@@ -95,9 +95,9 @@ class BoxGeometry {
   /** Flags for all three dimensions whether pbc are applied (default). */
   std::bitset<3> m_periodic = 0b111;
   /** Side lengths of the box */
-  Utils::Vector3d m_length = {1, 1, 1};
+  Utils::Vector3d m_length = {1., 1., 1.};
   /** Inverse side lengths of the box */
-  Utils::Vector3d m_length_inv = {1, 1, 1};
+  Utils::Vector3d m_length_inv = {1., 1., 1.};
   /** Half side lengths of the box */
   Utils::Vector3d m_length_half = {0.5, 0.5, 0.5};
 

src/core/CMakeLists.txt

@@ -1,6 +1,5 @@
 set(EspressoCore_SRC
     accumulators.cpp
-    bond_breakage.cpp
     bond_error.cpp
     cells.cpp
     collision.cpp
@@ -90,6 +89,7 @@ target_compile_definitions(EspressoCore PUBLIC $<$<BOOL:${H5MD}>:H5XX_USE_MPI>)
 
 add_subdirectory(accumulators)
 add_subdirectory(actor)
+add_subdirectory(bond_breakage)
 add_subdirectory(bonded_interactions)
 add_subdirectory(cluster_analysis)
 add_subdirectory(constraints)

src/core/CellStructure.cpp

@@ -25,7 +25,7 @@
 #include "CellStructureType.hpp"
 #include "RegularDecomposition.hpp"
 #include "grid.hpp"
-#include "lees_edwards.hpp"
+#include "lees_edwards/lees_edwards.hpp"
 
 #include <utils/contains.hpp>
 

src/core/CellStructure.hpp

@@ -375,8 +375,9 @@ struct CellStructure {
                         Utils::Vector3d const &additional_offset = {}) const {
     auto const lim = Utils::sqr(skin / 2.) - additional_offset.norm2();
     return std::any_of(
-        particles.begin(), particles.end(),
-        [lim](const auto &p) { return ((p.r.p - p.l.p_old).norm2() > lim); });
+        particles.begin(), particles.end(), [lim](const auto &p) {
+          return ((p.pos() - p.pos_at_last_verlet_update()).norm2() > lim);
+        });
   }
 
   auto get_le_pos_offset_at_last_resort() const {
@@ -656,7 +657,7 @@ struct CellStructure {
   Cell *find_current_cell(const Particle &p) {
     assert(not get_resort_particles());
 
-    if (p.l.ghost) {
+    if (p.is_ghost()) {
       return nullptr;
     }
 

src/core/ComFixed.hpp

@@ -41,9 +41,9 @@ template <typename ParticleRange> class ComFixed {
 
     for (auto const &p : particles) {
       /* Check if type is of interest */
-      auto it = m_type_index.find(p.p.type);
+      auto it = m_type_index.find(p.type());
       if (it != m_type_index.end()) {
-        ret[it->second] += Utils::Vector3d{p.f.f};
+        ret[it->second] += p.force();
       }
     }
 
@@ -55,9 +55,9 @@ template <typename ParticleRange> class ComFixed {
 
     for (auto const &p : particles) {
       /* Check if type is of interest */
-      auto it = m_type_index.find(p.p.type);
+      auto it = m_type_index.find(p.type());
       if (it != m_type_index.end()) {
-        ret[it->second] += p.p.mass;
+        ret[it->second] += p.mass();
       }
     }
 
@@ -96,12 +96,12 @@ template <typename ParticleRange> class ComFixed {
 
     for (auto &p : particles) {
       /* Check if type is of interest */
-      auto it = m_type_index.find(p.p.type);
+      auto it = m_type_index.find(p.type());
       if (it != m_type_index.end()) {
-        auto const mass_frac = p.p.mass / masses[it->second];
+        auto const mass_frac = p.mass() / masses[it->second];
         auto const &type_force = forces[it->second];
         for (int i = 0; i < 3; i++) {
-          p.f.f[i] -= mass_frac * type_force[i];
+          p.force()[i] -= mass_frac * type_force[i];
         }
       }
     }

src/core/MpiCallbacks.hpp

@@ -46,6 +46,7 @@
 
 #include <boost/mpi/collectives/all_reduce.hpp>
 #include <boost/mpi/collectives/broadcast.hpp>
+#include <boost/mpi/collectives/reduce.hpp>
 #include <boost/mpi/communicator.hpp>
 #include <boost/optional.hpp>
 #include <boost/range/algorithm/remove_if.hpp>

src/core/Particle.hpp

@@ -29,7 +29,6 @@
 
 #include <boost/serialization/vector.hpp>
 
-#include <bitset>
 #include <cstdint>
 #include <stdexcept>
 #include <vector>
@@ -44,26 +43,6 @@ inline bool get_nth_bit(uint8_t const bitfield, int const bit_idx) {
 }
 } // namespace detail
 
-enum : uint8_t {
-  ROTATION_FIXED = 0u,
-  ROTATION_X = 1u,
-  ROTATION_Y = 2u,
-  ROTATION_Z = 4u
-};
-
-#ifdef EXTERNAL_FORCES
-/** \ref ParticleProperties::ext_flag "ext_flag" value for fixed coordinate
- *  @c coord.
- */
-#define COORD_FIXED(coord) (2u << (coord))
-/** \ref ParticleProperties::ext_flag "ext_flag" mask to check whether any of
- *  the coordinates is fixed.
- */
-#define COORDS_FIX_MASK (COORD_FIXED(0) | COORD_FIXED(1) | COORD_FIXED(2))
-#else // EXTERNAL_FORCES
-#define COORD_FIXED(coord) (0)
-#endif // EXTERNAL_FORCES
-
 /** Properties of a self-propelled particle. */
 struct ParticleParametersSwimming {
   /** Is the particle a swimmer. */
@@ -117,11 +96,18 @@ struct ParticleProperties {
   static constexpr Utils::Vector3d rinertia = {1., 1., 1.};
 #endif
 
-  /** bitfield for the particle axes of rotation */
 #ifdef ROTATION
-  uint8_t rotation = ROTATION_FIXED;
+  /** Bitfield for the particle axes of rotation.
+   *  Values:
+   *  - 1: allow rotation around the x axis
+   *  - 2: allow rotation around the y axis
+   *  - 4: allow rotation around the z axis
+   *  By default, the particle cannot rotate.
+   */
+  uint8_t rotation = static_cast<uint8_t>(0b000u);
 #else
-  static constexpr uint8_t rotation = ROTATION_FIXED;
+  /** Bitfield for the particle axes of rotation. Particle cannot rotate. */
+  static constexpr uint8_t rotation = static_cast<uint8_t>(0b000u);
 #endif
 
   /** charge. */
@@ -190,22 +176,22 @@ struct ParticleProperties {
 #endif // THERMOSTAT_PER_PARTICLE
 
 #ifdef EXTERNAL_FORCES
-  /** Flag for fixed particle coordinates.
+  /** Bitfield for fixed particle coordinates.
    *  Values:
-   *  - 0: no fixed coordinates
-   *  - 2: fix translation along the x axis
-   *  - 4: fix translation along the y axis
-   *  - 8: fix translation along the z axis
+   *  - 1: fix translation along the x axis
+   *  - 2: fix translation along the y axis
+   *  - 4: fix translation along the z axis
    */
-  uint8_t ext_flag = 0;
+  uint8_t ext_flag = static_cast<uint8_t>(0b000u);
   /** External force. */
   Utils::Vector3d ext_force = {0, 0, 0};
 #ifdef ROTATION
   /** External torque. */
   Utils::Vector3d ext_torque = {0, 0, 0};
 #endif
 #else  // EXTERNAL_FORCES
-  static constexpr const uint8_t ext_flag = 0; // no fixed coordinates
+  /** Bitfield for fixed particle coordinates. Coordinates cannot be fixed. */
+  static constexpr const uint8_t ext_flag = static_cast<uint8_t>(0b000u);
 #endif // EXTERNAL_FORCES
 
 #ifdef ENGINE
@@ -481,9 +467,7 @@ struct Particle { // NOLINT(bugprone-exception-escape)
   constexpr auto &mass() const { return p.mass; }
 #endif
 #ifdef ROTATION
-  bool can_rotate() const {
-    return can_rotate_around(0) or can_rotate_around(1) or can_rotate_around(2);
-  }
+  bool can_rotate() const { return static_cast<bool>(p.rotation); }
   bool can_rotate_around(int const axis) const {
     detail::check_axis_idx_valid(axis);
     return detail::get_nth_bit(p.rotation, axis);
@@ -496,10 +480,9 @@ struct Particle { // NOLINT(bugprone-exception-escape)
       p.rotation &= static_cast<uint8_t>(~(1u << axis));
     }
   }
-  void set_can_rotate_all_axes() {
-    for (int axis = 0; axis <= 2; axis++) {
-      set_can_rotate_around(axis, true);
-    }
+  void set_can_rotate_all_axes() { p.rotation = static_cast<uint8_t>(0b111u); }
+  void set_cannot_rotate_all_axes() {
+    p.rotation = static_cast<uint8_t>(0b000u);
   }
   auto const &quat() const { return r.quat; }
   auto &quat() { return r.quat; }
@@ -553,33 +536,25 @@ struct Particle { // NOLINT(bugprone-exception-escape)
 #endif // ROTATION
 #endif // THERMOSTAT_PER_PARTICLE
 #ifdef EXTERNAL_FORCES
-  bool has_fixed_coordinates() const {
-    return detail::get_nth_bit(p.ext_flag, 0);
-  }
+  bool has_fixed_coordinates() const { return static_cast<bool>(p.ext_flag); }
   bool is_fixed_along(int const axis) const {
     detail::check_axis_idx_valid(axis);
-    return detail::get_nth_bit(p.ext_flag, axis + 1);
+    return detail::get_nth_bit(p.ext_flag, axis);
   }
   void set_fixed_along(int const axis, bool const fixed_flag) {
     // set new flag
     if (fixed_flag) {
-      p.ext_flag |= static_cast<uint8_t>(1u << (axis + 1));
-    } else {
-      p.ext_flag &= static_cast<uint8_t>(~(1u << (axis + 1)));
-    }
-    // check if any flag is set and store that in the 0th bit
-    if (p.ext_flag >> 1) {
-      p.ext_flag |= static_cast<uint8_t>(1u);
+      p.ext_flag |= static_cast<uint8_t>(1u << axis);
     } else {
-      p.ext_flag &= static_cast<uint8_t>(~1u);
+      p.ext_flag &= static_cast<uint8_t>(~(1u << axis));
     }
   }
   auto const &ext_force() const { return p.ext_force; }
   auto &ext_force() { return p.ext_force; }
 
 #else  // EXTERNAL_FORCES
   constexpr bool has_fixed_coordinates() const { return false; }
-  constexpr bool is_fixed_along(int const axis) const { return false; }
+  constexpr bool is_fixed_along(int const) const { return false; }
 #endif // EXTERNAL_FORCES
 #ifdef ENGINE
   auto const &swimming() const { return p.swim; }

src/core/RegularDecomposition.cpp

@@ -80,7 +80,7 @@ void RegularDecomposition::move_if_local(
     ParticleList &src, ParticleList &rest,
     std::vector<ParticleChange> &modified_cells) {
   for (auto &part : src) {
-    auto target_cell = position_to_cell(part.r.p);
+    auto target_cell = position_to_cell(part.pos());
 
     if (target_cell) {
       target_cell->particles().insert(std::move(part));
@@ -162,9 +162,9 @@ namespace {
  * @brief Fold coordinates to box and reset the old position.
  */
 void fold_and_reset(Particle &p, BoxGeometry const &box_geo) {
-  fold_position(p.r.p, p.l.i, box_geo);
+  fold_position(p.pos(), p.image_box(), box_geo);
 
-  p.l.p_old = p.r.p;
+  p.pos_at_last_verlet_update() = p.pos();
 }
 } // namespace
 

src/core/RegularDecomposition.hpp

@@ -106,7 +106,7 @@ struct RegularDecomposition : public ParticleDecomposition {
   }
 
   Cell *particle_to_cell(Particle const &p) override {
-    return position_to_cell(p.r.p);
+    return position_to_cell(p.pos());
   }
 
   void resort(bool global, std::vector<ParticleChange> &diff) override;

src/core/bond_breakage/CMakeLists.txt

@@ -0,0 +1,2 @@
+target_sources(EspressoCore
+               PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/bond_breakage.cpp)

src/core/bond_breakage/actions.hpp

@@ -0,0 +1,78 @@
+/*
+ * Copyright (C) 2022 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/>.
+ */
+
+#ifndef CORE_BOND_BREAKAGE_ACTIONS_HPP
+#define CORE_BOND_BREAKAGE_ACTIONS_HPP
+
+#include <boost/functional/hash.hpp>
+#include <boost/variant.hpp>
+
+#include <cstddef>
+
+namespace BondBreakage {
+
+// Delete Actions
+struct DeleteBond {
+  int particle_id;
+  int bond_partner_id;
+  int bond_type;
+  std::size_t hash_value() const {
+    std::size_t seed = 3875;
+    boost::hash_combine(seed, particle_id);
+    boost::hash_combine(seed, bond_partner_id);
+    boost::hash_combine(seed, bond_type);
+    return seed;
+  }
+  bool operator==(DeleteBond const &rhs) const {
+    return rhs.particle_id == particle_id and
+           rhs.bond_partner_id == bond_partner_id and
+           rhs.bond_type == bond_type;
+  }
+};
+
+struct DeleteAllBonds {
+  int particle_id_1;
+  int particle_id_2;
+  std::size_t hash_value() const {
+    std::size_t seed = 75;
+    boost::hash_combine(seed, particle_id_1);
+    boost::hash_combine(seed, particle_id_2);
+    return seed;
+  }
+  bool operator==(DeleteAllBonds const &rhs) const {
+    return rhs.particle_id_1 == particle_id_1 and
+           rhs.particle_id_2 == particle_id_2;
+  }
+};
+} // namespace BondBreakage
+
+// Hash support for std::unordered_set
+namespace boost {
+template <> struct hash<BondBreakage::DeleteBond> {
+  std::size_t operator()(BondBreakage::DeleteBond const &t) const noexcept {
+    return t.hash_value();
+  }
+};
+template <> struct hash<BondBreakage::DeleteAllBonds> {
+  std::size_t operator()(BondBreakage::DeleteAllBonds const &t) const noexcept {
+    return t.hash_value();
+  }
+};
+} // namespace boost
+#endif