Fix charge balancing error and dependency housekeeping

CHANGELOG.md

@@ -5,14 +5,21 @@ 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
+## [1.0.1] - 2024-06-17
+
+### Added
+
+- Added `NPY201` ruleset to `ruff` configuration to check support for `numpy` 2.0
 
 ### Fixed
 
+- `NativeEOS` and `PHREEQCEOS` `equilibrate`: Fixed a charge balancing bug that cause repeated calls to `equlibrate` to
+  severely distort the pH and/or the composition. (#141)[https://github.com/KingsburyLab/pyEQL/issues/141]
 - `UnicodeDecodeError` when trying to connect to ion database on non-English platforms ([#122](https://github.com/KingsburyLab/pyEQL/issues/122))
 
 ### Changed
 
+- Replace `math` with `numpy` functions throughout. This mainly changed calls to `log` and `exp`.
 - Docs: Fix missing close parentheses in docstring (@Andrew S. Rosen)
 
 ## [1.0.0] - 2024-03-17

pyproject.toml

@@ -47,6 +47,7 @@ select = [
   "UP",   # pyupgrade
   "W",    # pycodestyle warning
   "YTT",  # flake8-2020
+  "NPY201", # Numpy 2.0 migration
 ]
 ignore = [
   "B008",    # Do not perform function call `unit` in argument defaults

setup.cfg

@@ -50,9 +50,9 @@ python_requires = >=3.9
 # For more information, check out https://semver.org/.
 install_requires =
     pint>=0.19
-    numpy
+    numpy<2
     scipy
-    pymatgen>=2023.10.11
+    pymatgen==2024.5.1
     iapws
     monty
     maggma>=0.67.0

src/pyEQL/activity_correction.py

@@ -14,8 +14,8 @@
 """
 
 import logging
-import math
 
+import numpy as np
 from pint import Quantity
 
 from pyEQL import ureg
@@ -114,8 +114,8 @@ def _debye_parameter_activity(temperature: str = "25 degC") -> "Quantity":
 
     debyeparam = (
         ureg.elementary_charge**3
-        * (2 * math.pi * ureg.N_A * ureg.Quantity(water_substance.rho, "g/L")) ** 0.5
-        / (4 * math.pi * ureg.epsilon_0 * water_substance.epsilon * ureg.boltzmann_constant * T) ** 1.5
+        * (2 * np.pi * ureg.N_A * ureg.Quantity(water_substance.rho, "g/L")) ** 0.5
+        / (4 * np.pi * ureg.epsilon_0 * water_substance.epsilon * ureg.boltzmann_constant * T) ** 1.5
     )
 
     logger.debug(rf"Computed Debye-Huckel Limiting Law Constant A^{{\gamma}} = {debyeparam} at {temperature}")
@@ -248,7 +248,7 @@ def get_activity_coefficient_debyehuckel(ionic_strength, z=1, temperature="25 de
 
     log_f = -_debye_parameter_activity(temperature) * z**2 * ionic_strength**0.5
 
-    return math.exp(log_f) * ureg.Quantity(1, "dimensionless")
+    return np.exp(log_f) * ureg.Quantity(1, "dimensionless")
 
 
 def get_activity_coefficient_guntelberg(ionic_strength, z=1, temperature="25 degC"):
@@ -285,7 +285,7 @@ def get_activity_coefficient_guntelberg(ionic_strength, z=1, temperature="25 deg
 
     log_f = -_debye_parameter_activity(temperature) * z**2 * ionic_strength**0.5 / (1 + ionic_strength.magnitude**0.5)
 
-    return math.exp(log_f) * ureg.Quantity(1, "dimensionless")
+    return np.exp(log_f) * ureg.Quantity(1, "dimensionless")
 
 
 def get_activity_coefficient_davies(ionic_strength, z=1, temperature="25 degC"):
@@ -327,7 +327,7 @@ def get_activity_coefficient_davies(ionic_strength, z=1, temperature="25 degC"):
         * (ionic_strength.magnitude**0.5 / (1 + ionic_strength.magnitude**0.5) - 0.2 * ionic_strength.magnitude)
     )
 
-    return math.exp(log_f) * ureg.Quantity(1, "dimensionless")
+    return np.exp(log_f) * ureg.Quantity(1, "dimensionless")
 
 
 def get_activity_coefficient_pitzer(
@@ -437,7 +437,7 @@ def get_activity_coefficient_pitzer(
         b,
     )
 
-    return math.exp(loggamma) * ureg.Quantity(1, "dimensionless")
+    return np.exp(loggamma) * ureg.Quantity(1, "dimensionless")
 
 
 def get_apparent_volume_pitzer(
@@ -533,7 +533,7 @@ def get_apparent_volume_pitzer(
         (nu_cation + nu_anion)
         * abs(z_cation * z_anion)
         * (_debye_parameter_volume(temperature) / 2 / b)
-        * math.log(1 + b * ionic_strength**0.5)
+        * np.log(1 + b * ionic_strength**0.5)
     )
 
     third_term = (
@@ -566,7 +566,7 @@ def _pitzer_f1(x):
     # return 0 if the input is 0
     if x == 0:
         return 0
-    return 2 * (1 - (1 + x) * math.exp(-x)) / x**2
+    return 2 * (1 - (1 + x) * np.exp(-x)) / x**2
 
 
 def _pitzer_f2(x):
@@ -586,7 +586,7 @@ def _pitzer_f2(x):
     # return 0 if the input is 0
     if x == 0:
         return 0
-    return -2 * (1 - (1 + x + x**2 / 2) * math.exp(-x)) / x**2
+    return -2 * (1 - (1 + x + x**2 / 2) * np.exp(-x)) / x**2
 
 
 def _pitzer_B_MX(ionic_strength, alpha1, alpha2, beta0, beta1, beta2):
@@ -692,6 +692,10 @@ def _pitzer_B_phi(ionic_strength, alpha1, alpha2, beta0, beta1, beta2):
         and Representation with an Ion Interaction (Pitzer) Model.
         Journal of Chemical & Engineering Data, 55(2), 830-838. doi:10.1021/je900487a
     """
+    import math
+
+    # TODO - for some reason this specific method requires the use of math.exp rather than np.exp. Using np.exp raises
+    # a dimensionalityerror.
     return beta0 + beta1 * math.exp(-alpha1 * ionic_strength**0.5) + beta2 * math.exp(-alpha2 * ionic_strength**0.5)
 
 
@@ -772,7 +776,7 @@ def _pitzer_log_gamma(
         -1
         * abs(z_cation * z_anion)
         * _debye_parameter_osmotic(temperature)
-        * (ionic_strength**0.5 / (1 + b * ionic_strength**0.5) + 2 / b * math.log(1 + b * ionic_strength**0.5))
+        * (ionic_strength**0.5 / (1 + b * ionic_strength**0.5) + 2 / b * np.log(1 + b * ionic_strength**0.5))
     )
     second_term = 2 * molality * nu_cation * nu_anion / (nu_cation + nu_anion) * (B_MX + B_phi)
     third_term = 3 * molality**2 * (nu_cation * nu_anion) ** 1.5 / (nu_cation + nu_anion) * C_phi

src/pyEQL/engines.py

@@ -47,7 +47,7 @@ def get_activity_coefficient(self, solution, solute):
             solution: pyEQL Solution object
             solute: str identifying the solute of interest
 
-        Returns
+        Returns:
             Quantity: dimensionless quantity object
 
         Raises:
@@ -62,7 +62,7 @@ def get_osmotic_coefficient(self, solution):
         Args:
             solution: pyEQL Solution object
 
-        Returns
+        Returns:
             Quantity: dimensionless molal scale osmotic coefficient
 
         Raises:
@@ -77,7 +77,7 @@ def get_solute_volume(self):
         Args:
             solution: pyEQL Solution object
 
-        Returns
+        Returns:
             Quantity: solute volume in L
 
         Raises:
@@ -94,7 +94,7 @@ def equilibrate(self, solution):
         Args:
             solution: pyEQL Solution object
 
-        Returns
+        Returns:
             Nothing. The speciation of the Solution is modified in-place.
 
         Raises:
@@ -173,9 +173,7 @@ def __init__(
         self._stored_comp = None
 
     def _setup_ppsol(self, solution):
-        """
-        Helper method to set up a PhreeqPython solution for subsequent analysis.
-        """
+        """Helper method to set up a PhreeqPython solution for subsequent analysis."""
         self._stored_comp = solution.components.copy()
         solv_mass = solution.solvent_mass.to("kg").magnitude
         # inherit bulk solution properties
@@ -197,6 +195,16 @@ def _setup_ppsol(self, solution):
         # add the composition to the dict
         # also, skip H and O
         for el, mol in solution.get_el_amt_dict().items():
+            # CAUTION - care must be taken to avoid unintended behavior here. get_el_amt_dict() will return
+            # all distinct oxi states of each element present. If there are elements present whose oxi states
+            # are NOT recognized by PHREEQC (via SPECIAL_ELEMENTS) then the amount of only 1 oxi state will be
+            # entered into the composition dict. This can especially cause problems after equilibrate() has already
+            # been called once. For example, equilibrating a simple NaCl solution generates Cl species that are assigned
+            # various oxidations states, -1 mostly, but also 1, 2, and 3. Since the concentrations of everything
+            # except the -1 oxi state are tiny, this can result in Cl "disappearing" from the solution if
+            # equlibrate is called again. It also causes non-determinism, because the amount is taken from whatever
+            # oxi state happens to be iterated through last.
+
             # strip off the oxi state
             bare_el = el.split("(")[0]
             if bare_el in SPECIAL_ELEMENTS:
@@ -208,7 +216,12 @@ def _setup_ppsol(self, solution):
             else:
                 key = bare_el
 
-            d[key] = str(mol / solv_mass)
+            if key in d:
+                # when multiple oxi states for the same (non-SPECIAL) element are present, make sure to
+                # add all their amounts together
+                d[key] += str(mol / solv_mass)
+            else:
+                d[key] = str(mol / solv_mass)
 
             # tell PHREEQC which species to use for charge balance
             if (
@@ -660,7 +673,6 @@ def equilibrate(self, solution):
             solution.components[s] = mol
 
         # make sure all species are accounted for
-        charge_adjust = 0
         assert set(self._stored_comp.keys()) - set(solution.components.keys()) == set()
 
         # log a message if any components were not touched by PHREEQC
@@ -671,6 +683,11 @@ def equilibrate(self, solution):
                 f"After equilibration, the amounts of species {missing_species} were not modified "
                 "by PHREEQC. These species are likely absent from its database."
             )
+
+        # re-adjust charge balance for any missing species
+        # note that if balance_charge is set, it will have been passed to PHREEQC, so we only need to adjust
+        # for any missing species here.
+        charge_adjust = 0
         for s in missing_species:
             charge_adjust += -1 * solution.get_amount(s, "eq").magnitude
         if charge_adjust != 0:
@@ -679,11 +696,10 @@ def equilibrate(self, solution):
                 f" {charge_adjust} eq of charge were added via {solution.balance_charge}"
             )
 
-            # re-adjust charge balance
             if solution.balance_charge is None:
                 pass
             elif solution.balance_charge == "pH":
-                solution.components["H+"] += charge_adjust
+                solution.components["H+"] += charge_adjust.magnitude
             elif solution.balance_charge == "pE":
                 raise NotImplementedError
             else:

src/pyEQL/equilibrium.py

@@ -11,7 +11,8 @@
 
 # import libraries for scientific functions
 import logging
-import math
+
+import numpy as np
 
 from pyEQL import ureg
 
@@ -51,7 +52,7 @@ def adjust_temp_pitzer(c1, c2, c3, c4, c5, temp, temp_ref=ureg.Quantity("298.15
     return (
         c1
         + c2 * (1 / temp + 1 / temp_ref)
-        + c2 * math.log(temp / temp_ref)
+        + c2 * np.log(temp / temp_ref)
         + c3 * (temp - temp_ref)
         + c4 * (temp**2 - temp_ref**2)
         + c5 * (temp**-2 - temp_ref**-2)
@@ -91,7 +92,7 @@ def adjust_temp_vanthoff(equilibrium_constant, enthalpy, temperature, reference_
         >>> adjust_temp_vanthoff(0.15,ureg.Quantity('-197.6 kJ/mol'),ureg.Quantity('42 degC')) #doctest: +ELLIPSIS
         0.00203566...
     """
-    output = equilibrium_constant * math.exp(
+    output = equilibrium_constant * np.exp(
         enthalpy / ureg.R * (1 / reference_temperature.to("K") - 1 / temperature.to("K"))
     )
 
@@ -137,7 +138,7 @@ def adjust_temp_arrhenius(
         >>> adjust_temp_arrhenius(7,900*ureg.Quantity('kJ/mol'),37*ureg.Quantity('degC'),97*ureg.Quantity('degC')) #doctest: +ELLIPSIS
         1.8867225...e-24
     """
-    output = rate_constant * math.exp(
+    output = rate_constant * np.exp(
         activation_energy / ureg.R * (1 / reference_temperature.to("K") - 1 / temperature.to("K"))
     )
 

src/pyEQL/functions.py

@@ -7,7 +7,8 @@
 """
 
 import logging
-import math
+
+import numpy as np
 
 from pyEQL import Solution, ureg
 
@@ -55,7 +56,7 @@ def gibbs_mix(solution1: Solution, solution2: Solution):
     for solution in term_list:
         for solute in solution.components:
             if solution.get_amount(solute, "fraction") != 0:
-                term_list[solution] += solution.get_amount(solute, "mol") * math.log(solution.get_activity(solute))
+                term_list[solution] += solution.get_amount(solute, "mol") * np.log(solution.get_activity(solute))
 
     return (ureg.R * blend.temperature.to("K") * (term_list[blend] - term_list[concentrate] - term_list[dilute])).to(
         "J"
@@ -102,7 +103,7 @@ def entropy_mix(solution1: Solution, solution2: Solution):
     for solution in term_list:
         for solute in solution.components:
             if solution.get_amount(solute, "fraction") != 0:
-                term_list[solution] += solution.get_amount(solute, "mol") * math.log(
+                term_list[solution] += solution.get_amount(solute, "mol") * np.log(
                     solution.get_amount(solute, "fraction")
                 )
 
@@ -242,7 +243,7 @@ def donnan_solve(x):
 
         return (act_cation_mem / act_cation_soln) ** (1 / z_cation) * (act_anion_soln / act_anion_mem) ** (
             1 / z_anion
-        ) - math.exp(delta_pi * exp_term)
+        ) - np.exp(delta_pi * exp_term)
 
     # solve the function above using one of scipy's nonlinear solvers