Adding tests for solid-liquid equilibrium in equilibrate() #295

SuixiongTay · 2025-11-24T22:32:13Z

Summary

Added pytest test cases
Addressing enable solid-liquid and gas-liquid equilibrium calculations in equilibrate() #292 .

This PR revises the new test functions to tests/test_phreeqc.py:

package, installed with the "phreeqc" extra.

…hon incompatibility)

…cture

SuixiongTay · 2025-11-24T22:52:10Z

Hi @vineetbansal , I ran the pyEQL pytest based on the changes above and didn't encounter any errors. However, I noticed there were errors related to the Cl for [testing / docs (pull_request)] and [testing / test (3.12, ubuntu-latest) (pull_request)] checks in #294. Is there a way to run these Cl tests locally before merging into vb/issue186a?

SuixiongTay · 2025-11-25T21:02:42Z

The main issue is related to the saturation index should default to 0 (SI = 0). This is ready to merge after review. Thanks!

Example: Calcite 0 100 instead of Calcite 1 100

rkingsbury · 2025-11-25T22:11:54Z

tests/test_phreeqc.py

+    total_alk = HCO3 + 2 * CO3 + OH - H
+    assert alk.to("mg/L").magnitude == pytest.approx(total_alk, abs=0.01)


Here it appears you're comparing alkalinity calculated in mol/L with the value returned from the wrapper converted to mg/L. Why do you expect these to be equal?

Also the tolerance here is 0.01, but the total amount of CO2 in the starting solution is only 0.001 mol/L, so I think the test might be passing just because the tolerance is so high.

When running the following test:

def test_alkalinity(): solution = Solution({"CO2(aq)": "0.001 mol/L"}, pH=7, volume="1 L", engine="phreeqc") solution.equilibrate() alk = solution.alkalinity # Total alkalinity HCO3 = solution.get_amount("HCO3-", "mg/L").magnitude CO3 = solution.get_amount("CO3-2", "mg/L").magnitude OH = solution.get_amount("OH-", "mg/L").magnitude H = solution.get_amount("H+", "mg/L").magnitude # Alkalinity calculated from the excess of negative charges from weak acids total_alk = HCO3 + 2 * CO3 + OH - H assert alk.to("mg/L").magnitude == pytest.approx(total_alk, abs=0.001)

I encountered the following error:

============= short test summary info ========================= FAILED test_phreeqc.py::test_alkalinity - assert 0.0 == 50.05633916820964 ± 1.0e-03

The alkalinity units appears to be handling properly in mg/L but yields an incorrect value.

Upon inspecting def alkalinity(self), it looks like the carbonate species (HCO3-, CO3-2) are not included in the calculation, which explains why the method returns 0 mg/L. Would it make sense to create a separate PR to add carbonate alkalinity before we proceed with the test?

See my comment on this below - please mark this test an expected failure for now, and open an Issue with the content of this post. This results from a discrepancy in how pyEQL defines alkalinity, so it doesn't really indicate a problem with equilibrate, but it is something we should look into.

rkingsbury · 2025-11-25T22:13:20Z

tests/test_phreeqc.py

+def test_equilibrate_liquid():
    solution = Solution([["Cu+2", "4 mol/L"], ["O-2", "4 mol/L"]], volume="2 L", engine="phreeqc")
-    # Specifying an unrecognized solid raises an Exception
-    with pytest.raises(Exception):  # noqa: B017, PT011
-        solution.equilibrate(solids=["Ferroxite"])
+    solution.equilibrate()
+    assert solution.get_total_amount("Cu", "mol").magnitude == pytest.approx(1.654340558452752)
+


Please change this test to Ca2+ instead of copper, and use a lower concentration like 0.05 M so that the Ca doesn't precipitate. That way the expected result (0.05 M * 2 L) will be much clearer

tests/test_phreeqc.py

…nto pytest_vb186

Phreeqc wrapper

This reverts commit 9f0bc8e, reversing changes made to 0c01935.

rkingsbury

@SuixiongTay and @YitongPan1 , thank you for your work on these tests. With these recent changes they are getting close to ready - please both of you read the comments I've made, ask any technical questions especially if the rationale doesn't make sense to you, and do your best to address.

rkingsbury · 2025-12-01T18:10:18Z

tests/test_phreeqc.py

