Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Add unit tests for reaction methods #85

Merged
merged 3 commits into from
Aug 20, 2024
Merged

Add unit tests for reaction methods #85

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
34ac21f
Add unit tests for reaction methods
davidbbeyer Aug 19, 2024
be48338
Fixed formatting error
davidbbeyer Aug 19, 2024
9328333
Full coverage of reaction method setup
davidbbeyer Aug 20, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 2 additions & 0 deletions Makefile
View file
Open in desktop
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
Open in desktop
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
51 changes: 51 additions & 0 deletions testsuite/determine_reservoir_concentrations_unit_test.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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 ***")
Loading
Loading