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/reaction_methods_unit_tests.py b/testsuite/reaction_methods_unit_tests.py index cd8a562..824a333 100644 --- a/testsuite/reaction_methods_unit_tests.py +++ b/testsuite/reaction_methods_unit_tests.py @@ -1,4 +1,4 @@ -# + # Copyright (C) 2024 pyMBE-dev team # # This file is part of pyMBE. @@ -47,31 +47,40 @@ def reaction_method_test_template(parameters): # Define the ions pmb.define_particle( name="Na", - z=1, + z=parameters["z_Na"], sigma = 1*pmb.units('reduced_length'), epsilon = 1*pmb.units('reduced_energy')) - + pmb.define_particle( name="Cl", - z=-1, + z=parameters["z_Cl"], sigma = 1*pmb.units('reduced_length'), epsilon = 1*pmb.units('reduced_energy')) pmb.define_particle( name="H", - z=1, + z=parameters["z_H"], sigma = 1*pmb.units('reduced_length'), epsilon = 1*pmb.units('reduced_energy')) pmb.define_particle( name="OH", - z=-1, + 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 - cpH, _ = pmb.setup_cpH(counter_ion="H", constant_pH=parameters["pH"]) + 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) @@ -85,12 +94,23 @@ def reaction_method_test_template(parameters): 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( - 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) + gcmc = pmb.setup_gcmc(**input_parameters) # Check the number of reactions np.testing.assert_equal(len(gcmc.reactions), 2) @@ -101,14 +121,35 @@ def reaction_method_test_template(parameters): np.testing.assert_allclose(gcmc.reactions[1].gamma, 1/K_NaCl) elif parameters["method"] == "grxmc": - grxmc, *_ = pmb.setup_grxmc_reactions( - 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) + 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) @@ -156,12 +197,27 @@ def reaction_method_test_template(parameters): np.testing.assert_allclose(grxmc.reactions[27].gamma, cH_res*cCl_res/Ka_base) elif parameters["method"] == "grxmc_unified": - grxmc, *_ = pmb.setup_grxmc_unified( - 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) + 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) @@ -197,44 +253,110 @@ def reaction_method_test_template(parameters): # cpH test print("*** Unit test: check that reactions are correctly set up in the cpH method. ***") -parameters = { - "method": "cpH", - "pK_acid": 4.0, - "pK_base": 8.0, - "pH": 7.0 - } -reaction_method_test_template(parameters) +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. ***") -parameters = { - "method": "gcmc", - "c_salt_res": 1, - } -reaction_method_test_template(parameters) +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. ***") -parameters = { - "method": "grxmc", - "pK_acid": 4.0, - "pK_base": 9.0, - "c_salt_res": 1, - "pH_res": 5.0 - } -reaction_method_test_template(parameters) +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. ***") -parameters = { - "method": "grxmc_unified", - "pK_acid": 4.0, - "pK_base": 9.0, - "c_salt_res": 1, - "pH_res": 5.0 - } -reaction_method_test_template(parameters) +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 ***")