diff --git a/src/pyEQL/engines.py b/src/pyEQL/engines.py index 5eecead5..34a4b162 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: @@ -899,21 +899,232 @@ 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 PHRQSol # 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" + + try: + ppsol = self.pp.add_solution(PHRQSol(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: - 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() + return ureg.Quantity(osmotic, "dimensionless") def get_solute_volume(self, solution: "solution.Solution") -> ureg.Quantity: """Return the volume of the solutes.""" # 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/__init__.py b/src/pyEQL/pyEQL-phreeqc/pyEQL_phreeqc/__init__.py index db8b7556..7b8ad1e3 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 PHRQSol -__all__ = ["Phreeqc"] +__all__ = ["PHRQSol", "Phreeqc"] diff --git a/src/pyEQL/pyEQL-phreeqc/pyEQL_phreeqc/core.py b/src/pyEQL/pyEQL-phreeqc/pyEQL_phreeqc/core.py index 4f0c26eb..763d4082 100644 --- a/src/pyEQL/pyEQL-phreeqc/pyEQL_phreeqc/core.py +++ b/src/pyEQL/pyEQL-phreeqc/pyEQL_phreeqc/core.py @@ -1,10 +1,21 @@ +from inspect import cleandoc 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.solution import PHRQSol 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): @@ -14,7 +25,8 @@ 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._solutions: list[Solution] = [] + self._str = "" + 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 @@ -26,30 +38,147 @@ 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): return getattr(self._ext, item) raise AttributeError(f"Phreeqc has no attribute '{item}'") - def add_solution(self, solution_dict: dict) -> None: - solution = Solution(solution_dict) - index = len(self) - self.run_string(f""" - SOLUTION {index} - {solution} - SAVE SOLUTION {index} - END - """) - self._solutions.append(solution) - - def remove_solution(self, index: int) -> Solution: - self.run_string(f""" + def __call__(self, *args, **kwargs): + self.run_string(self._str) + + def __str__(self): + return cleandoc(self._str) + + def clear(self): + self._str = "" + + def accumulate(self, s: str) -> None: + self._str += dedent(s) + + 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: + index = len(self) + solution._phreeqc = ref(self) + solution._number = index + + # TODO: This should go in the PHRQSol 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) -> PHRQSol: + _str = ( + cleandoc(f""" DELETE -solution {index} """) + + "\n" + ) + self.accumulate(_str) + self() return self._solutions.pop(index) + 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]) + + self.clear() + self.accumulate(f""" + SELECTED_OUTPUT + -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), {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 line in output: + solution_i_props = {"species": {}, "eq_species": {}} + n_tokens = len(line) + j = 0 + + 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: + 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) + + 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: def __init__(self, phreeqc: Phreeqc): diff --git a/src/pyEQL/pyEQL-phreeqc/pyEQL_phreeqc/solution.py b/src/pyEQL/pyEQL-phreeqc/pyEQL_phreeqc/solution.py index 683c5141..6a8d02ed 100644 --- a/src/pyEQL/pyEQL-phreeqc/pyEQL_phreeqc/solution.py +++ b/src/pyEQL/pyEQL-phreeqc/pyEQL_phreeqc/solution.py @@ -1,5 +1,69 @@ -class Solution(dict): - # A solution is nothing but a dict of str: Any mapping +from copy import deepcopy +from typing import Any + + +class PHRQSol: + def __init__(self, props, phreeqc=None): + self._phreeqc = phreeqc + self._number = -1 + self._input_props = 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, 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) -> float: + return self._get_calculated_prop("ACT", species=species) + + def get_molality(self, species) -> float: + return self._get_calculated_prop("MOL", species=species) + + 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]: + return {species: self.get_moles(species) for species in self.species} + + def forget(self) -> None: + if self._phreeqc is not None: + self._phreeqc().remove_solution(self._number) 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..92b3b522 100644 --- a/src/pyEQL/solution.py +++ b/src/pyEQL/solution.py @@ -237,6 +237,17 @@ def __init__( self.database.connect() self.logger.debug(f"Connected to property database {self.database!s}") + if engine == "native": + warnings.warn( + '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, + ) + # set the equation of state engine self._engine = engine # self.engine: Optional[EOS] = None diff --git a/tests/phreeqc/test_phreeqc.py b/tests/phreeqc/test_phreeqc.py index 1938353e..df2da2f2 100644 --- a/tests/phreeqc/test_phreeqc.py +++ b/tests/phreeqc/test_phreeqc.py @@ -1,8 +1,9 @@ -import textwrap +from inspect import cleandoc from pathlib import Path import numpy as np -from pyEQL_phreeqc import Phreeqc +from pyEQL_phreeqc import Phreeqc, PHRQSol +from pytest import approx def test_load_database_internal(): @@ -132,17 +133,19 @@ def test_run_simple(): def test_run_add_solution(): phreeqc = Phreeqc() - phreeqc.add_solution( + solution = PHRQSol( {"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 = PHRQSol( {"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 @@ -150,8 +153,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 +163,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 +215,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 +231,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 +243,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 +293,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 +313,286 @@ 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_add_solution_input(): + phreeqc = Phreeqc() + + phreeqc.add_solution( + [ + 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} + ), + ] + ) + + expected = cleandoc(""" + SELECTED_OUTPUT + -reset false + + USER_PUNCH + 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 + 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""") + + assert str(phreeqc) == expected + + +def test_add_solution_output(): + phreeqc = Phreeqc() + + phreeqc.add_solution( + [ + 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} + ), + ] + ) + + # 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 + # +"\n" + [, ] repeated for each (3) equilibrium species + assert phreeqc.get_selected_output_column_count() == 23 + + +def test_kgw(): + phreeqc = Phreeqc() + phreeqc.add_solution( + 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) + + +def test_speciate(): + phreeqc = Phreeqc() + phreeqc.add_solution( + 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 = PHRQSol( + {"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 = PHRQSol( + {"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 = PHRQSol( + {"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 = [ + 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) + + expected = { + 0: { + "CELL_NO": 0, + "TOT['water']": 0.9970480319717386, + "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}, + }, + "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}, + "H2": {"ACT": 5.626700758118202e-41, "MOL": 0.0}, + "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}, + }, + }, + } + + 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 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]) + + 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_get_activity(): + solution = PHRQSol( + {"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_get_osmotic_coefficient(): + solution = PHRQSol( + {"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 89030606..26c42a06 100644 --- a/tests/phreeqc/test_phreeqc_solution.py +++ b/tests/phreeqc/test_phreeqc_solution.py @@ -1,9 +1,8 @@ -from pyEQL_phreeqc import Phreeqc +from pyEQL_phreeqc import PHRQSol -def test_add_solution(): - phreeqc = Phreeqc() - phreeqc.add_solution( +def test_create_solution(): + PHRQSol( { "Cl": "4.011842831773806", "Na": "4.011842831773806", diff --git a/tests/test_pyeql_phreeqc.py b/tests/test_pyeql_phreeqc.py new file mode 100644 index 00000000..9c4bcb81 --- /dev/null +++ b/tests/test_pyeql_phreeqc.py @@ -0,0 +1,416 @@ +""" +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 == 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 + 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) + # 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) + 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