Skip to content

Commit

Permalink
Full coverage of reaction method setup
Browse files Browse the repository at this point in the history
  • Loading branch information
@davidbbeyer
davidbbeyer committed Aug 20, 2024
1 parent be48338 commit 9328333
Show file tree
Hide file tree
Showing 3 changed files with 191 additions and 61 deletions.
2 changes: 2 additions & 0 deletions Makefile
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
20 changes: 13 additions & 7 deletions pyMBE.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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}.
Expand All @@ -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):
Expand Down Expand Up @@ -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
Expand All @@ -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]
Expand Down Expand Up @@ -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
)
Expand Down Expand Up @@ -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
)
Expand Down Expand Up @@ -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
)
Expand Down Expand Up @@ -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
Expand Down
230 changes: 176 additions & 54 deletions testsuite/reaction_methods_unit_tests.py
View file
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
#

# Copyright (C) 2024 pyMBE-dev team
#
# This file is part of pyMBE.
Expand Down Expand Up @@ -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)
Expand All @@ -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)
Expand All @@ -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)
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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 ***")

0 comments on commit 9328333

Please sign in to comment.