reaction_methods: new optional keyword exclusion_radius_per_type (#4469)

Fixes #4437

Description of changes:
- The  keyword `exclusion_radius`  has been renamed to `exclusion_range`.
- If the keyword `exclusion_radius` is provided, a keyerror warns the user that `exclusion_radius` is obsolete and the current equivalent keyword is `exclusion_range`
- The reaction methods now take `exclusion_radius_per_type` as an optional argument.
- The exclusion range between two particles is calculated as their sum of their radii in `exclusion_radius_per_type`, following the Lorentz-Berthelot combining rule. 
- 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.
- A new section have been added in `doc/sphinx/advanced_methods.rst` to explain how to couple the reaction methods to molecular dynamics, including the new features of this PR.
- New unit tests have been included in `src/core/reaction_methods/tests/ReactionAlgorithm_test.cpp` to check the behavior of `exclusion_radius_per_type`
-  The loop over all particles in `ReactionAlgorithm::check_exclusion_range`  now breaks if one particle is found within the exclusion range of the inserted particle. 
- All benchmarks, samples, tests and tutorials have been changed accordingly. 

I checked the computational efficiency using the benchmark in `maintainer/benchmarks/mc_acid_base_reservoir.py` and this PR improves the MC  computational by a 15%  compared to the current  python branch (MC timing: 1.268e-04 vs 1.395e-04, this PR vs current python, respectively). Probably the gain could be even higher for dense systems.

doc/sphinx/advanced_methods.rst

@@ -1772,6 +1772,25 @@ be converted to :math:`K_c` as
 where :math:`p^{\ominus}=1\,\mathrm{atm}` is the standard pressure.
 Consider using the python module pint for unit conversion.
 
+Coupling Reaction Methods to Molecular Dynamics
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+The Monte Carlo (MC) sampling of the reaction can  be coupled with a configurational sampling using Molecular Dynamics (MD).
+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.
+
 .. _Grand canonical ensemble simulation using the Reaction Ensemble:
 
 Grand canonical ensemble simulation
@@ -1813,7 +1832,7 @@ As before in the Reaction Ensemble one can define multiple reactions (e.g. for a
 .. code-block:: python
 
     cpH=reaction_ensemble.ConstantpHEnsemble(
-        temperature=1, exclusion_radius=1, seed=77)
+        temperature=1, exclusion_range=1, seed=77)
     cpH.add_reaction(gamma=K_diss, reactant_types=[0], reactant_coefficients=[1],
                     product_types=[1, 2], product_coefficients=[1, 1],
                     default_charges={0: 0, 1: -1, 2: +1})

doc/tutorials/constant_pH/constant_pH.ipynb

@@ -561,7 +561,7 @@
     "After the particles have been added to the system we initialize the `espressomd.reaction_ensemble`. The parameters to set are:\n",
     "\n",
     "* `kT` specifies the $k_\\mathrm{B}T$ value which is used as the inverse-temperature in the Boltzmann-factor to calculate the probabilities for the insertion.\n",
-    "* `exclusion_radius` specifies the distance around particles already existing in the system, where new particles will not be inserted. This ensures that there are no big overlaps with the newly inserted particles, so that MD integration does not crash. The value of the exclusion radius does not affect the limiting result, however, it affects the convergence and the stability of the integration. For interacting systems, exclusion radius should be comparable to the particle size. For non-interacting systems, we can set the exclusion radius to $0.0$.\n",
+    "* `exclusion_range` specifies the distance around particles already existing in the system, where new particles will not be inserted. This ensures that there are no big overlaps with the newly inserted particles, so that MD integration does not crash. The value of the exclusion range does not affect the limiting result, however, it affects the convergence and the stability of the integration. For interacting systems, exclusion range should be comparable to the particle size. For non-interacting systems, we can set the exclusion range to $0.0$.\n",
     "\n",
     "* `seed` for the random number generator (RNG) is an arbitrary number that defines the starting point in the RNG sequence. Using the same seed should yield reproducible simualtion results if the same RNG implementation is used. When running multiple copies of the same simulation, each of them should use a different seed.\n",
     "\n",
@@ -577,10 +577,10 @@
    },
    "source": [
     "```python\n",
-    "exclusion_radius = PARTICLE_SIZE_REDUCED if USE_WCA else 0.0\n",
+    "exclusion_range = PARTICLE_SIZE_REDUCED if USE_WCA else 0.0\n",
     "RE = espressomd.reaction_ensemble.ConstantpHEnsemble(\n",
     "    kT=KT_REDUCED,\n",
-    "    exclusion_radius=exclusion_radius,\n",
+    "    exclusion_range=exclusion_range,\n",
     "    seed=77,\n",
     "    constant_pH=2  # temporary value\n",
     ")\n",

maintainer/benchmarks/mc_acid_base_reservoir.py

@@ -227,10 +227,10 @@ def calc_donnan_coefficient(c_acid, I_res, charge=-1):
 system.actors.add(coulomb)
 
 # ### Set up the constant pH ensemble using the reaction ensemble module
-exclusion_radius = PARTICLE_SIZE_REDUCED
+exclusion_range = PARTICLE_SIZE_REDUCED
 RE = espressomd.reaction_ensemble.ReactionEnsemble(
     kT=KT_REDUCED,
-    exclusion_radius=exclusion_radius,
+    exclusion_range=exclusion_range,
     seed=77
 )
 # this parameter helps speed up the calculation in an interacting system

samples/grand_canonical.py

@@ -89,7 +89,7 @@
             epsilon=wca_eps, sigma=wca_sig)
 
 RE = espressomd.reaction_ensemble.ReactionEnsemble(
-    kT=temperature, exclusion_radius=wca_sig, seed=3)
+    kT=temperature, exclusion_range=wca_sig, seed=3)
 RE.add_reaction(
     gamma=cs_bulk**2 * np.exp(excess_chemical_potential_pair / temperature),
     reactant_types=[], reactant_coefficients=[], product_types=[1, 2],

samples/reaction_ensemble.py

@@ -80,7 +80,7 @@
 if args.mode == "reaction_ensemble":
     RE = espressomd.reaction_ensemble.ReactionEnsemble(
         kT=1,
-        exclusion_radius=1,
+        exclusion_range=1,
         seed=77)
     RE.add_reaction(gamma=K_diss,
                     reactant_types=[types["HA"]],
@@ -90,7 +90,7 @@
                     default_charges=charge_dict)
 elif args.mode == "constant_pH_ensemble":
     RE = espressomd.reaction_ensemble.ConstantpHEnsemble(
-        kT=1, exclusion_radius=1, seed=77, constant_pH=2)
+        kT=1, exclusion_range=1, seed=77, constant_pH=2)
     RE.add_reaction(gamma=K_diss, reactant_types=[types["HA"]],
                     product_types=[types["A-"], types["H+"]],
                     default_charges=charge_dict)
