diff --git a/src/pyEQL/engines.py b/src/pyEQL/engines.py index 34a4b162..461604fd 100644 --- a/src/pyEQL/engines.py +++ b/src/pyEQL/engines.py @@ -1128,3 +1128,23 @@ 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) + 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/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..b42c4d43 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(eq_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: @@ -31,6 +33,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..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 = ( @@ -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,8 +2222,20 @@ def _get_diffusion_coefficient(self, solute: str, activity_correction: bool = Tr pyEQL.activity_correction._debye_parameter_activity """ - D = self.get_property(solute, "transport.diffusion_coefficient") rform = standardize_formula(solute) + + if use_engine: + assert self._engine == "pyeql", "Only supported for pyeql engine" + 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") 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 df2da2f2..2fc716f0 100644 --- a/tests/phreeqc/test_phreeqc.py +++ b/tests/phreeqc/test_phreeqc.py @@ -406,7 +406,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 +454,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(): 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)