CHNOSZ
Thermodynamic Calculations and Diagrams for Geochemistry

Package index
Search the CHNOSZ package
Vignettes
Functions
314
Source code
137
Man pages
46
Home
/
R-Forge
/
CHNOSZ
/
An Introduction to CHNOSZ

An Introduction to CHNOSZ
In CHNOSZ: Thermodynamic Calculations and Diagrams for Geochemistry

```{css, echo=FALSE} html { font-size: 14px; } / zero margin around pre blocks (looks more like R console output) / pre { margin: 0; padding: 0; }

```{js, echo=FALSE}
function ToggleDiv(ID) {
  var D = document.getElementById("D-" + ID);
  var B = document.getElementById("B-" + ID);
  if (D.style.display === "none") {
    // open the div and change button text
    D.style.display = "block";
    B.innerText = "Hide code";
  } else {
    // close the div and change button text
    D.style.display = "none";
    B.innerText = "Show code";
  }
}

library(knitr)

## From "Tufte Handout" example dated 2016-12-27
# Invalidate cache when the tufte version changes
opts_chunk$set(tidy = FALSE, cache.extra = packageVersion('tufte'))
options(htmltools.dir.version = FALSE)

## 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 

## Colorize messages 20171031
## Adapted from https://gist.github.com/yihui/2629886#file-knitr-color-msg-rnw
color_block = function(color) {
  function(x, options) sprintf('<pre style="color:%s">%s</pre>', color, x)
}
knit_hooks$set(warning = color_block('magenta'), error = color_block('red'), message = color_block('blue'))

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.

Getting Started

After installing CHNOSZ from CRAN, load the package:

library(CHNOSZ)

This makes the thermodynamic database and functions available in your R session. To restore default settings at any point, use reset().

Core Functionality

CHNOSZ offers several primary functions for thermodynamic analysis:

Functions Without Side Effects (Return Values)

Functions With Side Effects (Modify System State)

Querying the Thermodynamic Database

The info() function provides access to the OBIGT thermodynamic database.

# Get database index for aqueous methane
info("CH4")

When searching by formula, aqueous species are returned if they are available. Use a species name or add the state to get a particular physical state - aq, cr, gas, or liq:

# Two ways to lookup methane gas
info("methane")
info("CH4", "gas")

Use info() recursively to return thermodynamic parameters:

# Get thermodynamic properties for aqueous methane
info(info("CH4"))

You can access fuzzy search functionality by using partial names:

# Search for ribose-related entries
info("ribose+")

Calculating Thermodynamic Properties

The subcrt() function [named after SUPCRT; @JOH92] calculates standard thermodynamic properties:

# Properties of aqueous methane at default T and P
subcrt("CH4")


# Custom T,P grid for water in supercritical region
subcrt("H2O", T = c(400, 500), P = c(250, 300))

Unit conversions are handled by T.units(), P.units(), and E.units():

# Change energy units from Joules to calories
E.units("cal")
subcrt("CH4", T = 25)
reset()  # Restore defaults

Working with Reactions

Define reactions with species names, states (optional), and coefficients:

# CO2 dissolution reaction
subcrt(c("CO2", "CO2"), c("gas", "aq"), c(-1, 1), T = 25)

Reactions can be automatically balanced using basis species:

basis(c("CO2", "H2O", "H+", "e-"))
species(c("CH4", "acetate"))
subcrt(c("CH4", "acetate"), c("aq", "aq"), c(1, -1), T = 25)

There are some keywords (e.g. CHNOS+, CHNOSe and QEC) for loading predefined sets of basis species. See the help page (?basis) for more information.

Chemical Affinity and Stability Diagrams

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
)

The affinity() function calculates chemical affinities over ranges of T, P, and activities:

# 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, col = 4, col.names = 4, italic = TRUE)

Note that diagram() automatically adds shading to represent where water is not stable with respect to O2 or H2.

For more sophisticated diagrams involving speciation of basis species, use the mosaic() function:

