Merge pull request #34 from rkingsbury/sanitize_formula

Solution: sanitize all formulas

Author: rkingsbury

CHANGELOG.md

@@ -5,16 +5,24 @@ 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).
 
-## [Unreleased]
+## [0.7.0] - 2023-08-22
+
+### Changed
+
+- `Solution` now more robustly converts any user-supplied formulas into unique values using `pymatgen.core.ion.Ion.reduced_formula`. This means that the `.components` or `solvent` attributes may now differ slightly from whatever is entered during `__init__`. For example, `Solution(solvent='H2O').solvent` gives `H2O(aq)`. This behavior resolved a small bug that could occur when mixing solutions. User supplied formulas passed to `get_amount` or `Solution.components[xxx]` can still be any valid formula. E.g., `Solution.components["Na+"]`, `Solution.components["Na+1"]`, and `Solution.components["Na[+]"]` will all return the same thing.
+
+## [0.6.1] - 2023-08-22
 
 ### Added
 
+- `Solution`: enable passing an `EOS` instance to the `engine` kwarg.
 - `Solution`: new properties `total_dissolved_solids` and alias `TDS`
 - `Solution`: support new units in `get_amount` - ppm, ppb, eq/L, etc.
 - `Solution`: implemented arithmetic operations `+` (for mixing two solutions), `*` and `\` for scaling their amounts
 
 ### Changed
 
+
 - `pyEQL.unit` was renamed to `pyEQL.ureg` (short for `UnitRegistry`) for consistency with the `pint` documentation and tutorials.
 
 ## [v0.6.0] - 2023-08-15

docs/chemistry.md

@@ -10,47 +10,58 @@ calculations.
 ## How to Enter Valid Chemical Formulas
 
 Generally speaking, type the chemical formula of your solute the "normal" way
-and `pyEQL` should be able to inerpret it. Internally, `pyEQL` uses the [`pymatgen.core.ion.Ion`](https://pymatgen.org/pymatgen.core.html#pymatgen.core.ion.Ion)
- class to "translate" chemical formulas into a consistent format. Anything that the `Ion` class can understand will
- be processed into a valid formula by `pyEQL`.
+and `pyEQL` should be able to inerpret it. Internally, `pyEQL` uses a utility function `pyEQL.utils.standardize_formula`
+to process all formulas into a standard form. At present, this is done by passing the formula through the
+[`pymatgen.core.ion.Ion`](https://pymatgen.org/pymatgen.core.html#pymatgen.core.ion.Ion) class. Anything that the `Ion`
+class can understand will be processed into a valid formula by `pyEQL`.
 
 Here are some examples:
 
 | Substance | You enter | `pyEQL` understands |
 | :--- | :---: | :---: |
-| Sodium Chloride | "NaCl", "NaCl(aq)", or "ClNa" | "NaCl(aq)" | 
-| Sodium Sulfate | "Na(SO4)2" or "NaS2O8" | "Na(SO4)2(aq)" | 
-| Sodium Ion | "Na+", "Na+1", "Na1+", or "Na[+]" | "Na[+1]" | 
-| Magnesium Ion | "Mg+2", "Mg++", or "Mg[++]" | "Mg[+2]" | 
-| Methanol | "CH3OH", "CH4O" | "'CH3OH(aq)'" | 
+| Sodium Chloride | "NaCl", "NaCl(aq)", or "ClNa" | "NaCl(aq)" |
+| Sodium Sulfate | "Na(SO4)2" or "NaS2O8" | "Na(SO4)2(aq)" |
+| Sodium Ion | "Na+", "Na+1", "Na1+", or "Na[+]" | "Na[+1]" |
+| Magnesium Ion | "Mg+2", "Mg++", or "Mg[++]" | "Mg[+2]" |
+| Methanol | "CH3OH", "CH4O" | "'CH3OH(aq)'" |
 
-Specifically, `pyEQL` uses `Ion.from_formula(<formula>).reduced_formla` (shown in the right hand column of the table) to
-identify solutes. Notice that for charged species, the charges are always placed inside square brackets (e.g., `Na[+1]`)
-and always include the charge number (even for monovalent ions). Uncharged species are always suffixed by `(aq)` to
-disambiguate them from solids.
+Specifically, `standardize_formula` uses `Ion.from_formula(<formula>).reduced_formla` (shown in the right hand column
+of the table) to identify solutes. Notice that for charged species, the charges are always placed inside square brackets
+(e.g., `Na[+1]`) and always include the charge number (even for monovalent ions). Uncharged species are always suffixed
+by `(aq)` to disambiguate them from solids.
 
 :::{important}
-**When writing multivalent ion formulas, it is strongly recommended that you put the charge number AFTER the + or - 
-sign** (e.g., type "Mg+2" NOT "Mg2+"). The latter formula is ambiguous - it could mean $`Mg_2^+`$ or $`Mg^{+2}`$
+**When writing multivalent ion formulas, it is strongly recommended that you put the charge number AFTER the + or -
+sign** (e.g., type "Mg+2" NOT "Mg2+"). The latter formula is ambiguous - it could mean $`Mg_2^+`$ or $`Mg^{+2}`$ and
+it will be processed incorrectly into `Mg[+0.5]`
 :::
 
 (manual-testing)=
 ## Manually testing a formula
 
-If you want to make sure `pyEQL` is understanding your formula correctly, you can manually test it via `pymatgen` as
-follows:
+If you want to make sure `pyEQL` is understanding your formula correctly, you can manually test it as follows:
 
 ```
