[CHNOSZ-commits] r876 - in pkg/CHNOSZ: . inst vignettes
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Apr 7 06:20:54 CEST 2025
Author: jedick
Date: 2025-04-07 06:20:53 +0200 (Mon, 07 Apr 2025)
New Revision: 876
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/inst/NEWS.Rd
pkg/CHNOSZ/vignettes/anintro.Rmd
pkg/CHNOSZ/vignettes/equilibrium.Rmd
Log:
Revise and shorten anintro.Rmd with AI assistance
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2025-03-08 12:15:09 UTC (rev 875)
+++ pkg/CHNOSZ/DESCRIPTION 2025-04-07 04:20:53 UTC (rev 876)
@@ -1,6 +1,6 @@
-Date: 2025-03-08
+Date: 2025-04-07
Package: CHNOSZ
-Version: 2.1.0-47
+Version: 2.1.0-48
Title: Thermodynamic Calculations and Diagrams for Geochemistry
Authors at R: c(
person("Jeffrey", "Dick", , "j3ffdick at gmail.com", role = c("aut", "cre"),
Modified: pkg/CHNOSZ/inst/NEWS.Rd
===================================================================
--- pkg/CHNOSZ/inst/NEWS.Rd 2025-03-08 12:15:09 UTC (rev 875)
+++ pkg/CHNOSZ/inst/NEWS.Rd 2025-04-07 04:20:53 UTC (rev 876)
@@ -15,7 +15,7 @@
\newcommand{\Cp}{\ifelse{latex}{\eqn{C_P}}{\ifelse{html}{\out{<I>C<sub>P</sub></I>}}{Cp}}}
\newcommand{\DG0}{\ifelse{latex}{\eqn{{\Delta}G^{\circ}}}{\ifelse{html}{\out{Δ<I>G</I>°}}{ΔG°}}}
-\section{Changes in CHNOSZ version 2.1.0-45 (2025-02-15)}{
+\section{Changes in CHNOSZ version 2.1.0-48 (2025-04-07)}{
\subsection{OBIGT DEFAULT DATA}{
\itemize{
@@ -87,6 +87,8 @@
\subsection{DOCUMENTATION}{
\itemize{
+ \item Revise and shorten \file{anintro.Rmd} with AI assistance.
+
\item Add FAQ question: Why are mineral stability boundaries curved on
mosaic diagrams?
Modified: pkg/CHNOSZ/vignettes/anintro.Rmd
===================================================================
--- pkg/CHNOSZ/vignettes/anintro.Rmd 2025-03-08 12:15:09 UTC (rev 875)
+++ pkg/CHNOSZ/vignettes/anintro.Rmd 2025-04-07 04:20:53 UTC (rev 876)
@@ -7,54 +7,26 @@
tufte_features: ["background"]
toc: true
mathjax: null
- 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}
%\VignetteEncoding{UTF-8}
bibliography: vig.bib
-link-citations: yes
csl: elementa.csl
---
-<style>
+```{css, echo=FALSE}
html {
font-size: 14px;
}
-body {
- font-family: ‘Times New Roman’, Times, serif;
-}
-li {
- padding: 0.25rem 0;
-}
/* zero margin around pre blocks (looks more like R console output) */
pre {
- margin-top: 0;
- margin-bottom: 0;
+ margin: 0;
+ padding: 0;
}
</style>
-
-```{r options, include=FALSE}
-options(width = 80)
-options(digits = 6)
```
-```{r HTML, include=FALSE}
-## Some frequently used HTML expressions
-logfO2 <- "log<i>f</i><sub>O<sub>2</sub></sub>"
-# Use lowercase here because these tend to be variable names in the examples
-zc <- "<i>Z</i><sub>C</sub>"
-o2 <- "O<sub>2</sub>"
-h2o <- "H<sub>2</sub>O"
-sio2 <- "SiO<sub>2</sub>"
-ch4 <- "CH<sub>4</sub>"
-```
-
```{r setup, include=FALSE}
library(knitr)
@@ -63,52 +35,12 @@
opts_chunk$set(tidy = FALSE, cache.extra = packageVersion('tufte'))
options(htmltools.dir.version = FALSE)
-## Adjust plot margins
-## First one from https://yihui.name/knitr/hooks/
-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
-})
-knit_hooks$set(tiny.mar = function(before, options, envir) {
- if (before) par(mar = c(.1, .1, .1, .1)) # tiny margin all around
-})
-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
-})
-
-## Use pngquant to optimize PNG images
+## Setup pngquant to optimize PNG images
knit_hooks$set(pngquant = hook_pngquant)
pngquant <- "--speed=1 --quality=0-25"
# pngquant isn't available on R-Forge ...
if (!nzchar(Sys.which("pngquant"))) pngquant <- NULL
-# Set dpi 20231129
-knitr::opts_chunk$set(
- dpi = if(nzchar(Sys.getenv("CHNOSZ_BUILD_LARGE_VIGNETTES"))) 72 else 50
-)
-hidpi = if(nzchar(Sys.getenv("CHNOSZ_BUILD_LARGE_VIGNETTES"))) 100 else 50
-
-## 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
-})
-
-## http://stackoverflow.com/questions/30530008/hook-to-time-knitr-chunks
-now = Sys.time()
-knit_hooks$set(timeit = function(before) {
- if (before) { now <<- Sys.time() }
- else {
- paste("%", sprintf("Chunk rendering time: %s seconds.\n", round(Sys.time() - now, digits = 3)))
- }
-})
-timeit <- NULL
-
## Colorize messages 20171031
## Adapted from https://gist.github.com/yihui/2629886#file-knitr-color-msg-rnw
color_block = function(color) {
@@ -117,1499 +49,225 @@
knit_hooks$set(warning = color_block('magenta'), error = color_block('red'), message = color_block('blue'))
```
-# Overview
+# CHNOSZ: An Introduction to Thermodynamic Calculations in R
-This vignette was made for version `r sessionInfo()$otherPkgs$CHNOSZ$Version` of CHNOSZ, a package for the [R software environment](https://www.r-project.org/).
-For more information on R, see "[An Introduction to R](https://cran.r-project.org/manuals.html)" and the [contributed documentation](https://cran.r-project.org/other-docs.html) for R.
+This vignette introduces CHNOSZ, an R package for thermodynamic calculations relevant to geochemistry and geobiochemistry. CHNOSZ provides functions and a thermodynamic database for calculating properties of reactions involving minerals, aqueous species, and gases across a range of temperatures and pressures.
-CHNOSZ has been developed since 2006 as a tool for thermodynamic calculations in geochemistry and geochemical biology.
-The package provides functions and a thermodynamic database that can be used to calculate the stoichiometric and energetic properties of reactions involving minerals and inorganic and/or organic aqueous species.
-These functions also enable calculations of chemical affinities and metastable equilibrium distributions of proteins.
-A major feature of the package is the production of diagrams to visualize the effects of changing temperature, pressure, and activities of basis species on the potential for reactions among various species.
+## Getting Started
-## Installing and loading CHNOSZ
+After installing CHNOSZ from CRAN, load the package:
-After starting R, install CHNOSZ by selecting the "Install packages from CRAN" or similar menu item in the R GUI or by using the following command:
-```{r install_CHNOSZ, eval=FALSE}
-install.packages("CHNOSZ")
-```
-
-Then load the CHNOSZ package to make its data and functions available in your R session:
```{r library_CHNOSZ}
library(CHNOSZ)
```
-CHNOSZ is now ready to go with the default thermodynamic database and an empty system definition.
-After running some calculations, you may want to "start over" with the default values.
-To clear the system settings and restore the default thermodynamic database, use <span style="color:red">`reset()`</span>.
-```{r reset}
-reset()
-```
+This makes the thermodynamic database and functions available in your R session. To restore default settings at any point, use `reset()`.
-Note: Throughout this document, syntax highlighting is applied to the *input* of the code chunks.
-Double hash marks (`##`) precede the *output*, where black text denotes *results* and blue text is used for *messages*.
+## Core Functionality
-## Getting help
+CHNOSZ offers several primary functions for thermodynamic analysis:
-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 pages) 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; this document is a vignette).
+### Functions Without Side Effects (Return Values)
-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}
-"`?`" is a shortcut to R's `help()` function.
-The command here is equivalent to <span style="color:blue">`help("CHNOSZ-package")`</span>.
-```
+* `info()`: Search the thermodynamic database
+* `subcrt()`: Calculate thermodynamic properties of species and reactions
+* `affinity()`: Calculate affinities of formation reactions
+* `equilibrate()`: Calculate equilibrium chemical activities
+* `diagram()`: Plot the results
-## Organization of major functions
+### Functions With Side Effects (Modify System State)
-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 diagram in @Dic08.
-Rectangles and ellipses represent functions and datasets; bold text indicates primary functions.
+* `basis()`: Set basis species and their chemical activities
+* `species()`: Set species of interest and their activities
+* `reset()`: Reset database and system settings to defaults
-<!-- https://stackoverflow.com/questions/14675913/how-to-change-image-size-markdown -->
-{ width=75% }
+## Querying the Thermodynamic Database
-Many functions in CHNOSZ have no side effects.
-That is, the function only returns a result; to use the result elsewhere, it can be assigned to a variable with `<-`.
-In this document, the names of these functions are shown in <span style="color:green">green text</span> (not applicable to the code chunks).
-```{marginfigure}
-When they are mentioned, names of functions in the base and recommended packages of R are said to belong to R.
-Example: Use R's `plot()` to plot the data.
-```
-Major functions without side effects in CHNOSZ are:
+The `info()` function provides access to the OBIGT thermodynamic database:
-* <span style="color:green">`info()`</span>: search for species in the thermodynamic database;
-* <span style="color:green">`subcrt()`</span>: calculate the thermodynamic properties of species and reactions;
-* <span style="color:green">`affinity()`</span>: calculate the affinities of formation reactions using given chemical activities;
-* <span style="color:green">`equilibrate()`</span>: calculate the equilibrium chemical activities of the species of interest;
-* <span style="color:green">`diagram()`</span>: plot the results.
-
-Some functions in CHNOSZ do have side effects: they modify the `thermo` data object in the current R session.
-In this document, the names of these functions are shown in <span style="color:red">red text</span> (but not in the code chunks).
-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">`reset()`</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 species and variables:
-```{r pseudocode, eval=FALSE}
-basis(...)
-species(...)
-a <- affinity(...)
-e <- equilibrate(a) ## optional
-diagram(e) ## or diagram(a)
-reset() ## clear settings for next calculation
-```
-
-# Querying the thermodynamic database
-
-## The <span style="color:green">`info()`</span> function
-
-<span style="color:green">`info()`</span> provides an interface to the OBIGT thermodynamic database that is packaged with CHNOSZ.
-Suppose you are interested in the thermodynamic properties of aqueous methane.
-Because the database is assembled with aqueous species first, they take precedence over other states.
-Searching by chemical formula alone gives the first matching species, in this case aqueous methane:
```{r info_CH4}
+# Get database index for aqueous methane
info("CH4")
```
-The number that is returned is the species index in the database.
-A second argument can be used to specify a physical state with lower precedence:
-```{r info_CH4_gas}
-info("CH4", "gas")
+```{r info_info_CH4}
+# Get thermodynamic properties for methane
+info(info("CH4"))
```
-While some species are identified only by chemical formula, others have distinct names (in English) listed in the database.
-For `r ch4` and inorganic substances that are represented by both gaseous and aqueous forms, the name is applied only to the gas.
-However, the names of organic substances other than methane are applied to aqueous species, which have precedence, and those in other states.
-The following commands get the species indices for some common gases:
-```{r info_names_gas}
-info("methane")
-info("oxygen")
-info("carbon dioxide")
-```
+You can access fuzzy search functionality by using partial names:
-A special case is sulfur; the name refers to both the native mineral, which has precedence, and the gas.
-These two phases can be identifed with the formulas S<sub>2</sub> and S, respectively.
-```{r info_S_S2}
-info("S")
-info("S2")
+```{r info_fuzzy}
+# Search for ribose-related entries
+info("ribose+")
```
-Taking the species number of aqueous methane returned by <span style="color:green">`info()`</span>, use the function again to retrieve the set of standard molal thermodynamic properties and equations of state parameters:
-```{r iCH4, message=FALSE}
-iCH4 <- info("CH4")
-info(iCH4)
-```
+## Calculating Thermodynamic Properties
-Liquid water is species number 1; it has NA entries in the database because dedicated functions are used to compute its properties:
-```{r info_info_water}
-info(info("water"))
-```
+The `subcrt()` function (named after SUPCRT) calculates standard thermodynamic properties:
-## Fuzzy searches
-
-Calling <span style="color:green">`info()`</span> 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 subcrt_CH4}
+# Properties of aqueous methane at default T and P
+subcrt("CH4")
```
-```{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.
-Since "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")
+<p></p>
+```{r subcrt_TP}
+# Custom T,P grid for water in supercritical region
+subcrt("H2O", T = c(400, 500), P = c(250, 300))
```
-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.
+Unit conversions are handled by `T.units()`, `P.units()`, and `E.units()`:
-## Counting elements, chemical formulas, <span style="color:green">`ZC()`</span>
-
-Continuing with the example of aqueous methane, let's look at its chemical formula:
-```{r info_CH4_formula, message=FALSE}
-info(iCH4)$formula
-```
-
-We can use <span style="color:green">`makeup()`</span> to count the elements in the formula, followed by <span style="color:green">`as.chemical.formula()`</span> to rewrite the formula on one line:
-```{r makeup_iCH4}
-makeup(iCH4)
-as.chemical.formula(makeup(iCH4))
-```
-
-For organic species, a calculation of the average oxidation state of carbon (`r zc`) is possible given the species index, chemical formula, or elemental count:
-```{r ZC_iCH4, message=FALSE}
-ZC(iCH4)
-ZC(info(iCH4)$formula)
-ZC(makeup(iCH4))
-```
-
-# Calculating thermodynamic properties
-
-To calculate the standard molal properties of species and reactions, use <span style="color:green">`subcrt()`</span>.
-The name of this function is derived from the SUPCRT package [@JOH92], which is also the source of the Fortran subroutine used to calculate the thermodynamic properties of H<sub>2</sub>O.
-
-If no reaction coefficients are given, <span style="color:green">`subcrt()`</span> calculates the standard molal properties of individual species:
-```{r subcrt_water}
-subcrt("water")
-```
-
-That uses the default temperature and pressure settings, i.e. equally spaced temperature intervals from 0 to 350 °C at *P*<sub>sat</sub>, i.e. 1 bar below 100 °C, or the pressure of liquid-vapor saturation (i.e. boiling) 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](#properties-of-reactions)), 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` and `Cp`), and cm<sup>3</sup> mol<sup>-1</sup> (`V`).
-```
-
-A custom temperature-pressure grid can be specified.
-Here, we calculate the properties of `r h2o` on a *T*, *P* grid in the supercritical region, with conditions grouped by pressure:
-```{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, out.width="100%", echo=FALSE, message=FALSE, fig.cap="Isothermal contours of density (g cm<sup>-3</sup>) and pressure (bar) of water.", cache=TRUE, pngquant=pngquant, timeit=timeit}
-sres <- subcrt("water", T=seq(0,1000,100), P=c(NA, seq(1,500,1)), grid="T")
-water <- sres$out$water
-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}
-```
-
-## Changing units
-
-The default units of temperature, pressure, and energy in the input and output of <span style="color:green">`subcrt()`</span> are °C, bar, and Joules.
-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 <span style="color:green">`subcrt()`</span> and some other functions in CHNOSZ.
-What is the Gibbs energy in J/mol (the default) and cal/mol of aqueous methane at 298.15 K and 0.1 MPa?
-```{r units_CH4}
-T.units("K")
-P.units("MPa")
-subcrt("CH4", T = 298.15, P = 0.1)$out$CH4$G
+```{r E.units}
+# Convert energy units from Joules to calories
E.units("cal")
-subcrt("CH4", T = 298.15, P = 0.1)$out$CH4$G
+subcrt("CH4", T = 25)
+reset() # Restore defaults
```
-The parameters for each species in the OBIGT database have the units indicated by the `E_units` column there, which is independent of the setting of <span style="color:red">`E.units()`</span>.
-Unlike <span style="color:green">`subcrt()`</span>, <span style="color:green">`info()`</span> always displays the database values in the units given in the database.
-A separate function, <span style="color:green">`convert()`</span>, can be used to convert values to selected units.
-The following code shows that the database parameters for CH<sub>4</sub>(aq) are in calories, then converts the standard Gibbs energy from cal/mol to J/mol.
-```{r convert_G, message=FALSE}
-(CH4dat <- info(info("CH4")))
-convert(CH4dat$G, "J")
-```
+## Working with Reactions
-Note that converting the database value to J/mol gives the same result as running `subcrt("CH4", T = 25)` with the default <span style="color:red">`E.units()`</span> setting of J.
+Define reactions with species names, states (optional), and coefficients:
-Use <span style="color:red">`reset()`</span> to restore the units and all other settings for CHNOSZ to their defaults:
-```{r reset}
+```{r subcrt_CO2_reaction}
+# CO2 dissolution reaction
+subcrt(c("CO2", "CO2"), c("gas", "aq"), c(-1, 1), T = 25)
```
-# Properties of reactions
+Reactions can be automatically balanced using basis species:
-## 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 <span style="color:green">`subcrt()`</span>.
-Here we calculate properties for the dissolution of CO<sub>2</sub>.
-This doesn't correspond to the _solubility_ of CO<sub>2</sub>; see the [solubility section](#complete-equilibrium-solubility) in this vignette and [<span style="color:blue">`demo(solubility)`</span>](../demo) for examples of solubility calculations.
-```{r subcrt_CO2}
-subcrt(c("CO2", "CO2"), c("gas", "aq"), c(-1, 1), T = seq(0, 250, 50))
+```{r basis}
+basis(c("CO2", "H2O", "H+", "e-"))
+species(c("CH4", "acetate"))
+subcrt(c("CH4", "acetate"), c("aq", "aq"), c(1, -1), T = 25)
```
-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:
+## Chemical Affinity and Stability Diagrams
-```{r CO2_logK, 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 figure-setup, include = FALSE}
+knitr::opts_chunk$set(
+ fig.margin=TRUE, fig.width=6, fig.height=4, out.width="100%", results="hide",
+ message=FALSE, cache=TRUE, pngquant=pngquant,
+ # Set dpi 20231129
+ dpi = if(nzchar(Sys.getenv("CHNOSZ_BUILD_LARGE_VIGNETTES"))) 72 else 50
+)
```
-```{r CO2_plot, fig.margin=TRUE, fig.width=4, fig.height=4, small.mar=TRUE, out.width="100%", echo=FALSE, message=FALSE, fig.cap="Calculated equilibrium constants for dissolution of CO<sub>2</sub>, CO, and CH<sub>4</sub>.", cache=TRUE, pngquant=pngquant, timeit=timeit}
-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(240, -2.37, expr.species("CO"))
-text(300, -2.57, expr.species("CH4"))
-```
-```{r CO2_logK, eval=FALSE}
-```
-Now we can make the plot, using R's `matplot()`.
-Here, <span style="color:green">`axis.label()`</span> and <span style="color:green">`expr.species()`</span> are used to create formatted axis labels and chemical formulas:
-```{r CO2_plot, eval=FALSE}
-```
+The `affinity()` function calculates chemical affinities over ranges of T, P, and activities:
-## Unbalanced reactions
-
-A balanced chemical reaction conserves mass.
-<span style="color:green">`subcrt()`</span> 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 linearly combine to give the composition of any species in the system.
-The basis species are similar to thermodynamic components, but can include charged species.
-Basis species are used in CHNOSZ to automatically balance reactions; they are 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.
-
-## Automatically balancing reactions
-
-Methanogenic metabolism in reducing environments may proceed by 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", "CH4"), c(-1, 1))$reaction
-```
-
-We can similarly consider reactions for hydrogenotrophic methanogenesis as well as acetate oxidation (without production of methane):
-```{r subcrt_methanogenesis, message=FALSE}
-acetate_oxidation <- subcrt("acetate", -1)
-hydrogenotrophic <- subcrt("CH4", 1)
-acetoclastic <- subcrt(c("acetate", "CH4"), c(-1, 1))
-```
-
-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, out.width="100%", pngquant=pngquant, timeit=timeit}
-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
-
-Usually, <span style="color:green">`subcrt()`</span> 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 <span style="color:green">`subcrt()`</span>.
-
-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 change the states of CO<sub>2</sub> and H<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 <span style="color:green">`subcrt()`</span>:
-```{r basis_mayumi, message=FALSE, results="hide"}
-basis(c("CO2", "H2", "H2O", "H+"))
-basis(c("CO2", "H2"), "gas")
-basis(c("H2", "pH"), c(-3.92, 7.3))
-```
-```{r affinity_acetoclastic, message=FALSE}
-subcrt(c("acetate", "CH4"), 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("CH4", 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", "CH4"), coeffs,
- c("aq", "gas"), logact = c(-3.4, -0.18), T = 55, P = 50)$out
-}
-```
-
-Now we're ready to calculate and plot the affinities.
-Here, we use R's `lapply()` 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.
-R's `do.call()` and `rbind()` are used to turn the list into a data frame that can be plotted with R's `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, out.width="100%", echo=FALSE, message=FALSE, fig.cap="Gibbs energies of acetate oxidation and methanogenesis (after Mayumi et al., 2013).", cache=TRUE, pngquant=pngquant, timeit=timeit}
-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, eval=FALSE}
-```
-
-Let's not forget to clear the system settings, which were modified by <span style="color:red">`basis()`</span>, before running other calculations:
-```{r reset, message=FALSE}
-```
-
-# Using <span style="color:green">`affinity()`</span>
-
-<span style="color:green">`affinity()`</span> offers calculations of chemical affinity of formation reactions over a configurable range of *T*, *P*, and activities of basis species.
-
-## From affinity to diagrams
-
-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 *formed species*.
-Let's consider the stoichiometry of some aqueous sulfur-bearing species.
-Here we use <span style="color:red">`basis()`</span> with a keyword to load a preset basis definition.
-To use Eh as a variable, the electron (*e*<sup>-</sup>) should be in the basis, so we use the keyword (<span style="color:red">`basis("CHNOSe")`</span>).
-```{marginfigure}
-Some available keywords are `CHNOS` (including CO<sub>2</sub>, H<sub>2</sub>O, NH<sub>3</sub>, H<sub>2</sub>S, and O<sub>2</sub>), `CHNOS+` (also including H<sup>+</sup>), and `CHNOSe` (including H<sup>+</sup>, and *e*<sup>-</sup> instead of O<sub>2</sub>).
-See <span style="color:blue">`?basis`</span> for more options.
-```
-```{marginfigure}
-What is `SO42-`? Is it 1 S, 4 O, and 2 negative charges, or 1 S, 42 O, and 1 negative charge?
-The ambiguity of a digit that could belong to the coefficient for the following charge or to that for the preceding element is why formulas in CHNOSZ are written with the number of charges after the + or - symbol.
-`SO4-2` is unambiguously parsed as 1 S, 4 O and 2 negative charges.
-```
-```{r basis_CHNOSZ, results="hide"}
+```{r diagram, fig.cap = "Eh-pH (Pourbaix) diagram for S"}
+# Set up the C-H-N-O-S basis system with electron
basis("CHNOSe")
+# Define aqueous sulfur species
+species(c("SO4-2", "HSO4-", "HS-", "H2S"))
+# Calculate affinities on an Eh-pH grid
+a <- affinity(pH = c(0, 14), Eh = c(-0.8, 1.2))
+# Create an Eh-pH (Pourbaix) diagram
+diagram(a, fill = "terrain")
```
-```{r species_sulfur}
-species(c("H2S", "HS-", "HSO4-", "SO4-2"))
-```
-Aqueous species are assigned default activities of 10<sup>-3</sup> (`logact` is -3).
-Now, we can use <span style="color:green">`affinity()`</span> to calculate the affinities of the formation reactions of each of the species.
-R's `unlist()` is used here to turn the list of values of affinity into a numeric object that can be printed in a couple of lines (note that the names correspond to `ispecies` above):
-```{r affinity}
-unlist(affinity()$values)
-```
+For more sophisticated diagrams involving speciation of basis species, use the `mosaic()` function:
-The values returned by <span style="color:green">`affinity()`</span> are dimensionless, i.e. *A*/(2.303*RT*).
-Although <span style="color:green">`subcrt()`</span> can also calculate affinities (in units of cal/mol or J/mol), the main advantage of <span style="color:green">`affinity()`</span> is that it can perform calculations on an arbitrary grid of *T*, *P*, or activities of basis species.
-This is the foundation for making many types of diagrams that are useful in geochemistry.
-
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/chnosz -r 876
More information about the CHNOSZ-commits
mailing list