From 9bfe5bcc1ee0c0aff20c9a43a8ccc828af4a5ef4 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= <jgrad@icp.uni-stuttgart.de>
Date: Fri, 28 Jul 2023 23:46:08 +0200
Subject: [PATCH 1/2] core: Extract exclusion range neighbor search

---
 doc/sphinx/reaction_methods.rst               |  38 +++--
 src/core/reaction_methods/CMakeLists.txt      |   6 +-
 .../reaction_methods/ConstantpHEnsemble.hpp   |  13 +-
 src/core/reaction_methods/ExclusionRadius.cpp | 154 ++++++++++++++++++
 src/core/reaction_methods/ExclusionRadius.hpp |  43 +++++
 .../reaction_methods/ReactionAlgorithm.cpp    |  64 +-------
 .../reaction_methods/ReactionAlgorithm.hpp    |  36 +---
 .../reaction_methods/ReactionEnsemble.hpp     |   9 +-
 src/core/reaction_methods/WidomInsertion.hpp  |  12 +-
 .../tests/ReactionAlgorithm_test.cpp          |  31 ++--
 src/python/espressomd/reaction_methods.py     |  60 ++++++-
 .../reaction_methods/ConstantpHEnsemble.hpp   |  11 +-
 .../reaction_methods/ExclusionRadius.hpp      | 109 +++++++++++++
 .../reaction_methods/ReactionAlgorithm.cpp    |  45 +----
 .../reaction_methods/ReactionAlgorithm.hpp    |  10 ++
 .../reaction_methods/ReactionEnsemble.hpp     |   9 +-
 .../reaction_methods/WidomInsertion.hpp       |  34 ++--
 .../reaction_methods/initialize.cpp           |   2 +
 .../tests/ConstantpHEnsemble_test.cpp         |   4 +
 .../tests/ReactionEnsemble_test.cpp           |   3 +
 .../python/reaction_methods_interface.py      |  15 +-
 21 files changed, 498 insertions(+), 210 deletions(-)
 create mode 100644 src/core/reaction_methods/ExclusionRadius.cpp
 create mode 100644 src/core/reaction_methods/ExclusionRadius.hpp
 create mode 100644 src/script_interface/reaction_methods/ExclusionRadius.hpp

diff --git a/doc/sphinx/reaction_methods.rst b/doc/sphinx/reaction_methods.rst
index 5b2a1b012f8..9a0d7140eb7 100644
--- a/doc/sphinx/reaction_methods.rst
+++ b/doc/sphinx/reaction_methods.rst
@@ -353,14 +353,30 @@ The Monte Carlo (MC) sampling of the reaction can  be coupled with a configurati
 For non-interacting systems this coupling is not an issue, but for interacting systems the insertion of new particles
 can lead to instabilities in the MD integration ultimately leading to a crash of the simulation.
 
-This integration instabilities can be avoided by defining a distance around the particles which already exist in the system
-where new particles will not be inserted, which is defined by the required keyword ``exclusion_range``.
-This prevents big overlaps with the newly inserted particles, avoiding too big forces between particles, which prevents the MD integration from crashing.
-The value of the exclusion range does not affect the limiting result and it only affects the convergence and the stability of the integration.  For interacting systems,
-it is usually a good practice to choose the exclusion range such that it is comparable to the diameter of the particles.
-
-If particles with significantly different sizes are present, it is desired to define a different exclusion range for each pair of particle types. This can be done by
-defining an exclusion radius per particle type by using the optional argument ``exclusion_radius_per_type``. Then, their exclusion range is calculated using
-the Lorentz-Berthelot combination rule, *i.e.* ``exclusion_range = exclusion_radius_per_type[particle_type_1] + exclusion_radius_per_type[particle_type_2]``.
-If the exclusion radius of one particle type is not defined, the value of the parameter provided in ``exclusion_range`` is used by default.
-If the value in ``exclusion_radius_per_type`` is equal to 0, then the exclusion range of that particle type with any other particle is 0.
+This integration instabilities can be avoided by defining a distance around
+the particles which already exist in the system where new particles will not
+be inserted, which is defined by the required keyword ``exclusion_range``.
+This prevents big overlaps with the newly inserted particles, which would
+otherwise cause large forces between particles and crash the MD integrator.
+The value of the exclusion range does not affect the limiting result
+and it only affects the convergence and the stability of the integration.
+For interacting systems, it is usually a good practice to choose the exclusion
+range such that it is comparable to the diameter of the particles.
+
+If particles with significantly different sizes are present, it is desirable
+to define a different exclusion range for each pair of particle types.
+This can be done by defining an exclusion radius per particle type via
+the optional argument ``exclusion_radius_per_type``.
+Then, their exclusion range is calculated using the Lorentz-Berthelot
+combination rule, *i.e.*
+
+.. code-block:: python
+
+    exclusion_range = exclusion_radius_per_type[particle_type_1] + \
+                      exclusion_radius_per_type[particle_type_2]
+
+If the exclusion radius of one particle type is not defined, the value
+of the parameter provided in ``exclusion_range`` is used by default.
+If the value in ``exclusion_radius_per_type`` is equal to 0,
+then the exclusion range of that particle type with any other particle is 0.
+For more detail, see :class:`espressomd.cell_system.ExclusionRadius`.
diff --git a/src/core/reaction_methods/CMakeLists.txt b/src/core/reaction_methods/CMakeLists.txt
index 88a046e9b79..fe0c8971060 100644
--- a/src/core/reaction_methods/CMakeLists.txt
+++ b/src/core/reaction_methods/CMakeLists.txt
@@ -18,8 +18,10 @@
 #
 
 target_sources(
-  espresso_core PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/ReactionAlgorithm.cpp
-                        ${CMAKE_CURRENT_SOURCE_DIR}/utils.cpp)
+  espresso_core
+  PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/ReactionAlgorithm.cpp
+          ${CMAKE_CURRENT_SOURCE_DIR}/ExclusionRadius.cpp
+          ${CMAKE_CURRENT_SOURCE_DIR}/utils.cpp)
 
 if(ESPRESSO_BUILD_TESTS)
   add_subdirectory(tests)
diff --git a/src/core/reaction_methods/ConstantpHEnsemble.hpp b/src/core/reaction_methods/ConstantpHEnsemble.hpp
index ea526aa0613..5bf33b7048c 100644
--- a/src/core/reaction_methods/ConstantpHEnsemble.hpp
+++ b/src/core/reaction_methods/ConstantpHEnsemble.hpp
@@ -19,9 +19,12 @@
 #ifndef REACTION_METHODS_CONSTANT_PH_ENSEMBLE_HPP
 #define REACTION_METHODS_CONSTANT_PH_ENSEMBLE_HPP
 
+#include "reaction_methods/ExclusionRadius.hpp"
 #include "reaction_methods/ReactionAlgorithm.hpp"
 
 #include <map>
