Reaction ensemble sample with complex reaction and analytical solution

samples/reaction_ensemble.py

@@ -17,13 +17,15 @@
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
 #
 """
-Guide for the reaction ensemble and the constant pH ensemble.
+Guide for the reaction ensemble and the constant pH ensemble. The modeled
+reaction is :math:`\\mathrm{AH} \\leftrightarrow \\mathrm{A}^- + \\mathrm{H}^+`.
 """
-epilog = """You can choose in which ensemble you want to simulate via either
+epilog = __doc__.split(":math:")[0] + "AH <-> A- + H+. " + \
+    """You can choose in which ensemble you want to simulate via either
 providing --reaction_ensemble or --constant_pH_ensemble as command line
 argument to the script. Be aware that in the case of the reaction ensemble,
 the dissociation constant gamma is not the thermodynamic reaction constant K,
-but rather K * 1 mol/l and therefore carries a unit! In the case of the of the
+but rather K * 1 mol/l and therefore carries a unit! In the case of the
 constant pH method, gamma is the thermodynamic reaction constant!
 """
 import numpy as np
@@ -32,7 +34,7 @@
 import espressomd
 from espressomd import reaction_ensemble
 
-parser = argparse.ArgumentParser(epilog=__doc__ + epilog)
+parser = argparse.ArgumentParser(epilog=epilog)
 group = parser.add_mutually_exclusive_group()
 group.add_argument('--reaction_ensemble', action='store_const', dest='mode',
                    const='reaction_ensemble')

samples/reaction_ensemble_complex_reaction.py

@@ -0,0 +1,162 @@
+#
+# Copyright (C) 2013-2018 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/>.
+#
+"""
+Guide for the reaction ensemble. The modeled reaction is
+:math:`2\\mathrm{A} + 3\\mathrm{B} \\leftrightarrow 4\\mathrm{C} + 1\\mathrm{D} + 3\\mathrm{E}`.
+"""
+
+import pprint
+import numpy as np
+import scipy.optimize
+
+import espressomd
+from espressomd import reaction_ensemble
+
+# System parameters
+#############################################################
+box_l = 35.0
+volume = box_l**3
+
+# Integration parameters
+#############################################################
+system = espressomd.System(box_l=[box_l] * 3)
+np.random.seed(seed=343)
+
+system.time_step = 0.02
+system.cell_system.skin = 0.4
+
+# Particle setup
+#############################################################
+# 2A +3B <-> 4C + 1D +3 E
+N0 = 50  # number of reactant pairs
+
+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]
+types_name = {type_A: 'A', type_B: 'B', type_C: 'C', type_D: 'D', type_E: 'E'}
+
+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]
+nu_bar = np.sum(nu_s)
+
+# reaction constant
+K = 0.001
+# reference concentration for which the reaction constant is reported, 1 mol/l
+c_ref_in_mol_per_l = 1.0
+# simulation units: 1 sigma = 3.55 Angstrom
+conversion_inv_sigma_cube_to_mol_per_l = 37.1
+c_ref_in_1_div_sigma_cubed = c_ref_in_mol_per_l / \
+    conversion_inv_sigma_cube_to_mol_per_l
+Gamma = K * (c_ref_in_1_div_sigma_cubed)**nu_bar
+
+# setup N0 pairs of reactants
+for i in range(N0):
+    for j in range(abs(nu_A)):
+        system.part.add(pos=np.random.random(3) * system.box_l, type=type_A)
+    for j in range(abs(nu_B)):
+        system.part.add(pos=np.random.random(3) * system.box_l, type=type_B)
+
+# use an exclusion radius of 0 to simulate an ideal gas
+RE = reaction_ensemble.ReactionEnsemble(
+    temperature=1, exclusion_radius=0, seed=4)
+
+
+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})
+pprint.pprint(RE.get_status())
+
+numbers = {type_A: [], type_B: [], type_C: [], type_D: [], type_E: []}
+
+# warmup
+RE.reaction(200)
+
+for i in range(200):
+    RE.reaction(10)
+    for _type in types:
+        numbers[_type].append(system.number_of_particles(type=_type))
+
+concentrations = {}
+concentrations_95ci = {}
+for ptype in types:
+    concentrations[ptype] = np.mean(
+        numbers[ptype]) / volume * conversion_inv_sigma_cube_to_mol_per_l
+    concentrations_95ci[ptype] = 1.96 * np.std(numbers[ptype], ddof=1) / np.sqrt(len(
+        numbers[ptype])) / volume * conversion_inv_sigma_cube_to_mol_per_l
+
+
+def equations(variables):
+    c_A, c_B, c_C, c_D, c_E = variables
+    # K = c_C**nu_C * c_D**nu_D * c_E**nu_E * c_A**nu_A * c_B**nu_B
+    eq1 = K - ((c_C / c_ref_in_mol_per_l)**nu_C * (c_D / c_ref_in_mol_per_l)**nu_D
+               * (c_E / c_ref_in_mol_per_l)**nu_E * (c_A / c_ref_in_mol_per_l)**nu_A
+               * (c_B / c_ref_in_mol_per_l)**nu_B)
+    eq2 = N0 - (1.0 / abs(nu_A) * c_A / conversion_inv_sigma_cube_to_mol_per_l +
+                1.0 / nu_D * c_D / conversion_inv_sigma_cube_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)
+
+
+initial_guess = [
+    concentrations[type_A],
+    concentrations[type_B],
+    concentrations[type_C],
+    concentrations[type_D],
+    concentrations[type_E]]
+
+c_A, c_B, c_C, c_D, c_E = scipy.optimize.fsolve(equations, initial_guess)
+
+concentrations_numerical = {
+    type_A: c_A,
+    type_B: c_B,
+    type_C: c_C,
+    type_D: c_D,
+    type_E: c_E}
+
+print("concentrations sampled with the reaction ensemble vs. analytical solutions:")
+for ptype in types:
+    print("  type {}: {:.4f} +/- {:.4f} mol/l (95% CI), expected: {:.4f} mol/l"
+          .format(types_name[ptype],
+                  concentrations[ptype],
+                  concentrations_95ci[ptype],
+                  concentrations_numerical[ptype]))
+
+K_sim = ((concentrations[type_C] / c_ref_in_mol_per_l)**nu_C
+         * (concentrations[type_D] / c_ref_in_mol_per_l)**nu_D
+         * (concentrations[type_E] / c_ref_in_mol_per_l)**nu_E
+         * (concentrations[type_A] / c_ref_in_mol_per_l)**nu_A
+         * (concentrations[type_B] / c_ref_in_mol_per_l)**nu_B)
+N0_sim = (1.0 / abs(nu_A) * concentrations[type_A] + 1.0 / nu_D *
+          concentrations[type_D]) / conversion_inv_sigma_cube_to_mol_per_l * volume
+print("properties of the simulated ensemble:")
+print("  K_sim = {:.1e} mol/l, expected: {:.1e} mol/l".format(K_sim, K))
+print("  N0_sim = {:.1f}, expected: {}".format(N0_sim, N0))

testsuite/scripts/samples/CMakeLists.txt

@@ -54,6 +54,7 @@ sample_test(FILE test_p3m.py SUFFIX cpu)
 sample_test(FILE test_p3m.py SUFFIX gpu LABELS "gpu")
 sample_test(FILE test_reaction_ensemble.py SUFFIX constant_pH_ensemble)
 sample_test(FILE test_reaction_ensemble.py SUFFIX reaction_ensemble)
+sample_test(FILE test_reaction_ensemble_complex_reaction.py)
 sample_test(FILE test_rigid_body.py)
 sample_test(FILE test_save_checkpoint.py)
 set_tests_properties(sample_save_checkpoint PROPERTIES FIXTURES_SETUP

testsuite/scripts/samples/test_reaction_ensemble_complex_reaction.py

@@ -0,0 +1,43 @@
+# Copyright (C) 2020 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
+
+sample, skipIfMissingFeatures = importlib_wrapper.configure_and_import(
+    "@SAMPLES_DIR@/reaction_ensemble_complex_reaction.py", random_seeds=False)
+
+
+@skipIfMissingFeatures
+class Sample(ut.TestCase):
+    system = sample.system
+
+    def test_concentrations(self):
+        err_msg = "Concentration of species {} doesn't match analytical result"
+        for ptype in sample.types:
+            self.assertAlmostEqual(sample.concentrations[ptype],
+                                   sample.concentrations_numerical[ptype],
+                                   delta=1e-3,
+                                   msg=err_msg.format(sample.types_name[ptype]))
+            self.assertLess(sample.concentrations_95ci[ptype], 1e-3,
+                            msg="95% confidence interval too large")
+        self.assertAlmostEqual(sample.K_sim, sample.K, delta=1e-3)
+        self.assertAlmostEqual(sample.N0_sim, sample.N0, delta=1e-3)
+
+
+if __name__ == "__main__":
+    ut.main()