# Create a mosaic diagram for Cu-S-Cl-H2O system
basis(c("Cu", "H2S", "Cl-", "H2O", "H+", "e-"))
basis("H2S", -6)
basis("Cl-", -1)
species(c("CuCl", "CuCl2-", "CuCl3-2", "CuCl+", "CuCl2", "CuCl3-", "CuCl4-2"))
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")
diagram(m$A.bases, add = TRUE, col = 4, col.names = 4, italic = TRUE)
water.lines(m$A.species, lty = 3, col = 8)

Here we've added dotted lines to help visualize the water stability limits.

Equilibrium Calculations

Calculate equilibrium distributions of species:

# Carbonate speciation vs pH
basis(c("CO2", "H2O", "H+", "e-"))
species(c("CO2", "HCO3-", "CO3-2"))
# 25 degrees C
a <- affinity(pH = c(0, 14))
e <- equilibrate(a)
diagram(e, alpha = TRUE)
# 100 degrees C
a <- affinity(pH = c(4, 12), T = 100)
e <- equilibrate(a)
diagram(e, alpha = TRUE, add = TRUE, col = 2, names = NA)

Calculate solubility of minerals or gases:

# Corundum solubility vs pH
basis(c("Al+3", "H2O", "H+", "e-"))
species("corundum")
iaq <- info(c("Al(OH)4-", "Al(OH)3", "Al(OH)2+", "AlOH+2", "Al+3"))
s <- solubility(iaq, pH = c(2, 10), in.terms.of = "Al+3")
diagram(s, col = 3, lwd = 2, ylim = c(-10, -2))
diagram(s, type = "loga.equil", add = TRUE)
legend("topright", c("25 °C", "1 bar"), bty = "n")

Activity Coefficients

Incorporate non-ideal behavior using the extended Debye-Hückel equation by setting the ionic strength parameter IS:

# Corundum solubility again
basis(c("Al+3", "H2O", "H+", "e-"))
species("corundum")
iaq <- info(c("Al(OH)4-", "Al(OH)3", "Al(OH)2+", "AlOH+2", "Al+3"))
# Calculate with ionic strength of 0 and 1 molal
s0 <- solubility(iaq, pH = c(2, 10))
s1 <- solubility(iaq, pH = c(2, 10), IS = 1)
diagram(s1, col = 4, lwd = 3, ylim = c(-8, -2))
diagram(s0, col = 3, lwd = 2, add = TRUE)
legend("topleft", legend = c(1, 0), lwd = c(3, 2), col = c(4, 3), title = "I (mol/kg)", bty = "n")
legend("topright", c("25 °C", "1 bar"), bty = "n")

Functions that have the IS parameter are subcrt(), affinity(), mosaic(), and solubility(). When a value for IS is specified, inputs and outputs labeled as activities are actually interpreted or reported as molalities.

Building Up

Having seen basic examples of the main features of CHNOSZ, here are some ideas for building more complex calculations and diagrams.

1. Use helper functions to create formatted logK and other labels for diagrams

Labeling diagrams is an importan part of creating publication-ready figures, and chemical formulas and reactions can be diffcult for beginners and experienced R coders alike. See the documentation for R's plotmath() for formatting mathematical expressions. My go-to function for building expressions programmatically is bquote(), which allows substituting variables into a formula:

CHNOSZ has several helper functions for creating labels. axis.label() and expr.species() are used 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.

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)
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"))

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.

2. Use retrieve() to search species by elements

Want to find all the Al complexes in the database? List them by calling retrieve() with the main element optionally followed by the ligand elements and state.

# List aqueous Al species in the default database
iaq <- retrieve("Al", ligands = c("F", "H", "O"), state = "aq")
# Print the first few rows and columns
info(iaq)[1:3, 1:5]
# Use the species index or name in subcrt()
subcrt(iaq[1:3], T = 25)
#subcrt(names(iaq)[1:3], T = 25)  # same as above

The result of retrieve() can also be used to add species to a diagram; see the example in #3 below.

basis(c("Al+3", "H2O", "F-", "H+", "e-"))
species(iaq)
species(names(iaq))  # same as above

3. Load optional data with add.OBIGT()

