Skip to content

PhreeqcEOS inherits from Phreeqc2026EOS #351

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
rkingsbury merged 7 commits into from
Feb 14, 2026
Merged

PhreeqcEOS inherits from Phreeqc2026EOS #351

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
0695e6c
renamed 2 parameters in wrapper equalize(); copied methods that Phree…
vineetbansal Feb 13, 2026
d81298f
PhreeqcEOS now inherits directly from EOS; has its own constructor
vineetbansal Feb 13, 2026
fd04313
changed order of class declaration
vineetbansal Feb 13, 2026
0980c9a
PhreeqcEOS now inherits from Phreeqc2026EOS
vineetbansal Feb 13, 2026
93a2c4e
equilibrate() in PhreeqcEOS gone (identical to parent)
vineetbansal Feb 13, 2026
9d67777
some methods from PhreeqcEOS gone (identical to its parent)
vineetbansal Feb 13, 2026
806b876
Phreeqc2026EOS now descends from NativeEOS (their equilibrate methods…
vineetbansal Feb 13, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
317 changes: 148 additions & 169 deletions src/pyEQL/engines.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -808,72 +808,7 @@ def __deepcopy__(self, memo) -> "NativeEOS":
return result


class PhreeqcEOS(NativeEOS):
"""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.
"""
super().__init__(phreeqc_db=phreeqc_db)

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.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:
"""
Return the *molal scale* osmotic coefficient of solute, given a Solution
object.

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."""
# TODO - phreeqc seems to have no concept of volume, but it does calculate density
return ureg.Quantity(0, "L")


class Phreeqc2026EOS(EOS):
class Phreeqc2026EOS(NativeEOS):
"""Engine based on the PhreeqC model, as implemented in the pyphreeqc
module of pyEQL."""

Expand Down Expand Up @@ -1028,125 +963,169 @@ 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(

class PhreeqcEOS(Phreeqc2026EOS):
"""Engine based on the PhreeqC model, as implemented via the phreeqpython package."""

def __init__(
self,
solution: "solution.Solution",
atmosphere: bool = False,
solids: list[str] | None = None,
gases: dict[str, str | float] | None = None,
phreeqc_db: Literal[
"phreeqc.dat", "vitens.dat", "wateq4f_PWN.dat", "pitzer.dat", "llnl.dat", "geothermal.dat"
] = "phreeqc.dat",
) -> 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}
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.
"""
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]
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

# 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)
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"

# Use supplied gases, merged with atmospheric gases
# if atmosphere == True
gases = (ATMOSPHERE if atmosphere else {}) | (gases or {})
# 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.

# Mapping from phase name to:
# (<saturation_index>, <amount_in_moles>) tuples (for solids).
# (<log_partial_pressure>, <amount_in_moles>) 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))
# 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

for k, v in gases.items():
v_quantity = ureg.Quantity(v)
if v_quantity.dimensionless:
log_partial_pressure = v_quantity.magnitude
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:
log_partial_pressure = math.log10(v_quantity.to("atm").magnitude)
phases |= {f"{k}(g)": (log_partial_pressure, EQUILIBRIUM_PHASE_AMOUNT)}
d[key] = str(mol / solv_mass)

if phases:
phase_names = list(phases.keys())
saturation_indices = [v[0] for v in phases.values()]
amounts = [v[1] for v in phases.values()]
# tell PHREEQC which species to use for charge balance
if solution.balance_charge is not None and solution._cb_species in solution.get_components_by_element()[el]:
d[key] += " charge"

try:
self.ppsol.equalize(phases=phase_names, saturation_indices=saturation_indices, amounts=amounts)
except Exception as e:
# Re-raise exception due to unrecognized phases as ValueError.
if "Phase not found in database" in str(e):
raise ValueError(str(e))
raise e
# 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."
)

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
self.ppsol = ppsol

# 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."
)
def _destroy_ppsol(self) -> None:
"""Remove the PhreeqPython solution from memory."""
if self.ppsol is not None:
self.ppsol.forget()
self.ppsol = None

# 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
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)

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, [])])
# 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

# 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
}
)
# 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)

# 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()
return ureg.Quantity(act, "dimensionless")

# 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 get_osmotic_coefficient(self, solution: "solution.Solution") -> ureg.Quantity:
"""
Return the *molal scale* osmotic coefficient of solute, given a Solution
object.

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 __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
4 changes: 2 additions & 2 deletions src/pyEQL/phreeqc/solution.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
Loading