diff --git a/src/core/Particle.hpp b/src/core/Particle.hpp
index a4c007617c4..425e518da27 100644
--- a/src/core/Particle.hpp
+++ b/src/core/Particle.hpp
@@ -467,8 +467,8 @@ struct Particle { // NOLINT(bugprone-exception-escape)
   constexpr auto &mass() const { return p.mass; }
 #endif
 #ifdef ROTATION
+  auto const &rotation() const { return p.rotation; }
   auto &rotation() { return p.rotation; }
-  auto rotation() const { return p.rotation; }
   bool can_rotate() const { return static_cast<bool>(p.rotation); }
   bool can_rotate_around(int const axis) const {
     detail::check_axis_idx_valid(axis);
@@ -498,8 +498,8 @@ struct Particle { // NOLINT(bugprone-exception-escape)
 #endif // EXTERNAL_FORCES
   auto calc_director() const { return r.calc_director(); }
 #else  // ROTATION
-  bool can_rotate() const { return false; }
-  bool can_rotate_around(int const axis) const { return false; }
+  auto can_rotate() const { return false; }
+  auto can_rotate_around(int const axis) const { return false; }
 #endif // ROTATION
 #ifdef DIPOLES
   auto const &dipm() const { return p.dipm; }
@@ -523,8 +523,9 @@ struct Particle { // NOLINT(bugprone-exception-escape)
   auto &mu_E() { return p.mu_E; }
 #endif
 #ifdef VIRTUAL_SITES
-  bool &virtual_flag() { return p.is_virtual; }
-  bool is_virtual() const { return p.is_virtual; }
+  auto &virtual_flag() { return p.is_virtual; }
+  auto const &virtual_flag() const { return p.is_virtual; }
+  auto is_virtual() const { return p.is_virtual; }
   void set_virtual(bool const virt_flag) { p.is_virtual = virt_flag; }
 #ifdef VIRTUAL_SITES_RELATIVE
   auto const &vs_relative() const { return p.vs_relative; }
@@ -542,6 +543,7 @@ struct Particle { // NOLINT(bugprone-exception-escape)
 #endif // ROTATION
 #endif // THERMOSTAT_PER_PARTICLE
 #ifdef EXTERNAL_FORCES
+  auto const &fixed() const { return p.ext_flag; }
   auto &fixed() { return p.ext_flag; }
   bool has_fixed_coordinates() const { return static_cast<bool>(p.ext_flag); }
   bool is_fixed_along(int const axis) const {
@@ -578,9 +580,6 @@ struct Particle { // NOLINT(bugprone-exception-escape)
 #ifdef EXCLUSIONS
   Utils::compact_vector<int> &exclusions() { return el; }
   Utils::compact_vector<int> const &exclusions() const { return el; }
-  std::vector<int> const exclusions_as_vector() const {
-    return {el.begin(), el.end()};
-  }
   bool has_exclusion(int pid) const {
     return std::find(el.begin(), el.end(), pid) != el.end();
   }
diff --git a/src/core/bond_breakage/bond_breakage.cpp b/src/core/bond_breakage/bond_breakage.cpp
index 22d6782bf87..caee1750393 100644
--- a/src/core/bond_breakage/bond_breakage.cpp
+++ b/src/core/bond_breakage/bond_breakage.cpp
@@ -22,6 +22,7 @@
 #include "cells.hpp"
 #include "communication.hpp"
 #include "errorhandling.hpp"
+#include "event.hpp"
 #include "particle_data.hpp"
 
 #include <utils/mpi/gather_buffer.hpp>
@@ -153,20 +154,45 @@ ActionSet actions_for_breakage(QueueEntry const &e) {
   return {};
 }
 
+/**
+ * @brief Delete specific bond.
+ */
+static void remove_bond(Particle &p, BondView const &view) {
+  auto &bond_list = p.bonds();
+  auto it = std::find(bond_list.begin(), bond_list.end(), view);
+  if (it != bond_list.end()) {
+    bond_list.erase(it);
+  }
+}
+
+/**
+ * @brief Delete pair bonds to a specific partner
+ */
+static void remove_pair_bonds_to(Particle &p, int other_pid) {
+  std::vector<std::pair<int, int>> to_delete;
+  for (auto b : p.bonds()) {
+    if (b.partner_ids().size() == 1 and b.partner_ids()[0] == other_pid)
+      to_delete.emplace_back(std::pair<int, int>{b.bond_id(), other_pid});
+  }
+  for (auto const &b : to_delete) {
+    remove_bond(p, BondView(b.first, {&b.second, 1}));
+  }
+}
+
 // Handler for the different delete events
 class execute : public boost::static_visitor<> {
 public:
   void operator()(DeleteBond const &d) const {
-    auto p = cell_structure.get_local_particle(d.particle_id);
-    if (!p)
-      return;
-    local_remove_bond(*p, {d.bond_type, d.bond_partner_id});
+    if (auto p = ::cell_structure.get_local_particle(d.particle_id)) {
+      remove_bond(*p, BondView(d.bond_type, {&d.bond_partner_id, 1}));
+    }
+    on_particle_change();
   }
   void operator()(DeleteAllBonds const &d) const {
-    auto p = cell_structure.get_local_particle(d.particle_id_1);
-    if (!p)
-      return;
-    local_remove_pair_bonds_to(*p, d.particle_id_2);
+    if (auto p = ::cell_structure.get_local_particle(d.particle_id_1)) {
+      remove_pair_bonds_to(*p, d.particle_id_2);
+    }
+    on_particle_change();
   }
 };
 
diff --git a/src/core/collision.cpp b/src/core/collision.cpp
index 14a2f8a267e..0c199d7d940 100644
--- a/src/core/collision.cpp
+++ b/src/core/collision.cpp
@@ -344,9 +344,7 @@ void place_vs_and_relate_to_particle(const int current_vs_pid,
   new_part.id() = current_vs_pid;
   new_part.pos() = pos;
   auto p_vs = cell_structure.add_particle(std::move(new_part));
-
-  local_vs_relate_to(*p_vs, get_part(relate_to));
-
+  vs_relate_to(*p_vs, get_part(relate_to));
   p_vs->set_virtual(true);
   p_vs->type() = collision_params.vs_particle_type;
 }
diff --git a/src/core/particle_data.cpp b/src/core/particle_data.cpp
index a55e449dcd2..f61db4ba881 100644
--- a/src/core/particle_data.cpp
+++ b/src/core/particle_data.cpp
@@ -25,28 +25,17 @@
 #include "cells.hpp"
 #include "communication.hpp"
 #include "event.hpp"
-#include "exclusions.hpp"
 #include "nonbonded_interactions/nonbonded_interaction_data.hpp"
-#include "partCfg_global.hpp"
 #include "particle_node.hpp"
-#include "rotation.hpp"
 
 #include "config/config.hpp"
 
-#include <utils/Span.hpp>
 #include <utils/Vector.hpp>
-#include <utils/quaternion.hpp>
 
-#include <boost/range/algorithm/copy.hpp>
-#include <boost/range/algorithm/find.hpp>
 #include <boost/serialization/variant.hpp>
-#include <boost/serialization/vector.hpp>
 #include <boost/variant.hpp>
 
-#include <iterator>
-#include <stdexcept>
-#include <unordered_map>
-#include <vector>
+#include <cassert>
 
 constexpr auto some_tag = 42;
 
@@ -90,168 +79,23 @@ using Prop = ParticleProperties;
 // clang-format off
 using UpdatePropertyMessage = boost::variant
         < UpdateProperty<int, &Prop::type>
-        , UpdateProperty<int, &Prop::mol_id>
-#ifdef MASS
-        , UpdateProperty<double, &Prop::mass>
-#endif
-#ifdef ROTATIONAL_INERTIA
-        , UpdateProperty<Utils::Vector3d, &Prop::rinertia>
-#endif
-#ifdef ROTATION
-        , UpdateProperty<uint8_t, &Prop::rotation>
-#endif
 #ifdef ELECTROSTATICS
         , UpdateProperty<double, &Prop::q>
 #endif
-#ifdef LB_ELECTROHYDRODYNAMICS
-        , UpdateProperty<Utils::Vector3d, &Prop::mu_E>
-#endif
-#ifdef ENGINE
-        , UpdateProperty<ParticleParametersSwimming, &Prop::swim>
-#endif
-#ifdef DIPOLES
-        , UpdateProperty<double, &Prop::dipm>
-#endif
-#ifdef VIRTUAL_SITES
-        , UpdateProperty<bool, &Prop::is_virtual>
-#ifdef VIRTUAL_SITES_RELATIVE
-        , UpdateProperty<ParticleProperties::VirtualSitesRelativeParameters,
-                         &Prop::vs_relative>
-#endif
-#endif
-#ifdef THERMOSTAT_PER_PARTICLE
-#ifndef PARTICLE_ANISOTROPY
-        , UpdateProperty<double, &Prop::gamma>
-#else
-        , UpdateProperty<Utils::Vector3d, &Prop::gamma>
-#endif // PARTICLE_ANISOTROPY
-#ifdef ROTATION
-#ifndef PARTICLE_ANISOTROPY
-        , UpdateProperty<double, &Prop::gamma_rot>
-#else
-        , UpdateProperty<Utils::Vector3d, &Prop::gamma_rot>
-#endif // PARTICLE_ANISOTROPY
-#endif // ROTATION
-#endif // THERMOSTAT_PER_PARTICLE
 #ifdef EXTERNAL_FORCES
-        , UpdateProperty<uint8_t, &Prop::ext_flag>
         , UpdateProperty<Utils::Vector3d, &Prop::ext_force>
-#ifdef ROTATION
-        , UpdateProperty<Utils::Vector3d, &Prop::ext_torque>
-#endif
 #endif
         >;
 
-using UpdateLocalPropertyMessage = boost::variant
-        < UpdateLocalProperty<double, &ParticleLocal::lees_edwards_offset>>;
-
 using UpdatePositionMessage = boost::variant
-        < UpdatePosition<Utils::Vector3d, &ParticlePosition::p>
-#ifdef ROTATION
-        , UpdatePosition<Utils::Quaternion<double>, &ParticlePosition::quat>
-#endif
-        >;
+        < UpdatePosition<Utils::Vector3d, &ParticlePosition::p> >;
 
 using UpdateMomentumMessage = boost::variant
-      < UpdateMomentum<Utils::Vector3d, &ParticleMomentum::v>
-#ifdef ROTATION
-      , UpdateMomentum<Utils::Vector3d, &ParticleMomentum::omega>
-#endif
-      >;
+      < UpdateMomentum<Utils::Vector3d, &ParticleMomentum::v> >;
 
 using UpdateForceMessage = boost::variant
-      < UpdateForce<Utils::Vector3d, &ParticleForce::f>
-#ifdef ROTATION
-      , UpdateForce<Utils::Vector3d, &ParticleForce::torque>
-#endif
-      >;
-// clang-format on
-
-/**
- * @brief Delete specific bond.
- */
-struct RemoveBond {
-  std::vector<int> bond;
-
-  void operator()(Particle &p) const {
-    assert(not bond.empty());
-    auto const view =
-        BondView(bond.front(), {bond.data() + 1, bond.size() - 1});
-    auto it = boost::find(p.bonds(), view);
-
-    if (it != p.bonds().end()) {
-      p.bonds().erase(it);
-    }
-  }
-
-  template <class Archive> void serialize(Archive &ar, long int) { ar &bond; }
-};
-
-/**
- * @brief Delete pair bonds to a specific partner
- */
-struct RemovePairBondsTo {
-  int other_pid;
-
-  void operator()(Particle &p) const {
-    using Bond = std::vector<int>;
-    std::vector<Bond> to_delete;
-    for (auto b : p.bonds()) {
-      if (b.partner_ids().size() == 1 and b.partner_ids()[0] == other_pid)
-        to_delete.push_back(Bond{b.bond_id(), other_pid});
-    }
-    for (auto b : to_delete) {
-      RemoveBond{b}(p);
-    }
-  }
-  template <class Archive> void serialize(Archive &ar, long int) {
-    ar &other_pid;
-  }
-};
-
-/**
- * @brief Delete all bonds.
- */
-struct RemoveBonds {
-  void operator()(Particle &p) const { p.bonds().clear(); }
-
-  template <class Archive> void serialize(Archive &, long int) {}
-};
-
-struct AddBond {
-  std::vector<int> bond;
-
-  void operator()(Particle &p) const {
-    auto const view = BondView(bond.at(0), {bond.data() + 1, bond.size() - 1});
-
-    p.bonds().insert(view);
-  }
-
-  template <class Archive> void serialize(Archive &ar, long int) { ar &bond; }
-};
-
-// clang-format off
-using UpdateBondMessage = boost::variant
-        < RemoveBond
-        , RemoveBonds
-        , AddBond
-        >;
-// clang-format on
-
-#ifdef ROTATION
-struct UpdateOrientation {
-  Utils::Vector3d axis;
-  double angle;
+      < UpdateForce<Utils::Vector3d, &ParticleForce::f> >;
 
-  void operator()(Particle &p) const { local_rotate_particle(p, axis, angle); }
-
-  template <class Archive> void serialize(Archive &ar, long int) {
-    ar &axis &angle;
-  }
-};
-#endif
-
-// clang-format off
 /**
  * @brief Top-level message.
  *
@@ -265,15 +109,10 @@ struct UpdateOrientation {
  * variants with leaves that have such an <tt>operator()</tt> member.
  */
 using UpdateMessage = boost::variant
-        < UpdateLocalPropertyMessage
-        , UpdatePropertyMessage
+        < UpdatePropertyMessage
         , UpdatePositionMessage
         , UpdateMomentumMessage
         , UpdateForceMessage
-        , UpdateBondMessage
-#ifdef ROTATION
-        , UpdateOrientation
-#endif
         >;
 // clang-format on
 
@@ -287,10 +126,6 @@ template <> struct message_type<ParticleProperties, &Particle::p> {
   using type = UpdatePropertyMessage;
 };
 
-template <> struct message_type<ParticleLocal, &Particle::l> {
-  using type = UpdateLocalPropertyMessage;
-};
-
 template <> struct message_type<ParticlePosition, &Particle::r> {
   using type = UpdatePositionMessage;
 };
@@ -335,14 +170,6 @@ struct UpdateVisitor : public boost::static_visitor<void> {
 };
 } // namespace
 
-void local_remove_bond(Particle &p, std::vector<int> const &bond) {
-  RemoveBond{bond}(p);
-}
-
-void local_remove_pair_bonds_to(Particle &p, int other_pid) {
-  RemovePairBondsTo{other_pid}(p);
-}
-
 static void mpi_send_update_message_local(int node, int id) {
   if (node == comm_cart.rank()) {
     UpdateMessage msg{};
@@ -407,102 +234,6 @@ void set_particle_v(int part, Utils::Vector3d const &v) {
                       &ParticleMomentum::v>(part, v);
 }
 
-void set_particle_lees_edwards_offset(int part, const double v) {
-  mpi_update_particle<ParticleLocal, &Particle::l, double,
-                      &ParticleLocal::lees_edwards_offset>(part, v);
-}
-
-#ifdef ENGINE
-void set_particle_swimming(int part, ParticleParametersSwimming swim) {
-  mpi_update_particle_property<ParticleParametersSwimming,
-                               &ParticleProperties::swim>(part, swim);
-}
-#endif
-
-void set_particle_f(int part, const Utils::Vector3d &f) {
-  mpi_update_particle<ParticleForce, &Particle::f, Utils::Vector3d,
-                      &ParticleForce::f>(part, f);
-}
-
-#ifdef MASS
-void set_particle_mass(int part, double mass) {
-  mpi_update_particle_property<double, &ParticleProperties::mass>(part, mass);
-}
-#else
-const constexpr double ParticleProperties::mass;
-#endif
-
-#ifdef ROTATIONAL_INERTIA
-void set_particle_rotational_inertia(int part,
-                                     Utils::Vector3d const &rinertia) {
-  mpi_update_particle_property<Utils::Vector3d, &ParticleProperties::rinertia>(
-      part, rinertia);
-}
-#else
-constexpr Utils::Vector3d ParticleProperties::rinertia;
-#endif
-
-#ifdef ROTATION
-void set_particle_rotation(int part, Utils::Vector3i const &flag) {
-  auto rot_flag = static_cast<uint8_t>(0u);
-  if (flag[0])
-    rot_flag |= static_cast<uint8_t>(1u);
-  if (flag[1])
-    rot_flag |= static_cast<uint8_t>(2u);
-  if (flag[2])
-    rot_flag |= static_cast<uint8_t>(4u);
-  mpi_update_particle_property<uint8_t, &ParticleProperties::rotation>(
-      part, rot_flag);
-}
-
-void rotate_particle(int part, const Utils::Vector3d &axis, double angle) {
-  mpi_send_update_message(part, UpdateOrientation{axis, angle});
-}
-#endif
-
-#ifdef DIPOLES
-void set_particle_dipm(int part, double dipm) {
-  mpi_update_particle_property<double, &ParticleProperties::dipm>(part, dipm);
-}
-
-void set_particle_dip(int part, Utils::Vector3d const &dip) {
-  auto const [quat, dipm] = convert_dip_to_quat(dip);
-  set_particle_dipm(part, dipm);
-  set_particle_quat(part, quat);
-}
-#endif
-
-#ifdef VIRTUAL_SITES
-void set_particle_virtual(int part, bool is_virtual) {
-  mpi_update_particle_property<bool, &ParticleProperties::is_virtual>(
-      part, is_virtual);
-}
-#endif
-
-#ifdef VIRTUAL_SITES_RELATIVE
-void set_particle_vs_quat(int part,
-                          Utils::Quaternion<double> const &vs_relative_quat) {
-  auto vs_relative = get_particle_data(part).vs_relative();
-  vs_relative.quat = vs_relative_quat;
-
-  mpi_update_particle_property<
-      ParticleProperties::VirtualSitesRelativeParameters,
-      &ParticleProperties::vs_relative>(part, vs_relative);
-}
-
-void set_particle_vs_relative(int part, int vs_relative_to, double vs_distance,
-                              Utils::Quaternion<double> const &rel_ori) {
-  ParticleProperties::VirtualSitesRelativeParameters vs_relative;
-  vs_relative.distance = vs_distance;
-  vs_relative.to_particle_id = vs_relative_to;
-  vs_relative.rel_orientation = rel_ori;
-
-  mpi_update_particle_property<
-      ParticleProperties::VirtualSitesRelativeParameters,
-      &ParticleProperties::vs_relative>(part, vs_relative);
-}
-#endif
-
 #ifdef ELECTROSTATICS
 void set_particle_q(int part, double q) {
   mpi_update_particle_property<double, &ParticleProperties::q>(part, q);
@@ -511,285 +242,8 @@ void set_particle_q(int part, double q) {
 const constexpr double ParticleProperties::q;
 #endif
 
-#ifdef LB_ELECTROHYDRODYNAMICS
-void set_particle_mu_E(int part, Utils::Vector3d const &mu_E) {
-  mpi_update_particle_property<Utils::Vector3d, &ParticleProperties::mu_E>(
-      part, mu_E);
-}
-#endif
-
 void set_particle_type(int p_id, int type) {
   make_particle_type_exist(type);
   on_particle_type_change(p_id, type);
   mpi_update_particle_property<int, &ParticleProperties::type>(p_id, type);
 }
-
-void set_particle_mol_id(int part, int mid) {
-  mpi_update_particle_property<int, &ParticleProperties::mol_id>(part, mid);
-}
-
-#ifdef ROTATION
-void set_particle_quat(int part, Utils::Quaternion<double> const &quat) {
-  mpi_update_particle<ParticlePosition, &Particle::r, Utils::Quaternion<double>,
-                      &ParticlePosition::quat>(part, quat);
-}
-
-void set_particle_director(int part, const Utils::Vector3d &director) {
-  Utils::Quaternion<double> quat =
-      convert_director_to_quaternion(director.normalized());
-  set_particle_quat(part, quat);
-}
-
-void set_particle_omega_lab(int part, const Utils::Vector3d &omega_lab) {
-  auto const &particle = get_particle_data(part);
-
-  mpi_update_particle<ParticleMomentum, &Particle::m, Utils::Vector3d,
-                      &ParticleMomentum::omega>(
-      part, convert_vector_space_to_body(particle, omega_lab));
-}
-
-void set_particle_omega_body(int part, const Utils::Vector3d &omega) {
-  mpi_update_particle<ParticleMomentum, &Particle::m, Utils::Vector3d,
-                      &ParticleMomentum::omega>(part, omega);
-}
-
-void set_particle_torque_lab(int part, const Utils::Vector3d &torque_lab) {
-  auto const &particle = get_particle_data(part);
-
-  mpi_update_particle<ParticleForce, &Particle::f, Utils::Vector3d,
-                      &ParticleForce::torque>(
-      part, convert_vector_space_to_body(particle, torque_lab));
-}
-#endif // ROTATION
-
-#ifdef THERMOSTAT_PER_PARTICLE
-#ifndef PARTICLE_ANISOTROPY
-void set_particle_gamma(int part, double gamma) {
-  mpi_update_particle_property<double, &ParticleProperties::gamma>(part, gamma);
-}
-#else
-void set_particle_gamma(int part, Utils::Vector3d const &gamma) {
-  mpi_update_particle_property<Utils::Vector3d, &ParticleProperties::gamma>(
-      part, gamma);
-}
-#endif // PARTICLE_ANISOTROPY
-
-#ifdef ROTATION
-#ifndef PARTICLE_ANISOTROPY
-void set_particle_gamma_rot(int part, double gamma_rot) {
-  mpi_update_particle_property<double, &ParticleProperties::gamma_rot>(
-      part, gamma_rot);
-}
-#else
-void set_particle_gamma_rot(int part, Utils::Vector3d const &gamma_rot) {
-  mpi_update_particle_property<Utils::Vector3d, &ParticleProperties::gamma_rot>(
-      part, gamma_rot);
-}
-#endif // PARTICLE_ANISOTROPY
-#endif // ROTATION
-#endif // THERMOSTAT_PER_PARTICLE
-
-#ifdef EXTERNAL_FORCES
-#ifdef ROTATION
-void set_particle_ext_torque(int part, const Utils::Vector3d &torque) {
-  mpi_update_particle_property<Utils::Vector3d,
-                               &ParticleProperties::ext_torque>(part, torque);
-}
-#endif // ROTATION
-
-void set_particle_ext_force(int part, const Utils::Vector3d &force) {
-  mpi_update_particle_property<Utils::Vector3d, &ParticleProperties::ext_force>(
-      part, force);
-}
-
-void set_particle_fix(int part, Utils::Vector3i const &flag) {
-  auto ext_flag = static_cast<uint8_t>(0u);
-  if (flag[0])
-    ext_flag |= static_cast<uint8_t>(1u);
-  if (flag[1])
-    ext_flag |= static_cast<uint8_t>(2u);
-  if (flag[2])
-    ext_flag |= static_cast<uint8_t>(4u);
-  mpi_update_particle_property<uint8_t, &ParticleProperties::ext_flag>(
-      part, ext_flag);
-}
-#endif // EXTERNAL_FORCES
-
-void delete_particle_bond(int part, Utils::Span<const int> bond) {
-  mpi_send_update_message(
-      part, UpdateBondMessage{RemoveBond{{bond.begin(), bond.end()}}});
-}
-
-void delete_particle_bonds(int part) {
-  mpi_send_update_message(part, UpdateBondMessage{RemoveBonds{}});
-}
-
-void add_particle_bond(int part, Utils::Span<const int> bond) {
-  mpi_send_update_message(
-      part, UpdateBondMessage{AddBond{{bond.begin(), bond.end()}}});
-}
-
-const std::vector<BondView> &get_particle_bonds(int part) {
-  static std::vector<BondView> ret;
-  ret.clear();
-
-  boost::copy(get_particle_data(part).bonds(), std::back_inserter(ret));
-
-  return ret;
-}
-
-void rescale_particles(int dir, double scale) {
-  for (auto &p : cell_structure.local_particles()) {
-    if (dir < 3)
-      p.pos()[dir] *= scale;
-    else {
-      p.pos() *= scale;
-    }
-  }
-  on_particle_change();
-}
-
-#ifdef EXCLUSIONS
-/**
- * @brief Locally remove an exclusion to a particle.
- * @param part1 the identity of the first exclusion partner
- * @param part2 the identity of the second exclusion partner
- */
-static void local_remove_exclusion(int part1, int part2) {
-  auto *p1 = cell_structure.get_local_particle(part1);
-  if (p1) {
-    delete_exclusion(*p1, part2);
-  }
-  auto *p2 = cell_structure.get_local_particle(part2);
-  if (p2) {
-    delete_exclusion(*p2, part1);
-  }
-}
-
-/**
- * @brief Locally add an exclusion to a particle.
- * @param part1 the identity of the first exclusion partner
- * @param part2 the identity of the second exclusion partner
- */
-static void local_add_exclusion(int part1, int part2) {
-  auto *p1 = cell_structure.get_local_particle(part1);
-  if (p1) {
-    add_exclusion(*p1, part2);
-  }
-  auto *p2 = cell_structure.get_local_particle(part2);
-  if (p2) {
-    add_exclusion(*p2, part1);
-  }
-}
-
-/* keep a unique list for particle i. Particle j is only added if it is not i
-   and not already in the list. */
-static void add_partner(std::vector<int> &il, int i, int j, int distance) {
-  if (j == i)
-    return;
-  for (int k = 0; k < il.size(); k += 2)
-    if (il[k] == j)
-      return;
-
-  il.push_back(j);
-  il.push_back(distance);
-}
-
-static void mpi_remove_exclusion_local(int part1, int part2) {
-  local_remove_exclusion(part1, part2);
-  on_particle_change();
-}
-
-REGISTER_CALLBACK(mpi_remove_exclusion_local)
-
-static void mpi_add_exclusion_local(int part1, int part2) {
-  local_add_exclusion(part1, part2);
-  on_particle_change();
-}
-
-REGISTER_CALLBACK(mpi_add_exclusion_local)
-
-static void check_particle_exists(int p_id) {
-  if (not particle_exists(p_id)) {
-    throw std::runtime_error("Particle with id " + std::to_string(p_id) +
-                             " not found");
-  }
-}
-
-static void particle_exclusion_sanity_checks(int part1, int part2) {
-  if (part1 == part2) {
-    throw std::runtime_error("Particles cannot exclude themselves (id " +
-                             std::to_string(part1) + ")");
-  }
-  check_particle_exists(part1);
-  check_particle_exists(part2);
-}
-
-void remove_particle_exclusion(int part1, int part2) {
-  particle_exclusion_sanity_checks(part1, part2);
-  mpi_call_all(mpi_remove_exclusion_local, part1, part2);
-}
-
-void add_particle_exclusion(int part1, int part2) {
-  particle_exclusion_sanity_checks(part1, part2);
-  mpi_call_all(mpi_add_exclusion_local, part1, part2);
-}
-
-void auto_exclusions(int distance) {
-  /* partners is a list containing the currently found excluded particles for
-     each particle, and their distance, as an interleaved list */
-  std::unordered_map<int, std::vector<int>> partners;
-
-  /* determine initial connectivity */
-  for (auto const &part1 : partCfg()) {
-    auto const p1 = part1.id();
-    for (auto const bond : part1.bonds()) {
-      if ((bond.partner_ids().size() == 1) and (bond.partner_ids()[0] != p1)) {
-        auto const p2 = bond.partner_ids()[0];
-        add_partner(partners[p1], p1, p2, 1);
-        add_partner(partners[p2], p2, p1, 1);
-      }
-    }
-  }
-
-  /* calculate transient connectivity. For each of the current neighbors,
-     also exclude their close enough neighbors.
-  */
-  for (int count = 1; count < distance; count++) {
-    for (auto const &p : partCfg()) {
-      auto const p1 = p.id();
-      for (int i = 0; i < partners[p1].size(); i += 2) {
-        auto const p2 = partners[p1][i];
-        auto const dist1 = partners[p1][i + 1];
-        if (dist1 > distance)
-          continue;
-        /* loop over all partners of the partner */
-        for (int j = 0; j < partners[p2].size(); j += 2) {
-          auto const p3 = partners[p2][j];
-          auto const dist2 = dist1 + partners[p2][j + 1];
-          if (dist2 > distance)
-            continue;
-          add_partner(partners[p1], p1, p3, dist2);
-          add_partner(partners[p3], p3, p1, dist2);
-        }
-      }
-    }
-  }
-
-  /* setup the exclusions and clear the arrays. We do not setup the exclusions
-     up there, since on_part_change clears the partCfg, so that we would have to
-     restore it continuously. Of course this could be optimized by bundling the
-     exclusions, but this is only done once and the overhead is as much as for
-     setting the bonds, which the user apparently accepted.
-  */
-  for (auto &kv : partners) {
-    auto const pid1 = kv.first;
-    auto const partner_list = kv.second;
-    for (int j = 0; j < partner_list.size(); j += 2) {
-      auto const pid2 = partner_list[j];
-      if (pid1 < pid2)
-        add_particle_exclusion(pid1, pid2);
-    }
-  }
-}
-#endif // EXCLUSIONS
diff --git a/src/core/particle_data.hpp b/src/core/particle_data.hpp
index 0418ce63340..0799203b009 100644
--- a/src/core/particle_data.hpp
+++ b/src/core/particle_data.hpp
@@ -28,13 +28,7 @@
 
 #include "config/config.hpp"
 
-#include "Particle.hpp"
-
-#include <utils/Span.hpp>
 #include <utils/Vector.hpp>
-#include <utils/quaternion.hpp>
-
-#include <vector>
 
 /** Call only on the head node: set particle velocity.
  *  @param part the particle.
@@ -42,60 +36,6 @@
  */
 void set_particle_v(int part, Utils::Vector3d const &v);
 
-/** Call only on the head node: set particle Lees-Edwards offset.
- *  @param part the particle.
- *  @param v new value for Lees-Edwards offset
- */
-void set_particle_lees_edwards_offset(int part, double v);
-
-#ifdef ENGINE
-/** Call only on the head node: set particle velocity.
- *  @param part the particle.
- *  @param swim struct containing swimming parameters
- */
-void set_particle_swimming(int part, ParticleParametersSwimming swim);
-#endif
-
-/** Call only on the head node: set particle force.
- *  @param part the particle.
- *  @param F its new force.
- */
-void set_particle_f(int part, const Utils::Vector3d &F);
-
-#ifdef MASS
-/** Call only on the head node: set particle mass.
- *  @param part the particle.
- *  @param mass its new mass.
- */
-void set_particle_mass(int part, double mass);
-#endif
-
-#ifdef ROTATIONAL_INERTIA
-/** Call only on the head node: set particle rotational inertia.
- *  @param part the particle.
- *  @param rinertia its new inertia.
- */
-void set_particle_rotational_inertia(int part, Utils::Vector3d const &rinertia);
-#endif
-
-#ifdef ROTATION
-/** Call only on the head node: Specifies whether a particle's rotational
- *  degrees of freedom are integrated or not. If set to zero, the content of
- *  the torque and omega variables are meaningless
- *  @param part the particle.
- *  @param flag the degrees of freedom flag.
- */
-void set_particle_rotation(int part, Utils::Vector3i const &flag);
-
-/** @brief rotate a particle around an axis
- *
- *  @param part particle id
- *  @param axis rotation axis
- *  @param angle rotation angle
- */
-void rotate_particle(int part, const Utils::Vector3d &axis, double angle);
-#endif
-
 #ifdef ELECTROSTATICS
 /** Call only on the head node: set particle charge.
  *  @param part the particle.
@@ -104,191 +44,10 @@ void rotate_particle(int part, const Utils::Vector3d &axis, double angle);
 void set_particle_q(int part, double q);
 #endif
 
-#ifdef LB_ELECTROHYDRODYNAMICS
-/** Call only on the head node: set particle electrophoretic mobility.
- *  @param part the particle.
- *  @param mu_E its new mobility.
- */
-void set_particle_mu_E(int part, Utils::Vector3d const &mu_E);
-#endif
-
 /** Call only on the head node: set particle type.
  *  @param p_id the particle.
  *  @param type its new type.
  */
 void set_particle_type(int p_id, int type);
 
-/** Call only on the head node: set particle's molecule id.
- *  @param part the particle.
- *  @param mid  its new mol id.
- */
-void set_particle_mol_id(int part, int mid);
-
-#ifdef ROTATION
-/** Call only on the head node: set particle orientation using quaternions.
- *  @param part the particle.
- *  @param quat its new value for quaternions.
- */
-void set_particle_quat(int part, Utils::Quaternion<double> const &quat);
-
-/** Call only on the head node: set particle orientation using director.
- *  The particle director defines the z-axis in the body-fixed frame.
- *  @param part the particle.
- *  @param director its new director vector (will be normalized if necessary)
- */
-void set_particle_director(int part, const Utils::Vector3d &director);
-
-/** Call only on the head node: set particle angular velocity from lab frame.
- *  @param part the particle.
- *  @param omega_lab its new angular velocity.
- */
-void set_particle_omega_lab(int part, const Utils::Vector3d &omega_lab);
-
-/** Call only on the head node: set particle angular velocity in body frame.
- *  @param part the particle.
- *  @param omega its new angular velocity.
- */
-void set_particle_omega_body(int part, const Utils::Vector3d &omega);
-
-/** Call only on the head node: set particle torque from lab frame.
- *  @param part the particle.
- *  @param torque_lab its new torque.
- */
-void set_particle_torque_lab(int part, const Utils::Vector3d &torque_lab);
-#endif
-
-#ifdef DIPOLES
-/** Call only on the head node: set particle dipole orientation.
- *  @param part the particle.
- *  @param dip its new dipole orientation.
- */
-void set_particle_dip(int part, Utils::Vector3d const &dip);
-
-/** Call only on the head node: set particle dipole moment (absolute value).
- *  @param part the particle.
- *  @param dipm its new dipole moment.
- */
-void set_particle_dipm(int part, double dipm);
-#endif
-
-#ifdef VIRTUAL_SITES
-/** Call only on the head node: set particle virtual flag.
- *  @param part the particle.
- *  @param is_virtual new @ref ParticleProperties::is_virtual "is_virtual" flag.
- */
-void set_particle_virtual(int part, bool is_virtual);
-#endif
-#ifdef VIRTUAL_SITES_RELATIVE
-void set_particle_vs_quat(int part,
-                          Utils::Quaternion<double> const &vs_relative_quat);
-void set_particle_vs_relative(int part, int vs_relative_to, double vs_distance,
-                              Utils::Quaternion<double> const &rel_ori);
-#endif
-
-#ifdef THERMOSTAT_PER_PARTICLE
-/** Call only on the head node: set particle frictional coefficient.
- *  @param part the particle.
- *  @param gamma its new frictional coefficient.
- */
-#ifndef PARTICLE_ANISOTROPY
-void set_particle_gamma(int part, double gamma);
-#else
-void set_particle_gamma(int part, Utils::Vector3d const &gamma);
-#endif
-#ifdef ROTATION
-#ifndef PARTICLE_ANISOTROPY
-void set_particle_gamma_rot(int part, double gamma);
-#else
-void set_particle_gamma_rot(int part, Utils::Vector3d const &gamma_rot);
-#endif
-#endif
-#endif // THERMOSTAT_PER_PARTICLE
-
-#ifdef EXTERNAL_FORCES
-#ifdef ROTATION
-/** Call only on the head node: set particle external torque.
- *  @param part  the particle.
- *  @param torque new value for ext_torque.
- */
-void set_particle_ext_torque(int part, const Utils::Vector3d &torque);
-#endif
-/** Call only on the head node: set particle external force.
- *  @param part  the particle.
- *  @param force new value for ext_force.
- */
-void set_particle_ext_force(int part, const Utils::Vector3d &force);
-/** Call only on the head node: set coordinate axes for which the particles
- *  motion is fixed.
- *  @param part  the particle.
- *  @param flag  coordinates to be fixed.
- */
-void set_particle_fix(int part, Utils::Vector3i const &flag);
-#endif // EXTERNAL_FORCES
-
-/** Call only on the head node: remove bond from particle.
- *  @param part     identity of principal atom of the bond.
- *  @param bond     field containing the bond type number and the identity
- *                  of all bond partners (secondary atoms of the bond).
- */
-void delete_particle_bond(int part, Utils::Span<const int> bond);
-
-/** Call only on the head node: remove all bonds from particle.
- *  @param part     identity of principal atom of the bond.
- */
-void delete_particle_bonds(int part);
-
-/** @brief Remove the specified bond from the particle
- *  @param p        The particle to update
- *  @param bond     The bond in the form
- *                  <tt>{bond_id, partner_1, partner_2, ...}</tt>
- */
-void local_remove_bond(Particle &p, std::vector<int> const &bond);
-
-/** @brief Remove all pair bonds on the particle which have the specified
- *  particle id as partner.
- *  @param p         The particle to update
- *  @param other_pid The particle id to filter for
- */
-void local_remove_pair_bonds_to(Particle &p, int other_pid);
-
-/** Call only on the head node: Add bond to particle.
- *  @param part     identity of principal atom of the bond.
- *  @param bond     field containing the bond type number and the identity
- *                  of all bond partners (secondary atoms of the bond).
- */
-void add_particle_bond(int part, Utils::Span<const int> bond);
-
-const std::vector<BondView> &get_particle_bonds(int part);
-
-#ifdef EXCLUSIONS
-/** @brief Remove particle exclusion.
- *  Call only on the head node.
- *  @param part1    identity of particle for which the exclusion is set.
- *  @param part2    identity of particle for which the exclusion is set.
- */
-void remove_particle_exclusion(int part1, int part2);
-
-/** @brief Add particle exclusion.
- *  Call only on the head node.
- *  @param part1    identity of particle for which the exclusion is set.
- *  @param part2    identity of particle for which the exclusion is set.
- */
-void add_particle_exclusion(int part1, int part2);
-
-/** Automatically add the next \<distance\> neighbors in each molecule to the
- *  exclusion list.
- *  This uses the bond topology obtained directly from the particles.
- *  To easily setup the bonds, all data should be on a single node,
- *  therefore the \ref partCfg array is used. With large amounts of
- *  particles, you should avoid this function and setup exclusions manually.
- */
-void auto_exclusions(int distance);
-#endif
-
-/** Rescale all particle positions in direction @p dir by a factor @p scale.
- *  @param dir   direction to scale (0/1/2 = x/y/z, 3 = x+y+z isotropically)
- *  @param scale factor by which to rescale (>1: stretch, <1: contract)
- */
-void rescale_particles(int dir, double scale);
-
 #endif
diff --git a/src/core/particle_node.cpp b/src/core/particle_node.cpp
index cd30c597787..00b99fb73ba 100644
--- a/src/core/particle_node.cpp
+++ b/src/core/particle_node.cpp
@@ -73,6 +73,18 @@ static std::unordered_map<int, int> particle_node;
  */
 static int max_seen_pid = -1;
 
+static auto rebuild_needed() {
+  auto is_rebuild_needed = ::particle_node.empty();
+  boost::mpi::broadcast(::comm_cart, is_rebuild_needed, 0);
+  return is_rebuild_needed;
+}
+
+static void mpi_synchronize_max_seen_pid_local() {
+  boost::mpi::broadcast(::comm_cart, ::max_seen_pid, 0);
+}
+
+REGISTER_CALLBACK(mpi_synchronize_max_seen_pid_local)
+
 void init_type_map(int type) {
   type_list_enable = true;
   if (type < 0)
@@ -98,8 +110,19 @@ static void add_id_to_type_map(int p_id, int type) {
     it->second.insert(p_id);
 }
 
+void on_particle_type_change_head(int p_id, int old_type, int new_type) {
+  if (type_list_enable) {
+    assert(::this_node == 0);
+    if (old_type != new_type) {
+      remove_id_from_map(p_id, old_type);
+    }
+    add_id_to_type_map(p_id, new_type);
+  }
+}
+
 void on_particle_type_change(int p_id, int type) {
   if (type_list_enable) {
+    assert(::this_node == 0);
     // check if the particle exists already and the type is changed, then remove
     // it from the list which contains it
     auto const &cur_par = get_particle_data(p_id);
@@ -254,8 +277,10 @@ static void mpi_who_has_local() {
   auto const n_part = static_cast<int>(local_particles.size());
   boost::mpi::gather(comm_cart, n_part, 0);
 
-  if (n_part == 0)
+  if (n_part == 0) {
+    mpi_synchronize_max_seen_pid_local();
     return;
+  }
 
   sendbuf.resize(n_part);
 
@@ -263,13 +288,12 @@ static void mpi_who_has_local() {
                  sendbuf.begin(), [](Particle const &p) { return p.id(); });
 
   comm_cart.send(0, some_tag, sendbuf);
+  mpi_synchronize_max_seen_pid_local();
 }
 
 REGISTER_CALLBACK(mpi_who_has_local)
 
-static void mpi_who_has() {
-  mpi_call(mpi_who_has_local);
-
+static void mpi_who_has_head() {
   auto local_particles = cell_structure.local_particles();
 
   static std::vector<int> n_parts;
@@ -295,12 +319,27 @@ static void mpi_who_has() {
       }
     }
   }
+  mpi_synchronize_max_seen_pid_local();
 }
 
 /**
  * @brief Rebuild the particle index.
  */
-static void build_particle_node() { mpi_who_has(); }
+static void build_particle_node() {
+  mpi_call(mpi_who_has_local);
+  mpi_who_has_head();
+}
+
+/**
+ * @brief Rebuild the particle index.
+ */
+static void build_particle_node_parallel() {
+  if (this_node == 0) {
+    mpi_who_has_head();
+  } else {
+    mpi_who_has_local();
+  }
+}
 
 int get_particle_node(int p_id) {
   if (p_id < 0) {
@@ -320,7 +359,10 @@ int get_particle_node(int p_id) {
   return needle->second;
 }
 
-void clear_particle_node() { particle_node.clear(); }
+void clear_particle_node() {
+  ::max_seen_pid = -1;
+  particle_node.clear();
+}
 
 static void clear_particle_type_map() {
   for (auto &kv : ::particle_type_map) {
@@ -381,7 +423,7 @@ static bool maybe_move_particle(int p_id, Utils::Vector3d const &pos) {
   return true;
 }
 
-static void particle_checks(int p_id, Utils::Vector3d const &pos) {
+void particle_checks(int p_id, Utils::Vector3d const &pos) {
   if (p_id < 0) {
     throw std::domain_error("Invalid particle id: " + std::to_string(p_id));
   }
@@ -396,6 +438,7 @@ static void particle_checks(int p_id, Utils::Vector3d const &pos) {
 static void mpi_remove_particle_local(int p_id) {
   cell_structure.remove_particle(p_id);
   on_particle_change();
+  mpi_synchronize_max_seen_pid_local();
 }
 
 REGISTER_CALLBACK(mpi_remove_particle_local)
@@ -403,15 +446,13 @@ REGISTER_CALLBACK(mpi_remove_particle_local)
 static void mpi_remove_all_particles_local() {
   cell_structure.remove_all_particles();
   on_particle_change();
+  clear_particle_node();
+  clear_particle_type_map();
 }
 
 REGISTER_CALLBACK(mpi_remove_all_particles_local)
 
-void remove_all_particles() {
-  mpi_call_all(mpi_remove_all_particles_local);
-  clear_particle_node();
-  clear_particle_type_map();
-}
+void remove_all_particles() { mpi_call_all(mpi_remove_all_particles_local); }
 
 void remove_particle(int p_id) {
   auto const &cur_par = get_particle_data(p_id);
@@ -433,9 +474,55 @@ void remove_particle(int p_id) {
       max_seen_pid = calculate_max_seen_id();
     }
   }
+  mpi_call_all(mpi_synchronize_max_seen_pid_local);
+}
+
+void remove_particle_parallel(int p_id) {
+  {
+    auto track_particle_types = ::type_list_enable;
+    boost::mpi::broadcast(::comm_cart, track_particle_types, 0);
+    if (track_particle_types) {
+      auto p = ::cell_structure.get_local_particle(p_id);
+      auto p_type = -1;
+      if (p != nullptr and not p->is_ghost()) {
+        if (this_node == 0) {
+          p_type = p->type();
+        } else {
+          ::comm_cart.send(0, 42, p->type());
+        }
+      } else if (this_node == 0) {
+        ::comm_cart.recv(boost::mpi::any_source, 42, p_type);
+      }
+      if (this_node == 0) {
+        assert(p_type != -1);
+        remove_id_from_map(p_id, p_type);
+      }
+    }
+  }
+
+  if (this_node == 0) {
+    particle_node[p_id] = -1;
+  }
+  mpi_remove_particle_local(p_id);
+  if (this_node == 0) {
+    particle_node.erase(p_id);
+    if (p_id == ::max_seen_pid) {
+      --::max_seen_pid;
+      // if there is a gap (i.e. there is no particle with id max_seen_pid - 1,
+      // then the cached value is invalidated and has to be recomputed (slow)
+      if (particle_node.count(::max_seen_pid) == 0 or
+          particle_node[::max_seen_pid] == -1) {
+        ::max_seen_pid = calculate_max_seen_id();
+      }
+    }
+  }
+  mpi_synchronize_max_seen_pid_local();
 }
 
-static void mpi_make_new_particle_local(int p_id, Utils::Vector3d const &pos) {
+void mpi_make_new_particle_local(int p_id, Utils::Vector3d const &pos) {
+  if (rebuild_needed()) {
+    build_particle_node_parallel();
+  }
   auto const has_created = maybe_insert_particle(p_id, pos);
   on_particle_change();
 
@@ -447,6 +534,7 @@ static void mpi_make_new_particle_local(int p_id, Utils::Vector3d const &pos) {
     max_seen_pid = std::max(max_seen_pid, p_id);
     assert(not has_created or node == 0);
   }
+  mpi_synchronize_max_seen_pid_local();
 }
 
 REGISTER_CALLBACK(mpi_make_new_particle_local)
@@ -456,7 +544,7 @@ void mpi_make_new_particle(int p_id, Utils::Vector3d const &pos) {
   mpi_call_all(mpi_make_new_particle_local, p_id, pos);
 }
 
-static void mpi_set_particle_pos_local(int p_id, Utils::Vector3d const &pos) {
+void mpi_set_particle_pos_local(int p_id, Utils::Vector3d const &pos) {
   auto const has_moved = maybe_move_particle(p_id, pos);
   ::cell_structure.set_resort_particles(Cells::RESORT_GLOBAL);
   on_particle_change();
@@ -524,6 +612,14 @@ int get_maximal_particle_id() {
   return max_seen_pid;
 }
 
+int get_maximal_particle_id_parallel() {
+  if (rebuild_needed()) {
+    build_particle_node_parallel();
+  }
+
+  return max_seen_pid;
+}
+
 int get_n_part() {
   if (particle_node.empty())
     build_particle_node();
diff --git a/src/core/particle_node.hpp b/src/core/particle_node.hpp
index c6fe1b79909..a519d4e5d78 100644
--- a/src/core/particle_node.hpp
+++ b/src/core/particle_node.hpp
@@ -70,6 +70,8 @@ std::size_t fetch_cache_max_size();
  */
 void clear_particle_node();
 
+void mpi_set_particle_pos_local(int p_id, Utils::Vector3d const &pos);
+
 /**
  * @brief Create a new particle and attach it to a cell.
  * Also call @ref on_particle_change.
@@ -77,6 +79,8 @@ void clear_particle_node();
  * @param pos   The particle position.
  */
 void mpi_make_new_particle(int p_id, Utils::Vector3d const &pos);
+void mpi_make_new_particle_local(int p_id, Utils::Vector3d const &pos);
+void particle_checks(int p_id, Utils::Vector3d const &pos);
 
 /**
  * @brief Move particle to a new position.
@@ -91,12 +95,14 @@ void mpi_set_particle_pos(int p_id, Utils::Vector3d const &pos);
  *  @param p_id     identity of the particle to remove
  */
 void remove_particle(int p_id);
+void remove_particle_parallel(int p_id);
 
 /** Remove all particles. */
 void remove_all_particles();
 
 void init_type_map(int type);
 void on_particle_type_change(int p_id, int type);
+void on_particle_type_change_head(int p_id, int old_type, int new_type);
 
 /** Find a particle of given type and return its id */
 int get_random_p_id(int type, int random_index_in_type_map);
@@ -129,6 +135,7 @@ std::vector<int> get_particle_ids();
  * @brief Get maximal particle id.
  */
 int get_maximal_particle_id();
+int get_maximal_particle_id_parallel();
 
 /**
  * @brief Get number of particles.
diff --git a/src/core/unit_tests/EspressoSystemStandAlone_test.cpp b/src/core/unit_tests/EspressoSystemStandAlone_test.cpp
index f4a5bcfb065..06e303e624f 100644
--- a/src/core/unit_tests/EspressoSystemStandAlone_test.cpp
+++ b/src/core/unit_tests/EspressoSystemStandAlone_test.cpp
@@ -33,6 +33,7 @@ namespace utf = boost::unit_test;
 #include "bonded_interactions/bonded_interaction_utils.hpp"
 #include "bonded_interactions/fene.hpp"
 #include "bonded_interactions/harmonic.hpp"
+#include "cells.hpp"
 #include "communication.hpp"
 #include "electrostatics/p3m.hpp"
 #include "electrostatics/registration.hpp"
@@ -95,6 +96,22 @@ static void mpi_create_bonds(int harm_bond_id, int fene_bond_id) {
   mpi_call_all(mpi_create_bonds_local, harm_bond_id, fene_bond_id);
 }
 
+static void mpi_add_bond_local(int p_id, int bond_id,
+                               std::vector<int> const &partner_ids) {
+  auto p = ::cell_structure.get_local_particle(p_id);
+  if (p != nullptr and not p->is_ghost()) {
+    p->bonds().insert(BondView(bond_id, partner_ids));
+  }
+  on_particle_change();
+}
+
+REGISTER_CALLBACK(mpi_add_bond_local)
+
+static void mpi_add_bond(int p_id, int bond_id,
+                         std::vector<int> const &partner_ids) {
+  mpi_call_all(mpi_add_bond_local, p_id, bond_id, partner_ids);
+}
+
 static void mpi_remove_translational_motion_local() {
   Galilei{}.kill_particle_motion(false);
 }
@@ -274,8 +291,8 @@ BOOST_FIXTURE_TEST_CASE(espresso_system_stand_alone, ParticleFactory,
         *boost::get<HarmonicBond>(bonded_ia_params.at(harm_bond_id).get());
     auto const &fene_bond =
         *boost::get<FeneBond>(bonded_ia_params.at(fene_bond_id).get());
-    add_particle_bond(pid2, std::vector<int>{harm_bond_id, pid1});
-    add_particle_bond(pid2, std::vector<int>{fene_bond_id, pid3});
+    mpi_add_bond(pid2, harm_bond_id, {pid1});
+    mpi_add_bond(pid2, fene_bond_id, {pid3});
 
     // measure energies
     auto const obs_energy = mpi_calculate_energy();
diff --git a/src/core/unit_tests/Particle_test.cpp b/src/core/unit_tests/Particle_test.cpp
index 76214e2f7f1..0611b5bf414 100644
--- a/src/core/unit_tests/Particle_test.cpp
+++ b/src/core/unit_tests/Particle_test.cpp
@@ -83,10 +83,6 @@ BOOST_AUTO_TEST_CASE(serialization) {
 
   BOOST_CHECK(q.id() == p.id());
   BOOST_CHECK((*q.bonds().begin() == BondView{bond_id, bond_partners}));
-
-#ifdef EXCLUSIONS
-  BOOST_CHECK(q.exclusions_as_vector() == el);
-#endif
 }
 
 namespace Utils {
diff --git a/src/core/unit_tests/Verlet_list_test.cpp b/src/core/unit_tests/Verlet_list_test.cpp
index 9749a70f42d..660463ba75f 100644
--- a/src/core/unit_tests/Verlet_list_test.cpp
+++ b/src/core/unit_tests/Verlet_list_test.cpp
@@ -36,6 +36,7 @@ namespace bdata = boost::unit_test::data;
 #include "MpiCallbacks.hpp"
 #include "Particle.hpp"
 #include "ParticleFactory.hpp"
+#include "cells.hpp"
 #include "communication.hpp"
 #include "event.hpp"
 #include "integrate.hpp"
@@ -70,6 +71,18 @@ struct if_head_node {
   boost::mpi::communicator world;
 };
 
+#ifdef EXTERNAL_FORCES
+static void mpi_set_particle_ext_force_local(int p_id,
+                                             Utils::Vector3d const &ext_force) {
+  if (auto p = ::cell_structure.get_local_particle(p_id)) {
+    p->ext_force() = ext_force;
+  }
+  on_particle_change();
+}
+
+REGISTER_CALLBACK(mpi_set_particle_ext_force_local)
+#endif // EXTERNAL_FORCES
+
 namespace Testing {
 /**
  * Helper class to setup an integrator and particle properties such that the
@@ -103,7 +116,8 @@ struct : public IntegratorHelper {
     mpi_call_all(mpi_set_integrator_sd_local);
   }
   void set_particle_properties(int pid) const override {
-    set_particle_ext_force(pid, {20., 0., 0.});
+    mpi_call_all(mpi_set_particle_ext_force_local, pid,
+                 Utils::Vector3d{20., 0., 0.});
   }
   char const *name() const override { return "SteepestDescent"; }
 } steepest_descent;
diff --git a/src/core/virtual_sites.cpp b/src/core/virtual_sites.cpp
index e50ab2cfe85..2c66a768e6f 100644
--- a/src/core/virtual_sites.cpp
+++ b/src/core/virtual_sites.cpp
@@ -31,15 +31,13 @@
 #include "grid.hpp"
 #include "integrate.hpp"
 #include "nonbonded_interactions/nonbonded_interaction_data.hpp"
-#include "particle_data.hpp"
-#include "particle_node.hpp"
 
 #include <utils/Vector.hpp>
 #include <utils/math/quaternion.hpp>
+#include <utils/quaternion.hpp>
 
 #include <cmath>
 #include <cstdio>
-#include <stdexcept>
 #include <tuple>
 
 namespace {
@@ -56,11 +54,11 @@ void set_virtual_sites(std::shared_ptr<VirtualSites> const &v) {
 #ifdef VIRTUAL_SITES_RELATIVE
 
 /** Calculate the rotation quaternion and distance between two particles */
-static std::tuple<Utils::Quaternion<double>, double>
-calculate_vs_relate_to_params(Particle const &p_current,
+std::tuple<Utils::Quaternion<double>, double>
+calculate_vs_relate_to_params(Particle const &p_vs,
                               Particle const &p_relate_to) {
   // get the distance between the particles
-  Utils::Vector3d d = box_geo.get_mi_vector(p_current.pos(), p_relate_to.pos());
+  auto d = ::box_geo.get_mi_vector(p_vs.pos(), p_relate_to.pos());
 
   // Check if the distance between virtual and non-virtual particles is larger
   // than minimum global cutoff. If so, warn user.
@@ -76,6 +74,13 @@ calculate_vs_relate_to_params(Particle const &p_current,
         << "increase the minimum cutoff.";
   }
 
+  // If the distance between real & virtual particle is 0 we just set the
+  // relative orientation to {1 0 0 0}, as it is irrelevant but needs to be
+  // a valid quaternion
+  if (dist == 0.) {
+    return std::make_tuple(Utils::Quaternion<double>::identity(), dist);
+  }
+
   // Now, calculate the quaternions which specify the angle between
   // the director of the particle we relate to and the vector
   // (particle_we_relate_to - this_particle)
@@ -87,73 +92,40 @@ calculate_vs_relate_to_params(Particle const &p_current,
   //   = quat_(obtained from desired director)
   // Resolving this for the quat_(virtual particle)
 
-  Utils::Quaternion<double> quat;
-  // If the distance between real & virtual particle is 0 we just set the
-  // relative orientation to {1 0 0 0}, as it is irrelevant but needs to be
-  // a valid quaternion
-  if (dist == 0) {
-    quat = Utils::Quaternion<double>::identity();
-  } else {
-    d.normalize();
-
-    // Obtain quaternion from desired director
-    Utils::Quaternion<double> quat_director =
-        Utils::convert_director_to_quaternion(d);
-
-    // Define quaternion as described above
-    auto relate_to_quat = p_relate_to.quat();
-    quat = Utils::Quaternion<double>{
-        {{{Utils::dot(relate_to_quat, quat_director),
-           -quat_director[0] * relate_to_quat[1] +
-               quat_director[1] * relate_to_quat[0] +
-               quat_director[2] * relate_to_quat[3] -
-               quat_director[3] * relate_to_quat[2],
-           relate_to_quat[1] * quat_director[3] +
-               relate_to_quat[0] * quat_director[2] -
-               relate_to_quat[3] * quat_director[1] -
-               relate_to_quat[2] * quat_director[0],
-           quat_director[3] * relate_to_quat[0] -
-               relate_to_quat[3] * quat_director[0] +
-               relate_to_quat[2] * quat_director[1] -
-               relate_to_quat[1] * quat_director[2]}}}};
-    quat /= relate_to_quat.norm2();
-
-    // Verify result
-    Utils::Quaternion<double> qtemp = relate_to_quat * quat;
-    for (int i = 0; i < 4; i++)
-      if (fabs(qtemp[i] - quat_director[i]) > 1E-9)
-        fprintf(stderr, "vs_relate_to: component %d: %f instead of %f\n", i,
-                qtemp[i], quat_director[i]);
+  d.normalize();
+
+  // Obtain quaternion from desired director
+  Utils::Quaternion<double> quat_director =
+      Utils::convert_director_to_quaternion(d);
+
+  // Define quaternion as described above
+  auto relate_to_quat = p_relate_to.quat();
+  auto quat =
+      Utils::Quaternion<double>{{{{Utils::dot(relate_to_quat, quat_director),
+                                   -quat_director[0] * relate_to_quat[1] +
+                                       quat_director[1] * relate_to_quat[0] +
+                                       quat_director[2] * relate_to_quat[3] -
+                                       quat_director[3] * relate_to_quat[2],
+                                   relate_to_quat[1] * quat_director[3] +
+                                       relate_to_quat[0] * quat_director[2] -
+                                       relate_to_quat[3] * quat_director[1] -
+                                       relate_to_quat[2] * quat_director[0],
+                                   quat_director[3] * relate_to_quat[0] -
+                                       relate_to_quat[3] * quat_director[0] +
+                                       relate_to_quat[2] * quat_director[1] -
+                                       relate_to_quat[1] * quat_director[2]}}}};
+  quat /= relate_to_quat.norm2();
+
+  // Verify result
+  Utils::Quaternion<double> qtemp = relate_to_quat * quat;
+  for (int i = 0; i < 4; i++) {
+    if (fabs(qtemp[i] - quat_director[i]) > 1E-9) {
+      fprintf(stderr, "vs_relate_to: component %d: %f instead of %f\n", i,
+              qtemp[i], quat_director[i]);
+    }
   }
   return std::make_tuple(quat, dist);
 }
 
-void local_vs_relate_to(Particle &p_current, Particle const &p_relate_to) {
-  // Set the particle id of the particle we want to relate to, the distance
-  // and the relative orientation
-  p_current.vs_relative().to_particle_id = p_relate_to.id();
-  std::tie(p_current.vs_relative().rel_orientation,
-           p_current.vs_relative().distance) =
-      calculate_vs_relate_to_params(p_current, p_relate_to);
-}
-
-void vs_relate_to(int part_num, int relate_to) {
-  if (part_num == relate_to) {
-    throw std::invalid_argument("A virtual site cannot relate to itself");
-  }
-  // Get the data for the particle we act on and the one we want to relate
-  // it to.
-  auto const &p_current = get_particle_data(part_num);
-  auto const &p_relate_to = get_particle_data(relate_to);
-
-  auto const [quat, dist] =
-      calculate_vs_relate_to_params(p_current, p_relate_to);
-
-  // Set the particle id of the particle we want to relate to, the distance
-  // and the relative orientation
-  set_particle_vs_relative(part_num, relate_to, dist, quat);
-  set_particle_virtual(part_num, true);
-}
-
 #endif // VIRTUAL_SITES_RELATIVE
 #endif // VIRTUAL_SITES
diff --git a/src/core/virtual_sites.hpp b/src/core/virtual_sites.hpp
index 18a282a1723..bf540b3beef 100644
--- a/src/core/virtual_sites.hpp
+++ b/src/core/virtual_sites.hpp
@@ -36,16 +36,28 @@ void set_virtual_sites(std::shared_ptr<VirtualSites> const &v);
 
 #ifdef VIRTUAL_SITES_RELATIVE
 
-/** Setup the @ref ParticleProperties::vs_relative "vs_relative" of a particle
- *  so that the given virtual particle will follow the given real particle.
- */
-void vs_relate_to(int part_num, int relate_to);
+#include <utils/quaternion.hpp>
 
-/** Setup the @ref ParticleProperties::vs_relative "vs_relative" of a particle
- *  so that the given virtual particle will follow the given real particle.
- */
-void local_vs_relate_to(Particle &p_current, Particle const &p_relate_to);
+#include <tuple>
 
-#endif
-#endif
+std::tuple<Utils::Quaternion<double>, double>
+calculate_vs_relate_to_params(Particle const &p_current,
+                              Particle const &p_relate_to);
+
+/**
+ * @brief Setup a virtual site to track a real particle.
+ * @param[in,out] p_vs         Virtual site.
+ * @param[in]     p_relate_to  Real particle to follow.
+ */
+inline void vs_relate_to(Particle &p_vs, Particle const &p_relate_to) {
+  // Set the particle id of the particle we want to relate to, the distance
+  // and the relative orientation
+  auto &vs_relative = p_vs.vs_relative();
+  vs_relative.to_particle_id = p_relate_to.id();
+  std::tie(vs_relative.rel_orientation, vs_relative.distance) =
+      calculate_vs_relate_to_params(p_vs, p_relate_to);
+}
+
+#endif // VIRTUAL_SITES_RELATIVE
+#endif // VIRTUAL_SITES
 #endif
diff --git a/src/python/espressomd/particle_data.py b/src/python/espressomd/particle_data.py
index 9467dcb54c4..c1c0ac18fbe 100644
--- a/src/python/espressomd/particle_data.py
+++ b/src/python/espressomd/particle_data.py
@@ -375,11 +375,20 @@ class ParticleHandle(ScriptInterfaceHelper):
     """
 
     _so_name = "Particles::ParticleHandle"
-    _so_creation_policy = "LOCAL"
+    _so_creation_policy = "GLOBAL"
     _so_bind_methods = (
         "delete_all_bonds",
     )
 
+    # here we must redefine the script interface setters
+
+    def set_params(self, **kwargs):
+        for name, value in kwargs.items():
+            self.set_parameter(name, value)
+
+    def set_parameter(self, name, value):
+        return self.call_method("set_param_parallel", name=name, value=value)
+
     def remove(self):
         """
         Delete the particle.
@@ -971,11 +980,32 @@ class ParticleList(ScriptInterfaceHelper):
         --------
         :meth:`espressomd.particle_data.ParticleHandle.remove`
 
+    auto_exclusions()
+        Add exclusions between particles that are connected by pair bonds,
+        including virtual bonds. Angle and diheral bonds are ignored. The most
+        common use case for this method is to auto-exclude virtual sites.
+
+        Another use case is to exclude 1-2, 1-3 and optionally 1-4 non-nonded
+        interactions on polymer chains. This technique is commonly used in
+        atomistic molecular dynamics engines such as NAMD, AMBER or GROMACS,
+        where the short-range part of the potential energy surface is better
+        approximated with Fourier sums (using dihedral bonds) than with pair
+        potentials. Linear, branched and circular topologies are supported.
+
+        Requires feature ``EXCLUSIONS``.
+
+        Parameters
+        ----------
+        distance : :obj:`int`
+            Maximal length of a chain in unit of bonds. The topology
+            will be traversed recursively until the bond chain either
+            terminates or reaches that distance.
+
     """
     _so_name = "Particles::ParticleList"
-    _so_creation_policy = "LOCAL"
+    _so_creation_policy = "GLOBAL"
     _so_bind_methods = (
-        "clear",
+        "clear", "auto_exclusions"
     )
 
     def by_id(self, p_id):
diff --git a/src/python/espressomd/system.py b/src/python/espressomd/system.py
index f2124133fa9..0d39573cc6d 100644
--- a/src/python/espressomd/system.py
+++ b/src/python/espressomd/system.py
@@ -485,4 +485,4 @@ def auto_exclusions(self, distance):
 
         """
         assert_features("EXCLUSIONS")
-        self.call_method("auto_exclusions", distance=distance)
+        self.part.auto_exclusions(distance=distance)
diff --git a/src/script_interface/cell_system/CellSystem.cpp b/src/script_interface/cell_system/CellSystem.cpp
index 5d942bafd32..cfd6c7c1465 100644
--- a/src/script_interface/cell_system/CellSystem.cpp
+++ b/src/script_interface/cell_system/CellSystem.cpp
@@ -230,9 +230,7 @@ Variant CellSystem::do_call_method(std::string const &name,
 
 std::vector<int> CellSystem::mpi_resort_particles(bool global_flag) const {
   ::cell_structure.resort_particles(global_flag, box_geo);
-  if (context()->is_head_node()) {
-    clear_particle_node();
-  }
+  clear_particle_node();
   auto const size = static_cast<int>(::cell_structure.local_particles().size());
   std::vector<int> n_part_per_node;
   boost::mpi::gather(context()->get_comm(), size, n_part_per_node, 0);
diff --git a/src/script_interface/particle_data/ParticleHandle.cpp b/src/script_interface/particle_data/ParticleHandle.cpp
index a1b82b381fd..f9248d13a60 100644
--- a/src/script_interface/particle_data/ParticleHandle.cpp
+++ b/src/script_interface/particle_data/ParticleHandle.cpp
@@ -21,11 +21,16 @@
 
 #include "ParticleHandle.hpp"
 
+#include "script_interface/Variant.hpp"
+#include "script_interface/communication.hpp"
 #include "script_interface/get_value.hpp"
 
 #include "core/bonded_interactions/bonded_interaction_data.hpp"
+#include "core/cells.hpp"
+#include "core/event.hpp"
+#include "core/exclusions.hpp"
 #include "core/grid.hpp"
-#include "core/particle_data.hpp"
+#include "core/nonbonded_interactions/nonbonded_interaction_data.hpp"
 #include "core/particle_node.hpp"
 #include "core/rotation.hpp"
 #include "core/virtual_sites.hpp"
@@ -33,18 +38,37 @@
 #include <utils/Vector.hpp>
 
 #include <boost/format.hpp>
+#include <boost/mpi/collectives/all_reduce.hpp>
+#include <boost/mpi/collectives/broadcast.hpp>
+#include <boost/mpi/communicator.hpp>
 
 #include <algorithm>
+#include <cassert>
 #include <cmath>
+#include <cstddef>
 #include <memory>
+#include <set>
 #include <sstream>
 #include <stdexcept>
 #include <string>
+#include <tuple>
+#include <type_traits>
 #include <vector>
 
 namespace ScriptInterface {
 namespace Particles {
 
+static uint8_t bitfield_from_flag(Utils::Vector3i const &flag) {
+  auto bitfield = static_cast<uint8_t>(0u);
+  if (flag[0])
+    bitfield |= static_cast<uint8_t>(1u);
+  if (flag[1])
+    bitfield |= static_cast<uint8_t>(2u);
+  if (flag[2])
+    bitfield |= static_cast<uint8_t>(4u);
+  return bitfield;
+}
+
 static auto error_msg(std::string const &name, std::string const &reason) {
   std::stringstream msg;
   msg << "attribute '" << name << "' of 'ParticleHandle' " << reason;
@@ -77,15 +101,58 @@ static auto get_gamma_safe(Variant const &value) {
 }
 #endif // THERMOSTAT_PER_PARTICLE
 
-static auto get_bond_vector(VariantMap const &params) {
-  auto const bond_id = get_value<int>(params, "bond_id");
-  auto const part_id = get_value<std::vector<int>>(params, "part_id");
-  std::vector<int> bond_view;
-  bond_view.emplace_back(bond_id);
-  for (auto const pid : part_id) {
-    bond_view.emplace_back(pid);
+static auto get_real_particle(boost::mpi::communicator const &comm, int p_id) {
+  if (p_id < 0) {
+    throw std::domain_error("Invalid particle id: " + std::to_string(p_id));
+  }
+  auto ptr = ::cell_structure.get_local_particle(p_id);
+  if (ptr != nullptr and ptr->is_ghost()) {
+    ptr = nullptr;
+  }
+  auto const n_found = boost::mpi::all_reduce(
+      comm, static_cast<int>(ptr != nullptr), std::plus<>());
+  if (n_found == 0) {
+    throw std::runtime_error("Particle with id " + std::to_string(p_id) +
+                             " not found");
+  }
+  return ptr;
+}
+
+template <typename T, class F>
+T ParticleHandle::get_particle_property(F const &fun) const {
+  auto const &comm = context()->get_comm();
+  auto const ptr = const_cast<Particle const *>(get_real_particle(comm, m_pid));
+  boost::optional<T> ret;
+  if (ptr == nullptr) {
+    ret = {};
+  } else {
+    ret = {fun(*ptr)};
+  }
+  return mpi_reduce_optional(comm, ret);
+}
+
+template <typename T>
+T ParticleHandle::get_particle_property(T const &(Particle::*getter)()
+                                            const) const {
+  return get_particle_property<T>(
+      [getter](Particle const &p) { return (p.*getter)(); });
+}
+
+template <class F>
+void ParticleHandle::set_particle_property(F const &fun) const {
+  auto const &comm = context()->get_comm();
+  auto const ptr = get_real_particle(comm, m_pid);
+  if (ptr != nullptr) {
+    fun(*ptr);
   }
-  return bond_view;
+  on_particle_change();
+}
+
+template <typename T>
+void ParticleHandle::set_particle_property(T &(Particle::*setter)(),
+                                           Variant const &value) const {
+  set_particle_property(
+      [&value, setter](Particle &p) { (p.*setter)() = get_value<T>(value); });
 }
 
 ParticleHandle::ParticleHandle() {
@@ -93,49 +160,57 @@ ParticleHandle::ParticleHandle() {
       {"id", AutoParameter::read_only, [this]() { return m_pid; }},
       {"type",
        [this](Variant const &value) {
-         auto const p_type = get_value<int>(value);
-         if (p_type < 0) {
+         auto const old_type = get_particle_property(&Particle::type);
+         auto const new_type = get_value<int>(value);
+         if (new_type < 0) {
            throw std::domain_error(
                error_msg("type", "must be an integer >= 0"));
          }
-         set_particle_type(m_pid, p_type);
+         make_particle_type_exist_local(new_type);
+         if (context()->is_head_node()) {
+           on_particle_type_change_head(m_pid, old_type, new_type);
+         }
+         set_particle_property(&Particle::type, value);
        },
-       [this]() { return particle().type(); }},
+       [this]() { return get_particle_data(m_pid).type(); }},
       {"pos",
        [this](Variant const &value) {
-         mpi_set_particle_pos(m_pid, get_value<Utils::Vector3d>(value));
+         mpi_set_particle_pos_local(m_pid, get_value<Utils::Vector3d>(value));
        },
        [this]() {
-         auto const &p = particle();
-         return unfolded_position(p.pos(), p.image_box(), ::box_geo.length());
+         auto const p = get_particle_data(m_pid);
+         auto const pos = p.pos();
+         auto const image_box = p.image_box();
+         return unfolded_position(pos, image_box, ::box_geo.length());
        }},
       {"v",
        [this](Variant const &value) {
-         set_particle_v(m_pid, get_value<Utils::Vector3d>(value));
+         set_particle_property(&Particle::v, value);
        },
-       [this]() { return particle().v(); }},
+       [this]() { return get_particle_data(m_pid).v(); }},
       {"f",
        [this](Variant const &value) {
-         set_particle_f(m_pid, get_value<Utils::Vector3d>(value));
+         set_particle_property(&Particle::force, value);
        },
-       [this]() { return particle().force(); }},
+       [this]() { return get_particle_data(m_pid).force(); }},
       {"mass",
 #ifdef MASS
        [this](Variant const &value) {
-         set_particle_mass(m_pid, get_value<double>(value));
+         set_particle_property(&Particle::mass, value);
        },
 #else  // MASS
        [](Variant const &value) {
-         if (std::abs(get_value<double>(value) - 1.) > 1e-10) {
+         auto const default_mass = Particle().mass();
+         if (std::abs(get_value<double>(value) - default_mass) > 1e-10) {
            throw std::runtime_error("Feature MASS not compiled in");
          }
        },
 #endif // MASS
-       [this]() { return particle().mass(); }},
+       [this]() { return get_particle_data(m_pid).mass(); }},
       {"q",
 #ifdef ELECTROSTATICS
        [this](Variant const &value) {
-         set_particle_q(m_pid, get_value<double>(value));
+         set_particle_property(&Particle::q, value);
        },
 #else  // ELECTROSTATICS
        [](Variant const &value) {
@@ -144,11 +219,11 @@ ParticleHandle::ParticleHandle() {
          }
        },
 #endif // ELECTROSTATICS
-       [this]() { return particle().q(); }},
+       [this]() { return get_particle_data(m_pid).q(); }},
       {"virtual",
 #ifdef VIRTUAL_SITES
        [this](Variant const &value) {
-         set_particle_virtual(m_pid, get_value<bool>(value));
+         set_particle_property(&Particle::virtual_flag, value);
        },
 #else  // VIRTUAL_SITES
        [](Variant const &value) {
@@ -157,129 +232,158 @@ ParticleHandle::ParticleHandle() {
          }
        },
 #endif // VIRTUAL_SITES
-       [this]() { return particle().is_virtual(); }},
+       [this]() { return get_particle_data(m_pid).is_virtual(); }},
 #ifdef ROTATION
       {"director",
        [this](Variant const &value) {
-         set_particle_director(m_pid, get_value<Utils::Vector3d>(value));
+         set_particle_property([&value](Particle &p) {
+           auto const director = get_value<Utils::Vector3d>(value).normalized();
+           p.quat() = Utils::convert_director_to_quaternion(director);
+         });
        },
-       [this]() { return particle().calc_director(); }},
+       [this]() {
+         auto const quat = get_particle_data(m_pid).quat();
+         return Utils::convert_quaternion_to_director(quat);
+       }},
       {"quat",
        [this](Variant const &value) {
-         set_particle_quat(m_pid, get_quaternion_safe("quat", value));
+         auto const quat = get_quaternion_safe("quat", value);
+         set_particle_property([&quat](Particle &p) { p.quat() = quat; });
        },
-       [this]() { return quat2vector(particle().quat()); }},
+       [this]() { return quat2vector(get_particle_data(m_pid).quat()); }},
       {"omega_body",
        [this](Variant const &value) {
-         set_particle_omega_body(m_pid, get_value<Utils::Vector3d>(value));
+         set_particle_property(&Particle::omega, value);
        },
-       [this]() { return particle().omega(); }},
+       [this]() { return get_particle_data(m_pid).omega(); }},
       {"rotation",
        [this](Variant const &value) {
-         set_particle_rotation(
-             m_pid, Utils::Vector3i{get_value<Utils::Vector3b>(value)});
+         set_particle_property([&value](Particle &p) {
+           auto const rotation_flag =
+               Utils::Vector3i{(get_value<Utils::Vector3b>(value))};
+           p.rotation() = bitfield_from_flag(rotation_flag);
+         });
        },
        [this]() {
-         auto const &p = particle();
-         return Utils::Vector3b{{p.can_rotate_around(0), p.can_rotate_around(1),
-                                 p.can_rotate_around(2)}};
+         auto const rotation_bits = get_particle_data(m_pid).rotation();
+         return Utils::Vector3b{{::detail::get_nth_bit(rotation_bits, 0),
+                                 ::detail::get_nth_bit(rotation_bits, 1),
+                                 ::detail::get_nth_bit(rotation_bits, 2)}};
        }},
       {"omega_lab",
        [this](Variant const &value) {
-         set_particle_omega_lab(m_pid, get_value<Utils::Vector3d>(value));
+         set_particle_property([&value](Particle &p) {
+           auto const omega = get_value<Utils::Vector3d>(value);
+           p.omega() = convert_vector_space_to_body(p, omega);
+         });
        },
        [this]() {
-         auto const &p = particle();
+         auto &p = get_particle_data(m_pid);
          return convert_vector_body_to_space(p, p.omega());
        }},
       {"torque_lab",
        [this](Variant const &value) {
-         set_particle_torque_lab(m_pid, get_value<Utils::Vector3d>(value));
+         set_particle_property([&value](Particle &p) {
+           auto const torque = get_value<Utils::Vector3d>(value);
+           p.torque() = convert_vector_space_to_body(p, torque);
+         });
        },
        [this]() {
-         auto const &p = particle();
+         auto &p = get_particle_data(m_pid);
          return convert_vector_body_to_space(p, p.torque());
        }},
 #endif // ROTATION
 #ifdef DIPOLES
       {"dip",
        [this](Variant const &value) {
-         set_particle_dip(m_pid, get_value<Utils::Vector3d>(value));
+         set_particle_property([&value](Particle &p) {
+           auto const dip = get_value<Utils::Vector3d>(value);
+           std::tie(p.quat(), p.dipm()) = convert_dip_to_quat(dip);
+         });
        },
-       [this]() { return particle().calc_dip(); }},
+       [this]() { return get_particle_data(m_pid).calc_dip(); }},
       {"dipm",
        [this](Variant const &value) {
-         set_particle_dipm(m_pid, get_value<double>(value));
+         set_particle_property(&Particle::dipm, value);
        },
-       [this]() { return particle().dipm(); }},
+       [this]() { return get_particle_data(m_pid).dipm(); }},
 #endif // DIPOLES
 #ifdef ROTATIONAL_INERTIA
       {"rinertia",
        [this](Variant const &value) {
-         set_particle_rotational_inertia(m_pid,
-                                         get_value<Utils::Vector3d>(value));
+         set_particle_property(&Particle::rinertia, value);
        },
-       [this]() { return particle().rinertia(); }},
+       [this]() { return get_particle_data(m_pid).rinertia(); }},
 #endif // ROTATIONAL_INERTIA
 #ifdef LB_ELECTROHYDRODYNAMICS
       {"mu_E",
        [this](Variant const &value) {
-         set_particle_mu_E(m_pid, get_value<Utils::Vector3d>(value));
+         set_particle_property(&Particle::mu_E, value);
        },
-       [this]() { return particle().mu_E(); }},
+       [this]() { return get_particle_data(m_pid).mu_E(); }},
 #endif // LB_ELECTROHYDRODYNAMICS
 #ifdef EXTERNAL_FORCES
       {"fix",
        [this](Variant const &value) {
-         set_particle_fix(m_pid,
-                          Utils::Vector3i{get_value<Utils::Vector3b>(value)});
+         set_particle_property([&value](Particle &p) {
+           auto const fix_flag =
+               Utils::Vector3i{(get_value<Utils::Vector3b>(value))};
+           p.fixed() = bitfield_from_flag(fix_flag);
+         });
        },
        [this]() {
-         auto const &p = particle();
-         return Utils::Vector3b{
-             {p.is_fixed_along(0), p.is_fixed_along(1), p.is_fixed_along(2)}};
+         auto const fixed = get_particle_data(m_pid).fixed();
+         return Utils::Vector3b{{::detail::get_nth_bit(fixed, 0),
+                                 ::detail::get_nth_bit(fixed, 1),
+                                 ::detail::get_nth_bit(fixed, 2)}};
        }},
       {"ext_force",
        [this](Variant const &value) {
-         set_particle_ext_force(m_pid, get_value<Utils::Vector3d>(value));
+         set_particle_property(&Particle::ext_force, value);
        },
-       [this]() { return particle().ext_force(); }},
+       [this]() { return get_particle_data(m_pid).ext_force(); }},
 #ifdef ROTATION
       {"ext_torque",
        [this](Variant const &value) {
-         set_particle_ext_torque(m_pid, get_value<Utils::Vector3d>(value));
+         set_particle_property(&Particle::ext_torque, value);
        },
-       [this]() { return particle().ext_torque(); }},
+       [this]() { return get_particle_data(m_pid).ext_torque(); }},
 #endif // ROTATION
 #endif // EXTERNAL_FORCES
 #ifdef THERMOSTAT_PER_PARTICLE
       {"gamma",
        [this](Variant const &value) {
-         set_particle_gamma(m_pid, get_gamma_safe(value));
+         set_particle_property(&Particle::gamma,
+                               Variant{get_gamma_safe(value)});
        },
-       [this]() { return particle().gamma(); }},
+       [this]() { return get_particle_data(m_pid).gamma(); }},
 #ifdef ROTATION
       {"gamma_rot",
        [this](Variant const &value) {
-         set_particle_gamma_rot(m_pid, get_gamma_safe(value));
+         set_particle_property(&Particle::gamma_rot,
+                               Variant{get_gamma_safe(value)});
        },
-       [this]() { return particle().gamma_rot(); }},
+       [this]() { return get_particle_data(m_pid).gamma_rot(); }},
 #endif // ROTATION
 #endif // THERMOSTAT_PER_PARTICLE
       {"pos_folded", AutoParameter::read_only,
-       [this]() { return folded_position(particle().pos(), ::box_geo); }},
+       [this]() {
+         return folded_position(get_particle_data(m_pid).pos(), ::box_geo);
+       }},
 
       {"lees_edwards_offset",
        [this](Variant const &value) {
-         set_particle_lees_edwards_offset(m_pid, get_value<double>(value));
+         set_particle_property(&Particle::lees_edwards_offset, value);
        },
-       [this]() { return particle().lees_edwards_offset(); }},
+       [this]() { return get_particle_data(m_pid).lees_edwards_offset(); }},
       {"lees_edwards_flag", AutoParameter::read_only,
-       [this]() { return particle().lees_edwards_flag(); }},
+       [this]() { return get_particle_data(m_pid).lees_edwards_flag(); }},
       {"image_box", AutoParameter::read_only,
-       [this]() { return particle().image_box(); }},
+       [this]() { return get_particle_data(m_pid).image_box(); }},
       {"node", AutoParameter::read_only,
-       [this]() { return get_particle_node(m_pid); }},
+       [this]() {
+         return (context()->is_head_node()) ? get_particle_node(m_pid) : -1;
+       }},
       {"mol_id",
        [this](Variant const &value) {
          auto const mol_id = get_value<int>(value);
@@ -287,15 +391,19 @@ ParticleHandle::ParticleHandle() {
            throw std::domain_error(
                error_msg("mol_id", "must be an integer >= 0"));
          }
-         set_particle_mol_id(m_pid, mol_id);
+         set_particle_property(&Particle::mol_id, Variant{mol_id});
        },
-       [this]() { return particle().mol_id(); }},
+       [this]() { return get_particle_data(m_pid).mol_id(); }},
 #ifdef VIRTUAL_SITES_RELATIVE
       {"vs_quat",
        [this](Variant const &value) {
-         set_particle_vs_quat(m_pid, get_quaternion_safe("vs_quat", value));
+         auto const quat = get_quaternion_safe("vs_quat", value);
+         set_particle_property(
+             [&quat](Particle &p) { p.vs_relative().quat = quat; });
        },
-       [this]() { return quat2vector(particle().vs_relative().quat); }},
+       [this]() {
+         return quat2vector(get_particle_data(m_pid).vs_relative().quat);
+       }},
       {"vs_relative",
        [this](Variant const &value) {
          ParticleProperties::VirtualSitesRelativeParameters vs_relative{};
@@ -312,56 +420,55 @@ ParticleHandle::ParticleHandle() {
            throw std::invalid_argument(error_msg(
                "vs_relative", "must take the form [id, distance, quaternion]"));
          }
-         set_particle_vs_relative(m_pid, vs_relative.to_particle_id,
-                                  vs_relative.distance,
-                                  vs_relative.rel_orientation);
+         set_particle_property(
+             [&vs_relative](Particle &p) { p.vs_relative() = vs_relative; });
        },
        [this]() {
-         auto const &p = particle();
-         return std::vector<Variant>{
-             {p.vs_relative().to_particle_id, p.vs_relative().distance,
-              quat2vector(p.vs_relative().rel_orientation)}};
+         auto const vs_rel = get_particle_data(m_pid).vs_relative();
+         return std::vector<Variant>{{vs_rel.to_particle_id, vs_rel.distance,
+                                      quat2vector(vs_rel.rel_orientation)}};
        }},
 #endif // VIRTUAL_SITES_RELATIVE
 #ifdef ENGINE
       {"swimming",
        [this](Variant const &value) {
-         ParticleParametersSwimming swim{};
-         swim.swimming = true;
-         auto const dict = get_value<VariantMap>(value);
-         if (dict.count("f_swim") != 0) {
-           swim.f_swim = get_value<double>(dict.at("f_swim"));
-         }
-         if (dict.count("v_swim") != 0) {
-           swim.v_swim = get_value<double>(dict.at("v_swim"));
-         }
-         if (swim.f_swim != 0. and swim.v_swim != 0.) {
-           throw std::invalid_argument(error_msg(
-               "swimming",
-               "cannot be set with 'v_swim' and 'f_swim' at the same time"));
-         }
-         if (dict.count("mode") != 0) {
-           auto const mode = get_value<std::string>(dict.at("mode"));
-           if (mode == "pusher") {
-             swim.push_pull = -1;
-           } else if (mode == "puller") {
-             swim.push_pull = +1;
-           } else if (mode == "N/A") {
-             swim.push_pull = 0;
-           } else {
+         set_particle_property([&value](Particle &p) {
+           ParticleParametersSwimming swim{};
+           swim.swimming = true;
+           auto const dict = get_value<VariantMap>(value);
+           if (dict.count("f_swim") != 0) {
+             swim.f_swim = get_value<double>(dict.at("f_swim"));
+           }
+           if (dict.count("v_swim") != 0) {
+             swim.v_swim = get_value<double>(dict.at("v_swim"));
+           }
+           if (swim.f_swim != 0. and swim.v_swim != 0.) {
              throw std::invalid_argument(
-                 error_msg("swimming.mode",
-                           "has to be either 'pusher', 'puller' or 'N/A'"));
+                 error_msg("swimming", "cannot be set with 'v_swim' and "
+                                       "'f_swim' at the same time"));
            }
-         }
-         if (dict.count("dipole_length") != 0) {
-           swim.dipole_length = get_value<double>(dict.at("dipole_length"));
-         }
-         set_particle_swimming(m_pid, swim);
+           if (dict.count("mode") != 0) {
+             auto const mode = get_value<std::string>(dict.at("mode"));
+             if (mode == "pusher") {
+               swim.push_pull = -1;
+             } else if (mode == "puller") {
+               swim.push_pull = +1;
+             } else if (mode == "N/A") {
+               swim.push_pull = 0;
+             } else {
+               throw std::invalid_argument(
+                   error_msg("swimming.mode",
+                             "has to be either 'pusher', 'puller' or 'N/A'"));
+             }
+           }
+           if (dict.count("dipole_length") != 0) {
+             swim.dipole_length = get_value<double>(dict.at("dipole_length"));
+           }
+           p.swimming() = swim;
+         });
        },
        [this]() {
-         auto const &p = particle();
-         auto const &swim = p.swimming();
+         auto const swim = get_particle_data(m_pid).swimming();
          std::string mode;
          if (swim.push_pull == -1) {
            mode = "pusher";
@@ -379,11 +486,65 @@ ParticleHandle::ParticleHandle() {
   });
 }
 
+#ifdef EXCLUSIONS
+/**
+ * @brief Locally add an exclusion to a particle.
+ * @param pid1 the identity of the first exclusion partner
+ * @param pid2 the identity of the second exclusion partner
+ */
+static void local_add_exclusion(int pid1, int pid2) {
+  if (auto p1 = ::cell_structure.get_local_particle(pid1)) {
+    add_exclusion(*p1, pid2);
+  }
+  if (auto p2 = ::cell_structure.get_local_particle(pid2)) {
+    add_exclusion(*p2, pid1);
+  }
+}
+
+/**
+ * @brief Locally remove an exclusion to a particle.
+ * @param pid1 the identity of the first exclusion partner
+ * @param pid2 the identity of the second exclusion partner
+ */
+static void local_remove_exclusion(int pid1, int pid2) {
+  if (auto p1 = ::cell_structure.get_local_particle(pid1)) {
+    delete_exclusion(*p1, pid2);
+  }
+  if (auto p2 = ::cell_structure.get_local_particle(pid2)) {
+    delete_exclusion(*p2, pid1);
+  }
+}
+
+void ParticleHandle::particle_exclusion_sanity_checks(int pid1,
+                                                      int pid2) const {
+  if (pid1 == pid2) {
+    throw std::runtime_error("Particles cannot exclude themselves (id " +
+                             std::to_string(pid1) + ")");
+  }
+  std::ignore = get_real_particle(context()->get_comm(), pid1);
+  std::ignore = get_real_particle(context()->get_comm(), pid2);
+}
+#endif // EXCLUSIONS
+
 Variant ParticleHandle::do_call_method(std::string const &name,
                                        VariantMap const &params) {
+  if (name == "set_param_parallel") {
+    auto const param_name = get_value<std::string>(params, "name");
+    if (params.count("value") == 0) {
+      throw Exception("Parameter '" + param_name + "' is missing.");
+    }
+    auto const &value = params.at("value");
+    context()->parallel_try_catch(
+        [&]() { do_set_parameter(param_name, value); });
+    return {};
+  }
   if (name == "get_bonds_view") {
+    if (not context()->is_head_node()) {
+      return {};
+    }
+    auto const bond_list = get_particle_data(m_pid).bonds();
     std::vector<std::vector<int>> bonds_flat;
-    for (auto const &bond_view : get_particle_bonds(m_pid)) {
+    for (auto const &&bond_view : bond_list) {
       std::vector<int> bond_flat;
       bond_flat.emplace_back(bond_view.bond_id());
       for (auto const pid : bond_view.partner_ids()) {
@@ -394,35 +555,87 @@ Variant ParticleHandle::do_call_method(std::string const &name,
     return make_vector_of_variants(bonds_flat);
   }
   if (name == "add_bond") {
-    add_particle_bond(m_pid, get_bond_vector(params));
+    set_particle_property([&params](Particle &p) {
+      auto const bond_id = get_value<int>(params, "bond_id");
+      auto const part_id = get_value<std::vector<int>>(params, "part_id");
+      auto const bond_view =
+          BondView(bond_id, {part_id.data(), part_id.size()});
+      p.bonds().insert(bond_view);
+    });
   } else if (name == "del_bond") {
-    delete_particle_bond(m_pid, get_bond_vector(params));
+    set_particle_property([&params](Particle &p) {
+      auto const bond_id = get_value<int>(params, "bond_id");
+      auto const part_id = get_value<std::vector<int>>(params, "part_id");
+      auto const bond_view =
+          BondView(bond_id, {part_id.data(), part_id.size()});
+      auto &bond_list = p.bonds();
+      auto it = std::find(bond_list.begin(), bond_list.end(), bond_view);
+      if (it != bond_list.end()) {
+        bond_list.erase(it);
+      }
+    });
   } else if (name == "delete_all_bonds") {
-    delete_particle_bonds(m_pid);
+    set_particle_property([&](Particle &p) { p.bonds().clear(); });
   } else if (name == "is_valid_bond_id") {
     auto const bond_id = get_value<int>(params, "bond_id");
     return ::bonded_ia_params.get_zero_based_type(bond_id) != 0;
   }
   if (name == "remove_particle") {
-    remove_particle(m_pid);
+    context()->parallel_try_catch([&]() {
+      std::ignore = get_real_particle(context()->get_comm(), m_pid);
+      remove_particle_parallel(m_pid);
+    });
 #ifdef VIRTUAL_SITES_RELATIVE
   } else if (name == "vs_relate_to") {
-    vs_relate_to(m_pid, get_value<int>(params, "pid"));
+    if (not context()->is_head_node()) {
+      return {};
+    }
+    auto const other_pid = get_value<int>(params, "pid");
+    if (m_pid == other_pid) {
+      throw std::invalid_argument("A virtual site cannot relate to itself");
+    }
+    if (other_pid < 0) {
+      throw std::domain_error("Invalid particle id: " +
+                              std::to_string(other_pid));
+    }
+    /* note: this code can be rewritten as parallel code, but only with a call
+     * to `cells_update_ghosts(DATA_PART_POSITION | DATA_PART_PROPERTIES)`, as
+     * there is no guarantee the virtual site has visibility of the relative
+     * particle through the ghost layer during particle creation. However,
+     * ghost updates can scramble the particle ordering in the local cells,
+     * which is an issue for checkpointing: the H5MD writer will use the
+     * scrambled ordering before writing to a checkpoint file and the
+     * non-scrambled ordering after reloading from a checkpoint file.
+     */
+    auto const &p_current = get_particle_data(m_pid);
+    auto const &p_relate_to = get_particle_data(other_pid);
+    auto const [quat, dist] =
+        calculate_vs_relate_to_params(p_current, p_relate_to);
+    set_parameter("vs_relative", Variant{std::vector<Variant>{
+                                     {other_pid, dist, quat2vector(quat)}}});
+    set_parameter("virtual", true);
 #endif // VIRTUAL_SITES_RELATIVE
 #ifdef EXCLUSIONS
   } else if (name == "has_exclusion") {
-    auto const &p = get_particle_data(m_pid);
-    return p.has_exclusion(get_value<int>(params, "pid"));
+    auto const other_pid = get_value<int>(params, "pid");
+    auto const p = get_real_particle(context()->get_comm(), m_pid);
+    if (p != nullptr) {
+      return p->has_exclusion(other_pid);
+    }
   }
   if (name == "add_exclusion") {
-    add_particle_exclusion(m_pid, get_value<int>(params, "pid"));
+    auto const other_pid = get_value<int>(params, "pid");
+    context()->parallel_try_catch(
+        [&]() { particle_exclusion_sanity_checks(m_pid, other_pid); });
+    local_add_exclusion(m_pid, other_pid);
+    on_particle_change();
   } else if (name == "del_exclusion") {
-    remove_particle_exclusion(m_pid, get_value<int>(params, "pid"));
+    auto const other_pid = get_value<int>(params, "pid");
+    context()->parallel_try_catch(
+        [&]() { particle_exclusion_sanity_checks(m_pid, other_pid); });
+    local_remove_exclusion(m_pid, other_pid);
+    on_particle_change();
   } else if (name == "set_exclusions") {
-    auto const &p = particle();
-    for (auto const pid : p.exclusions_as_vector()) {
-      remove_particle_exclusion(m_pid, pid);
-    }
     std::vector<int> exclusion_list;
     try {
       auto const pid = get_value<int>(params, "p_ids");
@@ -430,31 +643,50 @@ Variant ParticleHandle::do_call_method(std::string const &name,
     } catch (...) {
       exclusion_list = get_value<std::vector<int>>(params, "p_ids");
     }
-    for (auto const pid : exclusion_list) {
-      if (!p.has_exclusion(pid)) {
-        add_particle_exclusion(m_pid, pid);
+    context()->parallel_try_catch([&]() {
+      for (auto const pid : exclusion_list) {
+        particle_exclusion_sanity_checks(m_pid, pid);
       }
-    }
+    });
+    set_particle_property([this, &exclusion_list](Particle &p) {
+      for (auto const pid : p.exclusions()) {
+        local_remove_exclusion(m_pid, pid);
+      }
+      for (auto const pid : exclusion_list) {
+        if (!p.has_exclusion(pid)) {
+          local_add_exclusion(m_pid, pid);
+        }
+      }
+    });
   } else if (name == "get_exclusions") {
-    return particle().exclusions_as_vector();
+    if (not context()->is_head_node()) {
+      return {};
+    }
+    auto const excl_list = get_particle_data(m_pid).exclusions();
+    return Variant{std::vector<int>{excl_list.begin(), excl_list.end()}};
 #endif // EXCLUSIONS
 #ifdef ROTATION
   }
   if (name == "rotate_particle") {
-    rotate_particle(m_pid, get_value<Utils::Vector3d>(params, "axis"),
-                    get_value<double>(params, "angle"));
+    set_particle_property([&params](Particle &p) {
+      auto const axis = get_value<Utils::Vector3d>(params, "axis");
+      auto const angle = get_value<double>(params, "angle");
+      local_rotate_particle(p, axis, angle);
+    });
   }
   if (name == "convert_vector_body_to_space") {
-    auto const &p = get_particle_data(m_pid);
-    return convert_vector_body_to_space(
-               p, get_value<Utils::Vector3d>(params, "vec"))
-        .as_vector();
+    return get_particle_property<std::vector<double>>(
+        [&params](Particle const &p) {
+          auto const vec = get_value<Utils::Vector3d>(params, "vec");
+          return convert_vector_body_to_space(p, vec).as_vector();
+        });
   }
   if (name == "convert_vector_space_to_body") {
-    auto const &p = get_particle_data(m_pid);
-    return convert_vector_space_to_body(
-               p, get_value<Utils::Vector3d>(params, "vec"))
-        .as_vector();
+    return get_particle_property<std::vector<double>>(
+        [&params](Particle const &p) {
+          auto const vec = get_value<Utils::Vector3d>(params, "vec");
+          return convert_vector_space_to_body(p, vec).as_vector();
+        });
 #endif // ROTATION
   }
   return {};
@@ -481,15 +713,33 @@ void ParticleHandle::do_construct(VariantMap const &params) {
     return params.count(key) == 1;
   };
   m_pid = (has_param("id")) ? get_value<int>(params, "id")
-                            : get_maximal_particle_id() + 1;
+                            : get_maximal_particle_id_parallel() + 1;
+
+#ifndef NDEBUG
+  if (!has_param("id")) {
+    auto head_node_reference = m_pid;
+    boost::mpi::broadcast(context()->get_comm(), head_node_reference, 0);
+    assert(m_pid == head_node_reference && "global max_seen_pid has diverged");
+  }
+#endif
 
   // create a new particle if extra arguments were passed
-  if (n_extra_args > 0) {
-    if (particle_exists(m_pid)) {
+  if (n_extra_args == 0) {
+    return;
+  }
+
+  auto const pos = get_value<Utils::Vector3d>(params, "pos");
+  context()->parallel_try_catch([&]() {
+    particle_checks(m_pid, pos);
+    auto ptr = ::cell_structure.get_local_particle(m_pid);
+    if (ptr != nullptr) {
       throw std::invalid_argument("Particle " + std::to_string(m_pid) +
                                   " already exists");
     }
+  });
+
 #ifdef ROTATION
+  context()->parallel_try_catch([&]() {
     // if we are not constructing a particle from a checkpoint file,
     // check the quaternion is not accidentally set twice by the user
     if (not has_param("__cpt_sentinel")) {
@@ -502,14 +752,15 @@ void ParticleHandle::do_construct(VariantMap const &params) {
         }
       }
     }
+  });
 #endif // ROTATION
 
-    // create a default-constructed particle
-    auto const pos = get_value<Utils::Vector3d>(params, "pos");
-    mpi_make_new_particle(m_pid, pos);
+  // create a default-constructed particle
+  mpi_make_new_particle_local(m_pid, pos);
 
+  context()->parallel_try_catch([&]() {
     // set particle properties (filter out read-only and deferred properties)
-    std::vector<std::string> skip = {
+    std::set<std::string> const skip = {
         "pos_folded", "pos", "quat", "director",  "id",    "lees_edwards_flag",
         "exclusions", "dip", "node", "image_box", "bonds", "__cpt_sentinel",
     };
@@ -518,23 +769,27 @@ void ParticleHandle::do_construct(VariantMap const &params) {
     // can be allowed to; these conditionals are required to handle a reload
     // from a checkpoint file, where all properties exist (avoids accidentally
     // overwriting the quaternion by the default-constructed dipole moment)
-    if (has_param("quat")) {
-      do_set_parameter("quat", params.at("quat"));
-    } else if (has_param("director")) {
-      do_set_parameter("director", params.at("director"));
-    } else if (has_param("dip")) {
-      do_set_parameter("dip", params.at("dip"));
+    for (std::string name : {"quat", "director", "dip"}) {
+      if (has_param(name)) {
+        do_set_parameter(name, params.at(name));
+        break;
+      }
     }
 #endif // ROTATION
-    for (auto const &kv : params) {
-      if (std::find(skip.begin(), skip.end(), kv.first) == skip.end()) {
-        do_set_parameter(kv.first, kv.second);
+    std::vector<std::string> sorted_param_names = {};
+    std::for_each(params.begin(), params.end(), [&](auto const &kv) {
+      if (skip.count(kv.first) == 0) {
+        sorted_param_names.push_back(kv.first);
       }
+    });
+    std::sort(sorted_param_names.begin(), sorted_param_names.end());
+    for (auto const &name : sorted_param_names) {
+      do_set_parameter(name, params.at(name));
     }
     if (not has_param("type")) {
       do_set_parameter("type", 0);
     }
-  }
+  });
 }
 
 } // namespace Particles
diff --git a/src/script_interface/particle_data/ParticleHandle.hpp b/src/script_interface/particle_data/ParticleHandle.hpp
index 833a9dd9490..eee2c9f69ea 100644
--- a/src/script_interface/particle_data/ParticleHandle.hpp
+++ b/src/script_interface/particle_data/ParticleHandle.hpp
@@ -27,15 +27,24 @@
 
 #include <string>
 
-Particle const &get_particle_data(int p_id);
-
 namespace ScriptInterface {
 namespace Particles {
 
 class ParticleHandle : public AutoParameters<ParticleHandle> {
   int m_pid;
 
-  Particle const &particle() const { return get_particle_data(m_pid); }
+  template <typename T>
+  T get_particle_property(T const &(Particle::*getter)() const) const;
+  template <typename T, class F> T get_particle_property(F const &fun) const;
+
+  template <typename T>
+  void set_particle_property(T &(Particle::*setter)(),
+                             Variant const &value) const;
+
+  template <class F> void set_particle_property(F const &fun) const;
+#ifdef EXCLUSIONS
+  void particle_exclusion_sanity_checks(int pid1, int pid2) const;
+#endif // EXCLUSIONS
 
 public:
   ParticleHandle();
diff --git a/src/script_interface/particle_data/ParticleList.cpp b/src/script_interface/particle_data/ParticleList.cpp
index 8304ea4b7fe..9552ec16a9f 100644
--- a/src/script_interface/particle_data/ParticleList.cpp
+++ b/src/script_interface/particle_data/ParticleList.cpp
@@ -23,11 +23,19 @@
 #include "script_interface/ObjectState.hpp"
 #include "script_interface/ScriptInterface.hpp"
 
+#include "core/cells.hpp"
+#include "core/event.hpp"
+#include "core/exclusions.hpp"
 #include "core/particle_node.hpp"
 
 #include <utils/Vector.hpp>
+#include <utils/mpi/gather_buffer.hpp>
 #include <utils/serialization/pack.hpp>
 
+#include <boost/mpi/collectives.hpp>
+#include <boost/mpi/communicator.hpp>
+#include <boost/range/algorithm.hpp>
+
 #include <memory>
 #include <stdexcept>
 #include <string>
@@ -40,7 +48,7 @@ namespace Particles {
 
 #ifdef EXCLUSIONS
 static void set_exclusions(ParticleHandle &p, Variant const &exclusions) {
-  p.do_call_method("set_exclusions", {{"p_ids", exclusions}});
+  p.call_method("set_exclusions", {{"p_ids", exclusions}});
 }
 #endif // EXCLUSIONS
 
@@ -50,8 +58,8 @@ static void set_bonds(ParticleHandle &p, Variant const &bonds) {
     auto const bond_id = bond_flat[0];
     auto const part_id =
         std::vector<int>{bond_flat.begin() + 1, bond_flat.end()};
-    p.do_call_method("add_bond",
-                     {{"bond_id", bond_id}, {"part_id", std::move(part_id)}});
+    p.call_method("add_bond",
+                  {{"bond_id", bond_id}, {"part_id", std::move(part_id)}});
   }
 }
 
@@ -59,9 +67,10 @@ std::string ParticleList::get_internal_state() const {
   auto const p_ids = get_particle_ids();
   std::vector<std::string> object_states(p_ids.size());
 
-  boost::transform(p_ids, object_states.begin(), [](auto const p_id) {
-    ParticleHandle p_handle{};
-    p_handle.do_construct({{"id", p_id}});
+  boost::transform(p_ids, object_states.begin(), [this](auto const p_id) {
+    auto p_obj =
+        context()->make_shared("Particles::ParticleHandle", {{"id", p_id}});
+    auto &p_handle = dynamic_cast<ParticleHandle &>(*p_obj);
     auto const packed_state = p_handle.serialize();
     // custom particle serialization
     auto state = Utils::unpack<ObjectState>(packed_state);
@@ -90,25 +99,23 @@ void ParticleList::set_internal_state(std::string const &state) {
   std::unordered_map<int, Variant> bonds = {};
 
   for (auto const &packed_object : object_states) {
-    auto const state = Utils::unpack<ObjectState>(packed_object);
-    auto o = std::dynamic_pointer_cast<ParticleHandle>(
-        ObjectHandle::deserialize(packed_object, *ObjectHandle::context()));
-    auto const p_id = get_value<int>(o->get_parameter("id"));
+    auto state = Utils::unpack<ObjectState>(packed_object);
+    VariantMap params = {};
     for (auto const &kv : state.params) {
-      if (kv.first == "bonds") {
-        bonds[p_id] = unpack(kv.second, {});
-      }
+      params[kv.first] = unpack(kv.second, {});
+    }
+    auto const p_id = get_value<int>(params.at("id"));
+    bonds[p_id] = params.extract("bonds").mapped();
 #ifdef EXCLUSIONS
-      else if (kv.first == "exclusions") {
-        exclusions[p_id] = unpack(kv.second, {});
-      }
+    exclusions[p_id] = params.extract("exclusions").mapped();
 #endif // EXCLUSIONS
-    }
+    context()->make_shared("Particles::ParticleHandle", params);
   }
 
   for (auto const p_id : get_particle_ids()) {
-    ParticleHandle p_handle{};
-    p_handle.do_construct({{"id", p_id}});
+    auto p_obj =
+        context()->make_shared("Particles::ParticleHandle", {{"id", p_id}});
+    auto &p_handle = dynamic_cast<ParticleHandle &>(*p_obj);
     set_bonds(p_handle, bonds[p_id]);
 #ifdef EXCLUSIONS
     set_exclusions(p_handle, exclusions[p_id]);
@@ -116,8 +123,106 @@ void ParticleList::set_internal_state(std::string const &state) {
   }
 }
 
+#ifdef EXCLUSIONS
+/**
+ * @brief Use the bond topology to automatically add exclusions between
+ * particles that are up to @c n_bonds_max bonds apart in a chain.
+ */
+static void auto_exclusions(boost::mpi::communicator const &comm,
+                            int const n_bonds_max) {
+  // bookkeeping of particle exclusions, with their n-bond distance
+  std::unordered_map<int, std::vector<std::pair<int, int>>> partners;
+  std::vector<int> bonded_pairs;
+
+  // determine initial connectivity
+  for (auto const &p : ::cell_structure.local_particles()) {
+    auto const pid1 = p.id();
+    for (auto const bond : p.bonds()) {
+      if (bond.partner_ids().size() == 1) {
+        auto const pid2 = bond.partner_ids()[0];
+        if (pid1 != pid2) {
+          bonded_pairs.emplace_back(pid1);
+          bonded_pairs.emplace_back(pid2);
+        }
+      }
+    }
+  }
+
+  Utils::Mpi::gather_buffer(bonded_pairs, comm);
+
+  if (comm.rank() == 0) {
+    auto const add_partner = [&partners](int pid1, int pid2, int n_bonds) {
+      if (pid2 == pid1)
+        return;
+      for (auto const &partner : partners[pid1])
+        if (partner.first == pid2)
+          return;
+      partners[pid1].emplace_back(std::pair{pid2, n_bonds});
+    };
+
+    for (auto it = bonded_pairs.begin(); it != bonded_pairs.end(); it += 2) {
+      add_partner(it[0], it[1], 1);
+      add_partner(it[1], it[0], 1);
+    }
+
+    // determine transient connectivity
+    for (int iteration = 1; iteration < n_bonds_max; iteration++) {
+      std::vector<int> pids;
+      for (auto const &kv : partners) {
+        pids.emplace_back(kv.first);
+      }
+      for (auto const pid1 : pids) {
+        // loop over partners (counter-based loops due to iterator invalidation)
+        // NOLINTNEXTLINE(modernize-loop-convert)
+        for (int i = 0; i < partners[pid1].size(); ++i) {
+          auto const [pid2, dist21] = partners[pid1][i];
+          if (dist21 > n_bonds_max)
+            continue;
+          // loop over all partners of the partner
+          // NOLINTNEXTLINE(modernize-loop-convert)
+          for (int j = 0; j < partners[pid2].size(); ++j) {
+            auto const [pid3, dist32] = partners[pid2][j];
+            auto const dist31 = dist32 + dist21;
+            if (dist31 > n_bonds_max)
+              continue;
+            add_partner(pid1, pid3, dist31);
+            add_partner(pid3, pid1, dist31);
+          }
+        }
+      }
+    }
+  }
+
+  boost::mpi::broadcast(comm, partners, 0);
+  for (auto const &kv : partners) {
+    auto const pid1 = kv.first;
+    auto const &partner_list = kv.second;
+    for (auto const &partner : partner_list) {
+      auto const pid2 = partner.first;
+      if (auto p1 = ::cell_structure.get_local_particle(pid1)) {
+        add_exclusion(*p1, pid2);
+      }
+      if (auto p2 = ::cell_structure.get_local_particle(pid2)) {
+        add_exclusion(*p2, pid1);
+      }
+    }
+  }
+  on_particle_change();
+}
+#endif // EXCLUSIONS
+
 Variant ParticleList::do_call_method(std::string const &name,
                                      VariantMap const &params) {
+#ifdef EXCLUSIONS
+  if (name == "auto_exclusions") {
+    auto const distance = get_value<int>(params, "distance");
+    auto_exclusions(context()->get_comm(), distance);
+    return {};
+  }
+#endif // EXCLUSIONS
+  if (not context()->is_head_node()) {
+    return {};
+  }
   if (name == "get_n_part") {
     return get_n_part();
   }
@@ -131,21 +236,8 @@ Variant ParticleList::do_call_method(std::string const &name,
     remove_all_particles();
   } else if (name == "add_particle") {
     assert(params.count("bonds") == 0);
-    // sanitize particle properties
-#ifdef DIPOLES
-    if (params.count("dip") and params.count("dipm")) {
-      throw std::invalid_argument("Contradicting attributes: 'dip' and 'dipm'. \
-Setting 'dip' is sufficient as the length of the vector defines the scalar \
-dipole moment.");
-    }
-    if (params.count("dip") and params.count("quat")) {
-      throw std::invalid_argument("Contradicting attributes: 'dip' and 'quat'. \
-Setting 'dip' overwrites the rotation of the particle around the dipole axis. \
-Set attribute 'quat' together with 'dipm' (scalar dipole moment) instead.");
-    }
-#endif // DIPOLES
-    ParticleHandle p_handle{};
-    p_handle.do_construct(params);
+    auto obj = context()->make_shared("Particles::ParticleHandle", params);
+    auto &p_handle = dynamic_cast<ParticleHandle &>(*obj);
 #ifdef EXCLUSIONS
     if (params.count("exclusions")) {
       set_exclusions(p_handle, params.at("exclusions"));
diff --git a/src/script_interface/system/System.cpp b/src/script_interface/system/System.cpp
index 62f74339601..5b00ffb1b49 100644
--- a/src/script_interface/system/System.cpp
+++ b/src/script_interface/system/System.cpp
@@ -21,9 +21,10 @@
 
 #include "config/config.hpp"
 
+#include "core/cells.hpp"
+#include "core/event.hpp"
 #include "core/grid.hpp"
 #include "core/object-in-fluid/oif_global_forces.hpp"
-#include "core/particle_data.hpp"
 #include "core/particle_node.hpp"
 #include "core/rotate_system.hpp"
 
@@ -43,6 +44,21 @@ System::System() {
   });
 }
 
+/** Rescale all particle positions in direction @p dir by a factor @p scale.
+ *  @param dir   direction to scale (0/1/2 = x/y/z, 3 = x+y+z isotropically)
+ *  @param scale factor by which to rescale (>1: stretch, <1: contract)
+ */
+static void rescale_particles(int dir, double scale) {
+  for (auto &p : ::cell_structure.local_particles()) {
+    if (dir < 3)
+      p.pos()[dir] *= scale;
+    else {
+      p.pos() *= scale;
+    }
+  }
+  on_particle_change();
+}
+
 Variant System::do_call_method(std::string const &name,
                                VariantMap const &parameters) {
   if (name == "is_system_created") {
@@ -79,15 +95,6 @@ Variant System::do_call_method(std::string const &name,
     }
     return {};
   }
-#ifdef EXCLUSIONS
-  if (name == "auto_exclusions") {
-    if (context()->is_head_node()) {
-      auto const distance = get_value<int>(parameters, "distance");
-      auto_exclusions(distance);
-    }
-    return {};
-  }
-#endif
   if (name == "setup_type_map") {
     if (context()->is_head_node()) {
       auto const types = get_value<std::vector<int>>(parameters, "type_list");
diff --git a/src/utils/include/utils/Bag.hpp b/src/utils/include/utils/Bag.hpp
index 83c8553b904..86f47a22f2e 100644
--- a/src/utils/include/utils/Bag.hpp
+++ b/src/utils/include/utils/Bag.hpp
@@ -26,6 +26,7 @@
 #include <boost/serialization/vector.hpp>
 
 #include <cstddef>
+#include <type_traits>
 #include <utility>
 #include <vector>
 
@@ -40,12 +41,14 @@ namespace Utils {
  *
  * Elements in the container do not have a stable position and
  * removing elements can change the order of the other elements.
- * The loser contract (compared to a vector) allows removing any
+ * The looser contract (compared to a vector) allows removing any
  * element in the container in constant time.
  *
- * @tparam T Element Type, needs to be Swappable.
+ * @tparam T Element type, needs to be Swappable.
  */
 template <class T> class Bag {
+  static_assert(std::is_swappable_v<T>);
+
   /** Storage backend */
   using storage_type = std::vector<T>;
 
@@ -63,7 +66,7 @@ template <class T> class Bag {
 
 private:
   /** Underlying storage of the container */
-  std::vector<T> m_storage;
+  storage_type m_storage;
 
   friend boost::serialization::access;
   /**
@@ -110,7 +113,7 @@ template <class T> class Bag {
    *
    * Increase capacity to at least the specified value.
    *
-   *     @param new_capacity New minimum capacity.
+   * @param new_capacity New minimum capacity.
    */
   void reserve(std::size_t new_capacity) { m_storage.reserve(new_capacity); }
 
@@ -121,7 +124,7 @@ template <class T> class Bag {
    * If the new size is larger than the capacity, all
    * iterators into the container are invalidated.
    *
-   *     @param new_size Size to resize to.
+   * @param new_size Size to resize to.
    */
   void resize(std::size_t new_size) { m_storage.resize(new_size); }
 
diff --git a/testsuite/python/CMakeLists.txt b/testsuite/python/CMakeLists.txt
index 05c43ca0fbf..e878840e18b 100644
--- a/testsuite/python/CMakeLists.txt
+++ b/testsuite/python/CMakeLists.txt
@@ -172,7 +172,7 @@ python_test(FILE rotational_inertia.py MAX_NUM_PROC 4)
 python_test(FILE rotational-diffusion-aniso.py MAX_NUM_PROC 1 LABELS long)
 python_test(FILE rotational_dynamics.py MAX_NUM_PROC 1)
 python_test(FILE script_interface.py MAX_NUM_PROC 4)
-python_test(FILE reaction_methods_interface.py MAX_NUM_PROC 1)
+python_test(FILE reaction_methods_interface.py MAX_NUM_PROC 2)
 python_test(FILE reaction_ensemble.py MAX_NUM_PROC 4)
 python_test(FILE reaction_complex.py MAX_NUM_PROC 1)
 python_test(FILE reaction_bookkeeping.py MAX_NUM_PROC 1)
diff --git a/testsuite/python/auto_exclusions.py b/testsuite/python/auto_exclusions.py
index 37f5dd46866..2efbed05a42 100644
--- a/testsuite/python/auto_exclusions.py
+++ b/testsuite/python/auto_exclusions.py
@@ -112,7 +112,7 @@ def check():
             length = len(system.part)
             for pid in range(length):
                 system.part.by_id(pid).add_bond((bond, (pid + 1) % length))
-            system.auto_exclusions(distance=2)
+            system.part.auto_exclusions(distance=2)
 
             for pid in range(length):
                 excl = sorted(list(system.part.by_id(pid).exclusions))
@@ -145,7 +145,7 @@ def test_branched(self):
         p1.add_bond((bond, p0))
         p2.add_bond((bond, p1))
         p3.add_bond((bond, p1))
-        system.auto_exclusions(distance=2)
+        system.part.auto_exclusions(distance=2)
 
         self.assertEqual(sorted(list(p0.exclusions)), [1, 2, 3])
         self.assertEqual(sorted(list(p1.exclusions)), [0, 2, 3])
@@ -172,7 +172,7 @@ def test_diamond(self):
         p2.add_bond((bond, p4))
         p3.add_bond((bond, p4))
         p5.add_bond((bond, p4))
-        system.auto_exclusions(distance=3)
+        system.part.auto_exclusions(distance=3)
 
         self.assertEqual(sorted(list(p0.exclusions)), [1, 2, 3, 4])
         self.assertEqual(sorted(list(p1.exclusions)), [0, 2, 3, 4, 5])
@@ -191,7 +191,7 @@ def test_id_gaps(self):
 
         # topology: 0-4
         p0.add_bond((bond, p4))
-        system.auto_exclusions(distance=1)
+        system.part.auto_exclusions(distance=1)
 
         self.assertEqual(list(p0.exclusions), [4])
         self.assertEqual(list(p4.exclusions), [0])
@@ -219,7 +219,7 @@ def test_non_pairwise(self):
         if espressomd.has_features(["VIRTUAL_SITES_INERTIALESS_TRACERS"]):
             system.part.by_id(0).add_bond((volcons,))
 
-        system.auto_exclusions(distance=2)
+        system.part.auto_exclusions(distance=2)
 
         for p in system.part.all():
             self.assertEqual(list(p.exclusions), [])
diff --git a/testsuite/python/nsquare.py b/testsuite/python/nsquare.py
index ec57cb304ef..7351132dc14 100644
--- a/testsuite/python/nsquare.py
+++ b/testsuite/python/nsquare.py
@@ -44,8 +44,7 @@ def test_load_balancing(self):
             pos=n_part * [(0, 0, 0)], type=n_part * [1])
 
         # And now change their positions
-        partcls.pos = self.system.box_l * \
-            np.random.random((n_part, 3))
+        partcls.pos = self.system.box_l * np.random.random((n_part, 3))
 
         # Add an interacting particle in a corner of the box
         self.system.part.add(pos=[(0.01, 0.01, 0.01)], type=[0])
diff --git a/testsuite/python/particle.py b/testsuite/python/particle.py
index 097159c5997..b547f647823 100644
--- a/testsuite/python/particle.py
+++ b/testsuite/python/particle.py
@@ -446,8 +446,9 @@ def test_zz_remove_all(self):
         np.random.shuffle(ids)
         for pid in ids:
             self.system.part.by_id(pid).remove()
-        with self.assertRaises(Exception):
-            self.system.part.by_id(17).remove()
+        p = self.system.part.by_id(17)
+        with self.assertRaises(RuntimeError):
+            p.remove()
 
     def test_coord_fold_corner_cases(self):
         system = self.system
@@ -563,6 +564,8 @@ def test_update(self):
         new_pdict = p.to_dict()
         del new_pdict['id']
         self.assertEqual(str(new_pdict), str(pdict))
+        with self.assertRaisesRegex(RuntimeError, "Parameter 'test' is missing"):
+            p.call_method("set_param_parallel", name="test")
 
 
 if __name__ == "__main__":
diff --git a/testsuite/python/reaction_methods_interface.py b/testsuite/python/reaction_methods_interface.py
index d9072b734df..6e1f07be643 100644
--- a/testsuite/python/reaction_methods_interface.py
+++ b/testsuite/python/reaction_methods_interface.py
@@ -152,13 +152,15 @@ def count_by_type(types):
             reaction_backward['gamma'],
             delta=1e-10)
 
-        # check particle deletion
+        # check particle deletion on a worker node
         p1, _, p3 = self.system.part.add(
-            pos=3 * [(0., 0., 0.)], type=[5, 2, 3])
+            pos=3 * [(-1., -1., -1.)], type=[5, 2, 3])
         if isinstance(method, espressomd.reaction_methods.WidomInsertion):
             potential_energy = method.calculate_particle_insertion_potential_energy(
                 reaction_id=0)
             self.assertEqual(potential_energy, 0.)
+        if self.system.cell_system.get_state()["n_nodes"] > 1:
+            assert set(self.system.part.all().node) != {0}
         self.assertEqual(count_by_type([5, 2, 3, 0]), [1, 1, 1, 0])
         method.delete_particle(p_id=p3.id)
         self.assertEqual(count_by_type([5, 2, 3, 0]), [1, 1, 0, 0])
diff --git a/testsuite/python/regular_decomposition.py b/testsuite/python/regular_decomposition.py
index cf96f8caa08..b05746f5423 100644
--- a/testsuite/python/regular_decomposition.py
+++ b/testsuite/python/regular_decomposition.py
@@ -44,8 +44,10 @@ def check_resort(self):
             pos=n_part * [(0, 0, 0)], type=n_part * [1])
 
         # And now change their positions
-        particles.pos = self.system.box_l * \
-            np.random.random((n_part, 3))
+        particles.pos = self.system.box_l * np.random.random((n_part, 3))
+
+        # All particles should still be on node 0
+        np.testing.assert_array_equal(np.copy(particles.node), 0)
 
         # Add an interacting particle in a corner of the box
         self.system.part.add(pos=(0.01, 0.01, 0.01), type=0)
@@ -63,7 +65,7 @@ def check_resort(self):
 
         # Check that we can still access all the particles
         # This basically checks if part_node and local_particles
-        # is still in a valid state after the particle exchange
+        # are still in a valid state after the particle exchange
         self.assertEqual(sum(self.system.part.all().type), n_part)
 
         # Check that the system is still valid
@@ -74,6 +76,15 @@ def check_resort(self):
         # force calculation
         self.system.integrator.run(0, recalc_forces=True)
 
+        # Check particle transfer back to node 0
+        old_nodes = np.copy(particles.node)
+        particles.pos = n_part * [(0., 0., 0.)]
+        new_nodes = np.copy(particles.node)
+        np.testing.assert_array_equal(new_nodes, old_nodes)
+        self.system.cell_system.resort()
+        new_nodes = np.copy(particles.node)
+        np.testing.assert_array_equal(new_nodes, 0)
+
     def test_resort(self):
         self.check_resort()
 
diff --git a/testsuite/python/virtual_sites_relative.py b/testsuite/python/virtual_sites_relative.py
index 0b5ead8cae6..eda97ebafe0 100644
--- a/testsuite/python/virtual_sites_relative.py
+++ b/testsuite/python/virtual_sites_relative.py
@@ -159,6 +159,8 @@ def test_vs_exceptions(self):
         # relating to anything else other than a particle or id is not allowed
         with self.assertRaisesRegex(ValueError, "Argument of 'vs_auto_relate_to' has to be of type ParticleHandle or int"):
             p2.vs_auto_relate_to('0')
+        with self.assertRaisesRegex(ValueError, "Invalid particle id: -2"):
+            p2.vs_auto_relate_to(-2)
         # relating to itself is not allowed
         with self.assertRaisesRegex(ValueError, "A virtual site cannot relate to itself"):
             p2.vs_auto_relate_to(p2)