Merge pull request #36 from rkingsbury/elements

improvements in element / chemistry handling

CHANGELOG.md

@@ -5,6 +5,17 @@ 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]
+
+### Added
+
+- `Solution`: new properties `elements` and `chemical_system`, new function `get_el_amt_dict` to compute the total
+  number of moles of each element present in the Solution.
+
+### Fixed
+
+- Two issues with the formatting of the `H2O(aq)` entry in the database, `pyeql_db.json`
+
 ## [0.7.0] - 2023-08-22
 
 ### Changed

setup.cfg

@@ -54,7 +54,7 @@ install_requires =
     numpy
     scipy
     pint
-    pymatgen>=2022.0.17
+    pymatgen>2022.8.10
     iapws
     monty
     maggma

src/pyEQL/database/pyeql_db.json

@@ -3,9 +3,9 @@
         "formula": "H2O(aq)",
         "charge": 0,
         "molecular_weight": "18.01528 g/mol",
-        "elements": "[Element H, Element O]",
+        "elements": ["H", "O"],
         "chemsys": "H-O",
-        "pmg_ion": "H2 O1 (aq)",
+        "pmg_ion": {"H": 2.0, "O": 1.0, "charge": 0.0},
         "formula_html": "H<sub>2</sub>O",
         "formula_latex": "H$_{2}$O",
         "formula_hill": "H2 O",

src/pyEQL/solute.py

@@ -19,6 +19,8 @@
 import numpy as np
 from pymatgen.core.ion import Ion
 
+from pyEQL.utils import standardize_formula
+
 
 @dataclass
 class Datum:
@@ -105,14 +107,15 @@ def from_formula(cls, formula: str):
         of the IonDoc.
         """
         pmg_ion = Ion.from_formula(formula)
-        f = pmg_ion.reduced_formula
+        f, factor = pmg_ion.get_reduced_formula_and_factor()
+        rform = standardize_formula(formula)
         charge = int(pmg_ion.charge)
         els = [str(el) for el in pmg_ion.elements]
-        mw = f"{float(pmg_ion.weight)} g/mol"  # weight is a FloatWithUnit
+        mw = f"{float(pmg_ion.weight / factor)} g/mol"  # weight is a FloatWithUnit
         chemsys = pmg_ion.chemical_system
 
         return cls(
-            f,
+            rform,
             charge=charge,
             molecular_weight=mw,
             elements=els,

src/pyEQL/solution.py

@@ -429,6 +429,25 @@ def dielectric_constant(self) -> Quantity:
 
         return ureg.Quantity(di_water / denominator, "dimensionless")
 
+    @property
+    def chemical_system(self) -> str:
+        """
+        Return the chemical system of the Solution as a "-" separated list of elements, sorted alphabetically. For
+        example, a solution containing CaCO3 would have a chemical system of "C-Ca-H-O".
+        """
+        return "-".join(self.elements)
+
+    @property
+    def elements(self) -> list:
+        """
+        Return a list of elements that are present in the solution. For example,
+        a solution containing CaCO3 would return ["C", "Ca", "H", "O"]
+        """
+        els = []
+        for s in self.components:
+            els.extend(self.get_property(s, "elements"))
+        return sorted(set(els))
+
     # TODO - need tests for viscosity
     @property
     def viscosity_dynamic(self) -> Quantity:
@@ -1020,6 +1039,31 @@ def get_amount(self, solute: str, units: str = "mol/L") -> Quantity:
 
         raise ValueError(f"Unsupported unit {units} specified for get_amount")
 
+    def get_el_amt_dict(self):
+        """
+        Return a dict of Element: amount in mol
+
+        Elements (keys) are suffixed with their oxidation state in parentheses,
+        e.g. "Fe(2)", "Cl(-1)".
+        """
+        d = {}
+        for s, mol in self.components.items():
+            elements = self.get_property(s, "elements")
+            pmg_ion_dict = self.get_property(s, "pmg_ion")
+            oxi_states = self.get_property(s, "oxi_state_guesses")[0]
+
+            for el in elements:
+                # stoichiometric coefficient, mol element per mol solute
+                stoich = pmg_ion_dict.get(el)
+                oxi_state = oxi_states.get(el)
+                key = f"{el}({oxi_state})"
+                if d.get(key):
+                    d[key] += stoich * mol
+                else:
+                    d[key] = stoich * mol
+
+        return d
+
     def get_total_amount(self, element: str, units) -> Quantity:
         """
         Return the total amount of 'element' (across all solutes) in the solution.
