From 0d304eb97243657602acc334a301f55fec6aa8c1 Mon Sep 17 00:00:00 2001 From: Ryan Kingsbury Date: Sat, 14 Feb 2026 11:48:41 -0500 Subject: [PATCH 1/3] make the warning into an error --- src/pyEQL/engines.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/pyEQL/engines.py b/src/pyEQL/engines.py index d62b1ebb..96b04828 100644 --- a/src/pyEQL/engines.py +++ b/src/pyEQL/engines.py @@ -784,7 +784,7 @@ def equilibrate( } ) elif abs(orig_el_amount - new_el_amount) / orig_el_amount > _rtol: - logger.warning( + logger.error( f"PHREEQC returned a total Element {el} concentration of {new_el_amount} mol, " f"which differs from the original concentration of {orig_el_amount}. This " "should never occur and indicates an error in the PHREEQC database or calculation." From 7164d71dcfcdad5b8d74583ac31c395732a56d58 Mon Sep 17 00:00:00 2001 From: Ryan Kingsbury Date: Sat, 14 Feb 2026 12:35:28 -0500 Subject: [PATCH 2/3] NativeEOS inherits from PHREEQC --- src/pyEQL/engines.py | 1392 ++++++++++++++++++++---------------------- 1 file changed, 651 insertions(+), 741 deletions(-) diff --git a/src/pyEQL/engines.py b/src/pyEQL/engines.py index 96b04828..ec9ccf49 100644 --- a/src/pyEQL/engines.py +++ b/src/pyEQL/engines.py @@ -13,7 +13,7 @@ import warnings from abc import ABC, abstractmethod from pathlib import Path -from typing import TYPE_CHECKING, Literal +from typing import TYPE_CHECKING, Literal, Self from phreeqpython import PhreeqPython @@ -154,18 +154,15 @@ def equilibrate(self, solution: "solution.Solution") -> None: return -class NativeEOS(EOS): - """ - pyEQL's native EOS. Uses the Pitzer model when possible, falls - back to other models (e.g. Debye-Huckel) based on ionic strength - if sufficient parameters are not available. - """ +class Phreeqc2026EOS(EOS): + """Engine based on the PhreeqC model, as implemented in the pyphreeqc + module of pyEQL.""" 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: @@ -179,21 +176,19 @@ def __init__( may offer improved prediction of LSI but currently these databases are not usable because they do not allow for conductivity calculations. """ + + from pyEQL.phreeqc import IS_AVAILABLE, Phreeqc # noqa: PLC0415 + + if not IS_AVAILABLE: + raise RuntimeError("pyEQL phreeqc support is not available in this installation") + self.phreeqc_db = phreeqc_db # database files in this list are not distributed with phreeqpython self.db_path = ( Path(os.path.dirname(__file__)) / "database" if self.phreeqc_db in ["llnl.dat", "geothermal.dat"] else None ) # create the PhreeqcPython instance - # try/except added to catch unsupported architectures, such as Apple Silicon - try: - self.pp = PhreeqPython(database=self.phreeqc_db, database_directory=self.db_path) - except OSError: - logger.error( - "OSError encountered when trying to instantiate phreeqpython. Most likely this means you" - " are running on an architecture that is not supported by PHREEQC, such as Apple M1/M2 chips." - " pyEQL will work, but equilibrate() will have no effect." - ) + self.pp = Phreeqc(database=self.phreeqc_db, database_directory=self.db_path) # attributes to hold the PhreeqPython solution. self.ppsol = None # store the solution composition to see whether we need to re-instantiate the solution @@ -201,6 +196,9 @@ 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 PHRQSol # noqa: PLC0415 + self._stored_comp = solution.components.copy() solv_mass = solution.solvent_mass.to("kg").magnitude # inherit bulk solution properties @@ -254,354 +252,483 @@ 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" - # create the PHREEQC solution object try: - ppsol = self.pp.add_solution(d) + ppsol = self.pp.add_solution(PHRQSol(d)) except Exception as e: - print(d) # catch problems with the input to phreeqc raise ValueError( "There is a problem with your input. The error message received from " - f" phreeqpython is:\n\n {e}\n Check your input arguments, especially " + 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: - """Remove the PhreeqPython solution from memory.""" if self.ppsol is not None: - self.ppsol.forget() + 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): - r""" - Whenever the appropriate parameters are available, the Pitzer model [may11]_ is used. - If no Pitzer parameters are available, then the appropriate equations are selected - according to the following logic: [stumm96]_. + 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. - I <= 0.0005: Debye-Huckel equation - 0.005 < I <= 0.1: Guntelberg approximation - 0.1 < I <= 0.5: Davies equation - I > 0.5: Raises a warning and returns activity coefficient = 1 + Args: + solution: a pyEQL Solution object to equilibrate. The composition of this solution will be modified in-place. + 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) - The ionic strength, activity coefficients, and activities are all - calculated based on the molal (mol/kg) concentration scale. If a different - scale is given as input, then the molal-scale activity coefficient :math:`\gamma_\pm` is - converted according to [rbs68]_ + # store the original solvent mass + orig_solvent_moles = solution.components[solution.solvent] - .. math:: f_\pm = \gamma_\pm * (1 + M_w \sum_i \nu_i m_i) + # 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) - .. math:: y_\pm = \frac{m \rho_w}{C \gamma_\pm} + # Use supplied gases, merged with atmospheric gases + # if atmosphere == True + gases = (ATMOSPHERE if atmosphere else {}) | (gases or {}) - where :math:`f_\pm` is the rational activity coefficient, :math:`M_w` is - the molecular weight of water, the summation represents the total molality of - all solute species, :math:`y_\pm` is the molar activity coefficient, - :math:`\rho_w` is the density of pure water, :math:`m` and :math:`C` are - the molal and molar concentrations of the chosen salt (not individual solute), respectively. + # 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)) - Args: - solute: String representing the name of the solute of interest + 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)} - Returns: - The mean ion activity coefficient of the solute in question. + if phases: + phase_names = list(phases.keys()) + saturation_indices = [v[0] for v in phases.values()] + amounts = [v[1] for v in phases.values()] - Notes: - For multicomponent mixtures, pyEQL implements the "effective Pitzer model" - presented by Mistry et al. [mistry13]_. In this model, the activity coefficient - of a salt in a multicomponent mixture is calculated using an "effective - molality," which is the molality that would result in a single-salt - mixture with the same total ionic strength as the multicomponent solution. + try: + self.ppsol.equalize(phases=phase_names, to_si=saturation_indices, in_phase=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 - .. math:: m_{effective} = \frac{2 I}{(\nu_{+} z_{+}^2 + \nu_{-}- z_{-}^2)} + 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 - References: - .. [may11] May, P. M., Rowland, D., Hefter, G., & Königsberger, E. (2011). - A Generic and Updatable Pitzer Characterization of Aqueous Binary Electrolyte Solutions at 1 bar - and 25 °C. *Journal of Chemical & Engineering Data*, 56(12), 5066-5077. doi:10.1021/je2009329 + # 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." + ) - .. [stumm96] Stumm, Werner and Morgan, James J. *Aquatic Chemistry*, 3rd ed, - pp 165. Wiley Interscience, 1996. + # tolerance (in moles) for detecting cases where an element amount + # is no longer balanced because of species that are not recognized + # by PHREEQC. + _rtol = 0.05 # differing by more than 5% - .. [rbs68] Robinson, R. A.; Stokes, R. H. Electrolyte Solutions: Second Revised - Edition; Butterworths: London, 1968, p.32. + 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, [])]) - See Also: - :attr:`pyEQL.solution.Solution.ionic_strength` - :func:`pyEQL.activity_correction.get_activity_coefficient_debyehuckel` - :func:`pyEQL.activity_correction.get_activity_coefficient_guntelberg` - :func:`pyEQL.activity_correction.get_activity_coefficient_davies` - :func:`pyEQL.activity_correction.get_activity_coefficient_pitzer` - """ - # identify the predominant salt that this ion is a member of - salt = None - rform = standardize_formula(solute) - for d in solution.get_salt_dict().values(): - if rform == d["salt"].cation or rform == d["salt"].anion: - salt = d["salt"] - break + # If this element went "missing", add back all components that + # contain this element (for any valence value) + if new_el_amount == 0 and orig_el_amount > 0: + 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 + } + ) + elif abs(orig_el_amount - new_el_amount) / orig_el_amount > _rtol: + logger.error( + f"PHREEQC returned a total Element {el} concentration of {new_el_amount} mol, " + f"which differs from the original concentration of {orig_el_amount}. This " + "should never occur and indicates an error in the PHREEQC database or calculation." + ) - # show an error if no salt can be found that contains the solute - if salt is None: - logger.error(f"No salts found that contain solute {solute}. Returning unit activity coefficient.") - return ureg.Quantity(1, "dimensionless") + # 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() - # use the Pitzer model for higher ionic strength, if the parameters are available - # search for Pitzer parameters - param = solution.get_property(salt.formula, "model_parameters.activity_pitzer") - if param is not None: - # TODO - consider re-enabling a log message recording what salt(s) are used as basis for activity calculation - logger.info(f"Calculating activity coefficient based on parent salt {salt.formula}") + # 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 - # determine alpha1 and alpha2 based on the type of salt - # see the May reference for the rules used to determine - # alpha1 and alpha2 based on charge - if salt.z_cation >= 2 and salt.z_anion <= -2: - if salt.z_cation >= 3 or salt.z_anion <= -3: - alpha1 = 2.0 - alpha2 = 50.0 - else: - alpha1 = 1.4 - alpha2 = 12 - else: - alpha1 = 2.0 - alpha2 = 0.0 + def get_activity_coefficient(self, solution: "solution.Solution", solute: str) -> ureg.Quantity: + """ + 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) - # determine the average molality of the salt - # this is necessary for solutions inside e.g. an ion exchange - # membrane, where the cation and anion concentrations may be - # unequal - # molality = (solution.get_amount(salt.cation,'mol/kg')/salt.nu_cation+solution.get_amount(salt.anion,'mol/kg')/salt.nu_anion)/2 + # 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 - # determine the effective molality of the salt in the solution - molality = salt.get_effective_molality(solution.ionic_strength) + # calculate the molal scale activity coefficient + act = self.ppsol.get_activity(k) / self.ppsol.get_molality(k) - activity_coefficient = ac.get_activity_coefficient_pitzer( - solution.ionic_strength, - molality, - alpha1, - alpha2, - ureg.Quantity(param["Beta0"]["value"]).magnitude, - ureg.Quantity(param["Beta1"]["value"]).magnitude, - ureg.Quantity(param["Beta2"]["value"]).magnitude, - ureg.Quantity(param["Cphi"]["value"]).magnitude, - salt.z_cation, - salt.z_anion, - salt.nu_cation, - salt.nu_anion, - str(solution.temperature), - ) + return ureg.Quantity(act, "dimensionless") - logger.debug( - f"Calculated activity coefficient of species {solute} as {activity_coefficient} based on salt" - f" {salt} using Pitzer model" - ) - molal = activity_coefficient + def get_osmotic_coefficient(self, solution: "solution.Solution") -> ureg.Quantity: + if self.ppsol is not None: + self.ppsol.forget() + self._setup_ppsol(solution) - # for very low ionic strength, use the Debye-Huckel limiting law - elif solution.ionic_strength.magnitude <= 0.005: - logger.debug( - f"Ionic strength = {solution.ionic_strength}. Using Debye-Huckel to calculate activity coefficient." - ) - molal = ac.get_activity_coefficient_debyehuckel( - solution.ionic_strength, - solution.get_property(solute, "charge"), - str(solution.temperature), - ) + osmotic = self.ppsol.get_osmotic_coefficient() + if osmotic == 0: + # PHREEQC returns 0 when it assumes a unit osmotic coefficient, so we need to catch that case and return 1 instead. + osmotic = 1 + return ureg.Quantity(osmotic, "dimensionless") - # use the Guntelberg approximation for 0.005 < I < 0.1 - elif solution.ionic_strength.magnitude <= 0.1: - logger.debug( - f"Ionic strength = {solution.ionic_strength}. Using Guntelberg to calculate activity coefficient." - ) - molal = ac.get_activity_coefficient_guntelberg( - solution.ionic_strength, - solution.get_property(solute, "charge"), - str(solution.temperature), - ) + def get_solute_volume(self, solution: "solution.Solution") -> ureg.Quantity: + """Return the volume of the solutes.""" + # TODO - see if we can access molar volume or solute volume via the pyEQL-phreeqc wrapper + return ureg.Quantity(0, "L") - # use the Davies equation for 0.1 < I < 0.5 - elif solution.ionic_strength.magnitude <= 0.5: - logger.debug( - f"Ionic strength = {solution.ionic_strength}. Using Davies equation to calculate activity coefficient." - ) - molal = ac.get_activity_coefficient_davies( - solution.ionic_strength, - solution.get_property(solute, "charge"), - str(solution.temperature), - ) + def __deepcopy__(self, memo) -> Self: + # custom deepcopy required because the PhreeqPython instance used by the Native and Phreeqc engines + # is not pickle-able. - else: + cls = self.__class__ + result = cls.__new__(cls) + memo[id(self)] = result + for k, v in self.__dict__.items(): + if k == "pp": + result.pp = PhreeqPython(database=self.phreeqc_db, database_directory=self.db_path) + continue + setattr(result, k, copy.deepcopy(v, memo)) + return result + + +class PhreeqcEOS(Phreeqc2026EOS): + """Engine based on the PhreeqC model, as implemented via the phreeqpython package.""" + + def __init__( + self, + phreeqc_db: Literal[ + "phreeqc.dat", "vitens.dat", "wateq4f_PWN.dat", "pitzer.dat", "llnl.dat", "geothermal.dat" + ] = "phreeqc.dat", + ) -> None: + """ + Args: + phreeqc_db: Name of the PHREEQC database file to use for solution thermodynamics + and speciation calculations. Generally speaking, `llnl.dat` is recommended + for moderate salinity water and prediction of mineral solubilities, + `wateq4f_PWN.dat` is recommended for low to moderate salinity waters. It is + similar to vitens.dat but has many more species. `pitzer.dat` is recommended + when accurate activity coefficients in solutions above 1 M TDS are desired, but + it has fewer species than the other databases. `llnl.dat` and `geothermal.dat` + may offer improved prediction of LSI but currently these databases are not + usable because they do not allow for conductivity calculations. + """ + self.phreeqc_db = phreeqc_db + # database files in this list are not distributed with phreeqpython + self.db_path = ( + Path(os.path.dirname(__file__)) / "database" if self.phreeqc_db in ["llnl.dat", "geothermal.dat"] else None + ) + # create the PhreeqcPython instance + # try/except added to catch unsupported architectures, such as Apple Silicon + try: + self.pp = PhreeqPython(database=self.phreeqc_db, database_directory=self.db_path) + except OSError: logger.error( - f"Ionic strength too high to estimate activity for species {solute}. Specify parameters for Pitzer " - "model. Returning unit activity coefficient" + "OSError encountered when trying to instantiate phreeqpython. Most likely this means you" + " are running on an architecture that is not supported by PHREEQC, such as Apple M1/M2 chips." + " pyEQL will work, but equilibrate() will have no effect." ) + # attributes to hold the PhreeqPython solution. + self.ppsol = None + # store the solution composition to see whether we need to re-instantiate the solution + self._stored_comp = None - molal = ureg.Quantity(1, "dimensionless") + def _setup_ppsol(self, solution: "solution.Solution") -> None: + """Helper method to set up a PhreeqPython solution for subsequent analysis.""" + 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" - return molal + # 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. - def get_osmotic_coefficient(self, solution: "solution.Solution") -> ureg.Quantity: - r""" - Return the *molal scale* osmotic coefficient of solute, given a Solution - object. + # 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 - Osmotic coefficient is calculated using the Pitzer model. [may11]_ If appropriate parameters for - the model are not available, then pyEQL raises a WARNING and returns an osmotic - coefficient of 1. + 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) - If the 'rational' scale is given as input, then the molal-scale osmotic - coefficient :math:`\phi` is converted according to [rbs68]_ + # 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" - .. math:: g = - \phi M_{w} \frac{\sum_{i} \nu_{i} m_{i}}{\ln x_{w}} + # create the PHREEQC solution object + try: + ppsol = self.pp.add_solution(d) + except Exception as e: + print(d) + # catch problems with the input to phreeqc + raise ValueError( + "There is a problem with your input. The error message received from " + f" phreeqpython is:\n\n {e}\n Check your input arguments, especially " + "the composition dictionary, and try again." + ) - where :math:`g` is the rational osmotic coefficient, :math:`M_{w}` is - the molecular weight of water, the summation represents the total molality of - all solute species, and :math:`x_{w}` is the mole fraction of water. + self.ppsol = ppsol - Returns: - Quantity: - The osmotic coefficient + def _destroy_ppsol(self) -> None: + """Remove the PhreeqPython solution from memory.""" + if self.ppsol is not None: + self.ppsol.forget() + self.ppsol = None - See Also: - :meth:`pyEQL.solution.Solution.get_water_activity` - :meth:`pyEQL.solution.Solution.get_salt` - :attr:`pyEQL.solution.Solution.ionic_strength` + def get_activity_coefficient(self, solution: "solution.Solution", solute: str) -> ureg.Quantity: + """ + 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) - Notes: - For multicomponent mixtures, pyEQL adopts the "effective Pitzer model" - presented by Mistry et al. [mistry13]_. In this approach, the osmotic coefficient of - each individual salt is calculated using the normal Pitzer model based - on its respective concentration. Then, an effective osmotic coefficient - is calculated as the concentration-weighted average of the individual - osmotic coefficients. + # 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 - For example, in a mixture of 0.5 M NaCl and 0.5 M KBr, one would calculate - the osmotic coefficient for each salt using a concentration of 0.5 M and - an ionic strength of 1 M. Then, one would average the two resulting - osmotic coefficients to obtain an effective osmotic coefficient for the - mixture. + # calculate the molal scale activity coefficient + # act = self.ppsol.activity(k, "mol") / self.ppsol.molality(k, "mol") + act = self.ppsol.pp.ip.get_activity(self.ppsol.number, k) / self.ppsol.pp.ip.get_molality(self.ppsol.number, k) - (Note: in the paper referenced below, the effective osmotic coefficient is determined by weighting - using the "effective molality" rather than the true molality. Subsequent checking and correspondence with - the author confirmed that the weight factor should be the true molality, and that is what is implemented - in pyEQL.) + return ureg.Quantity(act, "dimensionless") - Examples: - >>> s1 = pyEQL.Solution({'Na+': '0.2 mol/kg', 'Cl-': '0.2 mol/kg'}) - >>> s1.get_osmotic_coefficient() - + def get_osmotic_coefficient(self, solution: "solution.Solution") -> ureg.Quantity: + """ + Return the *molal scale* osmotic coefficient of solute, given a Solution + object. - >>> s1 = pyEQL.Solution({'Mg+2': '0.3 mol/kg', 'Cl-': '0.6 mol/kg'},temperature='30 degC') - >>> s1.get_osmotic_coefficient() - + PHREEQC appears to assume a unit osmotic coefficient unless the pitzer database + is used. Unfortunately, there is no easy way to access the osmotic coefficient + via phreeqcpython + """ + # TODO - find a way to access or calculate osmotic coefficient + return ureg.Quantity(1, "dimensionless") - References: - [may11] - May, P. M., Rowland, D., Hefter, G., & Königsberger, E. (2011). - A Generic and Updatable Pitzer Characterization of Aqueous Binary Electrolyte Solutions at 1 bar - and 25 °C. Journal of Chemical & Engineering Data, 56(12), 5066-5077. doi:10.1021/je2009329 +class NativeEOS(PhreeqcEOS): + """ + pyEQL's native EOS. Uses the Pitzer model when possible, falls + back to other models (e.g. Debye-Huckel) based on ionic strength + if sufficient parameters are not available. + """ - [rbs68] + def __init__( + self, + phreeqc_db: Literal[ + "phreeqc.dat", "vitens.dat", "wateq4f_PWN.dat", "pitzer.dat", "llnl.dat", "geothermal.dat" + ] = "llnl.dat", + ) -> None: + """ + Args: + phreeqc_db: Name of the PHREEQC database file to use for solution thermodynamics + and speciation calculations. Generally speaking, `llnl.dat` is recommended + for moderate salinity water and prediction of mineral solubilities, + `wateq4f_PWN.dat` is recommended for low to moderate salinity waters. It is + similar to vitens.dat but has many more species. `pitzer.dat` is recommended + when accurate activity coefficients in solutions above 1 M TDS are desired, but + it has fewer species than the other databases. `llnl.dat` and `geothermal.dat` + may offer improved prediction of LSI but currently these databases are not + usable because they do not allow for conductivity calculations. + """ + self.phreeqc_db = phreeqc_db + # database files in this list are not distributed with phreeqpython + self.db_path = ( + Path(os.path.dirname(__file__)) / "database" if self.phreeqc_db in ["llnl.dat", "geothermal.dat"] else None + ) + # create the PhreeqcPython instance + # try/except added to catch unsupported architectures, such as Apple Silicon + try: + self.pp = PhreeqPython(database=self.phreeqc_db, database_directory=self.db_path) + except OSError: + logger.error( + "OSError encountered when trying to instantiate phreeqpython. Most likely this means you" + " are running on an architecture that is not supported by PHREEQC, such as Apple M1/M2 chips." + " pyEQL will work, but equilibrate() will have no effect." + ) + # attributes to hold the PhreeqPython solution. + self.ppsol = None + # store the solution composition to see whether we need to re-instantiate the solution + self._stored_comp = None - Robinson, R. A.; Stokes, R. H. Electrolyte Solutions: Second Revised Edition; Butterworths: London, 1968, - p.32. + def get_activity_coefficient(self, solution: "solution.Solution", solute: str): + r""" + Whenever the appropriate parameters are available, the Pitzer model [may11]_ is used. + If no Pitzer parameters are available, then the appropriate equations are selected + according to the following logic: [stumm96]_. - [mistry13] + I <= 0.0005: Debye-Huckel equation + 0.005 < I <= 0.1: Guntelberg approximation + 0.1 < I <= 0.5: Davies equation + I > 0.5: Raises a warning and returns activity coefficient = 1 - Mistry, K. H.; Hunter, H. a.; Lienhard V, J. H. Effect of composition and nonideal solution - behavior on desalination calculations for mixed electrolyte solutions with comparison to - seawater. Desalination 2013, 318, 34-47. + The ionic strength, activity coefficients, and activities are all + calculated based on the molal (mol/kg) concentration scale. If a different + scale is given as input, then the molal-scale activity coefficient :math:`\gamma_\pm` is + converted according to [rbs68]_ - """ - ionic_strength = solution.ionic_strength + .. math:: f_\pm = \gamma_\pm * (1 + M_w \sum_i \nu_i m_i) - effective_osmotic_sum = 0 - molality_sum = 0 + .. math:: y_\pm = \frac{m \rho_w}{C \gamma_\pm} - # loop through all the salts in the solution, calculate the osmotic - # coefficint for each, and average them into an effective osmotic - # coefficient - for d in solution.get_salt_dict().values(): - item = d["salt"] - # determine alpha1 and alpha2 based on the type of salt - # see the May reference for the rules used to determine - # alpha1 and alpha2 based on charge - if item.z_cation >= 2 and item.z_anion <= -2: - if item.z_cation >= 3 or item.z_anion <= -3: - alpha1 = 2.0 - alpha2 = 50.0 - else: - alpha1 = 1.4 - alpha2 = 12.0 - else: - alpha1 = 2.0 - alpha2 = 0 + where :math:`f_\pm` is the rational activity coefficient, :math:`M_w` is + the molecular weight of water, the summation represents the total molality of + all solute species, :math:`y_\pm` is the molar activity coefficient, + :math:`\rho_w` is the density of pure water, :math:`m` and :math:`C` are + the molal and molar concentrations of the chosen salt (not individual solute), respectively. - # set the concentration as the average concentration of the cation and - # anion in the salt, accounting for stoichiometry - # concentration = (solution.get_amount(Salt.cation,'mol/kg')/Salt.nu_cation + \ - # solution.get_amount(Salt.anion,'mol/kg')/Salt.nu_anion)/2 + Args: + solute: String representing the name of the solute of interest - # get the effective molality of the salt - concentration = ureg.Quantity(d["mol"], "mol") / solution.solvent_mass + Returns: + The mean ion activity coefficient of the solute in question. - molality_sum += concentration + Notes: + For multicomponent mixtures, pyEQL implements the "effective Pitzer model" + presented by Mistry et al. [mistry13]_. In this model, the activity coefficient + of a salt in a multicomponent mixture is calculated using an "effective + molality," which is the molality that would result in a single-salt + mixture with the same total ionic strength as the multicomponent solution. - param = solution.get_property(item.formula, "model_parameters.activity_pitzer") - if param is not None: - osmotic_coefficient = ac.get_osmotic_coefficient_pitzer( - ionic_strength, - concentration, - alpha1, - alpha2, - ureg.Quantity(param["Beta0"]["value"]).magnitude, - ureg.Quantity(param["Beta1"]["value"]).magnitude, - ureg.Quantity(param["Beta2"]["value"]).magnitude, - ureg.Quantity(param["Cphi"]["value"]).magnitude, - item.z_cation, - item.z_anion, - item.nu_cation, - item.nu_anion, - str(solution.temperature), - ) + .. math:: m_{effective} = \frac{2 I}{(\nu_{+} z_{+}^2 + \nu_{-}- z_{-}^2)} - logger.debug( - f"Calculated osmotic coefficient of water as {osmotic_coefficient} based on salt " - f"{item.formula} using Pitzer model" - ) - effective_osmotic_sum += concentration * osmotic_coefficient + References: + .. [may11] May, P. M., Rowland, D., Hefter, G., & Königsberger, E. (2011). + A Generic and Updatable Pitzer Characterization of Aqueous Binary Electrolyte Solutions at 1 bar + and 25 °C. *Journal of Chemical & Engineering Data*, 56(12), 5066-5077. doi:10.1021/je2009329 - else: - logger.debug( - f"Returning unit osmotic coefficient for salt {item.formula} because Pitzer parameters are not" - "available in database." - ) - effective_osmotic_sum += concentration * 1 + .. [stumm96] Stumm, Werner and Morgan, James J. *Aquatic Chemistry*, 3rd ed, + pp 165. Wiley Interscience, 1996. - try: - return effective_osmotic_sum / molality_sum - except ZeroDivisionError: - # this means the solution is empty - return ureg.Quantity(1.0) + .. [rbs68] Robinson, R. A.; Stokes, R. H. Electrolyte Solutions: Second Revised + Edition; Butterworths: London, 1968, p.32. - def get_solute_volume(self, solution: "solution.Solution") -> ureg.Quantity: - """Return the volume of the solutes.""" - # identify the predominant salt in the solution - salt = solution.get_salt() - solute_vol = ureg.Quantity(0, "L") + See Also: + :attr:`pyEQL.solution.Solution.ionic_strength` + :func:`pyEQL.activity_correction.get_activity_coefficient_debyehuckel` + :func:`pyEQL.activity_correction.get_activity_coefficient_guntelberg` + :func:`pyEQL.activity_correction.get_activity_coefficient_davies` + :func:`pyEQL.activity_correction.get_activity_coefficient_pitzer` + """ + # identify the predominant salt that this ion is a member of + salt = None + rform = standardize_formula(solute) + for d in solution.get_salt_dict().values(): + if rform == d["salt"].cation or rform == d["salt"].anion: + salt = d["salt"] + break - # use the pitzer approach if parameters are available - pitzer_calc = False - param = None if salt is None else solution.get_property(salt.formula, "model_parameters.molar_volume_pitzer") + # show an error if no salt can be found that contains the solute + if salt is None: + logger.error(f"No salts found that contain solute {solute}. Returning unit activity coefficient.") + return ureg.Quantity(1, "dimensionless") + # use the Pitzer model for higher ionic strength, if the parameters are available + # search for Pitzer parameters + param = solution.get_property(salt.formula, "model_parameters.activity_pitzer") if param is not None: - # determine the average molality of the salt - # this is necessary for solutions inside e.g. an ion exchange - # membrane, where the cation and anion concentrations may be - # unequal - molality = (solution.get_amount(salt.cation, "mol/kg") + solution.get_amount(salt.anion, "mol/kg")) / 2 + # TODO - consider re-enabling a log message recording what salt(s) are used as basis for activity calculation + logger.info(f"Calculating activity coefficient based on parent salt {salt.formula}") # determine alpha1 and alpha2 based on the type of salt # see the May reference for the rules used to determine @@ -617,7 +744,16 @@ def get_solute_volume(self, solution: "solution.Solution") -> ureg.Quantity: alpha1 = 2.0 alpha2 = 0.0 - apparent_vol = ac.get_apparent_volume_pitzer( + # determine the average molality of the salt + # this is necessary for solutions inside e.g. an ion exchange + # membrane, where the cation and anion concentrations may be + # unequal + # molality = (solution.get_amount(salt.cation,'mol/kg')/salt.nu_cation+solution.get_amount(salt.anion,'mol/kg')/salt.nu_anion)/2 + + # determine the effective molality of the salt in the solution + molality = salt.get_effective_molality(solution.ionic_strength) + + activity_coefficient = ac.get_activity_coefficient_pitzer( solution.ionic_strength, molality, alpha1, @@ -626,7 +762,6 @@ def get_solute_volume(self, solution: "solution.Solution") -> ureg.Quantity: ureg.Quantity(param["Beta1"]["value"]).magnitude, ureg.Quantity(param["Beta2"]["value"]).magnitude, ureg.Quantity(param["Cphi"]["value"]).magnitude, - ureg.Quantity(param["V_o"]["value"]).magnitude, salt.z_cation, salt.z_anion, salt.nu_cation, @@ -634,505 +769,280 @@ def get_solute_volume(self, solution: "solution.Solution") -> ureg.Quantity: str(solution.temperature), ) - solute_vol += ( - apparent_vol - * ( - solution.get_amount(salt.cation, "mol") / salt.nu_cation - + solution.get_amount(salt.anion, "mol") / salt.nu_anion - ) - / 2 + logger.debug( + f"Calculated activity coefficient of species {solute} as {activity_coefficient} based on salt" + f" {salt} using Pitzer model" ) + molal = activity_coefficient - pitzer_calc = True - - logger.debug(f"Updated solution volume using Pitzer model for solute {salt.formula}") + # for very low ionic strength, use the Debye-Huckel limiting law + elif solution.ionic_strength.magnitude <= 0.005: + logger.debug( + f"Ionic strength = {solution.ionic_strength}. Using Debye-Huckel to calculate activity coefficient." + ) + molal = ac.get_activity_coefficient_debyehuckel( + solution.ionic_strength, + solution.get_property(solute, "charge"), + str(solution.temperature), + ) - # add the partial molar volume of any other solutes, except for water - # or the parent salt, which is already accounted for by the Pitzer parameters - for solute, mol in solution.components.items(): - # ignore water - if solute in ["H2O", "HOH", "H2O(aq)"]: - continue + # use the Guntelberg approximation for 0.005 < I < 0.1 + elif solution.ionic_strength.magnitude <= 0.1: + logger.debug( + f"Ionic strength = {solution.ionic_strength}. Using Guntelberg to calculate activity coefficient." + ) + molal = ac.get_activity_coefficient_guntelberg( + solution.ionic_strength, + solution.get_property(solute, "charge"), + str(solution.temperature), + ) - # ignore the salt cation and anion, if already accounted for by Pitzer - if pitzer_calc is True and solute in [salt.anion, salt.cation]: - continue + # use the Davies equation for 0.1 < I < 0.5 + elif solution.ionic_strength.magnitude <= 0.5: + logger.debug( + f"Ionic strength = {solution.ionic_strength}. Using Davies equation to calculate activity coefficient." + ) + molal = ac.get_activity_coefficient_davies( + solution.ionic_strength, + solution.get_property(solute, "charge"), + str(solution.temperature), + ) - part_vol = solution.get_property(solute, "size.molar_volume") - if part_vol is not None: - solute_vol += part_vol * ureg.Quantity(mol, "mol") - logger.debug(f"Updated solution volume using direct partial molar volume for solute {solute}") + else: + logger.error( + f"Ionic strength too high to estimate activity for species {solute}. Specify parameters for Pitzer " + "model. Returning unit activity coefficient" + ) - else: - logger.warning( - f"Volume of solute {solute} will be ignored because partial molar volume data are not available." - ) + molal = ureg.Quantity(1, "dimensionless") - return solute_vol.to("L") + return molal - 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. + def get_osmotic_coefficient(self, solution: "solution.Solution") -> ureg.Quantity: + r""" + Return the *molal scale* osmotic coefficient of solute, given a Solution + object. - 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) + Osmotic coefficient is calculated using the Pitzer model. [may11]_ If appropriate parameters for + the model are not available, then pyEQL raises a WARNING and returns an osmotic + coefficient of 1. - # store the original solvent mass - orig_solvent_moles = solution.components[solution.solvent] + If the 'rational' scale is given as input, then the molal-scale osmotic + coefficient :math:`\phi` is converted according to [rbs68]_ - # 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) + .. math:: g = - \phi M_{w} \frac{\sum_{i} \nu_{i} m_{i}}{\ln x_{w}} - # Use supplied gases, merged with atmospheric gases - # if atmosphere == True - gases = (ATMOSPHERE if atmosphere else {}) | (gases or {}) + where :math:`g` is the rational osmotic coefficient, :math:`M_{w}` is + the molecular weight of water, the summation represents the total molality of + all solute species, and :math:`x_{w}` is the mole fraction of water. - # 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)) + Returns: + Quantity: + The osmotic coefficient - 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)} + See Also: + :meth:`pyEQL.solution.Solution.get_water_activity` + :meth:`pyEQL.solution.Solution.get_salt` + :attr:`pyEQL.solution.Solution.ionic_strength` - if phases: - phase_names = list(phases.keys()) - saturation_indices = [v[0] for v in phases.values()] - amounts = [v[1] for v in phases.values()] + Notes: + For multicomponent mixtures, pyEQL adopts the "effective Pitzer model" + presented by Mistry et al. [mistry13]_. In this approach, the osmotic coefficient of + each individual salt is calculated using the normal Pitzer model based + on its respective concentration. Then, an effective osmotic coefficient + is calculated as the concentration-weighted average of the individual + osmotic coefficients. - try: - self.ppsol.equalize(phases=phase_names, to_si=saturation_indices, in_phase=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 + For example, in a mixture of 0.5 M NaCl and 0.5 M KBr, one would calculate + the osmotic coefficient for each salt using a concentration of 0.5 M and + an ionic strength of 1 M. Then, one would average the two resulting + osmotic coefficients to obtain an effective osmotic coefficient for the + mixture. - 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 + (Note: in the paper referenced below, the effective osmotic coefficient is determined by weighting + using the "effective molality" rather than the true molality. Subsequent checking and correspondence with + the author confirmed that the weight factor should be the true molality, and that is what is implemented + in pyEQL.) - # 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." - ) + Examples: + >>> s1 = pyEQL.Solution({'Na+': '0.2 mol/kg', 'Cl-': '0.2 mol/kg'}) + >>> s1.get_osmotic_coefficient() + - # tolerance (in moles) for detecting cases where an element amount - # is no longer balanced because of species that are not recognized - # by PHREEQC. - _rtol = 0.05 # differing by more than 5% + >>> s1 = pyEQL.Solution({'Mg+2': '0.3 mol/kg', 'Cl-': '0.6 mol/kg'},temperature='30 degC') + >>> s1.get_osmotic_coefficient() + - 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, [])]) + References: + [may11] - # 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: - if new_el_amount == 0 and orig_el_amount > 0: - 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 - } - ) - elif abs(orig_el_amount - new_el_amount) / orig_el_amount > _rtol: - logger.error( - f"PHREEQC returned a total Element {el} concentration of {new_el_amount} mol, " - f"which differs from the original concentration of {orig_el_amount}. This " - "should never occur and indicates an error in the PHREEQC database or calculation." - ) + May, P. M., Rowland, D., Hefter, G., & Königsberger, E. (2011). + A Generic and Updatable Pitzer Characterization of Aqueous Binary Electrolyte Solutions at 1 bar + and 25 °C. Journal of Chemical & Engineering Data, 56(12), 5066-5077. doi:10.1021/je2009329 - # 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() + [rbs68] - # 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 + Robinson, R. A.; Stokes, R. H. Electrolyte Solutions: Second Revised Edition; Butterworths: London, 1968, + p.32. - def __deepcopy__(self, memo) -> "NativeEOS": - # custom deepcopy required because the PhreeqPython instance used by the Native and Phreeqc engines - # is not pickle-able. + [mistry13] - cls = self.__class__ - result = cls.__new__(cls) - memo[id(self)] = result - for k, v in self.__dict__.items(): - if k == "pp": - result.pp = PhreeqPython(database=self.phreeqc_db, database_directory=self.db_path) - continue - setattr(result, k, copy.deepcopy(v, memo)) - return result + Mistry, K. H.; Hunter, H. a.; Lienhard V, J. H. Effect of composition and nonideal solution + behavior on desalination calculations for mixed electrolyte solutions with comparison to + seawater. Desalination 2013, 318, 34-47. + """ + ionic_strength = solution.ionic_strength -class Phreeqc2026EOS(NativeEOS): - """Engine based on the PhreeqC model, as implemented in the pyphreeqc - module of pyEQL.""" + effective_osmotic_sum = 0 + molality_sum = 0 - def __init__( - self, - phreeqc_db: Literal[ - "phreeqc.dat", "vitens.dat", "wateq4f_PWN.dat", "pitzer.dat", "llnl.dat", "geothermal.dat" - ] = "phreeqc.dat", - ) -> None: - """ - Args: - phreeqc_db: Name of the PHREEQC database file to use for solution thermodynamics - and speciation calculations. Generally speaking, `llnl.dat` is recommended - for moderate salinity water and prediction of mineral solubilities, - `wateq4f_PWN.dat` is recommended for low to moderate salinity waters. It is - similar to vitens.dat but has many more species. `pitzer.dat` is recommended - when accurate activity coefficients in solutions above 1 M TDS are desired, but - it has fewer species than the other databases. `llnl.dat` and `geothermal.dat` - may offer improved prediction of LSI but currently these databases are not - usable because they do not allow for conductivity calculations. - """ + # loop through all the salts in the solution, calculate the osmotic + # coefficint for each, and average them into an effective osmotic + # coefficient + for d in solution.get_salt_dict().values(): + item = d["salt"] + # determine alpha1 and alpha2 based on the type of salt + # see the May reference for the rules used to determine + # alpha1 and alpha2 based on charge + if item.z_cation >= 2 and item.z_anion <= -2: + if item.z_cation >= 3 or item.z_anion <= -3: + alpha1 = 2.0 + alpha2 = 50.0 + else: + alpha1 = 1.4 + alpha2 = 12.0 + else: + alpha1 = 2.0 + alpha2 = 0 - from pyEQL.phreeqc import IS_AVAILABLE, Phreeqc # noqa: PLC0415 + # set the concentration as the average concentration of the cation and + # anion in the salt, accounting for stoichiometry + # concentration = (solution.get_amount(Salt.cation,'mol/kg')/Salt.nu_cation + \ + # solution.get_amount(Salt.anion,'mol/kg')/Salt.nu_anion)/2 - if not IS_AVAILABLE: - raise RuntimeError("pyEQL phreeqc support is not available in this installation") + # get the effective molality of the salt + concentration = ureg.Quantity(d["mol"], "mol") / solution.solvent_mass - self.phreeqc_db = phreeqc_db - # database files in this list are not distributed with phreeqpython - self.db_path = ( - Path(os.path.dirname(__file__)) / "database" if self.phreeqc_db in ["llnl.dat", "geothermal.dat"] else None - ) - # create the PhreeqcPython instance - self.pp = Phreeqc(database=self.phreeqc_db, database_directory=self.db_path) - # attributes to hold the PhreeqPython solution. - self.ppsol = None - # store the solution composition to see whether we need to re-instantiate the solution - self._stored_comp = None + molality_sum += concentration - def _setup_ppsol(self, solution: "solution.Solution") -> None: - """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. + param = solution.get_property(item.formula, "model_parameters.activity_pitzer") + if param is not None: + osmotic_coefficient = ac.get_osmotic_coefficient_pitzer( + ionic_strength, + concentration, + alpha1, + alpha2, + ureg.Quantity(param["Beta0"]["value"]).magnitude, + ureg.Quantity(param["Beta1"]["value"]).magnitude, + ureg.Quantity(param["Beta2"]["value"]).magnitude, + ureg.Quantity(param["Cphi"]["value"]).magnitude, + item.z_cation, + item.z_anion, + item.nu_cation, + item.nu_anion, + str(solution.temperature), + ) - # 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 + logger.debug( + f"Calculated osmotic coefficient of water as {osmotic_coefficient} based on salt " + f"{item.formula} using Pitzer model" + ) + effective_osmotic_sum += concentration * osmotic_coefficient - 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" + logger.debug( + f"Returning unit osmotic coefficient for salt {item.formula} because Pitzer parameters are not" + "available in database." + ) + effective_osmotic_sum += concentration * 1 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: - 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: - """ - 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: - if self.ppsol is not None: - self.ppsol.forget() - self._setup_ppsol(solution) - - osmotic = self.ppsol.get_osmotic_coefficient() - if osmotic == 0: - # PHREEQC returns 0 when it assumes a unit osmotic coefficient, so we need to catch that case and return 1 instead. - osmotic = 1 - return ureg.Quantity(osmotic, "dimensionless") + return effective_osmotic_sum / molality_sum + except ZeroDivisionError: + # this means the solution is empty + return ureg.Quantity(1.0) 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") - - -class PhreeqcEOS(Phreeqc2026EOS): - """Engine based on the PhreeqC model, as implemented via the phreeqpython package.""" - - def __init__( - self, - phreeqc_db: Literal[ - "phreeqc.dat", "vitens.dat", "wateq4f_PWN.dat", "pitzer.dat", "llnl.dat", "geothermal.dat" - ] = "phreeqc.dat", - ) -> None: - """ - Args: - phreeqc_db: Name of the PHREEQC database file to use for solution thermodynamics - and speciation calculations. Generally speaking, `llnl.dat` is recommended - for moderate salinity water and prediction of mineral solubilities, - `wateq4f_PWN.dat` is recommended for low to moderate salinity waters. It is - similar to vitens.dat but has many more species. `pitzer.dat` is recommended - when accurate activity coefficients in solutions above 1 M TDS are desired, but - it has fewer species than the other databases. `llnl.dat` and `geothermal.dat` - may offer improved prediction of LSI but currently these databases are not - usable because they do not allow for conductivity calculations. - """ - self.phreeqc_db = phreeqc_db - # database files in this list are not distributed with phreeqpython - self.db_path = ( - Path(os.path.dirname(__file__)) / "database" if self.phreeqc_db in ["llnl.dat", "geothermal.dat"] else None - ) - # create the PhreeqcPython instance - # try/except added to catch unsupported architectures, such as Apple Silicon - try: - self.pp = PhreeqPython(database=self.phreeqc_db, database_directory=self.db_path) - except OSError: - logger.error( - "OSError encountered when trying to instantiate phreeqpython. Most likely this means you" - " are running on an architecture that is not supported by PHREEQC, such as Apple M1/M2 chips." - " pyEQL will work, but equilibrate() will have no effect." - ) - # attributes to hold the PhreeqPython solution. - self.ppsol = None - # store the solution composition to see whether we need to re-instantiate the solution - self._stored_comp = None - - def _setup_ppsol(self, solution: "solution.Solution") -> None: - """Helper method to set up a PhreeqPython solution for subsequent analysis.""" - 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" + # identify the predominant salt in the solution + salt = solution.get_salt() + solute_vol = ureg.Quantity(0, "L") - # 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. + # use the pitzer approach if parameters are available + pitzer_calc = False + param = None if salt is None else solution.get_property(salt.formula, "model_parameters.molar_volume_pitzer") - # 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 param is not None: + # determine the average molality of the salt + # this is necessary for solutions inside e.g. an ion exchange + # membrane, where the cation and anion concentrations may be + # unequal + molality = (solution.get_amount(salt.cation, "mol/kg") + solution.get_amount(salt.anion, "mol/kg")) / 2 - 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) + # determine alpha1 and alpha2 based on the type of salt + # see the May reference for the rules used to determine + # alpha1 and alpha2 based on charge + if salt.z_cation >= 2 and salt.z_anion <= -2: + if salt.z_cation >= 3 or salt.z_anion <= -3: + alpha1 = 2.0 + alpha2 = 50.0 + else: + alpha1 = 1.4 + alpha2 = 12 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" + alpha1 = 2.0 + alpha2 = 0.0 - # create the PHREEQC solution object - try: - ppsol = self.pp.add_solution(d) - except Exception as e: - print(d) - # catch problems with the input to phreeqc - raise ValueError( - "There is a problem with your input. The error message received from " - f" phreeqpython is:\n\n {e}\n Check your input arguments, especially " - "the composition dictionary, and try again." + apparent_vol = ac.get_apparent_volume_pitzer( + solution.ionic_strength, + molality, + alpha1, + alpha2, + ureg.Quantity(param["Beta0"]["value"]).magnitude, + ureg.Quantity(param["Beta1"]["value"]).magnitude, + ureg.Quantity(param["Beta2"]["value"]).magnitude, + ureg.Quantity(param["Cphi"]["value"]).magnitude, + ureg.Quantity(param["V_o"]["value"]).magnitude, + salt.z_cation, + salt.z_anion, + salt.nu_cation, + salt.nu_anion, + str(solution.temperature), ) - self.ppsol = ppsol - - def _destroy_ppsol(self) -> None: - """Remove the PhreeqPython solution from memory.""" - if self.ppsol is not None: - self.ppsol.forget() - self.ppsol = None - - def get_activity_coefficient(self, solution: "solution.Solution", solute: str) -> ureg.Quantity: - """ - 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) + solute_vol += ( + apparent_vol + * ( + solution.get_amount(salt.cation, "mol") / salt.nu_cation + + solution.get_amount(salt.anion, "mol") / salt.nu_anion + ) + / 2 + ) - # 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 + pitzer_calc = True - # calculate the molal scale activity coefficient - # act = self.ppsol.activity(k, "mol") / self.ppsol.molality(k, "mol") - act = self.ppsol.pp.ip.get_activity(self.ppsol.number, k) / self.ppsol.pp.ip.get_molality(self.ppsol.number, k) + logger.debug(f"Updated solution volume using Pitzer model for solute {salt.formula}") - return ureg.Quantity(act, "dimensionless") + # add the partial molar volume of any other solutes, except for water + # or the parent salt, which is already accounted for by the Pitzer parameters + for solute, mol in solution.components.items(): + # ignore water + if solute in ["H2O", "HOH", "H2O(aq)"]: + continue - def get_osmotic_coefficient(self, solution: "solution.Solution") -> ureg.Quantity: - """ - Return the *molal scale* osmotic coefficient of solute, given a Solution - object. + # ignore the salt cation and anion, if already accounted for by Pitzer + if pitzer_calc is True and solute in [salt.anion, salt.cation]: + continue - PHREEQC appears to assume a unit osmotic coefficient unless the pitzer database - is used. Unfortunately, there is no easy way to access the osmotic coefficient - via phreeqcpython - """ - # TODO - find a way to access or calculate osmotic coefficient - return ureg.Quantity(1, "dimensionless") + part_vol = solution.get_property(solute, "size.molar_volume") + if part_vol is not None: + solute_vol += part_vol * ureg.Quantity(mol, "mol") + logger.debug(f"Updated solution volume using direct partial molar volume for solute {solute}") - def __deepcopy__(self, memo) -> "NativeEOS": - # custom deepcopy required because the PhreeqPython instance used by the Native and Phreeqc engines - # is not pickle-able. + else: + logger.warning( + f"Volume of solute {solute} will be ignored because partial molar volume data are not available." + ) - cls = self.__class__ - result = cls.__new__(cls) - memo[id(self)] = result - for k, v in self.__dict__.items(): - if k == "pp": - result.pp = PhreeqPython(database=self.phreeqc_db, database_directory=self.db_path) - continue - setattr(result, k, copy.deepcopy(v, memo)) - return result + return solute_vol.to("L") From 32c3aaa204fcefd7f9cb19712df9079c1e272403 Mon Sep 17 00:00:00 2001 From: Ryan Kingsbury Date: Sat, 14 Feb 2026 12:45:17 -0500 Subject: [PATCH 3/3] change typing import for python 3.10 --- src/pyEQL/engines.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/pyEQL/engines.py b/src/pyEQL/engines.py index ec9ccf49..1b6fbc56 100644 --- a/src/pyEQL/engines.py +++ b/src/pyEQL/engines.py @@ -13,9 +13,10 @@ import warnings from abc import ABC, abstractmethod from pathlib import Path -from typing import TYPE_CHECKING, Literal, Self +from typing import TYPE_CHECKING, Literal from phreeqpython import PhreeqPython +from typing_extensions import Self import pyEQL.activity_correction as ac from pyEQL import ureg