Make reaction methods run in parallel (#4666)

Partial fix for #4614
Partial fix for #4617
Follow-up to #4637

Description of changes:
- reaction methods are now fully parallel and no longer rely on MpiCallbacks
- reaction methods are no longer entangled with the particle management code
- the Monte Carlo acceptance probability is now implemented in Python
- API changes:
   - it is no longer possible to change the reaction constant of an existing reaction with `gamma <= 0`
   - `ReactionAlgorithm.reaction()` now takes `steps` instead of `reaction_steps` as argument
   - when setting up a MC method with two or more reactions, a runtime error is raised if a reaction accidentally overwrites the default charge of a specific type with a different value
- under-the-hood changes:
   - the `particle_data.cpp` source file was removed
   - the `EspressoSystemStandAlone` class no longer relies on MpiCallbacks
   - frequency of cell system invalidation during Monte Carlo trial moves has been significantly reduced
   - the new Python interface of reaction methods is around 40% slower than the original C++ implementation

.pylintrc

@@ -67,6 +67,7 @@ disable=all
 # multiple time (only on the command line, not in the configuration file where
 # it should appear only once). See also the "--disable" option for examples.
 enable=dangerous-default-value, # W0102
+       duplicate-key, # W0109
        wildcard-import, # W0401
        assert-on-tuple, # W0199
        unused-import, # W0611

doc/tutorials/constant_pH/constant_pH.ipynb

@@ -686,7 +686,7 @@
    "source": [
     "```python\n",
     "def equilibrate_reaction(reaction_steps=1):\n",
-    "    RE.reaction(reaction_steps=reaction_steps)\n",
+    "    RE.reaction(steps=reaction_steps)\n",
     "```"
    ]
   },
@@ -743,7 +743,7 @@
     "        if USE_WCA and np.random.random() < prob_integration:\n",
     "            system.integrator.run(integration_steps)\n",
     "        # we should do at least one reaction attempt per reactive particle\n",
-    "        RE.reaction(reaction_steps=reaction_steps)        \n",
+    "        RE.reaction(steps=reaction_steps)\n",
     "        num_As[i] = system.number_of_particles(type=type_A)\n",
     "```"
    ]

maintainer/benchmarks/mc_acid_base_reservoir.py

@@ -259,7 +259,7 @@ def calc_donnan_coefficient(c_acid, I_res, charge=-1):
 
 
 def equilibrate_reaction(reaction_steps=1):
-    RE.reaction(reaction_steps=reaction_steps)
+    RE.reaction(steps=reaction_steps)
 
 
 def report_progress(system, i, next_i):
@@ -295,7 +295,7 @@ def report_progress(system, i, next_i):
 
         if MC_STEPS_PER_SAMPLE > 0:
             tick_MC = time.time()
-            RE.reaction(reaction_steps=MC_STEPS_PER_SAMPLE)
+            RE.reaction(steps=MC_STEPS_PER_SAMPLE)
             tock_MC = time.time()
 
         t_MC = (tock_MC - tick_MC) / MC_STEPS_PER_SAMPLE
@@ -332,7 +332,7 @@ def report_progress(system, i, next_i):
     for i in range(NUM_SAMPLES):
         if RUN_INTEGRATION:
             system.integrator.run(INTEGRATION_STEPS_PER_SAMPLE)
-        RE.reaction(reaction_steps=MC_STEPS_PER_SAMPLE)
+        RE.reaction(steps=MC_STEPS_PER_SAMPLE)
         n_A = system.number_of_particles(type=TYPES['A'])
         n_As.append(n_A)
         n_All = len(system.part)

samples/grand_canonical.py

@@ -101,7 +101,7 @@
 # up the simulation
 RE.set_non_interacting_type(type=max(types) + 1)
 
-RE.reaction(reaction_steps=10000)
+RE.reaction(steps=10000)
 
 p3m = espressomd.electrostatics.P3M(prefactor=2.0, accuracy=1e-3)
 system.actors.add(p3m)
@@ -134,14 +134,14 @@
 system.thermostat.set_langevin(kT=temperature, gamma=.5, seed=42)
 
 # MC warmup
-RE.reaction(reaction_steps=1000)
+RE.reaction(steps=1000)
 
 n_int_cycles = 10000
 n_int_steps = 600
 num_As = []
 deviation = None
 for i in range(n_int_cycles):
-    RE.reaction(reaction_steps=10)
+    RE.reaction(steps=10)
     system.integrator.run(steps=n_int_steps)
     num_As.append(system.number_of_particles(type=1))
     if i > 2 and i % 50 == 0:

samples/reaction_ensemble_complex_reaction.py

@@ -100,10 +100,10 @@
 RE.set_non_interacting_type(type=max(types) + 1)
 
 # warmup