@@ -1825,8 +1869,15 @@ def _get_property(self, solute: str, name: str) -> Any | None:
                 return doc["model_parameters"]["molar_volume_pitzer"]
             return None
 
+        if name == "molecular_weight":
+            return ureg.Quantity(doc.get(name))
+
         # for parameters not named above, just return the base value
-        val = doc.get(name) if not isinstance(doc.get(name), dict) else doc[name].get("value")
+        if name == "pmg_ion" or not isinstance(doc.get(name), dict):
+            # if the queried value is not a dict, it is a root level key and should be returned as is
+            return doc.get(name)
+
+        val = doc[name].get("value")
         # logger.warning("%s has not been corrected for solution conditions" % name)
         if val is not None:
             return ureg.Quantity(val)

tests/test_solute.py

@@ -17,3 +17,6 @@ def test_from_formula():
     assert s.n_elements == 1
     assert s.oxi_state_guesses == ({"Mg": 2.0},)
     assert s.molecular_weight == "24.305 g/mol"
+    s2 = Solute.from_formula("O6")
+    assert s2.formula == "O3(aq)"
+    assert s2.molecular_weight == "47.9982 g/mol"

tests/test_solution.py

@@ -171,6 +171,30 @@ def test_pressure_temperature(s5):
     assert s5.volume < intermediate_V
 
 
+def test_elements(s5, s6):
+    assert s6.elements == sorted({"Ag", "Br", "C", "Ca", "H", "Mg", "Na", "O", "S"})
+    assert s6.chemical_system == "-".join(s6.elements)
+    assert s5.chemical_system == "C-Ca-H-O"
+
+
+def test_get_el_amt_dict(s6):
+    """ """
+    water_mol = s6.components["H2O(aq)"]
+    # scale volume to 8L
+    s6 *= 8
+    d = s6.get_el_amt_dict()
+    for el, amt in zip(
+        ["H(1)", "O(-2)", "Ca(2)", "Mg(2)", "Na(1)", "Ag(1)", "C(4)", "S(6)", "Br(-1)"],
+        [water_mol * 2 * 8, (water_mol + 0.018 + 0.24) * 8, 0.008, 0.040, 0.08, 0.08, 0.048, 0.48, 0.16],
+    ):
+        assert np.isclose(d[el], amt, atol=1e-3)
+
+    s = Solution({"Fe+2": "1 mM", "Fe+3": "5 mM", "FeCl2": "1 mM", "FeCl3": "5 mM"})
+    d = s.get_el_amt_dict()
+    for el, amt in zip(["Fe(2)", "Fe(3)", "Cl(-1)"], [0.002, 0.01, 0.002 + 0.015]):
+        assert np.isclose(d[el], amt, atol=1e-3)
+
+
 def test_p(s2):
     assert np.isclose(s2.p("Na+"), -1 * np.log10(s2.get_activity("Na+")))
     assert np.isclose(s2.p("Na+", activity=False), -1 * np.log10(s2.get_amount("Na+", "M").magnitude))

tests/test_utils.py

@@ -15,6 +15,7 @@ def test_standardize_formula():
     assert standardize_formula("Na[+]") == "Na[+1]"
     assert standardize_formula("SO4--") == "SO4[-2]"
     assert standardize_formula("Mg+2") == "Mg[+2]"
+    assert standardize_formula("O2") == "O2(aq)"
 
 
 def test_formula_dict():