inst/doc/anintro.R

In CHNOSZ: Thermodynamic Calculations and Diagrams for Geochemistry

## ----setup, include=FALSE-----------------------------------------------------
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'))

## ----library_CHNOSZ-----------------------------------------------------------
library(CHNOSZ)

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

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

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

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

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

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

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

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

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

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

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

## ----mosaic, fig.cap = "Mosaic diagram showing effect of aqueous S speciation on the relative stabilities of Cu minerals and aqueous species"----
# 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)

## ----equilibrate, fig.cap="Carbonate speciation as a function of pH and temperature"----
# 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)

## ----corundum_solubility, fig.cap="Solubility of corundum (green line) and equilibrium concentrations of aqueous species (black lines)"----
# 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")

## ----corundum_solubility_IS, fig.cap="Solubility of corundum dependent on ionic strength"----
# 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")

## ----bquote-------------------------------------------------------------------

## ----dissolution_logK, fig.cap="Equilibrium constants of dissolution reactions"----
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"))

## ----retrieve_Al_default, results = "show"------------------------------------
# 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

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

## ----add_OBIGT_SLOP98, results = "show"---------------------------------------
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))

## ----corundum_solubility_2, fig.cap="Corundum solubility with species from SLOP98"----
# 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()

## ----basis_Al_F_OH_1, error=TRUE----------------------------------------------
try({
basis(c("Al+3", "F-", "OH-"))
})

## ----basis_Al_F_OH_2, error=TRUE----------------------------------------------
try({
basis(c("Al+3", "F-", "H+", "H2O"))
})

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

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

## ----Pourbaix_Mn, fig.cap = "Pourbaix diagram for Mn 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")

## ----NaCl, results = "show"---------------------------------------------------
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))

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

## ----solubility, eval=FALSE, cache=FALSE--------------------------------------
# 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)

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

## ----convert_ppm, fig.cap = "Solubility in units of ppm"----------------------
sppm <- convert(s, "ppm")
levels <- c(1e-6, 1e-3, 1e0, 1e3)
diagram(sppm, levels = levels)

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

## ----transect_results, results="show", cache=FALSE----------------------------
# Print pH and Eh values used for calculation
a$vals
# Print affinity values calculated for each species
a$values

## ----rainbow, echo=FALSE, fig.cap="Affinities of organic synthesis reactions per mole of C, H2, or formed species", fig.margin=FALSE, fig.fullwidth=TRUE, fig.width=12, fig.height=3, cache=TRUE----
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)

## ----rainbow, eval=FALSE, cache=FALSE-----------------------------------------
# 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)

## ----subcrt_nonideal, results = "show"----------------------------------------
# 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)

## ----subcrt_affinity, results = "show"----------------------------------------
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]

## ----refs, eval = FALSE-------------------------------------------------------
# # 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")))

## ----the_end------------------------------------------------------------------
   ######    ##   ##    ##   ##    ######     #####  #####
 ##         ##===##    ## \\##   ##    ##     \\       //
 ######    ##   ##    ##   ##    ######    #####      #####

## ----sessionInfo--------------------------------------------------------------
sessionInfo()