[CHNOSZ-commits] r897 - in pkg/CHNOSZ: . R man vignettes
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue May 13 12:58:51 CEST 2025
Author: jedick
Date: 2025-05-13 12:58:50 +0200 (Tue, 13 May 2025)
New Revision: 897
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/R/AD.R
pkg/CHNOSZ/R/water.R
pkg/CHNOSZ/man/ionize.aa.Rd
pkg/CHNOSZ/man/subcrt.Rd
pkg/CHNOSZ/man/water.Rd
pkg/CHNOSZ/vignettes/CHNOSZ.dia
pkg/CHNOSZ/vignettes/CHNOSZ.png
pkg/CHNOSZ/vignettes/anintro.Rmd
pkg/CHNOSZ/vignettes/postprocess.sh
Log:
Update definition of Psat
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2025-05-11 11:20:20 UTC (rev 896)
+++ pkg/CHNOSZ/DESCRIPTION 2025-05-13 10:58:50 UTC (rev 897)
@@ -1,6 +1,6 @@
-Date: 2025-05-11
+Date: 2025-05-13
Package: CHNOSZ
-Version: 2.1.0-68
+Version: 2.1.0-69
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/R/AD.R
===================================================================
--- pkg/CHNOSZ/R/AD.R 2025-05-11 11:20:20 UTC (rev 896)
+++ pkg/CHNOSZ/R/AD.R 2025-05-13 10:58:50 UTC (rev 897)
@@ -120,7 +120,7 @@
GH2O_P <- water("G", T = T, P = P)$G
GH2O_1 <- water("G", T = T, P = 1)$G
f1 <- exp ( (GH2O_P - GH2O_1) / (8.31441 * T) )
- # For Psat, calculate the real liquid-vapor curve (not 1 bar below 100 degC)
+ # For Psat, calculate the real liquid-vapor curve (not floored at 1 bar)
if(isPsat) {
P <- water("Psat", T = T, P = "Psat", P1 = FALSE)$Psat
f1[P < 1] <- P[P < 1]
Modified: pkg/CHNOSZ/R/water.R
===================================================================
--- pkg/CHNOSZ/R/water.R 2025-05-11 11:20:20 UTC (rev 896)
+++ pkg/CHNOSZ/R/water.R 2025-05-13 10:58:50 UTC (rev 897)
@@ -113,7 +113,8 @@
if(identical(P, "Psat")) {
w.P <- H2O[[2]][2]
w.P[w.P == 0] <- NA
- # Psat specifies P = 1 below 100 degC
+ # By default, the value of Psat is floored at 1 bar; P1 = FALSE means
+ # that we want the actual values of vapor-liquid saturation pressure
if(P1) w.P[w.P < 1] <- 1
P.out[i] <- w.P
}
Modified: pkg/CHNOSZ/man/ionize.aa.Rd
===================================================================
--- pkg/CHNOSZ/man/ionize.aa.Rd 2025-05-11 11:20:20 UTC (rev 896)
+++ pkg/CHNOSZ/man/ionize.aa.Rd 2025-05-13 10:58:50 UTC (rev 897)
@@ -15,7 +15,7 @@
\item{aa}{data frame, amino acid composition in the format of \code{thermo()$protein}}
\item{property}{character, property to calculate}
\item{T}{numeric, temperature in \degC}
- \item{P}{numeric, pressure in bar, or \samp{Psat} for vapor pressure of \H2O above 100 \degC}
+ \item{P}{numeric, pressure in bar, or \samp{Psat} for the greater of 1 bar or the vapor-liquid saturation pressure for \H2O}
\item{pH}{numeric, pH}
\item{ret.val}{character, return the indicated value from intermediate calculations}
\item{suppress.Cys}{logical, suppress (ignore) the ionization of the cysteine groups?}
Modified: pkg/CHNOSZ/man/subcrt.Rd
===================================================================
--- pkg/CHNOSZ/man/subcrt.Rd 2025-05-11 11:20:20 UTC (rev 896)
+++ pkg/CHNOSZ/man/subcrt.Rd 2025-05-13 10:58:50 UTC (rev 897)
@@ -169,16 +169,15 @@
# Compare calculated values of heat capacity of iron with
# values from Robie and Hemingway, 1995
T.units("K")
-# We set pressure here otherwise subcrt uses Psat (saturation
-# vapor pressure of H2O above 100 degrees C) which can't be
+# We need to set a finite value for P because Psat can't be
# calculated above the critical point of H2O (~647 K)
-s <- subcrt("Fe", T=seq(300, 1800, 20), P=1)
-plot(s$out[[1]]$T, s$out[[1]]$Cp, type="l",
- xlab=axis.label("T"), ylab=axis.label("Cp"))
+s <- subcrt("Fe", T = seq(300, 1800, 20), P = 1)
+plot(s$out[[1]]$T, s$out[[1]]$Cp, type = "l",
+ xlab = axis.label("T"), ylab = axis.label("Cp"))
# Add points from RH95
-RH95 <- read.csv(system.file("extdata/misc/RH95.csv", package="CHNOSZ"))
+RH95 <- read.csv(system.file("extdata/misc/RH95.csv", package = "CHNOSZ"))
points(RH95[,1], RH95[,2])
-title(main=paste("Heat capacity of Fe(cr)\n",
+title(main = paste("Heat capacity of Fe(cr)\n",
"(points - Robie and Hemingway, 1995)"))
# Reset the units to default values
T.units("C")
Modified: pkg/CHNOSZ/man/water.Rd
===================================================================
--- pkg/CHNOSZ/man/water.Rd 2025-05-11 11:20:20 UTC (rev 896)
+++ pkg/CHNOSZ/man/water.Rd 2025-05-13 10:58:50 UTC (rev 897)
@@ -19,8 +19,8 @@
\arguments{
\item{property}{character, computational setting or property(s) to calculate}
\item{T}{numeric, temperature (K)}
- \item{P}{numeric, pressure (bar), or \samp{Psat} for vapor-liquid saturation}
- \item{P1}{logical, output pressure of 1 bar below 100 \degC instead of calculated values of \samp{Psat}?}
+ \item{P}{numeric, pressure (bar), or \samp{Psat} for the greater of 1 bar or the vapor-liquid saturation pressure for \H2O (when P1 is TRUE)}
+ \item{P1}{logical, if TRUE, uses the greater of 1 bar or the vapor-liquid saturation pressure for \samp{Psat}; if FALSE, uses the vapor-liquid saturation pressure for \samp{Psat}}
}
\details{
@@ -92,7 +92,7 @@
\code{water.SUPCRT92} interfaces to the FORTRAN subroutine taken from the \acronym{SUPCRT92} package (H2O92D.F) for calculating properties of water.
These calculations are based on data and equations of Levelt-Sengers et al., 1983, Haar et al., 1984, and Johnson and Norton, 1991, among others (see Johnson et al., 1992).
-A value of \code{P} set to \samp{Psat} refers to one bar below 100 \degC, otherwise to the vapor-liquid saturation pressure at temperatures below the critical point (\samp{Psat} is not available at temperatures above the critical point).
+A value of \samp{Psat} for \code{P} refers to the greater of 1 bar or the vapor-liquid saturation pressure of H2O (\samp{Psat} is not defined at temperatures above the critical point).
\code{water.SUPCRT92} provides a limited interface to the FORTRAN subroutine; some functions provided there are not made available here (e.g., using variable density instead of pressure, or calculating the properties of steam).
The stated temperature limits of validity of calculations in \code{water.SUPCRT92} are from the greater of 0 \degC or the melting temperature at pressure, to 2250 \degC (Johnson et al., 1992).
Modified: pkg/CHNOSZ/vignettes/CHNOSZ.dia
===================================================================
(Binary files differ)
Modified: pkg/CHNOSZ/vignettes/CHNOSZ.png
===================================================================
(Binary files differ)
Modified: pkg/CHNOSZ/vignettes/anintro.Rmd
===================================================================
--- pkg/CHNOSZ/vignettes/anintro.Rmd 2025-05-11 11:20:20 UTC (rev 896)
+++ pkg/CHNOSZ/vignettes/anintro.Rmd 2025-05-13 10:58:50 UTC (rev 897)
@@ -151,16 +151,20 @@
### Calculating thermodynamic properties
-The `subcrt()` function [loosely named after SUPCRT; @JOH92] calculates standard thermodynamic properties:
+The `subcrt()` function [loosely named after SUPCRT; @JOH92] calculates standard thermodynamic properties.
+The default conditions are 0.01-350 °C along the Psat curve (here, Psat is the greater of 1 bar or the vapor-liquid saturation pressure for H2O):
+
```{r subcrt_CH4}
# Properties of aqueous methane at default T and P
subcrt("CH4")
```
+You can customize the T-P grid by passing the appropriate arguments:
+
<p></p>
```{r subcrt_TP}
-# Custom T,P grid for water in supercritical region
+# Custom T, P grid for water in supercritical region
subcrt("H2O", T = c(400, 500), P = c(250, 300))
```
@@ -223,6 +227,10 @@
a <- affinity(pH = c(0, 14), Eh = c(-0.8, 1.2))
# Create an Eh-pH (Pourbaix) diagram
diagram(a, col = 4, col.names = 4, italic = TRUE)
+# Add legend
+TC <- convert(a$T, "C")
+legend <- c(describe.property("T", TC), quote(italic(P) == "1 bar"))
+legend("topright", legend = legend, bty = "n")
```
NOTE: `diagram()` automatically adds shading to regions of water instability with respect to O2 or H2.
@@ -238,9 +246,13 @@
species(c("chalcocite", "tenorite", "cuprite", "copper"), add = TRUE)
bases <- c("H2S", "HS-", "HSO4-", "SO4-2")
m <- mosaic(bases, pH = c(0, 12), Eh = c(-1, 1), T = 200)
-diagram(m$A.species, lwd = 2, bold = species()$state == "cr")
+a <- m$A.species
+diagram(a, lwd = 2, bold = species()$state == "cr")
diagram(m$A.bases, add = TRUE, col = 4, col.names = 4, italic = TRUE)
water.lines(m$A.species, lty = 3, col = 8)
+# Add legend
+legend <- lTP(convert(a$T, "C"), a$P)
+legend("topright", legend = legend, bty = "n")
```
Here we've added dotted lines to help visualize the water stability limits.
@@ -261,6 +273,9 @@
a <- affinity(pH = c(4, 12), T = 100)
e <- equilibrate(a)
diagram(e, alpha = TRUE, add = TRUE, col = 2, names = NA)
+# Add legend
+legend <- as.expression(list(lT(25), lT(100)))
+legend("left", legend = legend, lty = 1, col = c(1, 2))
```
Calculate solubility of minerals or gases:
@@ -308,12 +323,12 @@
thermo.refs()
```
-For citing CHNOSZ itself, see "How should CHNOSZ be cited?" in the [FAQ](FAQ.html).
+For citing CHNOSZ itself, see "[How should CHNOSZ be cited?](FAQ.html#how-should-chnosz-be-cited)" in the FAQ.
## Interlude: Affinity, Formation Reactions, and Balancing
The idea for creating stability diagrams in CHNOSZ came from Prof. Harold Helgeson's geochemistry course at UC Berkeley.
-There, we constructed diagrams that were "balanced on" a metal.
+There, the students constructed diagrams that were "balanced on" a metal.
For instance, in a system balanced on aluminum, Al is only present in the minerals on both sides of the reaction and is not free as an ion.
The reaction-based method, used for making diagrams by hand, looks at reactions between pairs of species (let's call them transformation reactions),
@@ -322,15 +337,15 @@
then selects the most stable species according to their affinity values.
Affinity is just the opposite of non-standard Gibbs energy of reaction.
-"Standard Gibbs energy of reaction" and "Gibbs energy of reaction" - which are two different things - have unfortunately similar names except for
-an optional "overall" or "non-standard" in front of the latter [the word choice varies among authors, e.g. @AL19;@STK19].
-"Non-standard Gibbs energy of reaction" doesn't lend itself to a short, unambiguous function name, which is why its opposite, "affinity", is used in CHNOSZ.
+"Standard Gibbs energy of reaction" and "Gibbs energy of reaction" – which are different concepts – have unfortunately similar names except for
+an optional *overall* or *non-standard* in front of the latter [the word choice varies among authors, e.g. @AL19;@STK19].
+"Non-standard Gibbs energy of reaction" doesn't lend itself to a short, unambiguous function name, which is why its opposite, *affinity*, is used in CHNOSZ.
-In the reaction-based method, transformation reactions are said to be "balanced on" a metal.
-The grid-based method implements this balancing constraint by dividing the affinities of formation reactions by the coefficients of a basis species.
-CHNOSZ uses these normalized affinities for making relative stability diagrams, referred to as the *maximum affinity method*.
-Both the reaction- and grid-based methods have the same limitation: every candidate species must have non-zero stoichiometry of a given metal
-(or of a basis species with that metal).
+In the reaction-based method, transformation reactions are said to be *balanced on* a metal.
+The grid-based method implements this balancing constraint by dividing the affinities of formation reactions by the stoichiometric coefficients of a basis species.
+CHNOSZ uses these normalized affinities for making relative stability diagrams, a technique referred to as the maximum affinity method [@Dic19].
+The reaction- and grid-based methods both have the same limitation: every species considered in the relative stability calculation
+must have non-zero stoichiometry of the metal the transformation reactions are balanced on (or equivalently of the *conserved* basis species that has that metal).
## Advanced Uses
@@ -339,11 +354,8 @@
### 1. Use helper functions to create formatted labels for diagrams
Labeling diagrams is an important but often difficult step for creating publication-ready figures.
-For general information on formatting mathematical expressions in R, see the documentation for `plotmath()` as well as `bquote()`,
-which allows substituting variables into an equation.
-CHNOSZ has several helper functions for creating labels.
-`axis.label()` and `expr.species()` are used to create formatted axis labels and chemical formulas.
+CHNOSZ provides the `axis.label()` and `expr.species()` functions to create formatted axis labels and chemical formulas.
Let's revisit the CO2 dissolution example seen earlier and add two other gases (carbon monoxide and methane).
This plot is similar to Figure 18 of @MSS13.
@@ -359,11 +371,25 @@
text(80, -1.7, expr.species("CO2"))
text(240, -2.37, expr.species("CO"))
text(300, -2.57, expr.species("CH4"))
+# Add legend
+legend <- describe.property("P", "Psat")
+legend("topright", legend = legend, bty = "n")
```
-See the help pages in CHNOSZ for additional functions for labeling diagrams, including `describe.reaction()` to format a chemical reaction
-from the output of `subcrt()`, and `lT()` and related functions for compact representations of temperature and other variables for plot legends.
+CHNOSZ has some other helper functions for labeling diagrams:
+- `describe.reaction()` to format chemical reactions from the output of `subcrt()`
+- `describe.property()` to format property values as equations (e.g. "*T* = 25 °C")
+- `describe.basis()` to format the activities of basis species
+- `lT()`, `lP()`, `lTP()`, and related functions for more compact representations of conditions (e.g. "25 °C, 1 bar")
+
+There is no single best approach to formatting text for legends, and sometimes it's easiest to use basic R functions:
+
+- See the documentation for `plotmath()` for general information on formatting mathematical expressions in R.
+- `bquote()` is useful for putting values of variables into expressions.
+- Some of the examples in this vignette and in the demos use `as.expression()` on the combined output from other functions
+ to produce complex legends for the `legend()` function.
+
### 2. Use <span style="color:green;font-family:monospace;">retrieve()</span> to search species by elements
Want to find all the Al complexes in the database?
@@ -400,7 +426,7 @@
info(setdiff(iaq_all, iaq))
```
-The convention for SUPCRT-family databases is to use anhydrous species.
+NOTE: The convention for SUPCRT-family databases is to use anhydrous species.
For example, AlO2- in SLOP98 corresponds to Al(OH)4- in the default database (see output above).
They are effectively the same species, which is why only the latter [taken from a more extensive compilation for Al species properties; @TS01]
is used in the default database.
@@ -435,12 +461,12 @@
reset()
```
-### 5. Use <span style="color:red;font-family:monospace;">basis()</span> species to select compositional variables to plot
+### 5. Use <span style="color:red;font-family:monospace;">basis()</span> species to define the compositional space
A common question is: what are the basis species used for?
-The basis species define the compositional variables that can be added to a diagram.
-In more precise terms, they define the thermodynamic components of a chemical system.
-The composition of any possible species in that system can be represented by a linear combination of the basis species.
+The basis species define the compositional variables of a system.
+The composition of any species that can possibly be formed in that system can be represented by a linear combination of the basis species.
+The basis species may also be referred to as thermodynamic components, although a strict definition of the latter does not include charged species.
CHNOSZ requires that the number of basis species is equal to the number of different elements in the basis species (plus charge, if present).
If you were studying the relative stability of F- and OH-complexes with Al, you might be tempted to try this basis definition:
@@ -475,6 +501,7 @@
Selected basis species can be assigned to plot axes (with a range of values) in `affinity()`.
NOTE: `logact` is the logarithm of *activity* for aqueous species, solids, and liquids, or logarithm of *fugacity* for gases.
+All logarithms in CHNOSZ are common logarithms (R function: `log10()`) unless indicated otherwise.
How about the formed species in the system - that is, the species whose stability fields we want to visualize?
We both list the species and set their activities using `species()`.
@@ -535,7 +562,7 @@
The number of such arguments is the dimensionality of the final plot.
The grid resolution (`res`) defaults to 256 and can be different for each variable.
The names of the variables can be the formulas of any of the basis species, or `T`, `P`, or `IS` for temperature, pressure, or ionic strength.
-These last three default to 25 °C, `Psat` (1 bar below 100 °C and saturation pressure at higher temperatures), and 0 mol/kg.
+These last three default to 25 °C, `Psat` (the greater of 1 bar or the vapor-liquid saturation pressure for H2O), and 0 mol/kg.
I often start with a low grid resolution to quickly iterate a calculation, then switch to a higher resolution when I'm satisfied with the result.
@@ -580,10 +607,12 @@
The example here uses Fe; try changing it to Cu, Zn, Pb, Au, or something else!
CHNOSZ generates warning messages about being above the Cp limits for various iron oxyhydroxides.
-If you see warning messages like this, it's a good idea not to ignore them; instead, consider whether you might be pushing extrapolations of the Cp equation too far.
For the present calculation, the warnings are probably harmless because the predicted set of stable minerals (pyrite, pyrrhotite, magnetite, and hematite)
is consistent with many published diagrams.
+NOTE: If you see warning messages like this, it's a good idea not to ignore them,
+but to consider whether you might be pushing the extrapolations of the Cp equation too far.
+
```{r solubility, echo=FALSE, fig.cap="Mineral stability diagram; aqueous species predominance diagram; composite diagram with one solubility contour; diagram with multiple solubility contours in units of log *m*", fig.margin=FALSE, fig.fullwidth=TRUE, fig.width=16, fig.height=3, cache=TRUE, dpi=72}
par(mfrow = c(1, 4))
@@ -955,7 +984,7 @@
- `linesout`: x, y coordinates of boundary lines between stability fields
NOTE: If `diagram()` was passed the output of `equilibrate()` or `solubility()`,
-then its output contains logarithms of activities instead of dimensionless affinities.
+then its output contains base-10 logarithms of activities instead of dimensionless affinities.
### 17. Writing chemical formulas and counting and summing elements with <span style="color:green;font-family:monospace;">makeup()</span>
Modified: pkg/CHNOSZ/vignettes/postprocess.sh
===================================================================
--- pkg/CHNOSZ/vignettes/postprocess.sh 2025-05-11 11:20:20 UTC (rev 896)
+++ pkg/CHNOSZ/vignettes/postprocess.sh 2025-05-13 10:58:50 UTC (rev 897)
@@ -44,7 +44,11 @@
sed -i 's/<code>axis.label()/<code><a href="..\/html\/util.expression.html" style="background-image:none;color:green;">axis.label()<\/a>/g' anintro.html
sed -i 's/<code>expr.species()/<code><a href="..\/html\/util.expression.html" style="background-image:none;color:green;">expr.species()<\/a>/g' anintro.html
sed -i 's/<code>describe.reaction()/<code><a href="..\/html\/util.expression.html" style="background-image:none;color:green;">describe.reaction()<\/a>/g' anintro.html
+sed -i 's/<code>describe.property()/<code><a href="..\/html\/util.expression.html" style="background-image:none;color:green;">describe.property()<\/a>/g' anintro.html
+sed -i 's/<code>describe.basis()/<code><a href="..\/html\/util.expression.html" style="background-image:none;color:green;">describe.basis()<\/a>/g' anintro.html
sed -i 's/<code>lT()/<code><a href="..\/html\/util.legend.html" style="background-image:none;color:green;">lT()<\/a>/g' anintro.html
+sed -i 's/<code>lP()/<code><a href="..\/html\/util.legend.html" style="background-image:none;color:green;">lP()<\/a>/g' anintro.html
+sed -i 's/<code>lTP()/<code><a href="..\/html\/util.legend.html" style="background-image:none;color:green;">lTP()<\/a>/g' anintro.html
sed -i 's/<code>thermo.refs()/<code><a href="..\/html\/util.data.html" style="background-image:none;color:green;">thermo.refs()<\/a>/g' anintro.html
sed -i 's/<code>mod.OBIGT()/<code><a href="..\/html\/add.OBIGT.html" style="background-image:none;color:red;">mod.OBIGT()<\/a>/g' anintro.html
sed -i 's/<code>T.units()/<code><a href="..\/html\/util.units.html" style="background-image:none;color:red;">T.units()<\/a>/g' anintro.html
@@ -66,6 +70,9 @@
sed -i 's/<code>help(/<code><a href="..\/..\/utils\/html\/help.html">help(<\/a>/g' anintro.html
sed -i 's/<code>help.start/<code><a href="..\/..\/utils\/html\/help.start.html">help.start<\/a>/g' anintro.html
sed -i 's/<code>maintainer/<code><a href="..\/..\/utils\/html\/maintainer.html">maintainer<\/a>/g' anintro.html
+sed -i 's/<code>as.expression/<code><a href="..\/..\/base\/html\/expression.html">as.expression<\/a>/g' anintro.html
+sed -i 's/<code>legend/<code><a href="..\/..\/graphics\/html\/legend.html">legend<\/a>/g' anintro.html
+sed -i 's/<code>log10()/<code><a href="..\/..\/base\/html\/Log.html">log10()<\/a>/g' anintro.html
### End processing anintro.html
More information about the CHNOSZ-commits
mailing list