"""
pyEQL Solution Class.
:copyright: 2013-2024 by Ryan S. Kingsbury
:license: LGPL, see LICENSE for more details.
"""
from __future__ import annotations
import logging
import os
import warnings
from functools import lru_cache
from importlib.resources import files
from pathlib import Path
from typing import Any, Literal
import numpy as np
from maggma.stores import JSONStore, Store
from monty.dev import deprecated
from monty.json import MontyDecoder, MSONable
from monty.serialization import dumpfn, loadfn
from pint import DimensionalityError, Quantity
from pymatgen.core import Element
from pymatgen.core.ion import Ion
from pyEQL import IonDB, ureg
from pyEQL.activity_correction import _debye_parameter_activity, _debye_parameter_B
from pyEQL.engines import EOS, IdealEOS, NativeEOS, Phreeqc2026EOS, PhreeqcEOS
from pyEQL.salt_ion_match import Salt
from pyEQL.solute import Solute
from pyEQL.utils import (
FormulaDict,
_translate_pint_quantity,
create_water_substance,
standardize_formula,
translate_units,
)
EQUIV_WT_CACO3 = ureg.Quantity(100.09 / 2, "g/mol")
# string to denote unknown oxidation states
UNKNOWN_OXI_STATE = "unk"
K_W = 1e-14 # ion product of water at 25 degC
# the following properties will be pre-cached for all solutes in the solution, to speed up property retrieval during calculations.
CORE_PROPERTIES = [
"formula",
"name",
"charge",
"molecular_weight",
"size.molar_volume",
"oxi_state_guesses",
"elements",
"pmg_ion",
"model_parameters",
"transport.diffusion_coefficient",
]
[docs]
class Solution(MSONable):
"""
Class representing the properties of a solution. Instances of this class
contain information about the solutes, solvent, and bulk properties.
"""
def __init__(
self,
solutes: list[list[str]] | dict[str, str] | None = None,
volume: str | None = None,
temperature: str = "298.15 K",
pressure: str = "1 atm",
pH: float = 7,
pE: float = 8.5,
balance_charge: str | None = None,
solvent: str | list = "H2O",
engine: EOS | Literal["native", "ideal", "phreeqc", "phreeqc2026"] = "native",
database: str | Path | Store | None = None,
default_diffusion_coeff: float = 1.6106e-9,
log_level: Literal["DEBUG", "INFO", "WARNING", "ERROR", "CRITICAL"] | None = "ERROR",
) -> None:
r"""
Instantiate a Solution from a composition.
Args:
solutes: dict, optional. Keys must be the chemical formula, while values must be
str Quantity representing the amount. For example:
{"Na+": "0.1 mol/L", "Cl-": "0.1 mol/L"}
Note that an older "list of lists" syntax is also supported; however this
will be deprecated in the future and is no longer recommended. The equivalent
list syntax for the above example is
[["Na+", "0.1 mol/L"], ["Cl-", "0.1 mol/L"]]
Defaults to empty (pure solvent) if omitted
volume: str, optional
Volume of the solvent, including the unit. Defaults to '1 L' if omitted.
Note that the total solution volume will be computed using partial molar
volumes of the respective solutes as they are added to the solution.
temperature: str, optional
The solution temperature, including the ureg. Defaults to '25 degC' if omitted.
pressure: Quantity, optional
The ambient pressure of the solution, including the unit.
Defaults to '1 atm' if omitted.
pH: number, optional
Negative log of H+ activity. If omitted, the solution will be
initialized to pH 7 (neutral) with appropriate quantities of
H+ and OH- ions
pE: the :math:`pe` value of the solution. :math:`pe` measures the relative abundance of electrons
analogous to how pH measures the relative abundance of protons. Specifically, :math:`pe` is defined in
terms of the activity of electrons :math:`[e^{-}]`:
.. math:: pe = - \log [e^{-}]
The relationship between the redox potential :math:`Eh` and :math:`pe` can be illustrated by considering
the general redox reaction,
.. math::
\begin{gather*}
\text{A}^x \pm ne^{-} \longrightarrow \text{A}^{x \mp n} \quad\quad
K = \frac{[\text{A}^{x \mp n}]}{[\text{A}^x][e^{-}]^{\pm n}}
\end{gather*}
Writing :math:`pe` in terms of the equilibrium constant :math:`K` and the activities,
:math:`[\text{A}^{x}]` and :math:`[\text{A}^{x \mp n}]`, we have:
.. math::
\begin{gather*}
pe = -\log[e^{-}] = \mp \frac{1}{n} \log\left(\frac{1}{K} \frac{[\text{A}^{x \mp n}]}{[\text{A}^x]}\right)
= \mp \frac{\Delta G}{nRT \ln 10} = \frac{FEh}{RT \ln 10}
\end{gather*}
Thus, the redox potential :math:`Eh` is then related to :math:`pe` via:
.. math:: Eh = 2.303 \frac{RT}{F}pe
where :math:`F` is Faraday's constant. Note that lower values of ``pE`` (and thus :math:`Eh`)
correspond to more reducing environments, while higher values = more oxidizing. At pH 7, water is stable
between approximately -7 to +14. The default value corresponds to a :math:`pe` value typical of natural
waters in equilibrium with the atmosphere.
balance_charge: The strategy for balancing charge during init and equilibrium calculations. Valid options
are
- 'pH', which will adjust the solution pH to balance charge,
- 'auto' which will use the majority cation or anion (i.e., that with the largest concentration)
as needed,
- 'pE' (not currently implemented) which will adjust the redox equilibrium to balance charge, or
the name of a dissolved species e.g. 'Ca+2' or 'Cl-' that will be added/subtracted to balance
charge.
- None (default), in which case no charge balancing will be performed either on init or when
equilibrate() is called. Note that in this case, equilibrate() can distort the charge balance!
solvent: Formula of the solvent. Solvents other than water are not supported at this time.
engine: Electrolyte modeling engine to use. See documentation for details on the available engines.
database: path to a .json file (str or Path) or maggma Store instance that
contains serialized SoluteDocs. `None` (default) will use the built-in pyEQL database.
log_level: Log messages of this or higher severity will be printed to stdout. Defaults to 'ERROR', meaning
that ERROR and CRITICAL messages will be shown, while WARNING, INFO, and DEBUG messages are not. If set
to None, nothing will be printed.
default_diffusion_coeff: Diffusion coefficient value in m^2/s to use in
calculations when there is no diffusion coefficient for a species in the database. This affects several
important property calculations including conductivity and transference number, which are related to the
weighted sums of diffusion coefficients of all species. Setting this argument to zero will exclude any
species that does not have a tabulated diffusion coefficient from these calculations, possibly resulting
in underestimation of the conductivity and/or inaccurate transference numbers.
Missing diffusion coefficients are especially likely in complex electrolytes containing, for example,
complexes or paired species such as NaSO4[-1]. In such cases, setting default_diffusion_coeff to zero
is likely to result in the above errors.
By default, this argument is set to the diffusion coefficient of NaCl salt, 1.61x10^-9 m2/s.
Examples:
>>> s1 = pyEQL.Solution({'Na+': '1 mol/L','Cl-': '1 mol/L'},temperature='20 degC',volume='500 mL')
>>> print(s1)
Components:
Volume: 0.500 l
Pressure: 1.000 atm
Temperature: 293.150 K
Components: ['H2O(aq)', 'H[+1]', 'OH[-1]', 'Na[+1]', 'Cl[-1]']
"""
# create a logger and attach it to this class
self.log_level = log_level.upper()
self.logger = logging.getLogger("pyEQL")
if self.log_level is not None:
# set the level of the module logger
self.logger.setLevel(self.log_level)
# clear handlers and add a StreamHandler
self.logger.handlers.clear()
# use rich for pretty log formatting, if installed
try:
from rich.logging import RichHandler # noqa: PLC0415
sh = RichHandler(rich_tracebacks=True)
except ImportError:
sh = logging.StreamHandler()
# the formatter determines what our logs will look like
formatter = logging.Formatter("[%(asctime)s] [%(levelname)8s] --- %(message)s (%(filename)s:%(lineno)d)")
sh.setFormatter(formatter)
self.logger.addHandler(sh)
# per-instance cache of get_property and other calls that do not depend
# on composition
# see https://rednafi.com/python/lru_cache_on_methods/
self.get_property = lru_cache(maxsize=8192)(self._get_property)
self.get_molar_conductivity = lru_cache()(self._get_molar_conductivity)
self.get_mobility = lru_cache()(self._get_mobility)
self.default_diffusion_coeff = default_diffusion_coeff
self.get_diffusion_coefficient = lru_cache()(self._get_diffusion_coefficient)
# initialize the volume recalculation flag
self.volume_update_required = False
# initialize the volume
if volume is not None:
self._volume = ureg.Quantity(volume).to("L")
else:
self._volume = 1 * ureg.L
# store the initial conditions as private variables in case they are
# changed later
self._temperature = ureg.Quantity(temperature)
self._pressure = ureg.Quantity(pressure)
self._pE = pE
self._pH = pH
self.pE = self._pE
if isinstance(balance_charge, str) and balance_charge not in ["pH", "pE", "auto"]:
self.balance_charge = standardize_formula(balance_charge)
else:
self.balance_charge = balance_charge #: Standardized formula of the species used for charge balancing.
# instantiate a water substance for property retrieval
self.water_substance = create_water_substance(self.temperature, self.pressure)
"""IAPWS instance describing water properties."""
# create an empty dictionary of components. This dict comprises {formula: moles}
# where moles is the number of moles in the solution.
self.components = FormulaDict({})
"""Special dictionary where keys are standardized formula and values are the moles present in Solution."""
# connect to the desired property database
if database is None:
# load the default database, which is a JSONStore
db_store = IonDB
elif isinstance(database, str | Path):
db_store = JSONStore(str(database), key="formula")
self.logger.debug(f"Created maggma JSONStore from .json file {database}")
else:
db_store = database
self.database = db_store
"""`Store` instance containing the solute property database."""
self.database.connect()
self.logger.debug(f"Connected to property database {self.database!s}")
# set the equation of state engine
self._engine = engine
# self.engine: Optional[EOS] = None
if isinstance(self._engine, EOS):
self.engine: EOS = self._engine
elif self._engine == "ideal":
self.engine = IdealEOS()
elif self._engine == "native":
self.engine = NativeEOS()
elif self._engine == "phreeqc":
self.engine = PhreeqcEOS()
elif self._engine == "phreeqc2026":
self.engine = Phreeqc2026EOS()
else:
raise ValueError(f'{engine} is not a valid value for the "engine" kwarg!')
# define the solvent. Allow for list input to support future use of mixed solvents
if not isinstance(solvent, list):
solvent = [solvent]
if len(solvent) > 1:
raise ValueError("Multiple solvents are not yet supported!")
if solvent[0] not in ["H2O", "H2O(aq)", "water", "Water", "HOH"]:
raise ValueError("Non-aqueous solvent detected. These are not yet supported!")
self.solvent = standardize_formula(solvent[0])
"""Formula of the component that is set as the solvent (currently only H2O(aq) is supported)."""
# store the provided solutes as a dict
if isinstance(solutes, dict):
self._solutes = solutes
elif isinstance(solutes, list):
msg = (
'List input of solutes (e.g., [["Na+", "0.5 mol/L]]) is deprecated! Use dictionary formatted input '
'(e.g., {"Na+":"0.5 mol/L"} instead.)'
)
self.logger.warning(msg)
warnings.warn(msg, DeprecationWarning)
self._solutes = {item[0]: item[1] for item in solutes}
elif solutes is None:
self._solutes = {}
# # pre-cache all the relevant properties of the solutes provided by the user
self._solute_data = list(
self.database.query({"formula": {"$in": list(self._solutes.keys())}}, properties=CORE_PROPERTIES)
)
if "H2O(aq)" in self._solutes:
self.components["H2O"] = ureg.Quantity(self._solutes["H2O(aq)"]).to("mol").magnitude
else:
# calculate the moles of solvent (water) based on the density and solution volume
self.components["H2O"] = (
self.volume.magnitude * 1000 * self.water_substance.rho / 18.01528
) # moles = density / molar mass * volume
# set the pH with H+ and OH-
self.add_solute("H+", str(10 ** (-1 * pH)) + "mol/L")
self.add_solute("OH-", str(K_W / (10 ** (-1 * pH))) + "mol/L")
# populate remaining solutes
CHECK_H = False
for k, v in self._solutes.items():
self.add_solute(k, v)
# if user has specified H+ in solutes, check consistency with pH kwarg
if standardize_formula(k) == "H[+1]":
CHECK_H = True
if CHECK_H:
# if user has not specified pH (default value), override the pH argument
if self._pH == 7:
self.logger.warning(f"H[+1] = {v} found in solutes. Overriding default pH with this value.")
# if user specifies non-default pH that does not match the supplied H+, raise an error
elif not np.isclose(self.pH, self._pH, atol=1e-4):
warnings.warn(
f"After initialization, the calculated solution pH of {self.pH:.3f} does not match the "
f"specified pH of {self._pH:.3f}. This might be a result of erroneous input (e.g., specifying "
"both pH and H+), or it can happen during from_dict / from_preset if you use a different "
"engine than the one which generated the original dict."
)
# determine the species that will be used for charge balancing, when needed.
# this is necessary to do even if the composition is already electroneutral,
# because the appropriate species also needs to be passed to equilibrate
# to keep from distorting the charge balance.
cb = self.charge_balance
if self.balance_charge is None:
self._cb_species = None
elif self.balance_charge == "pH":
self._cb_species = "H[+1]"
elif self.balance_charge == "pE":
raise NotImplementedError("Balancing charge via redox (pE) is not yet implemented!")
elif self.balance_charge == "auto":
# add the most abundant ion of the opposite charge
if cb <= 0:
self._cb_species = max(self.cations, key=self.cations.get)
elif cb > 0:
self._cb_species = max(self.anions, key=self.anions.get)
else:
ions = set().union(*[self.cations, self.anions]) # all ions
self._cb_species = self.balance_charge
if self._cb_species not in ions:
raise ValueError(
f"Charge balancing species {self._cb_species} was not found in the solution!. "
f"Species {ions} were found."
)
# adjust charge balance, if necessary
self._adjust_charge_balance()
# validate water stability, if necessary
self._check_water_stability()
@property
def mass(self) -> Quantity:
"""
Return the total mass of the solution.
The mass is calculated each time this method is called.
Returns: The mass of the solution, in kg
"""
mass = np.sum([self.get_amount(item, "kg").magnitude for item in self.components])
return mass * ureg.kg
@property
def solvent_mass(self) -> Quantity:
"""
Return the mass of the solvent.
This property is used whenever mol/kg (or similar) concentrations
are requested by get_amount()
Returns:
The mass of the solvent, in kg
See Also:
:py:meth:`get_amount()`
"""
return self.get_amount(self.solvent, "kg")
@property
def volume(self) -> Quantity:
"""
Return the volume of the solution.
Returns:
Quantity: the volume of the solution, in L
"""
# if the composition has changed, recalculate the volume first
if self.volume_update_required is True:
self._update_volume()
self.volume_update_required = False
return self._volume.to("L")
@volume.setter
def volume(self, volume: str):
"""Change the total solution volume to volume, while preserving
all component concentrations.
Args:
volume : Total volume of the solution, including the unit, e.g. '1 L'
Examples:
>>> mysol = Solution([['Na+','2 mol/L'],['Cl-','0.01 mol/L']],volume='500 mL')
>>> print(mysol.volume)
0.5000883925072983 l
"""
# figure out the factor to multiply the old concentrations by
scale_factor = ureg.Quantity(volume) / self.volume
# scale down the amount of all the solutes according to the factor
for solute in self.components:
self.components[solute] *= scale_factor.magnitude
# update the solution volume
self._volume *= scale_factor.magnitude
@property
def temperature(self) -> Quantity:
"""Return the temperature of the solution in Kelvin."""
return self._temperature.to("K")
@temperature.setter
def temperature(self, temperature: str):
"""
Set the solution temperature.
Args:
temperature: pint-compatible string, e.g. '25 degC'
"""
self._temperature = ureg.Quantity(temperature)
# update the water substance
self.water_substance = create_water_substance(self.temperature, self.pressure)
# recalculate the volume
self.volume_update_required = True
# clear any cached solute properties that may depend on temperature
self.get_property.cache_clear()
self.get_molar_conductivity.cache_clear()
self.get_mobility.cache_clear()
self.get_diffusion_coefficient.cache_clear()
@property
def pressure(self) -> Quantity:
"""Return the hydrostatic pressure of the solution in atm."""
return self._pressure.to("atm")
@pressure.setter
def pressure(self, pressure: str):
"""
Set the solution pressure.
Args:
pressure: pint-compatible string, e.g. '1.2 atmC'
"""
self._pressure = ureg.Quantity(pressure)
# update the water substance
self.water_substance = create_water_substance(self.temperature, self.pressure)
# recalculate the volume
self.volume_update_required = True
@property
def pH(self) -> float:
"""Return the pH of the solution."""
return self.p("H+", activity=False)
[docs]
def p(self, solute: str, activity=True) -> float:
"""
Return the negative log of the activity of solute.
Generally used for expressing concentration of hydrogen ions (pH)
Args:
solute : str
String representing the formula of the solute
activity: bool, optional
If False, the function will use the molar concentration rather
than the activity to calculate p. Defaults to True.
Returns:
Quantity
The negative log10 of the activity (or molar concentration if
activity = False) of the solute. If the solute has zero concentration
then np.nan (not a number) is returned.
"""
try:
if activity is True:
amt = self.get_activity(solute).magnitude
else:
amt = self.get_amount(solute, "mol/L").magnitude
return float(-1 * np.log10(amt))
# if the solute has zero or negative concentration, np.log10 raises a RuntimeWarning
except RuntimeWarning:
return np.nan
@property
def density(self) -> Quantity:
"""
Return the density of the solution.
Density is calculated from the mass and volume each time this method is called.
Returns:
Quantity: The density of the solution.
"""
return self.mass / self.volume
@property
def dielectric_constant(self) -> Quantity:
r"""
Returns the dielectric constant of the solution.
Returns:
Quantity: the dielectric constant of the solution, dimensionless.
Notes:
Implements the following equation as given by Zuber et al.
.. math:: \epsilon = \epsilon_{solvent} \over 1 + \sum_i \alpha_i x_i
where :math:`\alpha_i` is a coefficient specific to the solvent and ion, and :math:`x_i`
is the mole fraction of the ion in solution.
References:
A. Zuber, L. Cardozo-Filho, V.F. Cabral, R.F. Checoni, M. Castier,
An empirical equation for the dielectric constant in aqueous and nonaqueous
electrolyte mixtures, Fluid Phase Equilib. 376 (2014) 116-123.
doi:10.1016/j.fluid.2014.05.037.
"""
di_water = self.water_substance.epsilon
denominator = 1
for item in self.components:
# ignore water
if item != "H2O(aq)":
# skip over solutes that don't have parameters
# try:
fraction = self.get_amount(item, "fraction")
coefficient = self.get_property(item, "model_parameters.dielectric_zuber")
if coefficient is not None:
denominator += coefficient * fraction
# except TypeError:
# self.logger.warning("No dielectric parameters found for species %s." % item)
# continue
return ureg.Quantity(di_water / denominator)
@property
def chemical_system(self) -> str:
"""
Return the chemical system of the Solution as a "-" separated list of elements, sorted alphabetically. For
example, a solution containing CaCO3 would have a chemical system of "C-Ca-H-O".
"""
return "-".join(self.elements)
@property
def elements(self) -> list:
"""
Return a list of elements that are present in the solution.
For example, a solution containing CaCO3 would return ["C", "Ca", "H", "O"]
"""
els = []
for s in self.components:
els.extend(self.get_property(s, "elements"))
return sorted(set(els))
@property
def cations(self) -> dict[str, float]:
"""
Returns the subset of `components` that are cations.
The returned dict is sorted by amount in descending order.
"""
return {k: v for k, v in self.components.items() if self.get_property(k, "charge") > 0}
@property
def anions(self) -> dict[str, float]:
"""
Returns the subset of `components` that are anions.
The returned dict is sorted by amount in descending order.
"""
return {k: v for k, v in self.components.items() if self.get_property(k, "charge") < 0}
@property
def neutrals(self) -> dict[str, float]:
"""
Returns the subset of `components` that are neutral (not charged).
The returned dict is sorted by amount in descending order.
"""
return {k: v for k, v in self.components.items() if self.get_property(k, "charge") == 0}
# TODO - need tests for viscosity
@property
def viscosity_dynamic(self) -> Quantity:
"""
Return the dynamic (absolute) viscosity of the solution.
Calculated from the kinematic viscosity
See Also:
:attr:`viscosity_kinematic`
"""
return self.viscosity_kinematic * self.density
# TODO - before deprecating get_viscosity_relative, consider whether the Jones-Dole
# model should be integrated here as a fallback, in case salt parameters for the
# other model are not available.
# if self.ionic_strength.magnitude > 0.2:
# self.logger.warning('Viscosity calculation has limited accuracy above 0.2m')
# viscosity_rel = 1
# for item in self.components:
# # ignore water
# if item != 'H2O':
# # skip over solutes that don't have parameters
# try:
# conc = self.get_amount(item,'mol/kg').magnitude
# coefficients= self.get_property(item, 'jones_dole_viscosity')
# viscosity_rel += coefficients[0] * conc ** 0.5 + coefficients[1] * conc + \
# coefficients[2] * conc ** 2
# except TypeError:
# continue
# return (
# self.viscosity_dynamic / self.water_substance.mu * ureg.Quantity("1 Pa*s")
# )
@property
def viscosity_kinematic(self) -> Quantity:
r"""
Return the kinematic viscosity of the solution.
Notes:
The calculation is based on a model derived from the Eyring equation
and presented in
.. math::
\ln \nu = \ln {\nu_w MW_w \over \sum_i x_i MW_i } +
15 x_+^2 + x_+^3 \delta G^*_{123} + 3 x_+ \delta G^*_{23} (1-0.05x_+)
Where:
.. math:: \delta G^*_{123} = a_o + a_1 (T)^{0.75}
.. math:: \delta G^*_{23} = b_o + b_1 (T)^{0.5}
In which :math:`\nu` is the kinematic viscosity, MW is the molecular weight,
:math:`x_{+}` is the mole fraction of cations, and :math:`T` is the temperature in degrees C.
The a and b fitting parameters for a variety of common salts are included in the
database.
References:
Vásquez-Castillo, G.; Iglesias-Silva, G. a.; Hall, K. R. An extension of the McAllister model to correlate
kinematic viscosity of electrolyte solutions. Fluid Phase Equilib. 2013, 358, 44-49.
See Also:
:py:meth:`viscosity_dynamic`
"""
# identify the main salt in the solution
salt = self.get_salt()
a0 = a1 = b0 = b1 = 0
# retrieve the parameters for the delta G equations
params = None if salt is None else self.get_property(salt.formula, "model_parameters.viscosity_eyring")
if params is not None:
a0 = ureg.Quantity(params["a0"]["value"]).magnitude
a1 = ureg.Quantity(params["a1"]["value"]).magnitude
b0 = ureg.Quantity(params["b0"]["value"]).magnitude
b1 = ureg.Quantity(params["b1"]["value"]).magnitude
# compute the delta G parameters
temperature = self.temperature.to("degC").magnitude
G_123 = a0 + a1 * (temperature) ** 0.75
G_23 = b0 + b1 * (temperature) ** 0.5
# calculate the cation mole fraction
# x_cat = self.get_amount(cation, "fraction")
x_cat = self.get_amount(salt.cation, "fraction").magnitude
else:
# TODO - fall back to the Jones-Dole model! There are currently no eyring parameters in the database!
# proceed with the coefficients equal to zero and log a warning
self.logger.warning("Appropriate viscosity coefficients were not found. Viscosity will be approximate.")
G_123 = G_23 = 0
x_cat = 0
# get the kinematic viscosity of water, returned by IAPWS in m2/s
nu_w = self.water_substance.nu
# compute the effective molar mass of the solution
total_moles = np.sum([v for k, v in self.components.items()])
MW = self.mass.to("g").magnitude / total_moles
# get the MW of water
MW_w = self.get_property(self.solvent, "molecular_weight").magnitude
# calculate the kinematic viscosity
nu = np.log(nu_w * MW_w / MW) + 15 * x_cat**2 + x_cat**3 * G_123 + 3 * x_cat * G_23 * (1 - 0.05 * x_cat)
return np.exp(nu) * ureg.m**2 / ureg.s
@property
def conductivity(self) -> Quantity:
r"""
Compute the electrical conductivity of the solution.
Returns:
The electrical conductivity of the solution in Siemens / meter.
Notes:
Conductivity is calculated by summing the molar conductivities of the respective
solutes.
.. math::
EC = {F^2 \over R T} \sum_i D_i z_i ^ 2 m_i = \sum_i \lambda_i m_i
Where :math:`D_i` is the diffusion coefficient, :math:`m_i` is the molal concentration,
:math:`z_i` is the charge, and the summation extends over all species in the solution.
Alternatively, :math:`\lambda_i` is the molar conductivity of solute i.
Diffusion coefficients :math:`D_i` (and molar conductivities :math:`\lambda_i`) are
adjusted for the effects of temperature and ionic strength using the method implemented
in PHREEQC >= 3.4. [aq]_ [hc]_ See `get_diffusion_coefficient for` further details.
References:
.. [aq] https://www.aqion.de/site/electrical-conductivity
.. [hc] https://www.hydrochemistry.eu/exmpls/sc.html
See Also:
:py:attr:`ionic_strength`
:py:meth:`get_diffusion_coefficient`
:py:meth:`get_molar_conductivity`
"""
EC = (
np.asarray(
[
self.get_molar_conductivity(i).to("S*L/mol/m").magnitude * self.get_amount(i, "mol/L").magnitude
for i in self.components
]
)
* ureg.S
/ ureg.m
)
return np.sum(EC)
@property
def ionic_strength(self) -> Quantity:
r"""
Return the ionic strength of the solution.
Return the ionic strength of the solution, calculated as 1/2 * sum ( molality * charge ^2) over all the ions.
Molal (mol/kg) scale concentrations are used for compatibility with the activity correction formulas.
Returns:
Quantity:
The ionic strength of the parent solution, mol/kg.
See Also:
:py:meth:`get_activity`
:py:meth:`get_water_activity`
Notes:
The ionic strength is calculated according to:
.. math:: I = \sum_i m_i z_i^2
Where :math:`m_i` is the molal concentration and :math:`z_i` is the charge on species i.
Examples:
>>> s1 = pyEQL.Solution([['Na+','0.2 mol/kg'],['Cl-','0.2 mol/kg']])
>>> s1.ionic_strength
<Quantity(0.20000010029672785, 'mole / kilogram')>
>>> s1 = pyEQL.Solution([['Mg+2','0.3 mol/kg'],['Na+','0.1 mol/kg'],['Cl-','0.7 mol/kg']],temperature='30 degC')
>>> s1.ionic_strength
<Quantity(1.0000001004383303, 'mole / kilogram')>
"""
# compute using magnitudes only, for performance reasons
ionic_strength = np.sum(
[mol * self.get_property(solute, "charge") ** 2 for solute, mol in self.components.items()]
)
ionic_strength /= self.solvent_mass.to("kg").magnitude # convert to mol/kg
ionic_strength *= 0.5
return ionic_strength * ureg.mol / ureg.kg
@property
def charge_balance(self) -> float:
r"""
Return the signed charge balance of the solution, positive or negative.
Return the signed charge balance of the solution, positive or negative. The charge balance represents the net electric charge
of the solution and SHOULD equal zero at all times, but due to numerical errors will usually have a small nonzero value.
Positive values indicate excess cationic charge, while negative values indivate excess anionic charge. It is calculated according to:
.. math:: CB = \sum_i C_i z_i
where :math:`C_i` is the molar concentration, and :math:`z_i` is the charge on species i.
Returns:
float :
The signed charge balance of the solution, in equivalents (mol of charge) per L.
"""
charge_balance = 0
for solute in self.components:
charge_balance += self.get_amount(solute, "eq/L").magnitude
return charge_balance
# TODO - consider adding guard statements to prevent alkalinity from being negative
@property
def alkalinity(self) -> Quantity:
r"""
Return the alkalinity or acid neutralizing capacity of a solution.
Returns:
Quantity: The alkalinity of the solution in mg/L as CaCO3
Notes:
The alkalinity is calculated according to [stm]_
.. math:: Alk = \sum_{i} z_{i} C_{B} + \sum_{i} z_{i} C_{A}
Where :math:`C_{B}` and :math:`C_{A}` are conservative cations and anions, respectively
(i.e. ions that do not participate in acid-base reactions), and :math:`z_{i}` is their signed charge.
When conservative cations (Group I and II cations) or strong base anions are present, the alkalinity is calculated according to[stm]_
.. math:: Alk = \sum_{i} z_{i} C_{B} + \sum_{i} z_{i} C_{A}
Where :math:`C_{B}` and :math:`C_{A}` are conservative cations and strong base anions, respectively (i.e. ions that do not participate in acid-base reactions), and :math:`z_{i}` is their signed charge.
Alternatively, if those species are not present, then alkalinity is calculated based on the concentrations of weak acid and base species according to [stm]_
.. math:: Alk = -\sum_{i} z_{i} C_{i}
Where :math:`C_i` is the molar concentration of species i, and :math:`z_i` is its charge.
The summation should extend over all weak inorganic species that can participate in acid-base reactions. In this method, we consider HCO3[-1], CO3[-2], H2PO4[-1], HPO4[-2], PO4[-3], HS[-1], S[-2], H3SiO4[-1], H2SiO4[-2], B(OH)4[-1], NH3(aq), OH[-1], and H[+1] as the relevant weak acid/base species, while organics are excluded.
References:
.. [stm] Stumm, Werner and Morgan, James J. Aquatic Chemistry, 3rd ed, pp 165. Wiley Interscience, 1996.
"""
alkalinity = 0 * ureg.mol / ureg.L
base_cations = {
"Li[+1]",
"Na[+1]",
"K[+1]",
"Rb[+1]",
"Cs[+1]",
"Fr[+1]",
"Be[+2]",
"Mg[+2]",
"Ca[+2]",
"Sr[+2]",
"Ba[+2]",
"Ra[+2]",
}
acid_anions = {
"Cl[-1]",
"Br[-1]",
"I[-1]",
"SO4[-2]",
"NO3[-1]",
"ClO4[-1]",
"ClO3[-1]",
}
weak_species = {
"HCO3[-1]",
"CO3[-2]",
"H2PO4[-1]",
"HPO4[-2]",
"PO4[-3]",
"HS[-1]",
"S[-2]",
"H3SiO4[-1]",
"H2SiO4[-2]",
"B(OH)4[-1]",
"NH3(aq)",
"OH[-1]",
"H[+1]",
} # Note that organics are excluded
conservative_species = base_cations.union(acid_anions)
# check presence of conservative cations or strong base anions
conservative_def = any(item in conservative_species for item in self.components)
for item in self.components:
if item in conservative_species:
# Conservative cations and strong base anions
alkalinity += self.get_amount(item, "eq/L")
elif item in weak_species and not conservative_def:
# Weak acid/base species, exclude organics
alkalinity += self.get_amount(item, "eq/L") * (-1)
return (alkalinity * EQUIV_WT_CACO3).to("mg/L")
@property
def hardness(self) -> Quantity:
"""
Return the hardness of a solution.
Hardness is defined as the sum of the equivalent concentrations
of multivalent cations as calcium carbonate.
NOTE: at present pyEQL cannot distinguish between mg/L as CaCO3
and mg/L units. Use with caution.
Returns:
Quantity:
The hardness of the solution in mg/L as CaCO3
"""
hardness = 0 * ureg.mol / ureg.L
for item in self.components:
z = self.get_property(item, "charge")
if z > 1:
hardness += z * self.get_amount(item, "mol/L")
# convert the hardness to mg/L as CaCO3
return (hardness * EQUIV_WT_CACO3).to("mg/L")
@property
def total_dissolved_solids(self) -> Quantity:
"""
Total dissolved solids in mg/L (equivalent to ppm) including both charged and uncharged species.
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 = 0 * ureg.mg / ureg.L
for s in self.components:
# ignore pure water and dissolved gases, but not CO2
if s in ["H2O(aq)", "H[+1]", "OH[-1]"]:
continue
tds += self.get_amount(s, "mg/L")
return tds
@property
def TDS(self) -> Quantity:
"""Alias of :py:meth:`total_dissolved_solids`."""
return self.total_dissolved_solids
@property
def debye_length(self) -> Quantity:
r"""
Return the Debye length of a solution.
Debye length is calculated as [wk3]_
.. math::
\kappa^{-1} = \sqrt({\epsilon_r \epsilon_o k_B T \over (2 N_A e^2 I)})
where :math:`I` is the ionic strength, :math:`\epsilon_r` and :math:`\epsilon_r`
are the relative permittivity and vacuum permittivity, :math:`k_B` is the
Boltzmann constant, and :math:`T` is the temperature, :math:`e` is the
elementary charge, and :math:`N_A` is Avogadro's number.
Returns The Debye length, in nanometers.
References:
.. [wk3] https://en.wikipedia.org/wiki/Debye_length#In_an_electrolyte_solution
See Also:
:attr:`ionic_strength`
:attr:`dielectric_constant`
"""
# to preserve dimensionality, convert the ionic strength into mol/L units
ionic_strength = self.ionic_strength.magnitude * ureg.mol / ureg.L
dielectric_constant = self.dielectric_constant
debye_length = (
dielectric_constant
* ureg.epsilon_0
* ureg.k
* self.temperature
/ (2 * ureg.N_A * ureg.e**2 * ionic_strength)
) ** 0.5
return debye_length.to("nm")
@property
def bjerrum_length(self) -> Quantity:
r"""
Return the Bjerrum length of a solution.
Bjerrum length represents the distance at which electrostatic
interactions between particles become comparable in magnitude
to the thermal energy.:math:`\lambda_B` is calculated as
.. math::
\lambda_B = {e^2 \over (4 \pi \epsilon_r \epsilon_o k_B T)}
where :math:`e` is the fundamental charge, :math:`\epsilon_r` and :math:`\epsilon_r`
are the relative permittivity and vacuum permittivity, :math:`k_B` is the
Boltzmann constant, and :math:`T` is the temperature.
Returns:
Quantity:
The Bjerrum length, in nanometers.
References:
https://en.wikipedia.org/wiki/Bjerrum_length
Examples:
>>> s1 = pyEQL.Solution()
>>> s1.bjerrum_length
<Quantity(0.7152793009386953, 'nanometer')>
See Also:
:attr:`dielectric_constant`
"""
bjerrum_length = ureg.e**2 / (4 * np.pi * self.dielectric_constant * ureg.epsilon_0 * ureg.k * self.temperature)
return bjerrum_length.to("nm")
@property
def osmotic_pressure(self) -> Quantity:
r"""
Return the osmotic pressure of the solution relative to pure water.
Returns:
The osmotic pressure of the solution relative to pure water in Pa
See Also:
:attr:`get_water_activity`
:attr:`get_osmotic_coefficient`
:attr:`get_salt`
Notes:
Osmotic pressure is calculated based on the water activity [sata]_ [wk]_
.. math:: \Pi = -\frac{RT}{V_{w}} \ln a_{w}
Where :math:`\Pi` is the osmotic pressure, :math:`V_{w}` is the partial
molar volume of water (18.2 cm**3/mol), and :math:`a_{w}` is the water
activity.
References:
.. [sata] Sata, Toshikatsu. Ion Exchange Membranes: Preparation, Characterization, and Modification.
Royal Society of Chemistry, 2004, p. 10.
.. [wk] https://en.wikipedia.org/wiki/Osmotic_pressure#Derivation_of_the_van_'t_Hoff_formula
Examples:
>>> s1=pyEQL.Solution()
>>> s1.osmotic_pressure
<Quantity(0.495791416, 'pascal')>
>>> s1 = pyEQL.Solution([['Na+','0.2 mol/kg'],['Cl-','0.2 mol/kg']])
>>> soln.osmotic_pressure
<Quantity(906516.7318131207, 'pascal')>
"""
partial_molar_volume_water = self.get_property(self.solvent, "size.molar_volume")
osmotic_pressure = (
-1 * ureg.R * self.temperature / partial_molar_volume_water * np.log(self.get_water_activity())
)
self.logger.debug(
f"Calculated osmotic pressure of solution as {osmotic_pressure} Pa at T= {self.temperature} degrees C"
)
return osmotic_pressure.to("Pa")
# Concentration Methods
[docs]
def get_amount(self, solute: str, units: str = "mol/L") -> Quantity:
"""
Return the amount of 'solute' in the parent solution.
The amount of a solute can be given in a variety of unit types.
1. substance per volume (e.g., 'mol/L', 'M')
2. equivalents (i.e., moles of charge) per volume (e.g., 'eq/L', 'meq/L')
3. substance per mass of solvent (e.g., 'mol/kg', 'm')
4. mass of substance (e.g., 'kg')
5. moles of substance ('mol')
6. mole fraction ('fraction')
7. percent by weight (%)
8. number of molecules ('count')
9. "parts-per-x" units, where ppm = mg/L, ppb = ug/L ppt = ng/L
Args:
solute : str
String representing the name of the solute of interest
units : str
Units desired for the output. Examples of valid units are
'mol/L','mol/kg','mol', 'kg', and 'g/L'
Use 'fraction' to return the mole fraction.
Use '%' to return the mass percent
Returns:
The amount of the solute in question, in the specified units
See Also:
:attr:`mass`
:meth:`add_amount`
:meth:`set_amount`
:meth:`get_total_amount`
:meth:`get_osmolarity`
:meth:`get_osmolality`
:meth:`get_total_moles_solute`
:func:`pyEQL.utils.translate_units`
"""
z = 1
# sanitized unit to be passed to pint
if "eq" in units:
_units = units.replace("eq", "mol")
z = self.get_property(solute, "charge")
if z == 0: # uncharged solutes have zero equiv concentration
return ureg.Quantity(0, _units)
else:
_units = translate_units(units)
# retrieve the number of moles of solute and its molecular weight
try:
moles = self.components[solute] * ureg.mol
# if the solute is not present in the solution, we'll get a KeyError
# In that case, the amount is zero
except KeyError:
try:
return ureg.Quantity(0, _units)
except DimensionalityError:
self.logger.error(
f"Unsupported unit {units} specified for zero-concentration solute {solute}. Returned 0."
)
return ureg.Quantity(0)
# with pint unit conversions enabled, we just pass the unit to pint
# the logic tests here ensure that only the required arguments are
# passed to pint for the unit conversion. This avoids unnecessary
# function calls.
if units == "count":
return round((moles * ureg.N_A).to("dimensionless"), 0)
if units == "fraction":
return moles / (self.get_moles_solvent() + self.get_total_moles_solute())
mw = self.get_property(solute, "molecular_weight").to("g/mol")
if units == "%":
return moles.to("kg", "chem", mw=mw) / self.mass.to("kg") * 100
qty = ureg.Quantity(_units)
if _units in ["eq", "mol", "moles"] or qty.check("[substance]"):
return z * moles.to(_units)
if (
_units in ["mol/L", "eq/L", "g/L", "mg/L", "ug/L"]
or qty.check("[substance]/[length]**3")
or qty.check("[mass]/[length]**3")
):
return z * moles.to(_units, "chem", mw=mw, volume=self.volume)
if _units in ["mol/kg"] or qty.check("[substance]/[mass]") or qty.check("[mass]/[mass]"):
return z * moles.to(_units, "chem", mw=mw, solvent_mass=self.solvent_mass)
if _units in ["kg", "g"] or qty.check("[mass]"):
return moles.to(_units, "chem", mw=mw)
raise ValueError(f"Unsupported unit {units} specified for get_amount")
[docs]
def get_components_by_element(
self, nested: bool = False
) -> dict[str, list[str]] | dict[str, dict[float | str, list[str]]]:
"""
Return a list of all species associated with a given element.
Args:
nested : bool
Whether to return a nested dictionary of <element>
to <valence> => <list of species> mapping. False by default.
Returns:
A mapping of element to a list of species in the solution.
If nested is False (default), elements (keys) are suffixed with
their oxidation state in parentheses, e.g.,
{"Na(1.0)":["Na[+1]", "NaOH(aq)"]}
If nested is True, the dictionary is nested, e.g.,
{"Na": [{1:["Na[+1]", "NaOH(aq)"]}]}.
Note that the valence may be a string, assuming the value "unk"
denoting an unknown oxidation state.
Species associated with each element are sorted in descending order of the amount
present (i.e., the first species listed is the most abundant).
"""
d = {}
# by sorting the components according to amount, we ensure that the species
# are sorted in descending order of concentration in the resulting dict
for s in self.components:
# determine the element and oxidation state
elements = self.get_property(s, "elements")
for el in elements:
try:
oxi_states = self.get_property(s, "oxi_state_guesses")
oxi_state = oxi_states.get(el, UNKNOWN_OXI_STATE)
except (TypeError, IndexError):
self.logger.error(f"No oxidation state found for element {el}. Assigning '{UNKNOWN_OXI_STATE}'")
oxi_state = UNKNOWN_OXI_STATE
if d.get(el):
if d[el].get(oxi_state):
d[el][oxi_state].append(s)
else:
d[el][oxi_state] = [s]
else:
d[el] = {oxi_state: [s]}
if nested:
return d
return {f"{el}({val})": species for el, val_dict in d.items() for val, species in val_dict.items()}
[docs]
def get_el_amt_dict(self, nested: bool = False) -> dict[str, float] | dict[str, dict[float | str, float]]:
"""
Return a dict of Element: amount in mol.
Args:
nested : bool
Whether to return a nested dictionary of <element>
to <valence> => amount mapping. False by default.
Returns:
A mapping of element to its amount in moles in the solution.
If nested is False (default), elements (keys) are suffixed with
their oxidation state in parentheses, e.g.,
{"Fe(2.0)": 0.354, "Cl(-1.0)": 0.708}
If nested is True, the dictionary is nested, e.g.,
{"Fe": {2.0: 0.354}, "Cl": {-1.0: 0.708}}.}
Note that the valence may be a string, assuming the value "unk"
denoting an unknown oxidation state.
"""
d = {}
for s, mol in self.components.items():
elements = self.get_property(s, "elements")
pmg_ion_dict = self.get_property(s, "pmg_ion")
oxi_states = self.get_property(s, "oxi_state_guesses")
for el in elements:
# stoichiometric coefficient, mol element per mol solute
stoich = pmg_ion_dict.get(el)
try:
oxi_states = self.get_property(s, "oxi_state_guesses")
oxi_state = oxi_states.get(el, UNKNOWN_OXI_STATE)
except (TypeError, IndexError):
self.logger.error(f"No oxidation state found for element {el}. Assigning '{UNKNOWN_OXI_STATE}'")
oxi_state = UNKNOWN_OXI_STATE
if d.get(el):
if d[el].get(oxi_state):
d[el][oxi_state] += stoich * mol
else:
d[el][oxi_state] = stoich * mol
else:
d[el] = {oxi_state: stoich * mol}
if nested:
return d
return {f"{el}({val})": amount for el, val_dict in d.items() for val, amount in val_dict.items()}
[docs]
def get_total_amount(self, element: str, units: str) -> Quantity:
"""
Return the total amount of 'element' (across all solutes) in the solution.
Args:
element: The symbol of the element of interest. The symbol can optionally be followed by the
oxidation state in parentheses, e.g., "Na(1.0)", "Fe(2.0)", or "O(0.0)". If no oxidation state
is given, the total concentration of the element (over all oxidation states) is returned.
units : str
Units desired for the output. Any unit understood by `get_amount` can be used. Examples of valid
units are 'mol/L','mol/kg','mol', 'kg', and 'g/L'.
Returns:
The total amount of the element in the solution, in the specified units
See Also:
:meth:`get_amount`
:func:`pyEQL.utils.translate_units`
"""
_units = translate_units(units)
TOT: Quantity = ureg.Quantity(0, _units)
# standardize the element formula and units
el = str(Element(element.split("(")[0]))
units = translate_units(units)
# enumerate the species whose concentrations we need
comp_by_element = self.get_components_by_element()
# compile list of species in different ways depending whether there is an oxidation state
if "(" in element and UNKNOWN_OXI_STATE not in element:
ox = float(element.split("(")[-1].split(")")[0])
key = f"{el}({ox})"
species = comp_by_element.get(key, [])
else:
species = []
for k, v in comp_by_element.items():
if k.split("(")[0] == el:
species.extend(v)
# loop through the species of interest, adding moles of element
for item, amt in self.components.items():
if item in species:
amt = self.get_amount(item, units)
ion = Ion.from_formula(item)
# convert the solute amount into the amount of element by
# either the mole / mole or weight ratio
if ureg.Quantity(units).dimensionality in (
"[substance]",
"[substance]/[length]**3",
"[substance]/[mass]",
):
TOT += amt * ion.get_el_amt_dict()[el] # returns {el: mol per formula unit}
elif ureg.Quantity(units).dimensionality in (
"[mass]",
"[mass]/[length]**3",
"[mass]/[mass]",
):
TOT += amt * ion.to_weight_dict[el] # returns {el: wt fraction}
return TOT
[docs]
def add_solute(self, formula: str, amount: str):
"""Primary method for adding substances to a pyEQL solution.
Args:
formula (str): Chemical formula for the solute. Charged species must contain a+ or - and
(for polyvalent solutes) a number representing the net charge (e.g. 'SO4-2').
amount (str): The amount of substance in the specified unit system. The string should
contain both a quantity and a pint-compatible representation of a ureg. e.g. '5 mol/kg' or '0.1 g/L'.
"""
Q = ureg.Quantity(*_translate_pint_quantity(amount))
# if units are given on a per-volume basis,
# iteratively solve for the amount of solute that will preserve the
# original volume and result in the desired concentration
if Q.dimensionality in (
"[substance]/[length]**3",
"[mass]/[length]**3",
):
# store the original volume for later
orig_volume = self.volume
# add the new solute
mw = self.get_property(formula, "molecular_weight") # returns a quantity
target_mol = Q.to("moles", "chem", mw=mw, volume=self.volume, solvent_mass=self.solvent_mass)
self.components[formula] = target_mol.to("moles").magnitude
# calculate the volume occupied by all the solutes
solute_vol = self._get_solute_volume()
# determine the volume of solvent that will preserve the original volume
target_vol = orig_volume - solute_vol
# adjust the amount of solvent
# density is returned in kg/m3 = g/L
target_mass = target_vol * self.water_substance.rho * ureg.g / ureg.L
mw = self.get_property(self.solvent, "molecular_weight")
if mw is None:
raise ValueError(f"Molecular weight for solvent {self.solvent} not found in database. Cannot proceed.")
target_mol = target_mass.to("g") / mw.to("g/mol")
self.components[self.solvent] = target_mol.magnitude
else:
if Q.dimensionality == "[substance]":
# add the new solute directly
self.components[formula] = Q.to("moles").magnitude
else:
# add the new solute after converting to moles
mw = self.get_property(formula, "molecular_weight") # returns a quantity
target_mol = Q.to("moles", "chem", mw=mw, volume=self.volume, solvent_mass=self.solvent_mass)
self.components[formula] = target_mol.to("moles").magnitude
# 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.magnitude <= 0:
self.logger.error("All solvent has been depleted from the solution")
return
# set the volume recalculation flag
self.volume_update_required = True
[docs]
def add_amount(self, solute: str, amount: str):
"""
Add the amount of 'solute' to the parent solution.
Args:
solute : str
String representing the name of the solute of interest
amount : str quantity
String representing the concentration desired, e.g. '1 mol/kg'
If the units are given on a per-volume basis, the solution
volume is not recalculated
If the units are given on a mass, substance, per-mass, or
per-substance basis, then the solution volume is recalculated
based on the new composition
Returns:
Nothing. The concentration of solute is modified.
"""
Q = ureg.Quantity(*_translate_pint_quantity(amount))
# Get the current amount of the solute
current_amt = self.get_amount(solute, amount.split(" ")[1])
if current_amt.magnitude == 0:
self.logger.warning(f"Add new solute {solute} to the solution")
new_amt = Q + current_amt
self.set_amount(solute, new_amt)
[docs]
def set_amount(self, solute: str, amount: str):
"""
Set the amount of 'solute' in the parent solution.
Args:
solute : str
String representing the name of the solute of interest
amount : str Quantity
String representing the concentration desired, e.g. '1 mol/kg'
If the units are given on a per-volume basis, the solution
volume is not recalculated and the molar concentrations of
other components in the solution are not altered, while the
molal concentrations are modified.
If the units are given on a mass, substance, per-mass, or
per-substance basis, then the solution volume is recalculated
based on the new composition and the molal concentrations of
other components are not altered, while the molar concentrations
are modified.
Returns:
Nothing. The concentration of solute is modified.
"""
Q = ureg.Quantity(*_translate_pint_quantity(amount))
# raise an error if a negative amount is specified
if Q.magnitude < 0:
raise ValueError(f"Negative amount specified for solute {solute}. Concentration not changed.")
# if units are given on a per-volume basis,
# iteratively solve for the amount of solute that will preserve the
# original volume and result in the desired concentration
if Q.dimensionality in (
"[substance]/[length]**3",
"[mass]/[length]**3",
):
# store the original volume for later
orig_volume = self.volume
# change the amount of the solute present to match the desired amount
self.components[solute] = Q.to(
"moles",
"chem",
mw=ureg.Quantity(self.get_property(solute, "molecular_weight")),
volume=self.volume,
solvent_mass=self.solvent_mass,
).magnitude
# calculate the volume occupied by all the solutes
solute_vol = self._get_solute_volume()
# determine the volume of solvent that will preserve the original volume
target_vol = orig_volume - solute_vol
# adjust the amount of solvent
target_mass = target_vol * self.water_substance.rho * ureg.g / ureg.L
mw = self.get_property(self.solvent, "molecular_weight")
target_mol = target_mass / mw
self.components[self.solvent] = target_mol.to("mol").magnitude
else:
# change the amount of the solute present
self.components[solute] = Q.to(
"moles",
"chem",
mw=ureg.Quantity(self.get_property(solute, "molecular_weight")),
volume=self.volume,
solvent_mass=self.solvent_mass,
).magnitude
# make sure that there is still solvent present in the first place
if self.solvent_mass.magnitude <= 0:
self.logger.critical("All solvent has been depleted from the solution")
return
# update the volume to account for the space occupied by all the solutes
self.volume_update_required = True
[docs]
def get_total_moles_solute(self) -> Quantity:
"""Return the total moles of all solute in the solution."""
tot_mol = 0
for item in self.components:
if item != self.solvent:
tot_mol += self.components[item]
return tot_mol * ureg.mol
[docs]
def get_moles_solvent(self) -> Quantity:
"""
Return the moles of solvent present in the solution.
Returns:
The moles of solvent in the solution.
"""
return self.get_amount(self.solvent, "mol")
[docs]
def get_osmolarity(self, activity_correction=False) -> Quantity:
"""Return the osmolarity of the solution in Osm/L.
Args:
activity_correction : bool
If TRUE, the osmotic coefficient is used to calculate the
osmolarity. This correction is appropriate when trying to predict
the osmolarity that would be measured from e.g. freezing point
depression. Defaults to FALSE if omitted.
"""
factor = self.get_osmotic_coefficient() if activity_correction is True else 1
return factor * self.get_total_moles_solute() / self.volume.to("L")
[docs]
def get_osmolality(self, activity_correction=False) -> Quantity:
"""Return the osmolality of the solution in Osm/kg.
Args:
activity_correction : bool
If TRUE, the osmotic coefficient is used to calculate the
osmolarity. This correction is appropriate when trying to predict
the osmolarity that would be measured from e.g. freezing point
depression. Defaults to FALSE if omitted.
"""
factor = self.get_osmotic_coefficient() if activity_correction is True else 1
return factor * self.get_total_moles_solute() / self.solvent_mass.to("kg")
[docs]
def get_salt(self) -> Salt:
"""
Determine the predominant salt in a solution of ions.
Many empirical equations for solution properties such as activity coefficient,
partial molar volume, or viscosity are based on the concentration of
single salts (e.g., NaCl). When multiple ions are present (e.g., a solution
containing Na+, Cl-, and Mg+2), it is generally not possible to directly model
these quantities. pyEQL works around this problem by treating such solutions
as single salt solutions.
The get_salt() method examines the ionic composition of a solution and returns
an object that identifies the single most predominant salt in the solution, defined
by the cation and anion with the highest mole fraction. The Salt object contains
information about the stoichiometry of the salt to enable its effective concentration
to be calculated (e.g., if a solution contains 0.5 mol/kg of Na+ and Cl-, plus traces
of H+ and OH-, the matched salt is 0.5 mol/kg NaCl).
Returns:
Salt object containing information about the parent salt.
See Also:
:py:meth:`get_activity`
:py:meth:`get_activity_coefficient`
:py:meth:`get_water_activity`
:py:meth:`get_osmotic_coefficient`
:py:attr:`osmotic_pressure`
:py:attr:`viscosity_kinematic`
Examples:
>>> s1 = Solution([['Na+','0.5 mol/kg'],['Cl-','0.5 mol/kg']])
>>> s1.get_salt()
<pyEQL.salt_ion_match.Salt object at 0x7fe6d3542048>
>>> s1.get_salt().formula
'NaCl'
>>> s1.get_salt().nu_cation
1
>>> s1.get_salt().z_anion
-1
>>> s2 = pyEQL.Solution([['Na+','0.1 mol/kg'],['Mg+2','0.2 mol/kg'],['Cl-','0.5 mol/kg']])
>>> s2.get_salt().formula
'MgCl2'
>>> s2.get_salt().nu_anion
2
>>> s2.get_salt().z_cation
2
"""
try:
salt: Salt = next(d["salt"] for d in self.get_salt_dict().values())
return salt
except StopIteration:
return None
# TODO - modify? deprecate? make a salts property?
[docs]
def get_salt_dict(self, cutoff: float = 1e-6, use_totals: bool = True) -> dict[str, dict[str, float | Salt]]:
"""
Returns a dict that represents the salts of the Solution by pairing anions and cations.
The ``get_salt_dict()`` method examines the ionic composition of a solution and approximates it as a set of
salts instead of individual ions. The method returns a dictionary of Salt objects where the keys are the salt
formulas (e.g., 'NaCl'). The Salt object contains information about the stoichiometry of the salt to
enable its effective concentration to be calculated (e.g., 1 M MgCl2 yields 1 M Mg+2 and 2 M Cl-).
Args:
cutoff: Lowest molal concentration to consider. No salts below this value will be included in the output.
Useful for excluding analysis of trace anions. Defaults to 1e-6 (1 part per million).
use_totals: Whether or not to base the analysis on the concentration of the predominant species of each
element. Note that species in which a given element assumes a different oxidation state are always
treated separately.
Returns:
dict
A dictionary of representing salts in the solution, keyed by the salt formula.
Notes:
The dict maps salt formulas to dictionaries containing their amounts and composition. The amount is stored
in moles under the key "mol", and a :class:`pyEQL.salt_ion_match.Salt` object stored under the "salt" key
represents the composition. Salts are identified by pairing the predominant cations and anions in the
solution, in descending order of their respective equivalent amounts.
Many empirical equations for solution properties such as activity coefficient, partial molar volume, or
viscosity are based on the concentration of single salts (e.g., NaCl). When multiple ions are present
(e.g., a solution containing Na+, Cl-, and Mg+2), it is generally not possible to directly model
these quantities.
Examples:
>>> from pyEQL import Solution
>>> from pyEQL.salt_ion_match import Salt
>>> s1 = Solution(
... solutes={
... 'Na[+1]': '1 mol/L',
... 'Cl[-1]': '1 mol/L',
... 'Ca[+2]': '0.01 mol/kg',
... 'HCO3[-1]': '0.007 mol/kg',
... 'CO3[-2]': '0.001 mol/kg',
... 'ClO[-1]': '0.001 mol/kg',
... }
... )
>>> salt_dict = s1.get_salt_dict()
>>> list(salt_dict) # Only returns salts with concentrations > 1e-3 m
['NaCl', 'Ca(HCO3)2']
>>> salt_dict['NaCl']['salt']
<pyEQL.salt_ion_match.Salt object at ...>
>>> salt_dict['NaCl']['mol']
1.0
>>> salt_dict = s1.get_salt_dict(cutoff=1e-4)
>>> list(salt_dict) # Returns 'Ca(ClO)2' because of reduced cutoff and Cl has different oxidation state
['NaCl', 'Ca(HCO3)2', 'Ca(ClO)2']
>>> salt_dict = s1.get_salt_dict(cutoff=1e-4, use_totals=False)
>>> list(salt_dict) # Returns salts with minor (same oxidation state) species since use_totals=False
['NaCl', 'Ca(HCO3)2', 'CaCO3', 'Ca(ClO)2']
See Also:
:attr:`components`
:attr:`cations`
:attr:`anions`
:class:`pyEQL.salt_ion_match.Salt`
:py:meth:`get_activity_coefficient`
:py:meth:`get_water_activity`
:py:meth:`get_osmotic_coefficient`
"""
salt_dict: dict[str, dict[str, float | Salt]] = {}
if use_totals:
# use only the predominant species for each element
components = {}
for el, lst in self.get_components_by_element(nested=False).items():
component = lst[0]
el_no_oxi_state = el.split("(")[0]
nu_el = self.get_property(component, "pmg_ion").get(el_no_oxi_state, 0)
components[component] = self.get_total_amount(el, "mol").magnitude / nu_el
# add H+ and OH-, which would otherwise be excluded
for k in ["H[+1]", "OH[-1]"]:
if self.components.get(k):
components[k] = self.components[k]
else:
components = self.components
components = dict(sorted(components.items(), key=lambda x: x[1], reverse=True))
# warn if something other than water is the predominant component
if next(iter(components)) != "H2O(aq)":
self.logger.warning("H2O(aq) is not the most prominent component in this Solution!")
# equivalents (charge-weighted moles) of cations and anions
cations = set(self.cations.keys()).intersection(components.keys())
anions = set(self.anions.keys()).intersection(components.keys())
# calculate the charge-weighted (equivalent) concentration of each ion
cation_equiv = {k: self.get_property(k, "charge") * components[k] for k in cations}
anion_equiv = {
k: self.get_property(k, "charge") * components[k] * -1 for k in anions
} # make sure amounts are positive
# sort in descending order of equivalent concentration
cation_equiv = dict(sorted(cation_equiv.items(), key=lambda x: x[1], reverse=True))
anion_equiv = dict(sorted(anion_equiv.items(), key=lambda x: x[1], reverse=True))
len_cat = len(cation_equiv)
len_an = len(anion_equiv)
# start with the first cation and anion
index_cat = 0
index_an = 0
# list(dict) returns a list of [[key, value],]
cation_list = [[k, v] for k, v in cation_equiv.items()]
anion_list = [[k, v] for k, v in anion_equiv.items()]
solvent_mass = self.solvent_mass.to("kg").m
# tolerance for detecting edge cases where equilibrate() slightly changes the
# total amount of a solute
_atol = 1e-16
while index_cat < len_cat and index_an < len_an:
c1 = cation_list[index_cat][-1]
a1 = anion_list[index_an][-1]
salt = Salt(cation_list[index_cat][0], anion_list[index_an][0])
# Use the smaller of the two amounts
equivs_consumed = min(c1, a1)
cation_list[index_cat][-1] -= equivs_consumed
anion_list[index_an][-1] -= equivs_consumed
index_an += 1 if a1 == equivs_consumed else 0
index_cat += 1 if c1 == equivs_consumed else 0
mol = equivs_consumed / (salt.z_cation * salt.nu_cation)
# filter out water and zero, effectively zero, and sub-cutoff salt amounts
if salt.formula != "HOH" and (mol / solvent_mass + _atol) >= cutoff:
salt_dict[salt.formula] = {"salt": salt, "mol": mol}
return dict(sorted(salt_dict.items(), key=lambda x: x[1]["mol"], reverse=True))
[docs]
def equilibrate(
self,
atmosphere: bool = False,
solids: list[str] | None = None,
gases: dict[str, str | float] | None = None,
**kwargs,
) -> None:
"""
This method follows the equilibrate logic used in the NativeEOS engine, adapted as the default behavior for this class.
Adjust the speciation of a Solution object to achieve chemical equilibrium.
Keyword Args:
atmosphere:
Boolean indicating whether to equilibrate the solution w.r.t
atmospheric gases. By default, this considers equilibrium with
atmospheric CO2 (420 ppm) and O2 (0.21 atm). N2 is
typically not considered due to its low solubility and limited
impact on aqueous speciation.
solids:
A list of solids used to achieve liquid-solid equilibrium. Each
solid in this list should be the name of a mineral phase present
in the Phreeqc database (e.g. "Calcite"). We assume a target
saturation index of 0 and an infinite amount of material.
gases:
A dictionary of gases used to achieve liquid-gas equilibrium.
Each key denotes the gas species, and the corresponding value
denotes its concentration, as a log partial pressure value or
other interpretable pressure units. For example, the following
are equivalent (log10(0.000316) = -3.5)
{"CO2": "0.000316 atm"}
{"CO2": -3.5}
**kwargs:
Additional engine-specific options passed to the underlying equilbrium
solver. These may include solver tolerances, or other advanced configuration
parameters introduced in v1.4.0.
"""
if self.engine == "native":
warnings.warn(
'In the next release, the default engine ("native") will '
"transition to a new version of the PHREEQC wrapper for "
"speciation calculations. No change in your script is "
"required, but if you call .equilibrate(), compare results "
"carefully between releases.",
DeprecationWarning,
stacklevel=2,
)
self.engine.equilibrate(self, atmosphere=atmosphere, solids=solids, gases=gases, **kwargs)
# Activity-related methods
[docs]
def get_activity_coefficient(
self,
solute: str,
scale: Literal["molal", "molar", "fugacity", "rational"] = "molal",
) -> Quantity:
"""
Return the activity coefficient of a solute in solution.
The model used to calculate the activity coefficient is determined by the Solution's equation of state
engine.
Args:
solute: The solute for which to retrieve the activity coefficient
scale: The activity coefficient concentration scale
Returns:
Quantity: the activity coefficient as a dimensionless pint Quantity
"""
# 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)
try:
# get the molal-scale activity coefficient from the EOS engine
molal = self.engine.get_activity_coefficient(solution=self, solute=solute)
except (ValueError, ZeroDivisionError):
self.logger.error("Calculation unsuccessful. Returning unit activity coefficient.", exc_info=True)
# return ureg.Quantity(1, "dimensionless")
return ureg.Quantity(1)
# if necessary, convert the activity coefficient to another scale, and return the result
if scale == "molal":
return molal
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.g / ureg.L * total_molality / total_molarity).to(
"dimensionless"
)
if scale == "rational":
return molal * 1 + 0.018015 * ureg.kg / ureg.mol * self.get_total_moles_solute() / self.solvent_mass
raise ValueError("Invalid scale argument. Pass 'molal', 'molar', or 'rational'.")
[docs]
def get_activity(
self,
solute: str,
scale: Literal["molal", "molar", "rational"] = "molal",
) -> Quantity:
"""
Return the thermodynamic activity of the solute in solution on the chosen concentration scale.
Args:
solute:
String representing the name of the solute of interest
scale:
The concentration scale for the returned activity.
Valid options are "molal", "molar", and "rational" (i.e., mole fraction).
By default, the molal scale activity is returned.
Returns:
The thermodynamic activity of the solute in question (dimensionless Quantity)
Notes:
The thermodynamic activity depends on the concentration scale used [rs]_ .
By default, the ionic strength, activity coefficients, and activities are all
calculated based on the molal (mol/kg) concentration scale.
References:
.. [rs] Robinson, R. A.; Stokes, R. H. Electrolyte Solutions: Second Revised
Edition; Butterworths: London, 1968, p.32.
See Also:
:attr:`ionic_strength`
:py:meth:`get_activity_coefficient`
:py:meth:`get_salt`
"""
# switch to the water activity function if the species is H2O
if solute in ["H2O(aq)", "water", "H2O", "HOH"]:
activity = self.get_water_activity()
else:
# determine the concentration units to use based on the desired scale
if scale == "molal":
units = "mol/kg"
elif scale == "molar":
units = "mol/L"
elif scale == "rational":
units = "fraction"
else:
raise ValueError("Invalid scale argument. Pass 'molal', 'molar', or 'rational'.")
activity = (self.get_activity_coefficient(solute, scale=scale) * self.get_amount(solute, units)).magnitude
self.logger.debug(f"Calculated {scale} scale activity of solute {solute} as {activity}")
return ureg.Quantity(activity)
# TODO - engine method
[docs]
def get_osmotic_coefficient(self, scale: Literal["molal", "molar", "rational"] = "molal") -> Quantity:
"""
Return the osmotic coefficient of an aqueous solution.
The method used depends on the Solution object's equation of state engine.
"""
molal_phi = self.engine.get_osmotic_coefficient(self)
if scale == "molal":
return molal_phi
if scale == "rational":
return (
-molal_phi
* 0.018015
* ureg.kg
/ ureg.mol
* self.get_total_moles_solute()
/ self.solvent_mass
/ np.log(self.get_amount(self.solvent, "fraction"))
)
if scale == "fugacity":
return np.exp(
-molal_phi * 0.018015 * ureg.kg / ureg.mol * self.get_total_moles_solute() / self.solvent_mass
- np.log(self.get_amount(self.solvent, "fraction"))
) * ureg.Quantity(1)
raise ValueError("Invalid scale argument. Pass 'molal', 'rational', or 'fugacity'.")
[docs]
def get_water_activity(self) -> Quantity:
r"""
Return the water activity.
Returns:
Quantity:
The thermodynamic activity of water in the solution.
See Also:
:attr:`ionic_strength`
:py:meth:`get_activity_coefficient`
:py:meth:`get_salt`
Notes:
Water activity is related to the osmotic coefficient in a solution containing i solutes by:
.. math:: \ln a_{w} = - \Phi M_{w} \sum_{i} m_{i}
Where :math:`M_{w}` is the molar mass of water (0.018015 kg/mol) and :math:`m_{i}` is the molal
concentration of each species.
If appropriate Pitzer model parameters are not available, the
water activity is assumed equal to the mole fraction of water.
References:
Blandamer, Mike J., Engberts, Jan B. F. N., Gleeson, Peter T., Reis, Joao Carlos R., 2005. "Activity of
water in aqueous systems: A frequently neglected property." *Chemical Society Review* 34, 440-458.
Examples:
>>> s1 = pyEQL.Solution([['Na+','0.3 mol/kg'],['Cl-','0.3 mol/kg']])
>>> s1.get_water_activity()
<Quantity(0.9900944932888518, 'dimensionless')>
"""
osmotic_coefficient = self.get_osmotic_coefficient()
if osmotic_coefficient == 1:
self.logger.warning("Pitzer parameters not found. Water activity set equal to mole fraction")
return self.get_amount("H2O", "fraction")
concentration_sum = np.sum([mol for item, mol in self.components.items() if item != "H2O(aq)"])
concentration_sum /= self.solvent_mass.to("kg").magnitude # converts to mol/kg
self.logger.debug("Calculated water activity using osmotic coefficient")
return ureg.Quantity(np.exp(-osmotic_coefficient * 0.018015 * concentration_sum))
[docs]
def get_chemical_potential_energy(self, activity_correction: bool = True) -> Quantity:
r"""
Return the total chemical potential energy of a solution (not including
pressure or electric effects).
Args:
activity_correction : bool, optional
If True, activities will be used to calculate the true chemical
potential. If False, mole fraction will be used, resulting in
a calculation of the ideal chemical potential.
Returns:
Quantity
The actual or ideal chemical potential energy of the solution, in Joules.
Notes:
The chemical potential energy (related to the Gibbs mixing energy) is
calculated as follows: [koga]_
.. math:: E = R T \sum_i n_i \ln a_i
or
.. math:: E = R T \sum_i n_i \ln x_i
Where :math:`n` is the number of moles of substance, :math:`T` is the temperature in kelvin,
:math:`R` the ideal gas constant, :math:`x` the mole fraction, and :math:`a` the activity of
each component.
Note that dissociated ions must be counted as separate components,
so a simple salt dissolved in water is a three component solution (cation,
anion, and water).
References:
.. [koga] Koga, Yoshikata, 2007. *Solution Thermodynamics and its Application to Aqueous Solutions:*
*A differential approach.* Elsevier, 2007, pp. 23-37.
"""
E = 0 * ureg.J
# loop through all the components and add their potential energy
for item in self.components:
try:
if activity_correction is True:
E += (
ureg.R
* self.temperature.to("K")
* self.get_amount(item, "mol")
* np.log(self.get_activity(item))
)
else:
E += (
ureg.R
* self.temperature.to("K")
* self.get_amount(item, "mol")
* np.log(self.get_amount(item, "fraction"))
)
# If we have a solute with zero concentration, we will get a ValueError
except ValueError:
continue
return E.to("J")
[docs]
def _get_property(self, solute: str, name: str) -> Any | None:
"""Retrieve a thermodynamic property (such as diffusion coefficient)
for solute, and adjust it from the reference conditions to the conditions
of the solution.
Args:
solute: str
String representing the chemical formula of the solute species
name: str
The name of the property needed, e.g.
'diffusion coefficient'
Returns:
Quantity: The desired parameter or None if not found
"""
base_temperature = 25 * ureg.degC
# base_pressure = ureg.Quantity("1 atm")
# query the database using the standardized formula
rform = standardize_formula(solute)
# TODO - add molar volume of water to database?
if name == "size.molar_volume" and rform == "H2O(aq)":
# calculate the partial molar volume for water since it isn't in the database
vol = ureg.Quantity(self.get_property("H2O", "molecular_weight")) / (
self.water_substance.rho * ureg.g / ureg.L
)
return vol.to("cm **3 / mol")
# update the cached data if necessary (only if solute is not already present)
keys = [d["formula"] for d in self._solute_data]
if name not in CORE_PROPERTIES or rform not in keys:
self._solute_data.extend(list(self.database.query({"formula": rform, name: {"$exists": True}})))
data = [d for d in self._solute_data if d["formula"] == rform]
# formulas should always be unique in the database. len==0 indicates no
# data. len>1 indicates duplicate data.
if len(data) == 1:
pass
elif len(data) > 1:
self.logger.warning(f"Duplicate database entries for solute {solute} found!")
elif len(data) == 0:
# try to determine basic properties using pymatgen
doc = Solute.from_formula(rform).as_dict()
data = [doc]
doc: dict = data[0]
try:
# perform temperature-corrections or other adjustments for certain
# parameter types
if name == "transport.diffusion_coefficient":
data = doc["transport"]["diffusion_coefficient"]
if data is not None:
return ureg.Quantity(data["value"]).to("m**2/s")
# just return the base-value molar volume for now; find a way to adjust for concentration later
if name == "size.molar_volume":
data = doc["size"]["molar_volume"]
if data is not None:
base_value = ureg.Quantity(doc["size"]["molar_volume"].get("value"))
if self.temperature != base_temperature:
self.logger.warning(f"Partial molar volume for species {solute} not corrected for temperature")
return base_value
return data
if name == "model_parameters.dielectric_zuber":
return ureg.Quantity(doc["model_parameters"]["dielectric_zuber"]["value"])
if name == "model_parameters.activity_pitzer":
# return a dict
if doc["model_parameters"]["activity_pitzer"].get("Beta0") is not None:
return doc["model_parameters"]["activity_pitzer"]
return None
if name == "model_parameters.molar_volume_pitzer":
# return a dict
if doc["model_parameters"]["molar_volume_pitzer"].get("Beta0") is not None:
return doc["model_parameters"]["molar_volume_pitzer"]
return None
if name == "molecular_weight":
# if the unit is g/mol, we can speed up the Quantity instantiation by about 10x by using built in multiplication instead of parsing the string
if doc.get(name).split(" ")[-1] == "g/mol":
return float(doc.get(name).split(" ")[0]) * ureg("g/mol")
return ureg.Quantity(doc.get(name))
if name == "elements":
return doc.get(name)
if name == "oxi_state_guesses":
# ensure that all oxi states are returned as floats
return {k: float(v) for k, v in doc.get(name).items()}
# for parameters not named above, just return the base value
if name == "pmg_ion" or not isinstance(doc.get(name), dict):
# if the queried value is not a dict, it is a root level key and should be returned as is
return doc.get(name)
val = doc[name].get("value")
# self.logger.warning("%s has not been corrected for solution conditions" % name)
if val is not None:
return ureg.Quantity(val)
except KeyError:
self.logger.error(f"Property {name} for solute {solute} not found in database. Returning None.")
return None
if name == "model_parameters.molar_volume_pitzer":
# return a dict
if doc["model_parameters"]["molar_volume_pitzer"].get("Beta0") is not None:
return doc["model_parameters"]["molar_volume_pitzer"]
return None
if name == "molecular_weight":
return ureg.Quantity(doc.get(name))
if name == "oxi_state_guesses":
return doc.get(name)
# for parameters not named above, just return the base value
if name == "pmg_ion" or not isinstance(doc.get(name), dict):
# if the queried value is not a dict, it is a root level key and should be returned as is
return doc.get(name)
val = doc[name].get("value")
# self.logger.warning("%s has not been corrected for solution conditions" % name)
if val is not None:
return ureg.Quantity(val)
return None
[docs]
def get_transference_number(self, solute: str) -> Quantity:
"""Alias of get_transport_number(). Note that the transference number is only equal to the transport number if there are no concentration or pressure gradients."""
return self.get_transport_number(solute)
[docs]
def get_transport_number(self, solute: str) -> Quantity:
r"""Calculate the transference number of a solute in the solution. Note that this is the
same as the _transport_ number if (and only if) there are no concentration or pressure gradients.
Args:
solute: Formula of the solute for which the transference number is to be calculated.
Returns:
The transference number of `solute`, as a dimensionless Quantity.
Notes:
Transference number is calculated according to :
.. math::
t_i = {D_i z_i^2 C_i \over \sum D_i z_i^2 C_i}
Where :math:`C_i` is the concentration in mol/L, :math:`D_i` is the diffusion
coefficient, and :math:`z_i` is the charge, and the summation extends
over all species in the solution.
Diffusion coefficients :math:`D_i` are adjusted for the effects of temperature
and ionic strength using the method implemented in PHREEQC >= 3.4.
See `get_diffusion_coefficient for` further details.
References:
Bieusheuvel, P.M.; Dykstra, J.E.; *Introduction to Physical Processes in Environmental
Technology*, Section 6.2. https://www.physicsofelectrochemicalprocesses.com/book.pdf.
Geise, G. M.; Cassady, H. J.; Paul, D. R.; Logan, E.; Hickner, M. A. "Specific
ion effects on membrane potential and the permselectivity of ion exchange membranes.""
*Phys. Chem. Chem. Phys.* 2014, 16, 21673-21681.
See Also:
:py:meth:`get_transference_number`
:py:meth:`get_diffusion_coefficient`
:py:meth:`get_molar_conductivity`
"""
solute = standardize_formula(solute)
denominator = numerator = 0
for item, mol in self.components.items():
# the molar conductivity of each species is F/RT D * z^2, and the F/RT factor
# 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
if item == solute:
numerator = term
denominator += term
return ureg.Quantity(numerator / denominator)
[docs]
def _get_molar_conductivity(self, solute: str) -> Quantity:
r"""
Calculate the molar (equivalent) conductivity for a solute.
Args:
solute: String identifying the solute for which the molar conductivity is
to be calculated.
Returns:
The molar or equivalent conductivity of the species in the solution.
Zero if the solute is not charged.
Notes:
Molar conductivity is calculated from the Nernst-Einstein relation [smed]_
.. math::
\lambda_i = \frac{F^2}{RT} D_i z_i^2
Diffusion coefficients :math:`D_i` are adjusted for the effects of temperature
and ionic strength using the method implemented in PHREEQC >= 3.4. See `get_diffusion_coefficient`
for further details.
References:
1. .. [smed] Smedley, Stuart. The Interpretation of Ionic Conductivity in Liquids, pp 1-9. Plenum Press, 1980.
2. https://www.hydrochemistry.eu/exmpls/sc.html
3. Appelo, C.A.J. Solute transport solved with the Nernst-Planck equation for concrete pores with 'free'
water and a double layer. Cement and Concrete Research 101, 2017.
https://dx.doi.org/10.1016/j.cemconres.2017.08.030
4. CRC Handbook of Chemistry and Physics
See Also:
:py:meth:`get_diffusion_coefficient`
"""
D = self.get_diffusion_coefficient(solute)
if D != 0:
molar_cond = (
D * (ureg.e * ureg.N_A) ** 2 * self.get_property(solute, "charge") ** 2 / (ureg.R * self.temperature)
)
else:
# molar_cond = ureg.Quantity(0, "mS / cm / (mol/L)")
molar_cond = 0 * ureg.mS / ureg.cm / (ureg.mol / ureg.L)
self.logger.debug(f"Calculated molar conductivity as {molar_cond} from D = {D!s} at T={self.temperature}")
return molar_cond.to("mS / cm / (mol/L)")
[docs]
def _get_diffusion_coefficient(self, solute: str, activity_correction: bool = True) -> Quantity:
r"""
Get the **temperature-adjusted** diffusion coefficient of a solute.
Args:
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.
Notes:
This method is equivalent to self.get_property(solute, "transport.diffusion_coefficient")
ONLY when the Solution temperature is the same as the reference temperature for the diffusion coefficient
in the database (usually 25 C).
Otherwise, the reference D value is adjusted based on the Solution temperature and (optionally),
ionic strength. The adjustments are
.. math::
D_T = D_{298} \exp(\frac{d}{T} - \frac{d}{298}) \frac{\nu_{298}}{\nu_T}
.. math::
D_{\gamma} = D^0 \exp(\frac{-a1 A |z_i| \sqrt{I}}{1+\kappa a})
.. math::
\kappa a = B \sqrt{I} \frac{a2}{1+I^{0.75}}
where a1, a2, and d are parameters from Ref. 2, A and B are the parameters used in the Debye Huckel
equation, and I is the ionic strength. If the model parameters for a particular solute are not available,
default values of d=0, a1=1.6, and a2=4.73 (as recommended in Ref. 2) are used instead.
References:
1. https://www.hydrochemistry.eu/exmpls/sc.html
2. Appelo, C.A.J. Solute transport solved with the Nernst-Planck equation for concrete pores with 'free'
water and a double layer. Cement and Concrete Research 101, 2017.
https://dx.doi.org/10.1016/j.cemconres.2017.08.030
3. CRC Handbook of Chemistry and Physics
See Also:
pyEQL.activity_correction._debye_parameter_B
pyEQL.activity_correction._debye_parameter_activity
"""
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 "
f"{self.default_diffusion_coeff} m**2/s."
)
D = self.default_diffusion_coeff * ureg.m**2 / ureg.s
# assume reference temperature is 298.15 K (this is the case for all current DB entries)
T_ref = 298.15
mu_ref = 0.0008900225512925807 # water viscosity from IAPWS97 at 298.15 K
T_sol = self.temperature.to("K").magnitude
mu = self.water_substance.mu
# skip temperature correction if within 1 degree
if abs(T_sol - T_ref) > 1 or activity_correction is True:
# get the a1, a2, and d parameters required by the PHREEQC model
try:
doc = self.database.query_one({"formula": rform})
d = doc["model_parameters"]["diffusion_temp_smolyakov"]["d"]["value"]
a1 = doc["model_parameters"]["diffusion_temp_smolyakov"]["a1"]["value"]
a2 = doc["model_parameters"]["diffusion_temp_smolyakov"]["a2"]["value"]
# values will be a str, e.g. "1 dimensionless"
d = float(d.split(" ")[0])
a1 = float(a1.split(" ")[0])
a2 = float(a2.split(" ")[0])
except TypeError:
# this means the database doesn't contain a d value.
# according to Ref 2, the following are recommended default parameters
self.logger.warning(
f"Temperature and ionic strength correction parameters for solute {rform} diffusion "
"coefficient not in database. Using recommended default values of a1=1.6, a2=4.73, and d=0."
)
d = 0
a1 = 1.6
a2 = 4.73
# use the PHREEQC model from Ref 2 to correct for temperature if more than 1 degree different from T_ref
if abs(T_sol - T_ref) > 1:
D *= np.exp(d / T_sol - d / T_ref) * mu_ref / mu
if activity_correction:
A = _debye_parameter_activity(str(self.temperature)).to("kg**0.5/mol**0.5").magnitude / 2.303
B = _debye_parameter_B(str(self.temperature)).to("1/angstrom * kg**0.5/mol**0.5").magnitude
z = self.get_property(solute, "charge")
IS = self.ionic_strength.magnitude
kappaa = B * IS**0.5 * a2 / (1 + IS**0.75)
# correct for ionic strength
D *= np.exp(-a1 * A * abs(z) * IS**0.5 / (1 + kappaa))
# else:
# # per CRC handbook, D increases by 2-3% per degree above 25 C
# return D * (1 + 0.025 * (T_sol - T_ref))
return D
[docs]
def _get_mobility(self, solute: str) -> Quantity:
r"""
Calculate the ionic mobility of the solute.
Args:
solute (str): String identifying the solute for which the mobility is to be calculated.
Returns:
float: The ionic mobility. Zero if the solute is not charged.
Note:
This function uses the Einstein relation to convert a diffusion coefficient into an ionic mobility [smed]_
.. math::
\mu_i = {F |z_i| D_i \over RT}
References:
Smedley, Stuart I. The Interpretation of Ionic Conductivity in Liquids. Plenum Press, 1980.
"""
D = self.get_diffusion_coefficient(solute)
mobility = ureg.N_A * ureg.e * abs(self.get_property(solute, "charge")) * D / (ureg.R * self.temperature)
self.logger.debug(f"Calculated ionic mobility as {mobility} from D = {D!s} at T={self.temperature}")
return mobility.to("m**2/V/s")
[docs]
def get_lattice_distance(self, solute: str) -> Quantity:
r"""
Calculate the average distance between molecules.
Calculate the average distance between molecules of the given solute,
assuming that the molecules are uniformly distributed throughout the
solution.
Args:
solute : str
String representing the name of the solute of interest
Returns:
Quantity: The average distance between solute molecules
Examples:
>>> soln = Solution([['Na+','0.5 mol/kg'],['Cl-','0.5 mol/kg']])
>>> soln.get_lattice_distance('Na+')
1.492964.... nanometer
Notes:
The lattice distance is related to the molar concentration as follows:
.. math:: d = ( C_i N_A ) ^ {-{1 \over 3}}
"""
# calculate the volume per particle as the reciprocal of the molar concentration
# (times avogadro's number). Take the cube root of the volume to get
# the average distance between molecules
distance = (self.get_amount(solute, "mol/L") * ureg.N_A) ** (-1 / 3)
return distance.to("nm")
def _adjust_charge_balance(self, atol=1e-8) -> None:
"""Helper method to adjust the charge balance of the Solution."""
cb = self.charge_balance
if not np.isclose(cb, 0, atol=atol):
self.logger.info(f"Solution is not electroneutral (C.B. = {cb} eq/L).")
if self.balance_charge is None:
# Nothing to do.
self.logger.info("balance_charge is None, so no charge balancing will be performed.")
return
self.logger.info(f"Adjusting {self._cb_species} to compensate.")
if self.balance_charge == "pH":
# the charge imbalance associated with the H+ / OH- system can be expressed
# as ([H+] - [OH-]) or ([H+] - K_W/[H+]). If we adjust H+, we also have to
# adjust OH- to maintain water equilibrium.
C_hplus = self.get_amount("H+", "mol/L").magnitude
start_imbalance = C_hplus - K_W / C_hplus
new_imbalance = start_imbalance - cb
# calculate the new concentration of H+ that will balance the charge
# solve H^2 - new_imbalance H - K_W = 0, so a=1, b=-new_imbalance, c=-K_W. Note that this is guaranteed to have real roots as
# b^2-4ac > 0
new_hplus = max(
[
(new_imbalance + np.sqrt(new_imbalance**2 + 4 * 1 * K_W)) / 2,
(new_imbalance - np.sqrt(new_imbalance**2 + 4 * 1 * K_W)) / 2,
]
)
self.set_amount("H+", f"{new_hplus} mol/L")
self.set_amount("OH-", f"{K_W / new_hplus} mol/L")
return
z = self.get_property(self._cb_species, "charge")
try:
self.add_amount(self._cb_species, f"{-1 * cb / z} mol")
return
except ValueError:
# if the concentration is negative, it must mean there is not enough present.
# remove everything that's present and log an error.
self.components[self._cb_species] = 0.0
self.logger.error(
f"There is not enough {self._cb_species} present to balance the charge. Try a different species."
)
return
def _check_water_stability(self, tol=1e-6) -> None:
"""Helper method to adjust the thermodynamic stability of the Solution."""
temp = self.temperature.to("K")
E0_O2 = 1.229 * ureg.V
lower_limit = -float(self.pH)
upper_limit = (ureg.faraday_constant * E0_O2 / (2.303 * ureg.R * temp)).to_base_units().magnitude - float(
self.pH
)
if self.pE < lower_limit - tol:
msg = (
f"Water is thermodynamically unstable under reducing conditions: "
f"pH={self.pH:.2f}, pE={self.pE:.2f}. "
f"Hydrogen evolution may occur."
)
self.logger.warning(msg)
elif self.pE > upper_limit + tol:
msg = (
f"Water is thermodynamically unstable under oxidizing conditions: "
f"pH={self.pH:.2f}, pE={self.pE:.2f}. "
f"Oxygen evolution may occur."
)
self.logger.warning(msg)
def _update_volume(self):
"""Recalculate the solution volume based on composition."""
self._volume = self._get_solvent_volume() + self._get_solute_volume()
def _get_solvent_volume(self):
"""Return the volume of the pure solvent."""
# calculate the volume of the pure solvent
solvent_vol = self.solvent_mass / (self.water_substance.rho * ureg.g / ureg.L)
return solvent_vol.to("L")
def _get_solute_volume(self):
"""Return the volume of only the solutes."""
return self.engine.get_solute_volume(self)
[docs]
def as_dict(self) -> dict:
"""Convert the Solution into a dict representation that can be serialized to .json or other format."""
# clear the volume update flag, if required
if self.volume_update_required:
self._update_volume()
d = super().as_dict()
for k, v in d.items():
# convert all Quantity to str
if isinstance(v, Quantity):
d[k] = str(v)
# replace solutes with the current composition
d["solutes"] = {k: f"{v} mol" for k, v in self.components.items()}
# replace the engine with the associated str
d["engine"] = self._engine
# d["logger"] = self.logger.__dict__
return d
[docs]
@classmethod
def from_dict(cls, d: dict, **kwargs) -> Solution:
"""Instantiate a Solution from a dictionary generated by as_dict().
Args:
d (dict): A dictionary representation of a Solution, generated by the as_dict() method.
Kwargs:
Any kwargs passed to this method will be passed to the Solution __init__ method, and will override
any values in the dict.
"""
# because of the automatic volume updating that takes place during the __init__ process,
# care must be taken here to recover the exact quantities of solute and volume
# first we store the volume and solvent mass of the serialized solution
orig_volume = ureg.Quantity(d["volume"])
orig_moles = ureg.Quantity(d["solutes"]["H2O(aq)"]).to("mol").magnitude if "H2O(aq)" in d["solutes"] else None
# then instantiate a new one
decoded = {k: MontyDecoder().process_decoded(v) for k, v in d.items() if not k.startswith("@")}
decoded.update(kwargs)
new_sol = cls(**decoded)
# now determine how different the new solution volume is from the original
scale_factor = (orig_volume / new_sol.volume).magnitude
# reset the new solution volume to that of the original. In the process of
# doing this, all the solute amounts are scaled by new_sol.volume / volume
new_sol.volume = str(orig_volume)
# undo the scaling by diving by that scale factor
for sol in new_sol.components:
new_sol.components[sol] /= scale_factor
if orig_moles is not None:
new_sol.components["H2O(aq)"] = orig_moles
# ensure that another volume update won't be triggered by these changes
# (this line should in principle be unnecessary, but it doesn't hurt anything)
new_sol.volume_update_required = False
return new_sol
[docs]
@classmethod
def from_preset(
cls,
preset: Literal[
"ash",
"batt_mfg",
"batt_recycling",
"coal_washing",
"CRL",
"drilling",
"excavation",
"FGD",
"flotation",
"gasification",
"geothermal",
"leachate",
"mine_drainage",
"mine_tailings",
"normal saline",
"plating",
"pw_conv",
"pw_unconv",
"rainwater",
"refining",
"Ringers lactate",
"seawater",
"semiconductor",
"smelting",
"tanning",
"urine",
"waste_gas",
"wastewater",
],
**kwargs,
) -> Solution:
r"""Instantiate a solution from a preset composition.
Args:
preset (str): String representing the desired solution.
Valid entries are 'seawater', 'rainwater', 'wastewater',
'urine', 'normal saline' and 'Ringers lactate'.
Kwargs:
Any kwargs passed to this method will be passed to the Solution __init__ method, and
will override any values in the preset file. This allows you to use a preset as a starting
point and then modify it as needed by, e.g., changing the modeling engine or database.
Returns:
A pyEQL Solution object.
Raises:
FileNotFoundError: If the given preset file doesn't exist on the file system.
Notes:
The following sections explain the different solution options:
- 'ash' - bottom ash transport wastewater from fossil fuel combustion [kwptr2026]_
- 'batt_mfg' - wastewater from lead-acid and legacy battery manufacturing processes [kwptr2026]_
- 'batt_recycling' - wastewater from lithium ion battery recycling operations [kwptr2026]_
- 'coal_washing' - wastewater generated from coal preparation and washing [kwptr2026]_
- 'CRL' - combustion residual leachate (CRL) from fossil fuel combustion landfills [kwptr2026]_
- 'drilling' - oil and gas drilling wastewater from drilling fluids and cuttings [kwptr2026]_
- 'excavation' - wastewater from excavation in metal ore mining operations [kwptr2026]_
- 'FGD' - flue gas desulfurization wastewater from SO2 removal from fossil fuel combustion [kwptr2026]_
- 'flotation' - milling and flotation wastewater from metal ore extraction [kwptr2026]_
- 'gasification' - wastewater from gasification of carbon-based feedstocks to syngas [kwptr2026]_
- 'geothermal' - geothermal produced water from geothermal power generation [kwptr2026]_
- 'leachate' - leachate from metal ore mining wastes [kwptr2026]_
- 'mine_drainage' - acid mine drainage wastewater from coal and metal ore mining operations [kwptr2026]_
- 'mine_tailings' - mine tailings pond water from collective metal ore and mining wastes [kwptr2026]_
- 'normal saline' or 'NS' - normal saline solution used in medicine [saline]_
- 'plating' - typical wastewater from metal electroplating operations [kwptr2026]_
- 'pw_conv' - produced water from conventional hydrocarbon production [kwptr2026]_
- 'pw_unconv' - produced water from unconventional hydrocarbon production [kwptr2026]_
- 'rainwater' - pure water in equilibrium with atmospheric CO2 at pH 6
- 'refining' - petroleum refining wastewater from crude oil refineries [kwptr2026]_
- 'Ringers lacatate' or 'RL' - Ringer's lactate solution used in medicine [lactate]_
- 'seawater' or 'SW'- Standard Seawater. See Table 4 of the Reference for Composition [mf08]_
- 'semiconductor' - semiconductor and electronics manufacturing wastewater [kwptr2026]_
- 'smelting' - wastewater from metal ore smelting and refining from pyrometallurgical slags [kwptr2026]_
- 'tanning' - wastewater from leather tanning and finishing operations [kwptr2026]_
- 'urine' - typical human urine. See Table 3-15 of [me13]_
- 'waste_gas' - wastewater from waste gas treatment during pyrometallurgical processing in metal ore and mining operations [kwptr2026]_
- 'wastewater' or 'WW' - medium strength domestic wastewater. See Table 3-18 of [me13]_
References:
.. [mf08] Millero, Frank J. "The composition of Standard Seawater and the definition of
the Reference-Composition Salinity Scale." *Deep-sea Research. Part I* 55(1), 2008, 50-72.
.. [me13] Metcalf & Eddy, Inc. et al. *Wastewater Engineering: Treatment and Resource Recovery*, 5th Ed.
McGraw-Hill, 2013.
.. [saline] https://en.wikipedia.org/w/index.php?title=Saline_(medicine)&oldid=1298292693
.. [lactate] https://en.wikipedia.org/wiki/Ringer%27s_lactate_solution
.. [kwptr2026] Ryan S. Kingsbury, Monong Wang, Jaebeom Park et al. Composition and Critical Mineral Content of Major Industrial Wastewaters: Implications for Treatment and Resource Recovery Technologies, 05 February 2026, PREPRINT (Version 2) available at Research Square [https://www.researchsquare.com/article/rs-8743330/v2]
"""
# preset_dir = files("pyEQL") / "presets"
# Path to the YAML and JSON files corresponding to the preset
yaml_path = files("pyEQL") / "presets" / f"{preset}.yaml"
json_path = files("pyEQL") / "presets" / f"{preset}.json"
# Check if the file exists
if yaml_path.exists():
preset_path = yaml_path
elif json_path.exists():
preset_path = json_path
else:
raise FileNotFoundError(f"Invalid preset! File '{yaml_path}' or '{json_path} not found!")
# Create and return a Solution object
return cls().from_file(preset_path, **kwargs)
[docs]
def to_file(self, filename: str | Path) -> None:
"""Saving to a .yaml or .json file.
Args:
filename (str | Path): The path to the file to save Solution.
Valid extensions are .json or .yaml.
"""
str_filename = str(filename)
if not ("yaml" in str_filename.lower() or "json" in str_filename.lower()):
self.logger.error("Invalid file extension entered - {str_filename}")
raise ValueError("File extension must be .json or .yaml")
if "yaml" in str_filename.lower():
solution_dict = self.as_dict()
solution_dict.pop("database")
dumpfn(solution_dict, filename)
else:
dumpfn(self, filename)
[docs]
@classmethod
def from_file(cls, filename: str | Path, **kwargs) -> Solution:
"""Loading from a .yaml or .json file.
Args:
filename (str | Path): Path to the .json or .yaml file (including extension) to load the Solution from.
Valid extensions are .json or .yaml.
Kwargs:
Any kwargs passed to this method will be passed to the Solution __init__ method, and will override any values in the file. This allows you to use a file as a starting point and then modify it as needed by, e.g., changing the modeling engine or database.
NOTE: CURRENTLY ONLY SUPPORTED FOR YAML FILES!
Returns:
A pyEQL Solution object.
Raises:
FileNotFoundError: If the given filename doesn't exist on the file system.
"""
if not os.path.exists(filename):
raise FileNotFoundError(f"File '{filename}' not found!")
str_filename = str(filename)
if "yaml" in str_filename.lower():
true_keys = [
"solutes",
"volume",
"temperature",
"pressure",
"pH",
"pE",
"balance_charge",
"solvent",
"engine",
# "database",
]
solution_dict = loadfn(filename)
keys_to_delete = [key for key in solution_dict if key not in true_keys]
for key in keys_to_delete:
solution_dict.pop(key)
for k, v in kwargs.items():
solution_dict[k] = v
return cls.from_dict(solution_dict)
return loadfn(filename)
# arithmetic operations
def __add__(self, other: Solution) -> Solution:
"""
Solution addition: mix two solutions together.
Args:
other: The Solutions to be mixed with this solution.
Returns:
A Solution object that represents the result of mixing this solution and other.
Notes:
The initial volume of the mixed solution is set as the sum of the volumes of this solution and other.
The pressure and temperature are volume-weighted averages. The pH and pE values are currently APPROXIMATE
because they are calculated assuming H+ and e- mix conservatively (i.e., the mixing process does not
incorporate any equilibration reactions or buffering). Such support is planned in a future release.
"""
# check to see if the two solutions have the same solvent
if self.solvent != other.solvent:
raise ValueError("Cannot add Solution with different solvents!")
if self._engine != other._engine:
raise ValueError("Cannot add Solution with different engines!")
if self.database != other.database:
raise ValueError("Cannot add Solution with different databases!")
# set the pressure for the new solution
p1 = self.pressure
t1 = self.temperature
v1 = self.volume
p2 = other.pressure
t2 = other.temperature
v2 = other.volume
# set the initial volume as the sum of the volumes
mix_vol = v1 + v2
# check to see if the solutions have the same temperature and pressure
if p1 != p2:
self.logger.info(
"Adding two solutions of different pressure. Pressures will be averaged (weighted by volume)"
)
mix_pressure = (p1 * v1 + p2 * v2) / (mix_vol)
if t1 != t2:
self.logger.info(
"Adding two solutions of different temperature. Temperatures will be averaged (weighted by volume)"
)
# do all temperature conversions in Kelvin to avoid ambiguity associated with "offset units". See pint docs.
mix_temperature = (t1.to("K") * v1 + t2.to("K") * v2) / (mix_vol)
# retrieve the amount of each component in the parent solution and
# store in a list.
mix_amounts = FormulaDict({})
for sol, amt in [*self.components.items(), *other.components.items()]:
mix_amounts[sol] = amt + mix_amounts.get(sol, 0.0)
# TODO - call equilibrate() here once the method is functional to get new pH and pE, instead of the below
warnings.warn(
"The pH and pE value of the mixed solution is approximate! More accurate addition (mixing) of"
"this property is planned for a future release."
)
# calculate the new pH and pE (before reactions) by mixing
# for pH, we make sure to conserve the mass of H+ and OH-. By not passing
# a kwarg for pH (i.e., by using the default value), the H+ concentration
# will override and determine the pH value of the mixed solution.
# pE = -log[e-], so calculate the moles of e- in each solution and mix them
mol_e_self = 10 ** (-1 * self.pE) * self.volume.to("L").magnitude
mol_e_other = 10 ** (-1 * other.pE) * other.volume.to("L").magnitude
mix_pE = -np.log10((mol_e_self + mol_e_other) / mix_vol.to("L").magnitude)
solutes = {sol: f"{amount} mol" for sol, amount in mix_amounts.items()}
# create a new solution
return Solution(
solutes=solutes,
volume=str(mix_vol),
pressure=str(mix_pressure),
temperature=str(mix_temperature.to("K")),
# pH=7, # leave at default value so that H+ concentration determines pH
pE=mix_pE,
engine=self._engine,
solvent=self.solvent,
database=self.database,
)
def __sub__(self, other: Solution) -> None:
raise NotImplementedError("Subtraction of solutions is not implemented.")
def __mul__(self, factor: float) -> None:
"""
Solution multiplication: scale all components by a factor. For example, Solution * 2 will double the moles of
every component (including solvent). No other properties will change.
"""
self.volume *= factor
return self
def __truediv__(self, factor: float) -> None:
"""
Solution division: scale all components by a factor. For example, Solution / 2 will remove half of the moles
of every compoonents (including solvent). No other properties will change.
"""
self.volume /= factor
return self
# informational methods
[docs]
def print(
self,
mode: Literal["all", "ions", "cations", "anions", "neutrals"] = "all",
units: Literal["ppm", "mol", "mol/kg", "mol/L", "%", "activity"] = "mol",
places=4,
):
"""
Print details about the Solution.
Args:
mode: Whether to list the amounts of all solutes, or only anions, cations, any ion, or any neutral solute.
units: The units to list solute amounts in. "activity" will list dimensionless activities instead of
concentrations.
places: The number of decimal places to round the solute amounts.
"""
print(self)
str1 = "Activities" if units == "activity" else "Concentrations"
str2 = f" ({units})" if units != "activity" else ""
header = f"\nComponent {str1}{str2}:"
print(header)
print("=" * (len(header) - 1))
for i in self.components:
if mode != "all":
z = self.get_property(i, "charge")
if (
(z != 0 and mode == "neutrals")
or (z >= 0 and mode == "anions")
or (z <= 0 and mode == "cations")
or (z == 0 and mode == "ions")
):
continue
amt = self.get_activity(i).magnitude if units == "activity" else self.get_amount(i, units).magnitude
print(f"{i:<12} {amt:0.{places}f}")
def __str__(self) -> str:
# set output of the print() statement for the solution
l1 = f"Volume: {self.volume:.3f~}"
l2 = f"Temperature: {self.temperature:.3f~}"
l3 = f"Pressure: {self.pressure:.3f~}"
l4 = f"pH: {self.pH:.1f}"
l5 = f"pE: {self.pE:.1f}"
l6 = f"Solvent: {self.solvent}"
l7 = f"Components: {self.components.keys():}"
return f"{l1}\n{l2}\n{l3}\n{l4}\n{l5}\n{l6}\n{l7}"
"""
Legacy methods to be deprecated in a future release.
"""
[docs]
@deprecated(message="add_solute() is deprecated. Use add_amount() instead.")
def add_solvent(self, formula: str, amount: str): # pragma: no cover
"""Same as add_solute but omits the need to pass solvent mass to pint."""
quantity = ureg.Quantity(amount)
mw = self.get_property(formula, "molecular_weight")
target_mol = quantity.to("moles", "chem", mw=mw, volume=self.volume, solvent_mass=self.solvent_mass)
self.components[formula] = target_mol.to("moles").magnitude
[docs]
def get_saturation_index(self, get_plot=None) -> dict:
r"""
Calculate the saturation index of a solute in the solution.
Notes:
The saturation index (:math:`\mathrm{SI}`) is defined as log10(IAP/Ksp), where IAP is the ion activity product and Ksp is the solubility product constant.
This method calculates the saturation index based on the active engine and database from `__init__`. The interpretation of the saturation index values is as follows:
- :math:`\mathrm{SI} < 0`: The solution is **undersaturated**. The solid tends to dissolve if present.
- :math:`\mathrm{SI} = 0`: The solution is **at saturation equilibrium**. Therefore, at the saturation limit, the SI is zero.
- :math:`\mathrm{SI} > 0`: The solution is **supersaturated**. Precipitation is thermodynamically favored, although kinetic factors may delay or prevent it.
Args:
get_plot (bool, optional):
If True, displays an interactive bar plot of saturation indices sorted from most oversaturated to least. Defaults to None (no plot).
Returns:
dict:
A dictionary with mineral phase names as keys and their saturation index values as values, sorted in descending order (most oversaturated to least oversaturated).
"""
engine = self.engine
if not hasattr(engine, "ppsol"):
raise NotImplementedError(f"Engine {type(engine).__name__} does not support saturation index calculations.")
# caching method from Phrqsol
if (engine.ppsol is None) or (self.components != engine._stored_comp):
engine._destroy_ppsol()
engine._setup_ppsol(self)
ppsol = engine.ppsol
phases = list(ppsol.phases.keys())
eq_species_dict = {phase: ppsol.si(phase) for phase in phases}
sorted_eq_species_dict = dict(sorted(eq_species_dict.items(), key=lambda item: item[1], reverse=True))
if get_plot:
import pandas as pd # noqa: PLC0415
import plotly.express as px # noqa: PLC0415
df = pd.DataFrame( # noqa: PD901
{"species": list(sorted_eq_species_dict.keys()), "si": list(sorted_eq_species_dict.values())}
)
fig = px.bar(
df,
x="species",
y="si",
labels={"species": "Mineral Phase", "si": "Saturation Index"},
color="si",
color_continuous_scale="Mint",
)
fig.update_layout(
xaxis_tickangle=-45,
template="plotly_white",
height=500,
)
fig.show()
return sorted_eq_species_dict