Perhaps you'd like to use data from older databases that have been superseded by later updates. The OBIGT vignette briefly summarizes the superseded data for SUPCRT92 and SLOP98 [@SLOP98]. Use add.OBIGT() to load these old data entries.

add.OBIGT("SLOP98")
iaq_all <- retrieve("Al", ligands = c("F", "H", "O"), state = "aq")
# Use difference between sets to get the added species
info(setdiff(iaq_all, iaq))

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. Unless you have a specific reason to compare them, redundant species should not be used in the same equilibrium calculation.

4. Use OBIGT() or reset() to restore the default database and settings

OBIGT() restores just the database, without affecting other settings (e.g. E.units(), basis() and species()). reset() resets all settings in CHNOSZ, including the database. These functions are useful for both interactive use and scripts that compare different versions of data or plots for different systems or conditions.

Let's put items #1-3 together to remake the corundum solubility plot using only species available in SLOP98. To do this, we use add.OBIGT() followed by retrieve() to gather the species indices for all Al species, then taken only those species sourced from @SSWS97.

# Add superseded species from SLOP98
add.OBIGT("SLOP98")
# List all aqueous Al species
iaq <- retrieve("Al", state = "aq")
# Keep only species from Shock et al. (1997)
iaq <- iaq[grepl("SSWS97", info(iaq)$ref1)]
# Plot corundum solubility vs pH
basis(c("Al+3", "H2O", "H+", "e-"))
species("corundum")
s <- solubility(iaq, pH = c(2, 10), in.terms.of = "Al+3")
## Alternatively, we could use the species names
#s <- solubility(names(iaq), pH = c(2, 10), in.terms.of = "Al+3")
diagram(s, col = 3, lwd = 2, ylim = c(-10, -2))
diagram(s, type = "loga.equil", add = TRUE)
legend("topright", c("25 °C", "1 bar"), bty = "n")
# Reset the database for subsequent calculations
reset()

5. Use basis() species to select compositional variables to plot

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.

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:

basis(c("Al+3", "F-", "OH-"))

According to the message, we don't have enough basis species for the number of elements. Since hydroxide (OH-) is just water minus a proton, we could try this instead:

basis(c("Al+3", "F-", "H+", "H2O"))

That's still not enough species. As is often the case, we need to include a basis species representing oxidation-reduction (redox) reactions, even if there are no redox reactions between the formed species. Here are two possible basis definitions that do not give an error.

# Use "oxygen" to get oxygen gas (for logfO2 diagrams)
basis(c("Al+3", "F-", "H+", "H2O", "oxygen"))
# Use "e-" to get aqueous electron (for Eh diagrams)
basis(c("Al+3", "F-", "H+", "H2O", "e-"))

6. Set activities of formed species() to define a single solubility contour

In order to make a diagram with stability fields for different species, CHNOSZ needs to know about the activities of all the species in the reaction. The activities of the basis species start with constant values as shown in the output above (logact column). Selected basis species can be assigned to plot axes (with a range of values) in affinity().

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(). The function defaults to activities of 1e-3 (logact of -3) for aqueous species and unit activity (logact = 0) for minerals, gases, and liquids. Let's change this to activities of 1e-6 for the formed species.

basis(c("Al+3", "F-", "H+", "H2O", "e-"))
iaq <- retrieve("Al", ligands = c("F", "H", "O"), state = "aq")
# Check that the data are from the same source
stopifnot(all(info(iaq)$ref1 == "TS01"))
species(iaq, -6)

This value for logact defines a solubility contour, as we'll see below.

7. When to use add = TRUE

There are two places where you might see add = TRUE. First, in species() to add species not already in the list. Without add = TRUE, any existing species are discarded. Second, in diagram() to add to an existing plot.

Let's put items #4-7 together to make a Pourbaix (Eh-pH) diagram for Al with two solubility contours.

