diff --git a/Makefile b/Makefile
index 8189ada..88d40e5 100644
--- a/Makefile
+++ b/Makefile
@@ -55,6 +55,8 @@ unit_tests:
${PYTHON} testsuite/analysis_tests.py
${PYTHON} testsuite/charge_number_map_tests.py
${PYTHON} testsuite/generate_coordinates_tests.py
+ ${PYTHON} testsuite/reaction_methods_unit_tests.py
+ ${PYTHON} testsuite/determine_reservoir_concentrations_unit_test.py
functional_tests:
${PYTHON} testsuite/cph_ideal_tests.py
diff --git a/pyMBE.py b/pyMBE.py
index b8ca2de..7080851 100644
--- a/pyMBE.py
+++ b/pyMBE.py
@@ -2127,9 +2127,12 @@ def get_pka_set(self):
pka_set[name] = {'pka_value':pka_value,'acidity':acidity}
return pka_set
- def get_radius_map(self):
+ def get_radius_map(self, dimensionless=True):
'''
Gets the effective radius of each `espresso type` in `pmb.df`.
+
+ Args:
+ dimensionless('bool'): controlls if the returned radii have a dimension. Defaults to False.
Returns:
radius_map(`dict`): {espresso_type: radius}.
@@ -2142,6 +2145,9 @@ def get_radius_map(self):
state_one = pd.Series((df_state_one.sigma.values+df_state_one.offset.values)/2.0,index=df_state_one.state_one.es_type.values)
state_two = pd.Series((df_state_two.sigma.values+df_state_two.offset.values)/2.0,index=df_state_two.state_two.es_type.values)
radius_map = pd.concat([state_one,state_two],axis=0).to_dict()
+ if dimensionless:
+ for key in radius_map:
+ radius_map[key] = radius_map[key].magnitude
return radius_map
def get_reduced_units(self):
@@ -2753,7 +2759,7 @@ def setup_cpH (self, counter_ion, constant_pH, exclusion_range=None, pka_set=Non
exclusion_radius_per_type = {}
RE = reaction_methods.ConstantpHEnsemble(kT=self.kT.to('reduced_energy').magnitude,
- exclusion_range=exclusion_range.magnitude,
+ exclusion_range=exclusion_range,
seed=self.seed,
constant_pH=constant_pH,
exclusion_radius_per_type = exclusion_radius_per_type
@@ -2762,7 +2768,7 @@ def setup_cpH (self, counter_ion, constant_pH, exclusion_range=None, pka_set=Non
charge_number_map = self.get_charge_number_map()
for name in pka_set.keys():
if self.check_if_name_is_defined_in_df(name,pmb_type_to_be_defined='particle') is False :
- print('WARNING: the acid-base reaction of ' + name +' has not been set up because its espresso type is not defined in sg.type_map')
+ print('WARNING: the acid-base reaction of ' + name +' has not been set up because its espresso type is not defined in the type map.')
continue
gamma=10**-pka_set[name]['pka_value']
state_one_type = self.df.loc[self.df['name']==name].state_one.es_type.values[0]
@@ -2802,7 +2808,7 @@ def setup_gcmc(self, c_salt_res, salt_cation_name, salt_anion_name, activity_coe
exclusion_radius_per_type = {}
RE = reaction_methods.ReactionEnsemble(kT=self.kT.to('reduced_energy').magnitude,
- exclusion_range=exclusion_range.magnitude,
+ exclusion_range=exclusion_range,
seed=self.seed,
exclusion_radius_per_type = exclusion_radius_per_type
)
@@ -2874,7 +2880,7 @@ def setup_grxmc_reactions(self, pH_res, c_salt_res, proton_name, hydroxide_name,
exclusion_radius_per_type = {}
RE = reaction_methods.ReactionEnsemble(kT=self.kT.to('reduced_energy').magnitude,
- exclusion_range=exclusion_range.magnitude,
+ exclusion_range=exclusion_range,
seed=self.seed,
exclusion_radius_per_type = exclusion_radius_per_type
)
@@ -3078,7 +3084,7 @@ def setup_grxmc_unified(self, pH_res, c_salt_res, cation_name, anion_name, activ
exclusion_radius_per_type = {}
RE = reaction_methods.ReactionEnsemble(kT=self.kT.to('reduced_energy').magnitude,
- exclusion_range=exclusion_range.magnitude,
+ exclusion_range=exclusion_range,
seed=self.seed,
exclusion_radius_per_type = exclusion_radius_per_type
)
@@ -3118,7 +3124,7 @@ def setup_grxmc_unified(self, pH_res, c_salt_res, cation_name, anion_name, activ
charge_number_map = self.get_charge_number_map()
for name in pka_set.keys():
if self.check_if_name_is_defined_in_df(name,pmb_type_to_be_defined='particle') is False :
- print('WARNING: the acid-base reaction of ' + name +' has not been set up because its espresso type is not defined in sg.type_map')
+ print('WARNING: the acid-base reaction of ' + name +' has not been set up because its espresso type is not defined in the type map.')
continue
Ka = 10**-pka_set[name]['pka_value'] * self.units.mol/self.units.l
diff --git a/testsuite/determine_reservoir_concentrations_unit_test.py b/testsuite/determine_reservoir_concentrations_unit_test.py
new file mode 100644
index 0000000..ed92e43
--- /dev/null
+++ b/testsuite/determine_reservoir_concentrations_unit_test.py
@@ -0,0 +1,51 @@
+#
+# Copyright (C) 2024 pyMBE-dev team
+#
+# This file is part of pyMBE.
+#
+# pyMBE 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.
+#
+# pyMBE 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 .
+
+# Import pyMBE and other libraries
+import pyMBE
+import numpy as np
+
+def determine_reservoir_concentrations_test(pH_res, c_salt_res):
+
+ # Create an instance of the pyMBE library
+ pmb = pyMBE.pymbe_library(seed=42)
+
+ # Determine the reservoir composition using pyMBE
+ cH_res_pyMBE, cOH_res_pyMBE, cNa_res_pyMBE, cCl_res_pyMBE = pmb.determine_reservoir_concentrations(
+ pH_res,
+ c_salt_res * pmb.units.mol/ pmb.units.L,
+ lambda x: 1.0)
+
+ # Determine the reservoir concentrations without pyMBE
+ cH_res = 10**(-pH_res) * pmb.units.mol/ pmb.units.L
+ cOH_res = 10**(pH_res-14) * pmb.units.mol/ pmb.units.L
+ c_salt_res = c_salt_res * pmb.units.mol/ pmb.units.L
+ cNa_res = max(c_salt_res, c_salt_res + cOH_res - cH_res)
+ cCl_res = max(c_salt_res, c_salt_res + cH_res - cOH_res)
+
+ # Check the determine equilibrium constants
+ np.testing.assert_allclose(cH_res, cH_res_pyMBE)
+ np.testing.assert_allclose(cOH_res, cOH_res_pyMBE)
+ np.testing.assert_allclose(cNa_res, cNa_res_pyMBE)
+ np.testing.assert_allclose(cCl_res, cCl_res_pyMBE)
+
+print("*** Unit test: check that pyMBE correctly calculates the reservoir composition when setting up the grand-reaction method. ***")
+for pH_res in np.linspace(1.0, 13.0, 5):
+ for c_salt_res in np.logspace(-5, 0, 5):
+ determine_reservoir_concentrations_test(pH_res, c_salt_res)
+print("*** Unit test passed ***")
diff --git a/testsuite/reaction_methods_unit_tests.py b/testsuite/reaction_methods_unit_tests.py
new file mode 100644
index 0000000..824a333
--- /dev/null
+++ b/testsuite/reaction_methods_unit_tests.py
@@ -0,0 +1,362 @@
+
+# Copyright (C) 2024 pyMBE-dev team
+#
+# This file is part of pyMBE.
+#
+# pyMBE 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.
+#
+# pyMBE 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 .
+
+# Import pyMBE and other libraries
+import pyMBE
+import numpy as np
+import espressomd
+
+def reaction_method_test_template(parameters):
+
+ # Create an instance of the pyMBE library
+ pmb = pyMBE.pymbe_library(seed=42)
+
+ if parameters["method"] in ["cpH", "grxmc", "grxmc_unified"]:
+ # Define the acidic particle
+ pmb.define_particle(
+ name = "A",
+ acidity = "acidic",
+ pka = parameters["pK_acid"],
+ sigma = 1*pmb.units('reduced_length'),
+ epsilon = 1*pmb.units('reduced_energy'))
+
+ # Define the basic particle
+ pmb.define_particle(
+ name = "B",
+ acidity = "basic",
+ pka = parameters["pK_base"],
+ sigma = 1*pmb.units('reduced_length'),
+ epsilon = 1*pmb.units('reduced_energy'))
+
+
+ # Define the ions
+ pmb.define_particle(
+ name="Na",
+ z=parameters["z_Na"],
+ sigma = 1*pmb.units('reduced_length'),
+ epsilon = 1*pmb.units('reduced_energy'))
+
+ pmb.define_particle(
+ name="Cl",
+ z=parameters["z_Cl"],
+ sigma = 1*pmb.units('reduced_length'),
+ epsilon = 1*pmb.units('reduced_energy'))
+
+ pmb.define_particle(
+ name="H",
+ z=parameters["z_H"],
+ sigma = 1*pmb.units('reduced_length'),
+ epsilon = 1*pmb.units('reduced_energy'))
+
+ pmb.define_particle(
+ name="OH",
+ z=parameters["z_OH"],
+ sigma = 1*pmb.units('reduced_length'),
+ epsilon = 1*pmb.units('reduced_energy'))
+
+ if parameters["method"] == "cpH":
+ # Add the reactions using pyMBE
+ if "pka_set" in parameters:
+ cpH, _ = pmb.setup_cpH(counter_ion="H",
+ constant_pH=parameters["pH"],
+ use_exclusion_radius_per_type=parameters["use_exclusion_radius_per_type"],
+ pka_set=parameters["pka_set"])
+ else:
+ cpH, _ = pmb.setup_cpH(counter_ion="H",
+ constant_pH=parameters["pH"],
+ use_exclusion_radius_per_type=parameters["use_exclusion_radius_per_type"])
+
+
+ # Check the number of reactions
+ np.testing.assert_equal(len(cpH.reactions), 4)
+
+ # Check the equilibrium constants
+ np.testing.assert_allclose(cpH.reactions[0].gamma, 10**(-parameters["pK_acid"]))
+ np.testing.assert_allclose(cpH.reactions[1].gamma, 10**parameters["pK_acid"])
+
+ np.testing.assert_allclose(cpH.reactions[2].gamma, 10**(-parameters["pK_base"]))
+ np.testing.assert_allclose(cpH.reactions[3].gamma, 10**parameters["pK_base"])
+
+
+ elif parameters["method"] == "gcmc":
+ input_parameters = {
+ "c_salt_res": parameters["c_salt_res"] * pmb.units.mol/ pmb.units.L,
+ "salt_cation_name": "Na",
+ "salt_anion_name": "Cl",
+ "activity_coefficient": lambda x: 1.0,
+ "use_exclusion_radius_per_type": parameters["use_exclusion_radius_per_type"]}
+
+ # Check that pyMBE raises an error if wrong charge signs are provided
+ if parameters["z_Na"]<0:
+ np.testing.assert_raises(ValueError, pmb.setup_gcmc, **input_parameters)
+ return
+ if parameters["z_Cl"]>0:
+ np.testing.assert_raises(ValueError, pmb.setup_gcmc, **input_parameters)
+ return
+
+ # Add the reactions using pyMBE
+ gcmc = pmb.setup_gcmc(**input_parameters)
+
+ # Check the number of reactions
+ np.testing.assert_equal(len(gcmc.reactions), 2)
+
+ # Check the equilibrium constants
+ K_NaCl = (parameters["c_salt_res"] * pmb.units.mol/ pmb.units.L).to('1/(N_A * reduced_length**3)').magnitude**2
+ np.testing.assert_allclose(gcmc.reactions[0].gamma, K_NaCl)
+ np.testing.assert_allclose(gcmc.reactions[1].gamma, 1/K_NaCl)
+
+ elif parameters["method"] == "grxmc":
+ input_parameters = {
+ "pH_res": parameters["pH_res"],
+ "c_salt_res": parameters["c_salt_res"] * pmb.units.mol/ pmb.units.L,
+ "proton_name": "H",
+ "hydroxide_name": "OH",
+ "salt_cation_name": "Na",
+ "salt_anion_name": "Cl",
+ "activity_coefficient": lambda x: 1.0,
+ "use_exclusion_radius_per_type": parameters["use_exclusion_radius_per_type"]}
+
+ # Check that pyMBE raises an error if wrong charge signs are provided
+ if parameters["z_H"]<0:
+ np.testing.assert_raises(ValueError, pmb.setup_grxmc_reactions, **input_parameters)
+ return
+ if parameters["z_Na"]<0:
+ np.testing.assert_raises(ValueError, pmb.setup_grxmc_reactions, **input_parameters)
+ return
+ if parameters["z_OH"]>0:
+ np.testing.assert_raises(ValueError, pmb.setup_grxmc_reactions, **input_parameters)
+ return
+ if parameters["z_Cl"]>0:
+ np.testing.assert_raises(ValueError, pmb.setup_grxmc_reactions, **input_parameters)
+ return
+
+ if "pka_set" in parameters:
+ input_parameters["pka_set"] = parameters["pka_set"]
+ grxmc, *_ = pmb.setup_grxmc_reactions(**input_parameters)
+ else:
+ grxmc, *_ = pmb.setup_grxmc_reactions(**input_parameters)
+
+ # Check the number of reactions
+ np.testing.assert_equal(len(grxmc.reactions), 28)
+
+ # Determine the reservoir concentrations independent from pyMBE
+ cH_res = (10**(-parameters["pH_res"]) * pmb.units.mol/ pmb.units.L).to('1/(N_A * reduced_length**3)').magnitude
+ cOH_res = (10**(parameters["pH_res"]-14) * pmb.units.mol/ pmb.units.L).to('1/(N_A * reduced_length**3)').magnitude
+ c_salt_res = (parameters["c_salt_res"] * pmb.units.mol/ pmb.units.L).to('1/(N_A * reduced_length**3)').magnitude
+ cNa_res = max(c_salt_res, c_salt_res + cOH_res - cH_res)
+ cCl_res = max(c_salt_res, c_salt_res + cH_res - cOH_res)
+ Ka_acid = (10**(-parameters["pK_acid"]) * pmb.units.mol/ pmb.units.L).to('1/(N_A * reduced_length**3)').magnitude
+ Ka_base = (10**(-parameters["pK_base"]) * pmb.units.mol/ pmb.units.L).to('1/(N_A * reduced_length**3)').magnitude
+
+ # Check the equilibrium constants
+ np.testing.assert_allclose(grxmc.reactions[0].gamma, cH_res*cOH_res)
+ np.testing.assert_allclose(grxmc.reactions[1].gamma, 1/(cH_res*cOH_res))
+ np.testing.assert_allclose(grxmc.reactions[2].gamma, cNa_res*cCl_res)
+ np.testing.assert_allclose(grxmc.reactions[3].gamma, 1/(cNa_res*cCl_res))
+ np.testing.assert_allclose(grxmc.reactions[4].gamma, cNa_res*cOH_res)
+ np.testing.assert_allclose(grxmc.reactions[5].gamma, 1/(cNa_res*cOH_res))
+ np.testing.assert_allclose(grxmc.reactions[6].gamma, cH_res*cCl_res)
+ np.testing.assert_allclose(grxmc.reactions[7].gamma, 1/(cH_res*cCl_res))
+
+ np.testing.assert_allclose(grxmc.reactions[8].gamma, cNa_res/cH_res)
+ np.testing.assert_allclose(grxmc.reactions[9].gamma, cH_res/cNa_res)
+ np.testing.assert_allclose(grxmc.reactions[10].gamma, cCl_res/cOH_res)
+ np.testing.assert_allclose(grxmc.reactions[11].gamma, cOH_res/cCl_res)
+
+ np.testing.assert_allclose(grxmc.reactions[12].gamma, Ka_acid)
+ np.testing.assert_allclose(grxmc.reactions[13].gamma, 1/Ka_acid)
+ np.testing.assert_allclose(grxmc.reactions[14].gamma, Ka_acid*cNa_res/cH_res)
+ np.testing.assert_allclose(grxmc.reactions[15].gamma, cH_res/(cNa_res*Ka_acid))
+ np.testing.assert_allclose(grxmc.reactions[16].gamma, Ka_acid/(cH_res*cOH_res))
+ np.testing.assert_allclose(grxmc.reactions[17].gamma, cH_res*cOH_res/Ka_acid)
+ np.testing.assert_allclose(grxmc.reactions[18].gamma, Ka_acid/(cH_res*cCl_res))
+ np.testing.assert_allclose(grxmc.reactions[19].gamma, cH_res*cCl_res/Ka_acid)
+
+ np.testing.assert_allclose(grxmc.reactions[20].gamma, Ka_base)
+ np.testing.assert_allclose(grxmc.reactions[21].gamma, 1/Ka_base)
+ np.testing.assert_allclose(grxmc.reactions[22].gamma, Ka_base*cNa_res/cH_res)
+ np.testing.assert_allclose(grxmc.reactions[23].gamma, cH_res/(cNa_res*Ka_base))
+ np.testing.assert_allclose(grxmc.reactions[24].gamma, Ka_base/(cH_res*cOH_res))
+ np.testing.assert_allclose(grxmc.reactions[25].gamma, cH_res*cOH_res/Ka_base)
+ np.testing.assert_allclose(grxmc.reactions[26].gamma, Ka_base/(cH_res*cCl_res))
+ np.testing.assert_allclose(grxmc.reactions[27].gamma, cH_res*cCl_res/Ka_base)
+
+ elif parameters["method"] == "grxmc_unified":
+ input_parameters = {
+ "pH_res": parameters["pH_res"],
+ "c_salt_res": parameters["c_salt_res"] * pmb.units.mol/ pmb.units.L,
+ "cation_name": "H",
+ "anion_name": "OH",
+ "activity_coefficient": lambda x: 1.0,
+ "use_exclusion_radius_per_type": parameters["use_exclusion_radius_per_type"]}
+
+ # Check that pyMBE raises an error if wrong charge signs are provided
+ if parameters["z_H"]<0:
+ np.testing.assert_raises(ValueError, pmb.setup_grxmc_unified, **input_parameters)
+ return
+ if parameters["z_OH"]>0:
+ np.testing.assert_raises(ValueError, pmb.setup_grxmc_unified, **input_parameters)
+ return
+
+ if "pka_set" in parameters:
+ input_parameters["pka_set"] = parameters["pka_set"]
+ grxmc, *_ = pmb.setup_grxmc_unified(**input_parameters)
+ else:
+ grxmc, *_ = pmb.setup_grxmc_unified(**input_parameters)
+
+ # Check the number of reactions
+ np.testing.assert_equal(len(grxmc.reactions), 10)
+
+ # Determine the reservoir concentrations independent from pyMBE
+ cH_res = (10**(-parameters["pH_res"]) * pmb.units.mol/ pmb.units.L).to('1/(N_A * reduced_length**3)').magnitude
+ cOH_res = (10**(parameters["pH_res"]-14) * pmb.units.mol/ pmb.units.L).to('1/(N_A * reduced_length**3)').magnitude
+ c_salt_res = (parameters["c_salt_res"] * pmb.units.mol/ pmb.units.L).to('1/(N_A * reduced_length**3)').magnitude
+ cNa_res = max(c_salt_res, c_salt_res + cOH_res - cH_res)
+ cCl_res = max(c_salt_res, c_salt_res + cH_res - cOH_res)
+ c_cation_res = cH_res + cNa_res
+ c_anion_res = cOH_res + cCl_res
+ Ka_acid = (10**(-parameters["pK_acid"]) * pmb.units.mol/ pmb.units.L).to('1/(N_A * reduced_length**3)').magnitude
+ Ka_base = (10**(-parameters["pK_base"]) * pmb.units.mol/ pmb.units.L).to('1/(N_A * reduced_length**3)').magnitude
+
+ # Check the equilibrium constants
+ np.testing.assert_allclose(grxmc.reactions[0].gamma, c_cation_res*c_anion_res)
+ np.testing.assert_allclose(grxmc.reactions[1].gamma, 1/(c_cation_res*c_anion_res))
+
+ np.testing.assert_allclose(grxmc.reactions[2].gamma, Ka_acid*c_cation_res/cH_res)
+ np.testing.assert_allclose(grxmc.reactions[3].gamma, cH_res/(Ka_acid*c_cation_res))
+ np.testing.assert_allclose(grxmc.reactions[4].gamma, Ka_acid/(cH_res*c_anion_res))
+ np.testing.assert_allclose(grxmc.reactions[5].gamma, cH_res*c_anion_res/Ka_acid)
+
+ np.testing.assert_allclose(grxmc.reactions[6].gamma, Ka_base*c_cation_res/cH_res)
+ np.testing.assert_allclose(grxmc.reactions[7].gamma, cH_res/(Ka_base*c_cation_res))
+ np.testing.assert_allclose(grxmc.reactions[8].gamma, Ka_base/(cH_res*c_anion_res))
+ np.testing.assert_allclose(grxmc.reactions[9].gamma, cH_res*c_anion_res/Ka_base)
+
+
+# Set up the espresso system
+espresso_system=espressomd.System(box_l = [10.0]*3)
+
+# cpH test
+print("*** Unit test: check that reactions are correctly set up in the cpH method. ***")
+for use_exclusion_radius_per_type in [False, True]:
+ parameters = {
+ "method": "cpH",
+ "pK_acid": 4.0,
+ "pK_base": 8.0,
+ "pH": 7.0,
+ "z_Na": 1,
+ "z_Cl": -1,
+ "z_H": 1,
+ "z_OH": -1,
+ "use_exclusion_radius_per_type": use_exclusion_radius_per_type
+ }
+ reaction_method_test_template(parameters)
+
+ parameters["pka_set"] = {
+ "A": {"pka_value": 4.0, "acidity": "acidic"},
+ "B": {"pka_value": 8.0, "acidity": "basic"},
+ "C": {"pka_value": 7.0, "acidity": "acidi"}}
+ reaction_method_test_template(parameters)
+print("*** Unit test passed ***")
+
+# gcmc test
+print("*** Unit test: check that reactions are correctly set up in the GCMC method. ***")
+for use_exclusion_radius_per_type in [False, True]:
+ parameters = {
+ "method": "gcmc",
+ "c_salt_res": 1,
+ "z_Na": 1,
+ "z_Cl": -1,
+ "z_H": 1,
+ "z_OH": -1,
+ "use_exclusion_radius_per_type": use_exclusion_radius_per_type
+ }
+ reaction_method_test_template(parameters)
+
+ parameters["z_Cl"] = 1
+ reaction_method_test_template(parameters)
+
+ parameters["z_Na"] = -1
+ reaction_method_test_template(parameters)
+print("*** Unit test passed ***")
+
+# grxmc test
+print("*** Unit test: check that reactions are correctly set up in the G-RxMC method. ***")
+for use_exclusion_radius_per_type in [False, True]:
+ parameters = {
+ "method": "grxmc",
+ "pK_acid": 4.0,
+ "pK_base": 9.0,
+ "c_salt_res": 1,
+ "pH_res": 5.0,
+ "z_Na": 1,
+ "z_Cl": -1,
+ "z_H": 1,
+ "z_OH": -1,
+ "use_exclusion_radius_per_type": use_exclusion_radius_per_type
+ }
+ reaction_method_test_template(parameters)
+
+ parameters["pka_set"] = {
+ "A": {"pka_value": 4.0, "acidity": "acidic"},
+ "B": {"pka_value": 9.0, "acidity": "basic"},
+ "C": {"pka_value": 7.0, "acidity": "acidi"}}
+ reaction_method_test_template(parameters)
+
+ parameters["z_Cl"] = 1
+ reaction_method_test_template(parameters)
+
+ parameters["z_OH"] = 1
+ reaction_method_test_template(parameters)
+
+ parameters["z_Na"] = -1
+ reaction_method_test_template(parameters)
+
+ parameters["z_H"] = -1
+ reaction_method_test_template(parameters)
+print("*** Unit test passed ***")
+
+# grxmc unified test
+print("*** Unit test: check that reactions are correctly set up in the unified G-RxMC method. ***")
+for use_exclusion_radius_per_type in [False, True]:
+ parameters = {
+ "method": "grxmc_unified",
+ "pK_acid": 4.0,
+ "pK_base": 9.0,
+ "c_salt_res": 1,
+ "pH_res": 5.0,
+ "z_Na": 1,
+ "z_Cl": -1,
+ "z_H": 1,
+ "z_OH": -1,
+ "use_exclusion_radius_per_type": use_exclusion_radius_per_type
+ }
+ reaction_method_test_template(parameters)
+
+ parameters["pka_set"] = {
+ "A": {"pka_value": 4.0, "acidity": "acidic"},
+ "B": {"pka_value": 9.0, "acidity": "basic"},
+ "C": {"pka_value": 7.0, "acidity": "acidi"}}
+ reaction_method_test_template(parameters)
+
+ parameters["z_OH"] = 1
+ reaction_method_test_template(parameters)
+
+ parameters["z_H"] = -1
+ reaction_method_test_template(parameters)
+print("*** Unit test passed ***")