Solution: tests+tweaks for alk, hardness, charge bal

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`

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,15 +745,15 @@ 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")
             if z > 1:
                 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:

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)