Merge pull request #79 from KingsburyLab/bugfix

Miscellanous small enhancements

CHANGELOG.md

@@ -5,16 +5,18 @@ All notable changes to this project will be documented in this file.
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
-## [Unreleased]
-
-### Changed
-
-- `Solution.water_substance` - use the IAPWS97 model instead of IAPWS95 whenever possible, for a substantial speedup.
+## [0.10.1] - 2023-11-12
 
 ### Added
 
 - utility function `create_water_substance` with caching to speed up access to IAPWS instances
 
+### Changed
+
+- `Solution.get_diffusion_coefficient`: the default diffusion coefficient (returned when D for a solute is not found in
+  the database) is now adjusted for temperature and ionic strength.
+- `Solution.water_substance` - use the IAPWS97 model instead of IAPWS95 whenever possible, for a substantial speedup.
+
 ## [0.10.0] - 2023-11-12
 
 ### Added

src/pyEQL/engines.py

@@ -112,18 +112,18 @@ def get_activity_coefficient(self, solution, solute):
         Return the *molal scale* activity coefficient of solute, given a Solution
         object.
         """
-        return ureg.Quantity("1 dimensionless")
+        return ureg.Quantity(1, "dimensionless")
 
     def get_osmotic_coefficient(self, solution):
         """
         Return the *molal scale* osmotic coefficient of solute, given a Solution
         object.
         """
-        return ureg.Quantity("1 dimensionless")
+        return ureg.Quantity(1, "dimensionless")
 
     def get_solute_volume(self, solution):
         """Return the volume of the solutes."""
-        return ureg.Quantity("0 L")
+        return ureg.Quantity(0, "L")
 
     def equilibrate(self, solution):
         """Adjust the speciation of a Solution object to achieve chemical equilibrium."""
@@ -219,7 +219,7 @@ def get_activity_coefficient(self, solution, solute):
         # show an error if no salt can be found that contains the solute
         if salt is None:
             logger.warning("No salts found that contain solute %s. Returning unit activity coefficient." % solute)
-            return ureg.Quantity("1 dimensionless")
+            return ureg.Quantity(1, "dimensionless")
 
         # use the Pitzer model for higher ionic strength, if the parameters are available
         # search for Pitzer parameters
@@ -314,7 +314,7 @@ def get_activity_coefficient(self, solution, solute):
                 % solute
             )
 
-            molal = ureg.Quantity("1 dimensionless")
+            molal = ureg.Quantity(1, "dimensionless")
 
         return molal
 
@@ -473,7 +473,7 @@ def get_solute_volume(self, solution):
         """Return the volume of the solutes."""
         # identify the predominant salt in the solution
         salt = solution.get_salt()
-        solute_vol = ureg.Quantity("0 L")
+        solute_vol = ureg.Quantity(0, "L")
 
         # use the pitzer approach if parameters are available
         pitzer_calc = False
@@ -695,12 +695,12 @@ def get_osmotic_coefficient(self, solution):
         via phreeqcpython
         """
         # TODO - find a way to access or calculate osmotic coefficient
-        return ureg.Quantity("1 dimensionless")
+        return ureg.Quantity(1, "dimensionless")
 
     def get_solute_volume(self, solution):
         """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")
+        return ureg.Quantity(0, "L")
 
     def equilibrate(self, solution):
         """Adjust the speciation of a Solution object to achieve chemical equilibrium."""

src/pyEQL/solution.py

@@ -127,7 +127,7 @@ def __init__(
             self._volume = ureg.Quantity(volume).to("L")
         else:
             # volume_set = False
-            self._volume = ureg.Quantity("1 L")
+            self._volume = ureg.Quantity(1, "L")
         # store the initial conditions as private variables in case they are
         # changed later
         self._temperature = ureg.Quantity(temperature)
@@ -193,8 +193,8 @@ def __init__(
         # self.add_solvent(self.solvent, kwargs["solvent"][1])
 
         # calculate the moles of solvent (water) on the density and the solution volume
-        moles = self.volume / ureg.Quantity("55.55 mol/L")
-        self.components["H2O"] = moles.magnitude
+        moles = self.volume.magnitude / 55.55  # molarity of pure water
+        self.components["H2O"] = moles
 
         # set the pH with H+ and OH-
         self.add_solute("H+", str(10 ** (-1 * pH)) + "mol/L")
@@ -759,7 +759,7 @@ def alkalinity(self) -> Quantity:
         .. [stm] Stumm, Werner and Morgan, James J. Aquatic Chemistry, 3rd ed, pp 165. Wiley Interscience, 1996.
 
         """