+def test_equilibrate_water_pH7():
+    solution = Solution({}, pH=7.00, temperature="25 degC", volume="1 L", engine="phreeqc")
+    solution.equilibrate()
+    # pH = -log10[H+]
+    assert solution.get_amount("H+", "mol").magnitude == pytest.approx(1.0e-07 * 0.997, rel=0.01)
+    # 14 - pH = log10[OH-]
+    assert solution.get_amount("OH-", "mol").magnitude == pytest.approx(1.01e-07 * 0.997, rel=0.01)
+    # small amount of H2 gas
+    assert solution.get_amount("H2", "mol").magnitude == pytest.approx(0 * 0.997, rel=0.01)
+    # 55.5 mol/kg * 0.997 kg = total moles of water
+    assert solution.get_amount("H2O", "mol").magnitude == pytest.approx(55.5 * 0.997, rel=0.01)


Good. Two comments here:

Since pH = -log10[H+], pH 7 equals 1e-7 mol/L, and you have set the volume of the solution as 1 L. Therefore, I think the mol of H+ should be exactly 1e-7; you don't need to adjust for the density (0.997).

I think your test tolerance is too large (too easy to pass). rel=0.01 means that the actual and expected results can differ by 1%. That's ok in some cases (like comparing a model to experiment), but here we know exactly what the answer should be. Because of that, I'd prefer to see an absolute rather than a relative tolerance, and I would think at least 10x smaller than your expected result (so, atol=1e-8) or lower would be appropriate.

If the test won't pass with that tolerance, we can look at why and decide whether it's numerical noise, something unexpected that equlibrate is doing, etc.

Hi @rkingsbury ,
For this case, the output of PHREEQC is in mol/kg rather than mol, so a unit conversion may still be needed to make this part work correctly.
We have also noticed a difference in the density values reported by pyEQL and PHREEQC: for this solution, pyEQL gives a density of 0.99704 kg/L, while PHREEQC reports 0.99757 kg/L. If the tolerance is set too low, this difference could cause the tests to fail.

@YitongPan1, @SuixiongTay - for the simple phreeqc script:

SOLUTION 1 temp 25

I'm seeing the output (among other lines):

----------------------------Description of solution---------------------------- pH = 7.000 pe = 4.000 Specific Conductance (uS/cm, 25oC) = 0 Density (g/cm3) = 0.99704

So phreeqc's output matches pyEQL's (as it should). Let me know where you're seeing 0.99757 kg/L.

Given that this is the case, I'm wondering if its okay to change this test case to its simplest form (which does pass):

def test_equilibrate_water_pH7(): solution = Solution({}, pH=7.00, temperature="25 degC", volume="1 L", engine="phreeqc") solution.equilibrate() # pH = -log10[H+] assert np.isclose(solution.get_amount("H+", "mol/kg").magnitude, 1e-07) # # 14 - pH = log10[OH-] assert np.isclose(solution.get_amount("OH-", "mol/kg").magnitude, 1e-07) # # small amount of H2 gas assert solution.get_amount("H2", "mol/kg").magnitude > 0 # Density of H2O in phreeqc = 0.99704 kg/L # For 0.99704 kg H2O # mol = 0.99704 kg * 1000 g/kg / (18.01528 g/mol) = 55.34413009400909 # mol/kg = 55.34413009400909 / 0.99704 = 55.50843506179199 assert np.isclose(solution.get_amount("H2O", "mol/kg").magnitude, 55.50843506179199)

The default tolerance for np.isclose is 1e-8 anyway so I think its better to leave it out.

Yes I'm good with your changes here @vineetbansal

tests/test_phreeqc.py

rkingsbury · 2025-12-01T18:17:26Z

tests/test_phreeqc.py

+def test_equilibrate_FeO3H3_ppt():
+    solution = Solution({"Fe+3": "0.01 mol/L", "OH-": "10**-7 mol/L"}, volume="1 L", engine="phreeqc")
+    solution.equilibrate()
+    # Under supersaturation, Fe(OH)3 precipitates out with Ksp = 10**-38.8
+    Fe3_act_coef = solution.get_activity_coefficient("Fe+3")
+    OH_act_coef = solution.get_activity_coefficient("OH-")
+    Fe3_conc = solution.get_amount("Fe+3", "mol/L").magnitude
+    OH_conc = solution.get_amount("OH-", "mol/L").magnitude
+    SI_FeO3H3 = np.log10((Fe3_act_coef.magnitude * Fe3_conc) * (OH_act_coef.magnitude * OH_conc) ** 3 / 10**-38.8)
+    assert solution.engine.ppsol.si("Fe(OH)3(a)") > 0
+    assert solution.engine.ppsol.si("Fe(OH)3(a)") == pytest.approx(SI_FeO3H3, rel=0.5)


I like the thought process here. To streamline this test:

You can use get_activity('Fe_3') to directly get the activity instead of having to separately get the coefficient and the concentration

In the assert statements, since you are already testing the si numerically (2nd test), there's no reason to also test that it's > 0. Please delete the first assert.