@@ -106,6 +106,7 @@
 # up the simulation
 RE.set_non_interacting_type(type=max(types.values()) + 1)
 
+
 for i in range(10000):
     RE.reaction(reaction_steps=1)
     if i % 100 == 0:

samples/reaction_ensemble_complex_reaction.py

@@ -82,7 +82,7 @@
 
 # use an exclusion radius of 0 to simulate an ideal gas
 RE = espressomd.reaction_ensemble.ReactionEnsemble(
-    kT=1, exclusion_radius=0, seed=4)
+    kT=1, exclusion_range=0, seed=4)
 
 
 RE.add_reaction(

src/core/reaction_methods/ConstantpHEnsemble.hpp

@@ -40,9 +40,10 @@ namespace ReactionMethods {
  */
 class ConstantpHEnsemble : public ReactionAlgorithm {
 public:
-  ConstantpHEnsemble(int seed, double kT, double exclusion_radius,
-                     double constant_pH)
-      : ReactionAlgorithm(seed, kT, exclusion_radius),
+  ConstantpHEnsemble(
+      int seed, double kT, double exclusion_range, double constant_pH,
+      const std::unordered_map<int, double> &exclusion_radius_per_type)
+      : ReactionAlgorithm(seed, kT, exclusion_range, exclusion_radius_per_type),
         m_constant_pH(constant_pH) {}
   double m_constant_pH;
 

src/core/reaction_methods/ReactionAlgorithm.cpp

@@ -185,7 +185,7 @@ ReactionAlgorithm::make_reaction_attempt(
                               current_reaction.reactant_coefficients[i];
            j++) {
         int p_id = create_particle(current_reaction.product_types[i]);
-        check_exclusion_radius(p_id);
+        check_exclusion_range(p_id);
         p_ids_created_particles.push_back(p_id);
       }
     } else if (current_reaction.reactant_coefficients[i] -
@@ -196,7 +196,7 @@ ReactionAlgorithm::make_reaction_attempt(
            j++) {
         append_particle_property_of_random_particle(
             current_reaction.reactant_types[i], hidden_particles_properties);
-        check_exclusion_radius(hidden_particles_properties.back().p_id);
+        check_exclusion_range(hidden_particles_properties.back().p_id);
         hide_particle(hidden_particles_properties.back().p_id);
       }
     }
@@ -213,14 +213,14 @@ ReactionAlgorithm::make_reaction_attempt(
       for (int j = 0; j < current_reaction.reactant_coefficients[i]; j++) {
         append_particle_property_of_random_particle(
             current_reaction.reactant_types[i], hidden_particles_properties);
-        check_exclusion_radius(hidden_particles_properties.back().p_id);
+        check_exclusion_range(hidden_particles_properties.back().p_id);
         hide_particle(hidden_particles_properties.back().p_id);
       }
     } else {
       // create additional product_types particles
       for (int j = 0; j < current_reaction.product_coefficients[i]; j++) {
         int p_id = create_particle(current_reaction.product_types[i]);
-        check_exclusion_radius(p_id);
+        check_exclusion_range(p_id);
         p_ids_created_particles.push_back(p_id);
       }
     }
@@ -270,7 +270,7 @@ void ReactionAlgorithm::generic_oneway_reaction(
     SingleReaction &current_reaction, double &E_pot_old) {
 
   current_reaction.tried_moves += 1;
-  particle_inside_exclusion_radius_touched = false;
+  particle_inside_exclusion_range_touched = false;
   if (!all_reactant_particles_exist(current_reaction)) {
     // makes sure, no incomplete reaction is performed -> only need to consider
     // rollback of complete reactions
@@ -292,7 +292,7 @@ void ReactionAlgorithm::generic_oneway_reaction(
            hidden_particles_properties) =
       make_reaction_attempt(current_reaction);
 
-  auto const E_pot_new = (particle_inside_exclusion_radius_touched)
+  auto const E_pot_new = (particle_inside_exclusion_range_touched)
                              ? std::numeric_limits<double>::max()
                              : calculate_current_potential_energy_of_system();
 
@@ -375,16 +375,53 @@ void ReactionAlgorithm::hide_particle(int p_id) const {
 }
 
 /**
- * Check if the modified particle is too close to neighboring particles.
+ * Check if the inserted particle is too close to neighboring particles.
  */
-void ReactionAlgorithm::check_exclusion_radius(int p_id) {
-  if (exclusion_radius == 0.) {
-    return;
+void ReactionAlgorithm::check_exclusion_range(int inserted_particle_id) {
+
+  auto const &inserted_particle = get_particle_data(inserted_particle_id);
+
+  /* Check the excluded radius of the inserted particle */
+
+  if (exclusion_radius_per_type.count(inserted_particle.type()) != 0) {
+    if (exclusion_radius_per_type[inserted_particle.type()] == 0) {
+      return;
+    }
+  }
+
+  auto particle_ids = get_particle_ids();
+  /* remove the inserted particle id*/
+  particle_ids.erase(std::remove(particle_ids.begin(), particle_ids.end(),
+                                 inserted_particle_id),
+                     particle_ids.end());
+
+  /* Check  if the inserted particle within the excluded_range of any other
+   * particle*/
+  double excluded_distance;
+  for (const auto &particle_id : particle_ids) {
+    auto const &already_present_particle = get_particle_data(particle_id);
+    if (exclusion_radius_per_type.count(inserted_particle.type()) == 0 ||
+        exclusion_radius_per_type.count(inserted_particle.type()) == 0) {
+      excluded_distance = exclusion_range;
+    } else if (exclusion_radius_per_type[already_present_particle.type()] ==
+               0.) {
+      continue;
+    } else {
+      excluded_distance =
+          exclusion_radius_per_type[inserted_particle.type()] +
+          exclusion_radius_per_type[already_present_particle.type()];
+    }
+
+    auto const d_min =
+        box_geo
+            .get_mi_vector(already_present_particle.r.p, inserted_particle.r.p)
+            .norm();
+
+    if (d_min < excluded_distance) {
+      particle_inside_exclusion_range_touched = true;
+      return;
+    }
   }
-  auto const &p = get_particle_data(p_id);
-  auto const d_min = distto(partCfg(), p.pos(), p_id);
-  if (d_min < exclusion_radius)
-    particle_inside_exclusion_radius_touched = true;
 }
 
 /**
@@ -532,7 +569,7 @@ ReactionAlgorithm::generate_new_particle_positions(int type, int n_particles) {
     auto const prefactor = std::sqrt(kT / p.mass());
     auto const new_pos = get_random_position_in_box();
     move_particle(p_id, new_pos, prefactor);
-    check_exclusion_radius(p_id);
+    check_exclusion_range(p_id);
   }
 
   return old_positions;
@@ -544,7 +581,7 @@ ReactionAlgorithm::generate_new_particle_positions(int type, int n_particles) {
 bool ReactionAlgorithm::do_global_mc_move_for_particles_of_type(
     int type, int particle_number_of_type_to_be_changed) {
   m_tried_configurational_MC_moves += 1;
-  particle_inside_exclusion_radius_touched = false;
+  particle_inside_exclusion_range_touched = false;
 
   int particle_number_of_type = number_of_particles_with_type(type);
   if (particle_number_of_type == 0 or
@@ -558,7 +595,7 @@ bool ReactionAlgorithm::do_global_mc_move_for_particles_of_type(
   auto const original_positions = generate_new_particle_positions(
       type, particle_number_of_type_to_be_changed);
 
-  auto const E_pot_new = (particle_inside_exclusion_radius_touched)
+  auto const E_pot_new = (particle_inside_exclusion_range_touched)
                              ? std::numeric_limits<double>::max()
                              : calculate_current_potential_energy_of_system();
 

src/core/reaction_methods/ReactionAlgorithm.hpp

@@ -32,6 +32,7 @@
 #include <random>
 #include <stdexcept>
 #include <tuple>
+#include <unordered_map>
 #include <utility>
 #include <vector>
 
@@ -47,17 +48,31 @@ struct StoredParticleProperty {
 class ReactionAlgorithm {
 
 public:
-  ReactionAlgorithm(int seed, double kT, double exclusion_radius)
+  ReactionAlgorithm(
+      int seed, double kT, double exclusion_range,
+      const std::unordered_map<int, double> &exclusion_radius_per_type)
       : 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_radius < 0.) {
-      throw std::domain_error("Invalid value for 'exclusion_radius'");
+    if (exclusion_range < 0.) {
+      throw std::domain_error("Invalid value for 'exclusion_range'");
     }
     this->kT = kT;
-    this->exclusion_radius = exclusion_radius;
+    this->exclusion_range = exclusion_range;
+    for (auto const &item : exclusion_radius_per_type) {
+      auto type = item.first;
+      auto exclusion_radius = item.second;
+      if (exclusion_radius < 0) {
+        std::stringstream ss;
+        ss << "Invalid excluded_radius value for type " << type
+           << " value: " << exclusion_radius;
+        std::string error_message = ss.str();
+        throw std::domain_error(error_message);
+      }
+    }
+    this->exclusion_radius_per_type = exclusion_radius_per_type;
     update_volume();
   }
 
@@ -72,7 +87,8 @@ class ReactionAlgorithm {
    * infinite, therefore these configurations do not contribute
    * to the partition function and ensemble averages.
    */
-  double exclusion_radius;
+  double exclusion_range;
+  std::unordered_map<int, double> exclusion_radius_per_type;
   double volume;
   int non_interacting_type = 100;
 
@@ -84,7 +100,7 @@ class ReactionAlgorithm {
   }
 
   auto get_kT() const { return kT; }
-  auto get_exclusion_radius() const { return exclusion_radius; }
+  auto get_exclusion_range() const { return exclusion_range; }
   auto get_volume() const { return volume; }
   void set_volume(double new_volume) {
     if (new_volume <= 0.) {
@@ -114,7 +130,7 @@ class ReactionAlgorithm {
   bool do_global_mc_move_for_particles_of_type(int type,
                                                int particle_number_of_type);
 
-  bool particle_inside_exclusion_radius_touched = false;
+  bool particle_inside_exclusion_range_touched = false;
 
 protected:
   std::vector<int> m_empty_p_ids_smaller_than_max_seen_particle;
@@ -146,7 +162,6 @@ class ReactionAlgorithm {
   generate_new_particle_positions(int type, int n_particles);
   void
   restore_properties(std::vector<StoredParticleProperty> const &property_list);
-
   /**
    * @brief draws a random integer from the uniform distribution in the range
    * [0,maxint-1]
@@ -178,7 +193,7 @@ class ReactionAlgorithm {
   void replace_particle(int p_id, int desired_type) const;
   int create_particle(int desired_type);
   void hide_particle(int p_id) const;
-  void check_exclusion_radius(int p_id);
+  void check_exclusion_range(int inserted_particle_id);
   void move_particle(int p_id, Utils::Vector3d const &new_pos,
                      double velocity_prefactor);
 

src/core/reaction_methods/ReactionEnsemble.hpp

@@ -36,8 +36,11 @@ namespace ReactionMethods {
  */
 class ReactionEnsemble : public ReactionAlgorithm {
 public:
-  ReactionEnsemble(int seed, double kT, double exclusion_radius)
-      : ReactionAlgorithm(seed, kT, exclusion_radius) {}
+  ReactionEnsemble(
+      int seed, double kT, double exclusion_radius,
+      const std::unordered_map<int, double> &exclusion_radius_per_type)
+      : ReactionAlgorithm(seed, kT, exclusion_radius,
+                          exclusion_radius_per_type) {}
 
 protected:
   double calculate_acceptance_probability(

src/core/reaction_methods/WidomInsertion.hpp

@@ -28,8 +28,11 @@ namespace ReactionMethods {
 /** Widom insertion method */
 class WidomInsertion : public ReactionAlgorithm {
 public:
-  WidomInsertion(int seed, double kT, double exclusion_radius)
-      : ReactionAlgorithm(seed, kT, exclusion_radius) {}
+  WidomInsertion(
+      int seed, double kT, double exclusion_range,
+      const std::unordered_map<int, double> &exclusion_radius_per_type)
+      : ReactionAlgorithm(seed, kT, exclusion_range,
+                          exclusion_radius_per_type) {}
   double calculate_particle_insertion_potential_energy(
       SingleReaction &current_reaction);
 };

src/core/reaction_methods/tests/ConstantpHEnsemble_test.cpp

@@ -50,7 +50,7 @@ BOOST_AUTO_TEST_CASE(ConstantpHEnsemble_test) {
   };
   constexpr double tol = 100 * std::numeric_limits<double>::epsilon();
 
-  ConstantpHEnsembleTest r_algo(42, 20., 0., 1.);
+  ConstantpHEnsembleTest r_algo(42, 20., 0., 1., {});
 
   // exception if no reaction was added
   BOOST_CHECK_THROW(r_algo.check_reaction_method(), std::runtime_error);