basis(c("Mn+2", "H+", "H2O", "e-"))
icr <- retrieve("Mn", ligands = c("H", "O"), state = "cr")
iaq <- retrieve("Mn", ligands = c("H", "O"), state = "aq")
# First layer, logact(aq) = -3
species(icr)
species(iaq, add = TRUE)
a <- affinity(pH = c(4, 14), Eh = c(-1, 1), T = 100)
# Use names = NA to avoid plotting labels twice
diagram(a, lty = 2, names = NA)
# Second layer, logact(aq) = -4
species(icr)
species(iaq, -4, add = TRUE)
a <- affinity(pH = c(4, 14), Eh = c(-1, 1), T = 100)
d <- diagram(a, bold = species()$state == "cr", add = TRUE)
# Add water stability limits
water.lines(d, lty = 3, col = 8)
# Add legends
legend("topright", legend = c(lT(100), lP("Psat")), bty = "n")
title = as.expression(quote(log~italic(a)*"Mn(aq)"))
legend("bottomleft", legend = c(-3, -4), lty = c(2, 1), title = title, bty = "n")

The shaded areas in the diagram represent water instability regions and are automatically added by diagram(). We use water.lines() here to plot the water stability limits with dotted lines.

8. Set grid resolution and constant T, P, or IS in affinity()

After defining the basis species and formed species (and their constant activities), you have some choices about what variables to put on the plot, the grid resolution, and values for a few other variables. affinity() accepts one or more named arguments that specifying ranges of variables using the default grid resolution of 256 (c(min, max)) or ranges and a custom grid resolution (c(min, max, res)). 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.

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.

9. Use NaCl() to estimate ionic strength from NaCl concentration

Sodium chloride (NaCl) solutions are commonly used reference points for geochemical models. The NaCl() function provides a quick-and-dirty way to estimate ionic strength and activity of chloride (Cl-) for a given total amount of NaCl added to 1 kg of H2O. These values can then be used in setting up a calculation that involves these variables.

This function does not use either the basis() or species() definitions. The following example runs a calculation for 0.8 mol/kg NaCl and given T, P, and pH. See demo('sum_S') for the fully worked-out example that uses this code [based on a diagram in @SW02].

T <- 300
P <- 1000
pH <- 5
m_NaCl = 0.8
NaCl <- NaCl(m_NaCl = m_NaCl, T = T, P = P, pH = pH)
print(paste("mol NaCl added to 1 kg H2O:", m_NaCl))
print(paste("Ionic strength (mol/kg):", NaCl$IS))
print(paste("Chloride concentration (mol/kg):", NaCl$m_Clminus))

10. Use solubility() to draw multiple solubility contours