The existing assert checks the underlying engine (.engine.ppsol.si) but not the state of the Solution instance itself. Add an assert that check the total Fe concentration in the solution is less than 0.01 mol/L (indicating that some has precipitated).

The saturation index calculation is somewhat tricky as it deviates from the typical SI = log(IAP/pKsp) calculation. We have included the assertion assert solution.engine.ppsol.si("Fe(OH)3(a)") > 0 to check the expected behavior for now.

tests/test_phreeqc.py

rkingsbury · 2025-12-01T18:29:05Z

tests/test_phreeqc.py

+    assert solution.get_total_amount("Cu", "mol").magnitude == pytest.approx(0.001 * 2, rel=0.01)
+    assert solution.get_total_amount("N(0)", "mol").magnitude == pytest.approx(0)


Please add asserts that also verify both Oxygen and CO2(g) are in solution. Preferably, compute the expected concentrations from the respective Henry's law constants and test each quantitatively.

…nto pytest_vb186

rkingsbury · 2025-12-11T03:34:37Z

Closing to continue discussion in #292

vineetbansal and others added 12 commits November 21, 2025 12:22

Integrated the pypyreeqc development done so far, as the pyEQL-phreeqc

5a791de

package, installed with the "phreeqc" extra.

Merge branch 'main' into vb/pyphreeqc

42c5995

bumped version on wrapper to trigger CI

13fdd22

skipping some tests that should be skipped on macos anyway (phreeqpyt…

75cacef

…hon incompatibility)

version bump

76157fc

linkcheck ignore

e3fa93a

not looking for specific floats in log/dump strings, but overall stru…

da2a76d

…cture

tox build_ext step, invoked in workflow whenever build is invoked

a7eed3a

relaxed test for log strings

5540bd1

revised pytest

c41606d

revised pytest

e24c95e

revised pytest

f71e3d2

revised pytest

b2b1fcc

SuixiongTay requested review from YitongPan1, rkingsbury and vineetbansal and removed request for YitongPan1 November 25, 2025 13:58

SuixiongTay changed the title ~~'PHREEQC` pytest revision~~ PHREEQC pytest revision Nov 25, 2025

rkingsbury requested changes Nov 25, 2025

View reviewed changes

rkingsbury and others added 5 commits November 25, 2025 17:26

Merge branch 'vb/issue186a' into pytest_vb186

dba6dfc

Add/fix tests

0a2fab0

Merge branch 'pytest_vb186' of https://github.com/SuixiongTay/pyEQL i…

877b20f

…nto pytest_vb186

Merge pull request KingsburyLab#291 from KingsburyLab/vb/pyphreeqc

4dbcfe8

Phreeqc wrapper

Merge branch 'vb/issue186a' into pytest_vb186

0c01935

rkingsbury changed the base branch from vb/issue186a to main December 1, 2025 17:52

Merge branch 'main' into pytest_vb186

9f0bc8e

rkingsbury changed the base branch from main to vb/issue186a December 1, 2025 17:53

Revert "Merge branch 'main' into pytest_vb186"

0711890

This reverts commit 9f0bc8e, reversing changes made to 0c01935.

rkingsbury requested changes Dec 1, 2025

View reviewed changes

rkingsbury changed the title ~~PHREEQC pytest revision~~ Adding tests for solid-liquid equilibrium in equilibrate() Dec 1, 2025

Sui Xiong Tay added 2 commits December 10, 2025 01:36

pytest revisions

8a99f97

Merge branch 'pytest_vb186' of https://github.com/SuixiongTay/pyEQL i…

1d66ee1

…nto pytest_vb186

rkingsbury closed this Dec 11, 2025

SuixiongTay mentioned this pull request Dec 11, 2025

pyEQL alkalinity incorrect for pure carbonate solutions #299

Open

		total_alk = HCO3 + 2 * CO3 + OH - H
		assert alk.to("mg/L").magnitude == pytest.approx(total_alk, abs=0.01)

		assert solution.get_total_amount("Cu", "mol").magnitude == pytest.approx(0.001 * 2, rel=0.01)
		assert solution.get_total_amount("N(0)", "mol").magnitude == pytest.approx(0)

Adding tests for solid-liquid equilibrium in equilibrate() #295

Adding tests for solid-liquid equilibrium in equilibrate() #295

Uh oh!

Conversation

SuixiongTay commented Nov 24, 2025

Summary

Uh oh!

SuixiongTay commented Nov 24, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

SuixiongTay commented Nov 25, 2025

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

rkingsbury left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

rkingsbury commented Dec 11, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

4 participants

SuixiongTay commented Nov 24, 2025 •

edited

Loading

rkingsbury commented Dec 11, 2025 •

edited

Loading