Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
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
61 changes: 61 additions & 0 deletions src/pyEQL/solution.py
Original file line number Diff line number Diff line change
Expand Up @@ -11,6 +11,7 @@
import logging
import os
import warnings
import re
from functools import lru_cache
from importlib.resources import files
from pathlib import Path
Expand Down Expand Up @@ -2523,6 +2524,66 @@
return Solution(**solution_dict)
return loadfn(filename)



@classmethod
def from_phreeqc(cls, path_to_pqi: str) -> "Solution":
"""
Parse a PHREEQC .pqi input file and pull out only the SOLUTION block,
extracting solute amounts (assumed mol/L), plus pH, pE, and temperature if present.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I believe the default units (if the user doesn't specify anything) are mmol/kgw, so at a minimum you should assume that instead of mol/L. But it would be even better to parse the units and water keywords as well so that this method can detect when users specify different units. See

https://wwwbrr.cr.usgs.gov/projects/GWC_coupled/phreeqc/html/final-56.html#89789

for a description of how they work.

"""
import re

sol_lines: list[str] = []
in_solution = False

with open(path_to_pqi, "r") as f:
for raw in f:
line = raw.strip()
if not in_solution and line.upper().startswith("SOLUTION"):
in_solution = True
continue
# stop on blank or a new ALL-CAPS keyword
if in_solution and (line == "" or re.match(r"^[A-Z]+\b", line)):
break
if in_solution:
sol_lines.append(line)

if not sol_lines:
raise ValueError(f"No SOLUTION block found in '{path_to_pqi}'")

composition: dict[str, str] = {}
pH = pE = None
temperature = None # degrees C

for l in sol_lines:
parts = re.split(r"\s+", l)
key = parts[0].upper()
vals = parts[1:]
if key == "PH":
pH = float(vals[0])
elif key == "PE":
pE = float(vals[0])

Check warning on line 2566 in src/pyEQL/solution.py

View check run for this annotation

Codecov / codecov/patch

src/pyEQL/solution.py#L2566

Added line #L2566 was not covered by tests
elif key in ("TEMP", "TEMPERATURE"):
temperature = float(vals[0])
else:
# solute line: e.g. "Ca+2 0.01"
try:
amt = float(vals[0])
except ValueError:
continue

Check warning on line 2574 in src/pyEQL/solution.py

View check run for this annotation

Codecov / codecov/patch

src/pyEQL/solution.py#L2573-L2574

Added lines #L2573 - L2574 were not covered by tests
# **here** we tag the unit so pyEQL knows mol/L
composition[parts[0]] = f"{amt} mol/L"

return cls(
solutes=composition,
pH=pH,
pE=pE,
temperature=f"{temperature} degC" if temperature is not None else None,
engine="ideal",
)


# arithmetic operations
def __add__(self, other: Solution) -> Solution:
"""
Expand Down
32 changes: 32 additions & 0 deletions tests/test_parser_from_phreeqc.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,32 @@
import warnings
warnings.filterwarnings("ignore", category=RuntimeWarning)

import textwrap
import pytest

from pyEQL import Solution

def test_from_phreeqc_minimal_parser(tmp_path):
"""Smoke-test our new from_phreeqc() on a tiny SOLUTION block."""
pqi = tmp_path / "sample.pqi"
pqi.write_text(textwrap.dedent("""
SOLUTION foo
pH 7.5
temp 30
Ca+2 0.01
SO4-2 0.01
Comment on lines +16 to +17
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

One challenge with this is that PHREEQC inputs are normally specified by elements rather than specific ionic species. I have not tested (yet), but I'm pretty sure that the simple input file you create for this test is not actually a valid PHREEQC input, because it lists individual species rather than elements.

A typical PHREEQC concentration input would like, for example:

Ca        80.
S(6)      96.     as SO4
N(5) N(3) 14.     as N

See https://wwwbrr.cr.usgs.gov/projects/GWC_coupled/phreeqc/html/final-56.html#89789 for other examples.

In the above, the number in parentheses indicates the oxidation state. Note that some elements list "as" which tells you what species they should be interpreted as, but not all. For example "Ca" is just "Ca". Would it be possible to detect "as" keywords and then use whatever species comes after the "as"?

For elements that do not have "as", PHREEQC usually assigns a reasonable default species such as interpreting Ca as Ca+2). However, there is a related bug - #229 that we need to address when a bare element like Ca is used as input to pyEQL.

Can you investigate how to fix that bug? It could be done in a separate PR if you want, but it's necessary in order to get PHREEQC input to parse correctly.

END
"""))
sol = Solution.from_phreeqc(str(pqi))
assert sol.pH == pytest.approx(7.5, rel=1e-6)
assert sol.temperature.magnitude == pytest.approx(303.15, rel=1e-6)
assert sol.get_amount("Ca+2", "mol/L").magnitude == pytest.approx(0.01, rel=1e-6)
assert sol.get_amount("SO4-2", "mol/L").magnitude == pytest.approx(0.01, rel=1e-6)

def test_from_phreeqc_missing_block_raises(tmp_path):
"""If there is no SOLUTION block, we should get a ValueError."""
pqi = tmp_path / "empty.pqi"
pqi.write_text("THIS FILE HAS NO SOLUTION BLOCK\n")
with pytest.raises(ValueError):
Solution.from_phreeqc(str(pqi))
Loading