Solution: add automatic charge balancing and other misc improvements

.pre-commit-config.yaml

@@ -12,26 +12,26 @@ ci:
 
 repos:
   - repo: https://github.com/astral-sh/ruff-pre-commit
-    rev: v0.4.8
+    rev: v0.5.5
     hooks:
       - id: ruff
         args: [--fix, --ignore, "D,E501", "--show-fixes"]
 
   - repo: https://github.com/psf/black-pre-commit-mirror
-    rev: 24.2.0
+    rev: 24.4.2
     hooks:
       - id: black
 
   - repo: https://github.com/codespell-project/codespell
-    rev: v2.2.6
+    rev: v2.3.0
     hooks:
       - id: codespell
         stages: [commit, commit-msg]
         exclude_types: [html, svg]
         additional_dependencies: [tomli] # needed to read pyproject.toml below py3.11
 
   - repo: https://github.com/pre-commit/pre-commit-hooks
-    rev: v4.5.0
+    rev: v4.6.0
     hooks:
       - id: check-case-conflict
       - id: check-symlinks

CHANGELOG.md

@@ -5,16 +5,29 @@ All notable changes to this project will be documented in this file.
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
-## [1.0.4] - 2024-07-26
+## [1.1.0] - 2024-07-27
 
 ### Fixed
 
 - `equilibrate`: Fixed a bug that could cause an `AttributeError` when pH was used for charge
   balancing.
 - Database: `size.radius_ionic` was missing units for `Ni[+2]` and `Cr[+3]`. Correct units have been added.
 
+### Added
+
+- `Solution`: New automatic charge balancing method will automatically identify the majority (highest concentration)
+  cation or anion as appropriate (depending on the charge balance) for charge balancing. To use this mode, set
+  `balance_charge='auto'` when instantiating a `Solution`.
+
 ### Changed
 