+#include <memory>
+#include <utility>
 
 namespace ReactionMethods {
 
@@ -40,12 +43,10 @@ namespace ReactionMethods {
  */
 class ConstantpHEnsemble : public ReactionAlgorithm {
 public:
-  ConstantpHEnsemble(
-      boost::mpi::communicator const &comm, int seed, double kT,
-      double exclusion_range, double constant_pH,
-      const std::unordered_map<int, double> &exclusion_radius_per_type)
-      : ReactionAlgorithm(comm, seed, kT, exclusion_range,
-                          exclusion_radius_per_type),
+  ConstantpHEnsemble(boost::mpi::communicator const &comm, int seed, double kT,
+                     std::shared_ptr<ExclusionRadius> exclusion,
+                     double constant_pH)
+      : ReactionAlgorithm(comm, seed, kT, std::move(exclusion)),
         m_constant_pH(constant_pH) {}
   double m_constant_pH;
 };
diff --git a/src/core/reaction_methods/ExclusionRadius.cpp b/src/core/reaction_methods/ExclusionRadius.cpp
new file mode 100644
index 00000000000..edd3d594ecf
--- /dev/null
+++ b/src/core/reaction_methods/ExclusionRadius.cpp
@@ -0,0 +1,154 @@
+/*
+ * Copyright (C) 2022-2023 The ESPResSo project
+ *
+ * This file is part of ESPResSo.
+ *
+ * ESPResSo is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * ESPResSo is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#include "reaction_methods/ExclusionRadius.hpp"
+
+#include "cells.hpp"
+#include "event.hpp"
+#include "grid.hpp"
+#include "particle_node.hpp"
+
+#include <boost/mpi/collectives/all_reduce.hpp>
+#include <boost/mpi/collectives/broadcast.hpp>
+#include <boost/mpi/communicator.hpp>
+
+#include <algorithm>
+#include <cassert>
+#include <functional>
+#include <stdexcept>
+#include <string>
+#include <unordered_map>
+#include <vector>
+
+static auto get_real_particle(boost::mpi::communicator const &comm, int p_id) {
+  assert(p_id >= 0);
+  auto ptr = ::cell_structure.get_local_particle(p_id);
+  if (ptr != nullptr and ptr->is_ghost()) {
+    ptr = nullptr;
+  }
+  assert(boost::mpi::all_reduce(comm, static_cast<int>(ptr != nullptr),
+                                std::plus<>()) == 1);
+  return ptr;
+}
+
+void ExclusionRadius::set_exclusion_range(double range) {
+  if (range < 0.) {
+    throw std::domain_error("Invalid value for exclusion range");
+  }
+  exclusion_range = range;
+  recalc_derived_parameters();
+}
+
+void ExclusionRadius::set_exclusion_radius_per_type(map_type const &map) {
+  for (auto const &[type, exclusion_radius] : map) {
+    if (exclusion_radius < 0.) {
+      throw std::domain_error("Invalid exclusion radius for type " +
+                              std::to_string(type) + ": radius " +
+                              std::to_string(exclusion_radius));
+    }
+  }
+  exclusion_radius_per_type = map;
+  recalc_derived_parameters();
+}
+
+void ExclusionRadius::recalc_derived_parameters() {
+  m_max_exclusion_range = exclusion_range;
+  for (auto const &item : exclusion_radius_per_type) {
+    auto const radius = item.second;
+    m_max_exclusion_range = std::max(m_max_exclusion_range, 2. * radius);
+  }
+}
+
+/**
+ * Check if an inserted particle is too close to neighboring particles.
+ */
+bool ExclusionRadius::check_exclusion_range(int p_id, int p_type) {
+
+  /* Check the exclusion radius of the inserted particle */
+  if (exclusion_radius_per_type.count(p_type) != 0) {
+    if (exclusion_radius_per_type[p_type] == 0.) {
+      return false;
+    }
+  }
+
+  auto p1_ptr = get_real_particle(m_comm, p_id);
+
+  std::vector<int> particle_ids;
+  if (neighbor_search_order_n) {
+    auto all_ids = get_particle_ids_parallel();
+    /* remove the inserted particle id */
+    all_ids.erase(std::remove(all_ids.begin(), all_ids.end(), p_id),
+                  all_ids.end());
+    particle_ids = all_ids;
+  } else {
+    on_observable_calc();
+    auto const local_ids =
+        get_short_range_neighbors(p_id, m_max_exclusion_range);
+    assert(p1_ptr == nullptr or !!local_ids);
+    if (local_ids) {
+      particle_ids = std::move(*local_ids);
+    }
+  }
+
+  auto within_exclusion_range = false;
+  if (p1_ptr != nullptr) {
+    auto &p1 = *p1_ptr;
+
+    /* Check if the inserted particle within any exclusion radius */
+    for (auto const p2_id : particle_ids) {
+      if (auto const p2_ptr = ::cell_structure.get_local_particle(p2_id)) {
+        auto const &p2 = *p2_ptr;
+        double excluded_distance;
+        if (exclusion_radius_per_type.count(p_type) == 0 or
+            exclusion_radius_per_type.count(p2.type()) == 0) {
+          excluded_distance = exclusion_range;
+        } else if (exclusion_radius_per_type[p2.type()] == 0.) {
+          continue;
+        } else {
+          excluded_distance = exclusion_radius_per_type[p_type] +
+                              exclusion_radius_per_type[p2.type()];
+        }
+
+        auto const d_min = ::box_geo.get_mi_vector(p2.pos(), p1.pos()).norm();
+
+        if (d_min < excluded_distance) {
+          within_exclusion_range = true;
+          break;
+        }
+      }
+    }
+    if (m_comm.rank() != 0) {
+      m_comm.send(0, 1, within_exclusion_range);
+    }
+  } else if (m_comm.rank() == 0) {
+    m_comm.recv(boost::mpi::any_source, 1, within_exclusion_range);
+  }
+  boost::mpi::broadcast(m_comm, within_exclusion_range, 0);
+  return within_exclusion_range;
+}
+
+bool ExclusionRadius::check_exclusion_range(int pid) {
+  int type_local = 0;
+  if (auto p = get_real_particle(m_comm, pid)) {
+    type_local = p->type();
+  }
+  auto const type =
+      boost::mpi::all_reduce(m_comm, type_local, std::plus<int>());
+  return check_exclusion_range(pid, type);
+}
diff --git a/src/core/reaction_methods/ExclusionRadius.hpp b/src/core/reaction_methods/ExclusionRadius.hpp
new file mode 100644
index 00000000000..70013120ef9
--- /dev/null
+++ b/src/core/reaction_methods/ExclusionRadius.hpp
@@ -0,0 +1,43 @@
+/*
+ * Copyright (C) 2022-2023 The ESPResSo project
+ *
+ * This file is part of ESPResSo.
+ *
+ * ESPResSo is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * ESPResSo is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#pragma once
+
+#include <boost/mpi/communicator.hpp>
+
+#include <unordered_map>
+
+class ExclusionRadius {
+  boost::mpi::communicator const &m_comm;
+  double m_max_exclusion_range = 0.;
+  void recalc_derived_parameters();
+
+public:
+  using map_type = std::unordered_map<int, double>;
+  explicit ExclusionRadius(boost::mpi::communicator const &comm)
+      : m_comm{comm} {}
+  void set_exclusion_range(double range);
+  void set_exclusion_radius_per_type(map_type const &map);
+  bool check_exclusion_range(int p_id, int p_type);
+  bool check_exclusion_range(int pid);
+
+  bool neighbor_search_order_n = true;
+  double exclusion_range = 0.;
+  map_type exclusion_radius_per_type{};
+};
diff --git a/src/core/reaction_methods/ReactionAlgorithm.cpp b/src/core/reaction_methods/ReactionAlgorithm.cpp
index f19f91b6d40..9e364577736 100644
--- a/src/core/reaction_methods/ReactionAlgorithm.cpp
+++ b/src/core/reaction_methods/ReactionAlgorithm.cpp
@@ -305,68 +305,8 @@ void ReactionAlgorithm::hide_particle(int p_id, int p_type) const {
  * Check if the inserted particle is too close to neighboring particles.
  */
 void ReactionAlgorithm::check_exclusion_range(int p_id, int p_type) {
-
-  /* Check the exclusion radius of the inserted particle */
-  if (exclusion_radius_per_type.count(p_type) != 0) {
-    if (exclusion_radius_per_type[p_type] == 0.) {
-      return;
-    }
-  }
-
-  auto p1_ptr = get_real_particle(p_id);
-
-  std::vector<int> particle_ids;
-  if (neighbor_search_order_n) {
-    auto all_ids = get_particle_ids_parallel();
-    /* remove the inserted particle id */
-    all_ids.erase(std::remove(all_ids.begin(), all_ids.end(), p_id),
-                  all_ids.end());
-    particle_ids = all_ids;
-  } else {
-    on_observable_calc();
-    auto const local_ids =
-        get_short_range_neighbors(p_id, m_max_exclusion_range);
-    assert(p1_ptr == nullptr or !!local_ids);
-    if (local_ids) {
-      particle_ids = std::move(*local_ids);
-    }
-  }
-
-  if (p1_ptr != nullptr) {
-    auto &p1 = *p1_ptr;
-
-    /* Check if the inserted particle within the exclusion radius of any other
-     * particle */
-    for (auto const p2_id : particle_ids) {
-      if (auto const p2_ptr = ::cell_structure.get_local_particle(p2_id)) {
-        auto const &p2 = *p2_ptr;
-        double excluded_distance;
-        if (exclusion_radius_per_type.count(p_type) == 0 ||
-            exclusion_radius_per_type.count(p2.type()) == 0) {
-          excluded_distance = exclusion_range;
-        } else if (exclusion_radius_per_type[p2.type()] == 0.) {
-          continue;
-        } else {
-          excluded_distance = exclusion_radius_per_type[p_type] +
-                              exclusion_radius_per_type[p2.type()];
-        }
-
-        auto const d_min = ::box_geo.get_mi_vector(p2.pos(), p1.pos()).norm();
-
-        if (d_min < excluded_distance) {
-          particle_inside_exclusion_range_touched = true;
-          break;
-        }
-      }
-    }
-    if (m_comm.rank() != 0) {
-      m_comm.send(0, 1, particle_inside_exclusion_range_touched);
-    }
-  } else if (m_comm.rank() == 0) {
-    m_comm.recv(boost::mpi::any_source, 1,
-                particle_inside_exclusion_range_touched);
-  }
-  boost::mpi::broadcast(m_comm, particle_inside_exclusion_range_touched, 0);
+  particle_inside_exclusion_range_touched |=
+      m_exclusion->check_exclusion_range(p_id, p_type);
 }
 
 /**
diff --git a/src/core/reaction_methods/ReactionAlgorithm.hpp b/src/core/reaction_methods/ReactionAlgorithm.hpp
index e617dbdebca..85d1003737c 100644
--- a/src/core/reaction_methods/ReactionAlgorithm.hpp
+++ b/src/core/reaction_methods/ReactionAlgorithm.hpp
@@ -21,6 +21,7 @@
 
 #include "config/config.hpp"
 
+#include "ExclusionRadius.hpp"
 #include "SingleReaction.hpp"
 
 #include "Particle.hpp"
@@ -50,20 +51,14 @@ class ReactionAlgorithm {
   boost::mpi::communicator const &m_comm;
 
 public:
-  ReactionAlgorithm(
-      boost::mpi::communicator const &comm, int seed, double kT,
-      double exclusion_range,
-      std::unordered_map<int, double> const &exclusion_radius_per_type)
-      : m_comm{comm}, kT{kT}, exclusion_range{exclusion_range},
+  ReactionAlgorithm(boost::mpi::communicator const &comm, int seed, double kT,
+                    std::shared_ptr<ExclusionRadius> exclusion)
+      : m_comm{comm}, kT{kT}, m_exclusion{std::move(exclusion)},
         m_generator(Random::mt19937(std::seed_seq({seed, seed, seed}))),
         m_normal_distribution(0.0, 1.0), m_uniform_real_distribution(0.0, 1.0) {
     if (kT < 0.) {
       throw std::domain_error("Invalid value for 'kT'");
     }
-    if (exclusion_range < 0.) {
-      throw std::domain_error("Invalid value for 'exclusion_range'");
-    }
-    set_exclusion_radius_per_type(exclusion_radius_per_type);
     update_volume();
   }
 
@@ -78,8 +73,7 @@ class ReactionAlgorithm {
    * infinite, therefore these configurations do not contribute
    * to the partition function and ensemble averages.
    */
-  double exclusion_range;
-  std::unordered_map<int, double> exclusion_radius_per_type;
+  std::shared_ptr<ExclusionRadius> m_exclusion;
   double volume;
   int non_interacting_type = 100;
 
@@ -91,7 +85,6 @@ class ReactionAlgorithm {
   }
 
   auto get_kT() const { return kT; }
-  auto get_exclusion_range() const { return exclusion_range; }
   auto get_volume() const { return volume; }
   void set_volume(double new_volume) {
     if (new_volume <= 0.) {
@@ -100,23 +93,6 @@ class ReactionAlgorithm {
     volume = new_volume;
   }
   void update_volume();
-  void
-  set_exclusion_radius_per_type(std::unordered_map<int, double> const &map) {
-    auto max_exclusion_range = exclusion_range;
-    for (auto const &item : map) {
-      auto const type = item.first;
-      auto const exclusion_radius = item.second;
-      if (exclusion_radius < 0.) {
-        throw std::domain_error("Invalid excluded_radius value for type " +
-                                std::to_string(type) + ": radius " +
-                                std::to_string(exclusion_radius));
-      }
-      max_exclusion_range =
-          std::max(max_exclusion_range, 2. * exclusion_radius);
-    }
-    exclusion_radius_per_type = map;
-    m_max_exclusion_range = max_exclusion_range;
-  }
 
   void remove_constraint() { m_reaction_constraint = ReactionConstraint::NONE; }
   void set_cyl_constraint(double center_x, double center_y, double radius);
@@ -136,7 +112,6 @@ class ReactionAlgorithm {
   }
 
   bool particle_inside_exclusion_range_touched = false;
-  bool neighbor_search_order_n = true;
 
 protected:
   std::vector<int> m_empty_p_ids_smaller_than_max_seen_particle;
@@ -261,7 +236,6 @@ class ReactionAlgorithm {
   double m_cyl_y = -10.0;
   double m_slab_start_z = -10.0;
   double m_slab_end_z = -10.0;
-  double m_max_exclusion_range = 0.;
 
   Particle *get_real_particle(int p_id) const;
   Particle *get_local_particle(int p_id) const;
diff --git a/src/core/reaction_methods/ReactionEnsemble.hpp b/src/core/reaction_methods/ReactionEnsemble.hpp
index f6ef81dd7c1..50d430fa919 100644
--- a/src/core/reaction_methods/ReactionEnsemble.hpp
+++ b/src/core/reaction_methods/ReactionEnsemble.hpp
@@ -21,8 +21,6 @@
 
 #include "reaction_methods/ReactionAlgorithm.hpp"
 
-#include <unordered_map>
-
 namespace ReactionMethods {
 
 /** Reaction ensemble method.
@@ -36,12 +34,7 @@ namespace ReactionMethods {
  */
 class ReactionEnsemble : public ReactionAlgorithm {
 public:
-  ReactionEnsemble(
-      boost::mpi::communicator const &comm, int seed, double kT,
-      double exclusion_radius,
-      const std::unordered_map<int, double> &exclusion_radius_per_type)
-      : ReactionAlgorithm(comm, seed, kT, exclusion_radius,
-                          exclusion_radius_per_type) {}
+  using ReactionAlgorithm::ReactionAlgorithm;
 };
 
 } // namespace ReactionMethods
diff --git a/src/core/reaction_methods/WidomInsertion.hpp b/src/core/reaction_methods/WidomInsertion.hpp
index d176b84a12b..bf9ae4dda1f 100644
--- a/src/core/reaction_methods/WidomInsertion.hpp
+++ b/src/core/reaction_methods/WidomInsertion.hpp
@@ -19,22 +19,16 @@
 #ifndef REACTION_METHODS_WIDOM_INSERTION_HPP
 #define REACTION_METHODS_WIDOM_INSERTION_HPP
 
-#include "ReactionAlgorithm.hpp"
+#include "reaction_methods/ReactionAlgorithm.hpp"
 
-#include <unordered_map>
-#include <utility>
+#include <stdexcept>
 
 namespace ReactionMethods {
 
 /** Widom insertion method */
 class WidomInsertion : public ReactionAlgorithm {
 public:
-  WidomInsertion(
-      boost::mpi::communicator const &comm, int seed, double kT,
-      double exclusion_range,
-      const std::unordered_map<int, double> &exclusion_radius_per_type)
-      : ReactionAlgorithm(comm, seed, kT, exclusion_range,
-                          exclusion_radius_per_type) {}
+  using ReactionAlgorithm::ReactionAlgorithm;
 
   double calculate_particle_insertion_potential_energy(int reaction_id) {
     auto &reaction = *reactions[reaction_id];
diff --git a/src/core/reaction_methods/tests/ReactionAlgorithm_test.cpp b/src/core/reaction_methods/tests/ReactionAlgorithm_test.cpp
index 180577afa6e..2827e74c5d9 100644
--- a/src/core/reaction_methods/tests/ReactionAlgorithm_test.cpp
+++ b/src/core/reaction_methods/tests/ReactionAlgorithm_test.cpp
@@ -72,7 +72,8 @@ BOOST_FIXTURE_TEST_CASE(ReactionAlgorithm_test, ParticleFactory) {
   auto const comm = boost::mpi::communicator();
 
   // check acceptance rate
-  auto r_algo = Testing::ReactionAlgorithm(comm, 42, 1., 0., {});
+  auto exclusion = std::make_shared<ExclusionRadius>(comm);
+  auto r_algo = Testing::ReactionAlgorithm(comm, 42, 1., exclusion);
   for (int tried_moves = 1; tried_moves < 5; ++tried_moves) {
     for (int accepted_moves = 0; accepted_moves < 5; ++accepted_moves) {
       r_algo.m_tried_configurational_MC_moves = tried_moves;
@@ -144,7 +145,7 @@ BOOST_FIXTURE_TEST_CASE(ReactionAlgorithm_test, ParticleFactory) {
     // update particle positions and velocities
     BOOST_CHECK(!r_algo.particle_inside_exclusion_range_touched);
     r_algo.particle_inside_exclusion_range_touched = false;
-    r_algo.exclusion_range = box_l;
+    r_algo.m_exclusion->exclusion_range = box_l;
     r_algo.displacement_mc_move(0, 2);
     auto const &bookkeeping = r_algo.get_old_system_state();
     BOOST_CHECK(r_algo.particle_inside_exclusion_range_touched);
@@ -191,7 +192,7 @@ BOOST_FIXTURE_TEST_CASE(ReactionAlgorithm_test, ParticleFactory) {
     BOOST_CHECK_THROW(displacement_move(-2, 1), std::domain_error);
     // force all MC moves to be rejected by picking particles inside
     // their exclusion radius
-    r_algo.exclusion_range = box_l;
+    r_algo.m_exclusion->exclusion_range = box_l;
     BOOST_REQUIRE(!displacement_move(type_A, 2));
     // check none of the particles moved
     for (auto const pid : {0, 1}) {
@@ -202,7 +203,7 @@ BOOST_FIXTURE_TEST_CASE(ReactionAlgorithm_test, ParticleFactory) {
       }
     }
     // force a MC move to be accepted by using a constant Hamiltonian
-    r_algo.exclusion_range = 0.;
+    r_algo.m_exclusion->exclusion_range = 0.;
     BOOST_REQUIRE(displacement_move(type_A, 1));
     std::vector<double> distances(2);
     // check that only one particle moved
@@ -279,16 +280,16 @@ BOOST_FIXTURE_TEST_CASE(ReactionAlgorithm_test, ParticleFactory) {
 
   {
     // domain error if negative exclusion_range is provided
-    BOOST_CHECK_THROW(Testing::ReactionAlgorithm(comm, 40, 1., -1, {}),
-                      std::domain_error);
+    auto exclusion = std::make_shared<ExclusionRadius>(comm);
+    BOOST_CHECK_THROW(exclusion->set_exclusion_range(-1.), std::domain_error);
 
     // domain error if a negative value is provided in exclusion_radius_per_type
     std::unordered_map<int, double> exclusion_radius_per_type;
-    exclusion_radius_per_type[type_A] = 1;
-    exclusion_radius_per_type[type_B] = -1;
+    exclusion_radius_per_type[type_A] = 1.;
+    exclusion_radius_per_type[type_B] = -1.;
 
     BOOST_CHECK_THROW(
-        Testing::ReactionAlgorithm(comm, 40, 1., 1, exclusion_radius_per_type),
+        exclusion->set_exclusion_radius_per_type(exclusion_radius_per_type),
         std::domain_error);
 
     espresso::system->set_box_l({1., 1., 1.});
@@ -299,9 +300,9 @@ BOOST_FIXTURE_TEST_CASE(ReactionAlgorithm_test, ParticleFactory) {
     set_particle_type(0, type_A);
     set_particle_type(1, type_B);
     exclusion_radius_per_type[type_A] = 0.1;
-    exclusion_radius_per_type[type_B] = 1;
-    auto r_algo =
-        Testing::ReactionAlgorithm(comm, 40, 1., 0, exclusion_radius_per_type);
+    exclusion_radius_per_type[type_B] = 1.;
+    exclusion->set_exclusion_radius_per_type(exclusion_radius_per_type);
+    auto r_algo = Testing::ReactionAlgorithm(comm, 40, 1., exclusion);
 
     // the new position will always be in the excluded range since the sum of
     // the radii of both particle types is larger than box length. The exclusion
@@ -312,7 +313,7 @@ BOOST_FIXTURE_TEST_CASE(ReactionAlgorithm_test, ParticleFactory) {
 
     // the new position will never be in the excluded range because the
     // exclusion_radius of the particle is 0
-    r_algo.exclusion_radius_per_type[type_B] = 0;
+    r_algo.m_exclusion->exclusion_radius_per_type[type_B] = 0.;
     r_algo.particle_inside_exclusion_range_touched = false;
     r_algo.displacement_mc_move(type_B, 1);
     BOOST_REQUIRE(!r_algo.particle_inside_exclusion_range_touched);
@@ -320,8 +321,8 @@ BOOST_FIXTURE_TEST_CASE(ReactionAlgorithm_test, ParticleFactory) {
 
     // the new position will never accepted since the value in exclusion_range
     // will be used if the particle does not have a defined excluded radius
-    r_algo.exclusion_range = 1;
-    r_algo.exclusion_radius_per_type = {{type_A, 0}};
+    r_algo.m_exclusion->exclusion_range = 1.;
+    r_algo.m_exclusion->exclusion_radius_per_type = {{type_A, 0.}};
     r_algo.displacement_mc_move(type_B, 1);
     BOOST_REQUIRE(r_algo.particle_inside_exclusion_range_touched);
     r_algo.clear_old_system_state();
diff --git a/src/python/espressomd/reaction_methods.py b/src/python/espressomd/reaction_methods.py
index 7cf5eb9ca17..768b1a2072e 100644
--- a/src/python/espressomd/reaction_methods.py
+++ b/src/python/espressomd/reaction_methods.py
@@ -52,6 +52,64 @@ def make_backward_reaction(self):
             product_coefficients=self.reactant_coefficients)
 
 
+@script_interface_register
+class ExclusionRadius(ScriptInterfaceHelper):
+    """
+    Neighbor search algorithm that detects when a particle enters the exclusion
+    zone of another particle. The exclusion radii are particle type-dependent.
+
+    During the neighbor search, the following cases can arise:
+
+    * the central particle per-type exclusion radius is zero: return ``False``
+    * the neighbor particle per-type exclusion radius is zero: return ``False``
+    * the central and neighbor particles per-type exclusion radii are non-zero:
+      return ``True`` if the inter-particle distance is smaller than the sum of
+      their respective exclusion radii, ``False`` otherwise
+    * either the central particle type or the neighbor particle type is not in
+      ``exclusion_radius_per_type``: return ``True`` if the inter-particle
+      distance is smaller than ``exclusion_range``, ``False`` otherwise
+
+    Attributes
+    ----------
+    exclusion_radius_per_type : :obj:`dict`, optional
+         Mapping of particle types to exclusion radii.
+    exclusion_range : :obj:`float`
+        Minimal distance from any particle whose type
+        is not in ``exclusion_radius_per_type``.
+    search_algorithm : :obj:`str`
+        Pair search algorithm. Default is ``"order_n"``, which evaluates the
+        distance between the queried particle and all other particles in the
+        system, and scales with O(N). For MPI-parallel simulations, the
+        ``"parallel"`` method is faster. The ``"parallel"`` method is not
+        recommended for simulations on 1 MPI rank, since it comes with the
+        overhead of a ghost particle update.
+
+    Methods
+    -------
+    check_exclusion_range()
+        Check the neighborhood of a central particle and detect if any neighbor
+        is too close.
+
+        Parameters
+        -----------
+        pid : :obj:`int`
+            Particle id.
+        ptype : :obj:`int`, optional
+            Particle type. If not provided, it will be read from the particle
+            and communicated to all MPI ranks.
+
+        Returns
+        -------
+        :obj:`bool` :
+            Whether the particle is within the exclusion radius
+            of another particle.
+
+    """
+    _so_name = "ReactionMethods::ExclusionRadius"
+    _so_creation_policy = "GLOBAL"
+    _so_bind_methods = ("check_exclusion_range",)
+
+
 class ReactionAlgorithm(ScriptInterfaceHelper):
     """
 
@@ -280,7 +338,7 @@ def __init__(self, **kwargs):
         if 'exclusion_radius' in kwargs:
             raise KeyError(
                 'the keyword `exclusion_radius` is obsolete. Currently, the equivalent keyword is `exclusion_range`')
-        super().__init__(**kwargs)
+        super().__init__(exclusion=ExclusionRadius(**kwargs), **kwargs)
         if not 'sip' in kwargs:
             utils.check_valid_keys(self.valid_keys(), kwargs.keys())
         self._rebuild_reaction_cache()
diff --git a/src/script_interface/reaction_methods/ConstantpHEnsemble.hpp b/src/script_interface/reaction_methods/ConstantpHEnsemble.hpp
index 4156a98ef89..77e91f9c9fa 100644
--- a/src/script_interface/reaction_methods/ConstantpHEnsemble.hpp
+++ b/src/script_interface/reaction_methods/ConstantpHEnsemble.hpp
@@ -55,18 +55,13 @@ class ConstantpHEnsemble : public ReactionAlgorithm {
   }
 
   void do_construct(VariantMap const &params) override {
+    setup_neighbor_search(params);
     context()->parallel_try_catch([&]() {
       m_re = std::make_shared<::ReactionMethods::ConstantpHEnsemble>(
           context()->get_comm(), get_value<int>(params, "seed"),
-          get_value<double>(params, "kT"),
-          get_value<double>(params, "exclusion_range"),
-          get_value<double>(params, "constant_pH"),
-          get_value_or<std::unordered_map<int, double>>(
-              params, "exclusion_radius_per_type", {}));
+          get_value<double>(params, "kT"), m_exclusion->get_instance(),
+          get_value<double>(params, "constant_pH"));
     });
-    do_set_parameter("search_algorithm",
-                     Variant{get_value_or<std::string>(
-                         params, "search_algorithm", "order_n")});
   }
 
 protected:
diff --git a/src/script_interface/reaction_methods/ExclusionRadius.hpp b/src/script_interface/reaction_methods/ExclusionRadius.hpp
new file mode 100644
index 00000000000..616f17c50e3
--- /dev/null
+++ b/src/script_interface/reaction_methods/ExclusionRadius.hpp
@@ -0,0 +1,109 @@
+/*
+ * Copyright (C) 2023 The ESPResSo project
+ *
+ * This file is part of ESPResSo.
+ *
+ * ESPResSo is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * ESPResSo is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#pragma once
+
+#include "script_interface/ScriptInterface.hpp"
+#include "script_interface/auto_parameters/AutoParameters.hpp"
+
+#include "core/reaction_methods/ExclusionRadius.hpp"
+
+#include <memory>
+#include <stdexcept>
+#include <string>
+
+namespace ScriptInterface {
+namespace ReactionMethods {
+
+class ExclusionRadius : public AutoParameters<ExclusionRadius> {
+  std::shared_ptr<::ExclusionRadius> m_obj;
+
+public:
+  ExclusionRadius() {
+    add_parameters({{"search_algorithm",
+                     [this](Variant const &v) {
+                       context()->parallel_try_catch([&]() {
+                         auto const key = get_value<std::string>(v);
+                         if (key == "order_n") {
+                           m_obj->neighbor_search_order_n = true;
+                         } else if (key == "parallel") {
+                           m_obj->neighbor_search_order_n = false;
+                         } else {
+                           throw std::invalid_argument(
+                               "Unknown search algorithm '" + key + "'");
+                         }
+                       });
+                     },
+                     [this]() {
+                       if (m_obj->neighbor_search_order_n) {
+                         return std::string("order_n");
+                       }
+                       return std::string("parallel");
+                     }},
+                    {"exclusion_range",
+                     [this](Variant const &v) {
+                       context()->parallel_try_catch([&]() {
+                         m_obj->set_exclusion_range(get_value<double>(v));
+                       });
+                     },
+                     [this]() { return m_obj->exclusion_range; }},
+                    {"exclusion_radius_per_type",
+                     [this](Variant const &v) {
+                       context()->parallel_try_catch([&]() {
+                         m_obj->set_exclusion_radius_per_type(
+                             get_value<::ExclusionRadius::map_type>(v));
+                       });
+                     },
+                     [this]() {
+                       return make_unordered_map_of_variants(
+                           m_obj->exclusion_radius_per_type);
+                     }}});
+  }
+
+  void do_construct(VariantMap const &params) override {
+    context()->parallel_try_catch([&]() {
+      m_obj = std::make_shared<::ExclusionRadius>(context()->get_comm());
+    });
+    if (params.count("exclusion_range")) {
+      do_set_parameter("exclusion_range", params.at("exclusion_range"));
+    }
+    if (params.count("exclusion_radius_per_type")) {
+      do_set_parameter("exclusion_radius_per_type",
+                       params.at("exclusion_radius_per_type"));
+    }
+  }
+
+  Variant do_call_method(std::string const &name,
+                         VariantMap const &params) override {
+    if (name == "check_exclusion_range") {
+      auto const pid = get_value<int>(params, "pid");
+      if (params.count("ptype")) {
+        auto const ptype = get_value<int>(params, "ptype");
+        return m_obj->check_exclusion_range(pid, ptype);
+      }
+      return m_obj->check_exclusion_range(pid);
+    }
+    return {};
+  }
+
+  auto get_instance() const { return m_obj; }
+};
+
+} // namespace ReactionMethods
+} // namespace ScriptInterface
diff --git a/src/script_interface/reaction_methods/ReactionAlgorithm.cpp b/src/script_interface/reaction_methods/ReactionAlgorithm.cpp
index d11c875a4bf..7723624e27b 100644
--- a/src/script_interface/reaction_methods/ReactionAlgorithm.cpp
+++ b/src/script_interface/reaction_methods/ReactionAlgorithm.cpp
@@ -60,24 +60,9 @@ ReactionAlgorithm::ReactionAlgorithm() {
        {"kT", AutoParameter::read_only, [this]() { return RE()->get_kT(); }},
        {"search_algorithm",
         [this](Variant const &v) {
-          context()->parallel_try_catch([&]() {
-            auto const key = get_value<std::string>(v);
-            if (key == "order_n") {
-              RE()->neighbor_search_order_n = true;
-            } else if (key == "parallel") {
-              RE()->neighbor_search_order_n = false;
-            } else {
-              throw std::invalid_argument("Unknown search algorithm '" + key +
-                                          "'");
-            }
-          });
+          m_exclusion->do_set_parameter("search_algorithm", v);
         },
-        [this]() {
-          if (RE()->neighbor_search_order_n) {
-            return std::string("order_n");
-          }
-          return std::string("parallel");
-        }},
+        [this]() { return m_exclusion->get_parameter("search_algorithm"); }},
        {"particle_inside_exclusion_range_touched",
         [this](Variant const &v) {
           RE()->particle_inside_exclusion_range_touched = get_value<bool>(v);
@@ -87,32 +72,20 @@ ReactionAlgorithm::ReactionAlgorithm() {
         [this]() {
           return make_unordered_map_of_variants(RE()->charges_of_types);
         }},
-       {"exclusion_range", AutoParameter::read_only,
-        [this]() { return RE()->get_exclusion_range(); }},
+       {"exclusion_range",
+        [this](Variant const &v) {
+          m_exclusion->do_set_parameter("exclusion_range", v);
+        },
+        [this]() { return m_exclusion->get_parameter("exclusion_range"); }},
        {"exclusion_radius_per_type",
         [this](Variant const &v) {
-          context()->parallel_try_catch([&]() {
-            RE()->set_exclusion_radius_per_type(
-                get_value<std::unordered_map<int, double>>(v));
-          });
+          m_exclusion->do_set_parameter("exclusion_radius_per_type", v);
         },
         [this]() {
-          return make_unordered_map_of_variants(
-              RE()->exclusion_radius_per_type);
+          return m_exclusion->get_parameter("exclusion_radius_per_type");
         }}});
 }
 
-static auto get_real_particle(boost::mpi::communicator const &comm, int p_id) {
-  assert(p_id >= 0);
-  auto ptr = ::cell_structure.get_local_particle(p_id);
-  if (ptr != nullptr and ptr->is_ghost()) {
-    ptr = nullptr;
-  }
-  assert(boost::mpi::all_reduce(comm, static_cast<int>(ptr != nullptr),
-                                std::plus<>()) == 1);
-  return ptr;
-}
-
 Variant ReactionAlgorithm::do_call_method(std::string const &name,
                                           VariantMap const &params) {
   if (name == "calculate_factorial_expression") {
diff --git a/src/script_interface/reaction_methods/ReactionAlgorithm.hpp b/src/script_interface/reaction_methods/ReactionAlgorithm.hpp
index ae3fa822fa3..f21e3a2bde2 100644
--- a/src/script_interface/reaction_methods/ReactionAlgorithm.hpp
+++ b/src/script_interface/reaction_methods/ReactionAlgorithm.hpp
@@ -20,6 +20,7 @@
 #ifndef SCRIPT_INTERFACE_REACTION_METHODS_REACTION_ALGORITHM_HPP
 #define SCRIPT_INTERFACE_REACTION_METHODS_REACTION_ALGORITHM_HPP
 
+#include "ExclusionRadius.hpp"
 #include "SingleReaction.hpp"
 
 #include "script_interface/ScriptInterface.hpp"
@@ -41,6 +42,7 @@ class ReactionAlgorithm : public AutoParameters<ReactionAlgorithm> {
 protected:
   /** Keep track of the script interface pointer of each reaction. */
   std::vector<std::shared_ptr<SingleReaction>> m_reactions;
+  std::shared_ptr<ExclusionRadius> m_exclusion;
   /**
    * Check reaction id is within the reaction container bounds.
    * Since each reaction has a corresponding backward reaction,
@@ -58,6 +60,14 @@ class ReactionAlgorithm : public AutoParameters<ReactionAlgorithm> {
     return index;
   }
 
+  void setup_neighbor_search(VariantMap const &params) {
+    auto so_ptr = get_value<ObjectRef>(params, "exclusion");
+    m_exclusion = std::dynamic_pointer_cast<ExclusionRadius>(so_ptr);
+    m_exclusion->do_set_parameter("search_algorithm",
+                                  Variant{get_value_or<std::string>(
+                                      params, "search_algorithm", "order_n")});
+  }
+
 public:
   virtual std::shared_ptr<::ReactionMethods::ReactionAlgorithm> RE() = 0;
   virtual std::shared_ptr<::ReactionMethods::ReactionAlgorithm> const
diff --git a/src/script_interface/reaction_methods/ReactionEnsemble.hpp b/src/script_interface/reaction_methods/ReactionEnsemble.hpp
index 6d2ea320dd2..82d4b4524ec 100644
--- a/src/script_interface/reaction_methods/ReactionEnsemble.hpp
+++ b/src/script_interface/reaction_methods/ReactionEnsemble.hpp
@@ -44,17 +44,12 @@ class ReactionEnsemble : public ReactionAlgorithm {
   }
 
   void do_construct(VariantMap const &params) override {
+    setup_neighbor_search(params);
     context()->parallel_try_catch([&]() {
       m_re = std::make_shared<::ReactionMethods::ReactionEnsemble>(
           context()->get_comm(), get_value<int>(params, "seed"),
-          get_value<double>(params, "kT"),
-          get_value<double>(params, "exclusion_range"),
-          get_value_or<std::unordered_map<int, double>>(
-              params, "exclusion_radius_per_type", {}));
+          get_value<double>(params, "kT"), m_exclusion->get_instance());
     });
-    do_set_parameter("search_algorithm",
-                     Variant{get_value_or<std::string>(
-                         params, "search_algorithm", "order_n")});
   }
 
 private:
diff --git a/src/script_interface/reaction_methods/WidomInsertion.hpp b/src/script_interface/reaction_methods/WidomInsertion.hpp
index bd660f64187..d03d11e6b07 100644
--- a/src/script_interface/reaction_methods/WidomInsertion.hpp
+++ b/src/script_interface/reaction_methods/WidomInsertion.hpp
@@ -45,22 +45,31 @@ class WidomInsertion : public ReactionAlgorithm {
   }
 
   WidomInsertion() {
-    add_parameters({{"search_algorithm",
-                     [this](Variant const &) {
-                       if (context()->is_head_node()) {
-                         throw std::runtime_error(
-                             "No search algorithm for WidomInsertion");
-                       }
-                     },
-                     []() { return none; }}});
+    add_parameters(
+        {{"search_algorithm",
+          [this](Variant const &) { throw_on_exclusion_change(); },
+          []() { return none; }},
+         {"exclusion_range",
+          [this](Variant const &) { throw_on_exclusion_change(); },
+          [this]() { return m_exclusion->get_parameter("exclusion_range"); }},
+         {"exclusion_radius_per_type",
+          [this](Variant const &) { throw_on_exclusion_change(); },
+          [this]() {
+            return m_exclusion->get_parameter("exclusion_radius_per_type");
+          }}});
   }
 
   void do_construct(VariantMap const &params) override {
+    setup_neighbor_search(params);
     context()->parallel_try_catch([&]() {
+      auto exclusion = m_exclusion->get_instance();
+      if (exclusion->exclusion_range != 0. or
+          not exclusion->exclusion_radius_per_type.empty()) {
+        throw std::runtime_error("No search algorithm for WidomInsertion");
+      }
       m_re = std::make_shared<::ReactionMethods::WidomInsertion>(
           context()->get_comm(), get_value<int>(params, "seed"),
-          get_value<double>(params, "kT"), 0.,
-          std::unordered_map<int, double>{});
+          get_value<double>(params, "kT"), m_exclusion->get_instance());
     });
   }
 
@@ -80,6 +89,11 @@ class WidomInsertion : public ReactionAlgorithm {
 
 private:
   std::shared_ptr<::ReactionMethods::WidomInsertion> m_re;
+  void throw_on_exclusion_change() const {
+    if (context()->is_head_node()) {
+      throw std::runtime_error("No search algorithm for WidomInsertion");
+    }
+  }
 };
 
 } /* namespace ReactionMethods */
diff --git a/src/script_interface/reaction_methods/initialize.cpp b/src/script_interface/reaction_methods/initialize.cpp
index dcf692d228b..ad289de5ad2 100644
--- a/src/script_interface/reaction_methods/initialize.cpp
+++ b/src/script_interface/reaction_methods/initialize.cpp
@@ -19,6 +19,7 @@
 
 #include "initialize.hpp"
 
+#include "ExclusionRadius.hpp"
 #include "SingleReaction.hpp"
 
 #include "ConstantpHEnsemble.hpp"
@@ -34,6 +35,7 @@ void initialize(Utils::Factory<ObjectHandle> *om) {
   om->register_new<WidomInsertion>("ReactionMethods::WidomInsertion");
   om->register_new<ReactionEnsemble>("ReactionMethods::ReactionEnsemble");
   om->register_new<ConstantpHEnsemble>("ReactionMethods::ConstantpHEnsemble");
+  om->register_new<ExclusionRadius>("ReactionMethods::ExclusionRadius");
 }
 } // namespace ReactionMethods
 } // namespace ScriptInterface
diff --git a/src/script_interface/tests/ConstantpHEnsemble_test.cpp b/src/script_interface/tests/ConstantpHEnsemble_test.cpp
index ec3d485f7a9..19e59e17a51 100644
--- a/src/script_interface/tests/ConstantpHEnsemble_test.cpp
+++ b/src/script_interface/tests/ConstantpHEnsemble_test.cpp
@@ -28,6 +28,7 @@
 #include "script_interface/LocalContext.hpp"
 #include "script_interface/Variant.hpp"
 #include "script_interface/reaction_methods/ConstantpHEnsemble.hpp"
+#include "script_interface/reaction_methods/ExclusionRadius.hpp"
 
 #include "core/EspressoSystemStandAlone.hpp"
 #include "core/Particle.hpp"
@@ -78,6 +79,8 @@ BOOST_FIXTURE_TEST_CASE(ConstantpHEnsemble_test, ParticleFactory) {
   Utils::Factory<ScriptInterface::ObjectHandle> factory;
   factory.register_new<Testing::ConstantpHEnsemble>(
       "Testing::ConstantpHEnsemble");
+  factory.register_new<ScriptInterface::ReactionMethods::ExclusionRadius>(
+      "ExclusionRadius");
 
   auto const comm = boost::mpi::communicator();
   auto const make_algo = [&factory,
@@ -91,6 +94,7 @@ BOOST_FIXTURE_TEST_CASE(ConstantpHEnsemble_test, ParticleFactory) {
     params["constant_pH"] = pH;
     params["exclusion_range"] = exclusion_range;
     params["exclusion_radius_per_type"] = make_unordered_map_of_variants(radii);
+    params["exclusion"] = ctx->make_shared_local("ExclusionRadius", params);
     auto &&sp = ctx->make_shared_local("Testing::ConstantpHEnsemble", params);
     return std::dynamic_pointer_cast<Testing::ConstantpHEnsemble>(sp);
   };
diff --git a/src/script_interface/tests/ReactionEnsemble_test.cpp b/src/script_interface/tests/ReactionEnsemble_test.cpp
index dc29e554ef5..d1c5bb80d65 100644
--- a/src/script_interface/tests/ReactionEnsemble_test.cpp
+++ b/src/script_interface/tests/ReactionEnsemble_test.cpp
@@ -81,6 +81,8 @@ BOOST_FIXTURE_TEST_CASE(ReactionEnsemble_test, ParticleFactory) {
   factory.register_new<Testing::ReactionEnsemble>("Testing::ReactionEnsemble");
   factory.register_new<ScriptInterface::ReactionMethods::SingleReaction>(
       "SingleReaction");
+  factory.register_new<ScriptInterface::ReactionMethods::ExclusionRadius>(
+      "ExclusionRadius");
 
   auto const comm = boost::mpi::communicator();
   auto const make_algo = [&factory,
@@ -93,6 +95,7 @@ BOOST_FIXTURE_TEST_CASE(ReactionEnsemble_test, ParticleFactory) {
     params["kT"] = kT;
     params["exclusion_range"] = exclusion_range;
     params["exclusion_radius_per_type"] = make_unordered_map_of_variants(radii);
+    params["exclusion"] = ctx->make_shared_local("ExclusionRadius", params);
     auto &&sp = ctx->make_shared_local("Testing::ReactionEnsemble", params);
     return std::dynamic_pointer_cast<Testing::ReactionEnsemble>(sp);
   };
diff --git a/testsuite/python/reaction_methods_interface.py b/testsuite/python/reaction_methods_interface.py
index 8f2bfcf4276..29e5409b2f1 100644
--- a/testsuite/python/reaction_methods_interface.py
+++ b/testsuite/python/reaction_methods_interface.py
@@ -94,6 +94,9 @@ def count_by_type(types):
             self.assertAlmostEqual(
                 method.exclusion_radius_per_type[2], 0.2, delta=1e-10)
             exclusion_radius_per_type = {2: 0.2}
+            method.exclusion_range = 0.8
+            self.assertAlmostEqual(method.exclusion_range, 0.8, delta=1e-10)
+            method.exclusion_range = 0.8
         self.assertAlmostEqual(
             method.get_volume(), self.system.volume(), delta=1e-10)
         method.set_volume(volume=1.)
@@ -229,7 +232,7 @@ def test_exceptions(self):
         }
         widom = espressomd.reaction_methods.WidomInsertion(kT=1., seed=12)
         method = espressomd.reaction_methods.ReactionEnsemble(
-            kT=1.5, exclusion_range=0.8, seed=12, exclusion_radius_per_type={1: 0.1})
+            kT=1.5, exclusion_range=0.1, seed=12, exclusion_radius_per_type={1: 0.1})
         method.add_reaction(**reaction_params)
         widom.add_reaction(**reaction_params)
 
@@ -310,6 +313,10 @@ def test_exceptions(self):
             widom.calculate_particle_insertion_potential_energy(reaction_id=0)
         with self.assertRaisesRegex(RuntimeError, "No search algorithm for WidomInsertion"):
             widom.search_algorithm = "order_n"
+        with self.assertRaisesRegex(RuntimeError, "No search algorithm for WidomInsertion"):
+            widom.exclusion_range = 1.
+        with self.assertRaisesRegex(RuntimeError, "No search algorithm for WidomInsertion"):
+            widom.exclusion_radius_per_type = {1: 2.}
 
         # check other exceptions
         with self.assertRaisesRegex(ValueError, "Invalid value for 'volume'"):
@@ -346,10 +353,10 @@ def test_exceptions(self):
             method.set_non_interacting_type(type=-1)
 
         # check invalid exclusion ranges and radii
-        with self.assertRaisesRegex(ValueError, "Invalid value for 'exclusion_range'"):
+        with self.assertRaisesRegex(ValueError, "Invalid value for exclusion range"):
             espressomd.reaction_methods.ReactionEnsemble(
                 kT=1., seed=12, exclusion_range=-1.)
-        with self.assertRaisesRegex(ValueError, "Invalid excluded_radius value for type 1: radius -0.10"):
+        with self.assertRaisesRegex(ValueError, "Invalid exclusion radius for type 1: radius -0.10"):
             espressomd.reaction_methods.ReactionEnsemble(
                 kT=1., seed=12, exclusion_range=1., exclusion_radius_per_type={1: -0.1})
         with self.assertRaisesRegex(ValueError, "Unknown search algorithm 'unknown'"):
@@ -357,7 +364,7 @@ def test_exceptions(self):
                 kT=1., seed=12, exclusion_range=1., search_algorithm="unknown")
         method = espressomd.reaction_methods.ReactionEnsemble(
             kT=1., exclusion_range=1., seed=12, exclusion_radius_per_type={1: 0.1})
-        with self.assertRaisesRegex(ValueError, "Invalid excluded_radius value for type 2: radius -0.10"):
+        with self.assertRaisesRegex(ValueError, "Invalid exclusion radius for type 2: radius -0.10"):
             method.exclusion_radius_per_type = {2: -0.1}
         self.assertEqual(list(method.exclusion_radius_per_type.keys()), [1])
 

From 642309aac3ed4c4554a10531cdce7901787b6e3c Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= <jgrad@icp.uni-stuttgart.de>
Date: Fri, 28 Jul 2023 23:59:18 +0200
Subject: [PATCH 2/2] Rapid prototyping of Monte Carlo methods
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

Provide an alternative implementation of SingleReaction,
ReactionAlgorithm and ConstantpHEnsemble in pure Python
for rapid prototyping of new Monte Carlo methods.

Co-authored-by: Pablo Miguel Blanco Andrés <blancoapa@natur.cuni.cz>
---
 doc/sphinx/reaction_methods.rst               |  39 +-
 samples/monte_carlo.py                        | 519 ++++++++++++++++++
 src/python/espressomd/analyze.py              |  10 +
 src/python/espressomd/reaction_methods.py     |   2 +
 src/script_interface/analysis/Analysis.cpp    |  15 +
 testsuite/scripts/samples/CMakeLists.txt      |   4 +
 testsuite/scripts/samples/test_monte_carlo.py |  47 ++
 7 files changed, 635 insertions(+), 1 deletion(-)
 create mode 100644 samples/monte_carlo.py
 create mode 100644 testsuite/scripts/samples/test_monte_carlo.py

diff --git a/doc/sphinx/reaction_methods.rst b/doc/sphinx/reaction_methods.rst
index 9a0d7140eb7..91a93f4dd9c 100644
--- a/doc/sphinx/reaction_methods.rst
+++ b/doc/sphinx/reaction_methods.rst
@@ -379,4 +379,41 @@ If the exclusion radius of one particle type is not defined, the value
 of the parameter provided in ``exclusion_range`` is used by default.
 If the value in ``exclusion_radius_per_type`` is equal to 0,
 then the exclusion range of that particle type with any other particle is 0.
-For more detail, see :class:`espressomd.cell_system.ExclusionRadius`.
+For more detail, see :class:`~espressomd.reaction_methods.ExclusionRadius`.
+
+.. _Writing new Monte Carlo methods:
+
+Writing new Monte Carlo methods
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+Most of the logic for reaction methods is implemented at the Python level.
+The C++ core is only used for performance-relevant operations on particles.
+Hence, one can prototype new reaction methods with relative ease.
+For example, the acceptance probability for a Monte Carlo trial move
+is exposed in :meth:`ReactionAlgorithm.calculate_acceptance_probability()
+<espressomd.reaction_methods.ReactionAlgorithm.calculate_acceptance_probability>`.
+Reaction method classes override this function with their custom expression
+for the acceptance probability.
+
+Alternatively, the sample script :file:`samples/monte_carlo.py` provides
+a re-implementation of the core functionality of reaction methods in Python,
+with a focus on the :ref:`constant pH <Constant pH>` and
+:ref:`reaction ensemble <Reaction Ensemble>` methods.
+More specifically, the :class:`~espressomd.reaction_methods.SingleReaction`,
+:class:`~espressomd.reaction_methods.ReactionAlgorithm`,
+:class:`~espressomd.reaction_methods.ReactionEnsemble`, and
+:class:`~espressomd.reaction_methods.ConstantpHEnsemble`
+classes are rewritten in Python.
+
+The goal of this sample is to assist in the rapid prototyping of new Monte Carlo
+methods. In particular, the sampling and move generation schemes are expressed
+in Python, and can be leveraged by users without C++ programming experience.
+The sample is designed to run with the :ref:`kernprof` profiler attached:
+
+.. code-block:: bash
+
+    pypresso --kernprof monte_carlo.py --mode=core
+    pypresso --kernprof monte_carlo.py --mode=python
+
+These Python implementations are roughly four times slower
+than their corresponding core implementations.
diff --git a/samples/monte_carlo.py b/samples/monte_carlo.py
new file mode 100644
index 00000000000..c4c2f1ff482
--- /dev/null
+++ b/samples/monte_carlo.py
@@ -0,0 +1,519 @@
+#
+# Copyright (C) 2023 The ESPResSo project
+#
+# This file is part of ESPResSo.
+#
+# ESPResSo is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# ESPResSo is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+#
+
+"""
+Rapid prototyping of new Monte Carlo methods in Python.
+
+This sample provides a re-implementation of the core functionality
+of :ref:`reaction methods <Reaction methods>` in Python,
+with a focus on the :ref:`constant pH <Constant pH>` and
+:ref:`reaction ensemble <Reaction Ensemble>` methods.
+See :ref:`Writing new Monte Carlo methods` for more details.
+The sample simulates the acid-based titration of polyelectrolyte chains.
+
+The sample is designed to run with the :ref:`kernprof` profiler attached:
+
+.. code-block:: bash
+
+    pypresso --kernprof monte_carlo.py --mode=core
+    pypresso --kernprof monte_carlo.py --mode=python
+
+"""
+
+import numpy as np
+import itertools
+import argparse
+import math
+import time
+import tqdm
+import os
+
+import espressomd
+import espressomd.polymer
+import espressomd.electrostatics
+import espressomd.reaction_methods
+
+required_features = ["P3M", "WCA"]
+espressomd.assert_features(required_features)
+
+parser = argparse.ArgumentParser(
+    prog=f"pypresso --kernprof {os.path.basename(__file__)}",
+    epilog=__doc__.lstrip().split("\n", 1)[0])
+parser.add_argument("--mode", choices=["core", "python"], default="python",
+                    help="use C++ (core) or Python (python) implementation")
+parser.add_argument("--method", choices=["cph", "re"], default="cph",
+                    help="use constant pH (cph) or reaction ensemble (re)")
+args = parser.parse_args()
+
+
+if "line_profiler" not in dir():
+    def profile(func):
+        def wrapper(*args, **kwargs):
+            return func(*args, **kwargs)
+        return wrapper
+
+
+def factorial_Ni0_by_factorial_Ni0_plus_nu_i(nu_i, N_i0):
+    value = 1.
+    if nu_i > 0:
+        value /= math.factorial(N_i0 + nu_i) // math.factorial(N_i0)
+    elif nu_i < 0:
+        value *= math.factorial(N_i0) // math.factorial(N_i0 + nu_i)
+    return value
+
+
+class SingleReaction:
+
+    def __init__(self, **kwargs):
+        self.reactant_types = kwargs["reactant_types"]
+        self.reactant_coefficients = kwargs["reactant_coefficients"]
+        self.product_types = kwargs["product_types"]
+        self.product_coefficients = kwargs["product_coefficients"]
+        self.gamma = kwargs["gamma"]
+        self.accepted_moves = 0
+        self.trial_moves = 0
+        self.accumulator_potential_energy_difference_exponential = []
+        self.nu_bar = sum(self.product_coefficients) - \
+            sum(self.reactant_coefficients)
+
+    def get_acceptance_rate(self):
+        return self.accepted_moves / self.trial_moves
+
+    def make_backward_reaction(self):
+        return SingleReaction(
+            gamma=1. / self.gamma, reactant_types=self.product_types,
+            reactant_coefficients=self.product_coefficients,
+            product_types=self.reactant_types,
+            product_coefficients=self.reactant_coefficients)
+
+
+class ReactionAlgorithm:
+
+    def __init__(self, **kwargs):
+        self.system = kwargs["system"]
+        self.kT = kwargs["kT"]
+        self.non_interacting_type = 100
+        self.reactions = []
+        self.particle_inside_exclusion_range_touched = False
+        self.default_charges = {}
+        self.m_empty_p_ids_smaller_than_max_seen_particle = []
+        self.rng = np.random.default_rng(seed=kwargs["seed"])
+        self.exclusion = espressomd.reaction_methods.ExclusionRadius(**kwargs)
+
+    def set_non_interacting_type(self, type):
+        self.non_interacting_type = type
+
+    def add_reaction(self, **kwargs):
+        """
+        Set up a reaction in the forward and backward directions.
+        """
+        self.default_charges.update(kwargs["default_charges"])
+        forward_reaction = SingleReaction(
+            reactant_types=kwargs["reactant_types"],
+            reactant_coefficients=kwargs["reactant_coefficients"],
+            product_types=kwargs["product_types"],
+            product_coefficients=kwargs["product_coefficients"],
+            gamma=kwargs["gamma"])
+        backward_reaction = forward_reaction.make_backward_reaction()
+        self.reactions.append(forward_reaction)
+        self.reactions.append(backward_reaction)
+
+    @profile
+    def reaction(self, steps):
+        """
+        Perform reaction steps. Chemical reactions are selected at random.
+        """
+        E_pot = self.system.analysis.potential_energy()
+        random = self.rng.choice(len(self.reactions), size=steps, replace=True)
+        for i in random:
+            E_pot = self.generic_oneway_reaction(self.reactions[i], E_pot)
+
+    @profile
+    def generic_oneway_reaction(self, reaction, E_pot_old):
+        """
+        Carry out a generic one-way chemical reaction of the type
+        `A + B + ... --> C + D + ...` and return the new potential
+        energy if the trial move is accepted.
+        """
+        reaction.trial_moves += 1
+        self.particle_inside_exclusion_range_touched = False
+        if not self.all_reactant_particles_exist(reaction):
+            return E_pot_old
+
+        old_particle_numbers = self.save_old_particle_numbers(reaction)
+        p_properties_old_state = self.make_reaction_attempt(reaction)
+
+        if self.particle_inside_exclusion_range_touched:  # reject trial move
+            self.restore_system(p_properties_old_state)
+            self.particle_inside_exclusion_range_touched = False
+            return E_pot_old
+
+        E_pot_new = self.system.analysis.potential_energy()
+        E_pot_diff = E_pot_new - E_pot_old
+        bf = self.calculate_acceptance_probability(
+            reaction, E_pot_diff, old_particle_numbers)
+        reaction.accumulator_potential_energy_difference_exponential.append(
+            math.exp(-E_pot_diff / self.kT))
+        if self.rng.uniform() < bf:  # accept trial move
+            self.delete_hidden_particles(p_properties_old_state)
+            reaction.accepted_moves += 1
+            return E_pot_new
+        else:  # reject trial move
+            self.restore_system(p_properties_old_state)
+            return E_pot_old
+
+    @profile
+    def make_reaction_attempt(self, reaction):
+        """
+        Carry out a chemical reaction and save the old system state.
+        """
+        minimum_number_of_types = min(len(reaction.reactant_types),
+                                      len(reaction.product_types))
+        maximum_number_of_types = max(len(reaction.reactant_types),
+                                      len(reaction.product_types))
+        p_properties_old_state = {"changed_particles": [],
+                                  "created_particles": [],
+                                  "hidden_particles": []}
+
+        for index in range(minimum_number_of_types):
+            r_type = reaction.reactant_types[index]
+            p_type = reaction.product_types[index]
+            r_charge = self.default_charges[r_type]
+            p_charge = self.default_charges[p_type]
+
+            # change reactant particles to product particles
+            size = min(reaction.reactant_coefficients[index],
+                       reaction.product_coefficients[index])
+            for random_pid in self.get_random_pids(r_type, size):
+                p = self.system.part.by_id(random_pid)
+                p.type = p_type
+                p.q = p_charge
+                p_properties_old_state["changed_particles"].append(
+                    {"pid": random_pid, "type": r_type, "charge": r_charge})
+
+            # measure stoichiometric excess
+            delta_n = reaction.product_coefficients[index] - \
+                reaction.reactant_coefficients[index]
+
+            if delta_n > 0:
+                # create product particles
+                for _ in range(delta_n):
+                    pid = self.create_particle(p_type)
+                    self.check_exclusion_range(pid)
+                    p_properties_old_state["created_particles"].append(
+                        {"pid": pid, "type": p_type, "charge": p_charge})
+            elif delta_n < 0:
+                # hide reactant particles
+                for random_pid in self.get_random_pids(r_type, -delta_n):
+                    p_properties_old_state["hidden_particles"].append(
+                        {"pid": random_pid, "type": r_type, "charge": r_charge})
+                    self.check_exclusion_range(random_pid)
+                    self.hide_particle(random_pid)
+
+        # create/hide particles with non-corresponding replacement types
+        for index in range(minimum_number_of_types, maximum_number_of_types):
+            if len(reaction.product_types) < len(reaction.reactant_types):
+                r_type = reaction.reactant_types[index]
+                r_charge = self.default_charges[r_type]
+                size = reaction.reactant_coefficients[index]
+                # hide superfluous reactant particles
+                for random_pid in self.get_random_pids(r_type, size):
+                    p_properties_old_state["hidden_particles"].append(
+                        {"pid": random_pid, "type": r_type, "charge": r_charge})
+                    self.check_exclusion_range(random_pid)
+                    self.hide_particle(random_pid)
+            else:
+                p_type = reaction.product_types[index]
+                p_charge = self.default_charges[p_type]
+                # create additional product particles
+                for _ in range(reaction.product_coefficients[index]):
+                    pid = self.create_particle(p_type)
+                    self.check_exclusion_range(pid)
+                    p_properties_old_state["created_particles"].append(
+                        {"pid": pid, "type": p_type, "charge": p_charge})
+
+        return p_properties_old_state
+
+    def calculate_acceptance_probability(
+            self, reaction, E_pot_diff, old_particle_numbers):
+        raise NotImplementedError("Derived classes must implement this method")
+
+    def check_exclusion_range(self, pid):
+        self.particle_inside_exclusion_range_touched |= self.exclusion.check_exclusion_range(
+            pid=pid)
+
+    def all_reactant_particles_exist(self, reaction):
+        for r_type in reaction.reactant_types:
+            r_index = reaction.reactant_types.index(r_type)
+            r_coef = reaction.reactant_coefficients[r_index]
+            if self.system.number_of_particles(type=r_type) < r_coef:
+                return False
+        return True
+
+    def get_random_pids(self, ptype, size):
+        pids = self.system.analysis.call_method(
+            "get_pids_of_type", ptype=ptype)
+        indices = self.rng.choice(len(pids), size=size, replace=False)
+        return [pids[i] for i in indices]
+
+    def save_old_particle_numbers(self, reaction):
+        old_particle_numbers = {}
+        for r_type in reaction.reactant_types + reaction.product_types:
+            old_particle_numbers[r_type] = self.system.number_of_particles(
+                type=r_type)
+        return old_particle_numbers
+
+    def delete_created_particles(self, p_properties_old_state):
+        for particle_info in p_properties_old_state["created_particles"]:
+            self.system.part.by_id(particle_info["pid"]).remove()
+
+    def delete_hidden_particles(self, p_properties_old_state):
+        for particle_info in p_properties_old_state["hidden_particles"]:
+            self.system.part.by_id(particle_info["pid"]).remove()
+
+    def restore_system(self, p_properties_old_state):
+        # restore properties of changed and hidden particles
+        for particle_info in p_properties_old_state["changed_particles"] + \
+                p_properties_old_state["hidden_particles"]:
+            p = self.system.part.by_id(particle_info["pid"])
+            p.type = particle_info["type"]
+            p.q = particle_info["charge"]
+        # destroy created particles
+        self.delete_created_particles(p_properties_old_state)
+
+    def hide_particle(self, pid):
+        p = self.system.part.by_id(pid)
+        p.type = self.non_interacting_type
+        p.q = 0.
+
+    def create_particle(self, ptype):
+        if len(self.m_empty_p_ids_smaller_than_max_seen_particle) == 0:
+            pid = self.system.part.highest_particle_id + 1
+        else:
+            pid = min(self.m_empty_p_ids_smaller_than_max_seen_particle)
+        self.system.part.add(id=pid, type=ptype, q=self.default_charges[ptype],
+                             pos=self.rng.random((3,)) * self.system.box_l,
+                             v=self.rng.normal(size=3) * math.sqrt(self.kT))
+        return pid
+
+
+class ReactionEnsemble(ReactionAlgorithm):
+
+    def calculate_acceptance_probability(
+            self, reaction, E_pot_diff, old_particle_numbers):
+        """
+        Calculate the acceptance probability of a Monte Carlo move.
+        """
+
+        volume = self.system.volume()
+        expr = math.exp(-E_pot_diff / self.kT)
+        expr *= volume**reaction.nu_bar * reaction.gamma
+
+        # factorial contribution of reactants
+        for i in range(len(reaction.reactant_types)):
+            nu_i = -reaction.reactant_coefficients[i]
+            N_i0 = old_particle_numbers[reaction.reactant_types[i]]
+            expr *= factorial_Ni0_by_factorial_Ni0_plus_nu_i(nu_i, N_i0)
+
+        # factorial contribution of products
+        for i in range(len(reaction.product_types)):
+            nu_i = reaction.product_coefficients[i]
+            N_i0 = old_particle_numbers[reaction.product_types[i]]
+            expr *= factorial_Ni0_by_factorial_Ni0_plus_nu_i(nu_i, N_i0)
+
+        return expr
+
+
+class ConstantpHEnsemble(ReactionAlgorithm):
+
+    def __init__(self, **kwargs):
+        super().__init__(**kwargs)
+        self.constant_pH = kwargs["constant_pH"]
+
+    def add_reaction(self, **kwargs):
+        kwargs["reactant_coefficients"] = [1]
+        kwargs["product_coefficients"] = [1, 1]
+        super().add_reaction(**kwargs)
+
+    def calculate_acceptance_probability(
+            self, reaction, E_pot_diff, old_particle_numbers):
+        """
+        Calculate the acceptance probability of a Monte Carlo move.
+        """
+
+        ln_bf = E_pot_diff - reaction.nu_bar * self.kT * math.log(10.) * (
+            self.constant_pH + reaction.nu_bar * math.log10(reaction.gamma))
+
+        factorial_expr = math.exp(-ln_bf / self.kT)
+
+        # factorial contribution of reactants
+        nu_i = -reaction.reactant_coefficients[0]
+        N_i0 = old_particle_numbers[reaction.reactant_types[0]]
+        factorial_expr *= factorial_Ni0_by_factorial_Ni0_plus_nu_i(nu_i, N_i0)
+
+        # factorial contribution of products
+        nu_i = reaction.product_coefficients[0]
+        N_i0 = old_particle_numbers[reaction.product_types[0]]
+        factorial_expr *= factorial_Ni0_by_factorial_Ni0_plus_nu_i(nu_i, N_i0)
+
+        return factorial_expr
+
+
+# System parameters
+#############################################################
+box_l = 35
+system = espressomd.System(box_l=[box_l] * 3)
+SEED = 42
+N_acid = 50
+N_chain = 5
+sigma = 1.
+epsilon = 1.
+system.time_step = 0.01
+system.cell_system.skin = 1.0
+N_steps_MD = 1000
+N_steps_MC = 50
+
+# Reaction parameters
+#############################################################
+pKa = 7.
+pH = 7.25
+kT = 1.
+friction = 1.
+types = {
+    "HA": 0,
+    "A-": 1,
+    "H+": 2,
+}
+charges = {
+    "HA": 0.,
+    "A-": -1.,
+    "H+": +1.,
+}
+params = {
+    "kT": kT,
+    "exclusion_range": sigma,
+    "seed": SEED,
+}
+if args.method == "cph":
+    params["constant_pH"] = pH
+
+# Setup
+#############################################################
+np.random.seed(seed=SEED)
+
+positions = espressomd.polymer.linear_polymer_positions(
+    n_polymers=N_chain, beads_per_chain=N_acid // N_chain,
+    bond_length=sigma + 0.1, seed=SEED)
+
+bond = espressomd.interactions.HarmonicBond(k=10., r_0=sigma + 0.1)
+system.bonded_inter.add(bond)
+for polymer_pos in positions:
+    bond_partner = None
+    for pos in polymer_pos:
+        p = system.part.add(pos=pos, type=types["A-"], q=charges["A-"])
+        if bond_partner:
+            p.add_bond((bond, bond_partner))
+        bond_partner = p
+
+for _ in range(N_acid):
+    pos = np.random.random(3) * system.box_l
+    system.part.add(pos=pos, type=types["H+"], q=charges["H+"])
+
+for type_pair in itertools.combinations_with_replacement(types.values(), 2):
+    system.non_bonded_inter[type_pair[0], type_pair[1]].wca.set_params(
+        epsilon=epsilon, sigma=sigma)
+
+p3m = espressomd.electrostatics.P3M(
+    prefactor=2., accuracy=1e-2, mesh=8, cao=3, verbose=False)
+dh = espressomd.electrostatics.DH(
+    prefactor=2., kappa=0., r_cut=0.2 * box_l)
+
+# energy minimize the system
+system.integrator.set_steepest_descent(
+    f_max=0., gamma=0.1, max_displacement=0.1 * sigma)
+system.integrator.run(200)
+system.electrostatics.solver = p3m
+system.integrator.run(1000)
+
+# thermalize the system
+system.integrator.set_vv()
+system.thermostat.set_langevin(kT=kT, gamma=friction, seed=SEED)
+system.integrator.run(1000)
+system.electrostatics.solver = dh
+
+if args.mode == "core":
+    if args.method == "cph":
+        ReactionMethod = espressomd.reaction_methods.ConstantpHEnsemble
+    elif args.method == "re":
+        ReactionMethod = espressomd.reaction_methods.ReactionEnsemble
+elif args.mode == "python":
+    if args.method == "cph":
+        ReactionMethod = ConstantpHEnsemble
+    elif args.method == "re":
+        ReactionMethod = ReactionEnsemble
+    params["system"] = system
+
+# set up reaction method
+RE = ReactionMethod(**params)
+RE.set_non_interacting_type(type=max(types.values()) + 1)
+if args.method == "cph":
+    RE.add_reaction(
+        gamma=10**-pKa,
+        reactant_types=[types["HA"]],
+        product_types=[types["A-"], types["H+"]],
+        default_charges={types[name]: charges[name] for name in types.keys()})
+elif args.method == "re":
+    RE.add_reaction(
+        gamma=1e-3,
+        reactant_types=[types["HA"]],
+        reactant_coefficients=[1],
+        product_types=[types["A-"], types["H+"]],
+        product_coefficients=[1, 1],
+        default_charges={types[name]: charges[name] for name in types.keys()})
+reaction = RE.reactions[0]
+system.setup_type_map(type_list=list(types.values()))
+
+# equilibrate the polyelectrolyte chains
+for i in range(5):
+    RE.reaction(steps=10 * N_steps_MC)
+    system.integrator.run(N_steps_MD)
+
+
+@profile
+def sample_alpha(length):
+    alpha_list = []
+    for _ in tqdm.tqdm(range(length)):
+        system.integrator.run(steps=N_steps_MD)
+        RE.reaction(steps=N_steps_MC)
+        alpha = system.number_of_particles(type=types["A-"]) / N_acid
+        alpha_list.append(alpha)
+    return alpha_list
+
+
+sample_size = 100
+tick = time.time()
+alphas = sample_alpha(sample_size)
+tock = time.time()
+
+alpha_avg = np.mean(alphas)
+alpha_err = 1.96 * np.sqrt(np.var(alphas) / len(alphas))
+acceptance_rate = reaction.get_acceptance_rate()
+print(f"acceptance rate = {100. * acceptance_rate:.0f}%")
+print(f"alpha = {alpha_avg:.2f} +/- {alpha_err:.2f}")
+print(f"runtime = {tock - tick:.2f}s")
diff --git a/src/python/espressomd/analyze.py b/src/python/espressomd/analyze.py
index bec8e00b49c..c36e98c38c3 100644
--- a/src/python/espressomd/analyze.py
+++ b/src/python/espressomd/analyze.py
@@ -380,6 +380,15 @@ class Analysis(ScriptInterfaceHelper):
             Where [0] contains the midpoints of the bins,
             and [1] contains the values of the rdf.
 
+    potential_energy()
+        Calculate the potential energy of the system, i.e. the total energy
+        minus the kinetic energy.
+
+        Returns
+        -------
+        :obj: `float`
+            Potential energy.
+
     """
     _so_name = "ScriptInterface::Analysis::Analysis"
     _so_creation_policy = "GLOBAL"
@@ -387,6 +396,7 @@ class Analysis(ScriptInterfaceHelper):
         "linear_momentum",
         "center_of_mass",
         "nbhood",
+        "potential_energy",
         "particle_neighbor_pids",
         "calc_re",
         "calc_rg",
diff --git a/src/python/espressomd/reaction_methods.py b/src/python/espressomd/reaction_methods.py
index 768b1a2072e..4bea1d8d44d 100644
--- a/src/python/espressomd/reaction_methods.py
+++ b/src/python/espressomd/reaction_methods.py
@@ -117,6 +117,8 @@ class ReactionAlgorithm(ScriptInterfaceHelper):
     the Reaction Ensemble algorithm and the constant pH method.
     Initialize the reaction algorithm by setting the
     standard pressure, temperature, and the exclusion range.
+    The exclusion range mechanism is explained in more detail
+    in :class:`~espressomd.reaction_methods.ExclusionRadius`.
 
     Note: When creating particles the velocities of the new particles are set
     according the Maxwell-Boltzmann distribution. In this step the mass of the
diff --git a/src/script_interface/analysis/Analysis.cpp b/src/script_interface/analysis/Analysis.cpp
index d7f9f20f530..6e1e3e282f9 100644
--- a/src/script_interface/analysis/Analysis.cpp
+++ b/src/script_interface/analysis/Analysis.cpp
@@ -96,6 +96,10 @@ Variant Analysis::do_call_method(std::string const &name,
     return {};
   }
 #endif
+  if (name == "potential_energy") {
+    auto const obs = calculate_energy();
+    return obs->accumulate(-obs->kinetic[0]);
+  }
   if (name == "particle_neighbor_pids") {
     on_observable_calc();
     std::unordered_map<int, std::vector<int>> dict;
@@ -109,6 +113,17 @@ Variant Analysis::do_call_method(std::string const &name,
     });
     return make_unordered_map_of_variants(dict);
   }
+  if (name == "get_pids_of_type") {
+    auto const type = get_value<int>(parameters, "ptype");
+    std::vector<int> pids;
+    for (auto const &p : ::cell_structure.local_particles()) {
+      if (p.type() == type) {
+        pids.push_back(p.id());
+      }
+    }
+    Utils::Mpi::gather_buffer(pids, context()->get_comm());
+    return pids;
+  }
 #ifdef DPD
   if (name == "dpd_stress") {
     auto const result = dpd_stress(context()->get_comm());
diff --git a/testsuite/scripts/samples/CMakeLists.txt b/testsuite/scripts/samples/CMakeLists.txt
index 78d4a9c2675..989304e1662 100644
--- a/testsuite/scripts/samples/CMakeLists.txt
+++ b/testsuite/scripts/samples/CMakeLists.txt
@@ -77,6 +77,10 @@ sample_test(FILE test_p3m.py SUFFIX gpu LABELS "gpu")
 sample_test(FILE test_reaction_methods.py SUFFIX constant_pH_ensemble)
 sample_test(FILE test_reaction_methods.py SUFFIX reaction_ensemble)
 sample_test(FILE test_reaction_ensemble_complex_reaction.py)
+sample_test(FILE test_monte_carlo.py SUFFIX core_cph)
+sample_test(FILE test_monte_carlo.py SUFFIX core_re)
+sample_test(FILE test_monte_carlo.py SUFFIX python_cph)
+sample_test(FILE test_monte_carlo.py SUFFIX python_re)
 sample_test(FILE test_rigid_body.py)
 sample_test(FILE test_save_checkpoint.py)
 set_tests_properties(sample_save_checkpoint PROPERTIES FIXTURES_SETUP
diff --git a/testsuite/scripts/samples/test_monte_carlo.py b/testsuite/scripts/samples/test_monte_carlo.py
new file mode 100644
index 00000000000..1653f43e937
--- /dev/null
+++ b/testsuite/scripts/samples/test_monte_carlo.py
@@ -0,0 +1,47 @@
+#
+# Copyright (C) 2023 The ESPResSo project
+#
+# This file is part of ESPResSo.
+#
+# ESPResSo is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# ESPResSo is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+#
+
+import unittest as ut
+import importlib_wrapper
+
+mode, method = "@TEST_SUFFIX@".split("_")
+assert method in ("cph", "re")
+
+sample, skipIfMissingFeatures = importlib_wrapper.configure_and_import(
+    "@SAMPLES_DIR@/monte_carlo.py", script_suffix="@TEST_SUFFIX@",
+    cmd_arguments=["--mode", mode, "--method", method], sample_size=150)
+
+
+@skipIfMissingFeatures
+class Sample(ut.TestCase):
+    system = sample.system
+
+    def test(self):
+        if method == "cph":
+            self.assertAlmostEqual(sample.alpha_avg, 0.29, delta=0.05)
+            self.assertAlmostEqual(sample.alpha_err, 0.01, delta=0.02)
+            self.assertAlmostEqual(sample.acceptance_rate, 0.56, delta=0.10)
+        else:
+            self.assertAlmostEqual(sample.alpha_avg, 0.33, delta=0.05)
+            self.assertAlmostEqual(sample.alpha_err, 0.01, delta=0.02)
+            self.assertAlmostEqual(sample.acceptance_rate, 0.55, delta=0.10)
+
+
+if __name__ == "__main__":
+    ut.main()