From 0695e6c159bb5c9e03b3532cf6df417119f64227 Mon Sep 17 00:00:00 2001 From: Vineet Bansal Date: Fri, 13 Feb 2026 14:01:53 -0500 Subject: [PATCH 1/7] renamed 2 parameters in wrapper equalize(); copied methods that PhreeqcEOS inherits from NativeEOS into PhreeqcEOS itself --- src/pyEQL/engines.py | 214 +++++++++++++++++++++++++++++++++- src/pyEQL/phreeqc/solution.py | 4 +- 2 files changed, 215 insertions(+), 3 deletions(-) diff --git a/src/pyEQL/engines.py b/src/pyEQL/engines.py index 59fdaf90..3d19c520 100644 --- a/src/pyEQL/engines.py +++ b/src/pyEQL/engines.py @@ -831,6 +831,81 @@ def __init__( """ super().__init__(phreeqc_db=phreeqc_db) + 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" + + # 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" + + # 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." + ) + + 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 @@ -872,6 +947,143 @@ def get_solute_volume(self, solution: "solution.Solution") -> ureg.Quantity: # TODO - phreeqc seems to have no concept of volume, but it does calculate density return ureg.Quantity(0, "L") + def equilibrate( + self, + solution: "solution.Solution", + 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, 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 + + 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 + + def __deepcopy__(self, memo) -> "NativeEOS": + # custom deepcopy required because the PhreeqPython instance used by the Native and Phreeqc engines + # is not pickle-able. + + 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 Phreeqc2026EOS(EOS): """Engine based on the PhreeqC model, as implemented in the pyphreeqc @@ -1094,7 +1306,7 @@ def equilibrate( amounts = [v[1] for v in phases.values()] try: - self.ppsol.equalize(phases=phase_names, saturation_indices=saturation_indices, amounts=amounts) + 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): diff --git a/src/pyEQL/phreeqc/solution.py b/src/pyEQL/phreeqc/solution.py index b42c4d43..fb74aee9 100644 --- a/src/pyEQL/phreeqc/solution.py +++ b/src/pyEQL/phreeqc/solution.py @@ -48,9 +48,9 @@ def get_kgw(self) -> float: 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: + def equalize(self, phases: list[str], to_si: list[float], in_phase: list[float]) -> None: if self._phreeqc is not None: - self._phreeqc().equalize(self._number, phases, saturation_indices, amounts) + self._phreeqc().equalize(self._number, phases, to_si, in_phase) def si(self, eq_species) -> float: return self._get_calculated_prop("SI", eq_species=eq_species) From d81298f2db28657a72ea377850c86134f7bb0b37 Mon Sep 17 00:00:00 2001 From: Vineet Bansal Date: Fri, 13 Feb 2026 14:18:47 -0500 Subject: [PATCH 2/7] PhreeqcEOS now inherits directly from EOS; has its own constructor --- src/pyEQL/engines.py | 40 +++++++++++++++++++++++++++++----------- 1 file changed, 29 insertions(+), 11 deletions(-) diff --git a/src/pyEQL/engines.py b/src/pyEQL/engines.py index 3d19c520..6d9e8fc0 100644 --- a/src/pyEQL/engines.py +++ b/src/pyEQL/engines.py @@ -808,7 +808,7 @@ def __deepcopy__(self, memo) -> "NativeEOS": return result -class PhreeqcEOS(NativeEOS): +class PhreeqcEOS(EOS): """Engine based on the PhreeqC model, as implemented via the phreeqpython package.""" def __init__( @@ -819,17 +819,35 @@ def __init__( ) -> 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. + 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. """ - super().__init__(phreeqc_db=phreeqc_db) + 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.""" From fd04313aaafbb313858d59761f3ee8af1dd6fdb8 Mon Sep 17 00:00:00 2001 From: Vineet Bansal Date: Fri, 13 Feb 2026 14:19:58 -0500 Subject: [PATCH 3/7] changed order of class declaration --- src/pyEQL/engines.py | 132 +++++++++++++++++++++---------------------- 1 file changed, 66 insertions(+), 66 deletions(-) diff --git a/src/pyEQL/engines.py b/src/pyEQL/engines.py index 6d9e8fc0..b8435e90 100644 --- a/src/pyEQL/engines.py +++ b/src/pyEQL/engines.py @@ -808,8 +808,9 @@ def __deepcopy__(self, memo) -> "NativeEOS": return result -class PhreeqcEOS(EOS): - """Engine based on the PhreeqC model, as implemented via the phreeqpython package.""" +class Phreeqc2026EOS(EOS): + """Engine based on the PhreeqC model, as implemented in the pyphreeqc + module of pyEQL.""" def __init__( self, @@ -829,21 +830,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 @@ -851,6 +850,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 @@ -904,24 +906,21 @@ 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) -> ureg.Quantity: @@ -944,21 +943,20 @@ def get_activity_coefficient(self, solution: "solution.Solution", solute: str) - # 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) + 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: - """ - Return the *molal scale* osmotic coefficient of solute, given a Solution - object. + if self.ppsol is not None: + self.ppsol.forget() + self._setup_ppsol(solution) - 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") + 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") def get_solute_volume(self, solution: "solution.Solution") -> ureg.Quantity: """Return the volume of the solutes.""" @@ -1088,24 +1086,9 @@ def equilibrate( # call to equilibrate can thus result in a slight change in the Solution mass. solution.components[solution.solvent] = orig_solvent_moles - def __deepcopy__(self, memo) -> "NativeEOS": - # custom deepcopy required because the PhreeqPython instance used by the Native and Phreeqc engines - # is not pickle-able. - - 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 Phreeqc2026EOS(EOS): - """Engine based on the PhreeqC model, as implemented in the pyphreeqc - module of pyEQL.""" +class PhreeqcEOS(EOS): + """Engine based on the PhreeqC model, as implemented via the phreeqpython package.""" def __init__( self, @@ -1125,19 +1108,21 @@ 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 - self.pp = Phreeqc(database=self.phreeqc_db, database_directory=self.db_path) + # 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 @@ -1145,9 +1130,6 @@ 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 @@ -1201,21 +1183,24 @@ 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(PHRQSol(d)) + 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" phreeqc is:\n\n {e}\n Check your input arguments, especially " + f" phreeqpython 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.pp.remove_solution(0) # TODO: Are we only expecting a single solution per wrapper? + self.ppsol.forget() self.ppsol = None def get_activity_coefficient(self, solution: "solution.Solution", solute: str) -> ureg.Quantity: @@ -1238,20 +1223,21 @@ def get_activity_coefficient(self, solution: "solution.Solution", solute: str) - # 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) + act = self.ppsol.pp.ip.get_activity(self.ppsol.number, k) / self.ppsol.pp.ip.get_molality(self.ppsol.number, 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) + """ + Return the *molal scale* osmotic coefficient of solute, given a Solution + object. - 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") + 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") def get_solute_volume(self, solution: "solution.Solution") -> ureg.Quantity: """Return the volume of the solutes.""" @@ -1380,3 +1366,17 @@ def equilibrate( # 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 + + def __deepcopy__(self, memo) -> "NativeEOS": + # custom deepcopy required because the PhreeqPython instance used by the Native and Phreeqc engines + # is not pickle-able. + + 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 From 0980c9a187741dd096954138881da816f29e74fb Mon Sep 17 00:00:00 2001 From: Vineet Bansal Date: Fri, 13 Feb 2026 14:20:38 -0500 Subject: [PATCH 4/7] PhreeqcEOS now inherits from Phreeqc2026EOS --- 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 b8435e90..b21f6450 100644 --- a/src/pyEQL/engines.py +++ b/src/pyEQL/engines.py @@ -1087,7 +1087,7 @@ def equilibrate( solution.components[solution.solvent] = orig_solvent_moles -class PhreeqcEOS(EOS): +class PhreeqcEOS(Phreeqc2026EOS): """Engine based on the PhreeqC model, as implemented via the phreeqpython package.""" def __init__( From 93a2c4e56c617e726f582ae877940cd745a4a856 Mon Sep 17 00:00:00 2001 From: Vineet Bansal Date: Fri, 13 Feb 2026 14:21:40 -0500 Subject: [PATCH 5/7] equilibrate() in PhreeqcEOS gone (identical to parent) --- src/pyEQL/engines.py | 123 ------------------------------------------- 1 file changed, 123 deletions(-) diff --git a/src/pyEQL/engines.py b/src/pyEQL/engines.py index b21f6450..8993401c 100644 --- a/src/pyEQL/engines.py +++ b/src/pyEQL/engines.py @@ -1244,129 +1244,6 @@ def get_solute_volume(self, solution: "solution.Solution") -> ureg.Quantity: # TODO - phreeqc seems to have no concept of volume, but it does calculate density return ureg.Quantity(0, "L") - def equilibrate( - self, - solution: "solution.Solution", - 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, 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 - - 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 - def __deepcopy__(self, memo) -> "NativeEOS": # custom deepcopy required because the PhreeqPython instance used by the Native and Phreeqc engines # is not pickle-able. From 9d677770555852d35a7f5de9ac50e3e2e103cc15 Mon Sep 17 00:00:00 2001 From: Vineet Bansal Date: Fri, 13 Feb 2026 14:30:42 -0500 Subject: [PATCH 6/7] some methods from PhreeqcEOS gone (identical to its parent) --- src/pyEQL/engines.py | 11 ----------- 1 file changed, 11 deletions(-) diff --git a/src/pyEQL/engines.py b/src/pyEQL/engines.py index 8993401c..6c18018d 100644 --- a/src/pyEQL/engines.py +++ b/src/pyEQL/engines.py @@ -1197,12 +1197,6 @@ def _setup_ppsol(self, solution: "solution.Solution") -> None: 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 @@ -1239,11 +1233,6 @@ def get_osmotic_coefficient(self, solution: "solution.Solution") -> ureg.Quantit # TODO - find a way to access or calculate osmotic coefficient return ureg.Quantity(1, "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 __deepcopy__(self, memo) -> "NativeEOS": # custom deepcopy required because the PhreeqPython instance used by the Native and Phreeqc engines # is not pickle-able. From 806b87695fe591e11b943390b1a2c30dbb269ddb Mon Sep 17 00:00:00 2001 From: Vineet Bansal Date: Fri, 13 Feb 2026 14:52:18 -0500 Subject: [PATCH 7/7] Phreeqc2026EOS now descends from NativeEOS (their equilibrate methods are identical). PhreeqcEOS has its own _destroy_ppsol (identical to NativeEOS but differs slightly from Phreeqc2026EOS). --- src/pyEQL/engines.py | 131 +++---------------------------------------- 1 file changed, 7 insertions(+), 124 deletions(-) diff --git a/src/pyEQL/engines.py b/src/pyEQL/engines.py index 6c18018d..9ba9c338 100644 --- a/src/pyEQL/engines.py +++ b/src/pyEQL/engines.py @@ -808,7 +808,7 @@ def __deepcopy__(self, memo) -> "NativeEOS": return result -class Phreeqc2026EOS(EOS): +class Phreeqc2026EOS(NativeEOS): """Engine based on the PhreeqC model, as implemented in the pyphreeqc module of pyEQL.""" @@ -963,129 +963,6 @@ def get_solute_volume(self, solution: "solution.Solution") -> ureg.Quantity: # TODO - phreeqc seems to have no concept of volume, but it does calculate density return ureg.Quantity(0, "L") - def equilibrate( - self, - solution: "solution.Solution", - 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, 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 - - 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 - class PhreeqcEOS(Phreeqc2026EOS): """Engine based on the PhreeqC model, as implemented via the phreeqpython package.""" @@ -1197,6 +1074,12 @@ def _setup_ppsol(self, solution: "solution.Solution") -> None: 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