[CHNOSZ-commits] r651 - in pkg/CHNOSZ: . R demo inst man tests/testthat
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Mar 19 08:15:54 CET 2021
Author: jedick
Date: 2021-03-19 08:15:53 +0100 (Fri, 19 Mar 2021)
New Revision: 651
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/NAMESPACE
pkg/CHNOSZ/R/diagram.R
pkg/CHNOSZ/R/solubility.R
pkg/CHNOSZ/demo/E_coli.R
pkg/CHNOSZ/demo/Pourbaix.R
pkg/CHNOSZ/demo/zinc.R
pkg/CHNOSZ/inst/NEWS.Rd
pkg/CHNOSZ/inst/TODO
pkg/CHNOSZ/man/diagram.Rd
pkg/CHNOSZ/man/solubility.Rd
pkg/CHNOSZ/tests/testthat/test-solubility.R
Log:
solubility(): revamp arguments and support multiple minerals
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2021-03-16 14:29:58 UTC (rev 650)
+++ pkg/CHNOSZ/DESCRIPTION 2021-03-19 07:15:53 UTC (rev 651)
@@ -1,6 +1,6 @@
-Date: 2021-03-16
+Date: 2021-03-19
Package: CHNOSZ
-Version: 1.4.0-20
+Version: 1.4.0-21
Title: Thermodynamic Calculations and Diagrams for Geochemistry
Authors at R: c(
person("Jeffrey", "Dick", , "j3ffdick at gmail.com", role = c("aut", "cre"),
Modified: pkg/CHNOSZ/NAMESPACE
===================================================================
--- pkg/CHNOSZ/NAMESPACE 2021-03-16 14:29:58 UTC (rev 650)
+++ pkg/CHNOSZ/NAMESPACE 2021-03-19 07:15:53 UTC (rev 651)
@@ -58,9 +58,7 @@
"CHNOSZ", "thermo", "reset", "OBIGT", "retrieve", "AkDi", "moles",
"lNaCl", "lS", "lT", "lP", "lTP", "lex",
# added 20200716 or later
- "mash", "mix", "rebalance",
-# added 20210303 or later
- "solubilities"
+ "mash", "mix", "rebalance"
)
# Load shared objects
Modified: pkg/CHNOSZ/R/diagram.R
===================================================================
--- pkg/CHNOSZ/R/diagram.R 2021-03-16 14:29:58 UTC (rev 650)
+++ pkg/CHNOSZ/R/diagram.R 2021-03-19 07:15:53 UTC (rev 651)
@@ -46,8 +46,8 @@
if(length(efun)==0) efun <- ""
if(!(efun %in% c("affinity", "equilibrate") | grepl("solubilit", efun)))
stop("'eout' is not the output from affinity(), equilibrate(), solubility(), or solubilities()")
- # For solubilities(), type is always loga.balance 20210303
- if(identical(efun, "solubilities")) type <- "loga.balance"
+ # For solubilities(), default type is loga.balance 20210303
+ if(identical(efun, "solubilities") & missing(type)) type <- "loga.balance"
## 'type' can be:
# 'auto' - property from affinity() (1D) or maximum affinity (affinity 2D) (aout) or loga.equil (eout)
Modified: pkg/CHNOSZ/R/solubility.R
===================================================================
--- pkg/CHNOSZ/R/solubility.R 2021-03-16 14:29:58 UTC (rev 650)
+++ pkg/CHNOSZ/R/solubility.R 2021-03-19 07:15:53 UTC (rev 651)
@@ -1,8 +1,9 @@
# solubility.R: vectorized solubility calculations without uniroot
-# 20181031 jmd
+# 20181031 first version jmd
# 20181106 work on the output from affinity(); no "equilibrate()" needed!
# 20190117 add find.IS and test for dissociation reaction
# 20190730 add codeanal; test all species for dissociation reaction
+# 20210319 use list of aqueous species as main argument (with back-compatibility for affinity output) and handle multiple minerals
## if this file is interactively sourced, the following are also needed to provide unexported functions:
#source("equilibrate.R")
@@ -11,7 +12,70 @@
#source("util.args.R")
#source("util.character.R")
-solubility <- function(aout, dissociation = NULL, find.IS = FALSE, in.terms.of = NULL, codeanal = FALSE) {
+# Function to calculate solubilities of multiple minerals 20210303
+# species() should be used first to load the minerals (all bearing the same metal)
+# 'iaq' lists aqueous species that can be produced by dissolution of the minerals
+# '...' contains arguments for affinity() or mosaic() (i.e. plotting variables)
+# FIXME: what to do about 'dissociation' argument?
+solubility <- function(iaq, ..., in.terms.of = NULL, dissociation = NULL, find.IS = FALSE) {
+
+ # If iaq is the output of affinity(), use old method 20210318
+ if(is.list(iaq)) return(solubility_calc(aout = iaq, in.terms.of = in.terms.of, dissociation = dissociation, find.IS = find.IS))
+ # Check whether to use affinity() or mosaic()
+ ddd <- list(...)
+ if(identical(names(ddd)[1], "bases")) is.mosaic <- TRUE else is.mosaic <- FALSE
+
+ # Save current thermodynamic system settings
+ thermo <- get("thermo", CHNOSZ)
+ # Use current basis species as a template for the solubility calculations
+ ispecies <- basis()$ispecies
+ logact <- basis()$logact
+ # The current formed species are the minerals to be dissolved
+ mineral <- species()
+
+ # Make a list to store the calculated solubilities for each mineral
+ slist <- list()
+ # Loop over minerals
+ for(i in seq_along(mineral$ispecies)) {
+ # Print message
+ message(paste("solubility: calculating for", mineral$name[i]))
+ # Define basis species with the mineral first (so it will be dissolved)
+ ispecies[1] <- mineral$ispecies[i]
+ logact[1] <- mineral$logact[i]
+ basis(ispecies, logact)
+ # Add aqueous species (no need to define activities here - they will be calculated by solubility_calc)
+ species(iaq)
+ if(is.mosaic) a <- suppressMessages(mosaic(...)$A.species) else a <- suppressMessages(affinity(...))
+ # Calculate solubility of this mineral
+ scalc <- solubility_calc(a, in.terms.of = in.terms.of, dissociation = dissociation, find.IS = find.IS)
+ # Store the solubilities in the list
+ slist[[i]] <- scalc$loga.balance
+ }
+
+ # Restore the original thermodynamic system settings
+ assign("thermo", thermo, CHNOSZ)
+
+ if(length(mineral$ispecies) == 1) {
+ # For one mineral, return the results of the solubility calculation
+ scalc
+ } else {
+ # For multiple minerals, the overall solubility is the *minimum* among all the minerals
+ smin <- do.call(pmin, slist)
+ # Put this into the last-computed 'solubility' object
+ scalc$loga.balance <- smin
+ scalc$loga.equil <- slist
+ scalc$species <- mineral
+ # Change the function name stored in the object so diagram() plots loga.balance automatically
+ scalc$fun <- "solubilities"
+ # Return the object
+ scalc
+ }
+}
+
+
+# The "nuts and bolts" of solubility calculations
+# Moved from solubility() to solubility_calc() 20210318
+solubility_calc <- function(aout, dissociation = NULL, find.IS = FALSE, in.terms.of = NULL, codeanal = FALSE) {
## concept: the logarithms of activities of species at equilibrium are equal to
## Astar, the affinities calculated for unit activities of species
@@ -47,7 +111,7 @@
} else if(codeanal) print(" at least one formation reaction doesn't involve the second basis species")
} else if(codeanal) print("-- no dissociation reactions are possible")
message("solubility: test for dissociation reaction returns ", dissociation)
- } else message("solubility: from argument, dissociation ratio is ", dissociation)
+ }
# get starting ionic strength (probably zero, but could be anything set by user)
IS <- aout$IS
@@ -182,49 +246,3 @@
loga.balance=loga.balance, Astar=loga.equil, loga.equil=loga.equil))
}
-# Calculate solubilities of multiple minerals 20210303
-# a_cr: affinities for minerals (all bearing the same metal)
-# i_aq: aqueous species that can be produced by dissolution of the minerals
-# FIXME: what to do about 'dissociation' argument?
-solubilities <- function(a_cr, i_aq, in.terms.of = NULL, dissociation = NULL) {
- # If a_cr is the output from mosaic(), just get the species' affinities
- is.mosaic <- FALSE
- m_cr <- a_cr
- if(identical(a_cr$fun, "mosaic")) {
- a_cr <- a_cr$A.species
- is.mosaic <- TRUE
- }
- # Find all stable minerals across diagram
- d_cr <- diagram(a_cr, plot.it = FALSE)
- d_cr.stable <- d_cr$species$ispecies[unique(as.vector(d_cr$predominant))]
- # Use basis species in a_cr as a template for the solubility calculations
- ispecies <- a_cr$basis$ispecies
- logact <- a_cr$basis$logact
-
- # Make a list to store the calculated solubilities for each mineral
- slist <- list()
- # Loop over stable minerals
- for(i in seq_along(d_cr.stable)) {
- # Define basis species with the mineral first (so it will be dissolved)
- ispecies[1] <- d_cr.stable[i]
- basis(ispecies, logact)
- # Add aqueous species (no need to define activities here - they will be calculated)
- species(i_aq)
- # Calculate affinities of formation reactions of this mineral at same conditions as a_cr (argument recall)
- if(is.mosaic) a <- mosaic(m_cr)$A.species else a <- affinity(a_cr)
- # Calculate solubility of this mineral
- s <- solubility(a, in.terms.of = in.terms.of, dissociation = dissociation)
- # Store the solubilities in the list
- slist[[i]] <- s$loga.balance
- }
-
- # The overall solubility is the *minimum* among all the minerals
- smin <- do.call(pmin, slist)
- # Put this into the last-computed 'solubility' object
- s$loga.balance <- smin
- # Change the function name stored in the object so diagram() plots loga.balance automatically
- s$fun <- "solubilities"
- # Return the object
- s
-}
-
Modified: pkg/CHNOSZ/demo/E_coli.R
===================================================================
--- pkg/CHNOSZ/demo/E_coli.R 2021-03-16 14:29:58 UTC (rev 650)
+++ pkg/CHNOSZ/demo/E_coli.R 2021-03-19 07:15:53 UTC (rev 651)
@@ -6,6 +6,7 @@
# Polymerization scheme from Amend et al. (2013): https://doi.org/10.1098/rstb.2012.0255
library(CHNOSZ)
+opar <- par(no.readonly = TRUE)
# Concentrations of biomolecules (mol (g cell)-1)
# (from Table 1 of LaRowe and Amend, 2016)
@@ -79,7 +80,7 @@
loga_N <- switch(N, "NO3-" = -5, "NH4+" = -6)
loga_S <- switch(S, "SO4-2" = -3, "HS-" = -6)
# Set basis species
- # (Note: we set activity of e- in affinity())
+ # (Note: we set activity of e- in affinity() below)
basis(c(C, N, S, "HPO4-2", "H2O", "H+", "e-"), c(loga_C, loga_N, loga_S, -5, 0, -7, 0))
# Load formed species
species(biomolecules, -9)
@@ -122,11 +123,12 @@
}
-# Make plots with different combinations of basis species
+# Set graphical parameters
par(mfrow = c(3, 2))
par(mar = c(2.5, 3, 1.5, 1), mgp = c(1.5, 0.3, 0))
par(tcl = 0.25)
par(cex = 1)
+# Make plots with different combinations of basis species
plot_G("CO2", "NO3-", "SO4-2")
plot_G("CO2", "NH4+", "HS-")
plot_G("CH3COO-", "NO3-", "SO4-2")
@@ -136,3 +138,5 @@
# Reset CHNOSZ settings (units and OBIGT database)
reset()
+# Reset plot settings
+par(opar)
Modified: pkg/CHNOSZ/demo/Pourbaix.R
===================================================================
--- pkg/CHNOSZ/demo/Pourbaix.R 2021-03-16 14:29:58 UTC (rev 650)
+++ pkg/CHNOSZ/demo/Pourbaix.R 2021-03-19 07:15:53 UTC (rev 651)
@@ -82,9 +82,8 @@
# Calculate affinities for minerals
species(i_cr)
-a_cr <- affinity(pH = c(pH, res), Eh = c(Eh, res), T = T, P = P, IS = IS)
# Calculate overall solubility (i.e. minimum solubility given all candidate minerals)
-s <- solubilities(a_cr, i_aq, in.terms.of = element)
+s <- solubility(i_aq, pH = c(pH, res), Eh = c(Eh, res), T = T, P = P, IS = IS, in.terms.of = element)
# Plot diagram (LAYER 2: equisolubility lines)
diagram(s, levels = levels, contour.method = "flattest", add = TRUE, lwd = 1.5)
Modified: pkg/CHNOSZ/demo/zinc.R
===================================================================
--- pkg/CHNOSZ/demo/zinc.R 2021-03-16 14:29:58 UTC (rev 650)
+++ pkg/CHNOSZ/demo/zinc.R 2021-03-19 07:15:53 UTC (rev 651)
@@ -70,7 +70,7 @@
# (i.e. saturation condition for the solution)
# Use solubilities() 20210303
# FIXME: why do we need to set dissociation = FALSE here?
-s <- solubilities(mcr, iaq, in.terms.of = "Zn", dissociation = FALSE)
+s <- solubility(iaq, bases = bases, pH = pH, O2 = O2, T = T, P = P, IS = sat$IS, in.terms.of = "Zn", dissociation = FALSE)
# Specify contour levels
levels <- seq(-12, 9, 3)
diagram(s, type = "loga.balance", levels = levels, contour.method = "flattest")
Modified: pkg/CHNOSZ/inst/NEWS.Rd
===================================================================
--- pkg/CHNOSZ/inst/NEWS.Rd 2021-03-16 14:29:58 UTC (rev 650)
+++ pkg/CHNOSZ/inst/NEWS.Rd 2021-03-19 07:15:53 UTC (rev 651)
@@ -10,7 +10,7 @@
\newcommand{\s}{\ifelse{latex}{\eqn{_{#1}}}{\ifelse{html}{\out{<sub>#1</sub>}}{#1}}}
\newcommand{\S}{\ifelse{latex}{\eqn{^{#1}}}{\ifelse{html}{\out{<sup>#1</sup>}}{^#1}}}
-\section{Changes in CHNOSZ version 1.4.0-20 (2021-03-16)}{
+\section{Changes in CHNOSZ version 1.4.0-21 (2021-03-19)}{
\subsection{NEW FEATURES}{
\itemize{
@@ -19,9 +19,9 @@
lines, after \href{https://www.worldcat.org/oclc/563921897}{Pourbaix,
1974}). This demo depends on the following two changes.
- \item Add \code{solubilities()} for calculating overall (i.e. minimum)
- solubility for multiple minerals. This function is also now used in
- \samp{demo/zinc.R}.
+ \item Revamp \code{solubility()} to enable calculating overall (i.e.
+ minimum) solubility for multiple minerals. Calculations for multiple
+ minerals are also now used in \samp{demo/zinc.R}.
\item Restore \samp{lty.aq} and \samp{lty.cr} arguments to
\code{diagram()} to control plotting of aq-aq and cr-cr field boundaries.
Modified: pkg/CHNOSZ/inst/TODO
===================================================================
--- pkg/CHNOSZ/inst/TODO 2021-03-16 14:29:58 UTC (rev 650)
+++ pkg/CHNOSZ/inst/TODO 2021-03-19 07:15:53 UTC (rev 651)
@@ -45,8 +45,6 @@
- diagram(): scale predominant.values in output by balancing coefficients?
-- add HKF parameters from LaRowe and Amend, 2016 (https://doi.org/10.1038/ismej.2015.227)
-
[20210115]
BUG: Why doesn't this work?
Modified: pkg/CHNOSZ/man/diagram.Rd
===================================================================
--- pkg/CHNOSZ/man/diagram.Rd 2021-03-16 14:29:58 UTC (rev 650)
+++ pkg/CHNOSZ/man/diagram.Rd 2021-03-19 07:15:53 UTC (rev 651)
@@ -93,7 +93,7 @@
\details{
This function displays diagrams representing either chemical affinities, or equilibrium chemical activities of species.
-The first argument is the output from \code{\link{affinity}}, \code{\link{equilibrate}}, \code{\link{solubility}}, or \code{\link{solubilities}}.
+The first argument is the output from \code{\link{affinity}}, \code{\link{equilibrate}}, or \code{\link{solubility}}.
0-D diagrams, at a single point, are shown as \code{\link{barplot}}s.
1-D diagrams, for a single variable on the x-axis, are plotted as lines.
2-D diagrams, for two variables, are plotted as predominance fields.
@@ -126,12 +126,11 @@
In the case of 2-D diagrams, both of these options use \code{\link{contour}} to draw the lines, with the method specified in \code{contour.method}.
This paragraph describes the effect of the \code{type} argument when the output from \code{solubility} is being used.
-For \code{type} set to \samp{auto}, the equilibrium activities (i.e. computed solubilities) of each species are plotted.
+For one mineral or gas, if \code{type} set to \samp{auto}, the equilibrium activities of each aqueous species are plotted.
If \code{type} is \samp{loga.balance}, the activity of the balancing basis species (i.e. total solubility) is plotted; this is represented by contours on a 2-D diagram.
+For two or more minerals or gases, if \code{type} set to \samp{auto}, the values of \samp{loga.balance} (overall minimum solubility) are plotted.
+If \code{type} is \samp{loga.equil}, the solubilities of the individual minerals and gases are plotted.
For examples that use these features, see \code{\link{solubility}} and various \code{\link{demos}}: \samp{DEW}, \samp{contour}, \samp{gold}, \samp{solubility}, \samp{sphalerite}.
-
-When the output from \code{\link{solubilities}} is provided, the only supported option is to plot \samp{loga.balance}.
-See \code{demo("Pourbaix")} and \code{demo("zinc")}.
}
\section{1-D diagrams}{
Modified: pkg/CHNOSZ/man/solubility.Rd
===================================================================
--- pkg/CHNOSZ/man/solubility.Rd 2021-03-16 14:29:58 UTC (rev 650)
+++ pkg/CHNOSZ/man/solubility.Rd 2021-03-19 07:15:53 UTC (rev 651)
@@ -1,7 +1,6 @@
\encoding{UTF-8}
\name{solubility}
\alias{solubility}
-\alias{solubilities}
\title{Equilibrium Chemical Activities of Species}
\description{
Calculate chemical activities of aqueous species in equilibrium with a mineral or gas.
@@ -8,54 +7,48 @@
}
\usage{
- solubility(aout, dissociation = NULL, find.IS = FALSE, in.terms.of = NULL,
- codeanal = FALSE)
- solubilities(a_cr, i_aq, in.terms.of = NULL, dissociation = NULL)
+ solubility(iaq, ..., in.terms.of = NULL, dissociation = NULL, find.IS = FALSE)
}
\arguments{
- \item{aout}{list, output from \code{\link{affinity}}}
+ \item{iaq}{character (names) or numeric (species indices) of aqueous species}
+ \item{...}{arguments for \code{\link{affinity}} or \code{\link{mosaic}} (i.e. plotting variables)}
+ \item{in.terms.of}{character, express the total solubility in terms of moles of this species}
\item{dissociation}{logical, does the mineral undergo a dissociation reaction?}
\item{find.IS}{logical, find the equilibrium ionic strength by iteration?}
- \item{in.terms.of}{character, express the total solubility in terms of moles of this species}
- \item{codeanal}{logical, print diagnostic messages and values of internal variables for code analysis}
- \item{a_cr}{list, output from \code{\link{affinity}} or \code{\link{mosaic}} for one or more minerals}
- \item{i_aq}{character or numeric, aqueous species that can be formed by dissolving the minerals}
}
\details{
-\code{solubility} performs a simple task: given the values of \code{\link{affinity}} of formation reactions of species at specified activity, it works backward to find the activities of species that make the affinities zero.
-This corresponds to complete equilibrium with all of the basis species.
-Usually, the basis species should be set up so that the first basis species represents the substance being dissolved (a mineral such as CaCO\s3 or gas such as CO\s2).
-Internally, this is treated as the conserved basis species, so it must be present in all of the formation reactions of the species.
+\code{solubility} calculates the activities of aqueous species in equilibrium with one or more minerals or gases.
+The minerals or gases should be loaded as the formed \code{\link{species}}, and the aqueous species (including ions and/or neutral complexes) that can be produced by dissolution should be listed in the \code{iaq} argument.
+The definitions of plotting variables should be provided in \code{...}, which are passed as arguments to \code{\link{affinity}}, or to \code{\link{mosaic}} if the first one is named \code{bases}.
-To generate \code{aout}, the \code{\link{species}} should be defined to represent one set of aqueous species (including ions and/or neutral complexes) formed in solution, all involving the conserved basis species.
-For a dissociation reaction, the second basis species should be used to represent the counterion (cation or anion).
-Other variables (pH, ionic strength, activities of other basis species) should be defined in the call to \code{\link{affinity}} to make \code{aout}.
+It must be possible to obtain a valid set of basis species by substituting each of the minerals or gases in the first position of the current \code{\link{basis}} defintion, and all of the aqueous species must include that basis species in their formation reactions.
+(This essentially means that all minerals, gases and aqueous species must share a common element, which is what the reactions are balanced on.)
-A more complex calculation involves setting \code{find.IS} to TRUE in order to determine the final ionic strength of dissolving a substance in pure water.
-This works by calculating the ionic strength from the equilibrium solubility calculation, then re-running \code{affinity} with those values.
-Note that for dissociation reactions, the ionic strength is calculated from both the ions present in the species definition and the counter ion, which should be the second basis species.
-The calculation is iterated until the ionic strength deviation at every point is lower than a preset tolerance (1e-4).
-Alternatively, speciation of counterions (e.g. ionized forms of carbonate or sulfate) can also be accomplished by using the \code{\link{mosaic}} function instead of \code{affinity}; this is used in the second example below.
+The output of the function is structured to facilitate making useful diagrams.
+For a single mineral or gas, the output of \code{solubility} can be used by \code{\link{diagram}} with \code{type = NULL} (the default) to plot the activities of the aqueous species or with \code{type = "loga.balance"} to plot the sum of activities of aqueous species, which corresponds to the solubility of the mineral or gas.
+This value corresponds to the total extent of dissolution of the mineral or gas; \code{in.terms.of} can be used to express this value in terms of another species or element.
+For example, for dissolution of gaseous S\s{2}, \code{in.terms.of = "S"} gives the total amount of S in solution, which is twice the amount of S\s{2} dissolved.
+Likewise, the solubility of corundum (Al\s{2}O\s{3}) can be expressed in terms of the moles of Al\S{+3} in solution (see the vignette \viglink{anintro}).
-The output of \code{solubility} has the same format as that of \code{equilibrate}, and can be used by \code{\link{diagram}} with \code{type = "loga.balance"} to plot the solubilities, or with \code{type = NULL} to plot the activities of species.
-The value of \code{loga.balance} reflects the activity (or molality) of the conserved basis species, i.e. the thing being dissolved.
-Use \code{in.terms.of} to express this value in terms of another species.
-For example, the solubility of corundum (Al\s{2}O\s{3}) can be expressed in terms of the moles of Al\S{+3} in solution (see the vignette \viglink{anintro}).
+For multiple minerals, the function calculates the solubilities for each of the minerals separately; these are stored in the \code{loga.equil} element of the output.
+The overall \emph{minimum} solubility among all the minerals at each point is stored in \code{loga.balance}.
+This corresponds to the total activity of dissolved species in equilibrium with the most stable mineral.
+In contrast to the situation for a single mineral or gas, \code{\link{diagram}} be default plots \code{loga.balance}, and \code{type = "loga.balance"} should be used to plot the solubilities for the individual minerals or gases.
+}
-Whereas \code{solubility} is limited to calculations for a single mineral at a time, \code{solubilities} combines the results of solubility calculations for multiple minerals.
-\code{a_cr} should contain the affinities of only minerals with any set of basis species, as long as the first one contains the main element.
-\code{i_aq} gives the names or species indices of all the aqueous species that all have that element in combination with the other basis species (possibly generated using \code{\link{retrieve}}).
-The function calculates the solubilities for each of the minerals present in code{a_cr}.
-Then, it takes the \emph{minimum} solubility among all the minerals at each grid location.
-This corresponds to the maximum activity of dissolved species.
-At this point the solution is saturated (not oversaturated) with respect to the most stable mineral.
+\section{Backward Compatibility}{
+For compatibility with previous versions of the function, the \code{iaq} argument can be the output of \code{\link{affinity}} or \code{\link{mosaic}} for aqueous species.
+The examples for ionic strength and dissociation reactions were designed for this calling style.
-The result of \code{solubilities} is stored in the \code{loga.balance} component of the returned list object.
-Only this component is used by \code{\link{diagram}} for plotting; the other components, while present, refer to the solubility calculations for only the last mineral.
-This function is used for generating isosolubility (aka. equisolubility) lines in \code{demo("Pourbaix")} and \code{demo{"zinc"}}.
-The latter demo combines the calculation of solubilities with a \code{\link{mosaic}} calculation to account for the speciation of aqueous sulfur.
+In this case the (single) mineral or gas being dissolved is taken from the current \code{\link{basis}} species.
+Usually, the basis species should be set up so that the first basis species represents the substance being dissolved (a mineral such as CaCO\s3 or gas such as CO\s2).
+This is treated as the conserved basis species, so it must be present in all of the formation reactions of the aqueous species.
+
+The \code{\link{species}} should be defined to represent one set of aqueous species (including ions and/or neutral complexes) formed in solution, all involving the conserved basis species.
+For a dissociation reaction, the second basis species should be used to represent the counterion (cation or anion).
+Other variables (pH, ionic strength, activities of other basis species) should be defined in the call to \code{\link{affinity}} to make \code{iaq}.
}
\section{Dissociation Reactions}{
@@ -76,6 +69,14 @@
That does not give a complete equilibrium in the system, but may be required to reproduce some published diagrams.
}
+\section{Ionic Strength}{
+Set \code{find.IS} to TRUE to determine the final ionic strength due to dissolution of a substance in pure water.
+This works by calculating the ionic strength from the amounts of aqueous species formed, then re-running \code{affinity} with the calculated \code{IS} value.
+Note that for dissociation reactions, the ionic strength is calculated from both the ions present in the species definition and the counter ion, which should be the second basis species.
+The calculation is iterated until the ionic strength deviation at every point is lower than a preset tolerance (1e-4).
+Alternatively, speciation of counterions (e.g. ionized forms of carbonate or sulfate) can also be accomplished by using the \code{\link{mosaic}} function instead of \code{affinity}; this is used in the second example below.
+}
+
\section{Warning}{
This function has not been tested for systems that may form dimers or higher-order complexes (such as Au\s{2}S\s{2}\S{2-}).
Except for relatively simple systems, even after careful refinement, the results from CHNOSZ, which considers chemical activities as the independent variables, will not match the results from speciation-solubility (or Gibbs energy minimization) codes, where the system is defined by its bulk composition.
@@ -82,12 +83,16 @@
}
\seealso{
+\code{\link{retrieve}} provides a way to list all of the aqueous species in the database that have the specified elements.
+
\code{demo("solubility")} shows \T-pH diagrams related to the CO\s{2} and calcite example here.
\code{demo("gold")} shows solubility calculations for Au in aqueous solutions with hydroxide, chloride, and hydrosulfide complexes.
In the latter demo, the conserved basis species is not the first one (this might have been done to prevent the incorrect automatic detection of dissociation reactions).
-\code{\link{retrieve}} provides a way to list all of the aqueous species in the database that have the specified elements.
-In contrast to \code{solubility}, \code{\link{equilibrate}} calculates equilibrium chemical activities of species given a constant value of \code{loga.balance} (the logarithm of total activity of the conserved basis species).
+Solubility calculations for multiple minerals are used for generating isosolubility (aka. equisolubility) lines in \code{demo("Pourbaix")} and \code{demo{"zinc"}}.
+The latter demo combines the calculation of solubilities with a \code{\link{mosaic}} calculation to account for the speciation of aqueous sulfur.
+
+Whereas \code{solubility} yields a stable equilibrium condition (the affinities of formation reactions of aqueous species are zero), \code{\link{equilibrate}} generates metastable equilibrium (the affinities of formation reactions are equal to each other, but not necessarily zero).
}
\examples{\dontshow{reset()}
@@ -160,7 +165,7 @@
# The two methods give nearly equivalent results
stopifnot(all.equal(sa0$loga.balance, sm0$loga.balance))
stopifnot(all.equal(saI$loga.balance, smI$loga.balance, tolerance = 0.003))
-## NOTE: the second method (using mosaic) takes longer, but is
+## NOTE: the second method (using mosaic) is slower, but is
## more flexible; e.g. complexes with Ca+2 could be included
}
Modified: pkg/CHNOSZ/tests/testthat/test-solubility.R
===================================================================
--- pkg/CHNOSZ/tests/testthat/test-solubility.R 2021-03-16 14:29:58 UTC (rev 650)
+++ pkg/CHNOSZ/tests/testthat/test-solubility.R 2021-03-19 07:15:53 UTC (rev 651)
@@ -52,3 +52,31 @@
# TODO: write a test to get this error 20190730
#expect_error(solubility(a), "Unsure whether the first formation reaction is a dissociation reaction.")
})
+
+test_that("backward compatible and new calling styles produce identical results", {
+ # Test added 20210319
+ # Calculate solubility of a single substance:
+ # Gaseous S2 with a given fugacity
+
+ # Define basis species (any S-bearing basis species is allowed)
+ basis(c("HS-", "oxygen", "H2O", "H+"))
+ basis("pH", 6)
+ # Load the substances (minerals or gases) to be dissolved
+ species("S2", -20)
+ # List the formed aqueous species
+ i_aq <- info(c("SO4-2", "HS-"))
+ # Place arguments for affinity() or mosaic() after the first argument of solubility()
+ s1 <- solubility(i_aq, O2 = c(-55, -40), T = 125, in.terms.of = "SO4-2")
+
+ # Backward-compatible method limited to dissolving one species:
+ # Include S2(g) in the basis species
+ basis(c("S2", "oxygen", "H2O", "H+"))
+ basis("pH", 6)
+ basis("S2", -20)
+ # Calculate affinities for formation of aqueous species
+ species(c("SO4-2", "HS-"))
+ a <- affinity(O2 = c(-55, -40), T = 125)
+ s_old <- solubility(a, in.terms.of = "SO4-2")
+
+ expect_identical(s1, s_old)
+})
More information about the CHNOSZ-commits
mailing list