-RE.reaction(reaction_steps=200)
+RE.reaction(steps=200)
 
 for i in range(200):
-    RE.reaction(reaction_steps=10)
+    RE.reaction(steps=10)
     for _type in types:
         numbers[_type].append(system.number_of_particles(type=_type))
 

samples/reaction_methods.py

@@ -107,7 +107,7 @@
 
 
 for i in range(10000):
-    RE.reaction(reaction_steps=1)
+    RE.reaction(steps=1)
     if i % 100 == 0:
         print(f"HA {system.number_of_particles(type=types['HA'])}",
               f"A- {system.number_of_particles(type=types['A-'])}",

src/core/CMakeLists.txt

@@ -41,7 +41,6 @@ add_library(
   integrate.cpp
   npt.cpp
   partCfg_global.cpp
-  particle_data.cpp
   particle_node.cpp
   polymer.cpp
   pressure.cpp

src/core/EspressoSystemStandAlone.cpp

@@ -20,7 +20,6 @@
 #include "config/config.hpp"
 
 #include "EspressoSystemStandAlone.hpp"
-#include "MpiCallbacks.hpp"
 #include "communication.hpp"
 #include "event.hpp"
 #include "grid.hpp"
@@ -47,48 +46,21 @@ EspressoSystemStandAlone::EspressoSystemStandAlone(int argc, char **argv) {
 #ifdef VIRTUAL_SITES
   set_virtual_sites(std::make_shared<VirtualSitesOff>());
 #endif
-
-  // initialize the MpiCallbacks loop (blocking on worker nodes)
-  mpi_loop();
-}
-
-static void mpi_set_box_length_local(Utils::Vector3d const &value) {
-  set_box_length(value);
-}
-
-static void mpi_set_node_grid_local(Utils::Vector3i const &value) {
-  set_node_grid(value);
-}
-
-static void mpi_set_time_step_local(double const &value) {
-  set_time_step(value);
 }
 
-REGISTER_CALLBACK(mpi_set_box_length_local)
-REGISTER_CALLBACK(mpi_set_node_grid_local)
-REGISTER_CALLBACK(mpi_set_time_step_local)
-
 void EspressoSystemStandAlone::set_box_l(Utils::Vector3d const &box_l) const {
-  if (!head_node)
-    return;
-  mpi_call_all(mpi_set_box_length_local, box_l);
+  set_box_length(box_l);
 }
 
 void EspressoSystemStandAlone::set_node_grid(
     Utils::Vector3i const &node_grid) const {
-  if (!head_node)
-    return;
-  mpi_call_all(mpi_set_node_grid_local, node_grid);
+  ::set_node_grid(node_grid);
 }
 
 void EspressoSystemStandAlone::set_time_step(double time_step) const {
-  if (!head_node)
-    return;
-  mpi_call_all(mpi_set_time_step_local, time_step);
+  ::set_time_step(time_step);
 }
 
 void EspressoSystemStandAlone::set_skin(double new_skin) const {
-  if (!head_node)
-    return;
-  mpi_call_all(mpi_set_skin_local, new_skin);
+  mpi_set_skin_local(new_skin);
 }

src/core/bond_breakage/bond_breakage.cpp

@@ -23,7 +23,6 @@
 #include "communication.hpp"
 #include "errorhandling.hpp"
 #include "event.hpp"
-#include "particle_data.hpp"
 
 #include <utils/mpi/gather_buffer.hpp>
 

src/core/collision.cpp

@@ -193,33 +193,33 @@ void Collision_parameters::initialize() {
       throw std::domain_error("Collision detection particle type for virtual "
                               "sites needs to be >=0");
     }
-    make_particle_type_exist_local(collision_params.vs_particle_type);
+    make_particle_type_exist(collision_params.vs_particle_type);
   }
 
   if (collision_params.mode == CollisionModeType::GLUE_TO_SURF) {
     if (collision_params.vs_particle_type < 0) {
       throw std::domain_error("Collision detection particle type for virtual "
                               "sites needs to be >=0");
     }
-    make_particle_type_exist_local(collision_params.vs_particle_type);
+    make_particle_type_exist(collision_params.vs_particle_type);
 
     if (collision_params.part_type_to_be_glued < 0) {
       throw std::domain_error("Collision detection particle type to be glued "
                               "needs to be >=0");
     }
-    make_particle_type_exist_local(collision_params.part_type_to_be_glued);
+    make_particle_type_exist(collision_params.part_type_to_be_glued);
 
     if (collision_params.part_type_to_attach_vs_to < 0) {
       throw std::domain_error("Collision detection particle type to attach "
                               "the virtual site to needs to be >=0");
     }
-    make_particle_type_exist_local(collision_params.part_type_to_attach_vs_to);
+    make_particle_type_exist(collision_params.part_type_to_attach_vs_to);
 
     if (collision_params.part_type_after_glueing < 0) {
       throw std::domain_error("Collision detection particle type after gluing "
                               "needs to be >=0");
     }
-    make_particle_type_exist_local(collision_params.part_type_after_glueing);
+    make_particle_type_exist(collision_params.part_type_after_glueing);
   }
 
   on_short_range_ia_change();

src/core/constraints/ShapeBasedConstraint.hpp

@@ -77,7 +77,7 @@ class ShapeBasedConstraint : public Constraint {
 
   void set_type(const int &type) {
     part_rep.type() = type;
-    make_particle_type_exist_local(type);
+    make_particle_type_exist(type);
   }
 
   Utils::Vector3d total_force() const;

src/core/electrostatics/scafacos_impl.cpp

@@ -31,7 +31,6 @@
 #include "event.hpp"
 #include "grid.hpp"
 #include "integrate.hpp"
-#include "particle_data.hpp"
 #include "tuning.hpp"
 
 #include <utils/Span.hpp>

src/core/event.cpp

@@ -84,7 +84,7 @@ void on_program_start() {
   cells_re_init(CellStructureType::CELL_STRUCTURE_REGULAR);
 
   /* make sure interaction 0<->0 always exists */
-  make_particle_type_exist_local(0);
+  make_particle_type_exist(0);
 }
 
 void on_integration_start(double time_step) {

src/core/grid.cpp

@@ -28,7 +28,6 @@
 
 #include "communication.hpp"
 #include "event.hpp"
-#include "particle_data.hpp"
 
 #include <utils/Vector.hpp>
 #include <utils/mpi/cart_comm.hpp>

src/core/magnetostatics/dlc.cpp

@@ -35,7 +35,6 @@
 #include "communication.hpp"
 #include "errorhandling.hpp"
 #include "grid.hpp"
-#include "particle_data.hpp"
 
 #include <utils/constants.hpp>
 #include <utils/math/sqr.hpp>

src/core/magnetostatics/scafacos_impl.cpp

@@ -29,7 +29,6 @@
 #include "cells.hpp"
 #include "communication.hpp"
 #include "grid.hpp"
-#include "particle_data.hpp"
 
 #include <utils/Span.hpp>
 #include <utils/Vector.hpp>

src/core/nonbonded_interactions/nonbonded_interaction_data.cpp

@@ -23,7 +23,6 @@
  */
 #include "nonbonded_interactions/nonbonded_interaction_data.hpp"
 
-#include "communication.hpp"
 #include "electrostatics/coulomb.hpp"
 #include "event.hpp"
 
@@ -50,7 +49,7 @@ static double min_global_cut = INACTIVE_CUTOFF;
  * general low-level functions
  *****************************************/
 
-void mpi_realloc_ia_params_local(int new_size) {
+static void realloc_ia_params(int new_size) {
   auto const old_size = ::max_seen_particle_type;
   if (new_size <= old_size)
     return;
@@ -76,8 +75,6 @@ void mpi_realloc_ia_params_local(int new_size) {
   ::nonbonded_ia_params = std::move(new_params);
 }
 
-REGISTER_CALLBACK(mpi_realloc_ia_params_local)
-
 static double recalc_maximal_cutoff(const IA_parameters &data) {
   auto max_cut_current = INACTIVE_CUTOFF;
 
@@ -165,19 +162,7 @@ double maximal_cutoff_nonbonded() {
   return max_cut_nonbonded;
 }
 
-bool is_new_particle_type(int type) {
-  return (type + 1) > max_seen_particle_type;
-}
-
-void make_particle_type_exist(int type) {
-  if (is_new_particle_type(type))
-    mpi_call_all(mpi_realloc_ia_params_local, type + 1);
-}
-
-void make_particle_type_exist_local(int type) {
-  if (is_new_particle_type(type))
-    mpi_realloc_ia_params_local(type + 1);
-}
+void make_particle_type_exist(int type) { realloc_ia_params(type + 1); }
 
 void set_min_global_cut(double min_global_cut) {
   ::min_global_cut = min_global_cut;

src/core/nonbonded_interactions/nonbonded_interaction_data.hpp

@@ -387,17 +387,12 @@ inline IA_parameters &get_ia_param(int i, int j) {
   return *::nonbonded_ia_params[get_ia_param_key(i, j)];
 }
 
-void mpi_realloc_ia_params_local(int new_size);
-
-bool is_new_particle_type(int type);
 /** Make sure that ia_params is large enough to cover interactions
  *  for this particle type. The interactions are initialized with values
  *  such that no physical interaction occurs.
  */
 void make_particle_type_exist(int type);
 
-void make_particle_type_exist_local(int type);
-
 /** Check if a non-bonded interaction is defined */
 inline bool checkIfInteraction(IA_parameters const &data) {
   return data.max_cut != INACTIVE_CUTOFF;