+- `Solution.add_amount`: This method will now add solutes that are absent from the Solution. Previously, calling, e.g.,
+  `add_amount('Na+', '1 mol')` on a `Solution` that did not contain any sodium would result in an error. A warning
+  is logged if the method has to add a new solute.
+- Units: use the upstream chemistry context from `pint` instead of the custom one from 2013.
+- `pre-commit autoupdate`
+- Misc. linting and code quality improvements.
+- Unit tests: update `tmpdir` to `tmp_path` text fixture.
 - CI: Small updates to pre-commit and GitHub actions per scientific python [repo review](https://scientific-python.github.io/repo-review/?repo=kingsburylab%2FpyEQL&branch=main).
 
 ## [1.0.3] - 2024-07-20

MANIFEST.in

@@ -2,6 +2,5 @@ include README COPYING CHANGELOG LICENSE README.md README.rst README.txt CHANGES
 recursive-include src/pyEQL/database/ *
 recursive-include src/pyEQL/presets/ *
 recursive-include docs/ *
-include src/pyEQL/pint_custom_units.txt
 prune docs/build
 global-exclude *.pyc *~ .DS_Store *__pycache__* *.pyo

pyproject.toml

@@ -2,7 +2,7 @@
 name = "pyEQL"
 readme = "README.md"
 dynamic = ["version"]
-description="A python interace for solution chemistry"
+description="A python interface for solution chemistry"
 authors =[
     {name = "Ryan Kingsbury", email = "[email protected]"}
 ]
@@ -25,7 +25,7 @@ dependencies = [
         "scipy",
         "pymatgen==2024.5.1",
         "iapws",
-        "monty",
+        "monty>=2024.7.12",
         "maggma>=0.67.0",
         "phreeqpython",
     ]

src/pyEQL/__init__.py

@@ -37,11 +37,8 @@
 # convert "offset units" so that, e.g. Quantity('25 degC') works without error
 # see https://pint.readthedocs.io/en/0.22/user/nonmult.html?highlight=offset#temperature-conversion
 ureg.autoconvert_offset_to_baseunit = True
-# append custom unit definitions and contexts
-fname = files("pyEQL") / "pint_custom_units.txt"
-ureg.load_definitions(fname)
 # activate the "chemistry" context globally
-ureg.enable_contexts("chem")
+ureg.enable_contexts("chemistry")
 # set the default string formatting for pint quantities
 ureg.default_format = "P~"
 

src/pyEQL/solution.py

@@ -92,11 +92,15 @@ def __init__(
                 -7 to +14. The default value corresponds to a pE value typical of natural
                 waters in equilibrium with the atmosphere.
             balance_charge: The strategy for balancing charge during init and equilibrium calculations. Valid options
-                are 'pH', which will adjust the solution pH to balance charge, 'pE' which will adjust the
-                redox equilibrium to balance charge, or the name of a dissolved species e.g. 'Ca+2' or 'Cl-'
-                that will be added/subtracted to balance charge. If set to None, no charge balancing will be
-                performed either on init or when equilibrate() is called. Note that in this case, equilibrate()
-                can distort the charge balance!
+                are
+                    - 'pH', which will adjust the solution pH to balance charge,
+                    - 'auto' which will use the majority cation or anion (i.e., that with the largest concentration)
+                    as needed,
+                    - 'pE' (not currently implemented) which will adjust the redox equilibrium to balance charge, or
+                    the name of a dissolved species e.g. 'Ca+2' or 'Cl-' that will be added/subtracted to balance
+                    charge.
+                    - None (default), in which case no charge balancing will be performed either on init or when
+                    equilibrate() is called. Note that in this case, equilibrate() can distort the charge balance!
             solvent: Formula of the solvent. Solvents other than water are not supported at this time.
             engine: Electrolyte modeling engine to use. See documentation for details on the available engines.
             database: path to a .json file (str or Path) or maggma Store instance that
@@ -171,7 +175,7 @@ def __init__(
         self._pE = pE
         self._pH = pH
         self.pE = self._pE
-        if isinstance(balance_charge, str) and balance_charge not in ["pH", "pE"]:
+        if isinstance(balance_charge, str) and balance_charge not in ["pH", "pE", "auto"]:
             self.balance_charge = standardize_formula(balance_charge)
         else:
             self.balance_charge = balance_charge  #: Standardized formula of the species used for charge balancing.
@@ -273,13 +277,19 @@ def __init__(
                 raise NotImplementedError("Balancing charge via redox (pE) is not yet implemented!")
             else:
                 ions = set().union(*[self.cations, self.anions])  # all ions
+                if self.balance_charge == "auto":
+                    # add the most abundant ion of the opposite charge
+                    if cb <= 0:
+                        self.balance_charge = max(self.cations, key=self.cations.get)
+                    elif cb > 0:
+                        self.balance_charge = max(self.anions, key=self.anions.get)
                 if self.balance_charge not in ions:
                     raise ValueError(
                         f"Charge balancing species {self.balance_charge} was not found in the solution!. "
                         f"Species {ions} were found."
                     )
-                z = self.get_property(balance_charge, "charge")
-                self.components[balance_charge] += -1 * cb / z * self.volume.to("L").magnitude
+                z = self.get_property(self.balance_charge, "charge")
+                self.components[self.balance_charge] += -1 * cb / z * self.volume.to("L").magnitude
                 balanced = True
 
             if not balanced:
@@ -1282,81 +1292,12 @@ def add_amount(self, solute: str, amount: str):
         Returns:
             Nothing. The concentration of solute is modified.
         """
-        # if units are given on a per-volume basis,
-        # iteratively solve for the amount of solute that will preserve the
-        # original volume and result in the desired concentration
-        if ureg.Quantity(amount).dimensionality in (
-            "[substance]/[length]**3",
-            "[mass]/[length]**3",
-        ):
-            # store the original volume for later
-            orig_volume = self.volume
-
-            # change the amount of the solute present to match the desired amount
-            self.components[solute] += (
-                ureg.Quantity(amount)
-                .to(
-                    "moles",
-                    "chem",
-                    mw=self.get_property(solute, "molecular_weight"),
-                    volume=self.volume,
-                    solvent_mass=self.solvent_mass,
-                )
-                .magnitude
-            )
-
-            # set the amount to zero and log a warning if the desired amount
-            # change would result in a negative concentration
-            if self.get_amount(solute, "mol").magnitude < 0:
-                self.logger.error(
-                    "Attempted to set a negative concentration for solute %s. Concentration set to 0" % solute
-                )
-                self.set_amount(solute, "0 mol")
-
-            # calculate the volume occupied by all the solutes
-            solute_vol = self._get_solute_volume()
-
-            # determine the volume of solvent that will preserve the original volume
-            target_vol = orig_volume - solute_vol
-
-            # adjust the amount of solvent
-            # volume in L, density in kg/m3 = g/L
-            target_mass = target_vol * ureg.Quantity(self.water_substance.rho, "g/L")
-
-            mw = self.get_property(self.solvent, "molecular_weight")
-            target_mol = target_mass / mw
-            self.components[self.solvent] = target_mol.magnitude
-
-        else:
-            # change the amount of the solute present
-            self.components[solute] += (
-                ureg.Quantity(amount)
-                .to(
-                    "moles",
-                    "chem",
-                    mw=self.get_property(solute, "molecular_weight"),
-                    volume=self.volume,
-                    solvent_mass=self.solvent_mass,
-                )
-                .magnitude
-            )
-
-            # set the amount to zero and log a warning if the desired amount
-            # change would result in a negative concentration
-            if self.get_amount(solute, "mol").magnitude < 0:
-                self.logger.error(
-                    "Attempted to set a negative concentration for solute %s. Concentration set to 0" % solute
-                )
-                self.set_amount(solute, "0 mol")
-
-            # update the volume to account for the space occupied by all the solutes
-            # make sure that there is still solvent present in the first place
-            if self.solvent_mass <= ureg.Quantity(0, "kg"):
-                self.logger.error("All solvent has been depleted from the solution")
-                return
-
-            # set the volume recalculation flag
-            self.volume_update_required = True
+        # Get the current amount of the solute
+        current_amt = self.get_amount(solute, amount.split(" ")[1])
+        if current_amt.magnitude == 0:
+            self.logger.warning(f"Add new solute {solute} to the solution")
+        new_amt = ureg.Quantity(amount) + current_amt
+        self.set_amount(solute, new_amt)
 
     def set_amount(self, solute: str, amount: str):
         """
@@ -2465,7 +2406,7 @@ def to_file(self, filename: str | Path) -> None:
         """
         str_filename = str(filename)
         if not ("yaml" in str_filename.lower() or "json" in str_filename.lower()):
-            self.logger.error("Invalid file extension entered - %s" % str_filename)
+            self.logger.error("Invalid file extension entered - {str_filename}")
             raise ValueError("File extension must be .json or .yaml")
         if "yaml" in str_filename.lower():
             solution_dict = self.as_dict()

tests/test_solution.py

@@ -9,7 +9,6 @@
 import copy
 import os
 import platform
-from pathlib import Path
 
 import numpy as np
 import pytest
@@ -255,6 +254,31 @@ def test_charge_balance(s3, s5, s5_pH, s6, s6_Ca):
     assert np.isclose(s6.charge_balance, -0.12)
     assert np.isclose(s6_Ca.charge_balance, 0, atol=1e-8)
 
+    # test auto charge balance
+    s = 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",
+        balance_charge="auto",
+    )
+    assert s.balance_charge == "Na[+1]"
+    assert np.isclose(s.charge_balance, 0, atol=1e-8)
+    s.equilibrate()
+    assert s.balance_charge == "Na[+1]"
+
+    s = Solution({"Na+": "2 mM", "Cl-": "1 mM"}, balance_charge="auto")
+    assert s.balance_charge == "Cl[-1]"
+    assert np.isclose(s.charge_balance, 0, atol=1e-8)
+    s.equilibrate()
+    assert s.balance_charge == "Cl[-1]"
+
 
 def test_alkalinity_hardness(s3, s5, s6):
     assert np.isclose(s3.hardness, 0)
@@ -621,11 +645,11 @@ def test_as_from_dict(s1, s2):
     # assert s2_new.database != s2.database
 
 
-def test_serialization(s1, s2, tmpdir):
+def test_serialization(s1, s2, tmp_path):
     from monty.serialization import dumpfn, loadfn
 
-    dumpfn(s1, str(tmpdir / "s1.json"))
-    s1_new = loadfn(str(tmpdir / "s1.json"))
+    dumpfn(s1, str(tmp_path / "s1.json"))
+    s1_new = loadfn(str(tmp_path / "s1.json"))
     assert s1_new.volume.magnitude == 2
     assert s1_new._solutes["H[+1]"] == "2e-07 mol"
     assert s1_new.get_total_moles_solute() == s1.get_total_moles_solute()
@@ -644,8 +668,8 @@ def test_serialization(s1, s2, tmpdir):
     # TODO currently this test will fail due to a bug in maggma's __eq__
     # assert s1_new.database != s1.database
 
-    dumpfn(s2, str(tmpdir / "s2.json"))
-    s2_new = loadfn(str(tmpdir / "s2.json"))
+    dumpfn(s2, str(tmp_path / "s2.json"))
+    s2_new = loadfn(str(tmp_path / "s2.json"))
     assert s2_new.volume == s2.volume
     # components concentrations should be the same
     assert s2_new.components == s2.components
@@ -667,7 +691,7 @@ def test_serialization(s1, s2, tmpdir):
     # assert s2_new.database != s2.database
 
 
-def test_from_preset(tmpdir):
+def test_from_preset(tmp_path):
     from monty.serialization import dumpfn
 
     preset_name = "seawater"
@@ -685,7 +709,6 @@ def test_from_preset(tmpdir):
     with pytest.raises(FileNotFoundError):
         Solution.from_preset("nonexistent_preset")
     # test json as preset
-    tmp_path = Path(tmpdir)
     json_preset = tmp_path / "test.json"
     dumpfn(solution, json_preset)
     solution_json = Solution.from_preset(tmp_path / "test")
@@ -695,8 +718,7 @@ def test_from_preset(tmpdir):
     assert np.isclose(solution_json.pH, data["pH"], atol=0.01)
 
 
-def test_test_to_from_file(tmpdir, s1):
-    tmp_path = Path(tmpdir)
+def test_to_from_file(tmp_path, s1):
     for f in ["test.json", "test.yaml"]:
         filename = tmp_path / f
         s1.to_file(filename)