[CHNOSZ-commits] r335 - in pkg/CHNOSZ: . R demo inst man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Oct 30 14:57:03 CET 2018
Author: jedick
Date: 2018-10-30 14:57:03 +0100 (Tue, 30 Oct 2018)
New Revision: 335
Added:
pkg/CHNOSZ/demo/oldsolub.R
pkg/CHNOSZ/demo/solubility.R
Removed:
pkg/CHNOSZ/demo/solubility.R
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/R/examples.R
pkg/CHNOSZ/demo/00Index
pkg/CHNOSZ/inst/NEWS
pkg/CHNOSZ/man/examples.Rd
Log:
revise demo/solubility.R for vectorized calculation
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2018-10-30 02:55:56 UTC (rev 334)
+++ pkg/CHNOSZ/DESCRIPTION 2018-10-30 13:57:03 UTC (rev 335)
@@ -1,6 +1,6 @@
Date: 2018-10-30
Package: CHNOSZ
-Version: 1.1.3-42
+Version: 1.1.3-43
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/R/examples.R
===================================================================
--- pkg/CHNOSZ/R/examples.R 2018-10-30 02:55:56 UTC (rev 334)
+++ pkg/CHNOSZ/R/examples.R 2018-10-30 13:57:03 UTC (rev 335)
@@ -28,7 +28,7 @@
demos <- function(which=c("sources", "protein.equil", "affinity", "NaCl", "density",
"ORP", "revisit", "findit", "ionize", "buffer", "protbuff", "yeastgfp", "mosaic",
- "copper", "solubility", "wjd", "bugstab", "Shh", "activity_ratios",
+ "copper", "oldsolub", "solubility", "wjd", "bugstab", "Shh", "activity_ratios",
"adenine", "DEW", "lambda", "TCA", "go-IU", "bison"), save.png=FALSE) {
# run one or more demos from CHNOSZ with ask=FALSE, and return the value of the last one
for(i in 1:length(which)) {
Modified: pkg/CHNOSZ/demo/00Index
===================================================================
--- pkg/CHNOSZ/demo/00Index 2018-10-30 02:55:56 UTC (rev 334)
+++ pkg/CHNOSZ/demo/00Index 2018-10-30 13:57:03 UTC (rev 335)
@@ -12,7 +12,8 @@
yeastgfp Subcellular locations: log fO2 - log aH2O and log a - log fO2 diagrams
mosaic Eh-pH diagram for iron oxides, sulfides and carbonate with two sets of changing basis species
copper Another example of mosaic(): complexation of copper with glycine species
-solubility Solubility of calcite or CO2(gas) as a function of pH
+oldsolub Old solubility calculations using uniroot()
+solubility Solubility of calcite and CO2(gas) as a function of pH
wjd Gibbs energy minimization: prebiological atmospheres and cell periphery of yeast
dehydration log K of dehydration reactions; SVG file contains tooltips and links
bugstab Formation potential of microbial proteins in colorectal cancer
Copied: pkg/CHNOSZ/demo/oldsolub.R (from rev 334, pkg/CHNOSZ/demo/solubility.R)
===================================================================
--- pkg/CHNOSZ/demo/oldsolub.R (rev 0)
+++ pkg/CHNOSZ/demo/oldsolub.R 2018-10-30 13:57:03 UTC (rev 335)
@@ -0,0 +1,87 @@
+# CHNOSZ/demo/oldsolub.R
+# 20150306 added to CHNOSZ as solubility.R
+# 20181030 renamed to oldsolub.R
+
+# demo showing how to calculate CO2(gas) or calcite solubility and aqueous carbonate speciation
+# the affinity() ... equilibrate() sequence in CHNOSZ gives *metastable equilibrium activities*
+# (activities for a total activity of the balanced component given in loga.balance) ...
+# here we are interested in finding the value of loga.balance itself
+# (the total activity of [CO2, HCO3-, CO3-2] species in the aqueous phase).
+# this total activity is the solubility of CO2(gas) or calcite if the affinities of the
+# aqueous species as formed from CO2(gas) or calcite are all equal to zero.
+# note that the affinities for species in metastable equilibrium are all equal.
+# Afun() calculates the metastable equilibrium affinities for a given loga.balance
+# and uniroot() finds the loga.balance where they are zero
+# additionally, if we are reacting calcite, the activity of Ca+2 should be set equal to loga.balance
+
+# 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
+
+# set this to CO2 or calcite
+what <- "calcite"
+#what <- "CO2"
+
+# function to return the affinity of the metastable equilibrium species
+Afun <- function(loga.balance=-3, T=25) {
+ if(what=="calcite") basis("Ca+2", loga.balance)
+ a <- affinity(T=T)
+ e <- equilibrate(a, loga.balance=loga.balance)
+ # set metastable activities and re-calculate the affinity
+ species(1:3, unlist(e$loga.equil))
+ a <- affinity(T=T)
+ # check they're actually equal
+ stopifnot(all(abs(unlist(a$values) - as.vector(a$values[[1]])) < 1e-10))
+ return(a$values[[1]])
+}
+
+# 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"))
+T <- 25
+# decrease this for higher resolution
+pHstep <- 1
+
+# where we'll store the results
+loga.tot <- numeric()
+loga.CO2 <- loga.HCO3 <- loga.CO3 <- numeric()
+
+# loop over pH range
+pHs <- seq(0, 14, pHstep)
+for(pH in pHs) {
+ print(paste("pH =", pH))
+ basis("pH", pH)
+ # this is for the solubility
+ loga.balance <- suppressMessages(uniroot(Afun, c(-10, 10), T=T)$root)
+ loga.tot <- c(loga.tot, loga.balance)
+ # this is for the speciation
+ if(what=="calcite") basis("Ca+2", loga.balance)
+ a <- affinity(T=T)
+ e <- equilibrate(a, loga.balance=loga.balance)
+ loga <- unlist(e$loga.equil)
+ loga.CO2 <- c(loga.CO2, loga[1])
+ loga.HCO3 <- c(loga.HCO3, loga[2])
+ loga.CO3 <- c(loga.CO3, loga[3])
+}
+
+# make plot
+ylim <- c(-10, 4)
+thermo.plot.new(xlim=range(pHs), ylim=ylim, xlab="pH", ylab="log a")
+lines(pHs, loga.tot, lwd=4, col="green2")
+lines(pHs, loga.CO2, lwd=2)
+lines(pHs, loga.HCO3, lty=2, lwd=2)
+lines(pHs, loga.CO3, 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)))
Deleted: pkg/CHNOSZ/demo/solubility.R
===================================================================
--- pkg/CHNOSZ/demo/solubility.R 2018-10-30 02:55:56 UTC (rev 334)
+++ pkg/CHNOSZ/demo/solubility.R 2018-10-30 13:57:03 UTC (rev 335)
@@ -1,83 +0,0 @@
-# demo showing how to calculate CO2(gas) or calcite solubility and aqueous carbonate speciation
-# the affinity() ... equilibrate() sequence in CHNOSZ gives *metastable equilibrium activities*
-# (activities for a total activity of the balanced component given in loga.balance) ...
-# here we are interested in finding the value of loga.balance itself
-# (the total activity of [CO2, HCO3-, CO3-2] species in the aqueous phase).
-# this total activity is the solubility of CO2(gas) or calcite if the affinities of the
-# aqueous species as formed from CO2(gas) or calcite are all equal to zero.
-# note that the affinities for species in metastable equilibrium are all equal.
-# Afun() calculates the metastable equilibrium affinities for a given loga.balance
-# and uniroot() finds the loga.balance where they are zero
-# additionally, if we are reacting calcite, the activity of Ca+2 should be set equal to loga.balance
-
-# 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
-
-# set this to CO2 or calcite
-what <- "calcite"
-#what <- "CO2"
-
-# function to return the affinity of the metastable equilibrium species
-Afun <- function(loga.balance=-3, T=25) {
- if(what=="calcite") basis("Ca+2", loga.balance)
- a <- affinity(T=T)
- e <- equilibrate(a, loga.balance=loga.balance)
- # set metastable activities and re-calculate the affinity
- species(1:3, unlist(e$loga.equil))
- a <- affinity(T=T)
- # check they're actually equal
- stopifnot(all(abs(unlist(a$values) - as.vector(a$values[[1]])) < 1e-10))
- return(a$values[[1]])
-}
-
-# 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"))
-T <- 25
-# decrease this for higher resolution
-pHstep <- 1
-
-# where we'll store the results
-loga.tot <- numeric()
-loga.CO2 <- loga.HCO3 <- loga.CO3 <- numeric()
-
-# loop over pH range
-pHs <- seq(0, 14, pHstep)
-for(pH in pHs) {
- print(paste("pH =", pH))
- basis("pH", pH)
- # this is for the solubility
- loga.balance <- suppressMessages(uniroot(Afun, c(-10, 10), T=T)$root)
- loga.tot <- c(loga.tot, loga.balance)
- # this is for the speciation
- if(what=="calcite") basis("Ca+2", loga.balance)
- a <- affinity(T=T)
- e <- equilibrate(a, loga.balance=loga.balance)
- loga <- unlist(e$loga.equil)
- loga.CO2 <- c(loga.CO2, loga[1])
- loga.HCO3 <- c(loga.HCO3, loga[2])
- loga.CO3 <- c(loga.CO3, loga[3])
-}
-
-# make plot
-ylim <- c(-10, 4)
-thermo.plot.new(xlim=range(pHs), ylim=ylim, xlab="pH", ylab="log a")
-lines(pHs, loga.tot, lwd=4, col="green2")
-lines(pHs, loga.CO2, lwd=2)
-lines(pHs, loga.HCO3, lty=2, lwd=2)
-lines(pHs, loga.CO3, 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)))
Added: pkg/CHNOSZ/demo/solubility.R
===================================================================
--- pkg/CHNOSZ/demo/solubility.R (rev 0)
+++ pkg/CHNOSZ/demo/solubility.R 2018-10-30 13:57:03 UTC (rev 335)
@@ -0,0 +1,86 @@
+# CHNOSZ/demo/solubility.R: vectorized solubility without uniroot
+# adapted from CHNOSZ/demo/oldsolub.R
+# 20181030 jmd
+
+# 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(what in c("CO2", "calcite")) {
+
+ # 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 range and resolution, constant temperature and ionic strength
+ pH <- c(0, 14)
+ res <- 100
+ T <- 25
+ IS <- 0
+
+ # 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)
+
+ # 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]]
+
+ # 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 <- 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)))
+
+}
Modified: pkg/CHNOSZ/inst/NEWS
===================================================================
--- pkg/CHNOSZ/inst/NEWS 2018-10-30 02:55:56 UTC (rev 334)
+++ pkg/CHNOSZ/inst/NEWS 2018-10-30 13:57:03 UTC (rev 335)
@@ -1,4 +1,4 @@
-CHANGES IN CHNOSZ 1.1.3-42 (2018-10-30)
+CHANGES IN CHNOSZ 1.1.3-43 (2018-10-30)
---------------------------------------
THERMODYNAMIC DATA
@@ -135,6 +135,10 @@
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/examples.Rd
===================================================================
--- pkg/CHNOSZ/man/examples.Rd 2018-10-30 02:55:56 UTC (rev 334)
+++ pkg/CHNOSZ/man/examples.Rd 2018-10-30 13:57:03 UTC (rev 335)
@@ -15,9 +15,9 @@
examples(save.png = FALSE)
demos(which = c("sources", "protein.equil", "affinity", "NaCl",
"density", "ORP", "revisit", "findit", "ionize", "buffer",
- "protbuff", "yeastgfp", "mosaic", "copper", "solubility",
- "wjd", "bugstab", "Shh", "activity_ratios", "adenine",
- "DEW", "lambda", "TCA", "go-IU", "bison"),
+ "protbuff", "yeastgfp", "mosaic", "copper", "oldsolub",
+ "solubility", "wjd", "bugstab", "Shh", "activity_ratios",
+ "adenine", "DEW", "lambda", "TCA", "go-IU", "bison"),
save.png=FALSE)
}
@@ -45,7 +45,8 @@
\code{yeastgfp} \tab * Subcellular locations: \logfO2 - \logaH2O and \loga - \logfO2 diagrams (Dick, 2009) \cr
\code{mosaic} \tab * Eh-pH diagram with two sets of changing basis species (Garrels and Christ, 1965) \cr
\code{copper} \tab * Another example of \code{\link{mosaic}}: complexation of Cu with glycine (Aksu and Doyle, 2001) \cr
- \code{solubility} \tab * Solubility of calcite (cf. Manning et al., 2013) or \CO2 (cf. Stumm and Morgan, 1996) \cr
+ \code{oldsolub} \tab Old solubility calculations using \code{\link{uniroot}} \cr
+ \code{solubility} \tab * Solubility of calcite (cf. Manning et al., 2013) and \CO2 (cf. Stumm and Morgan, 1996) \cr
\code{wjd} \tab * \eqn{G}{G} minimization: prebiological atmospheres (Dayhoff et al., 1964) and cell periphery of yeast \cr
\code{dehydration} \tab * \logK of dehydration reactions; SVG file contains tooltips and links \cr
\code{bugstab} \tab * Formation potential of microbial proteins in colorectal cancer (Dick, 2016) \cr
More information about the CHNOSZ-commits
mailing list