Add unit tests for reaction methods (#85)

* Add unit tests for reaction methods

* Fixed formatting error

* Full coverage of reaction method setup

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

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

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 <http://www.gnu.org/licenses/>.
+
+# 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 ***")