diff --git a/CHANGELOG.md b/CHANGELOG.md index 602b8d81..5ac7d2c2 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -9,6 +9,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Added +- `Solution`: add tests for `charge_balance`, `alkalinity`, `hardness` - `Solution`: add support for passing solutes as a `dict` - Implement extensible system for connecting `Solution` to various activity and speciation models. Models can be integrated into pyEQL by implementing an `EOS` class. The desired @@ -22,6 +23,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Changed +- `Solution.charge_balance` now returns in equivalents instead of Coulombs - Replace `water_properties.py` with [iapws](https://github.com/jjgomera/iapws) package - Replace `elements.py`` with `pymatgen.core.periodic_table` - Migrate all tests to `pytest` diff --git a/src/pyEQL/solution.py b/src/pyEQL/solution.py index 5ed6f943..4d06488a 100644 --- a/src/pyEQL/solution.py +++ b/src/pyEQL/solution.py @@ -25,6 +25,7 @@ from pyEQL.logging_system import logger from pyEQL.salt_ion_match import generate_salt_list, identify_salt +EQUIV_WT_CACO3 = 100.09 / 2 * unit.Quantity("g/mol") class Solution(MSONable): """ @@ -650,24 +651,23 @@ def charge_balance(self) -> float: on the solution and SHOULD equal zero at all times, but due to numerical errors will usually have a small nonzero value. It is calculated according to: - .. math:: CB = F \\sum_i n_i z_i + .. math:: CB = \\sum_i n_i z_i - where :math:`n_i` is the number of moles, :math:`z_i` is the charge on species i, and :math:`F` is the Faraday constant. + where :math:`n_i` is the number of moles, and :math:`z_i` is the charge on species i. Returns ------- float : - The charge balance of the solution, in equivalents. + The charge balance of the solution, in equivalents (mol of charge). """ charge_balance = 0 - F = (unit.e * unit.N_A).magnitude for solute in self.components: - charge_balance += self.get_amount(solute, "mol").magnitude * self.get_property(solute, "charge") * F + charge_balance += self.get_amount(solute, "mol").magnitude * self.get_property(solute, "charge") - return charge_balance + return charge_balance.magnitude - # TODO - need tests for alkalinity + # TODO - consider adding guard statements to prevent alkalinity from being negative @property def alkalinity(self): """ @@ -682,10 +682,10 @@ def alkalinity(self): ----- The alkalinity is calculated according to [stm]_ - .. math:: Alk = F \\sum_{i} z_{i} C_{B} - \\sum_{i} z_{i} C_{A} + .. math:: Alk = \\sum_{i} z_{i} C_{B} + \\sum_{i} z_{i} C_{A} Where :math:`C_{B}` and :math:`C_{A}` are conservative cations and anions, respectively - (i.e. ions that do not participate in acid-base reactions), and :math:`z_{i}` is their charge. + (i.e. ions that do not participate in acid-base reactions), and :math:`z_{i}` is their signed charge. In this method, the set of conservative cations is all Group I and Group II cations, and the conservative anions are all the anions of strong acids. @@ -696,34 +696,32 @@ def alkalinity(self): """ alkalinity = 0 * unit.Quantity("mol/L") - equiv_wt_CaCO3 = 100.09 / 2 * unit.Quantity("g/mol") - - base_cations = [ - "Li+", - "Na+", - "K+", - "Rb+", - "Cs+", - "Fr+", - "Be+2", - "Mg+2", - "Ca+2", - "Sr+2", - "Ba+2", - "Ra+2", - ] - acid_anions = ["Cl-", "Br-", "I-", "SO4-2", "NO3-", "ClO4-", "ClO3-"] + + base_cations = { + "Li[+1]", + "Na[+1]", + "K[+1]", + "Rb[+1]", + "Cs[+1]", + "Fr[+1]", + "Be[+2]", + "Mg[+2]", + "Ca[+2]", + "Sr[+2]", + "Ba[+2]", + "Ra[+2]", + } + acid_anions = {"Cl[-1]", "Br[-1]", "I[-1]", "SO4[-2]", "NO3[-1]", "ClO4[-1]", "ClO3[-1]"} for item in self.components: - if item in base_cations: + # sanitize the formulas + rform = Ion.from_formula(item).reduced_formula + if rform in base_cations.union(acid_anions): z = self.get_property(item, "charge") alkalinity += self.get_amount(item, "mol/L") * z - if item in acid_anions: - z = self.get_property(item, "charge") - alkalinity -= self.get_amount(item, "mol/L") * z # convert the alkalinity to mg/L as CaCO3 - return (alkalinity * equiv_wt_CaCO3).to("mg/L") + return (alkalinity * EQUIV_WT_CACO3).to("mg/L") @property def hardness(self): @@ -747,7 +745,7 @@ def hardness(self): """ hardness = 0 * unit.Quantity("mol/L") - equiv_wt_CaCO3 = 100.09 / 2 * unit.Quantity("g/mol") + for item in self.components: z = self.get_property(item, "charge") @@ -755,7 +753,7 @@ def hardness(self): hardness += z * self.get_amount(item, "mol/L") # convert the hardness to mg/L as CaCO3 - return (hardness * equiv_wt_CaCO3).to("mg/L") + return (hardness * EQUIV_WT_CACO3).to("mg/L") @property def debye_length(self) -> Quantity: diff --git a/tests/test_solution.py b/tests/test_solution.py index 7ac8b82d..f0bd1434 100644 --- a/tests/test_solution.py +++ b/tests/test_solution.py @@ -30,6 +30,24 @@ def s3(): def s4(): return Solution([["Na+", "8 mol"], ["Cl-", "8 mol"]], volume="2 L") +@pytest.fixture() +def s5(): + # 100 mg/L as CaCO3 + return Solution([["Ca+2", "40 mg/L"], ["CO3-2", "60 mg/L"]], volume="1 L") + +@pytest.fixture() +def s6(): + # non-electroneutral solution with lots of hardness + # alk = -118 meq/L * 50 = -5900 mg/L, hardness = 12*50 = 600 mg/L as CaCO3 + # charge balance = 2+10+10+10-120-20-12 = -120 meq/L + return Solution([["Ca+2", "1 mM"], # 2 meq/L + ["Mg+2", "5 mM"], # 10 meq/L + ["Na+1", "10 mM"], # 10 meq/L + ["Ag+1", "10 mM"], # no contribution to alk or hardness + ["CO3-2", "6 mM"], # no contribution to alk or hardness + ["SO4-2", "60 mM"], # -120 meq/L + ["Br-", "20 mM"]], # -20 meq/L + volume="1 L") def test_empty_solution_3(): # create an empty solution @@ -85,6 +103,18 @@ def test_solute_addition(s2, s3, s4): result_mol = s4.solvent_mass.to("kg").magnitude assert result_molL < result_mol +def test_alkalinity_hardness_chargebalance(s3, s5, s6): + assert np.isclose(s3.charge_balance, 0) + assert np.isclose(s3.hardness, 0) + assert np.isclose(s3.alkalinity, 0) + + assert np.isclose(s5.alkalinity.magnitude, 100, rtol=0.005) + assert np.isclose(s5.hardness.magnitude, 100, rtol=0.005) + assert np.isclose(s5.charge_balance, 0, atol=1e-5) + + assert np.isclose(s6.alkalinity.magnitude, -5900, rtol=0.005) + assert np.isclose(s6.hardness.magnitude, 600, rtol=0.005) + assert np.isclose(s6.charge_balance, -0.12) def test_serialization(s1, s2): assert isinstance(s1.as_dict(), dict)