[CHNOSZ-commits] r143 - in pkg/CHNOSZ: . inst inst/extdata/cpetc man vignettes
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sat Feb 11 19:37:23 CET 2017
Author: jedick
Date: 2017-02-11 19:37:23 +0100 (Sat, 11 Feb 2017)
New Revision: 143
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/inst/NEWS
pkg/CHNOSZ/inst/extdata/cpetc/rubisco.csv
pkg/CHNOSZ/man/extdata.Rd
pkg/CHNOSZ/vignettes/anintro.Rmd
Log:
anintro.Rmd: use RSVGTipsDevice for Rubisco plot
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2017-02-11 12:23:32 UTC (rev 142)
+++ pkg/CHNOSZ/DESCRIPTION 2017-02-11 18:37:23 UTC (rev 143)
@@ -1,11 +1,11 @@
Date: 2017-02-11
Package: CHNOSZ
-Version: 1.0.8-32
+Version: 1.0.8-33
Title: Chemical Thermodynamics and Activity Diagrams
Author: Jeffrey Dick
Maintainer: Jeffrey Dick <j3ffdick at gmail.com>
Depends: R (>= 3.1.0)
-Suggests: limSolve, testthat, knitr, rmarkdown, tufte
+Suggests: limSolve, testthat, knitr, rmarkdown, tufte, RSVGTipsDevice
Imports: grDevices, graphics, stats, utils, colorspace
Description: Functions and data sets to support chemical thermodynamic modeling in biochemistry
and low-temperature geochemistry. The features include calculation of the standard molal
Modified: pkg/CHNOSZ/inst/NEWS
===================================================================
--- pkg/CHNOSZ/inst/NEWS 2017-02-11 12:23:32 UTC (rev 142)
+++ pkg/CHNOSZ/inst/NEWS 2017-02-11 18:37:23 UTC (rev 143)
@@ -1,4 +1,4 @@
-CHANGES IN CHNOSZ 1.0.8-32 (2017-02-11)
+CHANGES IN CHNOSZ 1.0.8-33 (2017-02-11)
---------------------------------------
DOCUMENTATION:
Modified: pkg/CHNOSZ/inst/extdata/cpetc/rubisco.csv
===================================================================
--- pkg/CHNOSZ/inst/extdata/cpetc/rubisco.csv 2017-02-11 12:23:32 UTC (rev 142)
+++ pkg/CHNOSZ/inst/extdata/cpetc/rubisco.csv 2017-02-11 18:37:23 UTC (rev 143)
@@ -1,28 +1,28 @@
-ID,T1,T2,domain,species
-G9FID8,0,6,E,"Phaeocystis antarctica"
-M9R7V1,4,10,B,"Octadecabacter antarcticus 307"
-K4MAK9,18,18,A,"Methanolobus psychrophilus"
-Q12TQ0,23,23,A,"Methanococcoides burtonii"
-Q9ZI34,25,30,B,"Bradyrhizobium japonicum"
-P0C916,28,30,B,"Thiobacillus ferrooxidans"
-P00874,30,30,E,"Zea mays"
-Q0EX22,30,30,B,"Mariprofundus ferrooxydans"
-L0RHZ1,35,35,B,"Desulfovibrio hydrothermalis"
-Q8THG2,35,40,A,"Methanosarcina acetivorans"
-F9ZLP0,45,45,B,"Acidithiobacillus caldus"
-P37393,45,45,E,"Cyanidium caldarium"
-P72383,45,50,B,"Sulfobacillus acidophilus"
-Q51856,52,52,B,"Pseudomonas hydrogenothermophila"
-Q2JIP3,50,55,B,"Synechococcus sp. (strain JA-2-3B'a(2-13))"
-A0B9K9,55,60,A,"Methanosaeta thermophila"
-Q8DIS5,57,57,B,"Thermosynechococcus elongatus"
-G8LZL2,60,60,B,"Clostridium clariflavum"
-F8IID7,65,65,B,"Bacillus acidocaldarius"
-A8F7V4,65,65,B,"Thermotoga lettingae"
-B9KXE5,70,70,B,"Thermomicrobium roseum"
-O28635,76,76,A,"Archaeoglobus fulgidus"
-Q58632,85,85,A,"Methanocaldococcus jannaschii"
-A1RZJ5,85,90,A,"Thermofilum pendens"
-A3DND9,85,92,A,"Staphylothermus marinus"
-O58677,98,98,A,"Pyrococcus horikoshii"
-Q8U1P9,100,100,A,"Pyrococcus furiosus"
+ID,T1,T2,domain,species,URL
+G9FID8,0,6,E,"Phaeocystis antarctica",http://dx.doi.org/10.3354/ab00256
+M9R7V1,4,10,B,"Octadecabacter antarcticus 307",http://dx.doi.org/10.1016/S0723-2020(97)80003-3
+K4MAK9,18,18,A,"Methanolobus psychrophilus",http://dx.doi.org/10.1128/AEM.01146-08
+Q12TQ0,23,23,A,"Methanococcoides burtonii",http://dx.doi.org/10.1016/S0723-2020(11)80117-7
+Q9ZI34,25,30,B,"Bradyrhizobium japonicum",http://dx.doi.org/10.1099/00207713-32-1-136
+P0C916,28,30,B,"Thiobacillus ferrooxidans",http://dx.doi.org/10.1080/01490459209377920
+P00874,30,30,E,"Zea mays",http://dx.doi.org/10.1016/S0065-2113(08)60322-3
+Q0EX22,30,30,B,"Mariprofundus ferrooxydans",http://dx.doi.org/10.1371/journal.pone.0000667
+L0RHZ1,35,35,B,"Desulfovibrio hydrothermalis",http://dx.doi.org/10.1099/ijs.0.02323-0
+Q8THG2,35,40,A,"Methanosarcina acetivorans",http://aem.asm.org/content/47/5/971
+F9ZLP0,45,45,B,"Acidithiobacillus caldus",http://dx.doi.org/10.1099/13500872-140-12-3451
+P37393,45,45,E,"Cyanidium caldarium",http://dx.doi.org/10.1007/BF00409031
+P72383,45,50,B,"Sulfobacillus acidophilus",http://dx.doi.org/10.1099/00221287-142-4-775
+Q51856,52,52,B,"Pseudomonas hydrogenothermophila",http://dx.doi.org/10.1271/bbb1961.41.685
+Q2JIP3,50,55,B,"Synechococcus sp. (strain JA-2-3B'a(2-13))",http://dx.doi.org/10.1128/AEM.72.1.544-550.2006
+A0B9K9,55,60,A,"Methanosaeta thermophila",http://dx.doi.org/10.1099/00207713-42-3-463
+Q8DIS5,57,57,B,"Thermosynechococcus elongatus",http://dx.doi.org/10.1093/oxfordjournals.pcp.a075684
+G8LZL2,60,60,B,"Clostridium clariflavum",http://dx.doi.org/10.1099/ijs.0.003483-0
+F8IID7,65,65,B,"Bacillus acidocaldarius",http://dx.doi.org/10.1099/00221287-67-1-9
+A8F7V4,65,65,B,"Thermotoga lettingae",http://dx.doi.org/10.1099/ijs.0.02165-0
+B9KXE5,70,70,B,"Thermomicrobium roseum",http://dx.doi.org/10.1371/journal.pone.0004207
+O28635,76,76,A,"Archaeoglobus fulgidus",http://aem.asm.org/content/60/4/1227
+Q58632,85,85,A,"Methanocaldococcus jannaschii",http://dx.doi.org/10.1007/BF00425213
+A1RZJ5,85,90,A,"Thermofilum pendens",http://dx.doi.org/10.1016/S0723-2020(83)80035-6
+A3DND9,85,92,A,"Staphylothermus marinus",http://dx.doi.org/10.1016/S0723-2020(86)80157-6
+O58677,98,98,A,"Pyrococcus horikoshii",http://dx.doi.org/10.1007/s007920050051
+Q8U1P9,100,100,A,"Pyrococcus furiosus",http://dx.doi.org/10.1007/BF00413027
Modified: pkg/CHNOSZ/man/extdata.Rd
===================================================================
--- pkg/CHNOSZ/man/extdata.Rd 2017-02-11 12:23:32 UTC (rev 142)
+++ pkg/CHNOSZ/man/extdata.Rd 2017-02-11 18:37:23 UTC (rev 143)
@@ -36,7 +36,7 @@
\item \code{BKM60_Fig7.dat} Eh-pH values for normal, wet and waterlogged soils from Fig. 7 of Baas Becking et al., 1960. See the \sQuote{anintro} vignette for an example that uses this file.
\item \code{SC10_Rainbow.csv} Values of temperature (\eqn{^{\circ}}{°}C), pH and logarithms of activity of \eqn{\mathrm{CO_2}}{CO2}, \eqn{\mathrm{H_2}}{H2}, \eqn{\mathrm{NH_4^+}}{NH4+}, \eqn{\mathrm{H_2S}}{H2S} and \eqn{\mathrm{CH_4}}{CH4} for mixing of seawater and hydrothermal fluid at Rainbow field (Mid-Atlantic Ridge), taken from Shock and Canovas, 2010.
\item \code{SS98_Fig5a.csv}, \code{SS98_Fig5b.csv} Values of logarithm of fugacity of \eqn{\mathrm{O_2}}{O2} and pH as a function of temperature for mixing of seawater and hydrothermal fluid, digitized from Figs. 5a and b of Shock and Schulte, 1998.
- \item \code{rubisco.csv} UniProt IDs for Rubisco, ranges of optimal growth temperature of organisms, and domain and name of organisms, from Dick, 2014.
+ \item \code{rubisco.csv} UniProt IDs for Rubisco, ranges of optimal growth temperature of organisms, domain and name of organisms, and URL of reference for growth temperature, from Dick, 2014.
}
Files in \code{fasta} contain protein sequences:
Modified: pkg/CHNOSZ/vignettes/anintro.Rmd
===================================================================
--- pkg/CHNOSZ/vignettes/anintro.Rmd 2017-02-11 12:23:32 UTC (rev 142)
+++ pkg/CHNOSZ/vignettes/anintro.Rmd 2017-02-11 18:37:23 UTC (rev 143)
@@ -29,22 +29,23 @@
```
```{r setup, include=FALSE}
+library(knitr)
# invalidate cache when the tufte version changes
-knitr::opts_chunk$set(tidy = FALSE, cache.extra = packageVersion('tufte'))
+opts_chunk$set(tidy = FALSE, cache.extra = packageVersion('tufte'))
options(htmltools.dir.version = FALSE)
-knitr::knit_hooks$set(small.mar = function(before, options, envir) {
+knit_hooks$set(small.mar = function(before, options, envir) {
if (before) par(mar = c(4.2, 4.2, .1, .1)) # smaller margin on top and right
})
-knitr::knit_hooks$set(tiny.mar = function(before, options, envir) {
+knit_hooks$set(tiny.mar = function(before, options, envir) {
if (before) par(mar = c(.1, .1, .1, .1)) # tiny margin all around
})
-knitr::knit_hooks$set(smallish.mar = function(before, options, envir) {
+knit_hooks$set(smallish.mar = function(before, options, envir) {
if (before) par(mar = c(4.2, 4.2, 0.9, 0.9)) # smallish margins on top and right
})
# dpi setting
dpi <- 72
# use pngquant to optimize PNG images
-knitr::knit_hooks$set(pngquant = knitr::hook_pngquant)
+knit_hooks$set(pngquant = hook_pngquant)
pngquant <- "--speed=1 --quality=0-50"
# some frequently used HTML expressions
logfO2 <- "log<i>f</i><sub>O<sub>2</sub></sub>"
@@ -52,6 +53,16 @@
zc <- "<i>Z</i><sub>C</sub>"
o2 <- "O<sub>2</sub>"
h2o <- "H<sub>2</sub>O"
+# http://stackoverflow.com/questions/23852753/knitr-with-gridsvg
+# Set up a chunk hook for manually saved plots.
+knit_hooks$set(custom.plot = hook_plot_custom)
+# hook to change <img /> to <embed /> - required for interactive SVG
+hook_plot <- knit_hooks$get("plot")
+knit_hooks$set(plot = function(x, options) {
+ x <- hook_plot(x, options)
+ if (!is.null(options$embed.tag) && options$embed.tag) x <- gsub("<img ", "<embed ", x)
+ x
+})
```
# First steps
@@ -244,13 +255,13 @@
See also <span style="color:blue">`demo(density)`</span>.
```
```{r subcrt_water_grid}
-subcrt("water", T=c(400, 500, 600), P=c(200, 400, 600), grid="P")$out$water
+subcrt("water", T = c(400, 500, 600), P = c(200, 400, 600), grid = "P")$out$water
```
```{r subcrt_water_plot, fig.margin=TRUE, fig.width=4, fig.height=4, small.mar=TRUE, dpi=dpi, out.width="100%", echo=FALSE, message=FALSE, fig.cap="Isothermal contours of density (g cm<sup>-3</sup>) and pressure (bar) of H<sub>2</sub>O.", cache=TRUE, pngquant=pngquant}
substuff <- subcrt("water", T=seq(0,1000,100), P=c(NA, seq(1,500,1)), grid="T")
water <- substuff$out$water
-plot(water$P, water$rho, type="l")
+plot(water$P, water$rho, type = "l")
```
The additional operations (`$out$water`) are used to extract a specific part of the results; this can be used with e.g. R's `write.table()` or `plot()` for further processing:
```{r subcrt_water_plot, eval=FALSE}
@@ -265,7 +276,7 @@
T.units("K")
P.units("MPa")
E.units("J")
-subcrt("methane", T=298.15, P=0.1)$out$methane$G
+subcrt("methane", T = 298.15, P = 0.1)$out$methane$G
data(thermo) ## restore default settings
```
@@ -288,7 +299,7 @@
For an example of a solubility calculation, see <span style="color:blue">`demo(solubility)`</span>, which is based on a figure in Manning et al. (2013).
```
```{r CO2_dissolution}
-subcrt(c("CO2", "CO2"), c("gas", "aq"), c(-1, 1), T=seq(0, 250, 50))
+subcrt(c("CO2", "CO2"), c("gas", "aq"), c(-1, 1), T = seq(0, 250, 50))
```
In order to make a plot like Figure 18 of @MSS13, let's run more calculations and store the results.
@@ -296,14 +307,14 @@
```{r dissolution, echo=FALSE, message=FALSE}
T <- seq(0, 350, 10)
-CO2 <- subcrt(c("CO2", "CO2"), c("gas", "aq"), c(-1, 1), T=T)$out$logK
-CO <- subcrt(c("CO", "CO"), c("gas", "aq"), c(-1, 1), T=T)$out$logK
-CH4 <- subcrt(c("CH4", "CH4"), c("gas", "aq"), c(-1, 1), T=T)$out$logK
+CO2 <- subcrt(c("CO2", "CO2"), c("gas", "aq"), c(-1, 1), T = T)$out$logK
+CO <- subcrt(c("CO", "CO"), c("gas", "aq"), c(-1, 1), T = T)$out$logK
+CH4 <- subcrt(c("CH4", "CH4"), c("gas", "aq"), c(-1, 1), T = T)$out$logK
logK <- data.frame(T, CO2, CO, CH4)
```
```{r dissolution_plot, fig.margin=TRUE, fig.width=4, fig.height=4, small.mar=TRUE, dpi=dpi, out.width="100%", echo=FALSE, message=FALSE, fig.cap="Equilibrium constants calculated for dissolution of CO<sub>2</sub>, CO, and CH<sub>4</sub>.", cache=TRUE, pngquant=pngquant}
-matplot(logK[, 1], logK[, -1], type="l", col=1, lty=1,
- xlab=axis.label("T"), ylab=axis.label("logK"))
+matplot(logK[, 1], logK[, -1], type = "l", col = 1, lty = 1,
+ xlab = axis.label("T"), ylab = axis.label("logK"))
text(80, -1.7, expr.species("CO2"))
text(270, -2.2, expr.species("CO"))
text(310, -2.4, expr.species("CH4"))
@@ -372,13 +383,13 @@
Use <span style="color:green">`describe.reaction()`</span> to write the reactions on a plot:
```{r describe_reaction_plot, fig.margin=TRUE, fig.width=3.5, fig.height=1.8, tiny.mar=TRUE, dpi=dpi, out.width="100%", pngquant=pngquant}
-plot(0, 0, type="n", axes=FALSE, ann=FALSE, xlim=c(0, 5), ylim=c(5.2, -0.2))
-text(0, 0, "acetoclastic methanogenesis", adj=0)
-text(5, 1, describe.reaction(acetoclastic$reaction), adj=1)
-text(0, 2, "acetate oxidation", adj=0)
-text(5, 3, describe.reaction(acetate_oxidation$reaction), adj=1)
-text(0, 4, "hydrogenotrophic methanogenesis", adj=0)
-text(5, 5, describe.reaction(hydrogenotrophic$reaction), adj=1)
+plot(0, 0, type = "n", axes = FALSE, ann=FALSE, xlim=c(0, 5), ylim=c(5.2, -0.2))
+text(0, 0, "acetoclastic methanogenesis", adj = 0)
+text(5, 1, describe.reaction(acetoclastic$reaction), adj = 1)
+text(0, 2, "acetate oxidation", adj = 0)
+text(5, 3, describe.reaction(acetate_oxidation$reaction), adj = 1)
+text(0, 4, "hydrogenotrophic methanogenesis", adj = 0)
+text(5, 5, describe.reaction(hydrogenotrophic$reaction), adj = 1)
```
## Chemical affinity
@@ -405,13 +416,13 @@
```
```{r affinity_acetoclastic, message=FALSE}
subcrt(c("acetate", "methane"), c(-1, 1),
- c("aq", "gas"), logact=c(-3.4, -0.18), T=55, P=50)$out
+ c("aq", "gas"), logact = c(-3.4, -0.18), T = 55, P = 50)$out
```
The new `A` column shows the affinity; the other columns are unaffected and still show the standard-state properties.
Let's repeat the calculation for hydrogenotrophic methanogenesis.
```{r affinity_hydrogenotrophic, message=FALSE}
-subcrt("methane", 1, "gas", logact=-0.18, T=55, P=50)$out
+subcrt("methane", 1, "gas", logact = -0.18, T = 55, P = 50)$out
```
Under the specified conditions, the affinities of hydrogenotrophic and acetoclastic methanogenesis are somewhat greater than and less than 20 kJ, respectively.
@@ -425,7 +436,7 @@
```{r rxnfun, message=FALSE}
rxnfun <- function(coeffs) {
subcrt(c("acetate", "methane"), coeffs,
- c("aq", "gas"), logact=c(-3.4, -0.18), T=55, P=50)$out
+ c("aq", "gas"), logact = c(-3.4, -0.18), T = 55, P = 50)$out
}
```
@@ -446,10 +457,10 @@
)
})
Adat <- do.call(rbind, Adat)
-matplot(Adat[, 1], -Adat[, -1]/1000, type="l", lty=1, lwd=2,
- xlab=axis.label("CO2"), ylab=axis.label("DG", prefix="k"))
+matplot(Adat[, 1], -Adat[, -1]/1000, type = "l", lty = 1, lwd = 2,
+ xlab = axis.label("CO2"), ylab = axis.label("DG", prefix = "k"))
legend("topleft", c("acetate oxidation", "acetoclastic methanogenesis",
- "hydrogenotrophic methanogenesis"), lty=1, col=2:4)
+ "hydrogenotrophic methanogenesis"), lty = 1, col = 2:4)
```
```{r methanogenesis_plot, eval=FALSE}
```
@@ -491,7 +502,7 @@
The same result (in energetic units) could be obtained using <span style="color:green">`subcrt()`</span>, but <span style="color:green">`affinity()`</span> has the advantage of being able to perform calculations on a grid of *T*, *P*, or activities of basis species.
Let's choose a set of variables commonly used in aqueous speciation diagrams: Eh and pH.
-To use Eh as a variable, the electron (`e-`) should be in the basis.
+To use Eh as a variable, the electron (*e*<sup>-</sup>) should be in the basis.
To put the electron in there, we could use a different keyword (<span style="color:red">`basis("CHNOSe")`</span>), or swap oxygen out of the existing basis:
```{r swap_basis}
swap.basis("O2", "e-")
@@ -503,14 +514,14 @@
```
```{r EhpH_plot, fig.margin=TRUE, fig.width=4, fig.height=4, small.mar=TRUE, dpi=dpi, out.width="100%", echo=FALSE, message=FALSE, cache=TRUE, pngquant=pngquant}
-a <- affinity(pH=c(0, 12), Eh=c(-1, 1))
-diagram(a, fill="heat")
+a <- affinity(pH = c(0, 12), Eh = c(-1, 1))
+diagram(a, fill = "heat")
```
Now we can calculate the affinities on an Eh-pH grid:
```{r EhpH_plot, echo=1, eval=FALSE}
```
-Given values of affinity, the <span style="color:green">`diagram()`</span> function uses the maximum affinity method to make a potential diagram:
+Given values of affinity, the <span style="color:green">`diagram()`</span> function uses the maximum affinity method to make a potential diagram (i.e. a Pourbaix diagram):
```{r EhpH_plot_echo, eval=FALSE}
diagram(a)
```
@@ -520,9 +531,9 @@
In systems where equilibrium is attainable, it makes sense to call this a *predominance diagram*, showing regions of maximum activity.
```{r EhpH_plot_color, fig.margin=TRUE, fig.width=4, fig.height=4, smallish.mar=TRUE, dpi=dpi, out.width="100%", echo=FALSE, message=FALSE, cache=TRUE, pngquant=pngquant}
-diagram(a, fill="terrain", lwd=3, lty=3,
- names=c("hydrogen sulfide", "bisulfide", "bisulfate", "sulfate"),
- tplot=FALSE, main="sulfur species, 25 °C", bty="n")
+diagram(a, fill = "terrain", lwd = 3, lty = 3,
+ names = c("hydrogen sulfide", "bisulfide", "bisulfate", "sulfate"),
+ tplot = FALSE, main = "sulfur species, 25 °C", bty = "n")
```
The default colors for diagrams shown on the screen uses R's `heat.colors()` palette.
Some arguments in <span style="color:green">`diagram()`</span> can be used to control the color, labels, and lines, and title (`main`).
@@ -574,10 +585,10 @@
T <- 200
res <- 300
bases <- c("H2S", "HS-", "HSO4-", "SO4-2")
-m1 <- mosaic(bases, blend=TRUE, pH=c(0, 12, res), Eh=c(-1.2, 0.75, res), T=T)
-diagram(m1$A.species, lwd=2)
-diagram(m1$A.bases, add=TRUE, col="blue", col.names="blue", lty=2)
-water.lines("pH", "Eh", T=convert(T, "K"), col="red", lwd=2, lty=2)
+m1 <- mosaic(bases, blend = TRUE, pH = c(0, 12, res), Eh=c(-1.2, 0.75, res), T=T)
+diagram(m1$A.species, lwd = 2)
+diagram(m1$A.bases, add = TRUE, col = "blue", col.names = "blue", lty = 2)
+water.lines("pH", "Eh", T = convert(T, "K"), col = "red", lwd = 2, lty = 2)
```
The argument `blend=TRUE` is used to combine the diagrams according to the equilibrium activities of the basis species (see below).
@@ -590,19 +601,19 @@
To do that, let's write a function to swap those basis species and make a diagram.
We use R's `do.call()` to construct the argument list for <span style="color:green">`mosaic()`</span>; this way, the name of the `newvar` argument to our function indicates the chosen variable.
```{r mosaicfun, fig.fullwidth=TRUE, fig.width=9, fig.height=3, small.mar=TRUE, dpi=dpi, out.width="85%", message=FALSE, results="hide", cache=TRUE, fig.cap="Different projections (defined by the basis species) of the same thermodynamic system. The choice between them depends on convenience rather than correctness.", pngquant=pngquant}
-mosaicfun <- function(newvar, T=200, res=300) {
+mosaicfun <- function(newvar, T = 200, res = 300) {
swap.basis("e-", names(newvar))
if(names(newvar) == "O2") basis("O2", "gas")
mosaicargs <- c(list(bases), blend=TRUE, pH=list(c(-2, 12, res)), newvar, T=T)
m1 <- do.call(mosaic, mosaicargs)
- diagram(m1$A.species, lwd=2, fill=rev(topo.colors(10)))
- diagram(m1$A.bases, add=TRUE, col="blue", col.names="blue", lty=3)
+ diagram(m1$A.species, lwd = 2, fill = rev(topo.colors(10)))
+ diagram(m1$A.bases, add = TRUE, col = "blue", col.names = "blue", lty = 3)
swap.basis(names(newvar), "e-")
}
-par(mfrow=c(1, 3))
-mosaicfun(list(Eh=c(-1, 1, res)))
-mosaicfun(list(H2=c(-15, 5, res)))
-mosaicfun(list(O2=c(-70, 0, res)))
+par(mfrow = c(1, 3))
+mosaicfun(list(Eh = c(-1, 1, res)))
+mosaicfun(list(H2 = c(-15, 5, res)))
+mosaicfun(list(O2 = c(-70, 0, res)))
```
## *T*, *P*, activity transects
@@ -613,8 +624,8 @@
Let's use this feature to calculate affinities (negative Gibbs energies) of methanogenesis and biosynthetic reactions in a hydrothermal system.
Some results of mixing calculations for seawater and vent fluid from the Rainbow hydrothermal field, reported by @SC10, are included in a data file in CHNOSZ:
```{r rainbow_data}
-file <- system.file("extdata/cpetc/SC10_Rainbow.csv", package="CHNOSZ")
-rb <- read.csv(file, check.names=FALSE)
+file <- system.file("extdata/cpetc/SC10_Rainbow.csv", package = "CHNOSZ")
+rb <- read.csv(file, check.names = FALSE)
```
We take a selection of the species from Shock and Canovas (2010) with activities equal to 10<sup>-6</sup>; methane is assigned an activity of 10<sup>-3</sup>.
@@ -642,22 +653,22 @@
```
Using <span style="color:green">`convert()`</span>, we also convert the result from dimensionless values (*A*/2.303*RT*) to kcal/mol.
```{r rainbow_affinity, message=FALSE}
-a <- affinity(T=rb$T, CO2=rb$CO2, H2=rb$H2,
- `NH4+`=rb$`NH4+`, H2S=rb$H2S, pH=rb$pH)
+a <- affinity(T = rb$T, CO2 = rb$CO2, H2 = rb$H2,
+ `NH4+` = rb$`NH4+`, H2S = rb$H2S, pH = rb$pH)
T <- convert(a$vals[[1]], "K")
a$values <- lapply(a$values, convert, "G", T)
a$values <- lapply(a$values, `*`, -0.001)
```
```{r rainbow_diagram, fig.margin=TRUE, fig.width=4, fig.height=4, small.mar=TRUE, dpi=dpi, out.width="100%", echo=FALSE, message=FALSE, cache=TRUE, fig.cap="Affinities of organic synthesis in a hydrothermal system, after Shock and Canovas (2010).", pngquant=pngquant}
-diagram(a, balance=1, ylim=c(-100, 100), ylab=axis.label("A", prefix="k"),
- col=rainbow(8), lwd=2, legend.x=NA, bg="slategray3")
-abline(h=0, lty=2, lwd=2)
+diagram(a, balance = 1, ylim = c(-100, 100), ylab = axis.label("A", prefix="k"),
+ col = rainbow(8), lwd = 2, legend.x = NA, bg = "slategray3")
+abline(h = 0, lty = 2, lwd = 2)
```
Finally, we use <span style="color:green">`diagram()`</span> to plot the results.
Although only temperature is shown on the *x*-axis, pH and the activities of CO<sub>2</sub>, H<sub>2</sub>, NH<sub>4</sub><sup>+</sup>, and H<sub>2</sub>S are also varied according to the data in `rb`.
By default, <span style="color:green">`diagram()`</span> attempts to scale the affinities by dividing by the reaction coefficients of a shared basis species (in this case, CO<sub>2</sub>).
-To override that behavior, we set `balance=1` to plot the affinities of the formation reactions as written (per mole of the species being formed).
+To override that behavior, we set `balance = 1` to plot the affinities of the formation reactions as written (per mole of the species being formed).
Also, `legend.x=NA` is used to suppress making a legend (so the labels are placed next to the lines instead).
```{r rainbow_diagram, eval=FALSE}
```
@@ -693,7 +704,7 @@
The affinity of formation of pyrite happens to be zero because it is identical to one of the selected basis species.
```
```{r PPM_affinity, message=FALSE}
-unlist(affinity(T=300, P=100)$values)
+unlist(affinity(T = 300, P = 100)$values)
```
We use <span style="color:red">`mod.buffer()`</span> to choose the `cr2` phase of phyrrhotite, which is stable at this temperature.
@@ -703,11 +714,11 @@
basis(c("H2S", "O2"), c("PPM", "PPM"))
```
```{r PPM_activities, message=FALSE}
-unlist(affinity(T=300, P=100, return.buffer=TRUE)[1:3])
+unlist(affinity(T = 300, P = 100, return.buffer = TRUE)[1:3])
```
```{r demo_buffer, fig.margin=TRUE, fig.width=4, fig.height=4, small.mar=TRUE, dpi=dpi, out.width="100%", message=FALSE, echo=FALSE, cache=TRUE, pngquant=pngquant}
-demo(buffer, echo=FALSE)
+demo(buffer, echo = FALSE)
```
Et voilà! We have found log*a*<sub>H<sub>2</sub>S</sub> and `r logfO2` that are compatible with the coexistence of the three minerals.
Under these conditions, the affinities of formation reactions of the minerals in the buffer are all equal to zero:
@@ -753,49 +764,35 @@
The reaction-matrix method is slower, but can be applied to systems were the balanced basis species has non-unity coefficients.
```{r bjerrum_diagram, fig.margin=TRUE, fig.width=3, fig.height=6, small.mar=TRUE, dpi=dpi, out.width="100%", echo=FALSE, results="hide", message=FALSE, cache=TRUE, fig.cap="Three views of carbonate speciation: affinity, metastable equilibrium activity, degree of formation.", pngquant=pngquant}
-par(mfrow=c(3, 1))
+par(mfrow = c(3, 1))
basis("CHNOS+")
species(c("CO2", "HCO3-", "CO3-2"))
-# affinity at 25 °C and 150 °C
-a25 <- affinity(pH=c(4, 12))
-a150 <- affinity(pH=c(4, 12), T=150)
-diagram(a25, legend.x="topleft", bty="n")
-diagram(a150, add=TRUE, col="red")
-# equilibration
-e25 <- equilibrate(a25, loga.balance=-3)
-e150 <- equilibrate(a150, loga.balance=-3)
-diagram(e25, ylim=c(-6, 0), bty="n")
-diagram(e150, add=TRUE, col="red")
-# degree of formation
-diagram(e25, alpha=TRUE, legend.x="center", bty="n")
-diagram(e150, alpha=TRUE, add=TRUE, col="red")
+a25 <- affinity(pH = c(4, 12))
+a150 <- affinity(pH = c(4, 12), T = 150)
+diagram(a25, legend.x = "topleft", bty = "n")
+diagram(a150, add = TRUE, col = "red")
+e25 <- equilibrate(a25, loga.balance = -3)
+e150 <- equilibrate(a150, loga.balance = -3)
+diagram(e25, ylim = c(-6, 0), bty = "n")
+diagram(e150, add = TRUE, col = "red")
+diagram(e25, alpha = TRUE, legend.x = "center", bty = "n")
+diagram(e150, alpha = TRUE, add = TRUE, col = "red")
```
-The distribution of aqueous carbonate species is a classic example of an equilibrium calculation.
-We can begin by plotting the affinities for equal activities of the species.
+
+The distribution of aqueous carbonate species as a function of pH (Bjerrum plot) is a classic example of an equilibrium calculation.
+We can begin by plotting the affinities calculated at 25 °C and 150 °C with equal activities of the species.
Here, CO<sub>2</sub> is in the basis, so it has zero affinity, which is greater than the affinities of HCO<sub>3</sub><sup>-</sup> and CO<sub>3</sub><sup>-2</sup> at low pH:
-```{r bjerrum_1, eval=FALSE}
-basis("CHNOS+")
-species(c("CO2", "HCO3-", "CO3-2"))
-a25 <- affinity(pH=c(4, 12))
-a150 <- affinity(pH=c(4, 12), T=150)
-diagram(a25, legend.x="topleft", bty="n")
-diagram(a150, add=TRUE, col="red")
+```{r bjerrum_diagram, echo=1:7, eval=FALSE}
```
Now we use <span style="color:green">`equilibrate()`</span> to calculate the activities of species.
Our balancing constraint is that the total activity of C is 10<sup>-3</sup>.
This shows a hypothetical *metastable equilibrium*; we know that for true equilibrium the total activity of C is affected by pH.
-```{r bjerrum_2, eval=FALSE}
-e25 <- equilibrate(a25, loga.balance=-3)
-e150 <- equilibrate(a150, loga.balance=-3)
-diagram(e25, ylim=c(-6, 0), bty="n")
-diagram(e150, add=TRUE, col="red")
+```{r bjerrum_diagram, echo=8:11, eval=FALSE}
```
To display the species distribution, or degree of formation, use the `alpha=TRUE` argument:
-```{r bjerrum_3, eval=FALSE}
-diagram(e25, alpha=TRUE, legend.x="center", bty="n")
-diagram(e150, alpha=TRUE, add=TRUE, col="red")
+```{r bjerrum_diagram, echo=12:13, eval=FALSE}
```
The possible reactions between species are all balanced on 1 C.
@@ -812,22 +809,22 @@
To demonstrate this feature, let's consider the distribution of carbon among organic and inorganic species in the hydrothermal mixing scenario described by @SS98.
First we define the basis and add two inorganic species.
-The `index.return=FALSE` argument tells <span style="color:green">`info()`</span> to return the index (number) of the species in the current species definition; these indices are saved for use below:
+The `index.return = TRUE` argument tells <span style="color:green">`info()`</span> to return the index (number) of the species in the current species definition; these indices are saved for use below:
```{r groups_basis, results="hide", message=FALSE}
basis("CHNOS+")
-ii <- species(c("CO2", "HCO3-"), index.return=TRUE)
+ii <- species(c("CO2", "HCO3-"), index.return = TRUE)
```
Next, we add each group of organic species: C<sub>1</sub>--C<sub>8</sub> alcohols, C<sub>3</sub>--C<sub>8</sub> ketones, C<sub>2</sub>--C<sub>12</sub> carboxylic acids and their corresponding anions, and C<sub>2</sub>--C<sub>8</sub> alkenes.
We provide <span style="color:green">`info()`</span> with a set of `ispecies` values to select these species.
-The species in each group are ordered by carbon number in the database, so the set is made from the starting and ending indices using R's `seq()` function, wrapped by `seq2()` to make the code shorter.
+The species in each group are ordered by carbon number in the database, so the set is made from the starting and ending indices using R's `seq()` function, wrapped by a `seq2()` we write to make the code shorter.
```{r groups_species, message=FALSE}
seq2 <- function(x) seq(x[1], x[2])
-ia <- species(seq2(info(c("methanol", "octanol"))), index.return=TRUE)
-ik <- species(seq2(info(c("acetone", "2-octanone"))), index.return=TRUE)
+ia <- species(seq2(info(c("methanol", "octanol"))), index.return = TRUE)
+ik <- species(seq2(info(c("acetone", "2-octanone"))), index.return = TRUE)
ic <- species(seq2(info(c("acetic acid","n-dodecanoic acid"))),index.return=TRUE)
-ica <- species(seq2(info(c("acetate", "n-dodecanoate"))), index.return=TRUE)
-ie <- species(seq2(info(c("ethylene", "octene"))), index.return=TRUE)
+ica <- species(seq2(info(c("acetate", "n-dodecanoate"))), index.return = TRUE)
+ie <- species(seq2(info(c("ethylene", "octene"))), index.return = TRUE)
```
Now we read two data files that contain values of `r logfO2` and pH as a function of temperature, digitized from Figure 5 of Schulte and Shock (1998).
@@ -838,9 +835,9 @@
Because of the noise introduced by digitization of the figure, we smooth the data using R's `smooth.spline()`; the lower *T* limit reflects the absence of data below this temperature in the figure for log*f*<sub>O<sub>2</sub></sub>.
```{r groups_data}
O2dat <- read.csv(system.file(
- "extdata/cpetc/SS98_Fig5a.csv", package="CHNOSZ"))
+ "extdata/cpetc/SS98_Fig5a.csv", package = "CHNOSZ"))
pHdat <- read.csv(system.file(
- "extdata/cpetc/SS98_Fig5b.csv", package="CHNOSZ"))
+ "extdata/cpetc/SS98_Fig5b.csv", package = "CHNOSZ"))
T <- seq(8, 350)
O2 <- predict(smooth.spline(O2dat$T, O2dat$logfO2), T)$y
pH <- predict(smooth.spline(pHdat$T, pHdat$pH), T)$y
@@ -854,8 +851,8 @@
The ability to vary the activity of the balanced basis species is not yet implemented in CHNOSZ, so a single value is used here.
```
```{r groups_affinity, message=FALSE, cache=TRUE}
-a <- affinity(T=T, O2=O2, pH=pH)
-e <- equilibrate(a, loga.balance=-2.5)
+a <- affinity(T = T, O2 = O2, pH = pH)
+e <- equilibrate(a, loga.balance = -2.5)
```
At last we come to the diagram.
@@ -863,18 +860,18 @@
When summing the activities of species in the groups, each activity is multiplied first by the balance coefficient on that species.
Therefore, the total activity is that of a basis species (or of an element that only in that basis species, like carbon in this example).
```{r groups_diagram, echo=1:4, eval=FALSE}
-par(mfrow=c(1, 3))
-groups <- list(inorganic=ii, alcohols=ia, ketones=ik,
- `carboxylic acids`=c(ic, ica), alkenes=ie)
-diagram(e, alpha=TRUE, groups=groups, legend.x=NA)
+par(mfrow = c(1, 3))
+groups <- list(inorganic = ii, alcohols = ia, ketones = ik,
+ `carboxylic acids` = c(ic, ica), alkenes = ie)
+diagram(e, alpha = TRUE, groups = groups, legend.x = NA)
# plot only alcohols
names <- within(species(), name[-ia] <- "")$name
-lty <- ifelse(names=="", 0, 1)
-diagram(e, alpha=TRUE, ylim=c(0, 0.5), lty=lty, names=names)
+lty <- ifelse(names == "", 0, 1)
+diagram(e, alpha = TRUE, ylim = c(0, 0.5), lty = lty, names = names)
# plot only ketones
names <- within(species(), name[-ik] <- "")$name
-lty <- ifelse(names=="", 0, 1)
-diagram(e, alpha=TRUE, ylim=c(0, 0.13), lty=lty, names=names)
+lty <- ifelse(names == "", 0, 1)
+diagram(e, alpha = TRUE, ylim = c(0, 0.13), lty = lty, names = names)
```
That makes a diagram that is similar to Figure 6b of Shock and Schulte (1998).
@@ -893,28 +890,28 @@
How about the choice between balancing constraints?
Be default, <span style="color:green">`equilibrate()`</span> and <span style="color:green">`diagram()`</span> balance reactions on the first basis species that is present in each of the species of interest.
Let's look at some amino acids in a hypothetical metastable equilibrium.
-This calculation is based on one described by @Sho90b for five amino acids, but here we include 20 proteinogenic amino acids, whose names are returned by `aminoacids("")`.
+This calculation is based on one described by @Sho90b for five amino acids, but here we include 20 proteinogenic amino acids, whose names are returned by <span style="color:green">`aminoacids("")`</span>.
We use <span style="color:green">`ZC.col()`</span> to generate colors based on the average oxidation state of carbon of the amino acids (red and blue for relatively reduced and oxidized).
```{r aminoacids_setup, results="hide", message=FALSE}
basis("CHNOS")
basis("CO2", "gas")
swap.basis("NH3", "N2")
species(aminoacids(""))
-a <- affinity(O2=c(-50, -25, 300), CO2=c(-10, 15, 300), T=250, P=265)
+a <- affinity(O2 = c(-50, -25, 300), CO2 = c(-10, 15, 300), T = 250, P = 265)
aa.ZC <- ZC(info(aminoacids("")))
col <- ZC.col(aa.ZC)
```
To make plots using different balance constraints, let's write a simple function that sets the `balance` argument of <span style="color:green">`diagram()`</span> and adds a title to the plot.
The first plot is the most similar to Figure 4 of Shock (1990), except for the absence of alanine (probably due to different thermodynamic data) and the presence of some other amino acids.
-There, we set `balance=1`, which indicates that moles of species are conserved; this is equivalent to balancing on the amino acid backbone.
+There, we set `balance = 1`, which indicates that moles of species are conserved; this is equivalent to balancing on the amino acid backbone.
The remaining plots balance on each of the basis species (except for O<sub>2</sub>), then on volume (in the last plot).
```{r aafun, fig.fullwidth=TRUE, fig.width=12.5, fig.height=2.5, small.mar=TRUE, dpi=dpi, out.width="100%", message=FALSE, results="hide", cache=TRUE, pngquant=pngquant}
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/chnosz -r 143
More information about the CHNOSZ-commits
mailing list