tests: Check reaction algorithm more thoroughly

The test case re-uses the reaction_ensemble_complex_reaction.py sample.

Co-authored-by: Jonas Landsgesell <[email protected]>

testsuite/python/CMakeLists.txt

@@ -173,6 +173,7 @@ python_test(FILE rotational_dynamics.py MAX_NUM_PROC 1)
 python_test(FILE script_interface.py MAX_NUM_PROC 4)
 python_test(FILE reaction_methods_interface.py MAX_NUM_PROC 1)
 python_test(FILE reaction_ensemble.py MAX_NUM_PROC 4)
+python_test(FILE reaction_complex.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_complex.py

@@ -0,0 +1,130 @@
+#
+# 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 numpy as np
+import scipy.optimize
+import espressomd
+import espressomd.reaction_methods
+
+
+class Test(ut.TestCase):
+
+    """
+    Test the reaction ensemble for reaction 2A + 3B <-> 4C + 1D + 3E.
+    """
+
+    system = espressomd.System(box_l=3 * [35.])
+    system.time_step = 0.02
+    system.cell_system.skin = 0.4
+    np.random.seed(seed=42)
+
+    def test_equilibrium_concentrations(self):
+        # reactant and product types
+        type_A = 0
+        type_B = 1
+        type_C = 2
+        type_D = 3
+        type_E = 4
+        types = [type_A, type_B, type_C, type_D, type_E]
+        # reactant and product stoichiometric coefficients
+        nu_A = -2
+        nu_B = -3
+        nu_C = 4
+        nu_D = 1
+        nu_E = 3
+        nu_s = [nu_A, nu_B, nu_C, nu_D, nu_E]
+        # reaction constant
+        K = 0.001
+        # reference concentration: 1 mol/l
+        c_ref = 1.0
+        # simulation units: 1 sigma = 3.55 Angstrom
+        conv_sim_unit_to_mol_per_l = 37.1
+        Gamma = K * (c_ref / conv_sim_unit_to_mol_per_l)**np.sum(nu_s)
+        volume = self.system.volume()
+
+        # setup N0 batches of reactants
+        N0 = 30
+        self.system.part.add(
+            pos=np.random.random((-nu_A * N0, 3)) * self.system.box_l,
+            type=-nu_A * N0 * [type_A])
+        self.system.part.add(
+            pos=np.random.random((-nu_B * N0, 3)) * self.system.box_l,
+            type=-nu_B * N0 * [type_B])
+
+        # use an exclusion radius of 0 to simulate an ideal gas
+        RE = espressomd.reaction_methods.ReactionEnsemble(
+            kT=1., exclusion_range=0., seed=42)
+
+        RE.add_reaction(
+            gamma=Gamma, reactant_types=[type_A, type_B],
+            reactant_coefficients=[abs(nu_A), abs(nu_B)],
+            product_types=[type_C, type_D, type_E],
+            product_coefficients=[nu_C, nu_D, nu_E],
+            default_charges={type_A: 0, type_B: 0, type_C: 0, type_D: 0, type_E: 0})
+
+        # Set the hidden particle type to the lowest possible number to speed
+        # up the simulation
+        RE.set_non_interacting_type(type=max(types) + 1)
+
+        # warmup
+        RE.reaction(reaction_steps=50)
+
+        numbers = {type_A: [], type_B: [], type_C: [], type_D: [], type_E: []}
+        for _ in range(40):
+            RE.reaction(reaction_steps=6)
+            for key in types:
+                numbers[key].append(self.system.number_of_particles(type=key))
+
+        def calculate_K(c_s, nu_s):
+            return np.prod(
+                [(c_i / c_ref)**nu_i for c_i, nu_i in zip(c_s, nu_s)])
+
+        def equations(concentrations):
+            c_A, c_B, c_C, c_D, c_E = concentrations
+            eq1 = K - calculate_K(concentrations, nu_s)
+            eq2 = N0 - (1.0 / abs(nu_A) * c_A / conv_sim_unit_to_mol_per_l +
+                        1.0 / nu_D * c_D / conv_sim_unit_to_mol_per_l) * volume
+            eq3 = c_A / c_B - float(nu_A) / float(nu_B)
+            eq4 = c_C / c_D - float(nu_C) / float(nu_D)
+            eq5 = c_C / c_E - float(nu_C) / float(nu_E)
+            return (eq1, eq2, eq3, eq4, eq5)
+
+        concentrations_sim = {
+            type_i: np.mean(n_i) / volume * conv_sim_unit_to_mol_per_l
+            for type_i, n_i in numbers.items()}
+        concentrations_analytic = dict(zip(types, scipy.optimize.fsolve(
+            equations, list(concentrations_sim.values()))))
+
+        K_sim = calculate_K(concentrations_sim.values(), nu_s)
+        N0_sim = (1.0 / abs(nu_A) * concentrations_sim[type_A] + 1.0 / nu_D *
+                  concentrations_sim[type_D]) / conv_sim_unit_to_mol_per_l * volume
+
+        err_msg = "Concentration of species {} doesn't match analytical result"
+        for key in types:
+            self.assertAlmostEqual(concentrations_sim[key],
+                                   concentrations_analytic[key],
+                                   delta=1e-3,
+                                   msg=err_msg.format(key))
+        self.assertAlmostEqual(K_sim, K, delta=1e-3)
+        self.assertAlmostEqual(N0_sim, N0, delta=1e-3)
+
+
+if __name__ == "__main__":
+    ut.main()