There are many uses for "composite diagrams" [@GC65], where stability fields for minerals and predominance fields for aqueous species are both present. As mentioned above (#6), setting the activities of formed aqueous species defines a single solubility contour. This represents a concentration-dependent boundary between minerals and aqueous species on a composite diagram, a concept referred to either as "equisolubility" [@Pou74] or "isosolubility" [@Hel64 and Garrels and Christ, 1965].

Composite diagrams are often drawn with multiple solubility contours in order to show the dependence of solubility on pH, redox, or other variables. See examples of Eh-pH composite diagrams in demo("Pourbaix").

You could loop over constant activities to make a series of solubility contours (see the above example for Mn). An easier solution is to use solubility() to visualize multiple solubility contours in one go. The basic idea is to first load the mineral(s) containing a single metal as the formed species(). Then, list the aqueous species with that metal as the first argument to solubility(). The remaining arguments to the function define the plot variables, just as in affinity() and mosaic().

Let's put together #8-10 to make a set of diagrams for a single metal. The example here uses Fe; try changing it to Cu, Zn, Pb, Au, or something else!

par(mfrow = c(1, 4))

# System definition
metal <- "Fe"
# The concentration to be used for a single solubility contour
logm_metal <- -6
T <- 150
P <- "Psat"
Stot <- 1e-3
wt_percent_NaCl <- 10
# Plot variables
res <- 300
pH <- c(0, 14, res)
O2 <- c(-55, -35, res)

# Work out NaCl molality from weight percent
# Convert to permil to get parts out of 1000 g (1 kg) of solution
wt_permil_NaCl <- wt_percent_NaCl * 10
# Divide by molecular weight to get moles of NaCl in 1000 g of solution
moles_NaCl <- wt_permil_NaCl / mass("NaCl")
# Subtract from 1000 g to get mass of H2O
grams_H2O <- 1000 - wt_permil_NaCl
# This gives the moles of NaCl added to 1 kg of H2O
m_NaCl <- 1000 * moles_NaCl / grams_H2O
# Now calculate ionic strength and molality of Cl-
NaCl_res <- NaCl(m_NaCl, T = T, P = P)
IS = NaCl_res$IS
m_Clminus = NaCl_res$m_Clminus

# Set up basis species
basis(c(metal, "H2S", "Cl-", "oxygen", "H2O", "H+"))
basis("H2S", log10(Stot))
basis("Cl-", log10(m_Clminus))
# Define basis species to change for mosaic calculation
bases <- c("H2S", "HS-", "HSO4-", "SO4-2")

# Git minerals and aqueous species
icr <- retrieve(metal, c("Cl", "S", "O"), state = "cr")
iaq <- retrieve(metal, c("Cl", "S", "O"), state = "aq")

# Make diagram for minerals
species(icr)
mcr <- mosaic(bases, pH = pH, O2 = O2, T = T, P = P, IS = IS)
diagram(mcr$A.species, bold = TRUE, fill = "#FAEBD788")
diagram(mcr$A.bases, lty = 2, col = 4, col.names = 4, italic = TRUE, add = TRUE)
title(paste(metal, "minerals"), font.main = 1)
# Add a legend
leg_list <- c(
  lTP(T, P),
  lNaCl(m_NaCl),
  lS(Stot)
)
leg_expr <- lex(leg_list)
legend("topright", legend = leg_expr, bty = "n")

# Make diagram for aqueous species
species(iaq)
maq <- mosaic(bases, pH = pH, O2 = O2, T = T, P = P, IS = IS)
diagram(maq$A.species, fill = "#F0F8FF88")
title(paste(metal, "aqueous species"), font.main = 1)

# Make diagram for minerals and aqueous species
species(icr)
species(iaq, logm_metal, add = TRUE)
m <- mosaic(bases, pH = pH, O2 = O2, T = T, P = P, IS = IS)
diagram(m$A.species, bold = species()$state == "cr")
diagram(m$A.bases, lty = 2, col = 4, col.names = 4, italic = TRUE, add = TRUE)
main = bquote("One solubility contour for" ~ log ~ italic(m)*.(metal)*"(aq)" == .(logm_metal))
title(main, font.main = 1)

# Make solubility contours
species(icr)
s <- solubility(iaq, bases = bases, pH = pH, O2 = O2, T = T, P = P, IS = IS, in.terms.of = metal)
levels <- seq(-12, 9, 3)
diagram(s, levels = levels, contour.method = "flattest")
title("Multiple solubility contours", font.main = 1)

11. Use convert() for common unit conversions

The convert() function offers several unit conversions. It implements reciprocal conversion between pairs of units, so only the destination unit needs to be specified. Some simple uses are to convert temperature, pressure, or energy values:

# Convert 100 degrees C to value in Kelvin
convert(100, "K")
# Convert 273.15 K to value in degrees C
convert(273.15, "C")
# Convert 1 bar to value in MPa
convert(1, "MPa")
# Convert 100 MPa to value in bar
convert(100, "bar")
# Convert 1 cal/mol to value in J/mol
convert(1, "J")
# Convert 10 J/mol to value in cal/mol
convert(10, "cal")

Another use of convert() is to convert the output of solubility() from log activity to ppm, ppb, log ppm, or log ppb. The following code continues from the example above:

sppm <- convert(s, "ppm")
levels <- c(1e-6, 1e-3, 1e0, 1e3)
diagram(sppm, levels = levels)

12. Use the transect mode of affinity() for synchronized variables

Specify 4 or more values for one or more variables (each variable should have the same number of values, or be set to a constant) to activate the transect mode of affinity(). The transect mode allows defining an arbitrary path in multidimensional space.

Here's a simple example:

basis("CHNOSe")
species(c("NO3-", "NO2-", "NH4+", "NH3"))
a <- affinity(pH = c(6, 8, 6, 8), Eh = c(0.5, 0.5, -0.5, -0.5))

# Print pH and Eh values used for calculation
a$vals
# Print affinity values calculated for each species
a$values

Note that affinity() returns dimensionless values of affinity (i.e., A/2.303RT). Below we'll see how to convert these to energetic units.

13. Choose non-default balancing constraints in diagram()

diagram() looks for the first basis species that has non-zero coefficients in each of the formed species. This is called the conserved basis species in the documentation. The affinity values are then divided by the coefficients of the conserved basis species to give normalized affinities. This is how "balancing on a metal" is implemented.

Let's put together #11-13 to calculate affinities of organic synthesis reactions in mixed seawater and hydrothermal fluid from the Rainbow vent field using speciation results from @SC10:

basis(c("CO2", "H2", "NH4+", "H2O", "H2S", "H+"))
# Constant activity of CH4 is a simplification of the model
species("CH4", -3)
# Add other organic species with lower activity
species(c("propanoic acid", "adenine", "leucine", "deoxyribose"), -6, add = TRUE)
# Read file with variable conditions
file <- system.file("extdata/misc/SC10_Rainbow.csv", package = "CHNOSZ")
# Use `check.names = FALSE` to preserve the `NH4+` column name
df <- read.csv(file, check.names = FALSE)
# Reverse the rows to get increasing T
df <- df[rev(rownames(df)), ]
# Change six variables on a transect
a <- affinity(T = df$T, CO2 = df$CO2, H2 = df$H2,
              `NH4+` = df$`NH4+`, H2S = df$H2S, pH = df$pH)
# Get T in Kelvin to make unit conversions
T <- convert(a$vals[[1]], "K")
# Convert dimensionless affinity to Gibbs energy in units of J/mol
# nb. this is the same as converting logK to ΔG of reaction
a$values <- lapply(a$values, convert, "G", T)
# Convert J/mol to cal/mol
a$values <- lapply(a$values, convert, "cal")
# Convert cal/mol to kcal/mol
a$values <- lapply(a$values, `*`, -0.001)
# Setup figure
par(mfrow = c(1, 3))
par(cex = 0.9)
# First plot: balance on C
ylab = quote(italic(A)~"(kcal/mol) per C")
diagram(a, ylim = c(-20, 40), ylab  = ylab, col = 1:5)
abline(h = 0, lty = 2, lwd = 2)
# Second plot: balance on H2
ylab = quote(italic(A)~"(kcal/mol) per H2")
diagram(a, balance = "H2", ylim = c(-15, 10), ylab = ylab, col = 1:5)
abline(h = 0, lty = 2, lwd = 2)
# Third plot: no balancing
ylab <- quote(italic(A)~"(kcal/mol) per species")
diagram(a, balance = 1, ylim = c(-100, 100), ylab = ylab, col = 1:5)
abline(h = 0, lty = 2, lwd = 2)

Although affinity() uses all of the variables in the transect, diagram() treats a transect like a system of one variable, so we get a plot of affinity vs. the first variable (temperature). We obtain three plots:

  1. The reactions are balanced on the first basis species with non-zero coefficients, CO2. This is the same as normalizing on C, because no other basis species has C.
  2. Balance on a different species, H2.
  3. No balancing constraint (balance = 1). This just shows the affinity of each reaction as given (that is, per mole of formed species), which is how the results were presented by Shock and Canovas (2010).

14. Calculate adjusted and non-standard Gibbs energy with subcrt()

In normal use, subcrt() returns standard molal properties (G, H, S, Cp, and V) for species in their standard states, referenced as unit activity or fugacity. Two deviations from this default setting are possible: non-standard properties for specified activity, and adjusted properties for activity coefficients.

First let's look at how adjusted properties depend on activity coefficients. This example uses a particular nonideality model based on @Alb03:

# Set the nonideality model
old_ni <- nonideal("Alberty")
# Calculate standard and adjusted Gibbs energy at 25 °C
species <- c("MgH2ATP", "MgHATP-", "MgATP-2")
subcrt(species, T = 25, IS = c(0, 0.25), property = "G")$out
# Restore the original nonideality model
nonideal(old_ni)

Notice how logarithms of the activity coefficients (loggam) are more negative for the higher-charged species. The activity coefficients have a stabilizing effect: adjusted Gibbs energies (at I > 0) are less than the standard Gibbs energies (at I = 0).

Now let's change the activities to get non-standard properties.

species <- c("Mg+2", "ATP-4", "MgATP-2")
coeffs <- c(-1, -1, 1)
T <- c(25, 50, 100)
# Drop some columns for more compact output
idrop <- c(2:4, 6:9)
# Unit activity: affinity is opposite of standard Gibbs energy
subcrt(species, coeffs, T = T, logact = c(0, 0, 0))$out[, -idrop]
# Non-unit activity: affinity is opposite of non-standard Gibbs energy
subcrt(species, coeffs, T = T, logact = c(-5, -4, -3))$out[, -idrop]

Note that:

The first call above specifies unit activities of all the species in the reaction.

The second call specifies non-unit activities of the species.

15. Use the output of subcrt() to format reactions for diagrams

16. Extract lines and image data from the output of diagram()

17. Counting elements and summing and writing formulas with makeup() and as.chemical.formula()

18. A brief introduction to buffers

19. A brief introduction to proteins and groupwise relative stabilities

Further Resources

Help Pages

Demos

Explore demos with demo(package = "CHNOSZ"). You can also use demos() to run all the demos or just one (e.g. demos("mosaic")).

More use cases for mosaic()

Solubility contours with solubility()

Other contour plots

Calculations using the output of diagram()

Activity buffers

Other thermodynamic models

Calculations with proteins

Other Vignettes

Additional vignettes cover:

Frequently asked questions

The FAQ is a non-comprehensive collection of questions and answers about CHNOSZ. Please use the Discussions forum on GitHub for new questions.

OBIGT thermodynamic database

The OBIGT vignette is generated from reference information in the database and lists all literature citations for species arranged by default and optional data files.

Customizing the thermodynamic database

The custom_data vignette describes add.OBIGT() for adding data from files, mod.OBIGT() for updating or adding parameters of particular species, and logK.to.OBIGT() for generating parameters from logK values.

Fitting thermodynamic data

The EOSregress vignette shows how to fit experimental data (volume and heat capacity) using constructed equation-of-state models.

Creating multi-metal diagrams (aside: why affinity?)

The idea for creating stability diagrams in CHNOSZ came from Harold Helgeson's geochemistry course. There, we constructed diagrams that were "balanced on" a metal. For instance, in a system balanced on Al, 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), then draws a line between stability fields where the non-standard Gibbs energy of reaction is zero. The grid-based method, used in CHNOSZ, looks at reactions to compose individual species from the basis species (let's call them formation reactions), 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 [depending on the author, e.g. @AL19;@STK19] "overall" or "non-standard" in front of the latter. "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 method have the same limitation: every candidate species must have non-zero stoichiometry of a given metal (or of a basis species with that metal).

The multi-metal vignette has some techniques for overcoming this limitation of balancing reactions on a single metal.

References

The CHNOSZ package incorporates data and methods from various sources. To view citation information for data sources:

# View all references in a browser
thermo.refs()

# Return data frame with references for one or more species
thermo.refs(info("CH4", c("aq", "gas")))

For citing CHNOSZ itself, see "How should CHNOSZ be cited?" in the FAQ.

Document history

wzxhzdk:41



Try the CHNOSZ package in your browser

Any scripts or data that you put into this service are public.

CHNOSZ documentation built on April 23, 2025, 3 p.m.

R Package Documentation

rdrr.io home R language documentation Run R code online

Browse R Packages

CRAN packages Bioconductor packages R-Forge packages GitHub packages

We want your feedback!

Note that we can't provide technical support on individual packages. You should contact the package authors for that.
GitHub issue tracker
ian@mutexlabs.com
Personal blog

 
Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.

Close