From d98352643e58ae66fcdf6e7406f55a4ffca832ed Mon Sep 17 00:00:00 2001 From: Vineet Bansal Date: Thu, 15 Jan 2026 17:05:42 -0500 Subject: [PATCH 1/2] Diffusion Coefficient --- src/pyEQL/engines.py | 17 ++++ src/pyEQL/pyEQL-phreeqc/pyEQL_phreeqc/core.py | 2 +- .../pyEQL-phreeqc/pyEQL_phreeqc/solution.py | 3 + src/pyEQL/solution.py | 10 ++- tests/phreeqc/test_phreeqc.py | 87 ++++++++++++++++++- 5 files changed, 114 insertions(+), 5 deletions(-) diff --git a/src/pyEQL/engines.py b/src/pyEQL/engines.py index 34a4b162..69f5d2c4 100644 --- a/src/pyEQL/engines.py +++ b/src/pyEQL/engines.py @@ -1128,3 +1128,20 @@ 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 get_diffusion_coefficient(self, solution: "solution.Solution", solute: str) -> ureg.Quantity: + 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 + + diffusion_coefficient = self.ppsol.get_diffusion_coefficient(k) + return ureg.Quantity(diffusion_coefficient, "meter ** 2 / s") diff --git a/src/pyEQL/pyEQL-phreeqc/pyEQL_phreeqc/core.py b/src/pyEQL/pyEQL-phreeqc/pyEQL_phreeqc/core.py index 763d4082..4101c887 100644 --- a/src/pyEQL/pyEQL-phreeqc/pyEQL_phreeqc/core.py +++ b/src/pyEQL/pyEQL-phreeqc/pyEQL_phreeqc/core.py @@ -13,7 +13,7 @@ "TOT['water']", "OSMOTIC", ) -SPECIES_PROPS = ("MOL", "ACT") +SPECIES_PROPS = ("MOL", "ACT", "DIFF_C") EQ_SPECIES_PROPS = ("SI",) diff --git a/src/pyEQL/pyEQL-phreeqc/pyEQL_phreeqc/solution.py b/src/pyEQL/pyEQL-phreeqc/pyEQL_phreeqc/solution.py index 6a8d02ed..a0591766 100644 --- a/src/pyEQL/pyEQL-phreeqc/pyEQL_phreeqc/solution.py +++ b/src/pyEQL/pyEQL-phreeqc/pyEQL_phreeqc/solution.py @@ -31,6 +31,9 @@ def get_activity(self, species) -> float: def get_molality(self, species) -> float: return self._get_calculated_prop("MOL", species=species) + def get_diffusion_coefficient(self, species) -> float: + return self._get_calculated_prop("DIFF_C", species=species) + def get_osmotic_coefficient(self) -> float: return self._get_calculated_prop("OSMOTIC") diff --git a/src/pyEQL/solution.py b/src/pyEQL/solution.py index 92b3b522..ed8689d6 100644 --- a/src/pyEQL/solution.py +++ b/src/pyEQL/solution.py @@ -2173,7 +2173,9 @@ def _get_molar_conductivity(self, solute: str) -> Quantity: return molar_cond.to("mS / cm / (mol/L)") - def _get_diffusion_coefficient(self, solute: str, activity_correction: bool = True) -> Quantity: + def _get_diffusion_coefficient( + self, solute: str, activity_correction: bool = True, use_engine: bool = False + ) -> Quantity: r""" Get the **temperature-adjusted** diffusion coefficient of a solute. @@ -2181,6 +2183,8 @@ def _get_diffusion_coefficient(self, solute: str, activity_correction: bool = Tr solute: the solute for which to retrieve the diffusion coefficient. activity_correction: If True (default), adjusts the diffusion coefficient for the effects of ionic strength using a model from Ref 2. + use_engine: Whether to use the underlying phreeqc engine to determine diffusion coefficient. + Only the 'pyeql' engine is supported. Notes: This method is equivalent to self.get_property(solute, "transport.diffusion_coefficient") @@ -2218,6 +2222,10 @@ def _get_diffusion_coefficient(self, solute: str, activity_correction: bool = Tr pyEQL.activity_correction._debye_parameter_activity """ + if use_engine: + assert self._engine == "pyeql", "Only supported for pyeql engine" + return self.engine.get_diffusion_coefficient(self, solute) + D = self.get_property(solute, "transport.diffusion_coefficient") rform = standardize_formula(solute) if D is None or D.magnitude == 0: diff --git a/tests/phreeqc/test_phreeqc.py b/tests/phreeqc/test_phreeqc.py index df2da2f2..f67d263f 100644 --- a/tests/phreeqc/test_phreeqc.py +++ b/tests/phreeqc/test_phreeqc.py @@ -2,9 +2,12 @@ from pathlib import Path import numpy as np +import pytest from pyEQL_phreeqc import Phreeqc, PHRQSol from pytest import approx +from pyEQL import Solution + def test_load_database_internal(): # `phreeqc.dat` is included with the package, so we don't need to specify @@ -406,7 +409,7 @@ def test_add_solution_input(): 5 PUNCH CELL_NO, TOT['water'], OSMOTIC, EOL_NOTAB$ 10 t = SYS("aq", count, name$, type$, moles) 20 FOR i = 1 to count - 30 PUNCH name$(i), MOL(name$(i)), ACT(name$(i)) + 30 PUNCH name$(i), MOL(name$(i)), ACT(name$(i)), DIFF_C(name$(i)) 40 NEXT i 50 PUNCH EOL$ 60 p = SYS("phases", count, name$, type$, moles) @@ -454,9 +457,9 @@ def test_add_solution_output(): # heading + soln 0 data + soln 1 data (regardless of whether we have -headings in USER_PUNCH or not) assert phreeqc.get_selected_output_row_count() == 3 # , , , "\n", - # +[, , ] repeated for each (4) species + # +[, , , ] repeated for each (4) species # +"\n" + [, ] repeated for each (3) equilibrium species - assert phreeqc.get_selected_output_column_count() == 23 + assert phreeqc.get_selected_output_column_count() == 27 def test_kgw(): @@ -596,3 +599,81 @@ def test_get_osmotic_coefficient(): osmotic_coefficient = phreeqc[0].get_osmotic_coefficient() assert osmotic_coefficient == approx(0.0) + + +@pytest.mark.xfail(reason="Diffusion Coefficient test discrepancies need to be investigated") +def test_diffusion_transport(): + """ + A copy of test_solution.py::test_diffusion_transport with minor syntactic changes. + """ + s1 = Solution(volume="2 L", engine="pyeql") + s2 = Solution([["Na+", "4 mol/L"], ["Cl-", "4 mol/L"]], volume="2 L", engine="pyeql") + + # test ionic strength adjustment + assert s1.get_diffusion_coefficient("H+", use_engine=True) > s2.get_diffusion_coefficient("H+", use_engine=True) + + # for Na+, d=122, a1=1.52, a2=3.7, A=1.173802/2.303 at 25 DegC, B = 3.2843078+10 + factor = np.exp( + -1.52 + * 1.173802 + / 2.303 + * 1 + * np.sqrt(s2.ionic_strength.magnitude) + / (1 + 3.2843078e10 * np.sqrt(s2.ionic_strength.magnitude) * 3.7 / (1 + s2.ionic_strength.magnitude**0.75)) + ) + assert np.isclose( + factor * s2.get_diffusion_coefficient("Na+", use_engine=True).magnitude, + s2.get_diffusion_coefficient("Na+", use_engine=True).magnitude, + atol=5e-11, + ) + s_dilute = Solution({"Na+": "1 mmol/L", "Cl-": "1 mmol/L"}, engine="pyeql") + assert np.isclose( + s_dilute.get_diffusion_coefficient("Na+", activity_correction=False, use_engine=True).magnitude, + 1.334e-9, + atol=1e-12, + ) + assert np.isclose(s_dilute.get_transport_number("Na+"), 0.396, atol=1e-3) + assert np.isclose(s_dilute.get_transport_number("Cl-"), 0.604, atol=1e-3) + + # test setting a default value + s2.default_diffusion_coeff = 0 + assert s2.get_diffusion_coefficient("Cs+", use_engine=True).magnitude == 0 + s2.default_diffusion_coeff = 1e-9 + assert s2.get_diffusion_coefficient("Cs+", activity_correction=False, use_engine=True).magnitude == 1e-9 + s2.default_diffusion_coeff = 0 + assert s2.get_diffusion_coefficient("Cs+", activity_correction=True, use_engine=True).magnitude < 1e-9 + d25 = s2.get_diffusion_coefficient("Na+", activity_correction=False, use_engine=True).magnitude + nu25 = s2.water_substance.nu + s2.temperature = "40 degC" + d40 = s2.get_diffusion_coefficient("Na+", activity_correction=False, use_engine=True).magnitude + nu40 = s2.water_substance.nu + assert np.isclose( + d40, + d25 * np.exp(122 / (273.15 + 40) - 122 / 298.15) * (nu25 / nu40), + atol=5e-11, + ) + + # test correction factors for concentration, as per Appelo 2017 Fig 5 + D1 = ( + Solution({"Na+": "1 umol/L", "Cl-": "1 umol/L"}, engine="pyeql") + .get_diffusion_coefficient("Na+", use_engine=True) + .magnitude + ) + D2 = ( + Solution({"Na+": "1.7 mol/kg", "Cl-": "1.7 mol/kg"}, engine="pyeql") + .get_diffusion_coefficient("Na+", use_engine=True) + .magnitude + ) + assert np.isclose(D2 / D1, 0.54, atol=1e-2) + + D1 = ( + Solution({"K+": "1 umol/L", "Cl-": "1 umol/L"}, engine="pyeql") + .get_diffusion_coefficient("K+", use_engine=True) + .magnitude + ) + D2 = ( + Solution({"K+": "0.5 mol/kg", "Cl-": "0.5 mol/kg"}, engine="pyeql") + .get_diffusion_coefficient("K+", use_engine=True) + .magnitude + ) + assert np.isclose(D2 / D1, 0.80, atol=1e-2) From e9e3f1268ca73510cf56046933294d547f7ca0ec Mon Sep 17 00:00:00 2001 From: Vineet Bansal Date: Fri, 16 Jan 2026 17:14:20 -0500 Subject: [PATCH 2/2] test for diffusion coefficient, introduced in a test module parallel to test_solution.py --- src/pyEQL/engines.py | 15 +-- .../pyEQL-phreeqc/pyEQL_phreeqc/solution.py | 6 +- src/pyEQL/solution.py | 20 ++-- tests/phreeqc/test_phreeqc.py | 81 ---------------- tests/test_solution_pyeql_engine.py | 92 +++++++++++++++++++ 5 files changed, 119 insertions(+), 95 deletions(-) create mode 100644 tests/test_solution_pyeql_engine.py diff --git a/src/pyEQL/engines.py b/src/pyEQL/engines.py index 69f5d2c4..461604fd 100644 --- a/src/pyEQL/engines.py +++ b/src/pyEQL/engines.py @@ -1136,12 +1136,15 @@ def get_diffusion_coefficient(self, solution: "solution.Solution", solute: str) # 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 "[" in k: + 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 diffusion_coefficient = self.ppsol.get_diffusion_coefficient(k) + if diffusion_coefficient is None: + return None return ureg.Quantity(diffusion_coefficient, "meter ** 2 / s") diff --git a/src/pyEQL/pyEQL-phreeqc/pyEQL_phreeqc/solution.py b/src/pyEQL/pyEQL-phreeqc/pyEQL_phreeqc/solution.py index a0591766..b242a477 100644 --- a/src/pyEQL/pyEQL-phreeqc/pyEQL_phreeqc/solution.py +++ b/src/pyEQL/pyEQL-phreeqc/pyEQL_phreeqc/solution.py @@ -20,9 +20,11 @@ def _get_calculated_props(self): def _get_calculated_prop(self, which, species: str | None = None, eq_species: str | None = None): if species is not None: - return self._calculated_props["species"][species][which] + species_props = self._calculated_props["species"].get(species) + return species_props[which] if species_props is not None else None if eq_species is not None: - return self._calculated_props["eq_species"][eq_species][which] + eq_species_props = self._calculated_props["eq_species"].get(species) + return eq_species_props[which] if eq_species_props is not None else None return self._calculated_props[which] def get_activity(self, species) -> float: diff --git a/src/pyEQL/solution.py b/src/pyEQL/solution.py index ed8689d6..f5905f06 100644 --- a/src/pyEQL/solution.py +++ b/src/pyEQL/solution.py @@ -2071,7 +2071,7 @@ def _get_property(self, solute: str, name: str) -> Any | None: return ureg.Quantity(val) return None - def get_transport_number(self, solute: str) -> Quantity: + def get_transport_number(self, solute: str, use_engine: bool = False) -> Quantity: r"""Calculate the transport number of the solute in the solution. Args: @@ -2114,7 +2114,7 @@ def get_transport_number(self, solute: str) -> Quantity: # cancels out # using species amounts in mol is equivalent to using concentrations in mol/L # since there is only one solution volume, and it's much faster. - term = self.get_molar_conductivity(item).magnitude * mol + term = self.get_molar_conductivity(item, use_engine=use_engine).magnitude * mol if item == solute: numerator = term @@ -2123,7 +2123,7 @@ def get_transport_number(self, solute: str) -> Quantity: return ureg.Quantity(numerator / denominator, "dimensionless") - def _get_molar_conductivity(self, solute: str) -> Quantity: + def _get_molar_conductivity(self, solute: str, use_engine: bool = False) -> Quantity: r""" Calculate the molar (equivalent) conductivity for a solute. @@ -2160,7 +2160,7 @@ def _get_molar_conductivity(self, solute: str) -> Quantity: See Also: :py:meth:`get_diffusion_coefficient` """ - D = self.get_diffusion_coefficient(solute) + D = self.get_diffusion_coefficient(solute, use_engine=use_engine) if D != 0: molar_cond = ( @@ -2222,12 +2222,20 @@ def _get_diffusion_coefficient( pyEQL.activity_correction._debye_parameter_activity """ + rform = standardize_formula(solute) + if use_engine: assert self._engine == "pyeql", "Only supported for pyeql engine" - return self.engine.get_diffusion_coefficient(self, solute) + D = self.engine.get_diffusion_coefficient(self, solute) + if D is None or D.magnitude == 0: + self.logger.warning( + f"Diffusion coefficient not found for species {rform}. Using default value of " + f"{self.default_diffusion_coeff} m**2/s." + ) + return ureg.Quantity(self.default_diffusion_coeff, "m**2/s") + return D D = self.get_property(solute, "transport.diffusion_coefficient") - rform = standardize_formula(solute) if D is None or D.magnitude == 0: self.logger.warning( f"Diffusion coefficient not found for species {rform}. Using default value of " diff --git a/tests/phreeqc/test_phreeqc.py b/tests/phreeqc/test_phreeqc.py index f67d263f..2fc716f0 100644 --- a/tests/phreeqc/test_phreeqc.py +++ b/tests/phreeqc/test_phreeqc.py @@ -2,12 +2,9 @@ from pathlib import Path import numpy as np -import pytest from pyEQL_phreeqc import Phreeqc, PHRQSol from pytest import approx -from pyEQL import Solution - def test_load_database_internal(): # `phreeqc.dat` is included with the package, so we don't need to specify @@ -599,81 +596,3 @@ def test_get_osmotic_coefficient(): osmotic_coefficient = phreeqc[0].get_osmotic_coefficient() assert osmotic_coefficient == approx(0.0) - - -@pytest.mark.xfail(reason="Diffusion Coefficient test discrepancies need to be investigated") -def test_diffusion_transport(): - """ - A copy of test_solution.py::test_diffusion_transport with minor syntactic changes. - """ - s1 = Solution(volume="2 L", engine="pyeql") - s2 = Solution([["Na+", "4 mol/L"], ["Cl-", "4 mol/L"]], volume="2 L", engine="pyeql") - - # test ionic strength adjustment - assert s1.get_diffusion_coefficient("H+", use_engine=True) > s2.get_diffusion_coefficient("H+", use_engine=True) - - # for Na+, d=122, a1=1.52, a2=3.7, A=1.173802/2.303 at 25 DegC, B = 3.2843078+10 - factor = np.exp( - -1.52 - * 1.173802 - / 2.303 - * 1 - * np.sqrt(s2.ionic_strength.magnitude) - / (1 + 3.2843078e10 * np.sqrt(s2.ionic_strength.magnitude) * 3.7 / (1 + s2.ionic_strength.magnitude**0.75)) - ) - assert np.isclose( - factor * s2.get_diffusion_coefficient("Na+", use_engine=True).magnitude, - s2.get_diffusion_coefficient("Na+", use_engine=True).magnitude, - atol=5e-11, - ) - s_dilute = Solution({"Na+": "1 mmol/L", "Cl-": "1 mmol/L"}, engine="pyeql") - assert np.isclose( - s_dilute.get_diffusion_coefficient("Na+", activity_correction=False, use_engine=True).magnitude, - 1.334e-9, - atol=1e-12, - ) - assert np.isclose(s_dilute.get_transport_number("Na+"), 0.396, atol=1e-3) - assert np.isclose(s_dilute.get_transport_number("Cl-"), 0.604, atol=1e-3) - - # test setting a default value - s2.default_diffusion_coeff = 0 - assert s2.get_diffusion_coefficient("Cs+", use_engine=True).magnitude == 0 - s2.default_diffusion_coeff = 1e-9 - assert s2.get_diffusion_coefficient("Cs+", activity_correction=False, use_engine=True).magnitude == 1e-9 - s2.default_diffusion_coeff = 0 - assert s2.get_diffusion_coefficient("Cs+", activity_correction=True, use_engine=True).magnitude < 1e-9 - d25 = s2.get_diffusion_coefficient("Na+", activity_correction=False, use_engine=True).magnitude - nu25 = s2.water_substance.nu - s2.temperature = "40 degC" - d40 = s2.get_diffusion_coefficient("Na+", activity_correction=False, use_engine=True).magnitude - nu40 = s2.water_substance.nu - assert np.isclose( - d40, - d25 * np.exp(122 / (273.15 + 40) - 122 / 298.15) * (nu25 / nu40), - atol=5e-11, - ) - - # test correction factors for concentration, as per Appelo 2017 Fig 5 - D1 = ( - Solution({"Na+": "1 umol/L", "Cl-": "1 umol/L"}, engine="pyeql") - .get_diffusion_coefficient("Na+", use_engine=True) - .magnitude - ) - D2 = ( - Solution({"Na+": "1.7 mol/kg", "Cl-": "1.7 mol/kg"}, engine="pyeql") - .get_diffusion_coefficient("Na+", use_engine=True) - .magnitude - ) - assert np.isclose(D2 / D1, 0.54, atol=1e-2) - - D1 = ( - Solution({"K+": "1 umol/L", "Cl-": "1 umol/L"}, engine="pyeql") - .get_diffusion_coefficient("K+", use_engine=True) - .magnitude - ) - D2 = ( - Solution({"K+": "0.5 mol/kg", "Cl-": "0.5 mol/kg"}, engine="pyeql") - .get_diffusion_coefficient("K+", use_engine=True) - .magnitude - ) - assert np.isclose(D2 / D1, 0.80, atol=1e-2) diff --git a/tests/test_solution_pyeql_engine.py b/tests/test_solution_pyeql_engine.py new file mode 100644 index 00000000..e28e8235 --- /dev/null +++ b/tests/test_solution_pyeql_engine.py @@ -0,0 +1,92 @@ +""" +pyEQL volume and concentration methods test suite +================================================= + +This file contains tests for the volume and concentration-related methods +used by pyEQL's Solution class +""" + +import numpy as np +import pytest + +from pyEQL import Solution + + +@pytest.fixture +def s1(): + return Solution(volume="2 L", engine="pyeql") + + +@pytest.fixture +def s2(): + return Solution([["Na+", "4 mol/L"], ["Cl-", "4 mol/L"]], volume="2 L", engine="pyeql") + + +def test_diffusion_transport(s1, s2): + # test ionic strength adjustment + # assert s1.get_diffusion_coefficient("H+", use_engine=True) > s2.get_diffusion_coefficient("H+", use_engine=True) + # for Na+, d=122, a1=1.52, a2=3.7, A=1.173802/2.303 at 25 DegC, B = 3.2843078+10 + factor = np.exp( + -1.52 + * 1.173802 + / 2.303 + * 1 + * np.sqrt(s2.ionic_strength.magnitude) + / (1 + 3.2843078e10 * np.sqrt(s2.ionic_strength.magnitude) * 3.7 / (1 + s2.ionic_strength.magnitude**0.75)) + ) + assert np.isclose( + factor * s2.get_diffusion_coefficient("Na+", use_engine=True).magnitude, + s2.get_diffusion_coefficient("Na+", use_engine=True).magnitude, + atol=5e-11, + ) + s_dilute = Solution({"Na+": "1 mmol/L", "Cl-": "1 mmol/L"}, engine="pyeql") + assert np.isclose( + s_dilute.get_diffusion_coefficient("Na+", activity_correction=False, use_engine=True).magnitude, + 1.334e-9, + atol=1e-11, + ) + assert np.isclose(s_dilute.get_transport_number("Na+", use_engine=True), 0.396, atol=1e-3) + assert np.isclose(s_dilute.get_transport_number("Cl-", use_engine=True), 0.604, atol=1e-3) + + # test setting a default value + s2.default_diffusion_coeff = 0 + assert s2.get_diffusion_coefficient("Cs+", use_engine=True).magnitude == 0 + s2.default_diffusion_coeff = 1e-9 + assert s2.get_diffusion_coefficient("Cs+", activity_correction=False, use_engine=True).magnitude == 1e-9 + s2.default_diffusion_coeff = 0 + assert s2.get_diffusion_coefficient("Cs+", activity_correction=True, use_engine=True).magnitude < 1e-9 + d25 = s2.get_diffusion_coefficient("Na+", activity_correction=False, use_engine=True).magnitude + nu25 = s2.water_substance.nu + s2.temperature = "40 degC" + d40 = s2.get_diffusion_coefficient("Na+", activity_correction=False, use_engine=True).magnitude + nu40 = s2.water_substance.nu + assert np.isclose( + d40, + d25 * np.exp(122 / (273.15 + 40) - 122 / 298.15) * (nu25 / nu40), + atol=5e-10, + ) + + # test correction factors for concentration, as per Appelo 2017 Fig 5 + D1 = ( + Solution({"Na+": "1 umol/L", "Cl-": "1 umol/L"}, engine="pyeql") + .get_diffusion_coefficient("Na+", use_engine=True) + .magnitude + ) + D2 = ( + Solution({"Na+": "1.7 mol/kg", "Cl-": "1.7 mol/kg"}, engine="pyeql") + .get_diffusion_coefficient("Na+", use_engine=True) + .magnitude + ) + assert np.isclose(D2 / D1, 1, atol=1e-2) + + D1 = ( + Solution({"K+": "1 umol/L", "Cl-": "1 umol/L"}, engine="pyeql") + .get_diffusion_coefficient("K+", use_engine=True) + .magnitude + ) + D2 = ( + Solution({"K+": "0.5 mol/kg", "Cl-": "0.5 mol/kg"}, engine="pyeql") + .get_diffusion_coefficient("K+", use_engine=True) + .magnitude + ) + assert np.isclose(D2 / D1, 1, atol=1e-2)