diff --git a/CHANGELOG.md b/CHANGELOG.md index d8c57e45..3bfe65cc 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,8 +2,12 @@ ## 0.6.0 (in progress) -- `Solution`: new properties `pressure`, `temperature`, `pE`, `pH`, `mass`, `density`, - `viscosity_dynamic`, `viscosity_kinematic` +- **DEPRECATION NOTICE** `Solution`: new properties `pressure`, `temperature`, `pE`, +- `pH`, `mass`, `density`, `viscosity_dynamic`, `viscosity_kinematic`, `ionic_strength`, +- `conductivity`, `debye_length`, `bjerrum_length`, `alkalinity`, `hardness`, +- `dielectric_constant`, `charge_balance` have replaced the corresponding get_XXX and + set_XXX (for temperature and pressure) methods, which will be removed in a future + release. `get_viscosity_relative` will be removed entirely. - `Solution`: add support for passing solutes as a `dict` - Implement extensible system for connecting `Solution` to various activity and speciation models. Models can be integrated into pyEQL by implementing an `EOS` class. The desired @@ -12,8 +16,6 @@ parameters are not available) or `ideal` for a dummy engine that returns unit activity coefficients. Support for additional external engines such as [`phreeqpython`](https://github.com/Vitens/phreeqpython) is planned. - **BREAKING CHANGE** disable 'verbose' kwarg in `get_activity` and `get_activity_coefficient` -- **DEPRECATION NOTICE** - `get_temperature()`, `set_temperature()`, `get_pressure`, `set_pressure()`, - `get_mass()`, `get_viscosity_dynamic`, `get_viscosity_kinematic`, and `get_density` will be removed in the next release. Use direct access via property (e.g. `Solution.pressure`) instead. `get_viscosity_relative` will be removed entirely. - Add more comprehensive platform testing via `tox` - Replace `water_properties.py` with [iapws](https://github.com/jjgomera/iapws) package - Replace elements.py with `pymatgen.core.periodic_table` diff --git a/src/pyEQL/solution.py b/src/pyEQL/solution.py index d04856de..02320e33 100644 --- a/src/pyEQL/solution.py +++ b/src/pyEQL/solution.py @@ -426,7 +426,8 @@ def density(self) -> Quantity: """ return self.mass / self.get_volume() - def get_dielectric_constant(self): + @property + def dielectric_constant(self) -> Quantity: """ Returns the dielectric constant of the solution. @@ -675,952 +676,812 @@ def conductivity(self): return EC.to("S/m") - def get_osmotic_pressure(self): + @property + def ionic_strength(self) -> Quantity: """ - Return the osmotic pressure of the solution relative to pure water. + Return the ionic strength of the solution. - Parameters - ---------- - None + 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 osmotic pressure of the solution relative to pure water in Pa + Quantity : + The ionic strength of the parent solution, mol/kg. See Also -------- + get_activity get_water_activity - get_osmotic_coefficient - get_salt Notes ----- - Osmotic pressure is calculated based on the water activity [#]_ [#]_ : - - .. math:: \\Pi = {RT \\over 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. - + The ionic strength is calculated according to: - References - ---------- - .. [#] Sata, Toshikatsu. Ion Exchange Membranes: Preparation, Characterization, and Modification. Royal Society of Chemistry, 2004, p. 10. + .. math:: I = \\sum_i m_i z_i^2 - .. [#] http://en.wikipedia.org/wiki/Osmotic_pressure#Derivation_of_osmotic_pressure + Where :math:`m_i` is the molal concentration and :math:`z_i` is the charge on species i. Examples -------- - >>> s1=pyEQL.Solution() - >>> s1.get_osmotic_pressure() - 0.0 - >>> s1 = pyEQL.Solution([['Na+','0.2 mol/kg'],['Cl-','0.2 mol/kg']]) - >>> soln.get_osmotic_pressure() - - """ - # TODO - tie this into parameter() and solvent() objects - partial_molar_volume_water = 1.82e-5 * unit("m ** 3/mol") + >>> s1.ionic_strength + - osmotic_pressure = ( - -1 - * unit.R - * self.temperature - / partial_molar_volume_water - * math.log(self.get_water_activity()) - ) - logger.info( - "Computed osmotic pressure of solution as %s Pa at T= %s degrees C" - % (osmotic_pressure, self.temperature) - ) - return osmotic_pressure.to("Pa") + >>> 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 + + """ + ionic_strength = 0 + for solute in self.components.keys(): + ionic_strength += ( + 0.5 + * self.get_amount(solute, "mol/kg") + * self.components[solute].get_formal_charge() ** 2 + ) - # Concentration Methods + return ionic_strength - def p(self, solute, activity=True): + @property + def charge_balance(self) -> float: """ - Return the negative log of the activity of solute. - - Generally used for expressing concentration of hydrogen ions (pH) + Return the charge balance of the solution. - Parameters - ---------- - 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. + Return the charge balance of the solution. The charge balance represents the net electric charge + on the solution and SHOULD equal zero at all times, but due to numerical errors will usually + have a small nonzero value. Returns ------- - Quantity - The negative log10 of the activity (or molar concentration if - activity = False) of the solute. + float : + The charge balance of the solution, in equivalents. - Examples - -------- - TODO + Notes + ----- + The charge balance is calculated according to: - """ - try: - if activity is True: - return -1 * math.log10(self.get_activity(solute)) - elif activity is False: - return -1 * math.log10(self.get_amount(solute, "mol/L").magnitude) - # if the solute has zero concentration, the log will generate a ValueError - except ValueError: - return 0 + .. math:: CB = F \\sum_i n_i z_i + + Where :math:`n_i` is the number of moles, :math:`z_i` is the charge on species i, and :math:`F` is the Faraday constant. - def get_amount(self, solute, units): """ - Return the amount of 'solute' in the parent solution. + charge_balance = 0 + F = (unit.e * unit.N_A).magnitude + for solute in self.components.keys(): + charge_balance += ( + self.get_amount(solute, "mol").magnitude + * self.components[solute].get_formal_charge() + * F + ) - The amount of a solute can be given in a variety of unit types. - 1. substance per volume (e.g., 'mol/L') - 1. substance per mass of solvent (e.g., 'mol/kg') - 1. mass of substance (e.g., 'kg') - 1. moles of substance ('mol') - 1. mole fraction ('fraction') - 1. percent by weight (%) + return charge_balance - Parameters - ---------- - 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 + # TODO - need tests for alkalinity + @property + def alkalinity(self): + """ + Return the alkalinity or acid neutralizing capacity of a solution Returns ------- - The amount of the solute in question, in the specified units + Quantity : + The alkalinity of the solution in mg/L as CaCO3 + Notes + ----- + The alkalinity is calculated according to: [#]_ - See Also - -------- - add_amount - set_amount - get_total_amount - get_osmolarity - get_osmolality - get_solvent_mass - get_mass - get_total_moles_solute + .. math:: Alk = F \\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 charge. + In this method, the set of conservative cations is all Group I and Group II cations, and the conservative anions are all the anions of strong acids. + + References + ---------- + .. [#] Stumm, Werner and Morgan, James J. Aquatic Chemistry, 3rd ed, + pp 165. Wiley Interscience, 1996. """ - # retrieve the number of moles of solute and its molecular weight - try: - moles = self.get_solute(solute).get_moles() - mw = self.get_solute(solute).get_molecular_weight() - # 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 0 * unit(units) - except DimensionalityError: - logger.error("Unsupported unit specified for get_amount") - return 0 + alkalinity = 0 * unit("mol/L") + equiv_wt_CaCO3 = 100.09 / 2 * unit("g/mol") - # 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 unecessary - # function calls. - if units == "fraction": - return moles / (self.get_moles_solvent() + self.get_total_moles_solute()) - elif units == "%": - return moles.to("kg", "chem", mw=mw) / self.mass.to("kg") * 100 - elif unit(units).dimensionality in ( - "[substance]/[length]**3", - "[mass]/[length]**3", - ): - return moles.to(units, "chem", mw=mw, volume=self.get_volume()) - elif unit(units).dimensionality in ("[substance]/[mass]", "[mass]/[mass]"): - return moles.to(units, "chem", mw=mw, solvent_mass=self.get_solvent_mass()) - elif unit(units).dimensionality == "[mass]": - return moles.to(units, "chem", mw=mw) - elif unit(units).dimensionality == "[substance]": - return moles.to(units) - else: - logger.error("Unsupported unit specified for get_amount") - return None + base_cations = [ + "Li+", + "Na+", + "K+", + "Rb+", + "Cs+", + "Fr+", + "Be+2", + "Mg+2", + "Ca+2", + "Sr+2", + "Ba+2", + "Ra+2", + ] + acid_anions = ["Cl-", "Br-", "I-", "SO4-2", "NO3-", "ClO4-", "ClO3-"] - def get_total_amount(self, element, units): + for item in self.components: + if item in base_cations: + z = self.get_solute(item).get_formal_charge() + alkalinity += self.get_amount(item, "mol/L") * z + if item in acid_anions: + z = self.get_solute(item).get_formal_charge() + alkalinity -= self.get_amount(item, "mol/L") * z + + # convert the alkalinity to mg/L as CaCO3 + return (alkalinity * equiv_wt_CaCO3).to("mg/L") + + @property + def hardness(self): """ - Return the total amount of 'element' (across all solutes) in the solution. + 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. Parameters ---------- - element : str - String representing the name of the element of interest - units : str - Units desired for the output. Examples of valid units are - 'mol/L','mol/kg','mol', 'kg', and 'g/L' + None Returns ------- - The total amount of the element in the solution, in the specified units - - Notes - ----- - There is currently no way to distinguish between different oxidation - states of the same element (e.g. TOTFe(II) vs. TOTFe(III)). This - is planned for a future release. (TODO) + Quantity + The hardness of the solution in mg/L as CaCO3 - See Also - -------- - get_amount """ - import pyEQL.chemical_formula as ch - - TOT = 0 * unit(units) + hardness = 0 * unit("mol/L") + equiv_wt_CaCO3 = 100.09 / 2 * unit("g/mol") - # loop through all the solutes, process each one containing element for item in self.components: - # check whether the solute contains the element - if ch.contains(item, element): - # start with the amount of the solute in the desired units - amt = self.get_amount(item, units) + z = self.get_solute(item).get_formal_charge() + if z > 1: + hardness += z * self.get_amount(item, "mol/L") - # convert the solute amount into the amount of element by - # either the mole / mole or weight ratio - if unit(units).dimensionality in ( - "[substance]", - "[substance]/[length]**3", - "[substance]/[mass]", - ): - TOT += amt * ch.get_element_mole_ratio(item, element) + # convert the hardness to mg/L as CaCO3 + return (hardness * equiv_wt_CaCO3).to("mg/L") - elif unit(units).dimensionality in ( - "[mass]", - "[mass]/[length]**3", - "[mass]/[mass]", - ): - TOT += amt * ch.get_element_weight_fraction(item, element) + @property + def debye_length(self) -> Quantity: + """ + Return the Debye length of a solution - return TOT + Debye length is calculated as [#]_ - def add_amount(self, solute, amount): - """ - Add the amount of 'solute' to the parent solution. + .. 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. Parameters ---------- - 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 + None Returns ------- - Nothing. The concentration of solute is modified. + Quantity + The Debye length, in nanometers. + References + ---------- + .. [#] https://en.wikipedia.org/wiki/Debye_length#Debye_length_in_an_electrolyte See Also -------- - Solute.add_moles() + ionic_strength + dielectric_constant + """ + # to preserve dimensionality, convert the ionic strength into mol/L units + ionic_strength = self.ionic_strength.magnitude * unit("mol/L") + dielectric_constant = self.dielectric_constant - # 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 unit(amount).dimensionality in ( - "[substance]/[length]**3", - "[mass]/[length]**3", - ): + debye_length = ( + dielectric_constant + * unit.epsilon_0 + * unit.k + * self.temperature + / (2 * unit.N_A * unit.e**2 * ionic_strength) + ) ** 0.5 - # store the original volume for later - orig_volume = self.get_volume() + return debye_length.to("nm") - # change the amount of the solute present to match the desired amount - self.get_solute(solute).add_moles( - amount, self.get_volume(), self.get_solvent_mass() - ) + @property + def bjerrum_length(self) -> Quantity: + """ + Return the Bjerrum length of a solution - # set the amount to zero and log a warning if the desired amount - # change would result in a negative concentration - if self.get_amount(solute, "mol").magnitude < 0: - logger.warning( - "Attempted to set a negative concentration for solute %s. Concentration set to 0" - % solute - ) - self.set_amount(solute, "0 mol") + Bjerrum length representes the distance at which electrostatic + interactions between particles become comparable in magnitude + to the thermal energy.:math:`\\lambda_B` is calculated as [#]_ - # calculate the volume occupied by all the solutes - solute_vol = self._get_solute_volume() + .. math:: - # determine the volume of solvent that will preserve the original volume - target_vol = orig_volume - solute_vol + \\lambda_B = {e^2 \\over (4 \\pi \\epsilon_r \\epsilon_o k_B T)} - # adjust the amount of solvent - # volume in L, density in kg/m3 = g/L - target_mass = ( - target_vol.magnitude * self.water_substance.rho * unit.Quantity("1 g") - ) + 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. - mw = self.get_solvent().get_molecular_weight() - target_mol = target_mass / mw - self.get_solvent().moles = target_mol + Parameters + ---------- + None - else: + Returns + ------- + Quantity + The Bjerrum length, in nanometers. - # change the amount of the solute present - self.get_solute(solute).add_moles( - amount, self.get_volume(), self.get_solvent_mass() - ) + References + ---------- + .. [#] https://en.wikipedia.org/wiki/Bjerrum_length - # set the amount to zero and log a warning if the desired amount - # change would result in a negative concentration - if self.get_amount(solute, "mol").magnitude < 0: - logger.warning( - "Attempted to set a negative concentration for solute %s. Concentration set to 0" - % solute - ) - self.set_amount(solute, "0 mol") + Examples + -------- + >>> s1 = pyEQL.Solution() + >>> s1.bjerrum_length + - # 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.get_solvent_mass() <= unit("0 kg"): - logger.error("All solvent has been depleted from the solution") - return None - else: - # set the volume recalculation flag - self.volume_update_required = True + See Also + -------- + dielectric_constant - def set_amount(self, solute, amount): """ - Set the amount of 'solute' in the parent solution. + bjerrum_length = unit.e**2 / ( + 4 + * math.pi + * self.dielectric_constant + * unit.epsilon_0 + * unit.k + * self.temperature + ) + return bjerrum_length.to("nm") + + def get_osmotic_pressure(self): + """ + Return the osmotic pressure of the solution relative to pure water. Parameters ---------- - 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. + None Returns ------- - Nothing. The concentration of solute is modified. - + Quantity + The osmotic pressure of the solution relative to pure water in Pa See Also -------- - Solute.set_moles() - """ - # raise an error if a negative amount is specified - if unit(amount).magnitude < 0: - logger.error( - "Negative amount specified for solute %s. Concentration not changed." - % solute - ) + get_water_activity + get_osmotic_coefficient + get_salt - # 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 - elif unit(amount).dimensionality in ( - "[substance]/[length]**3", - "[mass]/[length]**3", - ): + Notes + ----- + Osmotic pressure is calculated based on the water activity [#]_ [#]_ : - # store the original volume for later - orig_volume = self.get_volume() + .. math:: \\Pi = {RT \\over V_w} \\ln a_w - # change the amount of the solute present to match the desired amount - self.get_solute(solute).set_moles( - amount, self.get_volume(), self.get_solvent_mass() - ) + 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. - # 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 + References + ---------- + .. [#] Sata, Toshikatsu. Ion Exchange Membranes: Preparation, Characterization, and Modification. Royal Society of Chemistry, 2004, p. 10. - # adjust the amount of solvent - target_mass = ( - target_vol.magnitude - / 1000 - * self.water_substance.rho - * unit.Quantity("1 kg") - ) - mw = self.get_solvent().get_molecular_weight() - target_mol = target_mass / mw - self.get_solvent().moles = target_mol + .. [#] http://en.wikipedia.org/wiki/Osmotic_pressure#Derivation_of_osmotic_pressure - else: + Examples + -------- + >>> s1=pyEQL.Solution() + >>> s1.get_osmotic_pressure() + 0.0 - # change the amount of the solute present - self.get_solute(solute).set_moles( - amount, self.get_volume(), self.get_solvent_mass() - ) + >>> s1 = pyEQL.Solution([['Na+','0.2 mol/kg'],['Cl-','0.2 mol/kg']]) + >>> soln.get_osmotic_pressure() + + """ + # TODO - tie this into parameter() and solvent() objects + partial_molar_volume_water = 1.82e-5 * unit("m ** 3/mol") - # 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.get_solvent_mass() <= unit("0 kg"): - logger.error("All solvent has been depleted from the solution") - return None - else: - self._update_volume() - - def get_osmolarity(self, activity_correction=False): - """Return the osmolarity of the solution in Osm/L - - Parameters - ---------- - 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. - """ - if activity_correction is True: - factor = self.get_osmotic_coefficient() - else: - factor = 1 - return factor * self.get_total_moles_solute() / self.get_volume().to("L") + osmotic_pressure = ( + -1 + * unit.R + * self.temperature + / partial_molar_volume_water + * math.log(self.get_water_activity()) + ) + logger.info( + "Computed osmotic pressure of solution as %s Pa at T= %s degrees C" + % (osmotic_pressure, self.temperature) + ) + return osmotic_pressure.to("Pa") - def get_osmolality(self, activity_correction=False): - """Return the osmolality of the solution in Osm/kg + # Concentration Methods - Parameters - ---------- - 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. + def p(self, solute, activity=True): """ - if activity_correction is True: - factor = self.get_osmotic_coefficient() - else: - factor = 1 - return factor * self.get_total_moles_solute() / self.get_solvent_mass().to("kg") - - def get_total_moles_solute(self): - """Return the total moles of all solute in the solution""" - tot_mol = 0 - for item in self.components: - if item != self.solvent_name: - tot_mol += self.components[item].get_moles() - return tot_mol + Return the negative log of the activity of solute. - def get_moles_solvent(self): - """ - Return the moles of solvent present in the solution + Generally used for expressing concentration of hydrogen ions (pH) Parameters ---------- - None + 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 moles of solvent in the solution. - - """ + The negative log10 of the activity (or molar concentration if + activity = False) of the solute. - return self.get_amount(self.solvent_name, "mol") + Examples + -------- + TODO - def get_salt(self): """ - Determine the predominant salt in a solution of ions. + try: + if activity is True: + return -1 * math.log10(self.get_activity(solute)) + elif activity is False: + return -1 * math.log10(self.get_amount(solute, "mol/L").magnitude) + # if the solute has zero concentration, the log will generate a ValueError + except ValueError: + return 0 - 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 direclty model - these quantities. pyEQL works around this problem by treating such solutions - as single salt solutions. + def get_amount(self, solute, units): + """ + Return the amount of 'solute' in the parent solution. - 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., 1 M MgCl2 yields 1 M Mg+2 and 2 M Cl-). + The amount of a solute can be given in a variety of unit types. + 1. substance per volume (e.g., 'mol/L') + 1. substance per mass of solvent (e.g., 'mol/kg') + 1. mass of substance (e.g., 'kg') + 1. moles of substance ('mol') + 1. mole fraction ('fraction') + 1. percent by weight (%) Parameters ---------- - None + 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 ------- - Salt - Salt object containing information about the parent salt. + The amount of the solute in question, in the specified units - See Also - -------- - get_activity - get_activity_coefficient - get_water_activity - get_osmotic_coefficient - get_osmotic_pressure - get_viscosity_kinematic - Examples + See Also -------- - >>> s1 = Solution([['Na+','0.5 mol/kg'],['Cl-','0.5 mol/kg']]) - >>> s1.get_salt() - - >>> 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 - """ - # identify the predominant salt in the solution - return identify_salt(self) - - def get_salt_list(self): + add_amount + set_amount + get_total_amount + get_osmolarity + get_osmolality + get_solvent_mass + get_mass + get_total_moles_solute """ - Determine the predominant salt in a solution of ions. + # retrieve the number of moles of solute and its molecular weight + try: + moles = self.get_solute(solute).get_moles() + mw = self.get_solute(solute).get_molecular_weight() + # 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 0 * unit(units) + except DimensionalityError: + logger.error("Unsupported unit specified for get_amount") + return 0 - 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 direclty model - these quantities. + # 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 unecessary + # function calls. + if units == "fraction": + return moles / (self.get_moles_solvent() + self.get_total_moles_solute()) + elif units == "%": + return moles.to("kg", "chem", mw=mw) / self.mass.to("kg") * 100 + elif unit(units).dimensionality in ( + "[substance]/[length]**3", + "[mass]/[length]**3", + ): + return moles.to(units, "chem", mw=mw, volume=self.get_volume()) + elif unit(units).dimensionality in ("[substance]/[mass]", "[mass]/[mass]"): + return moles.to(units, "chem", mw=mw, solvent_mass=self.get_solvent_mass()) + elif unit(units).dimensionality == "[mass]": + return moles.to(units, "chem", mw=mw) + elif unit(units).dimensionality == "[substance]": + return moles.to(units) + else: + logger.error("Unsupported unit specified for get_amount") + return None - The get_salt_list() method examines the ionic composition of a solution and - simplifies it into a list of salts. The method retuns 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-). + def get_total_amount(self, element, units): + """ + Return the total amount of 'element' (across all solutes) in the solution. Parameters ---------- - None + element : str + String representing the name of the element of interest + units : str + Units desired for the output. Examples of valid units are + 'mol/L','mol/kg','mol', 'kg', and 'g/L' Returns ------- - dict - A dictionary of Salt objects, keyed to the salt formula + The total amount of the element in the solution, in the specified units + + Notes + ----- + There is currently no way to distinguish between different oxidation + states of the same element (e.g. TOTFe(II) vs. TOTFe(III)). This + is planned for a future release. (TODO) See Also -------- - get_activity - get_activity_coefficient - get_water_activity - get_osmotic_coefficient - get_osmotic_pressure - get_viscosity_kinematic - + get_amount """ - # identify the predominant salt in the solution - return generate_salt_list(self, unit="mol/kg") + import pyEQL.chemical_formula as ch - # Activity-related methods - def get_activity_coefficient( - self, - solute: str, - scale: Literal["molal", "molar", "fugacity", "rational"] = "molal", - verbose: bool = False, - ): - """ - Return the activity coefficient of a solute in solution. + TOT = 0 * unit(units) - The model used to calculte the activity coefficient is determined by the Solution's equation of state - engine. + # loop through all the solutes, process each one containing element + for item in self.components: + # check whether the solute contains the element + if ch.contains(item, element): + # start with the amount of the solute in the desired units + amt = self.get_amount(item, units) - Args: - solute: The solute for which to retrieve the activity coefficient - scale: The activity coefficient concentration scale - verbose: If True, pyEQL will print a message indicating the parent salt - that is being used for activity calculations. This option is - useful when modeling multicomponent solutions. False by default. + # convert the solute amount into the amount of element by + # either the mole / mole or weight ratio + if unit(units).dimensionality in ( + "[substance]", + "[substance]/[length]**3", + "[substance]/[mass]", + ): + TOT += amt * ch.get_element_mole_ratio(item, element) - 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 unit("1 dimensionless") - else: - try: - # get the molal-scale activity coefficient from the EOS engine - molal = self.engine.get_activity_coefficient( - solution=self, solute=solute - ) - except ValueError: - logger.warning( - "Calculation unsuccessful. Returning unit activity coefficient." - ) - return unit("1 dimensionless") + elif unit(units).dimensionality in ( + "[mass]", + "[mass]/[length]**3", + "[mass]/[mass]", + ): + TOT += amt * ch.get_element_weight_fraction(item, element) - # if necessary, convert the activity coefficient to another scale, and return the result - if scale == "molal": - return molal - elif scale == "molar": - total_molality = self.get_total_moles_solute() / self.get_solvent_mass() - total_molarity = self.get_total_moles_solute() / self.get_volume() - return ( - molal - * self.water_substance.rho - * unit.Quantity("1 g/L") - * total_molality - / total_molarity - ).to("dimensionless") - elif scale == "rational": - return molal * ( - 1 - + unit("0.018 kg/mol") - * self.get_total_moles_solute() - / self.get_solvent_mass() - ) - else: - logger.warning( - "Invalid scale argument. Returning molal-scale activity coefficient" - ) - return molal + return TOT - def get_activity( - self, - solute: str, - scale: Literal["molal", "molar", "rational"] = "molal", - verbose: bool = False, - ): + def add_amount(self, solute, amount): """ - Return the thermodynamic activity of the solute in solution on the chosen concentration scale. + Add the amount of 'solute' to the parent solution. - Args: - solute : str + Parameters + ---------- + solute : str String representing the name of the solute of interest - scale : str, optional - 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. - verbose : bool, optional - If True, pyEQL will print a message indicating the parent salt - that is being used for activity calculations. This option is - useful when modeling multicomponent solutions. False by default. + 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 ------- - The thermodynamic activity of the solute in question (dimensionless) + Nothing. The concentration of solute is modified. + See Also -------- - get_activity_coefficient - get_ionic_strength - get_salt - - Notes - ----- - The thermodynamic activity depends on the concentration scale used [#]. - By default, the ionic strength, activity coefficients, and activities are all - calculated based on the molal (mol/kg) concentration scale. + Solute.add_moles() + """ - References - ---------- - .. [#] Robinson, R. A.; Stokes, R. H. Electrolyte Solutions: Second Revised - Edition; Butterworths: London, 1968, p.32. + # 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 unit(amount).dimensionality in ( + "[substance]/[length]**3", + "[mass]/[length]**3", + ): - """ - # switch to the water activity function if the species is H2O - if solute == "H2O" or solute == "water": - activity = self.get_water_activity() - else: - # determine the concentration units to use based on the desired scale - if scale == "molal": - unit = "mol/kg" - elif scale == "molar": - unit = "mol/L" - elif scale == "rational": - unit = "fraction" - else: - logger.error("Invalid scale argument. Returning molal-scale activity.") - unit = "mol/kg" - scale = "molal" + # store the original volume for later + orig_volume = self.get_volume() - activity = ( - self.get_activity_coefficient(solute, scale=scale, verbose=verbose) - * self.get_amount(solute, unit).magnitude - ) - logger.info( - "Calculated %s scale activity of solute %s as %s" - % (scale, solute, activity) + # change the amount of the solute present to match the desired amount + self.get_solute(solute).add_moles( + amount, self.get_volume(), self.get_solvent_mass() ) - return activity - - # TODO - engine method - def get_osmotic_coefficient( - self, scale: Literal["molal", "molar", "rational"] = "molal" - ): - """ - Return the osmotic coefficient of an aqueous solution. + # set the amount to zero and log a warning if the desired amount + # change would result in a negative concentration + if self.get_amount(solute, "mol").magnitude < 0: + logger.warning( + "Attempted to set a negative concentration for solute %s. Concentration set to 0" + % solute + ) + self.set_amount(solute, "0 mol") - The method used depends on the Solution object's equation of state engine. + # calculate the volume occupied by all the solutes + solute_vol = self._get_solute_volume() - """ - molal_phi = self.engine.get_osmotic_coefficient(self) + # determine the volume of solvent that will preserve the original volume + target_vol = orig_volume - solute_vol - if scale == "molal": - return molal_phi - elif scale == "rational": - solvent = self.get_solvent().formula - return ( - -molal_phi - * unit("0.018 kg/mol") - * self.get_total_moles_solute() - / self.get_solvent_mass() - / math.log(self.get_amount(solvent, "fraction")) - ) - elif scale == "fugacity": - solvent = self.get_solvent().formula - return math.exp( - -molal_phi - * unit("0.018 kg/mol") - * self.get_total_moles_solute() - / self.get_solvent_mass() - - math.log(self.get_amount(solvent, "fraction")) - ) - else: - logger.warning( - "Invalid scale argument. Returning molal-scale osmotic coefficient" + # adjust the amount of solvent + # volume in L, density in kg/m3 = g/L + target_mass = ( + target_vol.magnitude * self.water_substance.rho * unit.Quantity("1 g") ) - return molal_phi - - def get_water_activity(self): - """ - Return the water activity. - Returns - ------- - Quantity : - The thermodynamic activity of water in the solution. + mw = self.get_solvent().get_molecular_weight() + target_mol = target_mass / mw + self.get_solvent().moles = target_mol - See Also - -------- - get_osmotic_coefficient - get_ionic_strength - get_salt + else: - Notes - ----- - Water activity is related to the osmotic coefficient in a solution containing i solutes by: [#]_ + # change the amount of the solute present + self.get_solute(solute).add_moles( + amount, self.get_volume(), self.get_solvent_mass() + ) - .. math:: \\ln a_w = - \\Phi M_w \\sum_i m_i + # set the amount to zero and log a warning if the desired amount + # change would result in a negative concentration + if self.get_amount(solute, "mol").magnitude < 0: + logger.warning( + "Attempted to set a negative concentration for solute %s. Concentration set to 0" + % solute + ) + self.set_amount(solute, "0 mol") - 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. + # 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.get_solvent_mass() <= unit("0 kg"): + logger.error("All solvent has been depleted from the solution") + return None + else: + # set the volume recalculation flag + self.volume_update_required = True - If appropriate Pitzer model parameters are not available, the - water activity is assumed equal to the mole fraction of water. + def set_amount(self, solute, amount): + """ + Set the amount of 'solute' in the parent solution. - References + Parameters ---------- - .. [#] 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. + 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. - Examples - -------- - >>> s1 = pyEQL.Solution([['Na+','0.3 mol/kg'],['Cl-','0.3 mol/kg']]) - >>> s1.get_water_activity() - - """ - """ - pseudo code + 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. - identify predominant salt for coefficients - check if coefficients exist for that salt - if so => calc osmotic coefficient and log an info message + Returns + ------- + Nothing. The concentration of solute is modified. - if not = > return mole fraction and log a warning message + See Also + -------- + Solute.set_moles() """ - osmotic_coefficient = self.get_osmotic_coefficient() - - if osmotic_coefficient == 1: - logger.warning( - "Pitzer parameters not found. Water activity set equal to mole fraction" + # raise an error if a negative amount is specified + if unit(amount).magnitude < 0: + logger.error( + "Negative amount specified for solute %s. Concentration not changed." + % solute ) - return self.get_amount("H2O", "fraction") - else: - concentration_sum = unit("0 mol/kg") - for item in self.components: - if item == "H2O": - pass - else: - concentration_sum += self.get_amount(item, "mol/kg") - logger.info("Calculated water activity using osmotic coefficient") + # 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 + elif unit(amount).dimensionality in ( + "[substance]/[length]**3", + "[mass]/[length]**3", + ): - return math.exp( - -osmotic_coefficient * 0.018015 * unit("kg/mol") * concentration_sum - ) * unit("1 dimensionless") + # store the original volume for later + orig_volume = self.get_volume() - @property - def ionic_strength(self) -> Quantity: - """ - Return the ionic strength of the solution. + # change the amount of the solute present to match the desired amount + self.get_solute(solute).set_moles( + amount, self.get_volume(), self.get_solvent_mass() + ) - Return the ionic strength of the solution, calculated as 1/2 * sum ( molality * charge ^2) over all the ions. + # calculate the volume occupied by all the solutes + solute_vol = self._get_solute_volume() - Molal (mol/kg) scale concentrations are used for compatibility with the activity correction formulas. + # determine the volume of solvent that will preserve the original volume + target_vol = orig_volume - solute_vol - Returns - ------- - Quantity : - The ionic strength of the parent solution, mol/kg. + # adjust the amount of solvent + target_mass = ( + target_vol.magnitude + / 1000 + * self.water_substance.rho + * unit.Quantity("1 kg") + ) + mw = self.get_solvent().get_molecular_weight() + target_mol = target_mass / mw + self.get_solvent().moles = target_mol - See Also - -------- - get_activity - get_water_activity + else: - Notes - ----- - The ionic strength is calculated according to: + # change the amount of the solute present + self.get_solute(solute).set_moles( + amount, self.get_volume(), self.get_solvent_mass() + ) - .. math:: I = \\sum_i m_i z_i^2 + # 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.get_solvent_mass() <= unit("0 kg"): + logger.error("All solvent has been depleted from the solution") + return None + else: + self._update_volume() - Where :math:`m_i` is the molal concentration and :math:`z_i` is the charge on species i. + def get_osmolarity(self, activity_correction=False): + """Return the osmolarity of the solution in Osm/L - Examples - -------- - >>> s1 = pyEQL.Solution([['Na+','0.2 mol/kg'],['Cl-','0.2 mol/kg']]) - >>> s1.ionic_strength - + Parameters + ---------- + 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. + """ + if activity_correction is True: + factor = self.get_osmotic_coefficient() + else: + factor = 1 + return factor * self.get_total_moles_solute() / self.get_volume().to("L") - >>> 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 - + def get_osmolality(self, activity_correction=False): + """Return the osmolality of the solution in Osm/kg + + Parameters + ---------- + 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. """ - ionic_strength = 0 - for solute in self.components.keys(): - ionic_strength += ( - 0.5 - * self.get_amount(solute, "mol/kg") - * self.components[solute].get_formal_charge() ** 2 - ) + if activity_correction is True: + factor = self.get_osmotic_coefficient() + else: + factor = 1 + return factor * self.get_total_moles_solute() / self.get_solvent_mass().to("kg") - return ionic_strength + def get_total_moles_solute(self): + """Return the total moles of all solute in the solution""" + tot_mol = 0 + for item in self.components: + if item != self.solvent_name: + tot_mol += self.components[item].get_moles() + return tot_mol - @property - def charge_balance(self) -> float: + def get_moles_solvent(self): """ - Return the charge balance of the solution. + Return the moles of solvent present in the solution - Return the charge balance of the solution. The charge balance represents the net electric charge - on the solution and SHOULD equal zero at all times, but due to numerical errors will usually - have a small nonzero value. + Parameters + ---------- + None Returns ------- - float : - The charge balance of the solution, in equivalents. - - Notes - ----- - The charge balance is calculated according to: - - .. math:: CB = F \\sum_i n_i z_i - - Where :math:`n_i` is the number of moles, :math:`z_i` is the charge on species i, and :math:`F` is the Faraday constant. + Quantity + The moles of solvent in the solution. """ - charge_balance = 0 - F = (unit.e * unit.N_A).magnitude - for solute in self.components.keys(): - charge_balance += ( - self.get_amount(solute, "mol").magnitude - * self.components[solute].get_formal_charge() - * F - ) - return charge_balance + return self.get_amount(self.solvent_name, "mol") - def get_alkalinity(self): + def get_salt(self): """ - 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: [#]_ + Determine the predominant salt in a solution of ions. - .. math:: Alk = F \\sum_i z_i C_B - \\sum_i z_i C_A + 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 direclty model + these quantities. pyEQL works around this problem by treating such solutions + as single salt solutions. - 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 charge. - In this method, the set of conservative cations is all Group I and Group II cations, and the conservative anions - are all the anions of strong acids. + 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., 1 M MgCl2 yields 1 M Mg+2 and 2 M Cl-). - References + Parameters ---------- - .. [#] Stumm, Werner and Morgan, James J. Aquatic Chemistry, 3rd ed, - pp 165. Wiley Interscience, 1996. - """ - alkalinity = 0 * unit("mol/L") - equiv_wt_CaCO3 = 100.09 / 2 * unit("g/mol") + None - base_cations = [ - "Li+", - "Na+", - "K+", - "Rb+", - "Cs+", - "Fr+", - "Be+2", - "Mg+2", - "Ca+2", - "Sr+2", - "Ba+2", - "Ra+2", - ] - acid_anions = ["Cl-", "Br-", "I-", "SO4-2", "NO3-", "ClO4-", "ClO3-"] + Returns + ------- + Salt + Salt object containing information about the parent salt. - for item in self.components: - if item in base_cations: - z = self.get_solute(item).get_formal_charge() - alkalinity += self.get_amount(item, "mol/L") * z - if item in acid_anions: - z = self.get_solute(item).get_formal_charge() - alkalinity -= self.get_amount(item, "mol/L") * z + See Also + -------- + get_activity + get_activity_coefficient + get_water_activity + get_osmotic_coefficient + get_osmotic_pressure + get_viscosity_kinematic - # convert the alkalinity to mg/L as CaCO3 - return (alkalinity * equiv_wt_CaCO3).to("mg/L") + Examples + -------- + >>> s1 = Solution([['Na+','0.5 mol/kg'],['Cl-','0.5 mol/kg']]) + >>> s1.get_salt() + + >>> s1.get_salt().formula + 'NaCl' + >>> s1.get_salt().nu_cation + 1 + >>> s1.get_salt().z_anion + -1 - def get_hardness(self): + >>> 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 """ - Return the hardness of a solution. + # identify the predominant salt in the solution + return identify_salt(self) - Hardness is defined as the sum of the equivalent concentrations - of multivalent cations as calcium carbonate. + def get_salt_list(self): + """ + Determine the predominant salt in a solution of ions. - NOTE: at present pyEQL cannot distinguish between mg/L as CaCO3 - and mg/L units. Use with caution. + 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 direclty model + these quantities. + + The get_salt_list() method examines the ionic composition of a solution and + simplifies it into a list of salts. The method retuns 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-). Parameters ---------- @@ -1628,120 +1489,262 @@ def get_hardness(self): Returns ------- - Quantity - The hardness of the solution in mg/L as CaCO3 - - """ - hardness = 0 * unit("mol/L") - equiv_wt_CaCO3 = 100.09 / 2 * unit("g/mol") + dict + A dictionary of Salt objects, keyed to the salt formula - for item in self.components: - z = self.get_solute(item).get_formal_charge() - if z > 1: - hardness += z * self.get_amount(item, "mol/L") + See Also + -------- + get_activity + get_activity_coefficient + get_water_activity + get_osmotic_coefficient + get_osmotic_pressure + get_viscosity_kinematic - # convert the hardness to mg/L as CaCO3 - return (hardness * equiv_wt_CaCO3).to("mg/L") + """ + # identify the predominant salt in the solution + return generate_salt_list(self, unit="mol/kg") - def get_debye_length(self): + # Activity-related methods + def get_activity_coefficient( + self, + solute: str, + scale: Literal["molal", "molar", "fugacity", "rational"] = "molal", + verbose: bool = False, + ): """ - Return the Debye length of a solution + Return the activity coefficient of a solute in solution. - Debye length is calculated as [#]_ + The model used to calculte the activity coefficient is determined by the Solution's equation of state + engine. - .. math:: + Args: + solute: The solute for which to retrieve the activity coefficient + scale: The activity coefficient concentration scale + verbose: If True, pyEQL will print a message indicating the parent salt + that is being used for activity calculations. This option is + useful when modeling multicomponent solutions. False by default. - \\kappa^{-1} = \\sqrt({\\epsilon_r \\epsilon_o k_B T \\over (2 N_A e^2 I)}) + 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 unit("1 dimensionless") + else: + try: + # get the molal-scale activity coefficient from the EOS engine + molal = self.engine.get_activity_coefficient( + solution=self, solute=solute + ) + except ValueError: + logger.warning( + "Calculation unsuccessful. Returning unit activity coefficient." + ) + return unit("1 dimensionless") - 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. + # if necessary, convert the activity coefficient to another scale, and return the result + if scale == "molal": + return molal + elif scale == "molar": + total_molality = self.get_total_moles_solute() / self.get_solvent_mass() + total_molarity = self.get_total_moles_solute() / self.get_volume() + return ( + molal + * self.water_substance.rho + * unit.Quantity("1 g/L") + * total_molality + / total_molarity + ).to("dimensionless") + elif scale == "rational": + return molal * ( + 1 + + unit("0.018 kg/mol") + * self.get_total_moles_solute() + / self.get_solvent_mass() + ) + else: + logger.warning( + "Invalid scale argument. Returning molal-scale activity coefficient" + ) + return molal - Parameters - ---------- - None + def get_activity( + self, + solute: str, + scale: Literal["molal", "molar", "rational"] = "molal", + verbose: bool = False, + ): + """ + Return the thermodynamic activity of the solute in solution on the chosen concentration scale. + + Args: + solute : str + String representing the name of the solute of interest + scale : str, optional + 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. + verbose : bool, optional + If True, pyEQL will print a message indicating the parent salt + that is being used for activity calculations. This option is + useful when modeling multicomponent solutions. False by default. Returns ------- - Quantity - The Debye length, in nanometers. + The thermodynamic activity of the solute in question (dimensionless) + + See Also + -------- + get_activity_coefficient + get_ionic_strength + get_salt + + Notes + ----- + The thermodynamic activity depends on the concentration scale used [#]. + By default, the ionic strength, activity coefficients, and activities are all + calculated based on the molal (mol/kg) concentration scale. References ---------- - .. [#] https://en.wikipedia.org/wiki/Debye_length#Debye_length_in_an_electrolyte + .. [#] Robinson, R. A.; Stokes, R. H. Electrolyte Solutions: Second Revised + Edition; Butterworths: London, 1968, p.32. - See Also - -------- - ionic_strength - get_dielectric_constant() + """ + # switch to the water activity function if the species is H2O + if solute == "H2O" or solute == "water": + activity = self.get_water_activity() + else: + # determine the concentration units to use based on the desired scale + if scale == "molal": + unit = "mol/kg" + elif scale == "molar": + unit = "mol/L" + elif scale == "rational": + unit = "fraction" + else: + logger.error("Invalid scale argument. Returning molal-scale activity.") + unit = "mol/kg" + scale = "molal" + + activity = ( + self.get_activity_coefficient(solute, scale=scale, verbose=verbose) + * self.get_amount(solute, unit).magnitude + ) + logger.info( + "Calculated %s scale activity of solute %s as %s" + % (scale, solute, activity) + ) + + return activity + # TODO - engine method + def get_osmotic_coefficient( + self, scale: Literal["molal", "molar", "rational"] = "molal" + ): """ - # to preserve dimensionality, convert the ionic strength into mol/L units - ionic_strength = self.ionic_strength.magnitude * unit("mol/L") - dielectric_constant = self.get_dielectric_constant() + Return the osmotic coefficient of an aqueous solution. - debye_length = ( - dielectric_constant - * unit.epsilon_0 - * unit.k - * self.temperature - / (2 * unit.N_A * unit.e**2 * ionic_strength) - ) ** 0.5 + The method used depends on the Solution object's equation of state engine. - return debye_length.to("nm") + """ + molal_phi = self.engine.get_osmotic_coefficient(self) + + if scale == "molal": + return molal_phi + elif scale == "rational": + solvent = self.get_solvent().formula + return ( + -molal_phi + * unit("0.018 kg/mol") + * self.get_total_moles_solute() + / self.get_solvent_mass() + / math.log(self.get_amount(solvent, "fraction")) + ) + elif scale == "fugacity": + solvent = self.get_solvent().formula + return math.exp( + -molal_phi + * unit("0.018 kg/mol") + * self.get_total_moles_solute() + / self.get_solvent_mass() + - math.log(self.get_amount(solvent, "fraction")) + ) + else: + logger.warning( + "Invalid scale argument. Returning molal-scale osmotic coefficient" + ) + return molal_phi - def get_bjerrum_length(self): + def get_water_activity(self): """ - Return the Bjerrum length of a solution + Return the water activity. - Bjerrum length representes the distance at which electrostatic - interactions between particles become comparable in magnitude - to the thermal energy.:math:`\\lambda_B` is calculated as [#]_ + Returns + ------- + Quantity : + The thermodynamic activity of water in the solution. - .. math:: + See Also + -------- + get_osmotic_coefficient + get_ionic_strength + get_salt - \\lambda_B = {e^2 \\over (4 \\pi \\epsilon_r \\epsilon_o k_B T)} + Notes + ----- + Water activity is related to the osmotic coefficient in a solution containing i solutes by: [#]_ - 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. + .. math:: \\ln a_w = - \\Phi M_w \\sum_i m_i - Parameters - ---------- - None + 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. - Returns - ------- - Quantity - The Bjerrum length, in nanometers. + If appropriate Pitzer model parameters are not available, the + water activity is assumed equal to the mole fraction of water. References ---------- - .. [#] https://en.wikipedia.org/wiki/Bjerrum_length + .. [#] 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() - >>> s1.get_bjerrum_length() - + >>> s1 = pyEQL.Solution([['Na+','0.3 mol/kg'],['Cl-','0.3 mol/kg']]) + >>> s1.get_water_activity() + + """ + """ + pseudo code - See Also - -------- - get_dielectric_constant() + identify predominant salt for coefficients + check if coefficients exist for that salt + if so => calc osmotic coefficient and log an info message + + if not = > return mole fraction and log a warning message """ - dielectric_constant = self.get_dielectric_constant() + osmotic_coefficient = self.get_osmotic_coefficient() - bjerrum_length = unit.e**2 / ( - 4 - * math.pi - * dielectric_constant - * unit.epsilon_0 - * unit.k - * self.temperature - ) - return bjerrum_length.to("nm") + if osmotic_coefficient == 1: + logger.warning( + "Pitzer parameters not found. Water activity set equal to mole fraction" + ) + return self.get_amount("H2O", "fraction") + else: + concentration_sum = unit("0 mol/kg") + for item in self.components: + if item == "H2O": + pass + else: + concentration_sum += self.get_amount(item, "mol/kg") + + logger.info("Calculated water activity using osmotic coefficient") + + return math.exp( + -osmotic_coefficient * 0.018015 * unit("kg/mol") * concentration_sum + ) * unit("1 dimensionless") def get_transport_number(self, solute, activity_correction=False): """Calculate the transport number of the solute in the solution @@ -2595,3 +2598,175 @@ def get_charge_balance(self): """ return self.charge_balance + + @deprecated( + message="get_alkalinity() will be removed in the next release. Access directly via the property Solution.alkalinity" + ) + def get_alkalinity(self): + """ + 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: [#]_ + + .. math:: Alk = F \\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 charge. + In this method, the set of conservative cations is all Group I and Group II cations, and the conservative anions are all the anions of strong acids. + + References + ---------- + .. [#] Stumm, Werner and Morgan, James J. Aquatic Chemistry, 3rd ed, + pp 165. Wiley Interscience, 1996. + """ + return self.alkalinity + + @deprecated( + message="get_hardness() will be removed in the next release. Access directly via the property Solution.hardness" + ) + def get_hardness(self): + """ + 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. + + Parameters + ---------- + None + + Returns + ------- + Quantity + The hardness of the solution in mg/L as CaCO3 + + """ + return self.hardness + + @deprecated( + message="get_debye_length() will be removed in the next release. Access directly via the property Solution.debye_length" + ) + def get_debye_length(self): + """ + Return the Debye length of a solution + + Debye length is calculated as [#]_ + + .. 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. + + Parameters + ---------- + None + + Returns + ------- + Quantity + The Debye length, in nanometers. + + References + ---------- + .. [#] https://en.wikipedia.org/wiki/Debye_length#Debye_length_in_an_electrolyte + + See Also + -------- + ionic_strength + get_dielectric_constant() + + """ + return self.debye_length + + @deprecated( + message="get_bjerrum_length() will be removed in the next release. Access directly via the property Solution.bjerrum_length" + ) + def get_bjerrum_length(self): + """ + Return the Bjerrum length of a solution + + Bjerrum length representes 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. + + Parameters + ---------- + None + + Returns + ------- + Quantity + The Bjerrum length, in nanometers. + + References + ---------- + .. [#] https://en.wikipedia.org/wiki/Bjerrum_length + + Examples + -------- + >>> s1 = pyEQL.Solution() + >>> s1.get_bjerrum_length() + + + See Also + -------- + get_dielectric_constant() + + """ + return self.bjerrum_length + + @deprecated( + message="get_dielectric_constant() will be removed in the next release. Access directly via the property Solution.dielectric_constant" + ) + def get_dielectric_constant(self): + """ + Returns the dielectric constant of the solution. + + Parameters + ---------- + None + + Returns + ------- + Quantity: the dielectric constant of the solution, dimensionless. + + Notes + ----- + Implements the following equation as given by [#]_ + + .. 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 + ---------- + .. [#] [1] 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. + """ + return self.dielectric_constant diff --git a/tests/test_bulk_properties.py b/tests/test_bulk_properties.py index 122565d4..4e97baa6 100644 --- a/tests/test_bulk_properties.py +++ b/tests/test_bulk_properties.py @@ -5,7 +5,7 @@ This file contains tests for the bulk properties calculated by Solution class methods. Currently included methods are: -- get_hardness() +- hardness """ @@ -17,14 +17,14 @@ class test_hardness: """ - test the get_hardness() method + test the hardness method ------------------------------ """ # an empty solution should have zero hardness def test_empty_solution(self): s1 = pyEQL.Solution() - result = s1.get_hardness().magnitude + result = s1.hardness.magnitude expected = 0 assert result == expected @@ -32,7 +32,7 @@ def test_empty_solution(self): # a solution with only monovalent ions should have zero hardness def test_hardness_1(self): s1 = pyEQL.Solution([["Na+", "0.2 mol/L"], ["Cl-", "0.2 mol/L"]]) - result = s1.get_hardness().magnitude + result = s1.hardness.magnitude expected = 0 assert result == expected @@ -40,7 +40,7 @@ def test_hardness_1(self): # a solution with only multivalent anions should have zero hardness def test_hardness_2(self): s1 = pyEQL.Solution([["Na+", "0.4 mol/L"], ["SO4-2", "0.2 mol/L"]]) - result = s1.get_hardness().magnitude + result = s1.hardness.magnitude expected = 0 assert result == expected @@ -49,7 +49,7 @@ def test_hardness_2(self): # molar concentration (e.g. multiply by the charge) def test_hardness_3(self): s1 = pyEQL.Solution([["Fe+3", "0.1 mol/L"], ["Cl-", "0.3 mol/L"]]) - result = s1.get_hardness().magnitude + result = s1.hardness.magnitude expected = 15013.5 assert round(abs(result - expected), 7) == 0 @@ -70,7 +70,7 @@ def test_hardness_4(self): ["PO4-3", "0.1 mol/L"], ] ) - result = s1.get_hardness().magnitude + result = s1.hardness.magnitude expected = 35031.5 assert round(abs(result - expected), 7) == 0 @@ -78,7 +78,7 @@ def test_hardness_4(self): # the hardness should return g/L units def test_hardness_5(self): s1 = pyEQL.Solution([["Fe+3", "0.1 mol/L"], ["Cl-", "0.3 mol/L"]]) - result = str(s1.get_hardness().dimensionality) + result = str(s1.hardness.dimensionality) expected = "[mass] / [length] ** 3" assert result == expected diff --git a/tests/test_debye_length.py b/tests/test_debye_length.py index fed45192..1ffca1d5 100644 --- a/tests/test_debye_length.py +++ b/tests/test_debye_length.py @@ -35,7 +35,7 @@ def test_debye_length_1(self): """ """ s1 = pyEQL.Solution([["Na+", "0.1 mmol/L"], ["Cl-", "0.1 mmol/L"]]) - result = s1.get_debye_length().magnitude + result = s1.debye_length.magnitude expected = 31 assert np.isclose(result, expected, rtol=RTOL) @@ -44,7 +44,7 @@ def test_debye_length_2(self): """ """ s1 = pyEQL.Solution([["Na+", "10 mmol/L"], ["Cl-", "10 mmol/L"]]) - result = s1.get_debye_length().magnitude + result = s1.debye_length.magnitude expected = 3.1 assert np.isclose(result, expected, rtol=RTOL) @@ -53,7 +53,7 @@ def test_debye_length_3(self): """ """ s1 = pyEQL.Solution([["Na+", "0.2 mmol/L"], ["SO4-2", "0.1 mmol/L"]]) - result = s1.get_debye_length().magnitude + result = s1.debye_length.magnitude expected = 18 assert np.isclose(result, expected, rtol=RTOL) @@ -62,7 +62,7 @@ def test_debye_length_4(self): """ """ s1 = pyEQL.Solution([["Na+", "20 mmol/L"], ["SO4-2", "10 mmol/L"]]) - result = s1.get_debye_length().magnitude + result = s1.debye_length.magnitude expected = 1.8 assert np.isclose(result, expected, rtol=RTOL) diff --git a/tests/test_dielectric.py b/tests/test_dielectric.py index d4232ef1..ad749eb5 100644 --- a/tests/test_dielectric.py +++ b/tests/test_dielectric.py @@ -32,7 +32,7 @@ def test_dielectric_constant(self): """ s1 = pyEQL.Solution([["Na+", "4.4 mol/kg"], ["Cl-", "4.4 mol/kg"]]) - result = s1.get_dielectric_constant().magnitude + result = s1.dielectric_constant.magnitude expected = 46 assert np.isclose(result, expected, rtol=RTOL) @@ -43,7 +43,7 @@ def test_dielectric_constant2(self): """ s1 = pyEQL.Solution([["Na+", "2 mol/kg"], ["Cl-", "2 mol/kg"]]) - result = s1.get_dielectric_constant().magnitude + result = s1.dielectric_constant.magnitude expected = 58 assert np.isclose(result, expected, rtol=RTOL) @@ -54,7 +54,7 @@ def test_dielectric_constant3(self): """ s1 = pyEQL.Solution([["Na+", "1 mol/kg"], ["Cl-", "1 mol/kg"]]) - result = s1.get_dielectric_constant().magnitude + result = s1.dielectric_constant.magnitude expected = 66 assert np.isclose(result, expected, rtol=RTOL) @@ -65,7 +65,7 @@ def test_dielectric_constant4(self): """ s1 = pyEQL.Solution([["K+", "1 mol/kg"], ["Br-", "1 mol/kg"]]) - result = s1.get_dielectric_constant().magnitude + result = s1.dielectric_constant.magnitude expected = 67 assert np.isclose(result, expected, rtol=RTOL) @@ -76,7 +76,7 @@ def test_dielectric_constant5(self): """ s1 = pyEQL.Solution([["K+", "3.4 mol/kg"], ["Br-", "3.4 mol/kg"]]) - result = s1.get_dielectric_constant().magnitude + result = s1.dielectric_constant.magnitude expected = 51 assert np.isclose(result, expected, rtol=RTOL) @@ -87,7 +87,7 @@ def test_dielectric_constant6(self): """ s1 = pyEQL.Solution([["Li+", "5 mol/kg"], ["Cl-", "5 mol/kg"]]) - result = s1.get_dielectric_constant().magnitude + result = s1.dielectric_constant.magnitude expected = 39 assert np.isclose(result, expected, rtol=RTOL) @@ -98,7 +98,7 @@ def test_dielectric_constant7(self): """ s1 = pyEQL.Solution([["Li+", "1 mol/kg"], ["Cl-", "1 mol/kg"]]) - result = s1.get_dielectric_constant().magnitude + result = s1.dielectric_constant.magnitude expected = 64 assert np.isclose(result, expected, rtol=RTOL) @@ -110,7 +110,7 @@ def test_dielectric_constant8(self): """ s1 = pyEQL.Solution([["Li+", "12 mol/kg"], ["Cl-", "12 mol/kg"]]) - result = s1.get_dielectric_constant().magnitude + result = s1.dielectric_constant.magnitude expected = 24 assert np.isclose(result, expected, rtol=RTOL) @@ -121,7 +121,7 @@ def test_dielectric_constant9(self): """ s1 = pyEQL.Solution([["Rb+", "6.5 mol/kg"], ["Cl-", "6.5 mol/kg"]]) - result = s1.get_dielectric_constant().magnitude + result = s1.dielectric_constant.magnitude expected = 43 assert np.isclose(result, expected, rtol=RTOL) @@ -132,7 +132,7 @@ def test_dielectric_constant9(self): """ s1 = pyEQL.Solution([["Rb+", "2.1 mol/kg"], ["Cl-", "2.1 mol/kg"]]) - result = s1.get_dielectric_constant().magnitude + result = s1.dielectric_constant.magnitude expected = 59 assert np.isclose(result, expected, rtol=RTOL) @@ -143,7 +143,7 @@ def test_dielectric_constant9(self): """ s1 = pyEQL.Solution([["Rb+", "0.5 mol/kg"], ["Cl-", "0.5 mol/kg"]]) - result = s1.get_dielectric_constant().magnitude + result = s1.dielectric_constant.magnitude expected = 73 assert np.isclose(result, expected, rtol=RTOL)