-        alkalinity = ureg.Quantity("0 mol/L")
+        alkalinity = ureg.Quantity(0, "mol/L")
 
         base_cations = {
             "Li[+1]",
@@ -806,7 +806,7 @@ def hardness(self) -> Quantity:
             The hardness of the solution in mg/L as CaCO3
 
         """
-        hardness = ureg.Quantity("0 mol/L")
+        hardness = ureg.Quantity(0, "mol/L")
 
         for item in self.components:
             z = self.get_property(item, "charge")
@@ -823,7 +823,7 @@ def total_dissolved_solids(self) -> Quantity:
 
         The TDS is defined as the sum of the concentrations of all aqueous solutes (not including the solvent), except for H[+1] and OH[-1]].
         """
-        tds = ureg.Quantity("0 mg/L")
+        tds = ureg.Quantity(0, "mg/L")
         for s in self.components:
             # ignore pure water and dissolved gases, but not CO2
             if s in ["H2O(aq)", "H[+1]", "OH[-1]"]:
@@ -1252,7 +1252,7 @@ def add_solute(self, formula: str, amount: str):
 
             # update the volume to account for the space occupied by all the solutes
             # make sure that there is still solvent present in the first place
-            if self.solvent_mass <= ureg.Quantity("0 kg"):
+            if self.solvent_mass <= ureg.Quantity(0, "kg"):
                 logger.error("All solvent has been depleted from the solution")
                 return
             # set the volume recalculation flag
@@ -1356,7 +1356,7 @@ def add_amount(self, solute: str, amount: str):
 
             # update the volume to account for the space occupied by all the solutes
             # make sure that there is still solvent present in the first place
-            if self.solvent_mass <= ureg.Quantity("0 kg"):
+            if self.solvent_mass <= ureg.Quantity(0, "kg"):
                 logger.error("All solvent has been depleted from the solution")
                 return
 
@@ -1444,7 +1444,7 @@ def set_amount(self, solute: str, amount: str):
 
             # update the volume to account for the space occupied by all the solutes
             # make sure that there is still solvent present in the first place
-            if self.solvent_mass <= ureg.Quantity("0 kg"):
+            if self.solvent_mass <= ureg.Quantity(0, "kg"):
                 logger.error("All solvent has been depleted from the solution")
                 return
 
@@ -1753,7 +1753,7 @@ def get_activity_coefficient(
         """
         # return unit activity coefficient if the concentration of the solute is zero
         if self.get_amount(solute, "mol").magnitude == 0:
-            return ureg.Quantity("1 dimensionless")
+            return ureg.Quantity(1, "dimensionless")
 
         try:
             # get the molal-scale activity coefficient from the EOS engine
@@ -1768,11 +1768,11 @@ def get_activity_coefficient(
         if scale == "molar":
             total_molality = self.get_total_moles_solute() / self.solvent_mass
             total_molarity = self.get_total_moles_solute() / self.volume
-            return (molal * self.water_substance.rho * ureg.Quantity("1 g/L") * total_molality / total_molarity).to(
+            return (molal * ureg.Quantity(self.water_substance.rho, "g/L") * total_molality / total_molarity).to(
                 "dimensionless"
             )
         if scale == "rational":
-            return molal * (1 + ureg.Quantity("0.018015 kg/mol") * self.get_total_moles_solute() / self.solvent_mass)
+            return molal * (1 + ureg.Quantity(0.018015, "kg/mol") * self.get_total_moles_solute() / self.solvent_mass)
 
         raise ValueError("Invalid scale argument. Pass 'molal', 'molar', or 'rational'.")
 
@@ -1953,7 +1953,7 @@ def get_chemical_potential_energy(self, activity_correction: bool = True) -> Qua
         .. [koga] Koga, Yoshikata, 2007. *Solution Thermodynamics and its Application to Aqueous Solutions: A differential approach.* Elsevier, 2007, pp. 23-37.
 
         """
-        E = ureg.Quantity("0 J")
+        E = ureg.Quantity(0, "J")
 
         # loop through all the components and add their potential energy
         for item in self.components:
@@ -2254,10 +2254,8 @@ def _get_diffusion_coefficient(self, solute: str, activity_correction: bool = Tr
         D = self.get_property(solute, "transport.diffusion_coefficient")
         rform = standardize_formula(solute)
         if D is None or D.magnitude == 0:
-            logger.info(
-                f"Diffusion coefficient not found for species {rform}. Return default value of {default} m**2/s."
-            )
-            return ureg.Quantity(default, "m**2/s")
+            logger.info(f"Diffusion coefficient not found for species {rform}. Use default value of {default} m**2/s.")
+            D = ureg.Quantity(default, "m**2/s")
 
         # assume reference temperature is 298.15 K (this is the case for all current DB entries)
         T_ref = 298.15

tests/test_solute_properties.py

@@ -130,15 +130,15 @@ def test_molar_conductivity_hydrogen(self):
     def test_molar_conductivity_neutral(self):
         s1 = Solution([["FeCl3", "0.001 mol/L"]], temperature="25 degC")
         result = s1.get_molar_conductivity("FeCl3").to("m**2*S/mol").magnitude
-        expected = ureg.Quantity("0 m**2 * S / mol").magnitude
+        expected = ureg.Quantity(0, "m**2 * S / mol").magnitude
 
         assert round(abs(result - expected), 5) == 0
 
     # molar conductivity of water should be zero
     def test_molar_conductivity_water(self):
         s1 = Solution(temperature="25 degC")
         result = s1.get_molar_conductivity("H2O").to("m**2*S/mol").magnitude
-        expected = ureg.Quantity("0 m**2 * S / mol").magnitude
+        expected = ureg.Quantity(0, "m**2 * S / mol").magnitude
 
         assert round(abs(result - expected), 5) == 0
 

tests/test_solution.py

@@ -139,7 +139,8 @@ def test_diffusion_transport(s1, s2):
 
     # test setting a default value
     assert s2.get_diffusion_coefficient("Cs+").magnitude == 0
-    assert s2.get_diffusion_coefficient("Cs+", default=1e-9).magnitude == 1e-9
+    assert s2.get_diffusion_coefficient("Cs+", default=1e-9, activity_correction=False).magnitude == 1e-9
+    assert s2.get_diffusion_coefficient("Cs+", default=1e-9, activity_correction=True).magnitude < 1e-9
     d25 = s2.get_diffusion_coefficient("Na+", activity_correction=False).magnitude
     nu25 = s2.water_substance.nu
     s2.temperature = "40 degC"

tox.ini

@@ -3,8 +3,8 @@
 # THIS SCRIPT IS SUPPOSED TO BE AN EXAMPLE. MODIFY IT ACCORDING TO YOUR NEEDS!
 
 [tox]
-minversion = 3.15
-envlist = {py38,py39,py310,py11}-{linux,windows,macos}
+minversion = 4
+envlist = {py39,py310,py11,py312}-{linux,windows,macos}
 
 
 [testenv]