->>> from pymatgen.core.ion import Ion
->>> Ion.from_formula(<your_formula>).reduced_formula
+>>> from pyEQL.utils import standardize_formula
+>>> standardize_formula(<your_formula>)
 ...
 ```
 
 ## Formulas you will see when using `Solution`
 
-When using the `Solution` class, 
+When using the `Solution` class,
 
-- When creating a `Solution`, you can enter chemical formulas in any format you prefer, as long as `pymatgen` can understand it (see [manual testing](#manually-testing-a-formula)).
-- The keys (solute formulas) in `Solution.components` are preserved in the same format the user enters them. So if you entered `Na+` for sodium ion, it will stay that way.
-- Arguments to `Solution.get_property` can be entered in any format you prefer. When `pyEQL` queries the database, it will automatically convert the formula to the canonical one from `pymatgen`
-- Property data in the database is uniquely identified by the canonical ion formula (output of `Ion.from_formula(<formula>).reduced_formla`, e.g. "Na[+1]" for sodium ion).
+- When creating a `Solution`, you can enter chemical formulas in any format you prefer, as long as `standardize_formula` can understand it (see [manual testing](#manually-testing-a-formula)).
+- The keys (solute formulas) in `Solution.components` are standardized. So if you entered `Na+` for sodium ion, it will appear in `components` as `Na[+1]`.
+- However, the `components` attribute is a special dictionary that automatically standardizes formulas when accessed. So, you can still enter the formula
+  however you want. For example, the following all access or modify the same element in `components`:
+
+  ```
+  Solution.components.get('Na+')
+  Solution.components["Na+1"]
+  Solution.components.update("Na[+]": 2)
+  Solution.components["Na[+1]"]
+  ```
+
+- Arguments to `Solution.get_property` can be entered in any format you prefer. When `pyEQL` queries the database, it will automatically standardize the formula.
+- Property data in the database is uniquely identified by the standardized ion formula (output of `Ion.from_formula(<formula>).reduced_formla`, e.g. "Na[+1]" for sodium ion).

src/pyEQL/engines.py

@@ -7,8 +7,6 @@
 """
 from abc import ABC, abstractmethod
 
-from pymatgen.core.ion import Ion
-
 # internal pyEQL imports
 import pyEQL.activity_correction as ac
 
@@ -17,6 +15,7 @@
 from pyEQL import ureg
 from pyEQL.logging_system import logger
 from pyEQL.salt_ion_match import generate_salt_list
+from pyEQL.utils import standardize_formula
 
 
 class EOS(ABC):
@@ -192,7 +191,7 @@ def get_activity_coefficient(self, solution, solute):
         """
         # identify the predominant salt that this ion is a member of
         Salt = None
-        rform = Ion.from_formula(solute).reduced_formula
+        rform = standardize_formula(solute)
         salt_list = generate_salt_list(solution, unit="mol/kg")
         for item in salt_list:
             if rform == item.cation or rform == item.anion:
@@ -456,18 +455,9 @@ def get_solute_volume(self, solution):
         """Return the volume of the solutes."""
         # identify the predominant salt in the solution
         salt = solution.get_salt()
-        # reverse-convert the sanitized formula back to whatever was in self.components
-        for i in solution.components:
-            rform = Ion.from_formula(i).reduced_formula
-            if rform == salt.cation:
-                cation = i
-            if rform == salt.anion:
-                anion = i
-
         solute_vol = ureg.Quantity("0 L")
 
         # use the pitzer approach if parameters are available
-
         pitzer_calc = False
 
         param = solution.get_property(salt.formula, "model_parameters.molar_volume_pitzer")
@@ -476,7 +466,7 @@ def get_solute_volume(self, solution):
             # this is necessary for solutions inside e.g. an ion exchange
             # membrane, where the cation and anion concentrations may be
             # unequal
-            molality = (solution.get_amount(cation, "mol/kg") + solution.get_amount(anion, "mol/kg")) / 2
+            molality = (solution.get_amount(salt.cation, "mol/kg") + solution.get_amount(salt.anion, "mol/kg")) / 2
 
             # determine alpha1 and alpha2 based on the type of salt
             # see the May reference for the rules used to determine
@@ -512,8 +502,8 @@ def get_solute_volume(self, solution):
             solute_vol += (
                 apparent_vol
                 * (
-                    solution.get_amount(cation, "mol") / salt.nu_cation
-                    + solution.get_amount(anion, "mol") / salt.nu_anion
+                    solution.get_amount(salt.cation, "mol") / salt.nu_cation
+                    + solution.get_amount(salt.anion, "mol") / salt.nu_anion
                 )
                 / 2
             )
@@ -530,7 +520,7 @@ def get_solute_volume(self, solution):
                 continue
 
             # ignore the salt cation and anion, if already accounted for by Pitzer
-            if pitzer_calc is True and solute in [anion, cation]:
+            if pitzer_calc is True and solute in [salt.anion, salt.cation]:
                 continue
 
             part_vol = solution.get_property(solute, "size.molar_volume")

src/pyEQL/salt_ion_match.py

@@ -14,6 +14,7 @@
 from pymatgen.core.ion import Ion
 
 from pyEQL.logging_system import logger
+from pyEQL.utils import standardize_formula
 
 
 class Salt:
@@ -44,9 +45,9 @@ def __init__(self, cation, anion):
         # create pymatgen Ion objects
         pmg_cat = Ion.from_formula(cation)
         pmg_an = Ion.from_formula(anion)
-        # sanitize the cation and anion formulas
-        self.cation = pmg_cat.reduced_formula
-        self.anion = pmg_an.reduced_formula
+        # standardize the cation and anion formulas
+        self.cation = standardize_formula(cation)
+        self.anion = standardize_formula(anion)
 
         # get the charges on cation and anion
         self.z_cation = pmg_cat.charge
@@ -166,12 +167,12 @@ def identify_salt(sol):
     anion = "OH-"
 
     # return water if there are no solutes
-    if len(sort_list) < 3 and sort_list[0] == "H2O":
+    if len(sort_list) < 3 and sort_list[0] == "H2O(aq)":
         logger.info("Salt matching aborted because there are not enough solutes.")
         return Salt(cation, anion)
 
     # warn if something other than water is the predominant component
-    if sort_list[0] != "H2O":
+    if sort_list[0] != "H2O(aq)":
         logger.warning("H2O is not the most prominent component")
 
     # take the dominant cation and anion and assemble a salt from them

src/pyEQL/solution.py

@@ -26,6 +26,7 @@
 # logging system
 from pyEQL.logging_system import logger
 from pyEQL.salt_ion_match import generate_salt_list, identify_salt
+from pyEQL.utils import FormulaDict, standardize_formula
 
 EQUIV_WT_CACO3 = 100.09 / 2 * ureg.Quantity("g/mol")
 
@@ -126,7 +127,7 @@ def __init__(
 
         # create an empty dictionary of components. This dict comprises {formula: moles}
         #  where moles is the number of moles in the solution.
-        self.components: dict = {}
+        self.components = FormulaDict({})
 
         # connect to the desired property database
         if database is None:
@@ -163,7 +164,7 @@ def __init__(
             raise ValueError("Multiple solvents are not yet supported!")
         if solvent[0] not in ["H2O", "H2O(aq)", "water", "Water", "HOH"]:
             raise ValueError("Non-aqueous solvent detected. These are not yet supported!")
-        self.solvent = solvent[0]
+        self.solvent = standardize_formula(solvent[0])
 
         # TODO - do I need the ability to specify the solvent mass?
         # # raise an error if the solvent volume has also been given
@@ -415,7 +416,7 @@ def dielectric_constant(self) -> Quantity:
         denominator = 1
         for item in self.components:
             # ignore water
-            if item != "H2O":
+            if item != "H2O(aq)":
                 # skip over solutes that don't have parameters
                 # try:
                 fraction = self.get_amount(item, "fraction")
@@ -499,10 +500,6 @@ def viscosity_kinematic(self) -> Quantity:
         """
         # identify the main salt in the solution
         salt = self.get_salt()
-        # reverse-convert the sanitized formula back to whatever was in self.components
-        for i in self.components:
-            if Ion.from_formula(i).reduced_formula == salt.cation:
-                cation = i
 
         a0 = a1 = b0 = b1 = 0
 
@@ -532,7 +529,8 @@ def viscosity_kinematic(self) -> Quantity:
         MW_w = ureg.Quantity(self.get_property(self.solvent, "molecular_weight"))
 
         # calculate the cation mole fraction
-        x_cat = self.get_amount(cation, "fraction")
+        # x_cat = self.get_amount(cation, "fraction")
+        x_cat = self.get_amount(salt.cation, "fraction")
 
         # calculate the kinematic viscosity
         nu = math.log(nu_w * MW_w / MW) + 15 * x_cat**2 + x_cat**3 * G_123 + 3 * x_cat * G_23 * (1 - 0.05 * x_cat)
@@ -730,9 +728,7 @@ def alkalinity(self) -> Quantity:
         acid_anions = {"Cl[-1]", "Br[-1]", "I[-1]", "SO4[-2]", "NO3[-1]", "ClO4[-1]", "ClO3[-1]"}
 
         for item in self.components:
-            # sanitize the formulas
-            rform = Ion.from_formula(item).reduced_formula
-            if rform in base_cations.union(acid_anions):
+            if item in base_cations.union(acid_anions):
                 z = self.get_property(item, "charge")
                 alkalinity += self.get_amount(item, "mol/L") * z
 
@@ -780,7 +776,7 @@ def total_dissolved_solids(self) -> Quantity:
         tds = ureg.Quantity("0 mg/L")
         for s in self.components:
             # ignore pure water and dissolved gases, but not CO2
-            if s in ["H2O", "H+", "OH-", "H2", "O2"]:
+            if s in ["H2O(aq)", "H[+1]", "OH[-1]"]:
                 continue
             tds += self.get_amount(s, "mg/L")
 
@@ -1565,7 +1561,7 @@ def get_activity(
 
         """
         # switch to the water activity function if the species is H2O
-        if solute == "H2O" or solute == "water":
+        if solute in ["H2O(aq)", "water", "H2O", "HOH"]:
             activity = self.get_water_activity()
         else:
             # determine the concentration units to use based on the desired scale
@@ -1659,7 +1655,7 @@ def get_water_activity(self) -> Quantity:
 
         concentration_sum = 0
         for item in self.components:
-            if item == "H2O":
+            if item == "H2O(aq)":
                 pass
             else:
                 concentration_sum += self.get_amount(item, "mol/kg").magnitude
@@ -1755,8 +1751,8 @@ def _get_property(self, solute: str, name: str) -> Any | None:
         base_temperature = ureg.Quantity("25 degC")
         # base_pressure = ureg.Quantity("1 atm")
 
-        # query the database using the sanitized formula
-        rform = Ion.from_formula(solute).reduced_formula
+        # query the database using the standardized formula
+        rform = standardize_formula(solute)
         # TODO - there seems to be a bug in mongomock / JSONStore wherein properties does
         # not properly return dot-notation fields, e.g. size.molar_volume will not be returned.
         # also $exists:True does not properly return dot notated fields.
@@ -2125,7 +2121,7 @@ def __add__(self, other: Solution):
 
         # retrieve the amount of each component in the parent solution and
         # store in a list.
-        mix_species = {}
+        mix_species = FormulaDict({})
         for sol, amt in self.components.items():
             mix_species.update({sol: f"{amt} mol"})
         for sol2, amt2 in other.components.items():
@@ -2150,7 +2146,7 @@ def __add__(self, other: Solution):
 
         # create a new solution
         return Solution(
-            mix_species,
+            mix_species.data,  # pass a regular dict instead of the FormulaDict
             volume=str(mix_vol),
             pressure=str(mix_pressure),
             temperature=str(mix_temperature.to("K")),

src/pyEQL/utils.py

@@ -0,0 +1,50 @@
+"""
+pyEQL utilities
+
+:copyright: 2023 by Ryan S. Kingsbury
+:license: LGPL, see LICENSE for more details.
+
+"""
+
+from collections import UserDict
+
+from pymatgen.core.ion import Ion
+
+
+def standardize_formula(formula: str):
+    """
+    Convert a chemical formula into standard form.
+
+    Args:
+        formula: the chemical formula to standardize.
+
+    Returns:
+        A standardized chemical formula
+
+    Raises:
+        ValueError if `formula` cannot be processed or is invalid.
+
+    Notes:
+        Currently this method standardizes formulae by passing them through pymatgen.core.ion.Ion.reduced_formula(). For ions, this means that 1) the
+        charge number will always be listed explicitly and 2) the charge number will be enclosed in square brackets to remove any ambiguity in the meaning of the formula. For example, 'Na+', 'Na+1', and 'Na[+]' will all
+        standardize to "Na[+1]"
+    """
+    return Ion.from_formula(formula).reduced_formula
+
+
+class FormulaDict(UserDict):
+    """
+    Automatically converts keys on get/set using pymatgen.core.Ion.from_formula(key).reduced_formula.
+
+    This allows getting/setting/updating of Solution.components using flexible
+    formula notation (e.g., "Na+", "Na+1", "Na[+]" all have the same effect)
+    """
+
+    def __getitem__(self, key):
+        return super().__getitem__(standardize_formula(key))
+
+    def __setitem__(self, key, value):
+        super().__setitem__(standardize_formula(key), value)
+
+    def __delitem__(self, key):
+        super().__delitem__(standardize_formula(key))