[CHNOSZ-commits] r127 - in pkg/CHNOSZ: . R inst man vignettes
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sat Feb 4 13:55:30 CET 2017
Author: jedick
Date: 2017-02-04 13:55:30 +0100 (Sat, 04 Feb 2017)
New Revision: 127
Added:
pkg/CHNOSZ/vignettes/CHNOSZ.dia
pkg/CHNOSZ/vignettes/anintro.Rmd
Removed:
pkg/CHNOSZ/vignettes/anintro.Rnw
pkg/CHNOSZ/vignettes/anintro.lyx
pkg/CHNOSZ/vignettes/flowchart.dia
pkg/CHNOSZ/vignettes/flowchart.pdf
pkg/CHNOSZ/vignettes/newintro.Rmd
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/R/basis.R
pkg/CHNOSZ/R/diagram.R
pkg/CHNOSZ/inst/NEWS
pkg/CHNOSZ/man/diagram.Rd
pkg/CHNOSZ/vignettes/vig.bib
Log:
replace anintro.Rnw with updated anintro.Rmd
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2017-02-04 02:51:27 UTC (rev 126)
+++ pkg/CHNOSZ/DESCRIPTION 2017-02-04 12:55:30 UTC (rev 127)
@@ -1,6 +1,6 @@
Date: 2017-02-04
Package: CHNOSZ
-Version: 1.0.8-16
+Version: 1.0.8-17
Title: Chemical Thermodynamics and Activity Diagrams
Author: Jeffrey Dick
Maintainer: Jeffrey Dick <j3ffdick at gmail.com>
Modified: pkg/CHNOSZ/R/basis.R
===================================================================
--- pkg/CHNOSZ/R/basis.R 2017-02-04 02:51:27 UTC (rev 126)
+++ pkg/CHNOSZ/R/basis.R 2017-02-04 12:55:30 UTC (rev 127)
@@ -174,15 +174,18 @@
# pH transformation
if("pH" %in% species) {
logact[species=="pH"] <- -logact[species=="pH"]
- species[species=="pH"] <- "H+"
+ if(!is.null(logact)) species[species=="pH"] <- "H+"
}
# Eh and pe transformations
- if(any(c("pe","Eh") %in% species)) {
+ if("pe" %in% species) {
logact[species=="pe"] <- -logact[species=="pe"]
+ if(!is.null(logact)) species[species=="pe"] <- "e-"
+ }
+ if("Eh" %in% species) {
# 20090209 should be careful with this conversion as it's only for 25 deg C
- # to be sure, just don"t call species("Eh")
- logact[species=="Eh"] <- -convert(logact[species=="Eh"],"pe")
- species[species %in% c("pe","Eh")] <- "e-"
+ # to be sure, just don't call species("Eh")
+ if(!is.null(logact)) logact[species=="Eh"] <- -convert(logact[species=="Eh"],"pe")
+ species[species=="Eh"] <- "e-"
}
## if all species are in the existing basis definition,
## *and* at least one of state or logact is not NULL
Modified: pkg/CHNOSZ/R/diagram.R
===================================================================
--- pkg/CHNOSZ/R/diagram.R 2017-02-04 02:51:27 UTC (rev 126)
+++ pkg/CHNOSZ/R/diagram.R 2017-02-04 12:55:30 UTC (rev 127)
@@ -421,7 +421,7 @@
}
if(is.null(fill)) fill <- "transparent"
else if(isTRUE(fill[1]=="rainbow")) fill <- rainbow(ngroups)
- else if(isTRUE(fill[1]=="heat")) fill <- heat.colors(ngroups)
+ else if(isTRUE(fill[1] %in% c("heat", "terrain", "topo", "cm"))) fill <- get(paste0(fill[1], ".colors"))(ngroups)
fill <- rep(fill, length.out=ngroups)
# the x and y values
xs <- eout$vals[[1]]
Modified: pkg/CHNOSZ/inst/NEWS
===================================================================
--- pkg/CHNOSZ/inst/NEWS 2017-02-04 02:51:27 UTC (rev 126)
+++ pkg/CHNOSZ/inst/NEWS 2017-02-04 12:55:30 UTC (rev 127)
@@ -1,4 +1,4 @@
-CHANGES IN CHNOSZ 1.0.8-16 (2017-02-04)
+CHANGES IN CHNOSZ 1.0.8-17 (2017-02-04)
---------------------------------------
- Add "AA" as a keyword for preset species in basis() (cysteine,
@@ -39,7 +39,7 @@
message() from base R. As a result, the messages don't appear in Sweave
vignettes, but can now be turned on or off easily in knitr vignettes.
-- Add vignette newintro.Rmd.
+- Replace anintro.Rnw (Sweave) with updated anintro.Rmd (knitr, tufte).
CHANGES IN CHNOSZ 1.0.8 (2016-05-28)
------------------------------------
Modified: pkg/CHNOSZ/man/diagram.Rd
===================================================================
--- pkg/CHNOSZ/man/diagram.Rd 2017-02-04 02:51:27 UTC (rev 126)
+++ pkg/CHNOSZ/man/diagram.Rd 2017-02-04 12:55:30 UTC (rev 127)
@@ -85,6 +85,7 @@
On 2-D diagrams, the fields represent the species with the highest equilibrium activity.
\code{fill} determines the color of the predominance fields, \code{col} that of the boundary lines.
By default, \code{\link{heat.colors}} are used to fill the predominance fields in diagrams on the screen plot device.
+\code{fill} can be any color specification, or the word \samp{rainbow}, \samp{heat}, \samp{terrain}, \samp{topo}, or \samp{cm}, indicating a palette from \pkg{grDevices}.
As of CHNOSZ 1.0.8-11, a new default line-drawing procedure has been implemented.
This uses \code{\link{contour}} to draw smooth-looking diagonal and curved lines, at the expense of not coinciding exactly with the rectangular grid (which is still used for drawing colors).
Copied: pkg/CHNOSZ/vignettes/CHNOSZ.dia (from rev 126, pkg/CHNOSZ/vignettes/flowchart.dia)
===================================================================
(Binary files differ)
Copied: pkg/CHNOSZ/vignettes/anintro.Rmd (from rev 126, pkg/CHNOSZ/vignettes/newintro.Rmd)
===================================================================
--- pkg/CHNOSZ/vignettes/anintro.Rmd (rev 0)
+++ pkg/CHNOSZ/vignettes/anintro.Rmd 2017-02-04 12:55:30 UTC (rev 127)
@@ -0,0 +1,550 @@
+---
+title: "An Introduction to [CHNOSZ](http://www.chnosz.net)"
+author: "Jeffrey M. Dick"
+date: "`r Sys.Date()`"
+output:
+ tufte::tufte_html:
+ tufte_features: ["background"]
+ css: "vig.css"
+ toc: true
+ tufte::tufte_handout:
+ citation_package: natbib
+ latex_engine: xelatex
+ tufte::tufte_book:
+ citation_package: natbib
+ latex_engine: xelatex
+vignette: >
+ %\VignetteIndexEntry{An Introduction to CHNOSZ}
+ %\VignetteEngine{knitr::rmarkdown}
+ \usepackage[utf8]{inputenc}
+bibliography: vig.bib
+link-citations: yes
+csl: elementa.csl
+---
+
+```{r options, include=FALSE}
+options(width = 80)
+options(digits = 6)
+```
+
+```{r setup, include=FALSE}
+# invalidate cache when the tufte version changes
+knitr::opts_chunk$set(tidy = FALSE, cache.extra = packageVersion('tufte'))
+options(htmltools.dir.version = FALSE)
+knitr::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) {
+ if (before) par(mar = c(.1, .1, .1, .1)) # tiny margin all around
+})
+knitr::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
+})
+```
+
+# About
+
+This document introduces the basic functionality of CHNOSZ, a package for the [R software environment](http://r-project.org).
+For more information on R, see "An Introduction to R" and other documents in [The R Manuals](http://cran.r-project.org/manuals.html).
+
+CHNOSZ has been developed since 2006 to support research projects in geochemistry and biochemistry.
+The package provides functions and an extensive thermodynamic database that can be used to calculate the stoichiometric and energetic properties of reactions among minerals and organic or inorganic aqueous species.
+These same functions enable calculations of chemical affinities and metastable equilibrium for systems of proteins.
+
+## Installing and loading CHNOSZ
+
+If you have just installed R for the first time, installing CHNOSZ should be as simple as selecting the "Install packages from CRAN" or similar menu item in the R GUI or using the following command to start the package installation process:
+```{marginfigure}
+Or, install the package from a local package file, which you can download from [CRAN](https://cran.r-project.org/package=CHNOSZ) or from the [CHNOSZ website](http://chnosz.net/download).
+```
+```{r install_CHNOSZ, eval=FALSE}
+install.packages("CHNOSZ")
+```
+
+Then load the CHNOSZ package to make its functions available in your R session.
+```{r library_CHNOSZ}
+library(CHNOSZ)
+```
+
+Then load the thermo object, which contains the thermodynamic database and system settings for CHNOSZ.
+```{r data_thermo}
+data(thermo)
+```
+
+## Getting help
+
+After CHNOSZ is installed, type <span style="color:blue">`help.start()`</span> to browse the R help documents, then choose "Packages" followed by "CHNOSZ".
+That shows an index of the manual (help page) for each function; many of the help pages include examples.
+There are also links to the *demos* (longer examples) and *vignettes* (more in-depth documentation).
+Suggestions for accessing the documentation are indicated here with <span style="color:blue">blue text</span>.
+For example, read <span style="color:blue">`?"CHNOSZ-package"`</span> to get an overview of the package and a list of features.
+```{marginfigure}
+That page identifies some features as experimental, i.e. not based on published algorithms or extensively compared with published results.
+```
+
+## Organization of major functions
+
+CHNOSZ is made up of a set of functions and supporting datasets.
+The major components of the package are shown in the figure below, which is an updated version of the flowchart from @Dic08.
+Rectangles and ellipses represent functions and datasets; bold text indicates primary functions.
+
+
+
+Many functions in CHNOSZ have no side effects.
+That is, the function only returns a result; to use the result in other functions, it can be assigned to a variable with `<-`.
+Major functions without side effects in CHNOSZ are:
+
+* `info()`: search for species in the thermodynamic database
+* `subcrt()`: calculate the thermodynamic properties of species and reactions
+* `affinity()`: calculate the affinities of formation reactions using given chemical activities
+* `equilibrate()`: calculate the equilibrium chemical activities of the species of interest
+* `diagram()`: plot the results
+
+Some functions in CHNOSZ do have side effects: they modify the `thermo` data object in the current R session.
+In the text (not code) of this document, the names of these functions are shown in <span style="color:red">red text</span>.
+The major functions with side effects are:
+
+* <span style="color:red">`basis()`</span>: set the basis species and their chemical activities
+* <span style="color:red">`species()`</span>: set the species of interest and their (non-equilibrium) chemical activities
+* <span style="color:red">`data(thermo)`</span>: reset the database, restoring all settings to their default values.
+
+The following pseudocode shows a common sequence of commands.
+In actual usage, the `...` are replaced by arguments that define the chemical makeup and range of conditions of the system:
+```{r, eval=FALSE}
+basis(...)
+species(...)
+a <- affinity(...)
+e <- equilibrate(a) ## optional
+diagram(e) ## or diagram(a)
+data(thermo) ## clear system settings
+```
+
+Some experimental functions are available:
+
+* (**experimental**) using `revisit()` to calculate/plot summary statistics of the chemical activities of the species of interest and `findit()` to search for combinations of activities of basis species, temperature and/or pressure that optimize those statistics.
+
+# Thermodynamic database and chemical formulas
+
+While an attempt has been made to provide a primary database (`OBIGT.csv`) that is generally internally consistent, all thermodynamic data, calculations, and examples are provided *as is*.
+```{marginfigure}
+For crucial problems, check not only the accuracy of the database, but also the *suitability of the data* for your problem.
+If there is any doubt about the suitability of data, please consult the primary sources (see <span style="color:blue">?`browse.refs`</span>).
+```
+Where possible, data with known or suspected inconsistencies have been placed into a secondary database (`OBIGT-2.csv`) that should be regarded as experimental.
+
+## The `info()` function
+
+The `info()` function provides an interface to the thermodynamic database packaged with CHNOSZ.
+Suppose you are interested in the thermodynamic properties aqueous ethylene.
+You can search for the species by name:
+```{r info_ethylene}
+info("ethylene")
+```
+
+Multiple entries exist for ethylene; the index of the `aq` (aqueous) species is returned by default.
+A second argument can be used to specify a different physical state:
+```{r info_ethylene_gas}
+info("ethylene", "gas")
+```
+
+Knowing that aqueous ethylene is species number 88 in the database, you can again use `info()` to retrieve the set of standard molal thermodynamic properties and equations of state parameters:
+```{r info_88}
+info(88)
+```
+
+This number can be used as an argument (`ispecies`) for other functions in CHNOSZ to uniquely identify any species; some commonly used functions also accept the species names.
+Liquid water is species number 1; it has NA entries in the database because specialized functions are used to compute its properties:
+```{r info_info_water}
+info(info("water"))
+```
+
+## Fuzzy searches
+
+Calling `info()` with a string that does not exactly match the name of any species invokes a fuzzy search of the database:
+```{r width180, include=FALSE}
+options(width = 180)
+```
+```{r info_acid}
+info("acid")
+```
+```{r width80, include=FALSE}
+options(width = 80)
+```
+
+The message includes e.g. "uracil" and "metacinnabar" because their names have some similarity to the search term.
+
+As "ribose" is the name of a species in the database, to find species with similar names, add an extra character to the search:
+```{r info_ribose}
+info(" ribose")
+```
+
+The messages may be useful for browsing the database, but owing to their ambiguous results, these fuzzy searches return an `NA` value for the species index.
+
+## Counting elements, chemical formulas, `ZC()`
+
+Continuing with the example of ethylene, let's look at its chemical formula:
+```{r info_88_formula}
+info(88)$formula
+```
+
+We can use the `makeup()` function to count the elements in the formula, followed by `as.chemical.formula()` to return to the formula:
+```{r makeup_88}
+makeup(88)
+as.chemical.formula(makeup(88))
+```
+
+For organic species, a simple calculation of the average oxidation state of carbon ($Z_C$) is possible given the species index, chemical formula, or elemental count:
+```{r ZC_88}
+ZC(88)
+ZC(info(88)$formula)
+ZC(makeup(88))
+```
+
+# Calculating thermodynamic properties
+
+To calculate the standard molal properties of species and reactions, use `subcrt()`.
+```{marginfigure}
+The inspiration for the name `subcrt()`, and the source of the FORTRAN subroutine used to calculate the thermodynamic properties of H<sub>2</sub>O, is SUPCRT (Johnson et al., 1992).
+```
+<sup>[- at JOH92]</sup>
+If no reaction coefficients are given, `subcrt()` calculates the standard molal properties of invididual species:
+```{r subcrt_water}
+subcrt("water")
+```
+
+That uses the default temperature and pressure settings: equally-spaced temperature grid from 0 to 350 °C at *P*<sub>sat</sub>.
+```{marginfigure}
+*P*<sub>sat</sub> is 1 bar below 100 °C, or the pressure of liquid-vapor coexistence at higher temperatures.
+```
+The columns in the output are temperature, pressure, density of water, logarithm of the equilibrium constant (only meaningful for reactions; see below), standard molal Gibbs energy and enthalpy of formation from the elements, standard molal entropy, volume, and heat capacity.
+```{marginfigure}
+The corresponding units are °C (`T`), bar (`P`), g cm<sup>-3</sup> (`rho`), cal mol<sup>-1</sup> (`G` and `H`), cal K<sup>-1</sup> mol<sup>-1</sup> (`S`), cm<sup>3</sup> mol<sup>-1</sup> (`V`), and cal K<sup>-1</sup> mol<sup>-1</sup> (`Cp`).
+```
+
+A custom temperature-pressure grid can be specified.
+Here, we calculate the properties of H<sub>2</sub>O on a *T*, *P* grid in the supercritical region, with conditions grouped by pressure:
+```{marginfigure}
+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
+```
+
+```{r subcrt_water_plot, fig.margin=TRUE, fig.width=4, fig.height=4, small.mar=TRUE, dpi=50, 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}
+H2O <- subcrt("water", T=seq(0, 1000, 100), P=c(NA, seq(1, 500, 1)), grid="T")
+H2O <- H2O$out$water
+plot(H2O$P, H2O$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. `write.table()` or `plot()` for further processing:
+```{r subcrt_water_plot, eval=FALSE}
+```
+
+## Changing units
+
+The default units of temperature, pressure, and energy are °C, bar, and calories.
+The functions <span style="color:red">`T.units()`</span>, <span style="color:red">`P.units()`</span>, and <span style="color:red">`E.units()`</span> can be used to change the units used by various functions in CHNOSZ.
+What is the Gibbs energy (J/mol) of aqueous methane at 298.15 K and 0.1 MPa?
+```{r methane_units, message=FALSE}
+T.units("K")
+P.units("MPa")
+E.units("J")
+subcrt("methane", T=298.15, P=0.1)$out$methane$G
+data(thermo) ## restore default settings
+```
+
+A related function, `convert()`, can be used to convert given values between units.
+Let's convert the standard Gibbs energy of aqueous methane listed in the database from cal/mol to J/mol:
+```{r methane_G, message=FALSE}
+convert(info(info("methane"))$G, "J")
+```
+
+As expected, we get the same result from both operations.
+
+# Properties of reactions
+
+## Reaction definitions
+
+To calculate the thermodynamic properties of reactions, give the names of species, the physical states (optional), and reaction coefficients as the arguments to `subcrt()`.
+Here we calculate properties for the dissolution of CO<sub>2</sub>:
+```{marginfigure}
+Because of aqueous speciation, this doesn't give the _solubility_ of CO<sub>2</sub>.
+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))
+```
+
+<p> <!--- needed to make the citation appear correction 20170202 -->
+In order to make a plot like Figure 18 of @MSS13, let's run more calculations and store the results.
+In addition to the reaction definition, we specify a greater number of temperature points than the default:
+```{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
+logK <- data.frame(T, CO2, CO, CH4)
+```
+```{r dissolution_plot, fig.margin=TRUE, fig.width=4, fig.height=4, small.mar=TRUE, dpi=50, 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}
+matplot(logK[, 1], logK[, -1], type="l",
+ xlab=axis.label("T"), ylab=axis.label("logK"))
+```
+```{r dissolution, eval=FALSE}
+```
+
+Now we can make the plot.
+Here, the `axis.label()` function of CHNOSZ is used to create formatted axis labels:
+```{r dissolution_plot, eval=FALSE}
+```
+
+## Unbalanced reactions
+
+A balanced chemical reaction conserves mass.
+`subcrt()` won't stop you from running an unbalanced reaction, but it will give you a warning:
+```{r subcrt_unbalanced, results="hide"}
+subcrt(c("CO2", "CH4"), c(-1, 1))
+```
+
+In other words, to balance the reaction, we should add 4 H to the left and 2 O to the right.
+That could be done manually be redefining the reaction with the appropriate species.
+There is another option: balancing the reaction automatically using basis species.
+
+## Setting the basis species
+
+_Basis species_ are a minimal number of chemical species that represent the compositional variation in a system.
+The basis species are similar to thermodynamic components, but can include charged species.
+You might want to use basis species to automatically balance reactions.
+A basis setting is also required for making chemical activity diagrams.
+
+Let's start with an example that doesn't work:
+```{r basis_singular, error=TRUE}
+basis(c("CO2", "H2", "H2CO2"))
+```
+
+That set of species has a singular (non-invertible) stoichiometric matrix.
+An error would also result from either an underdetermined or overdetermined system.
+A valid set of basis species has an invertible stoichiometric matrix and the same number of species as elements:
+```{r basis_CHO}
+basis(c("CO2", "H2", "H2O"))
+```
+The composition of any species made up of C, H, and O can be represented by a single linear combination of these basis species.
+
+## Auto-balancing reactions
+
+Methanogenic metabolism in reducing environments may take advantage of acetoclastic or hydrogenotrophic processes.
+To consider reactions involving a charged species (acetate), let's define a basis with H<sup>+</sup>:
+```{r basis_CHOZ}
+basis(c("CO2", "H2", "H2O", "H+"))
+```
+
+By identifying species *other than* the basis species, the reactions will be automatically balanced.
+This produces the balanced reaction for acetoclastic methanogenesis:
+```{r subcrt_acetoclastic, message=FALSE}
+subcrt(c("acetate", "methane"), c(-1, 1))$reaction
+```
+
+We can similarly consider reactions for hydrogenotrophic methanogenesis as well as acetate oxidation (no production of methane):
+```{r subcrt_methanogenesis, message=FALSE}
+acetate_oxidation <- subcrt("acetate", -1)
+hydrogenotrophic <- subcrt("methane", 1)
+acetoclastic <- subcrt(c("acetate", "methane"), c(-1, 1))
+```
+
+<p>
+Use `describe.reaction()` 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=50, out.width="100%"}
+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)
+```
+
+## Non-standard state affinities and Gibbs energies
+
+Usually, `subcrt()` returns only standard state thermodynamic properties.
+```{marginfigure}
+The standard state adopted for H<sub>2</sub>O is unit activity of the pure component at any *T* and *P*.
+The standard state for aqueous species is unit activity of a hypothetical one molal solution referenced to infinite dilution at any *T* and *P*.
+```
+Thermodynamic models often consider a non-standard state (i.e. non-unit activity).
+The activities of basis species can be modified with <span style="color:red">`basis()`</span>, and those of the other species using the `logact` argument in `subcrt()`.
+
+Let us calculate the chemical affinity of acetoclastic methanogenesis.
+```{marginfigure}
+The affinity is equal to the negative of the overall (non-standard) Gibbs energy change of the reaction.
+```
+We begin by changing the energy units to Joules.
+Then, we change the state of H<sub>2</sub> and CO<sub>2</sub> in the basis from `aq` (aqueous) to `gas`, and set the logarithm of fugacity of gaseous H<sub>2</sub> and the pH, using values from @MDS_13.
+The activity of acetate and fugacity of methane, as well as temperature and pressure, are set in the call to `subcrt()`:
+```{r basis_mayumi, message=FALSE, results="hide"}
+E.units("J")
+basis(c("CO2", "H2"), "gas")
+basis(c("H2", "pH"), c(-3.92, 7.3))
+```
+```{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
+```
+
+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
+```
+
+Under the specified conditions, the affinities of hydrogenotrophic and acetoclastic methanogenesis are somewhat greater than and less than 20 kJ, respectively.
+This result matches Figure 4b in Mayumi et al. (2013) at unit fugacity of CO<sub>2</sub>.
+
+We can go even further and reproduce their plot.
+```{marginfigure}
+The reproduction is not identical, owing to differences of thermodynamic data and of calculations of the effects of temperature and pressure.
+```
+To make the code neater, we write a function that can run any of the reactions:
+```{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
+}
+```
+
+<p>
+Now we're ready to calculate and plot the affinities.
+Here, we use the `lapply()` function of R to list the results at two values of logarithm of fugacity of CO<sub>2</sub>.
+We insert an empty reaction to get a line at zero affinity.
+`do.call(rbind, Adat)` turns the list into a data frame that can be plotted with `matplot()`.
+There, we plot the negative affinities, equal to Gibbs energy, as shown in the plot of Mayumi et al. (2013).
+```{r methanogenesis_plot, fig.margin=TRUE, fig.width=4.1, fig.height=4.1, small.mar=TRUE, dpi=50, out.width="100%", echo=FALSE, message=FALSE, fig.cap="Gibbs energies of acetate oxidation and methanogenesis (after Mayumi et al., 2013).", cache=TRUE}
+Adat <- lapply(c(-3, 3), function(logfCO2) {
+ basis("CO2", logfCO2)
+ data.frame(logfCO2,
+ rxnfun(c(0, 0))$A,
+ rxnfun(c(-1, 0))$A,
+ rxnfun(c(-1, 1))$A,
+ rxnfun(c(0, 1))$A
+ )
+})
+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"))
+legend("topleft", c("acetate oxidation", "acetoclastic methanogenesis",
+ "hydrogenotrophic methanogenesis"), lty=1, col=2:4)
+```
+```{r methanogenesis_plot, echo=TRUE, eval=FALSE}
+```
+
+Let't not forget to clear the system settings, which were modified by <span style="color:red">`basis()`</span> and <span style="color:red">`E.units()`</span>, before running other calculations:
+```{r data_thermo, message=FALSE}
+```
+
+# Affinity, species of interest, and potential diagrams
+
+`affinity()` offers calculations of chemical affinity of formation reactions over a configurable range of T, P, and activities of basis species.
+
+By *formation reaction* is meant the stoichiometric requirements for formation of one mole of any species from the basis species.
+The <span style="color:red">`species()`</span> function is used to set these *species of interest*.
+Let's consider the stoichiometry of some aqueous sulfur-bearing species.
+```{r basis_CHNOSZ, results="hide"}
+basis("CHNOS+")
+```
+```{r sulfur_species}
+species(c("SO4-2", "HS-", "H2S", "HSO4-"))
+```
+
+There we used <span style="color:red">`basis()`</span> with a keyword to identify a preset basis definition.
+Now, we can use `affinity()` to calculate the affinities of the formation reactions of each of the species:
+```{marginfigure}
+The values returned by `affinity()` are dimensionless, i.e. log<sub>10</sub>(*A*/*RT*).
+```
+```{r affinity}
+affinity()$values
+```
+
+The same result (in energetic units) could be obtained using `subcrt()`, but `affinity()` 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 get 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-")
+```
+
+The <span style="color:red">swap.basis()</span> changed the basis species and recalculated their activities, but preserved the species of interest.
+```{marginfigure}
+That is, running `affinity()$values` again would give the same result.
+```
+
+<p>
+```{r EhpH_plot, fig.margin=TRUE, fig.width=4, fig.height=4, small.mar=TRUE, dpi=50, out.width="100%", echo=FALSE, message=FALSE, cache=TRUE}
+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 `diagram()` function uses the maximum affinity method to make a potential diagram:
+```{r EhpH_plot_echo, eval=FALSE}
+diagram(a)
+```
+
+Note that the calculation of affinity implies a non-equilibrium reference state of equal activities of species (see above).
+Generally, then, `diagram()` gives a *potential diagram* because it shows regions of maximum affinity.
+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=50, out.width="100%", echo=FALSE, message=FALSE, cache=TRUE}
+diagram(a, fill="terrain", lwd=3, lty=3,
+ names=c("sulfate", "bisulfide", "hydrogen sulfide", "bisulfate"),
+ tplot=FALSE, main="sulfur species, 25 °C", bty="n")
+```
+The default colors for diagrams shown on screen are R's "heat colors".
+Some arguments in `diagram()` can be used to control the color, labels, and lines, and title (`main`).
+The `tplot` argument turns off plot customizations used in CHNOSZ.
+Additional arguments can be passed via `...`; here, we use that to remove the box around the plot:
+```{r EhpH_plot_color, echo=TRUE, eval=FALSE}
+```
+
+## Mosaic diagrams
+
+If sulfur is in the basis species, then we should consider that its speciation is sensitive to Eh and pH, as shown in the preceding diagram.
+Mosaic diagrams, which are often shown for metal oxide, sulfide, and carbonate minerals, account for speciation of the basis species.
+These diagrams are made by constructing individual diagrams for the possible basis species.
+The individual diagrams are then combined, each one contributing to the final diagram only in the range of stability of the corresponding basis species.
+
+Here, we use the `mosaic()` function in CHNOSZ to make a diagram for aqueous species and minerals in the Cu-S-Cl-H<sub>2</sub>O system, similar Figure 5a of @CPCC17.
+The key argument is `bases`, which identifies the candidate basis species (starting with the one in the current basis set).
+The other arguments, like those of `affinity()`, specify the ranges of the variables.
+
+# Equilibration
+
+## Getting from affinity to equilibrium
+
+## Setting the constraints
+
+# Other things you can do with affinity
+
+## Buffers
+
+## T-, P-, activity-transect
+
+# Other things you can do with equilibrate
+
+## Choosing different balancing constraints
+
+# Other things you can do with diagrams
+
+## Groups of species
+
+# Proteins
+
+## Group additivity
+
+## Sources of amino acid data
+
+## Ionization
+
+## Normalizing for different lengths
+
+# Functions outside of the main workflow
+
+transfer, wjd, eqdata, RH2obigt, EOSregress
Deleted: pkg/CHNOSZ/vignettes/anintro.Rnw
===================================================================
--- pkg/CHNOSZ/vignettes/anintro.Rnw 2017-02-04 02:51:27 UTC (rev 126)
+++ pkg/CHNOSZ/vignettes/anintro.Rnw 2017-02-04 12:55:30 UTC (rev 127)
@@ -1,1077 +0,0 @@
-%% LyX 2.2.2 created this file. For more info, see http://www.lyx.org/.
-%% Do not edit unless you really know what you are doing.
-\documentclass[english,noae,round]{article}
-\usepackage{mathpazo}
-\usepackage[T1]{fontenc}
-\usepackage[utf8]{inputenc}
-\usepackage[letterpaper]{geometry}
-\geometry{verbose,tmargin=2.5cm,bmargin=2.5cm,lmargin=2.5cm,rmargin=2.5cm}
-\usepackage{color}
-\usepackage{babel}
-\usepackage{amsbsy}
-\usepackage{amssymb}
-\usepackage{graphicx}
-\usepackage[authoryear]{natbib}
-\usepackage[unicode=true,
- bookmarks=true,bookmarksnumbered=false,bookmarksopen=false,
- breaklinks=false,pdfborder={0 0 1},backref=false,colorlinks=true]
- {hyperref}
-\hypersetup{pdftitle={An introduction to CHNOSZ},
- pdfauthor={Jeffrey M. Dick},
- citecolor=blue}
-\usepackage{breakurl}
-
-\makeatletter
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Textclass specific LaTeX commands.
-<<echo=F>>=
- if(exists(".orig.enc")) options(encoding = .orig.enc)
-@
-\newenvironment{lyxcode}
-{\par\begin{list}{}{
-\setlength{\rightmargin}{\leftmargin}
-\setlength{\listparindent}{0pt}% needed for AMS classes
-\raggedright
-\setlength{\itemsep}{0pt}
-\setlength{\parsep}{0pt}
-\normalfont\ttfamily}%
- \item[]}
-{\end{list}}
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% User specified LaTeX commands.
-%\VignetteIndexEntry{An introduction to CHNOSZ}
-
-% so DOIs in bibliography show up as hyperlinks
-\newcommand*{\doi}[1]{\href{http://dx.doi.org/#1}{doi: #1}}
-
-\makeatother
-
-\begin{document}
-
-\title{An introduction to CHNOSZ}
-
-\author{Jeffrey M. Dick}
-
-\maketitle
-<<echo=FALSE>>=
-options(width=90)
-@
-
-
-\section{About}
-
-This document will orient you to the basic functionality of CHNOSZ,
-a package for the R software environment. R is a powerful language
-and also very fun to use. Don't worry if you're new to it; just plow
-through the examples below and you'll start to get the hang of it.
-If you want a more structured approach to learning the language, there
-are some excellent guides in the Manuals section of the \href{http://www.r-project.org}{R Project page}.
-
-The package was developed since 2006 to support a research project
-on the thermodynamic properties of proteins. Since that time, the
-functions in the package have expanded to include calculation of the
-thermodynamic properties of reactions, and especially the construction
-of equilibrium chemical activity diagrams for both inorganic and organic
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/chnosz -r 127
More information about the CHNOSZ-commits
mailing list