From 4f01bd28e329e0b3715d60de08487eb2a2d88b93 Mon Sep 17 00:00:00 2001 From: Vineet Bansal Date: Wed, 17 Dec 2025 16:01:09 -0500 Subject: [PATCH 1/7] support for accumulating strings to run; some tests with USER_PUNCH in preparation for speciation --- src/pyEQL/pyEQL-phreeqc/pyEQL_phreeqc/core.py | 28 ++- tests/phreeqc/test_phreeqc.py | 182 ++++++++++++++++-- 2 files changed, 193 insertions(+), 17 deletions(-) diff --git a/src/pyEQL/pyEQL-phreeqc/pyEQL_phreeqc/core.py b/src/pyEQL/pyEQL-phreeqc/pyEQL_phreeqc/core.py index 4f0c26eb..03b6e7b6 100644 --- a/src/pyEQL/pyEQL-phreeqc/pyEQL_phreeqc/core.py +++ b/src/pyEQL/pyEQL-phreeqc/pyEQL_phreeqc/core.py @@ -1,4 +1,6 @@ +from inspect import cleandoc from pathlib import Path +from textwrap import dedent, indent from typing import Any from pyEQL_phreeqc._bindings import PyIPhreeqc @@ -14,6 +16,7 @@ def __init__(self, database: str = "phreeqc.dat", database_directory: Path | Non database_directory = Path(__file__).parent / "database" self._ext.load_database(str(database_directory / database)) + self._str = "" self._solutions: list[Solution] = [] # TODO: Is VAR the common denominator for most operations? @@ -32,22 +35,41 @@ def __getattr__(self, item) -> None: return getattr(self._ext, item) raise AttributeError(f"Phreeqc has no attribute '{item}'") + def __call__(self, *args, **kwargs): + self.run_string(self._str) + + def __str__(self): + return self._str + + def accumulate(self, s: str) -> None: + self._str += dedent(s) + def add_solution(self, solution_dict: dict) -> None: solution = Solution(solution_dict) index = len(self) - self.run_string(f""" + template = ( + "\n" + + cleandoc(f""" SOLUTION {index} - {solution} + {{solution}} SAVE SOLUTION {index} END """) + + "\n" + ) + _str = template.format(solution=indent(str(solution), " ")) + self.accumulate(_str) self._solutions.append(solution) def remove_solution(self, index: int) -> Solution: - self.run_string(f""" + _str = ( + cleandoc(f""" DELETE -solution {index} """) + + "\n" + ) + self.accumulate(_str) return self._solutions.pop(index) diff --git a/tests/phreeqc/test_phreeqc.py b/tests/phreeqc/test_phreeqc.py index 1938353e..620f9db8 100644 --- a/tests/phreeqc/test_phreeqc.py +++ b/tests/phreeqc/test_phreeqc.py @@ -1,4 +1,4 @@ -import textwrap +from inspect import cleandoc from pathlib import Path import numpy as np @@ -150,8 +150,7 @@ def test_run_add_delete_solution(): def test_run_dumpstring(): phreeqc = Phreeqc() - phreeqc.run_string( - textwrap.dedent(""" + phreeqc.run_string(""" SOLUTION 0 temp 25.0 REACTION 1 @@ -161,24 +160,21 @@ def test_run_dumpstring(): SAVE SOLUTION 0 END """) - ) phreeqc.set_dump_string_on(1) - phreeqc.run_string( - textwrap.dedent(""" + phreeqc.run_string(""" DUMP -solution 0 END """) - ) dump_string = phreeqc.get_dump_string() phreeqc.set_dump_string_on(0) # Due to platform specific differences, we only compare the first token # from each of the lines below, to the output. - expected = textwrap.dedent(""" + expected = cleandoc(""" SOLUTION_RAW 0 Solution after simulation 1. -temp 25 -pressure 1 @@ -216,7 +212,7 @@ def test_run_dumpstring(): USE reaction none USE reaction_temperature none USE reaction_pressure none - """).lstrip("\n") + """) dump_lines = dump_string.splitlines() expected_lines = expected.splitlines() @@ -232,8 +228,7 @@ def test_run_dumpstring(): def test_run_logstring(): phreeqc = Phreeqc() phreeqc.set_log_string_on(1) - phreeqc.run_string( - textwrap.dedent(""" + phreeqc.run_string(""" KNOBS -logfile true SOLUTION 0 @@ -245,13 +240,13 @@ def test_run_logstring(): SAVE SOLUTION 0 END """) - ) + log_string = phreeqc.get_log_string() phreeqc.set_log_string_on(0) # Due to platform specific differences, we only compare the first token # from each of the lines below, to the output. - expected = textwrap.dedent(""" + expected = cleandoc(""" ------------------------------------------- Beginning of initial solution calculations. ------------------------------------------- @@ -295,7 +290,7 @@ def test_run_logstring(): --------------------------------- End of Run after X Seconds. --------------------------------- - """).strip("\n") + """) log_lines = log_string.strip("\n").splitlines() expected_lines = expected.splitlines() @@ -315,3 +310,162 @@ def test_run_logstring(): continue assert got_first == exp_first + + +def test_speciation_one_solution(): + phreeqc = Phreeqc() + + phreeqc.run_string(""" + SELECTED_OUTPUT + -reset false + + USER_PUNCH + 10 t = SYS("aq", count, name$, type$, moles) + 20 FOR i = 1 to count + 30 PUNCH name$(i), MOL(name$(i)), ACT(name$(i)) + 40 NEXT i + + SOLUTION 0 + temp 25.0 + units mol/kgw + pH 7.0 + pe 8.5 + redox pe + water 0.9970480319717386 + END + """) + + # heading + soln 0 data (regardless of whether we have -headings in USER_PUNCH or not) + assert phreeqc.get_selected_output_row_count() == 2 + # [, , ], repeated for each (4) species + assert phreeqc.get_selected_output_column_count() == 12 + + +def test_speciation_two_solutions(): + phreeqc = Phreeqc() + + phreeqc.run_string(""" + SELECTED_OUTPUT + -reset false + + USER_PUNCH + 10 t = SYS("aq", count, name$, type$, moles) + 20 FOR i = 1 to count + 30 PUNCH name$(i), MOL(name$(i)), ACT(name$(i)) + 40 NEXT i + + SOLUTION 0 + temp 25.0 + units mol/kgw + pH 7.0 + pe 8.5 + redox pe + water 0.9970480319717386 + SAVE SOLUTION 0 + END + + SOLUTION 1 + temp 50.0 + units mol/kgw + pH 10.0 + pe 8.5 + redox pe + water 0.9970480319717386 + SAVE SOLUTION 1 + END + """) + + # heading + soln 0 data + soln 1 data (regardless of whether we have -headings in USER_PUNCH or not) + assert phreeqc.get_selected_output_row_count() == 3 + # [, , ], repeated for each (4) species + assert phreeqc.get_selected_output_column_count() == 12 + + +def test_accumulate(): + phreeqc = Phreeqc() + + phreeqc.accumulate( + cleandoc(""" + SELECTED_OUTPUT + -reset false + + USER_PUNCH + 10 t = SYS("aq", count, name$, type$, moles) + 20 FOR i = 1 to count + 30 PUNCH name$(i), MOL(name$(i)), ACT(name$(i)) + 40 NEXT i + """) + + "\n" + ) + + phreeqc.add_solution( + {"pH": 7.0, "pe": 8.5, "redox": "pe", "temp": 25.0, "units": "mol/kgw", "water": 0.9970480319717386} + ) + phreeqc.add_solution( + {"pH": 10.0, "pe": 8.5, "redox": "pe", "temp": 50.0, "units": "mol/kgw", "water": 0.9970480319717386} + ) + + expected = ( + cleandoc(""" + SELECTED_OUTPUT + -reset false + + USER_PUNCH + 10 t = SYS("aq", count, name$, type$, moles) + 20 FOR i = 1 to count + 30 PUNCH name$(i), MOL(name$(i)), ACT(name$(i)) + 40 NEXT i + + SOLUTION 0 + pH 7.0 + pe 8.5 + redox pe + temp 25.0 + units mol/kgw + water 0.9970480319717386 + SAVE SOLUTION 0 + END + + SOLUTION 1 + pH 10.0 + pe 8.5 + redox pe + temp 50.0 + units mol/kgw + water 0.9970480319717386 + SAVE SOLUTION 1 + END + """) + + "\n" + ) + + assert str(phreeqc) == expected + + +def test_speciation_add_solutions(): + phreeqc = Phreeqc() + + phreeqc.accumulate(""" + SELECTED_OUTPUT + -reset false + + USER_PUNCH + 10 t = SYS("aq", count, name$, type$, moles) + 20 FOR i = 1 to count + 30 PUNCH name$(i), MOL(name$(i)), ACT(name$(i)) + 40 NEXT i + """) + + phreeqc.add_solution( + {"pH": 7.0, "pe": 8.5, "redox": "pe", "temp": 25.0, "units": "mol/kgw", "water": 0.9970480319717386} + ) + phreeqc.add_solution( + {"pH": 10.0, "pe": 8.5, "redox": "pe", "temp": 50.0, "units": "mol/kgw", "water": 0.9970480319717386} + ) + # execute! + phreeqc() + + # heading + soln 0 data + soln 1 data (regardless of whether we have -headings in USER_PUNCH or not) + assert phreeqc.get_selected_output_row_count() == 3 + # [, , ], repeated for each (4) species + assert phreeqc.get_selected_output_column_count() == 12 From 17b827ec133059dea1253a1dbdbbe66b97c6b36d Mon Sep 17 00:00:00 2001 From: Vineet Bansal Date: Wed, 17 Dec 2025 17:05:54 -0500 Subject: [PATCH 2/7] a speciate method for the wrapper --- src/pyEQL/pyEQL-phreeqc/pyEQL_phreeqc/core.py | 52 ++++++++++++++++++- tests/phreeqc/test_phreeqc.py | 47 ++++++++++++++--- tests/phreeqc/test_phreeqc_solution.py | 4 +- 3 files changed, 94 insertions(+), 9 deletions(-) diff --git a/src/pyEQL/pyEQL-phreeqc/pyEQL_phreeqc/core.py b/src/pyEQL/pyEQL-phreeqc/pyEQL_phreeqc/core.py index 03b6e7b6..34304c37 100644 --- a/src/pyEQL/pyEQL-phreeqc/pyEQL_phreeqc/core.py +++ b/src/pyEQL/pyEQL-phreeqc/pyEQL_phreeqc/core.py @@ -41,11 +41,13 @@ def __call__(self, *args, **kwargs): def __str__(self): return self._str + def clear(self): + self._str = "" + def accumulate(self, s: str) -> None: self._str += dedent(s) - def add_solution(self, solution_dict: dict) -> None: - solution = Solution(solution_dict) + def add_solution(self, solution: Solution) -> None: index = len(self) template = ( "\n" @@ -72,6 +74,52 @@ def remove_solution(self, index: int) -> Solution: self.accumulate(_str) return self._solutions.pop(index) + def speciate( + self, solutions: Solution | list[Solution], props: tuple[str] | None = None + ) -> dict[int, dict[str, Any]]: + if props is None: + props = ("MOL", "ACT") + punch_line = ", ".join([f"{prop}(name$(i))" for prop in props]) + + self.clear() + self.accumulate(f""" + SELECTED_OUTPUT + -reset false + + USER_PUNCH + 10 t = SYS("aq", count, name$, type$, moles) + 20 FOR i = 1 to count + 30 PUNCH name$(i), {punch_line} + 40 NEXT i + """) + + if isinstance(solutions, Solution): + solutions = [solutions] + + for solution in solutions: + self.add_solution(solution) + + self() + + # first line is always the header, but ensure this so we can skip it. + assert self.output[0][0].startswith("no_heading") + + all_solution_species: dict[int, dict[str, Any]] = {} + for solution_i, line in enumerate(self.output[1:]): + this_solution_species = {} + n_tokens = len(line) + j = 0 + while j < n_tokens: + species = line[j] + this_solution_species[species] = {} + j += 1 + for prop in props: + this_solution_species[species][prop] = line[j] + j += 1 + all_solution_species[solution_i] = this_solution_species + + return all_solution_species + class PhreeqcOutput: def __init__(self, phreeqc: Phreeqc): diff --git a/tests/phreeqc/test_phreeqc.py b/tests/phreeqc/test_phreeqc.py index 620f9db8..e872f95d 100644 --- a/tests/phreeqc/test_phreeqc.py +++ b/tests/phreeqc/test_phreeqc.py @@ -3,6 +3,7 @@ import numpy as np from pyEQL_phreeqc import Phreeqc +from pyEQL_phreeqc.solution import Solution def test_load_database_internal(): @@ -132,17 +133,19 @@ def test_run_simple(): def test_run_add_solution(): phreeqc = Phreeqc() - phreeqc.add_solution( + solution = Solution( {"pH": 7.0, "pe": 8.5, "redox": "pe", "temp": 25.0, "units": "mol/kgw", "water": 0.9970480319717386} ) + phreeqc.add_solution(solution) assert len(phreeqc) == 1 def test_run_add_delete_solution(): phreeqc = Phreeqc() - phreeqc.add_solution( + solution = Solution( {"pH": 7.0, "pe": 8.5, "redox": "pe", "temp": 25.0, "units": "mol/kgw", "water": 0.9970480319717386} ) + phreeqc.add_solution(solution) phreeqc.remove_solution(0) assert len(phreeqc) == 0 @@ -399,10 +402,10 @@ def test_accumulate(): ) phreeqc.add_solution( - {"pH": 7.0, "pe": 8.5, "redox": "pe", "temp": 25.0, "units": "mol/kgw", "water": 0.9970480319717386} + Solution({"pH": 7.0, "pe": 8.5, "redox": "pe", "temp": 25.0, "units": "mol/kgw", "water": 0.9970480319717386}) ) phreeqc.add_solution( - {"pH": 10.0, "pe": 8.5, "redox": "pe", "temp": 50.0, "units": "mol/kgw", "water": 0.9970480319717386} + Solution({"pH": 10.0, "pe": 8.5, "redox": "pe", "temp": 50.0, "units": "mol/kgw", "water": 0.9970480319717386}) ) expected = ( @@ -457,10 +460,10 @@ def test_speciation_add_solutions(): """) phreeqc.add_solution( - {"pH": 7.0, "pe": 8.5, "redox": "pe", "temp": 25.0, "units": "mol/kgw", "water": 0.9970480319717386} + Solution({"pH": 7.0, "pe": 8.5, "redox": "pe", "temp": 25.0, "units": "mol/kgw", "water": 0.9970480319717386}) ) phreeqc.add_solution( - {"pH": 10.0, "pe": 8.5, "redox": "pe", "temp": 50.0, "units": "mol/kgw", "water": 0.9970480319717386} + Solution({"pH": 10.0, "pe": 8.5, "redox": "pe", "temp": 50.0, "units": "mol/kgw", "water": 0.9970480319717386}) ) # execute! phreeqc() @@ -469,3 +472,35 @@ def test_speciation_add_solutions(): assert phreeqc.get_selected_output_row_count() == 3 # [, , ], repeated for each (4) species assert phreeqc.get_selected_output_column_count() == 12 + + +def test_speciate(): + solutions = [ + Solution({"pH": 7.0, "pe": 8.5, "redox": "pe", "temp": 25.0, "units": "mol/kgw", "water": 0.9970480319717386}), + Solution({"pH": 10.0, "pe": 8.5, "redox": "pe", "temp": 50.0, "units": "mol/kgw", "water": 0.9970480319717386}), + ] + phreeqc = Phreeqc() + all_solution_species = phreeqc.speciate(solutions) + + expected_species = { + 0: { + "H+": {"ACT": 1.0001522689856982e-07, "MOL": 1.0005246407839175e-07}, + "H2": {"ACT": 7.079457681907915e-35, "MOL": 7.079457517820301e-35}, + "O2": {"ACT": 8.317637520771417e-25, "MOL": 8.317637327985348e-25}, + "OH-": {"ACT": 1.0123126727760366e-07, "MOL": 1.0126897880811404e-07}, + }, + 1: { + "H+": {"ACT": 1e-10, "MOL": 1.0197827284798617e-10}, + "H2": {"ACT": 5.626700758118202e-41, "MOL": 0.0}, + "O2": {"ACT": 3.659468532125681e-05, "MOL": 3.6592332322612637e-05}, + "OH-": {"ACT": 0.0005473549298676792, "MOL": 0.0005585111533405827}, + }, + } + + assert len(all_solution_species) == len(expected_species) + for solution_index, solution_species in all_solution_species.items(): + assert solution_index in expected_species + assert set(solution_species.keys()) == set(expected_species[solution_index].keys()) + for k in solution_species: + for prop in ("ACT", "MOL"): + assert np.isclose(solution_species[k][prop], expected_species[solution_index][k][prop]) diff --git a/tests/phreeqc/test_phreeqc_solution.py b/tests/phreeqc/test_phreeqc_solution.py index 89030606..937e34fb 100644 --- a/tests/phreeqc/test_phreeqc_solution.py +++ b/tests/phreeqc/test_phreeqc_solution.py @@ -1,9 +1,10 @@ from pyEQL_phreeqc import Phreeqc +from pyEQL_phreeqc.solution import Solution def test_add_solution(): phreeqc = Phreeqc() - phreeqc.add_solution( + solution = Solution( { "Cl": "4.011842831773806", "Na": "4.011842831773806", @@ -15,3 +16,4 @@ def test_add_solution(): "water": 0.9970480319717386, } ) + phreeqc.add_solution(solution) From 2e2f6d05707136af57858aa72a7117f365602d4d Mon Sep 17 00:00:00 2001 From: Vineet Bansal Date: Thu, 18 Dec 2025 14:32:58 -0500 Subject: [PATCH 3/7] support for molality/activity (species specific) and solution level properties (osmotic coefficient) --- src/pyEQL/engines.py | 91 ++++++++++- .../pyEQL-phreeqc/pyEQL_phreeqc/__init__.py | 3 +- src/pyEQL/pyEQL-phreeqc/pyEQL_phreeqc/core.py | 96 +++++++----- .../pyEQL-phreeqc/pyEQL_phreeqc/solution.py | 33 +++- tests/phreeqc/test_phreeqc.py | 144 ++++++++++-------- tests/phreeqc/test_phreeqc_solution.py | 9 +- 6 files changed, 265 insertions(+), 111 deletions(-) diff --git a/src/pyEQL/engines.py b/src/pyEQL/engines.py index 5eecead5..0a2164c0 100644 --- a/src/pyEQL/engines.py +++ b/src/pyEQL/engines.py @@ -899,16 +899,99 @@ def __init__( self._stored_comp = None def _setup_ppsol(self, solution: "solution.Solution") -> None: - raise NotImplementedError + """Helper method to set up a PhreeqPython solution for subsequent analysis.""" + + from pyEQL_phreeqc import Solution # noqa: PLC0415 + + self._stored_comp = solution.components.copy() + solv_mass = solution.solvent_mass.to("kg").magnitude + # inherit bulk solution properties + d = { + "temp": solution.temperature.to("degC").magnitude, + "units": "mol/kgw", # to avoid confusion about volume, use mol/kgw which seems more robust in PHREEQC + "pH": solution.pH, + "pe": solution.pE, + "redox": "pe", # hard-coded to use the pe + # PHREEQC will assume 1 kg if not specified, there is also no direct way to specify volume, so we + # really have to specify the solvent mass in 1 liter of solution + "water": solv_mass, + } + if solution.balance_charge == "pH": + d["pH"] = str(d["pH"]) + " charge" + if solution.balance_charge == "pE": + d["pe"] = str(d["pe"]) + " charge" + + # add the composition to the dict + # also, skip H and O + for el, mol in solution.get_el_amt_dict().items(): + # CAUTION - care must be taken to avoid unintended behavior here. get_el_amt_dict() will return + # all distinct oxi states of each element present. If there are elements present whose oxi states + # are NOT recognized by PHREEQC (via SPECIAL_ELEMENTS) then the amount of only 1 oxi state will be + # entered into the composition dict. This can especially cause problems after equilibrate() has already + # been called once. For example, equilibrating a simple NaCl solution generates Cl species that are assigned + # various oxidations states, -1 mostly, but also 1, 2, and 3. Since the concentrations of everything + # except the -1 oxi state are tiny, this can result in Cl "disappearing" from the solution if + # equlibrate is called again. It also causes non-determinism, because the amount is taken from whatever + # oxi state happens to be iterated through last. + + # strip off the oxi state + bare_el = el.split("(")[0] + if bare_el in SPECIAL_ELEMENTS: + # PHREEQC will ignore float-formatted oxi states. Need to make sure we are + # passing, e.g. 'C(4)' and not 'C(4.0)' + key = f"{bare_el}({int(float(el.split('(')[-1].split(')')[0]))})" + elif bare_el in ["H", "O"]: + continue + else: + key = bare_el + + if key in d: + # when multiple oxi states for the same (non-SPECIAL) element are present, make sure to + # add all their amounts together + d[key] += str(mol / solv_mass) + else: + d[key] = str(mol / solv_mass) + + # tell PHREEQC which species to use for charge balance + if solution.balance_charge is not None and solution._cb_species in solution.get_components_by_element()[el]: + d[key] += " charge" + + ppsol = self.pp.add_solution(Solution(d)) + self.ppsol = ppsol def _destroy_ppsol(self) -> None: - raise NotImplementedError + if self.ppsol is not None: + self.pp.remove_solution(0) # TODO: Are we only expecting a single solution per wrapper? + self.ppsol = None def get_activity_coefficient(self, solution: "solution.Solution", solute: str) -> ureg.Quantity: - raise NotImplementedError + """ + Return the *molal scale* activity coefficient of solute, given a Solution + object. + """ + if (self.ppsol is None) or (solution.components != self._stored_comp): + self._destroy_ppsol() + self._setup_ppsol(solution) + + # translate the species into keys that phreeqc will understand + k = standardize_formula(solute) + spl = k.split("[") + el = spl[0] + chg = spl[1].split("]")[0] + if chg[-1] == "1": + chg = chg[0] # just pass + or -, not +1 / -1 + k = el + chg + + # calculate the molal scale activity coefficient + # act = self.ppsol.activity(k, "mol") / self.ppsol.molality(k, "mol") + act = self.ppsol.get_activity(k) / self.ppsol.get_molality(k) + + return ureg.Quantity(act, "dimensionless") def get_osmotic_coefficient(self, solution: "solution.Solution") -> ureg.Quantity: - raise NotImplementedError + osmotic = self.ppsol.get_osmotic_coefficient() # noqa: F841 + # TODO: Why is the wrapper defaulting to 1 instead of 0 (default behavior for phreeqc) for phreeqpython? + return ureg.Quantity(1, "dimensionless") def get_solute_volume(self, solution: "solution.Solution") -> ureg.Quantity: """Return the volume of the solutes.""" diff --git a/src/pyEQL/pyEQL-phreeqc/pyEQL_phreeqc/__init__.py b/src/pyEQL/pyEQL-phreeqc/pyEQL_phreeqc/__init__.py index db8b7556..459a055d 100644 --- a/src/pyEQL/pyEQL-phreeqc/pyEQL_phreeqc/__init__.py +++ b/src/pyEQL/pyEQL-phreeqc/pyEQL_phreeqc/__init__.py @@ -1,3 +1,4 @@ from .core import Phreeqc +from .solution import Solution -__all__ = ["Phreeqc"] +__all__ = ["Phreeqc", "Solution"] diff --git a/src/pyEQL/pyEQL-phreeqc/pyEQL_phreeqc/core.py b/src/pyEQL/pyEQL-phreeqc/pyEQL_phreeqc/core.py index 34304c37..4ad47ff0 100644 --- a/src/pyEQL/pyEQL-phreeqc/pyEQL_phreeqc/core.py +++ b/src/pyEQL/pyEQL-phreeqc/pyEQL_phreeqc/core.py @@ -29,6 +29,9 @@ def __init__(self, database: str = "phreeqc.dat", database_directory: Path | Non def __len__(self): return len(self._solutions) + def __getitem__(self, item): + return self._solutions[item] + def __getattr__(self, item) -> None: """Delegate attribute access to the underlying PyIPhreeqc instance.""" if hasattr(self._ext, item): @@ -39,7 +42,7 @@ def __call__(self, *args, **kwargs): self.run_string(self._str) def __str__(self): - return self._str + return cleandoc(self._str) def clear(self): self._str = "" @@ -47,21 +50,29 @@ def clear(self): def accumulate(self, s: str) -> None: self._str += dedent(s) - def add_solution(self, solution: Solution) -> None: - index = len(self) - template = ( - "\n" - + cleandoc(f""" - SOLUTION {index} - {{solution}} - SAVE SOLUTION {index} - END - """) - + "\n" - ) - _str = template.format(solution=indent(str(solution), " ")) - self.accumulate(_str) - self._solutions.append(solution) + def _add_solution(self, solution: Solution | list[Solution]) -> Solution | list[Solution]: + singleton = isinstance(solution, Solution) + solutions = [solution] if singleton else solution + + for solution in solutions: + index = len(self) + solution._number = index + # TODO: This should go in the Solution class + template = ( + "\n" + + cleandoc(f""" + SOLUTION {index} + {{solution}} + SAVE SOLUTION {index} + END + """) + + "\n" + ) + _str = template.format(solution=indent(str(solution), " ")) + self.accumulate(_str) + self._solutions.append(solution) + + return self._solutions[-1] if singleton else self._solutions def remove_solution(self, index: int) -> Solution: _str = ( @@ -72,14 +83,22 @@ def remove_solution(self, index: int) -> Solution: + "\n" ) self.accumulate(_str) + self() return self._solutions.pop(index) - def speciate( - self, solutions: Solution | list[Solution], props: tuple[str] | None = None - ) -> dict[int, dict[str, Any]]: - if props is None: - props = ("MOL", "ACT") - punch_line = ", ".join([f"{prop}(name$(i))" for prop in props]) + def add_solution( + self, + solution: Solution | list[Solution], + solution_props: tuple[str] | None = None, + species_props: tuple[str] | None = None, + ) -> Solution | list[Solution]: + if solution_props is None: + solution_props = ("OSMOTIC",) + solution_punch_line = ", ".join(list(solution_props)) + + if species_props is None: + species_props = ("MOL", "ACT") + species_punch_line = ", ".join([f"{prop}(name$(i))" for prop in species_props]) self.clear() self.accumulate(f""" @@ -87,38 +106,43 @@ def speciate( -reset false USER_PUNCH + 5 PUNCH {solution_punch_line}, EOL_NOTAB$ 10 t = SYS("aq", count, name$, type$, moles) 20 FOR i = 1 to count - 30 PUNCH name$(i), {punch_line} + 30 PUNCH name$(i), {species_punch_line} 40 NEXT i """) - if isinstance(solutions, Solution): - solutions = [solutions] - - for solution in solutions: - self.add_solution(solution) + return_value = self._add_solution(solution) self() # first line is always the header, but ensure this so we can skip it. assert self.output[0][0].startswith("no_heading") - all_solution_species: dict[int, dict[str, Any]] = {} - for solution_i, line in enumerate(self.output[1:]): - this_solution_species = {} + output = self.output[:][1:] + for solution_i, line in enumerate(output): + solution_i_props = {"species": {}} n_tokens = len(line) j = 0 while j < n_tokens: + # everything before the first '\n' are solution_props + if not solution_i_props["species"]: + while (solution_prop := line[j]) != "\n": + solution_i_props[solution_props[j]] = solution_prop + j += 1 + j += 1 # skip "\n" + species = line[j] - this_solution_species[species] = {} + solution_i_props["species"][species] = {} j += 1 - for prop in props: - this_solution_species[species][prop] = line[j] + for prop in species_props: + solution_i_props["species"][species][prop] = line[j] j += 1 - all_solution_species[solution_i] = this_solution_species - return all_solution_species + self[solution_i]._set_calculated_props(solution_i_props) + + return return_value class PhreeqcOutput: diff --git a/src/pyEQL/pyEQL-phreeqc/pyEQL_phreeqc/solution.py b/src/pyEQL/pyEQL-phreeqc/pyEQL_phreeqc/solution.py index 683c5141..8c7e8f3e 100644 --- a/src/pyEQL/pyEQL-phreeqc/pyEQL_phreeqc/solution.py +++ b/src/pyEQL/pyEQL-phreeqc/pyEQL_phreeqc/solution.py @@ -1,5 +1,32 @@ -class Solution(dict): - # A solution is nothing but a dict of str: Any mapping +from copy import deepcopy +from typing import Any + + +class Solution: + def __init__(self, props): + self._number = -1 + self._input_props = deepcopy(props) + self._calculated_props = {} def __str__(self): - return "\n".join(f"{k} {v}" for k, v in self.items()) + return "\n".join(f"{k} {v}" for k, v in self._input_props.items()) + + def _set_calculated_props(self, props: dict[Any, Any]): + self._calculated_props = deepcopy(props) + + def _get_calculated_props(self): + return self._calculated_props + + def _get_calculated_prop(self, which, species: str | None = None): + if species is not None: + return self._calculated_props["species"][species][which] + return self._calculated_props[which] + + def get_activity(self, species): + return self._get_calculated_prop("ACT", species) + + def get_molality(self, species): + return self._get_calculated_prop("MOL", species) + + def get_osmotic_coefficient(self): + return self._get_calculated_prop("OSMOTIC") diff --git a/tests/phreeqc/test_phreeqc.py b/tests/phreeqc/test_phreeqc.py index e872f95d..5e8c1dc1 100644 --- a/tests/phreeqc/test_phreeqc.py +++ b/tests/phreeqc/test_phreeqc.py @@ -2,8 +2,8 @@ from pathlib import Path import numpy as np -from pyEQL_phreeqc import Phreeqc -from pyEQL_phreeqc.solution import Solution +from pyEQL_phreeqc import Phreeqc, Solution +from pytest import approx def test_load_database_internal(): @@ -384,36 +384,26 @@ def test_speciation_two_solutions(): assert phreeqc.get_selected_output_column_count() == 12 -def test_accumulate(): +def test_add_solution_input(): phreeqc = Phreeqc() - phreeqc.accumulate( - cleandoc(""" - SELECTED_OUTPUT - -reset false - - USER_PUNCH - 10 t = SYS("aq", count, name$, type$, moles) - 20 FOR i = 1 to count - 30 PUNCH name$(i), MOL(name$(i)), ACT(name$(i)) - 40 NEXT i - """) - + "\n" - ) - phreeqc.add_solution( - Solution({"pH": 7.0, "pe": 8.5, "redox": "pe", "temp": 25.0, "units": "mol/kgw", "water": 0.9970480319717386}) - ) - phreeqc.add_solution( - Solution({"pH": 10.0, "pe": 8.5, "redox": "pe", "temp": 50.0, "units": "mol/kgw", "water": 0.9970480319717386}) + [ + Solution( + {"pH": 7.0, "pe": 8.5, "redox": "pe", "temp": 25.0, "units": "mol/kgw", "water": 0.9970480319717386} + ), + Solution( + {"pH": 10.0, "pe": 8.5, "redox": "pe", "temp": 50.0, "units": "mol/kgw", "water": 0.9970480319717386} + ), + ] ) - expected = ( - cleandoc(""" + expected = cleandoc(""" SELECTED_OUTPUT -reset false USER_PUNCH + 5 PUNCH OSMOTIC, EOL_NOTAB$ 10 t = SYS("aq", count, name$, type$, moles) 20 FOR i = 1 to count 30 PUNCH name$(i), MOL(name$(i)), ACT(name$(i)) @@ -439,39 +429,28 @@ def test_accumulate(): SAVE SOLUTION 1 END """) - + "\n" - ) assert str(phreeqc) == expected -def test_speciation_add_solutions(): +def test_add_solution_output(): phreeqc = Phreeqc() - phreeqc.accumulate(""" - SELECTED_OUTPUT - -reset false - - USER_PUNCH - 10 t = SYS("aq", count, name$, type$, moles) - 20 FOR i = 1 to count - 30 PUNCH name$(i), MOL(name$(i)), ACT(name$(i)) - 40 NEXT i - """) - - phreeqc.add_solution( - Solution({"pH": 7.0, "pe": 8.5, "redox": "pe", "temp": 25.0, "units": "mol/kgw", "water": 0.9970480319717386}) - ) phreeqc.add_solution( - Solution({"pH": 10.0, "pe": 8.5, "redox": "pe", "temp": 50.0, "units": "mol/kgw", "water": 0.9970480319717386}) + [ + Solution( + {"pH": 7.0, "pe": 8.5, "redox": "pe", "temp": 25.0, "units": "mol/kgw", "water": 0.9970480319717386} + ), + Solution( + {"pH": 10.0, "pe": 8.5, "redox": "pe", "temp": 50.0, "units": "mol/kgw", "water": 0.9970480319717386} + ), + ] ) - # execute! - phreeqc() # heading + soln 0 data + soln 1 data (regardless of whether we have -headings in USER_PUNCH or not) assert phreeqc.get_selected_output_row_count() == 3 - # [, , ], repeated for each (4) species - assert phreeqc.get_selected_output_column_count() == 12 + # , "\n", +[, , ] repeated for each (4) species + assert phreeqc.get_selected_output_column_count() == 14 def test_speciate(): @@ -480,27 +459,70 @@ def test_speciate(): Solution({"pH": 10.0, "pe": 8.5, "redox": "pe", "temp": 50.0, "units": "mol/kgw", "water": 0.9970480319717386}), ] phreeqc = Phreeqc() - all_solution_species = phreeqc.speciate(solutions) + phreeqc.add_solution(solutions) - expected_species = { + expected = { 0: { - "H+": {"ACT": 1.0001522689856982e-07, "MOL": 1.0005246407839175e-07}, - "H2": {"ACT": 7.079457681907915e-35, "MOL": 7.079457517820301e-35}, - "O2": {"ACT": 8.317637520771417e-25, "MOL": 8.317637327985348e-25}, - "OH-": {"ACT": 1.0123126727760366e-07, "MOL": 1.0126897880811404e-07}, + "OSMOTIC": 0.0, + "species": { + "H+": {"ACT": 1.0001522689856982e-07, "MOL": 1.0005246407839175e-07}, + "H2": {"ACT": 7.079457681907915e-35, "MOL": 7.079457517820301e-35}, + "O2": {"ACT": 8.317637520771417e-25, "MOL": 8.317637327985348e-25}, + "OH-": {"ACT": 1.0123126727760366e-07, "MOL": 1.0126897880811404e-07}, + }, }, 1: { - "H+": {"ACT": 1e-10, "MOL": 1.0197827284798617e-10}, - "H2": {"ACT": 5.626700758118202e-41, "MOL": 0.0}, - "O2": {"ACT": 3.659468532125681e-05, "MOL": 3.6592332322612637e-05}, - "OH-": {"ACT": 0.0005473549298676792, "MOL": 0.0005585111533405827}, + "OSMOTIC": 0.0, + "species": { + "H+": {"ACT": 1e-10, "MOL": 1.0197827284798617e-10}, + "H2": {"ACT": 5.626700758118202e-41, "MOL": 0.0}, + "O2": {"ACT": 3.659468532125681e-05, "MOL": 3.6592332322612637e-05}, + "OH-": {"ACT": 0.0005473549298676792, "MOL": 0.0005585111533405827}, + }, }, } - assert len(all_solution_species) == len(expected_species) - for solution_index, solution_species in all_solution_species.items(): - assert solution_index in expected_species - assert set(solution_species.keys()) == set(expected_species[solution_index].keys()) - for k in solution_species: + assert len(phreeqc) == len(expected) + for solution_index, solution_props in expected.items(): + props = phreeqc[solution_index]._get_calculated_props() + assert set(props.keys()) == set(expected[solution_index].keys()) + for k in expected[solution_index]: + if k != "species": + assert props[k] == approx(expected[solution_index][k]) + + for k in solution_props["species"]: for prop in ("ACT", "MOL"): - assert np.isclose(solution_species[k][prop], expected_species[solution_index][k][prop]) + assert np.isclose(props["species"][k][prop], expected[solution_index]["species"][k][prop]) + + +def test_speciate_molality(): + solution = Solution( + {"pH": 7.0, "pe": 8.5, "redox": "pe", "temp": 25.0, "units": "mol/kgw", "water": 0.9970480319717386} + ) + phreeqc = Phreeqc() + phreeqc.add_solution(solution) + + activity = phreeqc[0].get_molality("H+") + assert activity == approx(1.0005246407839175e-07) + + +def test_speciate_activity(): + solution = Solution( + {"pH": 7.0, "pe": 8.5, "redox": "pe", "temp": 25.0, "units": "mol/kgw", "water": 0.9970480319717386} + ) + phreeqc = Phreeqc() + phreeqc.add_solution(solution) + + activity = phreeqc[0].get_activity("H+") + assert activity == approx(1.0001522689856982e-07) + + +def test_osmotic_coefficient(): + solution = Solution( + {"pH": 7.0, "pe": 8.5, "redox": "pe", "temp": 25.0, "units": "mol/kgw", "water": 0.9970480319717386} + ) + phreeqc = Phreeqc() + phreeqc.add_solution(solution) + + osmotic_coefficient = phreeqc[0].get_osmotic_coefficient() + assert osmotic_coefficient == approx(0.0) diff --git a/tests/phreeqc/test_phreeqc_solution.py b/tests/phreeqc/test_phreeqc_solution.py index 937e34fb..d1a9f27c 100644 --- a/tests/phreeqc/test_phreeqc_solution.py +++ b/tests/phreeqc/test_phreeqc_solution.py @@ -1,10 +1,8 @@ -from pyEQL_phreeqc import Phreeqc -from pyEQL_phreeqc.solution import Solution +from pyEQL_phreeqc import Solution -def test_add_solution(): - phreeqc = Phreeqc() - solution = Solution( +def test_create_solution(): + Solution( { "Cl": "4.011842831773806", "Na": "4.011842831773806", @@ -16,4 +14,3 @@ def test_add_solution(): "water": 0.9970480319717386, } ) - phreeqc.add_solution(solution) From a9241c990087634f359e56cd6077468c40b2a01d Mon Sep 17 00:00:00 2001 From: Vineet Bansal Date: Mon, 12 Jan 2026 18:08:09 -0500 Subject: [PATCH 4/7] all properties needed for pyEQL use; test_pyeql_phreeqc.py corresponding to test_phreeqc.py --- src/pyEQL/engines.py | 137 +++++- src/pyEQL/pyEQL-phreeqc/pyEQL_phreeqc/core.py | 85 ++-- .../pyEQL-phreeqc/pyEQL_phreeqc/solution.py | 55 ++- tests/phreeqc/test_phreeqc.py | 106 ++++- tests/test_pyeql_phreeqc.py | 414 ++++++++++++++++++ 5 files changed, 742 insertions(+), 55 deletions(-) create mode 100644 tests/test_pyeql_phreeqc.py diff --git a/src/pyEQL/engines.py b/src/pyEQL/engines.py index 0a2164c0..0329be5f 100644 --- a/src/pyEQL/engines.py +++ b/src/pyEQL/engines.py @@ -869,7 +869,7 @@ def __init__( self, phreeqc_db: Literal[ "phreeqc.dat", "vitens.dat", "wateq4f_PWN.dat", "pitzer.dat", "llnl.dat", "geothermal.dat" - ] = "llnl.dat", + ] = "phreeqc.dat", ) -> None: """ Args: @@ -956,7 +956,16 @@ def _setup_ppsol(self, solution: "solution.Solution") -> None: if solution.balance_charge is not None and solution._cb_species in solution.get_components_by_element()[el]: d[key] += " charge" - ppsol = self.pp.add_solution(Solution(d)) + try: + ppsol = self.pp.add_solution(Solution(d)) + except Exception as e: + # catch problems with the input to phreeqc + raise ValueError( + "There is a problem with your input. The error message received from " + f" phreeqc is:\n\n {e}\n Check your input arguments, especially " + "the composition dictionary, and try again." + ) + self.ppsol = ppsol def _destroy_ppsol(self) -> None: @@ -998,5 +1007,125 @@ def get_solute_volume(self, solution: "solution.Solution") -> ureg.Quantity: # TODO - phreeqc seems to have no concept of volume, but it does calculate density return ureg.Quantity(0, "L") - def equilibrate(self, solution: "solution.Solution") -> None: - raise NotImplementedError + def equilibrate( + self, + solution: "solution.Solution", + atmosphere: bool = False, + solids: list[str] | None = None, + gases: dict[str, str | float] | None = None, + ) -> None: + """ + Adjust the speciation of a Solution object to achieve chemical equilibrium. + + Args: + atmosphere: + Boolean indicating whether to equilibrate the solution + w.r.t atmospheric gases. + solids: + A list of solids used to achieve liquid-solid equilibrium. Each + solid in this list should be present in the Phreeqc database. + We assume a target saturation index of 0 and an infinite + amount of material. + gases: + A dictionary of gases used to achieve liquid-gas equilibrium. + Each key denotes the gas species, and the corresponding value + denotes its concentration, as a log partial pressure value or + other interpretable pressure units. For example, the following + are equivalent (log10(0.000316) = -3.5) + {"CO2": "0.000316 atm"} + {"CO2": -3.5} + """ + if self.ppsol is not None: + self.ppsol.forget() + self._setup_ppsol(solution) + + # store the original solvent mass + orig_solvent_moles = solution.components[solution.solvent] + + # store the original solvent elements and components with that element + # that we may need to add back later. + orig_el_dict = solution.get_el_amt_dict(nested=True) + orig_components_by_element = solution.get_components_by_element(nested=True) + + # Use supplied gases, merged with atmospheric gases + # if atmosphere == True + gases = (ATMOSPHERE if atmosphere else {}) | (gases or {}) + + # Mapping from phase name to: + # (, ) tuples (for solids). + # (, ) tuples (for gases). + phases = {} + if solids is not None: + # Assume saturation index of 0 for all solids. + phases |= dict.fromkeys(solids, (0, EQUILIBRIUM_PHASE_AMOUNT)) + + for k, v in gases.items(): + v_quantity = ureg.Quantity(v) + if v_quantity.dimensionless: + log_partial_pressure = v_quantity.magnitude + else: + log_partial_pressure = math.log10(v_quantity.to("atm").magnitude) + phases |= {f"{k}(g)": (log_partial_pressure, EQUILIBRIUM_PHASE_AMOUNT)} + + if phases: + phase_names = list(phases.keys()) + saturation_indices = [v[0] for v in phases.values()] + amounts = [v[1] for v in phases.values()] + + try: + self.ppsol.equalize(phases=phase_names, saturation_indices=saturation_indices, amounts=amounts) + except Exception as e: + # Re-raise exception due to unrecognized phases as ValueError. + if "Phase not found in database" in str(e): + raise ValueError(str(e)) + raise e + + solution.components = FormulaDict({}) + # use the output from PHREEQC to update the Solution composition + # the .species_moles attribute should return MOLES (not moles per ___) + for s, mol in self.ppsol.species_moles.items(): + solution.components[s] = mol + + # log a message if any components were not touched by PHREEQC + # if that was the case, re-adjust the charge balance to account for those species (since PHREEQC did not) + missing_species = set(self._stored_comp.keys()) - {standardize_formula(s) for s in self.ppsol.species} + if len(missing_species) > 0: + logger.warning( + f"After equilibration, the amounts of species {sorted(missing_species)} were not modified " + "by PHREEQC. These species are likely absent from its database." + ) + + # tolerance (in moles) for detecting cases where an element amount + # is no longer balanced because of species that are not recognized + # by PHREEQC. + _atol = 1e-16 + + new_el_dict = solution.get_el_amt_dict(nested=True) + for el in orig_el_dict: + orig_el_amount = sum([orig_el_dict[el][k] for k in orig_el_dict[el]]) + new_el_amount = sum([new_el_dict[el][k] for k in new_el_dict.get(el, [])]) + + # If this element went "missing", add back all components that + # contain this element (for any valence value) + if orig_el_amount - new_el_amount > _atol: + logger.info( + f"PHREEQC discarded element {el} during equilibration. Adding all components for this element." + ) + solution.components.update( + { + component: self._stored_comp[component] + for components in orig_components_by_element[el].values() + for component in components + if component not in solution.components + } + ) + + # re-adjust charge balance for any missing species + # note that if balance_charge is set, it will have been passed to PHREEQC, so the only reason to re-adjust charge balance here is to account for any missing species. + solution._adjust_charge_balance() + + # rescale the solvent mass to ensure the total mass of solution does not change + # this is important because PHREEQC and the pyEQL database may use slightly different molecular + # weights for water. Since water amount is passed to PHREEQC in kg but returned in moles, each + # call to equilibrate can thus result in a slight change in the Solution mass. + solution.components[solution.solvent] = orig_solvent_moles diff --git a/src/pyEQL/pyEQL-phreeqc/pyEQL_phreeqc/core.py b/src/pyEQL/pyEQL-phreeqc/pyEQL_phreeqc/core.py index 4ad47ff0..589803a0 100644 --- a/src/pyEQL/pyEQL-phreeqc/pyEQL_phreeqc/core.py +++ b/src/pyEQL/pyEQL-phreeqc/pyEQL_phreeqc/core.py @@ -2,11 +2,20 @@ from pathlib import Path from textwrap import dedent, indent from typing import Any +from weakref import ref from pyEQL_phreeqc._bindings import PyIPhreeqc from pyEQL_phreeqc.solution import Solution from pyEQL_phreeqc.var import Var +SOLUTION_PROPS = ( + "CELL_NO", + "TOT['water']", + "OSMOTIC", +) +SPECIES_PROPS = ("MOL", "ACT") +EQ_SPECIES_PROPS = ("SI",) + class Phreeqc: def __init__(self, database: str = "phreeqc.dat", database_directory: Path | None = None): @@ -56,7 +65,9 @@ def _add_solution(self, solution: Solution | list[Solution]) -> Solution | list[ for solution in solutions: index = len(self) + solution._phreeqc = ref(self) solution._number = index + # TODO: This should go in the Solution class template = ( "\n" @@ -86,19 +97,10 @@ def remove_solution(self, index: int) -> Solution: self() return self._solutions.pop(index) - def add_solution( - self, - solution: Solution | list[Solution], - solution_props: tuple[str] | None = None, - species_props: tuple[str] | None = None, - ) -> Solution | list[Solution]: - if solution_props is None: - solution_props = ("OSMOTIC",) - solution_punch_line = ", ".join(list(solution_props)) - - if species_props is None: - species_props = ("MOL", "ACT") - species_punch_line = ", ".join([f"{prop}(name$(i))" for prop in species_props]) + def add_solution(self, solution: Solution | list[Solution]) -> Solution | list[Solution]: + solution_punch_line = ", ".join(list(SOLUTION_PROPS)) + species_punch_line = ", ".join([f"{prop}(name$(i))" for prop in SPECIES_PROPS]) + eq_species_punch_line = ", ".join([f"{prop}(name$(j))" for prop in EQ_SPECIES_PROPS]) self.clear() self.accumulate(f""" @@ -111,38 +113,71 @@ def add_solution( 20 FOR i = 1 to count 30 PUNCH name$(i), {species_punch_line} 40 NEXT i + 50 PUNCH EOL$ + 60 p = SYS("phases", count, name$, type$, moles) + 70 FOR j = 1 TO count + 80 PUNCH name$(j), {eq_species_punch_line} + 90 NEXT j """) return_value = self._add_solution(solution) - self() + self._parse_output() + + return return_value + def _parse_output(self): # first line is always the header, but ensure this so we can skip it. assert self.output[0][0].startswith("no_heading") output = self.output[:][1:] - for solution_i, line in enumerate(output): - solution_i_props = {"species": {}} + for line in output: + solution_i_props = {"species": {}, "eq_species": {}} n_tokens = len(line) j = 0 - while j < n_tokens: - # everything before the first '\n' are solution_props - if not solution_i_props["species"]: - while (solution_prop := line[j]) != "\n": - solution_i_props[solution_props[j]] = solution_prop - j += 1 - j += 1 # skip "\n" + while j < len(SOLUTION_PROPS): + solution_i_props[SOLUTION_PROPS[j]] = line[j] + j += 1 + + j += 1 # skip "\n" + + while line[j] != "\n": species = line[j] solution_i_props["species"][species] = {} j += 1 - for prop in species_props: + for prop in SPECIES_PROPS: solution_i_props["species"][species][prop] = line[j] j += 1 + j += 1 # skip "\n" + while j < n_tokens: + # We encounter None values here, indicating end of valid + # entries. + if line[j] is None: + break + eq_species = line[j] + solution_i_props["eq_species"][eq_species] = {} + j += 1 + for prop in EQ_SPECIES_PROPS: + solution_i_props["eq_species"][eq_species][prop] = line[j] + j += 1 + + solution_i = int(solution_i_props["CELL_NO"]) self[solution_i]._set_calculated_props(solution_i_props) - return return_value + def equalize(self, index: int, phases: list[str], saturation_indices: list[float], amounts: list[float]) -> None: + phases_lines = "\n".join( + f"{p} {si} {amount}" for (p, si, amount) in zip(phases, saturation_indices, amounts, strict=False) + ) + self.accumulate(f""" + USE SOLUTION {index} + EQUILIBRIUM PHASES 1 + {phases_lines} + END + """) + self() + self._parse_output() class PhreeqcOutput: diff --git a/src/pyEQL/pyEQL-phreeqc/pyEQL_phreeqc/solution.py b/src/pyEQL/pyEQL-phreeqc/pyEQL_phreeqc/solution.py index 8c7e8f3e..0e8f06bd 100644 --- a/src/pyEQL/pyEQL-phreeqc/pyEQL_phreeqc/solution.py +++ b/src/pyEQL/pyEQL-phreeqc/pyEQL_phreeqc/solution.py @@ -3,9 +3,10 @@ class Solution: - def __init__(self, props): + def __init__(self, props, phreeqc=None): + self._phreeqc = phreeqc self._number = -1 - self._input_props = deepcopy(props) + self._input_props = props self._calculated_props = {} def __str__(self): @@ -17,16 +18,54 @@ def _set_calculated_props(self, props: dict[Any, Any]): def _get_calculated_props(self): return self._calculated_props - def _get_calculated_prop(self, which, species: str | None = None): + def _get_calculated_prop(self, which, species: str | None = None, eq_species: str | None = None): if species is not None: return self._calculated_props["species"][species][which] + if eq_species is not None: + return self._calculated_props["eq_species"][eq_species][which] return self._calculated_props[which] - def get_activity(self, species): - return self._get_calculated_prop("ACT", species) + def get_activity(self, species) -> float: + return self._get_calculated_prop("ACT", species=species) - def get_molality(self, species): - return self._get_calculated_prop("MOL", species) + def get_molality(self, species) -> float: + return self._get_calculated_prop("MOL", species=species) - def get_osmotic_coefficient(self): + def get_osmotic_coefficient(self) -> float: return self._get_calculated_prop("OSMOTIC") + + def get_species_list(self) -> list[str]: + return list(self._calculated_props["species"]) + + def get_kgw(self) -> float: + return self._get_calculated_prop("TOT['water']") + + def get_moles(self, species) -> float: + return self.get_molality(species) * self.get_kgw() + + def equalize(self, phases: list[str], saturation_indices: list[float], amounts: list[float]) -> None: + if self._phreeqc is not None: + self._phreeqc().equalize(self._number, phases, saturation_indices, amounts) + + def si(self, eq_species) -> float: + return self._get_calculated_prop("SI", eq_species=eq_species) + + """ + The following properties are somewhat redundant, but included in here + so we can act as a drop-in replacement for PhreeqPython as far as its + usage in pyEQL. + """ + + @property + def species(self) -> list[str]: + return self.get_species_list() + + @property + def species_moles(self) -> dict[str, float]: + amounts = {species: self.get_moles(species) for species in self.species} + # filtering only till issue (TODO) is resolved + return {k: v for k, v in amounts.items() if k != "(CO2)2"} + + def forget(self) -> None: + if self._phreeqc is not None: + self._phreeqc().remove_solution(self._number) diff --git a/tests/phreeqc/test_phreeqc.py b/tests/phreeqc/test_phreeqc.py index 5e8c1dc1..2cba7b06 100644 --- a/tests/phreeqc/test_phreeqc.py +++ b/tests/phreeqc/test_phreeqc.py @@ -403,11 +403,16 @@ def test_add_solution_input(): -reset false USER_PUNCH - 5 PUNCH OSMOTIC, EOL_NOTAB$ + 5 PUNCH CELL_NO, TOT['water'], OSMOTIC, EOL_NOTAB$ 10 t = SYS("aq", count, name$, type$, moles) 20 FOR i = 1 to count 30 PUNCH name$(i), MOL(name$(i)), ACT(name$(i)) 40 NEXT i + 50 PUNCH EOL$ + 60 p = SYS("phases", count, name$, type$, moles) + 70 FOR j = 1 TO count + 80 PUNCH name$(j), SI(name$(j)) + 90 NEXT j SOLUTION 0 pH 7.0 @@ -427,8 +432,7 @@ def test_add_solution_input(): units mol/kgw water 0.9970480319717386 SAVE SOLUTION 1 - END - """) + END""") assert str(phreeqc) == expected @@ -449,11 +453,70 @@ def test_add_solution_output(): # heading + soln 0 data + soln 1 data (regardless of whether we have -headings in USER_PUNCH or not) assert phreeqc.get_selected_output_row_count() == 3 - # , "\n", +[, , ] repeated for each (4) species - assert phreeqc.get_selected_output_column_count() == 14 + # , , , "\n", + # +[, , ] repeated for each (4) species + # +"\n" + [, ] repeated for each (3) equilibrium species + assert phreeqc.get_selected_output_column_count() == 23 + + +def test_kgw(): + phreeqc = Phreeqc() + phreeqc.add_solution( + Solution({"pH": 7.0, "pe": 8.5, "redox": "pe", "temp": 25.0, "units": "mol/kgw", "water": 0.9970480319717386}) + ) + assert phreeqc[0].get_kgw() == approx(0.9970480319717386) def test_speciate(): + phreeqc = Phreeqc() + phreeqc.add_solution( + Solution({"pH": 7.0, "pe": 8.5, "redox": "pe", "temp": 25.0, "units": "mol/kgw", "water": 0.9970480319717386}) + ) + assert phreeqc[0].get_species_list() == ["OH-", "H+", "O2", "H2"] + + +def test_speciate_get_molality(): + solution = Solution( + {"pH": 7.0, "pe": 8.5, "redox": "pe", "temp": 25.0, "units": "mol/kgw", "water": 0.9970480319717386} + ) + phreeqc = Phreeqc() + phreeqc.add_solution(solution) + + activity = phreeqc[0].get_molality("H+") + assert activity == approx(1.0005246407839175e-07) + + +def test_speciate_get_moles(): + solution = Solution( + {"pH": 7.0, "pe": 8.5, "redox": "pe", "temp": 25.0, "units": "mol/kgw", "water": 0.9970480319717386} + ) + phreeqc = Phreeqc() + phreeqc.add_solution(solution) + + activity = phreeqc[0].get_moles("H+") + assert activity == approx(1.0005246407839175e-07 * 0.9970480319717386) + + +def test_speciate_species_moles(): + solution = Solution( + {"pH": 7.0, "pe": 8.5, "redox": "pe", "temp": 25.0, "units": "mol/kgw", "water": 0.9970480319717386} + ) + phreeqc = Phreeqc() + phreeqc.add_solution(solution) + + species_moles = phreeqc[0].species_moles + expected = { + "H+": 9.975711240328358e-08, + "H2": 7.05855918557026e-35, + "O2": 8.293083928522462e-25, + "OH-": 1.009700360204178e-07, + } + assert set(species_moles.keys()) == set(expected.keys()) + for k, v in species_moles.items(): + assert v == approx(expected[k]) + + +def test_species_all_props(): solutions = [ Solution({"pH": 7.0, "pe": 8.5, "redox": "pe", "temp": 25.0, "units": "mol/kgw", "water": 0.9970480319717386}), Solution({"pH": 10.0, "pe": 8.5, "redox": "pe", "temp": 50.0, "units": "mol/kgw", "water": 0.9970480319717386}), @@ -463,6 +526,8 @@ def test_speciate(): expected = { 0: { + "CELL_NO": 0, + "TOT['water']": 0.9970480319717386, "OSMOTIC": 0.0, "species": { "H+": {"ACT": 1.0001522689856982e-07, "MOL": 1.0005246407839175e-07}, @@ -470,8 +535,15 @@ def test_speciate(): "O2": {"ACT": 8.317637520771417e-25, "MOL": 8.317637327985348e-25}, "OH-": {"ACT": 1.0123126727760366e-07, "MOL": 1.0126897880811404e-07}, }, + "eq_species": { + "H2(g)": {"SI": -31.048923033678367}, + "H2O(g)": {"SI": -1.5028233204748613}, + "O2(g)": {"SI": -21.18764110890733}, + }, }, 1: { + "CELL_NO": 1, + "TOT['water']": 0.9970480319717386, "OSMOTIC": 0.0, "species": { "H+": {"ACT": 1e-10, "MOL": 1.0197827284798617e-10}, @@ -479,6 +551,11 @@ def test_speciate(): "O2": {"ACT": 3.659468532125681e-05, "MOL": 3.6592332322612637e-05}, "OH-": {"ACT": 0.0005473549298676792, "MOL": 0.0005585111533405827}, }, + "eq_species": { + "H2(g)": {"SI": -37.113374354647966}, + "H2O(g)": {"SI": -0.9158658033784839}, + "O2(g)": {"SI": -1.4065718231105961}, + }, }, } @@ -487,26 +564,19 @@ def test_speciate(): props = phreeqc[solution_index]._get_calculated_props() assert set(props.keys()) == set(expected[solution_index].keys()) for k in expected[solution_index]: - if k != "species": + if k not in ("species", "eq_species"): assert props[k] == approx(expected[solution_index][k]) for k in solution_props["species"]: for prop in ("ACT", "MOL"): assert np.isclose(props["species"][k][prop], expected[solution_index]["species"][k][prop]) - -def test_speciate_molality(): - solution = Solution( - {"pH": 7.0, "pe": 8.5, "redox": "pe", "temp": 25.0, "units": "mol/kgw", "water": 0.9970480319717386} - ) - phreeqc = Phreeqc() - phreeqc.add_solution(solution) - - activity = phreeqc[0].get_molality("H+") - assert activity == approx(1.0005246407839175e-07) + for k in solution_props["eq_species"]: + for prop in ("SI",): + assert np.isclose(props["eq_species"][k][prop], expected[solution_index]["eq_species"][k][prop]) -def test_speciate_activity(): +def test_speciate_get_activity(): solution = Solution( {"pH": 7.0, "pe": 8.5, "redox": "pe", "temp": 25.0, "units": "mol/kgw", "water": 0.9970480319717386} ) @@ -517,7 +587,7 @@ def test_speciate_activity(): assert activity == approx(1.0001522689856982e-07) -def test_osmotic_coefficient(): +def test_get_osmotic_coefficient(): solution = Solution( {"pH": 7.0, "pe": 8.5, "redox": "pe", "temp": 25.0, "units": "mol/kgw", "water": 0.9970480319717386} ) diff --git a/tests/test_pyeql_phreeqc.py b/tests/test_pyeql_phreeqc.py new file mode 100644 index 00000000..a82af318 --- /dev/null +++ b/tests/test_pyeql_phreeqc.py @@ -0,0 +1,414 @@ +""" +pyEQL volume and concentration methods test suite +================================================= + +This file contains tests for the volume and concentration-related methods +used by pyEQL's Solution class +""" + +import logging + +import numpy as np +import pytest + +from pyEQL import Solution +from pyEQL.engines import PyEQLEOS + + +@pytest.fixture +def s1(): + return Solution(volume="2 L", engine="pyeql") + + +@pytest.fixture +def s2(): + return Solution([["Na+", "4 mol/L"], ["Cl-", "4 mol/L"]], volume="2 L", engine="pyeql") + + +@pytest.fixture +def s3(): + return Solution([["Na+", "4 mol/kg"], ["Cl-", "4 mol/kg"]], volume="2 L", engine="pyeql") + + +@pytest.fixture +def s5(): + # 100 mg/L as CaCO3 ~ 1 mM + return Solution([["Ca+2", "40.078 mg/L"], ["CO3-2", "60.0089 mg/L"]], volume="1 L", engine="pyeql") + + +@pytest.fixture +def s5_pH(): + # 100 mg/L as CaCO3 ~ 1 mM + return Solution( + [["Ca+2", "40.078 mg/L"], ["CO3-2", "60.0089 mg/L"]], volume="1 L", balance_charge="pH", engine="pyeql" + ) + + +@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", + engine="pyeql", + ) + + +@pytest.fixture +def s6_Ca(): + # 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", + balance_charge="Ca+2", + engine="pyeql", + ) + + +def test_empty_solution_3(): + # create an empty solution + s1 = Solution(database=None, engine="pyeql") + # It should return type Solution + assert isinstance(s1, Solution) + # It should have exactly 1L volume + assert s1.volume.to("L").magnitude == 1.0 + # the solvent should be water + assert s1.solvent == "H2O(aq)" + # It should have 0.997 kg water mass + assert np.isclose(s1.solvent_mass.to("kg").magnitude, 0.997, atol=1e-3) + # the temperature should be 25 degC + assert s1.temperature.to("degC").magnitude == 25 + # the pressure should be 1 atm + assert s1.pressure.to("atm").magnitude == 1 + # the pH should be 7.0 + assert np.isclose(s1.get_activity("H+"), 1e-7, atol=1e-9) + # assert np.isclose(s1.pH, 7.0, atol=0.01) + assert np.isclose(s1.pE, 8.5) + # it should contain H2O, H+, and OH- species + assert set(s1.components.keys()) == {"H2O(aq)", "OH[-1]", "H[+1]"} + + +def test_init_engines(): + """ + Test passing an EOS instance as well as the ideal and native EOS + """ + s = Solution([["Na+", "4 mol/L"], ["Cl-", "4 mol/L"]], engine="pyeql") + assert isinstance(s.engine, PyEQLEOS) + assert s.get_activity_coefficient("Na+").magnitude * s.get_activity_coefficient("Cl-").magnitude < 1 + assert s.get_osmotic_coefficient().magnitude == 1 + # with pytest.warns(match="Solute Mg+2 not found"): + assert s.get_activity_coefficient("Mg+2").magnitude == 1 + assert s.get_activity("Mg+2").magnitude == 0 + s.engine._destroy_ppsol() + assert s.engine.ppsol is None + + +def test_conductivity(s1): + # even an empty solution should have some conductivity + assert s1.conductivity > 0 + + for conc, cond in zip([0.001, 0.05, 0.1], [123.68, 111.01, 106.69], strict=False): + s1 = Solution({"Na+": f"{conc} mol/L", "Cl-": f"{conc} mol/L"}) + assert np.isclose(s1.conductivity.to("S/m").magnitude, conc * cond / 10, atol=0.5), ( + f"Conductivity test failed for NaCl at {conc} mol/L. Result = {s1.conductivity.to('S/m').magnitude}" + ) + + # higher concentration data points from Appelo, 2017 Figure 4. + s1 = Solution({"Na+": "2 mol/kg", "Cl-": "2 mol/kg"}) + assert np.isclose(s1.conductivity.to("mS/cm").magnitude, 145, atol=10) + + # MgCl2 + for conc, cond in zip([0.001, 0.05, 0.1], [124.15, 114.49, 97.05], strict=False): + s1 = Solution({"Mg+2": f"{conc} mol/L", "Cl-": f"{2 * conc} mol/L"}) + assert np.isclose(s1.conductivity.to("S/m").magnitude, 2 * conc * cond / 10, atol=1), ( + f"Conductivity test failed for MgCl2 at {conc} mol/L. Result = {s1.conductivity.to('S/m').magnitude}" + ) + + # per CRC handbook "standard KCl solutions for calibrating conductiVity cells", + # 0.1m KCl has a conductivity of 12.824 mS/cm at 25 C + s_kcl = Solution({"K+": "0.1 mol/kg", "Cl-": "0.1 mol/kg"}) + assert np.isclose(s_kcl.conductivity.magnitude, 1.2824, atol=0.02) # conductivity is in S/m + + s_kcl.temperature = "5 degC" + assert np.isclose(s_kcl.conductivity.magnitude, 0.81837, atol=0.06) + + s_kcl.temperature = "50 degC" + assert np.isclose(s_kcl.conductivity.magnitude, 1.91809, atol=0.18) + + # TODO - conductivity model not very accurate at high conc. + s_kcl = Solution({"K+": "1 mol/kg", "Cl-": "1 mol/kg"}) + assert np.isclose(s_kcl.conductivity.magnitude, 10.862, rtol=0.05) + + +def test_equilibrate(s1, s2, s5_pH, s6_Ca, caplog): + assert "H2(aq)" not in s1.components + orig_pH = s1.pH + orig_pE = s1.pE + orig_mass = s1.mass + orig_density = s2.density.magnitude + orig_solv_mass = s2.solvent_mass.magnitude + s1.equilibrate() + assert "H2(aq)" in s1.components + assert np.isclose(s1.charge_balance, 0, atol=1e-8) + assert np.isclose(s1.pH, orig_pH, atol=0.01) + assert np.isclose(s1.pE, orig_pE) + + assert "NaOH(aq)" not in s2.components + orig_density = s2.density.magnitude + orig_solv_mass = s2.solvent_mass.magnitude + orig_pH = s2.pH + orig_pE = s2.pE + orig_mass = s2.mass + s2.equilibrate() + assert np.isclose(s2.mass, orig_mass) + assert np.isclose(s2.density.magnitude, orig_density) + assert np.isclose(s2.solvent_mass.magnitude, orig_solv_mass) + # assert "NaOH(aq)" in s2.components + + # total element concentrations should be conserved after equilibrating + assert np.isclose(s2.get_total_amount("Na", "mol").magnitude, 8) + assert np.isclose(s2.get_total_amount("Cl", "mol").magnitude, 8) + assert np.isclose(s2.solvent_mass.magnitude, orig_solv_mass) + assert np.isclose(s2.density.magnitude, orig_density) + # the pH should drop due to hydrolysis of Ca+2 + assert s2.pH < orig_pH + assert np.isclose(s2.pE, orig_pE) + assert np.isclose(s2.mass, orig_mass) + # this solution has balance_charge=None, therefore, the charge balance + # may be off after equilibration + assert not np.isclose(s2.charge_balance, 0, atol=1e-8) + eq_Hplus = s2.components["H+"] + s2.balance_charge = "pH" + s2.equilibrate() + assert np.isclose(s2.charge_balance, 0, atol=1e-8) + assert s2.components["H+"] > eq_Hplus + + # test log message if there is a species not present in the phreeqc database + s_zr = Solution( + {"Zr+4": "0.05 mol/kg", "Na+": "0.05 mol/kg", "Cl-": "0.1 mol/kg"}, engine="pyeql", log_level="WARNING" + ) + totzr = s_zr.get_total_amount("Zr", "mol") + with caplog.at_level(logging.WARNING, "pyEQL.engines"): + s_zr.equilibrate() + assert "likely absent from its database" in caplog.text + assert "Zr[+4]" in s_zr.components + assert s_zr.get_total_amount("Zr", "mol") == totzr + + # this solution is the only one in the test that contains alkalinity + # and equilibrating it results in a shift in the pH + # the CO3-2 initially present reacts with the water to consume H+ and + # increase the pH by approximately 0.0006 M (b/c at pH 7 virtually all + # carbonate is present as HCO3-) -log10(0.001) = + assert "HCO3[-1]" not in s5_pH.components + assert np.isclose(s5_pH.charge_balance, 0) + orig_pH = s5_pH.pH + orig_pE = s5_pH.pE + orig_density = s5_pH.density.magnitude + orig_solv_mass = s5_pH.solvent_mass.magnitude + set(s5_pH.components.keys()) + s5_pH.equilibrate() + assert np.isclose(s5_pH.get_total_amount("Ca", "mol").magnitude, 0.001) + assert np.isclose(s5_pH.get_total_amount("C(4)", "mol").magnitude, 0.001, atol=1e-7) + # due to the large pH shift, water mass and density need not be perfectly conserved + assert np.isclose(s5_pH.solvent_mass.magnitude, orig_solv_mass, atol=1e-3) + assert np.isclose(s5_pH.density.magnitude, orig_density, atol=1e-3) + assert np.isclose(s5_pH.charge_balance, 0) + assert "CaOH[+1]" in s5_pH.components + assert "HCO3[-1]" in s5_pH.components + s5_pH_after = s5_pH.pH + assert s5_pH_after > orig_pH + assert np.isclose(s5_pH.pE, orig_pE) + + # repeated calls to equilibrate should not change the properties (much) + for i in range(10): + s5_pH.equilibrate() + assert np.isclose(s5_pH.charge_balance, 0, atol=1e-8), f"C.B. failed at iteration {i}" + assert np.isclose(s5_pH.pH, s5_pH_after, atol=0.01), f"pH failed at iteration {i}" + assert np.isclose(s5_pH.pE, orig_pE), f"pE failed at iteration {i}" + + # test equilibrate() with a non-pH balancing species + assert np.isclose(s6_Ca.charge_balance, 0, atol=1e-8) + initial_Ca = s6_Ca.get_total_amount("Ca", "mol").magnitude + assert s6_Ca.balance_charge == "Ca[+2]" + s6_Ca.equilibrate() + assert s6_Ca.get_total_amount("Ca", "mol").magnitude != initial_Ca + assert np.isclose(s6_Ca.charge_balance, 0, atol=1e-8) + + +def test_equilibrate_water_pH7(): + solution = Solution({}, pH=7.00, temperature="25 degC", volume="1 L", engine="pyeql") + solution.equilibrate() + # pH = -log10[H+] + assert np.isclose(solution.get_amount("H+", "mol/kg").magnitude, 1.001e-07) + # 14 - pH = log10[OH-] + assert np.isclose(solution.get_amount("OH-", "mol/kg").magnitude, 1.013e-07) + # small amount of H2 gas + assert solution.get_amount("H2", "mol/kg") > 0 + # Density of H2O in phreeqc = 0.99704 kg/L + # For 0.99704 kg H2O + # mol = 0.99704 kg * 1000 g/kg / (18.01528 g/mol) = 55.34413009400909 + # mol/kg = 55.34413009400909 / 0.99704 = 55.50843506179199 + assert np.isclose(solution.get_amount("H2O", "mol/kg").magnitude, 55.50843506179199) + + +def test_equilibrate_CO2_with_calcite(): + solution = Solution({}, pH=7.0, volume="1 L", engine="pyeql") + solution.equilibrate(atmosphere=False, gases={"CO2": -2.95, "O2": -0.6778}, solids=["Calcite"]) + # 5 rxns: I) CaCO3 dissolution, II) Ka1, III) Ka2, IV) water dissociation, V) CaHCO3+ rxn in PHREEQC + # 9 species, 5 components, 4 rxns exclude water dissociation + assert solution.get_amount("Na+", "mol") == 0 + assert np.isclose(solution.get_amount("CO2(aq)", "mol/kg").magnitude, 3.816e-05) + assert np.isclose( + solution.get_amount("HCO3-", "mol/kg").magnitude, 1.48e-03, atol=1e-5 + ) # slight tolerance adjustment + assert np.isclose( + solution.get_amount("Ca+2", "mol/kg").magnitude, 7.427e-04, atol=1e-5 + ) # slight tolerance adjustment + + +def test_equilibrate_FeO3H3_ppt(): + """Test an oversaturated solution""" + solution = Solution({"Fe+3": "0.01 mol/L", "OH-": "10**-7 mol/L"}, volume="1 L", engine="pyeql") + solution.equilibrate() + assert np.isclose(solution.get_amount("Fe+3", "mol/L").magnitude, 3.093e-11) + assert np.isclose(solution.get_amount("OH-", "mol/L").magnitude, 1.067e-07) + # The following assert passes, but we need to explain why phreeqc gives a + # Fe(OH)3(a) SI value of ~5.408 instead of the expected ~7.2 + # SI_FeO3H3 = np.log10((Fe_3) * (OH_) ** 3 / 10**-38.8) + assert solution.engine.ppsol.si("Fe(OH)3(a)") > 0 + + +def test_equilibrate_nophaseeq(): + # Test to see that equilibrating without any solids/gases has no effect + # on the concentrations. + solution0 = Solution({"CO2(aq)": "0.001 mol/L"}, pH=3.0, volume="1 L", engine="pyeql") + solution0.equilibrate() + solution1 = Solution({"CO2(aq)": "0.001 mol/L"}, pH=3.0, volume="1 L", engine="pyeql") + solution1.equilibrate(atmosphere=False, solids=None, gases=None) + + assert np.isclose(solution0.get_amount("CO2(aq)", "mol"), solution1.get_amount("CO2(aq)", "mol")) + assert np.isclose(solution0.get_amount("HCO3-", "mol"), solution1.get_amount("HCO3-", "mol")) + assert np.isclose(solution0.get_amount("CO3-2", "mol"), solution1.get_amount("CO3-2", "mol")) + + +def test_equilibrate_logC_pH_carbonate_8_3(): + solution = Solution({"CO2(aq)": "0.001 mol/L"}, pH=8.3, volume="1 L", engine="pyeql") + solution.equilibrate() + + # To evaluate CO2 equilibrium using Henry's law + # Note: In the asserts below, the results were calculated assuming Henry's Law constant for CO2 as 0.034, which + # is slightly different from what phreeqc uses (logK = -1.468, which gives us 0.0340408), + # Hence the relaxed atols. + + # CO2 + HCO3- + CO3-2 = total C = 0.001 mol with pKa1 = 6.3 and pKa2 = 10.3 + # [H2CO3] + 100[H2CO3] + 2.5119e-3[H2CO3] = total C = 0.001 mol + assert np.isclose(solution.get_amount("CO2(aq)", "mol/kg").magnitude, 1.079e-5, atol=1e-5) + # The atols below are more relaxed because of the error carried forward from CO2(aq) + # HCO3- approx CO2(aq) dominant species at pH 7 + # Note: The pKa1 value used in the calculation below is 6.3. Phreeqc uses 6.352 + assert np.isclose(solution.get_amount("HCO3-", "mol/kg").magnitude, 9.8219e-4, atol=1e-5) + + # CO2 + HCO3- + CO3-2 = total C = 0.001 mol with pKa1 = 6.3 and pKa2 = 10.3 + # [H2CO3] + 100[H2CO3] + 5.0119e-9[H2CO3] = total C = 0.001 mol + # Note: The pKa1 value used in the calculation below is 6.3. Phreeqc uses 6.352 + # Note: The pKa2 value used in the calculation below is 10.3. Phreeqc uses 10.329 + assert np.isclose(solution.get_amount("CO3-2", "mol/kg").magnitude, 9.923e-6, atol=1e-5) + + +def test_equilibrate_logC_pH_carbonate_13(): + solution = Solution({"CO2(aq)": "0.001 mol/L"}, pH=13.0, volume="1 L", engine="pyeql") + solution.equilibrate() + + # To evaluate CO2 equilibrium using Henry's law + # Note: In the asserts below, the results were calculated assuming Henry's Law constant for CO2 as 0.034, which + # is slightly different from what phreeqc uses (logK = -1.468, which gives us 0.0340408), + # Hence the relaxed atols. + + # H2CO3 approx CO2(aq) is negligible at pH 13, but still > 0 + assert solution.get_amount("CO2(aq)", "mol/kg") > 0 + + # CO2 + HCO3- + CO3-2 = total C = 0.001 mol with pKa1 = 6.3 and pKa2 = 10.3 + # 1.995e-7[HCO3-] + [HCO3-] + (10^-10.3 / 10^-13)[CO3-2] = total C = 0.001 mol + # Note: The pKa1 value used in the calculation below is 6.3. Phreeqc uses 6.352 + assert np.isclose(solution.get_amount("HCO3-", "mol/kg").magnitude, 1.152e-6, atol=1e-5) + + # CO2 + HCO3- + CO3-2 = total C = 0.001 mol with pKa1 = 6.3 and pKa2 = 10.3 + # (10^-3*2 / 10^-16.6)[CO3-2] + (10^-3 / 10^-10.3)[CO3-2] + [CO3-2] = total C = 0.001 mol + # Note: The pKa1 value used in the calculation below is 6.3. Phreeqc uses 6.352 + # Note: The pKa2 value used in the calculation below is 10.3. Phreeqc uses 10.329 + assert np.isclose(solution.get_amount("CO3-2", "mol/kg").magnitude, 9.979e-4, atol=1e-5) + + +@pytest.mark.xfail(strict=True, reason="alkalinity discrepancy needs to be investigated") +def test_alkalinity(): + solution = Solution({"CO2(aq)": "0.001 mol/L"}, pH=7, volume="1 L", engine="pyeql") + solution.equilibrate() + + HCO3 = solution.get_amount("HCO3-", "mg/L").magnitude + CO3 = solution.get_amount("CO3-2", "mg/L").magnitude + OH = solution.get_amount("OH-", "mg/L").magnitude + H = solution.get_amount("H+", "mg/L").magnitude + + # Alkalinity calculated from the excess of negative charges from weak acids + calculated_alk = HCO3 + 2 * CO3 + OH - H + assert solution.alkalinity.to("mg/L").magnitude == pytest.approx(calculated_alk, abs=0.001) + + +def test_equilibrate_2L(): + solution = Solution({"Cu+2": "1 umol/L", "O-2": "1 umol/L"}, volume="2 L", engine="pyeql") + solution.equilibrate(atmosphere=True) + assert np.isclose(solution.get_total_amount("Cu", "umol").magnitude, 1.9999998687384424) + + +def test_equilibrate_unrecognized_component(): + solution = Solution({}, engine="pyeql") + # Specifying an unrecognized solid raises a ValueError + with pytest.raises(ValueError, match="Phase not found in database, Ferroxite."): + solution.equilibrate(solids=["Ferroxite"]) + + +def test_equilibrate_OER_region(): + # The combination of pH and pE values don't fall within the water stability region. + solution = Solution({}, pH=12.0, pE=13, volume="1 L", engine="pyeql") + with pytest.raises(ValueError, match="Activity of water has not converged."): + solution.equilibrate() + + +def test_equilibrate_gas_units(): + # Specify CO2 partial pressure directly as log10 partial pressure, as well + # as an explicit pressure unit. + # Note: log10(0.000316) = -3.5 + s0 = Solution({}, pH=7.0, volume="1 L", engine="pyeql") + s0.equilibrate(atmosphere=True, gases={"CO2": "0.00031622776601683794 atm"}) + s1 = Solution({}, pH=7.0, volume="1 L", engine="pyeql") + s1.equilibrate(atmosphere=True, gases={"CO2": -3.5}) + assert s0.components == s1.components + + +def test_equilibrate_with_atm(): + s1 = Solution({}, pH=7.0, volume="1 L", engine="pyeql") + s1.equilibrate(atmosphere=True) + # PHREEQCUI final CO2, O2, and N2 concentrations were slightly adjusted for consistency with wrapper outputs + assert np.isclose(s1.get_amount("CO2(aq)", "mol/L").magnitude, 1.397209520185121e-05) # PHREQCUI - 1.389e-05 + assert np.isclose(s1.get_amount("O2(aq)", "mol/L").magnitude, 0.00025982718533575387) # PHREEQCUI - 2.615e-04 + assert np.isclose(s1.get_amount("N2(aq)", "mol/L").magnitude, 0.0005043306329272451) # PHREEQCUI - 5.064e-04 From 2706d7f3b4658de163ee6d55279df9037a6cbbc2 Mon Sep 17 00:00:00 2001 From: Vineet Bansal Date: Tue, 13 Jan 2026 16:22:03 -0500 Subject: [PATCH 5/7] Changes to prep and test a release of the wrappers. --- src/pyEQL/engines.py | 5 ++--- src/pyEQL/pyEQL-phreeqc/pyproject.toml | 2 +- src/pyEQL/solution.py | 12 +++++++++++- tests/test_pyeql_phreeqc.py | 6 ++++-- 4 files changed, 18 insertions(+), 7 deletions(-) diff --git a/src/pyEQL/engines.py b/src/pyEQL/engines.py index 0329be5f..fc2cd09e 100644 --- a/src/pyEQL/engines.py +++ b/src/pyEQL/engines.py @@ -998,9 +998,8 @@ def get_activity_coefficient(self, solution: "solution.Solution", solute: str) - return ureg.Quantity(act, "dimensionless") def get_osmotic_coefficient(self, solution: "solution.Solution") -> ureg.Quantity: - osmotic = self.ppsol.get_osmotic_coefficient() # noqa: F841 - # TODO: Why is the wrapper defaulting to 1 instead of 0 (default behavior for phreeqc) for phreeqpython? - return ureg.Quantity(1, "dimensionless") + osmotic = self.ppsol.get_osmotic_coefficient() + return ureg.Quantity(osmotic, "dimensionless") def get_solute_volume(self, solution: "solution.Solution") -> ureg.Quantity: """Return the volume of the solutes.""" diff --git a/src/pyEQL/pyEQL-phreeqc/pyproject.toml b/src/pyEQL/pyEQL-phreeqc/pyproject.toml index 0bb07478..6eaa3eed 100644 --- a/src/pyEQL/pyEQL-phreeqc/pyproject.toml +++ b/src/pyEQL/pyEQL-phreeqc/pyproject.toml @@ -4,7 +4,7 @@ build-backend = "scikit_build_core.build" [project] name = "pyEQL-phreeqc" -version = "0.0.3" +version = "3.8.6.17100a1" requires-python = ">=3.10" [project.optional-dependencies] diff --git a/src/pyEQL/solution.py b/src/pyEQL/solution.py index 952c19e7..3941d3fa 100644 --- a/src/pyEQL/solution.py +++ b/src/pyEQL/solution.py @@ -54,7 +54,7 @@ def __init__( pE: float = 8.5, balance_charge: str | None = None, solvent: str | list = "H2O", - engine: EOS | Literal["native", "ideal", "phreeqc", "pyeql"] = "native", + engine: EOS | Literal["native", "ideal", "phreeqc", "pyeql"] | None = None, database: str | Path | Store | None = None, default_diffusion_coeff: float = 1.6106e-9, log_level: Literal["DEBUG", "INFO", "WARNING", "ERROR", "CRITICAL"] | None = "ERROR", @@ -237,6 +237,16 @@ def __init__( self.database.connect() self.logger.debug(f"Connected to property database {self.database!s}") + if engine is None: + engine = "native" + warnings.warn( + "The default value of the 'engine' parameter will change from " + "'native' to 'pyeql' in a future release. " + "Please specify 'engine' explicitly to silence this warning.", + DeprecationWarning, + stacklevel=2, + ) + # set the equation of state engine self._engine = engine # self.engine: Optional[EOS] = None diff --git a/tests/test_pyeql_phreeqc.py b/tests/test_pyeql_phreeqc.py index a82af318..9c4bcb81 100644 --- a/tests/test_pyeql_phreeqc.py +++ b/tests/test_pyeql_phreeqc.py @@ -115,7 +115,7 @@ def test_init_engines(): s = Solution([["Na+", "4 mol/L"], ["Cl-", "4 mol/L"]], engine="pyeql") assert isinstance(s.engine, PyEQLEOS) assert s.get_activity_coefficient("Na+").magnitude * s.get_activity_coefficient("Cl-").magnitude < 1 - assert s.get_osmotic_coefficient().magnitude == 1 + assert s.get_osmotic_coefficient().magnitude == 0 # with pytest.warns(match="Solute Mg+2 not found"): assert s.get_activity_coefficient("Mg+2").magnitude == 1 assert s.get_activity("Mg+2").magnitude == 0 @@ -183,7 +183,9 @@ def test_equilibrate(s1, s2, s5_pH, s6_Ca, caplog): assert np.isclose(s2.mass, orig_mass) assert np.isclose(s2.density.magnitude, orig_density) assert np.isclose(s2.solvent_mass.magnitude, orig_solv_mass) - # assert "NaOH(aq)" in s2.components + # Phreeqc 3.8 does not include NaOH log_k values, so instead of checking + # NaOH(aq), we check HCl(aq) instead. + assert "HCl(aq)" in s2.components # total element concentrations should be conserved after equilibrating assert np.isclose(s2.get_total_amount("Na", "mol").magnitude, 8) From 892a67666ee1a0f17150ac41a91ec3089e3e7244 Mon Sep 17 00:00:00 2001 From: Vineet Bansal Date: Tue, 13 Jan 2026 16:50:10 -0500 Subject: [PATCH 6/7] removed temp fix for a bug --- src/pyEQL/pyEQL-phreeqc/pyEQL_phreeqc/solution.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/pyEQL/pyEQL-phreeqc/pyEQL_phreeqc/solution.py b/src/pyEQL/pyEQL-phreeqc/pyEQL_phreeqc/solution.py index 0e8f06bd..ed5a6119 100644 --- a/src/pyEQL/pyEQL-phreeqc/pyEQL_phreeqc/solution.py +++ b/src/pyEQL/pyEQL-phreeqc/pyEQL_phreeqc/solution.py @@ -62,9 +62,7 @@ def species(self) -> list[str]: @property def species_moles(self) -> dict[str, float]: - amounts = {species: self.get_moles(species) for species in self.species} - # filtering only till issue (TODO) is resolved - return {k: v for k, v in amounts.items() if k != "(CO2)2"} + return {species: self.get_moles(species) for species in self.species} def forget(self) -> None: if self._phreeqc is not None: From ce31cc1b71e0f7f400682bfdd1e44331e8dfe8d8 Mon Sep 17 00:00:00 2001 From: Vineet Bansal Date: Thu, 15 Jan 2026 10:05:26 -0500 Subject: [PATCH 7/7] renamed Solution->PHRQSol to avoid confusion with pyEQL's Solution class; DeprecationWarning on any use of the native engine --- src/pyEQL/engines.py | 4 +-- .../pyEQL-phreeqc/pyEQL_phreeqc/__init__.py | 4 +-- src/pyEQL/pyEQL-phreeqc/pyEQL_phreeqc/core.py | 14 ++++---- .../pyEQL-phreeqc/pyEQL_phreeqc/solution.py | 2 +- src/pyEQL/solution.py | 13 ++++---- tests/phreeqc/test_phreeqc.py | 32 +++++++++---------- tests/phreeqc/test_phreeqc_solution.py | 4 +-- 7 files changed, 37 insertions(+), 36 deletions(-) diff --git a/src/pyEQL/engines.py b/src/pyEQL/engines.py index fc2cd09e..34a4b162 100644 --- a/src/pyEQL/engines.py +++ b/src/pyEQL/engines.py @@ -901,7 +901,7 @@ def __init__( def _setup_ppsol(self, solution: "solution.Solution") -> None: """Helper method to set up a PhreeqPython solution for subsequent analysis.""" - from pyEQL_phreeqc import Solution # noqa: PLC0415 + from pyEQL_phreeqc import PHRQSol # noqa: PLC0415 self._stored_comp = solution.components.copy() solv_mass = solution.solvent_mass.to("kg").magnitude @@ -957,7 +957,7 @@ def _setup_ppsol(self, solution: "solution.Solution") -> None: d[key] += " charge" try: - ppsol = self.pp.add_solution(Solution(d)) + ppsol = self.pp.add_solution(PHRQSol(d)) except Exception as e: # catch problems with the input to phreeqc raise ValueError( diff --git a/src/pyEQL/pyEQL-phreeqc/pyEQL_phreeqc/__init__.py b/src/pyEQL/pyEQL-phreeqc/pyEQL_phreeqc/__init__.py index 459a055d..7b8ad1e3 100644 --- a/src/pyEQL/pyEQL-phreeqc/pyEQL_phreeqc/__init__.py +++ b/src/pyEQL/pyEQL-phreeqc/pyEQL_phreeqc/__init__.py @@ -1,4 +1,4 @@ from .core import Phreeqc -from .solution import Solution +from .solution import PHRQSol -__all__ = ["Phreeqc", "Solution"] +__all__ = ["PHRQSol", "Phreeqc"] diff --git a/src/pyEQL/pyEQL-phreeqc/pyEQL_phreeqc/core.py b/src/pyEQL/pyEQL-phreeqc/pyEQL_phreeqc/core.py index 589803a0..763d4082 100644 --- a/src/pyEQL/pyEQL-phreeqc/pyEQL_phreeqc/core.py +++ b/src/pyEQL/pyEQL-phreeqc/pyEQL_phreeqc/core.py @@ -5,7 +5,7 @@ from weakref import ref from pyEQL_phreeqc._bindings import PyIPhreeqc -from pyEQL_phreeqc.solution import Solution +from pyEQL_phreeqc.solution import PHRQSol from pyEQL_phreeqc.var import Var SOLUTION_PROPS = ( @@ -26,7 +26,7 @@ def __init__(self, database: str = "phreeqc.dat", database_directory: Path | Non self._ext.load_database(str(database_directory / database)) self._str = "" - self._solutions: list[Solution] = [] + self._solutions: list[PHRQSol] = [] # TODO: Is VAR the common denominator for most operations? # Here we create one and modify it in operations instead of having @@ -59,8 +59,8 @@ def clear(self): def accumulate(self, s: str) -> None: self._str += dedent(s) - def _add_solution(self, solution: Solution | list[Solution]) -> Solution | list[Solution]: - singleton = isinstance(solution, Solution) + def _add_solution(self, solution: PHRQSol | list[PHRQSol]) -> PHRQSol | list[PHRQSol]: + singleton = isinstance(solution, PHRQSol) solutions = [solution] if singleton else solution for solution in solutions: @@ -68,7 +68,7 @@ def _add_solution(self, solution: Solution | list[Solution]) -> Solution | list[ solution._phreeqc = ref(self) solution._number = index - # TODO: This should go in the Solution class + # TODO: This should go in the PHRQSol class template = ( "\n" + cleandoc(f""" @@ -85,7 +85,7 @@ def _add_solution(self, solution: Solution | list[Solution]) -> Solution | list[ return self._solutions[-1] if singleton else self._solutions - def remove_solution(self, index: int) -> Solution: + def remove_solution(self, index: int) -> PHRQSol: _str = ( cleandoc(f""" DELETE @@ -97,7 +97,7 @@ def remove_solution(self, index: int) -> Solution: self() return self._solutions.pop(index) - def add_solution(self, solution: Solution | list[Solution]) -> Solution | list[Solution]: + def add_solution(self, solution: PHRQSol | list[PHRQSol]) -> PHRQSol | list[PHRQSol]: solution_punch_line = ", ".join(list(SOLUTION_PROPS)) species_punch_line = ", ".join([f"{prop}(name$(i))" for prop in SPECIES_PROPS]) eq_species_punch_line = ", ".join([f"{prop}(name$(j))" for prop in EQ_SPECIES_PROPS]) diff --git a/src/pyEQL/pyEQL-phreeqc/pyEQL_phreeqc/solution.py b/src/pyEQL/pyEQL-phreeqc/pyEQL_phreeqc/solution.py index ed5a6119..6a8d02ed 100644 --- a/src/pyEQL/pyEQL-phreeqc/pyEQL_phreeqc/solution.py +++ b/src/pyEQL/pyEQL-phreeqc/pyEQL_phreeqc/solution.py @@ -2,7 +2,7 @@ from typing import Any -class Solution: +class PHRQSol: def __init__(self, props, phreeqc=None): self._phreeqc = phreeqc self._number = -1 diff --git a/src/pyEQL/solution.py b/src/pyEQL/solution.py index 3941d3fa..92b3b522 100644 --- a/src/pyEQL/solution.py +++ b/src/pyEQL/solution.py @@ -54,7 +54,7 @@ def __init__( pE: float = 8.5, balance_charge: str | None = None, solvent: str | list = "H2O", - engine: EOS | Literal["native", "ideal", "phreeqc", "pyeql"] | None = None, + engine: EOS | Literal["native", "ideal", "phreeqc", "pyeql"] = "native", database: str | Path | Store | None = None, default_diffusion_coeff: float = 1.6106e-9, log_level: Literal["DEBUG", "INFO", "WARNING", "ERROR", "CRITICAL"] | None = "ERROR", @@ -237,12 +237,13 @@ def __init__( self.database.connect() self.logger.debug(f"Connected to property database {self.database!s}") - if engine is None: - engine = "native" + if engine == "native": warnings.warn( - "The default value of the 'engine' parameter will change from " - "'native' to 'pyeql' in a future release. " - "Please specify 'engine' explicitly to silence this warning.", + 'In the next release, the default engine ("native") will' + "transition to a new version of the PHREEQC wrapper for" + "speciation calculations. No change in your script is" + "required, but if you call .equilibrate(), compare results" + "carefully between releases.", DeprecationWarning, stacklevel=2, ) diff --git a/tests/phreeqc/test_phreeqc.py b/tests/phreeqc/test_phreeqc.py index 2cba7b06..df2da2f2 100644 --- a/tests/phreeqc/test_phreeqc.py +++ b/tests/phreeqc/test_phreeqc.py @@ -2,7 +2,7 @@ from pathlib import Path import numpy as np -from pyEQL_phreeqc import Phreeqc, Solution +from pyEQL_phreeqc import Phreeqc, PHRQSol from pytest import approx @@ -133,7 +133,7 @@ def test_run_simple(): def test_run_add_solution(): phreeqc = Phreeqc() - solution = Solution( + solution = PHRQSol( {"pH": 7.0, "pe": 8.5, "redox": "pe", "temp": 25.0, "units": "mol/kgw", "water": 0.9970480319717386} ) phreeqc.add_solution(solution) @@ -142,7 +142,7 @@ def test_run_add_solution(): def test_run_add_delete_solution(): phreeqc = Phreeqc() - solution = Solution( + solution = PHRQSol( {"pH": 7.0, "pe": 8.5, "redox": "pe", "temp": 25.0, "units": "mol/kgw", "water": 0.9970480319717386} ) phreeqc.add_solution(solution) @@ -389,10 +389,10 @@ def test_add_solution_input(): phreeqc.add_solution( [ - Solution( + PHRQSol( {"pH": 7.0, "pe": 8.5, "redox": "pe", "temp": 25.0, "units": "mol/kgw", "water": 0.9970480319717386} ), - Solution( + PHRQSol( {"pH": 10.0, "pe": 8.5, "redox": "pe", "temp": 50.0, "units": "mol/kgw", "water": 0.9970480319717386} ), ] @@ -442,10 +442,10 @@ def test_add_solution_output(): phreeqc.add_solution( [ - Solution( + PHRQSol( {"pH": 7.0, "pe": 8.5, "redox": "pe", "temp": 25.0, "units": "mol/kgw", "water": 0.9970480319717386} ), - Solution( + PHRQSol( {"pH": 10.0, "pe": 8.5, "redox": "pe", "temp": 50.0, "units": "mol/kgw", "water": 0.9970480319717386} ), ] @@ -462,7 +462,7 @@ def test_add_solution_output(): def test_kgw(): phreeqc = Phreeqc() phreeqc.add_solution( - Solution({"pH": 7.0, "pe": 8.5, "redox": "pe", "temp": 25.0, "units": "mol/kgw", "water": 0.9970480319717386}) + PHRQSol({"pH": 7.0, "pe": 8.5, "redox": "pe", "temp": 25.0, "units": "mol/kgw", "water": 0.9970480319717386}) ) assert phreeqc[0].get_kgw() == approx(0.9970480319717386) @@ -470,13 +470,13 @@ def test_kgw(): def test_speciate(): phreeqc = Phreeqc() phreeqc.add_solution( - Solution({"pH": 7.0, "pe": 8.5, "redox": "pe", "temp": 25.0, "units": "mol/kgw", "water": 0.9970480319717386}) + PHRQSol({"pH": 7.0, "pe": 8.5, "redox": "pe", "temp": 25.0, "units": "mol/kgw", "water": 0.9970480319717386}) ) assert phreeqc[0].get_species_list() == ["OH-", "H+", "O2", "H2"] def test_speciate_get_molality(): - solution = Solution( + solution = PHRQSol( {"pH": 7.0, "pe": 8.5, "redox": "pe", "temp": 25.0, "units": "mol/kgw", "water": 0.9970480319717386} ) phreeqc = Phreeqc() @@ -487,7 +487,7 @@ def test_speciate_get_molality(): def test_speciate_get_moles(): - solution = Solution( + solution = PHRQSol( {"pH": 7.0, "pe": 8.5, "redox": "pe", "temp": 25.0, "units": "mol/kgw", "water": 0.9970480319717386} ) phreeqc = Phreeqc() @@ -498,7 +498,7 @@ def test_speciate_get_moles(): def test_speciate_species_moles(): - solution = Solution( + solution = PHRQSol( {"pH": 7.0, "pe": 8.5, "redox": "pe", "temp": 25.0, "units": "mol/kgw", "water": 0.9970480319717386} ) phreeqc = Phreeqc() @@ -518,8 +518,8 @@ def test_speciate_species_moles(): def test_species_all_props(): solutions = [ - Solution({"pH": 7.0, "pe": 8.5, "redox": "pe", "temp": 25.0, "units": "mol/kgw", "water": 0.9970480319717386}), - Solution({"pH": 10.0, "pe": 8.5, "redox": "pe", "temp": 50.0, "units": "mol/kgw", "water": 0.9970480319717386}), + PHRQSol({"pH": 7.0, "pe": 8.5, "redox": "pe", "temp": 25.0, "units": "mol/kgw", "water": 0.9970480319717386}), + PHRQSol({"pH": 10.0, "pe": 8.5, "redox": "pe", "temp": 50.0, "units": "mol/kgw", "water": 0.9970480319717386}), ] phreeqc = Phreeqc() phreeqc.add_solution(solutions) @@ -577,7 +577,7 @@ def test_species_all_props(): def test_speciate_get_activity(): - solution = Solution( + solution = PHRQSol( {"pH": 7.0, "pe": 8.5, "redox": "pe", "temp": 25.0, "units": "mol/kgw", "water": 0.9970480319717386} ) phreeqc = Phreeqc() @@ -588,7 +588,7 @@ def test_speciate_get_activity(): def test_get_osmotic_coefficient(): - solution = Solution( + solution = PHRQSol( {"pH": 7.0, "pe": 8.5, "redox": "pe", "temp": 25.0, "units": "mol/kgw", "water": 0.9970480319717386} ) phreeqc = Phreeqc() diff --git a/tests/phreeqc/test_phreeqc_solution.py b/tests/phreeqc/test_phreeqc_solution.py index d1a9f27c..26c42a06 100644 --- a/tests/phreeqc/test_phreeqc_solution.py +++ b/tests/phreeqc/test_phreeqc_solution.py @@ -1,8 +1,8 @@ -from pyEQL_phreeqc import Solution +from pyEQL_phreeqc import PHRQSol def test_create_solution(): - Solution( + PHRQSol( { "Cl": "4.011842831773806", "Na": "4.011842831773806",