[CHNOSZ-commits] r337 - in pkg/CHNOSZ: . R demo inst man tests/testthat vignettes
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Oct 31 07:25:36 CET 2018
Author: jedick
Date: 2018-10-31 07:25:35 +0100 (Wed, 31 Oct 2018)
New Revision: 337
Added:
pkg/CHNOSZ/R/solubility.R
pkg/CHNOSZ/man/solubility.Rd
pkg/CHNOSZ/tests/testthat/test-solubility.R
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/NAMESPACE
pkg/CHNOSZ/R/diagram.R
pkg/CHNOSZ/demo/solubility.R
pkg/CHNOSZ/inst/NEWS
pkg/CHNOSZ/man/nonideal.Rd
pkg/CHNOSZ/vignettes/anintro.Rmd
Log:
add solubility() function, test, and expanded demo
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2018-10-31 03:29:20 UTC (rev 336)
+++ pkg/CHNOSZ/DESCRIPTION 2018-10-31 06:25:35 UTC (rev 337)
@@ -1,6 +1,6 @@
Date: 2018-10-31
Package: CHNOSZ
-Version: 1.1.3-44
+Version: 1.1.3-45
Title: Thermodynamic Calculations and Diagrams for Geo(bio)chemistry
Authors at R: c(
person("Jeffrey", "Dick", , "j3ffdick at gmail.com", role = c("aut", "cre"),
Modified: pkg/CHNOSZ/NAMESPACE
===================================================================
--- pkg/CHNOSZ/NAMESPACE 2018-10-31 03:29:20 UTC (rev 336)
+++ pkg/CHNOSZ/NAMESPACE 2018-10-31 06:25:35 UTC (rev 337)
@@ -56,7 +56,7 @@
"calculateEpsilon", "calculateQ", "water.DEW", "berman",
"maxdiff", "expect_maxdiff", "Bdot",
# added 20171121 or later
- "dumpdata", "thermo.axis"
+ "dumpdata", "thermo.axis", "solubility"
)
# Load shared objects
Modified: pkg/CHNOSZ/R/diagram.R
===================================================================
--- pkg/CHNOSZ/R/diagram.R 2018-10-31 03:29:20 UTC (rev 336)
+++ pkg/CHNOSZ/R/diagram.R 2018-10-31 06:25:35 UTC (rev 337)
@@ -47,6 +47,7 @@
# 'loga.equil' - equilibrium activities of species of interest (eout)
# name of basis species - equilibrium activity of a basis species (aout)
# 'saturation' - affinity=0 line for each species (2D)
+ # 'loga.balance' - activity of balanced basis species (eout from solubility())
eout.is.aout <- FALSE
plot.loga.basis <- FALSE
if(type %in% c("auto", "saturation")) {
@@ -62,7 +63,7 @@
if(isTRUE(alpha) | is.character(alpha)) stop("equilibrium activities of basis species not available with alpha=TRUE")
plot.loga.basis <- TRUE
} else if(type=="loga.equil" & !"loga.equil" %in% names(eout)) stop("'eout' is not the output from equil()")
- else if(type!="loga.equil") stop(type, " is not a valid diagram type")
+ else if(!type %in% c("loga.equil", "loga.balance")) stop(type, " is not a valid diagram type")
## consider a different number of species if we're grouping them together
ngroups <- length(groups)
@@ -72,6 +73,12 @@
plotvals <- eout$loga.equil
plotvar <- "loga.equil"
+ ## handle loga.balance here (i.e. solubility calculations)
+ if(type=="loga.balance") {
+ plotvals <- list(eout$loga.balance)
+ plotvar <- "loga.balance"
+ }
+
## number of dimensions (T, P or chemical potentials that are varied)
# length(eout$vars) - the number of variables = the maximum number of dimensions
# length(dim(eout$values[[1]])) - nd=1 if it was a transect along multiple variables
@@ -209,7 +216,10 @@
### general plot parameters ###
## handle line type/width/color arguments
- if(is.null(lty)) lty <- 1:ngroups
+ if(is.null(lty)) {
+ if(type=="loga.balance") lty <- 1
+ else lty <- 1:ngroups
+ }
lty <- rep(lty, length.out=ngroups)
lwd <- rep(lwd, length.out=ngroups)
col <- rep(col, length.out=ngroups)
@@ -302,7 +312,7 @@
lines(spline.x, spline.y, col=col[i], lty=lty[i], lwd=lwd[i])
}
}
- if(!add & !is.null(legend.x)) {
+ if(type %in% c("auto", "loga.equil") & !is.null(legend.x)) {
# 20120521: use legend.x=NA to label lines rather than make legend
if(is.na(legend.x)) {
maxvals <- do.call(pmax, pv)
@@ -581,7 +591,7 @@
else contour(xs, ys, zs, add=TRUE, col=col, lty=lty, lwd=lwd, labcex=cex, levels=0, labels=names[i], method=contour.method)
}
} else {
- # otherwise, make contours of properties using first species only
+ # contour solubilities (loga.balance), or properties using first species only
if(length(plotvals) > 1) warning("showing only first species in 2-D property diagram")
zs <- plotvals[[1]]
contour(xs, ys, zs, add=TRUE, col=col, lty=lty, lwd=lwd, labcex=cex, method=contour.method)
Added: pkg/CHNOSZ/R/solubility.R
===================================================================
--- pkg/CHNOSZ/R/solubility.R (rev 0)
+++ pkg/CHNOSZ/R/solubility.R 2018-10-31 06:25:35 UTC (rev 337)
@@ -0,0 +1,32 @@
+# solubility.R: vectorized solubility calculations without uniroot
+# 20181031 jmd
+
+solubility <- function(eout, exp = 1) {
+ # exp = 1: e.g. dissolution of CO2
+ # exp = 2: e.g. dissolution (dissociation) of CaCO3
+
+ # bookkeeping: track any single species
+ itrack <- 1
+ # the log activity used to calculate the affinity
+ loga.species.track <- eout$species$logact[itrack]
+ # the affinities at the starting loga.balance
+ A.track <- eout$values[[itrack]]
+ # the loga.equil at the starting loga.balance
+ loga.equil.track <- eout$loga.equil[[itrack]]
+
+ # subjunctive: what would the affinities be if the
+ # activity of the tracked species was set to loga.equil?
+ A.whatif <- loga.species.track + A.track - loga.equil.track
+
+ # predictive: assuming the species distribution doesn't change,
+ # what is the total loga that would give zero affinity?
+ # TODO: modify this according to stoichiometry (species with > 1 of the balanced basis species)
+ loga.total <- (eout$loga.balance + A.whatif) / exp
+
+ message("solubility: calculated logarithm of total activity of ", eout$balance)
+
+ # use the predicted loga.total to re-calculate activities of species
+ aout <- eout[1:which(names(eout)=="values")]
+ equilibrate(aout, loga.balance = loga.total)
+}
+
Modified: pkg/CHNOSZ/demo/solubility.R
===================================================================
--- pkg/CHNOSZ/demo/solubility.R 2018-10-31 03:29:20 UTC (rev 336)
+++ pkg/CHNOSZ/demo/solubility.R 2018-10-31 06:25:35 UTC (rev 337)
@@ -1,86 +1,69 @@
# CHNOSZ/demo/solubility.R: vectorized solubility without uniroot
-# adapted from CHNOSZ/demo/oldsolub.R
-# 20181030 jmd
+# 20181030 adapted from CHNOSZ/demo/oldsolub.R
+# 20181031 use new solubility(); add T-pH plots
-# for comparison with published calcite solubility plot, see Fig. 4A in
-# Manning et al., 2013, Reviews in Mineralogy & Geochemistry, v. 75, pp. 109-148
-# (doi: 10.2138/rmg.2013.75.5)
-
# for comparison with published CO2 solubility plot, see Fig. 4.5 in
# Stumm and Morgan, 1996, Aquatic Chemistry: Chemical Equilibria and Rates in Natural Waters
# (New York: John Wiley & Sons), 3rd edition
-par(mfrow=c(1, 2))
+# for comparison with published calcite solubility plot, see Fig. 4A in
+# Manning et al., 2013, Reviews in Mineralogy & Geochemistry, v. 75, pp. 109-148
+# (doi: 10.2138/rmg.2013.75.5)
-for(what in c("CO2", "calcite")) {
+layout(matrix(1:4, nrow = 2))
- # set up system
- if(what=="CO2") {
- basis("CHNOS+")
- basis("CO2", "gas")
- # ca. atmospheric PCO2
- basis("CO2", -3.5)
- } else if(what=="calcite") {
- basis(c("calcite", "Ca+2", "H2O", "O2", "H+"))
- }
- species(c("CO2", "HCO3-", "CO3-2"))
+# set pH and T range and resolution, constant temperature and ionic strength
+pH <- c(0, 14)
+T <- c(0, 300)
+res <- 100
+T1 <- 25
+IS <- 0
- # set pH range and resolution, constant temperature and ionic strength
- pH <- c(0, 14)
- res <- 100
- T <- 25
- IS <- 0
+# start with CO2
+basis(c("carbon dioxide", "H2O", "O2", "H+"))
+# ca. atmospheric PCO2
+basis("CO2", -3.5)
+species(c("CO2", "HCO3-", "CO3-2"))
+a <- affinity(pH = c(pH, res), T = T1, IS = IS)
+e <- equilibrate(a)
+s <- solubility(e)
+# first plot total activity line
+diagram(s, ylim = c(-10, 4), type = "loga.balance", lwd = 4, col = "green2")
+# add activities of species
+diagram(s, ylim=c(-10, 4), add = TRUE, dy = 1)
+# add legend
+lexpr <- as.expression(c("total", expr.species("CO2", state = "aq"),
+ expr.species("HCO3-"), expr.species("CO3-2")))
+legend("topleft", lty = c(1, 1:3), lwd = c(4, 2, 2, 2),
+ col = c("green2", rep("black", 3)), legend = lexpr)
+title(main = substitute("Solubility of"~what~"at"~T~degree*"C",
+ list(what = expr.species("CO2"), T = T1)), line = 1.6)
+mtext("cf. Fig. 4.5 of Stumm and Morgan, 1996")
- # start with loga.balance = 0
- loga.balance <- 0
- a0 <- affinity(pH = c(pH, res), T = T, IS = IS)
- e0 <- equilibrate(a0, loga.balance = loga.balance)
+# CO2 T-pH plot
+a <- affinity(pH = c(pH, res), T = c(T, res), IS = IS)
+e <- equilibrate(a)
+s <- solubility(e)
+diagram(s, type = "loga.balance")
+title(main = substitute("Solubility of"~what, list(what = expr.species("CO2"))))
- # bookkeeping: track any single species
- itrack <- 1
- # the log activity used to calculate the affinity
- loga.species.track <- species(itrack)$logact
- # the affinities at the starting loga.balance
- A.track <- a0$values[[itrack]]
- # the loga.equil at the starting loga.balance
- loga.equil.track <- e0$loga.equil[[itrack]]
+# now do calcite
+basis(c("calcite", "Ca+2", "H2O", "O2", "H+"))
+species(c("CO2", "HCO3-", "CO3-2"))
+a <- affinity(pH = c(pH, res), T = T1, IS = IS)
+e <- equilibrate(a)
+s <- solubility(e, exp = 2)
+diagram(s, ylim = c(-10, 4), type = "loga.balance", lwd = 4, col = "green2")
+diagram(s, add = TRUE, dy = 1)
+legend("topright", lty = c(1, 1:3), lwd = c(4, 2, 2, 2),
+ col = c("green2", rep("black", 3)), legend = lexpr)
+title(main = substitute("Solubility of"~what~"at"~T~degree*"C",
+ list(what = "calcite", T = T1)), line = 1.6)
+mtext("cf. Fig. 4A of Manning et al., 2013")
- # subjunctive: what would the affinities be if the
- # activity of the tracked species was set to loga.equil?
- A.whatif <- loga.species.track + A.track - loga.equil.track
-
- # predictive: assuming the species distribution doesn't change,
- # what is the log(total activity) that gives zero affinity?
- # TODO: modify this according to stoichiometry (species with > 1 of the balanced basis species)
- loga.total <- loga.balance + A.whatif
-
- # a dissociation reaction makes two things, so the exponent is not 1
- if(what=="calcite") loga.total <- loga.total / 2
-
- # use the predicted loga.total to re-calculate activities of species
- e1 <- equilibrate(a0, loga.balance = loga.total)
-
-# # check that we got stable conditions
-# # (not easily vectorized because we change the activities of species)
-# print("checking for stable conditions (affinity = 0)")
-# for(i in 1:length(a0$vals[[1]])) {
-# basis(a0$vars[1], a0$vals[[1]][i])
-# if(what=="calcite") basis("Ca+2", loga.total[i])
-# species(1:3, sapply(e1$loga.equil, "[", i))
-# atest <- suppressMessages(affinity(T = T, IS = IS))
-# stopifnot(all(sapply(as.numeric(unlist(atest$values)), all.equal, 0)))
-# }
-
- # make plot
- ylim <- c(-10, 4)
- thermo.plot.new(xlim = range(pH), ylim = ylim, xlab = "pH", ylab = "log a")
- lines(a0$vals[[1]], loga.total, lwd = 4, col = "green2")
- lines(a0$vals[[1]], e1$loga.equil[[1]], lwd = 2)
- lines(a0$vals[[1]], e1$loga.equil[[2]], lty = 2, lwd = 2)
- lines(a0$vals[[1]], e1$loga.equil[[3]], lty = 3, lwd = 2)
-
- legend(ifelse(what=="calcite", "topright", "topleft"), lty = c(1, 1:3), lwd = c(4, 2, 2, 2), col = c("green2", rep("black", 3)),
- legend = as.expression(c("total", expr.species("CO2", state = "aq"), expr.species("HCO3-"), expr.species("CO3-2"))))
- title(main = substitute("Solubility of"~what~"at"~T~degree*"C", list(what = what, T = T)))
-
-}
+# calcite T-pH plot
+a <- affinity(pH = c(pH, res), T = c(T, res), IS = IS)
+e <- equilibrate(a)
+s <- solubility(e, exp = 2)
+diagram(s, type = "loga.balance")
+title(main = "Solubility of calcite", font.main = 1)
Modified: pkg/CHNOSZ/inst/NEWS
===================================================================
--- pkg/CHNOSZ/inst/NEWS 2018-10-31 03:29:20 UTC (rev 336)
+++ pkg/CHNOSZ/inst/NEWS 2018-10-31 06:25:35 UTC (rev 337)
@@ -1,6 +1,17 @@
-CHANGES IN CHNOSZ 1.1.3-43 (2018-10-30)
+CHANGES IN CHNOSZ 1.1.3-45 (2018-10-31)
---------------------------------------
+NEW FEATURES
+
+- Add solubility(). Run this after equilibrate() to calculate the
+ solubility (loga.balance) of the balanced basis species.
+
+- Revise demo/solubility.R to show solubility calculations for CO2(gas)
+ and calcite as a function of T and pH.
+
+- The old solubility demo, which uses uniroot() instead of the
+ vectorized calculations in solubility(), has been renamed oldsolub.R.
+
THERMODYNAMIC DATA
- The Berman data (Berman, 1988 and later additions) have replaced the
@@ -135,10 +146,6 @@
duplicated references (but not including any secondary references in
thermo$obigt$ref2). Thanks to Evgeniy Bastrakov for the suggestion.
-- Revise demo/solubility.R to perform vectorized, direct solubility
- calculations. The old demo that uses uniroot() has been moved to
- demo/oldsolub.R.
-
CHANGES IN CHNOSZ 1.1.3 (2017-11-13)
------------------------------------
Modified: pkg/CHNOSZ/man/nonideal.Rd
===================================================================
--- pkg/CHNOSZ/man/nonideal.Rd 2018-10-31 03:29:20 UTC (rev 336)
+++ pkg/CHNOSZ/man/nonideal.Rd 2018-10-31 06:25:35 UTC (rev 337)
@@ -117,7 +117,7 @@
a <- affinity(IS=c(0, 0.14), T=T)
e <- equilibrate(a)
if(T==25) diagram(e, ylim=c(-3.0, -2.6), legend.x=NULL)
- else diagram(e, ylim=c(-3.0, -2.6), add=TRUE, col="red")
+ else diagram(e, add=TRUE, names=NULL, col="red")
}
title(main="Non-ideality model for phosphate species")
dp <- describe.property(c("pH", "T", "T"), c(7, Ts))
Added: pkg/CHNOSZ/man/solubility.Rd
===================================================================
--- pkg/CHNOSZ/man/solubility.Rd (rev 0)
+++ pkg/CHNOSZ/man/solubility.Rd 2018-10-31 06:25:35 UTC (rev 337)
@@ -0,0 +1,86 @@
+\encoding{UTF-8}
+\name{solubility}
+\alias{solubility}
+\title{Equilibrium Chemical Activities of Species}
+\description{
+Calculate chemical activities of species in equilibrium with a soluble basis species.
+}
+
+\usage{
+ solubility(eout, exp = 1)
+}
+
+\arguments{
+ \item{eout}{list, output from \code{\link{equilibrate}}}
+ \item{exp}{numeric, exponent characterizing the stoichiometry of the dissociation reaction}
+}
+
+\details{
+Use this function to calculate solubilities of minerals (such as CaCO\s3) or gases (such as CO\s2).
+Start by using \code{\link{equilibrate}} to calculate equilibrium chemical activities of species given a constant value of \code{loga.balance} (the logarithm of total activity of the balanced basis species).
+Note that this produces affinities of formation reactions of species that are equal to each other, but are generally not equal to zero.
+\code{solubility} adjusts \code{loga.balance} such that the affinities of the formation reaction become zero.
+This corresponds to \dQuote{true} equilibrium for a solution in contact with the balanced basis species - i.e. the solubility of that species.
+
+Normally, the balance is automatically identified as the first basis species that is present in all of the species.
+If that is not adequate, it can be explicitly set via the \samp{balance} setting in \code{equilibrate}.
+
+The value of \code{exp} should be changed when calculating solubility of species that dissociate (not just dissolve).
+For example, set \code{exp} to 2 for calculating the solubility of calcite (CaCO\s3).
+
+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"}.
+}
+
+\examples{\dontshow{data(thermo)}
+# solubility of CO2 and calcite as a function of pH
+par(mfrow = c(1, 2))
+
+# set pH range and resolution, constant temperature and ionic strength
+pH <- c(0, 14)
+res <- 100
+T <- 25
+IS <- 0
+
+# start with CO2
+basis(c("carbon dioxide", "H2O", "O2", "H+"))
+# ca. atmospheric PCO2
+basis("CO2", -3.5)
+species(c("CO2", "HCO3-", "CO3-2"))
+a <- affinity(pH = c(pH, res), T = T, IS = IS)
+e <- equilibrate(a)
+s <- solubility(e)
+# first plot total activity line
+diagram(s, ylim = c(-10, 4), type = "loga.balance", lwd = 4, col = "green2")
+# add activities of species
+diagram(s, ylim=c(-10, 4), add = TRUE, dy = 1)
+# add legend
+lexpr <- as.expression(c("total", expr.species("CO2", state = "aq"),
+ expr.species("HCO3-"), expr.species("CO3-2")))
+legend("topleft", lty = c(1, 1:3), lwd = c(4, 2, 2, 2),
+ col = c("green2", rep("black", 3)), legend = lexpr)
+title(main = substitute("Solubility of"~what~"at"~T~degree*"C",
+ list(what = expr.species("CO2"), T = T)), line = 1.5)
+mtext("cf. Fig. 4.5 of Stumm and Morgan, 1996")
+
+# now do calcite
+basis(c("calcite", "Ca+2", "H2O", "O2", "H+"))
+species(c("CO2", "HCO3-", "CO3-2"))
+a <- affinity(pH = c(pH, res), T = T, IS = IS)
+e <- equilibrate(a)
+s <- solubility(e, exp = 2)
+diagram(s, ylim = c(-10, 4), type = "loga.balance", lwd = 4, col = "green2")
+diagram(s, add = TRUE, dy = 1)
+legend("topright", lty = c(1, 1:3), lwd = c(4, 2, 2, 2),
+ col = c("green2", rep("black", 3)), legend = lexpr)
+title(main = substitute("Solubility of"~what~"at"~T~degree*"C",
+ list(what = "calcite", T = T)))
+mtext("cf. Fig. 4A of Manning et al., 2013")
+}
+
+\references{
+Manning, C. E., Shock, E. L. and Sverjensky, D. A. (2013) The chemistry of carbon in aqueous fluids at crustal and upper-mantle conditions: Experimental and theoretical constraints. \emph{Rev. Mineral. Geochem.} \bold{75}, 109--148. \url{https://doi.org/10.2138/rmg.2013.75.5}
+
+Stumm, W. and Morgan, J. J. (1996) \emph{Aquatic Chemistry: Chemical Equilibria and Rates in Natural Waters}, John Wiley & Sons, New York, 1040 p. \url{http://www.worldcat.org/oclc/31754493}
+}
+
+\concept{Main workflow}
Added: pkg/CHNOSZ/tests/testthat/test-solubility.R
===================================================================
--- pkg/CHNOSZ/tests/testthat/test-solubility.R (rev 0)
+++ pkg/CHNOSZ/tests/testthat/test-solubility.R 2018-10-31 06:25:35 UTC (rev 337)
@@ -0,0 +1,35 @@
+context("solubility")
+
+test_that("solubility() produces stable conditions (affinity = 0)", {
+ # set pH range and resolution, constant temperature and ionic strength
+ pH <- c(0, 14)
+ res <- 100
+ T <- 25
+ IS <- 0
+
+ # start with CO2
+ basis(c("carbon dioxide", "H2O", "O2", "H+"))
+ # ca. atmospheric PCO2
+ basis("CO2", -3.5)
+ species(c("CO2", "HCO3-", "CO3-2"))
+ a <- affinity(pH = c(pH, res), T = T, IS = IS)
+ e <- equilibrate(a)
+ s <- solubility(e)
+
+ # check for stable conditions (affinity = 0)
+ species(1:3, 0)
+ atest <- affinity(pH = s$vals[[1]], T = T, IS = IS)
+ expect_true(all(sapply(unlist(atest$values) - unlist(s$loga.equil), all.equal, 0)))
+
+ # now do calcite
+ basis(c("calcite", "Ca+2", "H2O", "O2", "H+"))
+ species(c("CO2", "HCO3-", "CO3-2"))
+ a <- affinity(pH = c(pH, res), T = T, IS = IS)
+ e <- equilibrate(a)
+ s <- solubility(e, exp = 2)
+
+ # check for stable conditions (affinity = 0)
+ species(1:3, 0)
+ atest <- affinity(pH = s$vals[[1]], `Ca+2` = s$loga.balance, T = T, IS = IS)
+ expect_true(all(sapply(unlist(atest$values) - unlist(s$loga.equil), all.equal, 0)))
+})
Modified: pkg/CHNOSZ/vignettes/anintro.Rmd
===================================================================
--- pkg/CHNOSZ/vignettes/anintro.Rmd 2018-10-31 03:29:20 UTC (rev 336)
+++ pkg/CHNOSZ/vignettes/anintro.Rmd 2018-10-31 06:25:35 UTC (rev 337)
@@ -895,13 +895,13 @@
a25 <- affinity(pH = c(4, 13))
a150 <- affinity(pH = c(4, 13), T = 150)
diagram(a25, dy = 0.4)
-diagram(a150, add = TRUE, col = "red")
+diagram(a150, add = TRUE, names = NULL, col = "red")
e25 <- equilibrate(a25, loga.balance = -3)
e150 <- equilibrate(a150, loga.balance = -3)
diagram(e25, ylim = c(-6, 0), dy = 0.15)
-diagram(e150, add = TRUE, col = "red")
+diagram(e150, add = TRUE, names = NULL, col = "red")
diagram(e25, alpha = TRUE, dy = -0.25)
-diagram(e150, alpha = TRUE, add = TRUE, col = "red")
+diagram(e150, alpha = TRUE, add = TRUE, names = NULL, col = "red")
```
The distribution of aqueous carbonate species as a function of pH (a type of Bjerrum plot) is a classic example of an equilibrium calculation.
More information about the CHNOSZ-commits
mailing list