Allow multiple concurrent reaction methods (espressomd#4609)

Fixes espressomd#4595

Description of changes:
- Reaction methods now rebuild the list of free particle ids every time `ReactionAlgorithm::do_reaction()` and `WidomInsertion::calculate_particle_insertion_potential_energy()` are called
- Added a new test that checks that the reaction methods do the bookkeeping of empty ids by setting up two different instances of the reaction methods (constant pH and Widom insertion) with competing reactions.

I benchmarked my implementation against the current python branch using the `mc_acid_base_reservoir.py` script and I get the following timings: 1.564e-04 (my PR) vs  1.473e-04 (current python branch). That means that my implementation comes with a 6% performance loss in its current state.

src/core/reaction_methods/ReactionAlgorithm.cpp

@@ -55,6 +55,8 @@ void ReactionAlgorithm::do_reaction(int reaction_steps) {
   // assume that the kinetic part drops out in the process of calculating
   // ensemble averages (kinetic part may be separated and crossed out)
   auto current_E_pot = calculate_current_potential_energy_of_system();
+  // Setup the list of empty pids for bookkeeping
+  setup_bookkeeping_of_empty_pids();
   for (int i = 0; i < reaction_steps; i++) {
     int reaction_id = i_random(static_cast<int>(reactions.size()));
     generic_oneway_reaction(*reactions[reaction_id], current_E_pot);
@@ -649,4 +651,22 @@ bool ReactionAlgorithm::displacement_move_for_particles_of_type(int type,
   return false;
 }
 
+/**
+ * Cleans the list of empty pids and searches for empty pid in the system
+ */
+void ReactionAlgorithm::setup_bookkeeping_of_empty_pids() {
+  // Clean-up the list of empty pids
+  m_empty_p_ids_smaller_than_max_seen_particle.clear();
+
+  auto particle_ids = get_particle_ids();
+  std::sort(particle_ids.begin(), particle_ids.end());
+  auto pid1 = -1;
+  for (auto pid2 : particle_ids) {
+    for (int pid = pid1 + 1; pid < pid2; ++pid) {
+      m_empty_p_ids_smaller_than_max_seen_particle.push_back(pid);
+    }
+    pid1 = pid2;
+  }
+}
+
 } // namespace ReactionMethods

src/core/reaction_methods/ReactionAlgorithm.hpp

@@ -127,6 +127,7 @@ class ReactionAlgorithm {
     return {m_slab_start_z, m_slab_end_z};
   }
 
+  void setup_bookkeeping_of_empty_pids();
   void delete_particle(int p_id);
   void add_reaction(std::shared_ptr<SingleReaction> const &new_reaction);
   void delete_reaction(int reaction_id) {

src/core/reaction_methods/WidomInsertion.cpp

@@ -38,6 +38,9 @@ double WidomInsertion::calculate_particle_insertion_potential_energy(
 
   auto const E_pot_old = calculate_current_potential_energy_of_system();
 
+  // Setup the list of empty pids for bookkeeping
+  setup_bookkeeping_of_empty_pids();
+
   // make reaction attempt
   std::vector<int> p_ids_created_particles;
   std::vector<StoredParticleProperty> hidden_particles_properties;

testsuite/python/CMakeLists.txt

@@ -174,6 +174,7 @@ python_test(FILE script_interface.py MAX_NUM_PROC 4)
 python_test(FILE reaction_methods_interface.py MAX_NUM_PROC 1)
 python_test(FILE reaction_ensemble.py MAX_NUM_PROC 4)
 python_test(FILE reaction_complex.py MAX_NUM_PROC 1)
+python_test(FILE reaction_bookkeeping.py MAX_NUM_PROC 1)
 python_test(FILE widom_insertion.py MAX_NUM_PROC 1)
 python_test(FILE constant_pH.py MAX_NUM_PROC 1)
 python_test(FILE constant_pH_stats.py MAX_NUM_PROC 4 LABELS long)

testsuite/python/reaction_bookkeeping.py

@@ -0,0 +1,104 @@
+#
+# Copyright (C) 2013-2022 The ESPResSo project
+#
+# This file is part of ESPResSo.
+#
+# ESPResSo is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# ESPResSo is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+#
+
+import unittest as ut
+import unittest_decorators as utx
+
+import espressomd
+import espressomd.reaction_methods
+import numpy as np
+
+
+@utx.skipIfMissingFeatures(["ELECTROSTATICS"])
+class ReactionMethodsBookkeepingTest(ut.TestCase):
+    """
+    Test that two different instances of the reaction methods
+    do not break the particle id bookkeeping.
+    """
+    pH = 10.
+    pKa = 7.
+    exclusion_range = 1.
+    seed = 12345
+    kT = 1.
+    BOX_LENGTH = 100.
+    N_SALT = 10
+    N_acid = 10
+
+    types = {"H": 0, "Na": 1, "Cl": 2, "HA": 3, "A": 4}
+    charges = {"H": 1.0, "Na": 1.0, "Cl": -1.0, "HA": 0.0, "A": -1.0}
+
+    system = espressomd.System(box_l=[BOX_LENGTH, ] * 3)
+
+    cph = espressomd.reaction_methods.ConstantpHEnsemble(
+        constant_pH=pH,
+        kT=kT,
+        exclusion_range=exclusion_range,
+        seed=seed)
+    widom = espressomd.reaction_methods.WidomInsertion(
+        kT=kT,
+        seed=seed)
+
+    @classmethod
+    def setUpClass(cls):
+        cls.system.part.add(type=[cls.types["Na"]] * cls.N_SALT,
+                            pos=np.random.rand(cls.N_SALT, 3) * cls.BOX_LENGTH,
+                            q=[cls.charges["Na"]] * cls.N_SALT,
+                            id=list(range(20, 20 + cls.N_SALT)))
+        cls.system.part.add(type=[cls.types["Cl"]] * cls.N_SALT,
+                            pos=np.random.rand(cls.N_SALT, 3) * cls.BOX_LENGTH,
+                            q=[cls.charges["Cl"]] * cls.N_SALT)
+        cls.system.part.add(type=[cls.types["HA"]] * cls.N_acid,
+                            pos=np.random.rand(cls.N_acid, 3) * cls.BOX_LENGTH,
+                            q=[cls.charges["HA"]] * cls.N_acid)
+
+        cls.cph.add_reaction(
+            gamma=10**(-cls.pKa),
+            reactant_types=[cls.types["HA"]],
+            product_types=[cls.types["A"], cls.types["H"]],
+            default_charges={cls.types["HA"]: cls.charges["HA"],
+                             cls.types["A"]: cls.charges["A"],
+                             cls.types["H"]: cls.charges["H"]}
+        )
+
+        cls.widom.add_reaction(
+            reactant_types=[],
+            reactant_coefficients=[],
+            product_types=[cls.types["Na"], cls.types["Cl"]],
+            product_coefficients=[1, 1],
+            default_charges={cls.types["Na"]: cls.charges["Na"],
+                             cls.types["Cl"]: cls.charges["Cl"]}
+        )
+        cls.system.setup_type_map(type_list=list(cls.types.values()))
+
+    def test_reaction_bookeeping(self):
+        self.widom.calculate_particle_insertion_potential_energy(reaction_id=0)
+        self.cph.reaction(reaction_steps=100)
+
+        # Measure the chemical potential
+        for _ in range(50):
+            self.widom.calculate_particle_insertion_potential_energy(
+                reaction_id=0)
+            charge = 0.
+            for p in self.system.part:
+                charge += p.q
+            self.assertEqual(charge, 0.)
+
+
+if __name__ == "__main__":
+    ut.main()