Package 'CHNOSZ'

Thermodynamic calculations and diagrams for geochemistry

Description

CHNOSZ is a package for thermodynamic calculations, primarily with applications in geochemistry and compositional biology. It can be used to calculate the standard molal thermodynamic properties and chemical affinities of reactions relevant to geobiochemical processes, and to visualize the equilibrium activities of species on chemical speciation and predominance diagrams.

Warm Tips

• To view the manual, run help.start() then select ‘Packages’ and ‘CHNOSZ’. Examples in the function help pages can be run by pasting the code block into the R console.

• Also check out the vignette anintro (An Introduction to CHNOSZ).

• Run the command examples() to run all of the examples provided in CHNOSZ. This should take about a minute.

Getting Help

Each help page (other than this one) has been given one of the following “concept index entries”:

• Main workflow: info, subcrt, basis, species, affinity, equilibrate, diagram

• Extended workflow: swap.basis, buffer, mosaic, EOSregress

• Thermodynamic data: data, extdata, add.OBIGT, util.data

• Thermodynamic calculations: util.formula, makeup, util.units, Berman, nonideal, util.misc

• Water properties: water, util.water, DEW, IAPWS95

• Protein properties: protein.info, add.protein, util.protein, util.seq, ionize.aa

• Other tools: examples, taxonomy

• Utility functions: util.expression, util.plot, util.array, util.list, palply

These concept entries are visible to help.search (aka ??). For example, help pages related to thermodynamic data can be listed using ??"thermodynamic data".

Warning

All thermodynamic data and examples are provided on an as-is basis. It is up to you to check not only the accuracy of the data, but also the suitability of the data AND computational techniques for your problem. By combining data taken from different sources, it is possible to build an inconsistent and/or nonsensical calculation. An attempt has been made to provide a default database (OBIGT) that is internally consistent, but no guarantee can be made. If there is any doubt about the accuracy or suitability of data for a particular problem, please consult the primary sources (see thermo.refs).

Acknowledgements

This package would not exist without the scientific influence and friendship of the late Professor Harold C. Helgeson. The ‘src/H2O92D.f’ file with Fortran code for calculating the thermodynamic and electrostatic properties of H₂O is modified from the SUPCRT92 package (Johnson et al., 1992).

Work on CHNOSZ at U.C. Berkeley from ca. 2003 to 2008 was supported in part by research grants to HCH from the U.S. National Science Foundation and Department of Energy. In 2009–2011, development of this package was partially supported by NSF grant EAR-0847616 to JMD.

References

Johnson, J. W., Oelkers, E. H. and Helgeson, H. C. (1992) SUPCRT92: A software package for calculating the standard molal thermodynamic properties of minerals, gases, aqueous species, and reactions from 1 to 5000 bar and 0 to 1000°C. Comp. Geosci. 18, 899–947. doi:10.1016/0098-3004(92)90029-Q

---

Functions to work with the thermodynamic database

Description

Add or modify species in the thermodynamic database.

Usage

add.OBIGT(file, species = NULL, force = TRUE)
  mod.OBIGT(..., zap = FALSE)

Arguments

file	character, path to a file
species	character, names of species to load from file
force	logical, force replacement of already existing species?
...	character or numeric, properties of species to modify in the thermodynamic database
zap	logical, clear preexisting parameters?

Details

Note: change made to OBIGT are lost if you reload the database by calling reset or OBIGT or if you quit the R session without saving it.

add.OBIGT is used to update the thermodynamic database (thermo$OBIGT) in the running session. The format (column names) of the specified file must be the same as the extdata/OBIGT/*.csv files provided with CHNOSZ. Note that this includes both the E_units and model columns, which were added in versions 1.3.3 and 2.0.0.

file is first matched against the names of files in the extdata/OBIGT directory packaged with CHNOSZ. In this case, the file suffixes are removed, so ‘⁠DEW⁠’, ‘⁠organic_aq⁠’, and ‘⁠organic_cr⁠’ are valid names. If there are no matches to a system file, then file is interpreted as the path a user-supplied file.

If species is NULL (default), all species listed in the file are used. If species is given and matches the name(s) of species in the file, only those species are added to the database.

By default, species in the file replace any existing species having the same combination of name and state. Set force to FALSE to avoid replacing species that are present in (thermo()$OBIGT).

When adding (not replacing) species, there is no attempt made to keep the order of physical states in the database (aq-cr-liq-gas); the function simply adds new rows to the end of thermo$OBIGT. As a result, retrieving the properties of an added aqueous species using info requires an explicit state="aq" argument to that function if a species with the same name is present in one of the cr, liq or gas states.

mod.OBIGT changes one or more of the properties of species or adds species to the thermodynamic database. The name of the species to add or change must be supplied as the first argument of ... or as a named argument (named ‘⁠name⁠’). Additional arguments to mod.OBIGT refer to the name of the property(s) to be updated and are matched to any part of compound column names in thermo()$OBIGT. For instance, either ‘⁠z⁠’ or ‘⁠T⁠’ matches the ‘⁠z.T⁠’ column. The values provided should also include order-of-magnitude scaling of HKF and DEW model parameters (see thermo).

When adding new species, a chemical formula should be included along with the values of any of the thermodynamic properties. The formula is taken from the ‘⁠formula⁠’ argument, or if that is missing, is taken to be the same as the ‘⁠name⁠’ of the species. An error occurs if the formula is not valid (i.e. can not be parsed by makeup). For new species, properties that are not specified become NA, except for ‘⁠state⁠’ and ‘⁠E_units⁠’, which take default values from thermo()$opt. These defaults can be overridden by giving a value for ‘⁠state⁠’ or ‘⁠E_units⁠’ in the arguments.

‘⁠model⁠’, if missing, is set to ‘⁠HKF⁠’ for state == "aq" or ‘⁠CGL⁠’ otherwise. When modifying some existing minerals in OBIGT, model = "CGL" should be explicitly given in order to override the Berman model.

When modifying species, the parameters indicated by the named arguments of mod.OBIGT are updated. Use zap = TRUE to replace all prexisting parameters (except for state and model) with NA values.

Value

The values returned (invisible-y) are the indices of the added and/or modified species.

References

Apps, J. and Spycher, N. (2004) Data qualification for thermodynamic data used to support THC calculations. DOC.20041118.0004 ANL-NBS-HS-000043 REV 00. Bechtel SAIC Company, LLC.

Bazarkina, E. F., Zotov, A. V., and Akinfiev, N. N. (2010) Pressure-dependent stability of cadmium chloride complexes: Potentiometric measurements at 1-1000 bar and 25°C. Geology of Ore Deposits 52, 167–178. doi:10.1134/S1075701510020054

Kitadai, N. (2014) Thermodynamic prediction of glycine polymerization as a function of temperature and pH consistent with experimentally obtained results. J. Mol. Evol. 78, 171–187. doi:10.1007/s00239-014-9616-1

Shock, E. L., Helgeson, H. C. and Sverjensky, D. A. (1989) Calculation of the thermodynamic and transport properties of aqueous species at high pressures and temperatures: Standard partial molal properties of inorganic neutral species. Geochim. Cosmochim. Acta 53, 2157–2183. doi:10.1016/0016-7037(89)90341-4

Stefánsson, A. (2001) Dissolution of primary minerals of basalt in natural waters. I. Calculation of mineral solubilities from 0°C to 350°C. Chem. Geol. 172, 225–250. doi:10.1016/S0009-2541(00)00263-1

Sverjensky, D. A., Shock, E. L., and Helgeson, H. C. (1997) Prediction of the thermodynamic properties of aqueous metal complexes to 1000 °C and 5 kbar. Geochim. Cosmochim. Acta 61, 1359–1412. doi:10.1016/S0016-7037(97)00009-4

See Also

thermo (description of OBIGT), mod.buffer (modify buffer definitions), logB.to.OBIGT (fit thermodynamic parameters to formation constants)

Examples

## Modify an existing species (not real properties)
ialanine <- mod.OBIGT("alanine", state = "cr", G = 0, H = 0, S = 0)
# We have made the values of G, H, and S inconsistent
# with the elemental composition of alanine, so the following 
# now produces a message about that
info(ialanine)

## Add an aqueous species (default) with Gibbs energy given in J/mol 
## (the same as the default) and today's date
date <- as.character(Sys.Date())
iCl2O <- mod.OBIGT("Cl2O", date = date, E_units = "J", G = 87738)
info(iCl2O)

## Add a solid species with a name that is different from the formula
mod.OBIGT("lorem-ipsum", formula = "C123", state = "cr", G = -12345678)
# Retrieve the data for this species using either name or formula
info(info("lorem-ipsum"))
info(info("C123"))
# Reset database
OBIGT()

## Using add.OBIGT():
# Compare stepwise stability constants of cadmium chloride complexes
# using data from Sverjensky et al., 1997 and Bazarkina et al., 2010
Cdspecies <- c("Cd+2", "CdCl+", "CdCl2", "CdCl3-", "CdCl4-2")
P <- c(1, seq(25, 1000, 25))
SSH97 <- lapply(1:4, function(i) {
  subcrt(c(Cdspecies[i], "Cl-", Cdspecies[i+1]),
    c(-1, -1, 1), T=25, P=P)$out$logK
})
file <- system.file("extdata/adds/BZA10.csv", package="CHNOSZ")
add.OBIGT(file)
BZA10 <- lapply(1:4, function(i) {
  subcrt(c(Cdspecies[i], "Cl-", Cdspecies[i+1]),
    c(-1, -1, 1), T=25, P=P)$out$logK
})
matplot(P, do.call(cbind, SSH97), type="l")
matplot(P, do.call(cbind, BZA10), type="l", add=TRUE, lwd=2)
legend("topleft", legend=c("", "", "Sverjensky et al., 1997",
  "Bazarkina et al., 2010"), lwd=c(0, 0, 1, 2), bty="n")
# Make reaction labels
y <- c(1.8, 0.2, -0.5, -1)
invisible(lapply(1:4, function(i) {
  text(800, y[i], describe.reaction(subcrt(c(Cdspecies[i], "Cl-",
    Cdspecies[i+1]), c(-1, -1, 1), T=25, P=1)$reaction))
}))
# Restore default database
OBIGT()

# Another use of add.OBIGT()
# Compare Delta G of AABB = UPBB + H2O
# (Figure 9 of Kitadai, 2014)
# Default database has values from Kitadai, 2014
Kit14 <- subcrt(c("[AABB]", "[UPBB]", "H2O"), c(-1, 1, 1), T = seq(0, 300, 10))
# Load superseded parameters for [UPBB] from Dick et al., 2006
mod.OBIGT("[UPBB]", G = -21436, H = -45220, S = 1.62)
DLH06 <- subcrt(c("[AABB]", "[UPBB]", "H2O"), c(-1, 1, 1), T = seq(0, 300, 10))
xlab <- axis.label("T"); ylab <- axis.label("DG", prefix="k")
plot(Kit14$out$T, Kit14$out$G/1000, type = "l", ylim = c(10, 35),
     xlab = xlab, ylab = ylab)
lines(DLH06$out$T, DLH06$out$G/1000, lty = 2)
legend("topleft", c("Dick et al., 2006", "Kitadai, 2014"), lty = c(2, 1))
title(main = "AABB = UPBB + H2O; after Figure 9 of Kitadai, 2014")
# Restore default database
OBIGT()

# Another use of add.OBIGT(): calculate Delta G of
# H4SiO4 = SiO2 + 2H2O using various sources of data for SiO2.
# First, get H4SiO4 from Stefansson, 2001
add.OBIGT("AS04", "H4SiO4")
T <- seq(0, 350, 10)
s1 <- subcrt(c("H4SiO4", "SiO2", "H2O"), c(-1, 1, 2), T = T)
# Now, get SiO2 from Apps and Spycher, 2004
add.OBIGT("AS04", "SiO2")
s2 <- subcrt(c("H4SiO4", "SiO2", "H2O"), c(-1, 1, 2), T = T)
# Plot logK from the first and second calculations
plot(T, s1$out$G, type = "l", xlab = axis.label("T"),
  ylab = axis.label("DG"), ylim = c(-500, 2500))
lines(T, s2$out$G, lty = 2)
# Add title and legend
title(main = describe.reaction(s1$reaction))
stxt <- lapply(c("H4SiO4", "SiO2", "SiO2"), expr.species)
legend("top", c("Shock et al., 1989", "Apps and Spycher, 2004"),
  title = as.expression(expr.species("SiO2")), lty = c(1, 2))
legend("topright", "Stef\u00e1nsson, 2001",
  title = as.expression(expr.species("H4SiO4")))
abline(h = 0, lty = 3, col = 8)
# Take-home message: SiO2 from Ste01 is compatible with H4SiO4 from Ste01
# at low T, but SiO2 from Shock et al., 1989 (the default in OBIGT) isn't
OBIGT()

---

Amino acid compositions of proteins

Description

Functions to get amino acid compositions and add them to protein list for use by other functions.

Usage

add.protein(aa, as.residue = FALSE)

Arguments

aa	data frame, amino acid composition in the format of thermo()$protein
as.residue	logical, normalize by protein length?

Details

A ‘⁠protein⁠’ in CHNOSZ is defined by its identifying information and the amino acid composition, stored in thermo$protein. The names of proteins in CHNOSZ are distinguished from those of other chemical species by having an underscore character ("_") that separates two identifiers, referred to as the protein and organism. An example is ‘⁠LYSC_CHICK⁠’. The purpose of the functions described here is to identify proteins and work with their amino acid compositions. From the amino acid compositions, the thermodynamic properties of the proteins can be estimated by group additivity.

Given a data frame of amino acid compositions in the format of thermo()$protein, add.protein adds them to thermo()$protein for use by other functions in CHNOSZ. The amino acid compositions of proteins in aa with the same name as one in thermo()$protein are replaced. Set as.residue to TRUE to normalize by protein length; each input amino acid composition is divided by the corresponding number of residues, with the result that the sum of amino acid frequencies for each protein is 1.

Value

For add.protein, the rownumbers of thermo()$protein that are added and/or replaced.

See Also

read_fasta for reading amino acid compositions from FASTA files.

pinfo for protein-level functions (length, chemical formulas, reaction coefficients of basis species).

Examples

# Read a file with the amino acid compositions of several poliovirus protein subunits
file <- system.file("extdata/protein/POLG.csv", package = "CHNOSZ")
aa <- read.csv(file)

# Add the proteins to CHNOSZ
iprotein <- add.protein(aa)
# Calculate length and elemental formula
protein.length(iprotein)
protein.formula(iprotein)

---

Chemical affinities of formation reactions

Description

Calculate the chemical affinities of formation reactions of species.

Usage

affinity(..., property = NULL, sout = NULL, exceed.Ttr = FALSE,
    exceed.rhomin = FALSE, return.buffer = FALSE, return.sout = FALSE,
    balance = "PBB", iprotein = NULL, loga.protein = 0, transect = NULL)

Arguments

...	numeric, zero or more named arguments, used to identify the variables of interest in the calculations. For argument recall, pass the output from a previous calculation of affinity as an unnamed first argument.
property	character, the property to be calculated. Default is ‘⁠A⁠’, for chemical affinity of formation reactions of species of interest
sout	list, output from subcrt
exceed.Ttr	logical, allow subcrt to compute properties for phases beyond their transition temperature?
exceed.rhomin	logical, allow subcrt to compute properties of species in the HKF model below 0.35 g cm-3?
return.buffer	logical. If TRUE, and a buffer has been associated with one or more basis species in the system, return the values of the activities of the basis species calculated using the buffer. Default is FALSE.
return.sout	logical, return only the values calculated with subcrt?
balance	character. This argument is used to identify a conserved basis species (or ‘⁠PBB⁠’) in a chemical activity buffer. Default is ‘⁠PBB⁠’.
iprotein	numeric, indices of proteins in thermo$protein for which to calculate properties
loga.protein	numeric, logarithms of activities of proteins identified in iprotein
transect	logical, force a transect calculation, even for three or fewer values of the variables?

Details

affinity calculates the chemical affinities of reactions to form the species of interest from the basis species. The equation used to calculate chemical affinity (A), written for base-10 (decimal) logarithms, is $A/(2.303RT)$=$\log (K/Q)$, where $K$ is the equilibrium constant of the reaction, $Q$ is the activity product of the species in the reaction, and 2.303 is the conversion factor from natural to decimal logarithms. The calculation of chemical affinities relies on the current definitions of the basis species and species of interest. Calculations are possible at single values of temperature, pressure, ionic strength and chemical activities of the basis species, or as a function of one or more of these variables.

The argument property can be changed to calculate other thermodynamic properties of formation reactions. Valid properties are ‘⁠A⁠’ or NULL for chemical affinity, ‘⁠logK⁠’ or ‘⁠logQ⁠’ for logarithm of equilibrium constant and reaction activity product, or any of the properties available in subcrt except for ‘⁠rho⁠’. The properties returned are those of the formation reactions of the species of interest from the basis species. It is also possible to calculate the properties of the species of interest themselves (not their formation reactions) by setting the property to ‘⁠G.species⁠’, ‘⁠Cp.species⁠’, etc. Except for ‘⁠A⁠’, the properties of proteins or their reactions calculated in this manner are restricted to nonionized proteins.

Zero, one, or more leading arguments to the function identify which of the chemical activities of basis species, temperature, pressure and/or ionic strength to vary. The names of each of these arguments may be the formula of any of the basis species of the system, or ‘⁠T⁠’, ‘⁠P⁠’, ‘⁠pe⁠’, ‘⁠pH⁠’, ‘⁠Eh⁠’, or ‘⁠IS⁠’ (but names may not be repeated). The names of charged basis species such as ‘⁠K+⁠’ and ‘⁠SO4-2⁠’ should be quoted when used as arguments. The value of each argument is of the form c(min, max) or c(min, max, res) where min and max refer to the minimimum and maximum values of variable identified by the name of the argument, and res is the resolution, or number of points along which to do the calculations; res is assigned a default value of 256 if it is missing. For any arguments that refer to basis species, the numerical values are the logarithms of activity (or fugacity for gases) of that basis species.

If ‘⁠T⁠’, ‘⁠P⁠’, and/or ‘⁠IS⁠’ are not among the vars, their constant values can be supplied in T, P, or IS (in mol kg$^{-1}$). The units of ‘⁠T⁠’ and ‘⁠P⁠’ are those set by T.units and P.units (on program start-up these are °C and bar, respectively). sout, if provided, replaces the call to subcrt, which can greatly speed up the calculations if this intermediate result is stored by other functions. exceed.Ttr is passed to subcrt so that the properties of mineral phases beyond their transition temperatures can optionally be calculated.

If one or more buffers are assigned to the definition of basis species, the logarithms of activities of these basis species are taken from the buffer (see buffer).

The iprotein and loga.protein arguments can be used to compute the chemical affinities of formation reactions of proteins that are not in the current species definition. iprotein contains the indices (rownumbers) of desired proteins in thermo$protein. This uses some optimizations to calculate the properties of many proteins in a fraction of the time it would take to calculate them individually.

When the length(s) of the variables is(are) greater than 3, the function enters the ‘⁠transect⁠’ mode of operation. In this mode of operation, instead of performing the calculations on an $n$-dimensional grid, the affinities are calculated on a transect of changing T, P, and/or chemical activity of basis species.

Argument recall is invoked by passing a previous result of affinity as the first argument. The function then calls itself using the settings from the previous calculation, with additions or modifications indicated by the remaining arguments in the current function call.

Value

A list, elements of which are fun the name of the function (‘⁠affinity⁠’), args all of the arguments except for ‘⁠sout⁠’ (these are used for argument recall), sout output from subcrt, property name of the calculated property (‘⁠A⁠’ for chemical affinity), basis and species definition of basis species and species of interest in effect at runtime, T and P temperature and pressure, in the system units of Kelvin and bar, set to numeric() (length=0) if either one is a variable, vars the names of the variables, vals the values of the variables (a list, one element for each variable), values the result of the calculation (a list, one element for each species, with names taken from the species index in thermo$OBIGT). The elements of the lists in vals and values are arrays of $n$ dimensions, where $n$ is the number of variables. The values of chemical affinity of formation reactions of the species are returned in dimensionless units (for use with decimal logarithms, i.e., A/$2.303RT$).

Names other than ‘⁠T⁠’ or ‘⁠P⁠’ in vars generally refer to basis species, and the corresponding vals are the logarithms of activity or fugacity. However, if one or more of pe, Eh or pH is among the variables of interest, vals holds the values of the those variables as indicated.

References

Helgeson, H. C., Richard, L, McKenzie, W. F., Norton, D. L. and Schmitt, A. (2009) A chemical and thermodynamic model of oil generation in hydrocarbon source rocks. Geochim. Cosmochim. Acta 73, 594–695. doi:10.1016/j.gca.2008.03.004

See Also

ionize.aa, activated if proteins are among the species of interest, ‘⁠H+⁠’ is in the basis and thermo()$opt$ionize.aa is TRUE. equilibrate for using the results of affinity to calculate equilibrium activities of species, and diagram to plot the results. demo("saturation") for an example using the argument recall feature.

Examples

## Set up a system and calculate
## chemical affinities of formation reactions
basis(c("SiO2", "MgO", "H2O", "O2"), c(-5, -5, 0, 999))
species(c("quartz","enstatite","forsterite"))
# Chemical affinities (A/2.303RT) at 25 deg C and 1 bar
affinity()
# At higher temperature and pressure
affinity(T = 500, P = 2000)
# At 25 temperatures and pressures,
# some are in the low-density region so we suppress warnings
suppressWarnings(affinity(T = c(500, 1000, 5), P = c(1000, 5000, 5)))
# Equilibrium constants of formation reactions
affinity(property = "logK")
# Standard molal Gibbs energies of species,
# in units set by E.units() (default: J/mol)
affinity(property = "G.species")
# Standard molal Gibbs energies of reactions
affinity(property = "G")
# A T,P-transect
# (fluid pressure from Helgeson et al., 2009 Fig. 7)
affinity(T = c(25, 110, 115, 215), P = c(11, 335, 500, 1450))

---

Define basis species

Description

Define the basis species of a chemical system.

Usage

basis(species = NULL, state = NULL, logact = NULL,
    delete = FALSE, add = FALSE)

Arguments

species	character, names or formulas of species, or numeric, indices of species
state	character, physical states or names of buffers
logact	numeric, logarithms of activities or fugacities
delete	logical, delete the current basis definition?
add	logical, add species to the current basis definition?

Details

The basis species represent the possible range of chemical compositions for all the species of interest. As used here, a set of basis species is valid only if it satisifes two conditions: 1) the number of basis species is the same as the number of chemical elements (including charge) in those species and 2) the square matrix representing the elemental stoichiometries of the basis species has a real inverse.

To create a basis definition, call basis with the names or formulas of the basis species in the species argument, or all numeric values as species indices (rownumbers in thermo()$OBIGT). The special names ‘⁠pH⁠’, ‘⁠pe⁠’ and ‘⁠Eh⁠’ can also be used; they get translated into the names of the proton (‘⁠H+⁠’) and electron (‘⁠e-⁠’) as appropriate. If desired, include the state for the named species and the logarithms of activity (fugacity for gases) in logact. The latter defaults to zero (unit activity) if not specified.

To modify an existing basis definition, the physical states or logarithms of activities of species can be changed by calling basis with a species argument that has the formulas (not names) or indices of species in the existing basis. If either of the second or third arguments to basis is of type character, it refers to the physical state (if present in thermo()$OBIGT$state) or a chemical activity buffer (if present in thermo()$buffers$name). If either of these arguments is numeric it specifies the logarithms of activities (or fugacities for gases) of the basis species. In case ‘⁠pH⁠’, ‘⁠pe⁠’ or ‘⁠Eh⁠’ is named, the logarithm of activity of the basis species is converted from these values. For example, a value of 7 for pH is stored as a logarithm of activity of -7.

If add is TRUE, then the function attempts to add the indicated species to the basis definition. This only works if the enlarged set of species is a valid basis set as described above. If the formed species are currently defined, their formation reactions are modified accordingly (with zeroes for the newly added basis species).

If add is FALSE, and if basis is called with NULL values of both state and logact, the new set of species, if they are a valid basis set, completely replaces any existing basis definition. When this occurs, any existing species definition (created by the species function) is deleted. Call basis with delete set to TRUE or species set to ‘⁠""⁠’ to clear the basis definition and that of the species, if present.

If the value of basis is one of the keywords in the following table, the corresponding set of basis species is loaded, and their activities are given preset values. The basis species identified by these keywords are aqueous except for H₂O (liq), O₂ (gas) and Fe₂O₃ (hematite).

CHNOS	CO2, H2O, NH3, H2S, O2
CHNOS+	CO2, H2O, NH3, H2S, O2, H+
CHNOSe	CO2, H2O, NH3, H2S, e-, H+
CHNOPS+	CO2, H2O, NH3, H3PO4, H2S, O2, H+
CHNOPSe	CO2, H2O, NH3, H3PO4, H2S, e-, H+
MgCHNOPS+	Mg+2, CO2, H2O, NH3, H3PO4, H2S, O2, H+
MgCHNOPSe	Mg+2, CO2, H2O, NH3, H3PO4, H2S, e-, H+
FeCHNOS	Fe2O3, CO2, H2O, NH3, H2S, O2
FeCHNOS+	Fe2O3, CO2, H2O, NH3, H2S, O2, H+
QEC4	cysteine, glutamic acid, glutamine, H2O, O2
QEC	cysteine, glutamic acid, glutamine, H2O, O2
QEC+	cysteine, glutamic acid, glutamine, H2O, O2, H+
QCa	glutamine, cysteine, acetic acid, H2O, O2
QCa+	glutamine, cysteine, acetic acid, H2O, O2, H+

The logarithms of activities of amino acids in the ‘⁠QEC4⁠’ basis are -4 (i.e., basis II in Dick, 2016); those in ‘⁠QEC⁠’ and ‘⁠QEC+⁠’ are set to approximate concentrations in human plasma (see Dick, 2017).

Value

Returns the value of thermo()$basis after any modifications; or, if delete is TRUE, its value before deletion (invisibly).

References

Dick, J. M. (2016) Proteomic indicators of oxidation and hydration state in colorectal cancer. PeerJ 4:e2238. doi:10.7717/peerj.2238

Dick, J. M. (2017) Chemical composition and the potential for proteomic transformation in cancer, hypoxia, and hyperosmotic stress. PeerJ 5:e3421 doi:10.7717/peerj.3421

See Also

info to query the thermodynamic database in order to find what species are available. makeup is used by basis to generate the stoichiometric matrix from chemical formulas. swap.basis is used to change the chemical compounds (species formulas) used in the basis definition while keeping the chemical potentials of the elements unaltered. species for setting up the formation reactions from basis species.

Examples

## Define basis species
# with one, two or three elements
basis("O2")
basis(c("H2O", "O2"))
basis(c("H2O", "O2", "H+"))
## Clear the basis species
basis("")

## Not run: 
## Marked dontrun because they produce errors
# Fewer species than elements
basis(c("H2O", "H+"))
# More species than elements
basis(c("H2O", "O2", "H2", "H+"))
# Non-independent species
basis(c("CO2", "H2O", "HCl", "Cl-", "H+"))
## End(Not run)

## Specify activities and states
basis(c("H2O", "O2", "CO2"), c(-2, -78, -3), c("liq", "aq", "aq"))
# Change logarithms of activities/fugacities	
basis(c("H2O", "O2"), c(0, -80))	
# Change state of CO2
basis("CO2", "gas")

---

Thermodynamic properties of minerals

Description

Calculate thermodynamic properties of minerals using the equations of Berman (1988).

Usage

Berman(name, T = 298.15, P = 1, check.G = FALSE,
         calc.transition = TRUE, calc.disorder = TRUE)

Arguments

name	character, name of mineral
T	numeric, temperature(s) at which to calculate properties (K)
P	numeric, pressure(s) at which to calculate properties (bar)
check.G	logical, check consistency of G, H, and S?
calc.transition	logical, include calculation of polymorphic transition properties?
calc.disorder	logical, include calculation of disordering properties?

Details

This function calculates the thermodynamic properties of minerals at high P and T using equations given by Berman (1988). These minerals should be listed in thermo()$OBIGT with the state ‘⁠cr⁠’ and chemical formula, and optionally an abbreviation and references, but all other properties set to NA.

The standard state thermodynamic properties and parameters for the calculations are stored in data files under extdata/Berman, or can be read from a user-created file specified by thermo()$opt$Berman.

The equation used for heat capacity is C_P = k0 + k1*T^-0.5 + k2*T^-2 + k3*T^-3 + k4*T^-1 + k5*T + k6*T². This is an extended form Eq. 4 of Berman (1988) as used in the winTWQ program (Berman, 2007). The equation used for volume is V(P, T) / V(1 bar, 298.15 K) = 1 + v1 * (T - 298.15) + v2 * (T - 298.15)² + v3 * (P - 1) + v4 * (P - 1)² (Berman, 1988, Eq. 5, with terms reordered to follow winTWQ format). The equations used for lambda transitions follow Eqs. 8-14 of Berman (1988). The equation used for the disorder contribution between Tmin and Tmax is C_P[dis] = d0 + d1*T^-0.5 + d2*T^-2 + d3*T + d4*T² (Berman, 1988, Eq. 15). The parameters correspond to Tables 2 (GfPrTr, HfPrTr, SPrTr, VPrTr), 3a (k0 to k3), 4 (v1 to v4), 3b (transition parameters: Tlambda to dTH), and 5 (disorder parameters: Tmax, Tmin, d1 to d4 and Vad) of Berman (1988). Following the winTWQ data format, multipliers are applied to the volume parameters only (see below). Note that VPrTr is tabulated in J bar^-1 mol^-1, which is equal to 10 cm³ mol^-1.

A value for GfPrTr is not required and is only used for optional checks (see below). Numeric values (possibly 0) should be assigned for all of HfPrTr, SPrTr, VPrTr, k0 to k6 and v1 to v4. Missing (or NA) values are permitted for the transition and disorder parameters, for minerals where they are not used. The data files have the following 30 columns:

name	mineral name (must match an entry with a formula but NA properties in thermo()$OBIGT)
GfPrTr	standard Gibbs energy at 298.15 K and 1 bar (J mol-1) (Benson-Helgeson convention)
HfPrTr	standard enthalpy at 298.15 K and 1 bar (J mol-1)
SPrTr	standard entropy at 298.15 K and 1 bar (J mol-1 K-1)
VPrTr	standard volume at 298.15 K and 1 bar (J bar-1) [1 J bar-1 = 10 cm3]
k0 ... k6	k0 (J mol-1 K-1) to k6
v1	v1 (K-1) * 105
v2	v2 (K-2) * 105
v3	v3 (bar-1) * 105
v4	v4 (bar-2) * 108
Tlambda	Tλ (K)
Tref	Tref (K)
dTdP	dT / dP (K bar-1)
l1	l1 ((J/mol)0.5 K-1)
l2	l2 ((J/mol)0.5 K-2)
DtH	ΔTtH (J mol-1)
Tmax	temperature at which phase is fully disordered (TD in Berman, 1988) (K)
Tmin	reference temperature for onset of disordering (t in Berman, 1988) (K)
d0 ... d4	d0 (J mol-1 K-1) to d4
Vad	constant that scales the disordering enthalpy to volume of disorder (d5 in Berman, 1988)

The function outputs apparent Gibbs energies according to the Benson-Helgeson convention (ΔG = ΔH - TΔS) using the entropies of the elements in the chemical formula of the mineral to calculate ΔS (cf. Anderson, 2005). If check.G is TRUE, the tabulated value of GfTrPr (Benson-Helgeson) is compared with that calculated from HfPrTr - 298.15*DSPrTr (DSPrTr is the difference between the entropies of the elements in the formula and SPrTr in the table). A warning is produced if the absolute value of the difference between tabulated and calculated GfTrPr is greater than 1000 J/mol.

If the function is called with missing name, the parameters for all available minerals are returned.

Value

A data frame with T (K), P (bar), G, H, S, and Cp (energetic units in Joules), and V (cm³ mol^-1).

References

Anderson, G. M. (2005) Thermodynamics of Natural Systems, 2nd ed., Cambridge University Press, 648 p. https://www.worldcat.org/oclc/474880901

Berman, R. G. (1988) Internally-consistent thermodynamic data for minerals in the system Na₂O-K₂O-CaO-MgO-FeO-Fe₂O₃-Al₂O₃-SiO₂-TiO₂-H₂O-CO₂. J. Petrol. 29, 445-522. doi:10.1093/petrology/29.2.445

Berman, R. G. and Aranovich, L. Ya. (1996) Optimized standard state and solution properties of minerals. I. Model calibration for olivine, orthopyroxene, cordierite, garnet, and ilmenite in the system FeO-MgO-CaO-Al₂O₃-TiO₂-SiO₂. Contrib. Mineral. Petrol. 126, 1-24. doi:10.1007/s004100050233

Berman, R. G. (2007) winTWQ (version 2.3): A software package for performing internally-consistent thermobarometric calculations. Open File 5462, Geological Survey of Canada, 41 p. doi:10.4095/223425

Helgeson, H. C., Delany, J. M., Nesbitt, H. W. and Bird, D. K. (1978) Summary and critique of the thermodynamic properties of rock-forming minerals. Am. J. Sci. 278-A, 1–229. https://www.worldcat.org/oclc/13594862

Examples

# Other than the formula, the parameters aren't stored in
# thermo()$OBIGT, so this shows NAs
info(info("quartz", "cr"))
# Properties of alpha-quartz (aQz) at 298.15 K and 1 bar
Berman("quartz")
# Gibbs energies of aQz and coesite at higher T and P
T <- seq(200, 1300, 100)
P <- seq(22870, 31900, length.out = length(T))
G_aQz <- Berman("quartz", T = T, P = P)$G
G_Cs <- Berman("coesite", T = T, P = P)$G
# That is close to the univariant curve (Ber88 Fig. 4),
# so the difference in G is close to 0
DGrxn <- G_Cs - G_aQz
all(abs(DGrxn) < 100)  # TRUE

# Make a P-T diagram for SiO2 minerals (Ber88 Fig. 4)
basis(c("SiO2", "O2"), c("cr", "gas"))
species(c("quartz", "quartz,beta", "coesite"), "cr")
a <- affinity(T = c(200, 1700, 200), P = c(0, 50000, 200))
diagram(a)

## Getting data from a user-supplied file
## Ol-Opx exchange equilibrium, after Berman and Aranovich, 1996
species <- c("fayalite", "enstatite", "ferrosilite", "forsterite")
coeffs <- c(-1, -2, 2, 1)
T <- seq(600, 1500, 50)
Gex_Ber88 <- subcrt(species, coeffs, T = T, P = 1)$out$G
# Add data from BA96
datadir <- system.file("extdata/Berman/testing", package = "CHNOSZ")
add.OBIGT(file.path(datadir, "BA96_OBIGT.csv"))
thermo("opt$Berman" = file.path(datadir, "BA96_Berman.csv"))
Gex_BA96 <- subcrt(species, coeffs, T = seq(600, 1500, 50), P = 1)$out$G
# Ber88 is lower than BA96 at low T
(Gex_BA96 - Gex_Ber88)[1] > 0  # TRUE
# The curves cross at about 725 deg C (BA96 Fig. 8)
# (actually, in our calculation they cross closer to 800 deg C)
T[which.min(abs(Gex_BA96 - Gex_Ber88))]  # 800
# Reset the database (thermo()$OBIGT, and thermo()$opt$Berman)
reset()

---

Calculating buffered chemical activities

Description

Calculate values of activity or fugacity of basis species buffered by an assemblage of one or more species.

Usage

mod.buffer(name, species = NULL, state = "cr", logact = 0)

Arguments

name	character, name of buffer to add to or find in thermo()$buffer.
species	character, names or formulas of species in a buffer.
state	character, physical states of species in buffer.
logact	numeric, logarithms of activities of species in buffer.

Details

A buffer is treated here as assemblage of one or more species whose presence constrains values of the chemical activity (or fugacity) of one or more basis species. To perform calculations for buffers use basis to associate the name of the buffer with one or more basis species. After this, calls to affinity will invoke the required calculations. The calculated values of the buffered activites can be retrieved by setting return.buffer to TRUE (in affinity). The maximum number of buffered chemical activities possible for any buffer is equal to the number of species in the buffer; however, the user may then elect to work with the values for only one or some of the basis species calculated with the buffer.

The identification of a conserved basis species (or other reaction balancing rule) is required in calculations for buffers of more than one species. For example, in the pyrite-pyrrhotite-magnetite buffer ($\mathrm{FeS_2}$-$\mathrm{FeS}$-$\mathrm{Fe_3O_4}$) a basis species common to each species is one representing $Fe$. Therefore, when writing reactions between the species in this buffer $Fe$ is conserved while $\mathrm{H_2S}$ and $\mathrm{O2}$ are the variables of interest. The calculation for buffers attempts to determine which of the available basis species qualifies as a conserved quantity. This can be overriden with balance. The default value of balance is ‘⁠PBB⁠’, which instructs the function to use the protein backbone group as the conserved quantity in buffers consisting of proteins, but has no overriding effect on the computations for buffers without proteins.

To view the available buffers, print the thermo()$buffer object. Buffer definitions can be added to this dataframe with mod.buffer. The defaults for state and logact are intended for mineral buffers. If name identifies an already defined buffer, this function modifies the logarithms of activities or states of species in that buffer, optionally restricted to only those species given in species.

It is possible to assign different buffers to different basis species, in which case the order of their calculation depends on their order in thermo()$buffers. This function is compatible with systems of proteins, but note that for buffers made of proteins the buffer calculations presently use whole protein formulas (instead of residue equivalents) and consider nonionized proteins only.

References

Garrels, R. M. (1960) Mineral Equilibria. Harper & Brothers, New York, 254 p. https://www.worldcat.org/oclc/552690

See Also

diagram with type set to the name of a basis species solves for the activity of the basis species.

Examples

## List the buffers
thermo()$buffer
# Another way to do it, for a specific buffer
print(mod.buffer("PPM"))

## Buffer made of one species
# Calculate the activity of CO2 in equilibrium with
# (a buffer made of) acetic acid at a given activity
basis("CHNOS")
basis("CO2", "AC")
# What activity of acetic acid are we using?
print(mod.buffer("AC"))
# Return the activity of CO2
affinity(return.buffer = TRUE)$CO2  # -7.057521
# As a function of oxygen fugacity
affinity(O2 = c(-85, -70, 4), return.buffer = TRUE)
# As a function of logfO2 and temperature
affinity(O2 = c(-85, -70, 4), T = c(25, 100, 4), return.buffer = TRUE)
# Change the activity of species in the buffer
mod.buffer("AC", logact = -10)
affinity(O2 = c(-85,-70,4), T = c(25, 100, 4), return.buffer = TRUE)

## Buffer made of three species
## Pyrite-Pyrrhotite-Magnetite (PPM)
# Specify basis species and initial activities
basis(c("FeS2", "H2S", "O2", "H2O"), c(0, -10, -50, 0))
# Note that the affinity of formation of pyrite,
# which corresponds to FeS2 in the basis, is zero
species(c("pyrite", "pyrrhotite", "magnetite"))
affinity(T = c(200, 400, 11), P = 2000)$values
# Setup H2S and O2 to be buffered by PPM
basis(c("H2S", "O2"), c("PPM", "PPM"))
# Inspect values of H2S activity and O2 fugacity
affinity(T = c(200, 400, 11), P = 2000, return.buffer = TRUE, exceed.Ttr = TRUE)
# Calculate affinities of formation reactions of species in the buffer
a <- affinity(T = c(200, 400, 11), P = 2000, exceed.Ttr = TRUE)$values
# The affinities for species in the buffer are all equal to zero
all.equal(as.numeric(a[[1]]), rep(0, 11))  # TRUE
all.equal(as.numeric(a[[2]]), rep(0, 11))  # TRUE
all.equal(as.numeric(a[[3]]), rep(0, 11))  # TRUE

## Buffer made of one species: show values of logfO2 on an 
## Eh-pH diagram; after Garrels, 1960, Figure 6
basis("CHNOSe")
# Here we will buffer the activity of the electron by O2
mod.buffer("O2", "O2", "gas", 999)
basis("e-", "O2")
# Start our plot, then loop over values of logfO2
thermo.plot.new(xlim = c(0, 14), ylim = c(-0.8, 1.2),
  xlab = "pH",ylab = axis.label("Eh"))
# The upper and lower lines correspond to the upper
# and lower stability limits of water
logfO2 <- c(0, -20, -40, -60, -83.1)
for(i in 1:5) {
  # Update the logarithm of fugacity (logact) of O2 in the buffer
  mod.buffer("O2", "O2", "gas", logfO2[i])
  # Get the values of the logarithm of activity of the electron
  a <- affinity(pH = c(0, 14, 15), return.buffer = TRUE)
  # Convert values of pe (-logact of the electron) to Eh
  Eh <- convert(-as.numeric(a$`e-`), "Eh")
  lines(seq(0, 14, length.out = 15), Eh)
  # Add some labels
  text(seq(0, 14, length.out = 15)[i*2+2], Eh[i*2+2],
    paste("logfO2 =", logfO2[i]))
}
title(main = paste("Relation between logfO2(g), Eh and pH at\n",
  "25 degC and 1 bar. After Garrels, 1960"))

## Buffer made of two species
# Conditions for metastable equilibrium among 
# CO2 and acetic acid. note their starting activities:
print(mod.buffer("CO2-AC")) 
basis("CHNOS")
basis("O2", "CO2-AC")
affinity(return.buffer = TRUE)  # logfO2 = -75.94248
basis("CO2", 123)  # what the buffer reactions are balanced on
affinity(return.buffer = TRUE)  # unchanged
# Consider more oxidizing conditions
mod.buffer("CO2-AC", logact = c(0, -10))
affinity(return.buffer = TRUE)

---

Deep Earth Water (DEW) model

Description

Calculate thermodynamic properties of water using the Deep Earth Water (DEW) model.

Usage

calculateDensity(pressure, temperature, error = 0.01)
  calculateGibbsOfWater(pressure, temperature)
  calculateEpsilon(density, temperature)
  calculateQ(density, temperature)

Arguments

pressure	numeric, pressure (bar)
temperature	numeric, temperature (°C)
error	numeric, residual error for bisection calculation
density	numeric, density (g/cm^3)

Details

The Deep Earth Water (DEW) model, described by Sverjensky et al., 2014, extends the applicability of the revised HKF equations of state to 60 kbar. This implementation of DEW is based on the VBA macro code in the May, 2017 version of the DEW spreadsheet downloaded from http://www.dewcommunity.org/. The spreadsheet provides multiple options for some calculations; here the default equations for density of water (Zhang and Duan, 2005), dielectric constant (Sverjensky et al., 2014) and Gibbs energy of water (integral of volume, equation created by Brandon Harrison) are used.

Comments in the original code indicate that calculateGibbsOfWater is valid for 100 ≤ T ≤ 1000 °C and P ≥ 1000 bar. Likewise, the power function fit of the dielectric constant (epsilon) is valid for 100 ≤ T ≤ 1200 °C and P ≥ 1000 bar (Sverjensky et al., 2014).

Value

The calculated values of density, Gibbs energy, and the Q Born coefficient have units of g/cm^3, cal/mol, and bar^-1 (epsilon is dimensionless).

References

Sverjensky, D. A., Harrison, B. and Azzolini, D. (2014) Water in the deep Earth: The dielectric constant and the solubilities of quartz and corundum to 60 kb and 1,200 °C. Geochim. Cosmochim. Acta 129, 125–145. doi:10.1016/j.gca.2013.12.019

Zhang, Z. and Duan, Z. (2005) Prediction of the PVT properties of water over wide range of temperatures and pressures from molecular dynamics simulation. Phys. Earth Planet. Inter. 149, 335–354. doi:10.1016/j.pepi.2004.11.003

See Also

water.DEW; use water("DEW") to activate these equations for the main functions in CHNOSZ.

Examples

pressure <- c(1000, 60000)
temperature <- c(100, 1000)
calculateGibbsOfWater(pressure, temperature)
(density <- calculateDensity(pressure, temperature))
calculateEpsilon(density, temperature)
calculateQ(density, temperature)

---

Chemical activity diagrams

Description

Plot equilibrium chemical activity (1-D speciation) or equal-activity (2-D predominance) diagrams as a function of chemical activities of basis species, temperature and/or pressure.

Usage

diagram(
    # species affinities or activities
    eout,
    # type of plot
    type = "auto", alpha = FALSE, normalize = FALSE,
    as.residue = FALSE, balance = NULL, groups = as.list(1:length(eout$values)),
    # figure size and sides for axis tick marks
    xrange = NULL, mar = NULL, yline = par("mgp")[1]+0.3, side = 1:4,
    # axis limits and labels
    ylog = TRUE, xlim = NULL, ylim = NULL, xlab = NULL, ylab = NULL,
    # character sizes
    cex = par("cex"), cex.names = 1, cex.axis = par("cex"),
    # line styles
    lty = NULL, lty.cr = NULL, lty.aq = NULL, lwd = par("lwd"), dotted = NULL,
    spline.method = NULL, contour.method = "edge", levels = NULL,
    # colors
    col = par("col"), col.names = par("col"), fill = NULL,
    fill.NA = "gray80", limit.water = NULL,
    # field and line labels
    names = NULL, format.names = TRUE, bold = FALSE, italic = FALSE, 
    font = par("font"), family = par("family"), adj = 0.5,
    dx = 0, dy = 0, srt = 0, min.area = 0,
    # title and legend
    main = NULL, legend.x = NA,
    # plotting controls
    add = FALSE, plot.it = TRUE, tplot = TRUE, ...)
  find.tp(x)

Arguments

eout	list, object returned by affinity, equilibrate or related functions
type	character, type of plot, or name of basis species whose activity to plot
alpha	logical or character (‘⁠balance⁠’), for speciation diagrams, plot degree of formation instead of activities?
normalize	logical, divide chemical affinities by balance coefficients and rescale activities to whole formulas?
as.residue	logical, divide chemical affinities by balance coefficients (no rescaling)?
balance	character, balancing constraint; see equilibrate
groups	list of numeric, groups of species to consider as a single effective species
xrange	numeric, range of x-values between which predominance field boundaries are plotted
mar	numeric, margins of plot frame
yline	numeric, margin line on which to plot the y-axis name
side	numeric, which sides of plot to draw axes
xlim	numeric, limits of x-axis
ylim	numeric, limits of y-axis
xlab	character, label to use for x-axis
ylab	character, label to use for y-axis
ylog	logical, use a logarithmic y-axis (on 1D degree diagrams)?
cex	numeric, character expansion (scaling relative to current)
cex.names	numeric, character expansion factor to be used for names of species on plots
cex.axis	numeric, character expansion factor for names of axes
lty	numeric, line types to be used in plots
lty.cr	numeric, line types for cr-cr boundaries (between two minerals)
lty.aq	numeric, line types for aq-aq boundaries (between two aqueous species)
lwd	numeric, line width
dotted	numeric, how often to skip plotting points on predominance field boundaries (to gain the effect of dotted or dashed boundary lines)
spline.method	character, method used in splinefun
contour.method	character, labelling method used in contour (use NULL for no labels).
levels	numeric, levels at which to draw contour lines
col	character, color of activity lines (1D diagram) or predominance field boundaries (2D diagram)
col.names	character, colors for labels of species
fill	character, colors used to fill predominance fields
fill.NA	character, color for grid points with NA values
limit.water	NULL or logical, set NA values beyond water stability limits?
names	character, names of species for activity lines or predominance fields
format.names	logical, apply formatting to chemical formulas?
bold	logical, use bold formatting for names?
italic	logical, use italic formatting for names?
font	character, font type for names (has no effect if format.names is TRUE)
family	character, font family for names
adj	numeric, adjustment for line labels
dx	numeric, x offset for line or field labels
dy	numeric, y offset for line or field labels
srt	numeric, rotation for line labels
min.area	numeric, minimum area of fields that should be labeled, expressed as a fraction of the total plot area
main	character, a main title for the plot; NULL means to plot no title
legend.x	character, description of legend placement passed to legend
add	logical, add to current plot?
plot.it	logical, make a plot?
tplot	logical, set up plot with thermo.plot.new?
x	matrix, value of the predominant list element from diagram
...	additional arguments passed to plot or barplot

Details

This function displays diagrams representing either chemical affinities or equilibrium chemical activities of species. The first argument is the output from affinity, rank.affinity, equilibrate, or solubility. 0-D diagrams, at a single point, are shown as barplots. 1-D diagrams, for a single variable on the x-axis, are plotted as lines. 2-D diagrams, for two variables, are plotted as predominance fields. The allowed variables are any that affinity or the other functions accepts: temperature, pressure, or the chemical activities of the basis species.

A new plot is started unless add is TRUE. If plot.it is FALSE, no plot will be generated but all the intermediate computations will be performed and the results returned.

Line or field labels use the names of the species as provided in eout; formatting is applied to chemical formulas (see expr.species) unless format.names is FALSE. Set names to TRUE or NULL to plot the names, or FALSE, NA, or "" to prevent plotting the names, or a character argument to replace the default species names. Alternatively, supply a numeric value to names to indicate a subset of default names that should or shouldn't be plotted (positive and negative indices, respectively). Use col.names and cex.names to change the colors and size of the labels. Use cex and cex.axis to adjust the overall character expansion factors (see par) and those of the axis labels. The x- and y-axis labels are automatically generated unless they are supplied in xlab and ylab.

If groups is supplied, the activities of the species identified in each numeric element of this list are multiplied by the balance coefficients of the species, then summed together. The names of the list are used to label the lines or fields for the summed activities of the resulting groups.

Normalizing protein formulas by length gives “residue equivalents” (Dick and Shock, 2011) that are useful for equilibrium calculations with proteins. normalize and as.residue are only usable when eout is the output from affinity, and only one can be TRUE. If normalize is TRUE, formation reactions and their affinities are first divided by protein length, so equal activities of residue equivalents are considered; then, the residue activities are rescaled to whole proteins for making the plot. If as.residue is TRUE, no rescaling is performed, so the diagram represents activities of the residues, not the whole proteins.

type argument

This paragraph describes the effect of the type argument when the output from affinity is being used. For type set to ‘⁠auto⁠’, and with 0 or 1 variables defined in affinity, the property computed by affinity for each species is plotted. This is usually the affinity of the formation reactions, but can be set to some other property (using the property argument of affinity), such as the equilibrium constant (‘⁠logK⁠’). For two variables, equilibrium predominance (maximum affinity) fields are plotted. This “maximum affinity method” (Dick, 2019) uses balancing coefficients that are specified by the balance argument. If type is ‘⁠saturation⁠’, the function plots the line for each species where the affinity of formation equals zero (see demo("saturation") for an example). If for a given species no saturation line is possible or the range of the diagram does not include the saturation line, the function prints a message instead. If type is the name of a basis species, then the equilibrium activity of the selected basis species in each of the formation reactions is plotted (see the CO₂-acetic acid example in buffer). In the case of 2-D diagrams, both of these options use contour to draw the lines, with the method specified in contour.method.

This paragraph describes the effect of the type argument when the output from solubility is being used. For one mineral or gas, if type set to ‘⁠auto⁠’, the equilibrium activities of each aqueous species are plotted. If type is ‘⁠loga.balance⁠’, the activity of the balancing basis species (i.e. total solubility) is plotted; this is represented by contours on a 2-D diagram. For two or more minerals or gases, if type set to ‘⁠auto⁠’, the values of ‘⁠loga.balance⁠’ (overall minimum solubility) are plotted. If type is ‘⁠loga.equil⁠’, the solubilities of the individual minerals and gases are plotted. For examples that use these features, see solubility and various demos: ‘⁠DEW⁠’, ‘⁠contour⁠’, ‘⁠gold⁠’, ‘⁠solubility⁠’, ‘⁠sphalerite⁠’.

1-D diagrams

For 1-D diagrams, the default setting for the y-axis is a logarithmic scale (unless alpha is TRUE) with limits corresponding to the range of logarithms of activities (or 0,1 if alpha is TRUE); these actions can be overridden by ylog and ylim. If legend.x is NA (the default), the lines are labeled with the names of the species near the maximum value. Otherwise, a legend is placed at the location identified by legend.x, or omitted if legend.x is NULL.

If alpha is TRUE, the fractional degrees of formation (ratios of activities to total activity) are plotted. Or, setting alpha to ‘⁠balance⁠’ allows the activities to be multiplied by the number of the balancing component; this is useful for making “percent carbon” diagrams where the species differ in carbon number. The line type and line width can be controlled with lty and lwd, respectively. To connect the points with splines instead of lines, set spline.method to one of the methods in splinefun.

2-D diagrams

On 2-D diagrams, the fields represent the species with the highest equilibrium activity. fill determines the color of the predominance fields, col that of the boundary lines. The default of NULL for fill uses a light blue, light tan, and darker tan color for fields with aqueous species, one solid, or two solids. These correspond to the web colors "aliceblue", "antiquewhite", and "burlywood" with some transparency added; see multi-metal for an example with two solids produced using mix. If all the species in the diagram have the same state, or if the fill argument is NA or a 0-length value, the predominance fields are transparent, i.e. no fill color is used. Otherwise, fill can be any colors, or the word ‘⁠rainbow⁠’, ‘⁠heat⁠’, ‘⁠terrain⁠’, ‘⁠topo⁠’, or ‘⁠cm⁠’, indicating a palette from grDevices. Starting with R version 3.6.0, fill can be the name of any available HCL color palette, matched in the same way as the palette argument of hcl.colors.

fill.NA gives the color for empty fields, i.e. points at which NA values are present for any species. This may occur when there are missing thermodynamic data or the temperature or pressure are not in the range of the equations of state. To make overlay diagrams easier to construct, the default for fill.NA is automatically changed to ‘⁠transparent⁠’ when add is TRUE.

If limit.water is TRUE, the diagram is clipped to the the water stability region on Eh-pH (and some other) diagrams. That is, predominance fields are shown only where water is stable, and fill.NA is used for areas where H₂O is not stable. The default of NULL for limit.water does not clip the main diagram but instead overlays it on the water stability fields. Change limit.water to FALSE to not show the water stability regions at all; this is automatically done if limit.water is missing and add is TRUE.

The default line-drawing algorithm uses contourLines to obtain smooth-looking diagonal and curved lines, at the expense of not coinciding exactly with the rectangular grid that is used for drawing colors. lty, col, and lwd can be specified, but limiting the lines via xrange is not currently supported. Set lty.cr or lty.aq to 0 to suppress boundary lines between minerals or aqueous species.

To go back to the old behavior for drawing lines, set dotted to ‘⁠0⁠’. The old behavior does not respect lty; instead, the style of the boundary lines on 2-D diagrams can be altered by supplying one or more non-zero integers in dotted, which indicates the fraction of line segments to omit; a value of ‘⁠1⁠’ or NULL for dotted has the effect of not drawing the boundary lines.

Activity Coefficients

The wording in this page and names of variables in functions refer exclusively to ‘⁠activities⁠’ of aqueous species. However, if activity coefficients are calculated (using the IS argument in affinity), then these variables are effectively transformed to molalities (see inst/tinytest/test-logmolality.R). So that the labels on diagrams are adjusted accordingly, diagram sets the molality argument of axis.label to TRUE if IS was supplied as an argument to affinity. The labeling as molality takes effect even if IS is set to 0; this way, by including (or not) the IS = 0 argument to affinity, the user decides whether to label aqueous species variables as molality (or activity) for calculations at zero ionic strength (where molality = activity).

Other Functions

find.tp finds the locations in a matrix of integers that are surrounded by the greatest number of different values. The function counts the unique values in a 3x3 grid around each point and returns a matrix of indices (similar to which(..., arr.ind = TRUE)) for the maximum count (ties result in more than one pair of indices). It can be used with the output from diagram for calculations in 2 dimensions to approximately locate the triple points on the diagram.

Value

diagram returns an invisible list containing, first, the contents of eout, i.e. the output of affinity or equilibrate supplied in the function call. To this are added the names of the plotted variable in plotvar, the labels used for species (which may be plotmath expressions if format.names is TRUE) in names, and the values used for plotting in a list named plotvals. For 1-D diagrams, plotvals usually corresponds to the chemical activities of the species (i.e. eout$loga.equil), or, if alpha is TRUE, their mole fractions (degrees of formation). For 2-D diagrams, plotvals corresponds to the values of affinity (from eout$values) divided by the respective balancing coefficients for each species. For 2-D diagrams, the output also contains the matrices predominant, which identifies the predominant species in eout$species at each grid point, and predominant.values, which has the affinities of the predominant species divided by the balancing coefficients (if eout is the output of affinity) or the activities of the predominant species (if eout is the output of equilibrate). The rows and columns of these matrices correspond to the x and y variables, respectively. Finally, the output for 2-D diagrams contains a lines component, giving the x- and y-coordinates of the field boundaries computed using contourLines; the values are padded to equal length with NAs to faciliate exporting the results using write.csv.

References

Aksu, S. and Doyle, F. M. (2001) Electrochemistry of copper in aqueous glycine solutions. J. Electrochem. Soc. 148, B51–B57.

Dick, J. M. (2019) CHNOSZ: Thermodynamic calculations and diagrams for geochemistry. Front. Earth Sci. 7:180. doi:10.3389/feart.2019.00180

Dick, J. M. and Shock, E. L. (2011) Calculation of the relative chemical stabilities of proteins as a function of temperature and redox chemistry in a hot spring. PLOS One 6, e22782. doi:10.1371/journal.pone.0022782

Helgeson, H. C. (1970) A chemical and thermodynamic model of ore deposition in hydrothermal systems. Mineral. Soc. Amer. Spec. Pap. 3, 155–186. https://www.worldcat.org/oclc/583263

Helgeson, H. C., Delany, J. M., Nesbitt, H. W. and Bird, D. K. (1978) Summary and critique of the thermodynamic properties of rock-forming minerals. Am. J. Sci. 278-A, 1–229. https://www.worldcat.org/oclc/13594862

LaRowe, D. E. and Helgeson, H. C. (2007) Quantifying the energetics of metabolic reactions in diverse biogeochemical systems: electron flow and ATP synthesis. Geobiology 5, 153–168. doi:10.1111/j.1472-4669.2007.00099.x

Majzlan, J., Navrotsky, A., McClesky, R. B. and Alpers, C. N. (2006) Thermodynamic properties and crystal structure refinement of ferricopiapite, coquimbite, rhomboclase, and Fe₂(SO₄)₃(H₂O)₅. Eur. J. Mineral. 18, 175–186. doi:10.1127/0935-1221/2006/0018-0175

Tagirov, B. and Schott, J. (2001) Aluminum speciation in crustal fluids revisited. Geochim. Cosmochim. Acta 65, 3965–3992. doi:10.1016/S0016-7037(01)00705-0

See Also

Berman, mix, mosaic, nonideal, solubility, and util.plot are other help topics that use diagram in their examples. See the demos for even more examples.

Examples

## Calculate the equilibrium logarithm of activity of a 
## basis species in different reactions
basis("CHNOS")
species(c("ethanol", "lactic acid", "deoxyribose", "ribose"))
a <- affinity(T = c(0, 150))
diagram(a, type = "O2", legend.x = "topleft", col = rev(rainbow(4)), lwd = 2)
title(main = "Equilibrium logfO2 for 1e-3 mol/kg of CO2 and ... ")

### 1-D diagrams: logarithms of activities

## Degrees of formation of ionized forms of glycine
## After Fig. 1 of Aksu and Doyle, 2001 
basis("CHNOS+")
species(ispecies <- info(c("glycinium", "glycine", "glycinate")))
a <- affinity(pH = c(0, 14))
e <- equilibrate(a)
diagram(e, alpha = TRUE, lwd = 1)
title(main = paste("Degrees of formation of aqueous glycine species\n",
  "after Aksu and Doyle, 2001"))

## Degrees of formation of ATP species as a function of 
## temperature, after LaRowe and Helgeson, 2007, Fig. 10b
# to make a similar diagram, activity of Mg+2 here is set to
# 10^-4, which is different from LH07, who used 10^-3 total molality
basis(c("CO2", "NH3", "H2O", "H3PO4", "O2", "H+", "Mg+2"),
  c(999, 999, 999, 999, 999, -5, -4))
species(c("HATP-3", "H2ATP-2", "MgATP-2", "MgHATP-"))
a <- affinity(T = c(0, 120, 25))
e <- equilibrate(a)
diagram(e, alpha = TRUE)  
title(main = paste("Degrees of formation of ATP species,\n",
  "pH=5, log(aMg+2)=-3. After LaRowe and Helgeson, 2007"),
  cex.main = 0.9)

### 2-D diagrams: predominance diagrams
### These use the maximum affinity method

## Fe-S-O at 200 deg C, after Helgeson, 1970
basis(c("Fe", "oxygen", "S2"))
species(c("iron", "ferrous-oxide", "magnetite",
  "hematite", "pyrite", "pyrrhotite"))
# The calculations include the polymorphic transitions of
# pyrrhotite; no additional step is needed
a <- affinity(S2 = c(-50, 0), O2 = c(-90, -10), T=200)
diagram(a, fill = "heat")
title(main = paste("Fe-S-O, 200 degrees C, 1 bar",
  "After Helgeson, 1970", sep = "\n"))

## pe-pH diagram for hydrated iron sulfides,
## goethite and pyrite, after Majzlan et al., 2006
basis(c("Fe+2", "SO4-2", "H2O", "H+", "e-"), 
  c(0, log10(3), log10(0.75), 999, 999))
species(c("rhomboclase", "ferricopiapite", "hydronium jarosite",
  "goethite", "melanterite", "pyrite"))
a <- affinity(pH = c(-1, 4, 256), pe = c(-5, 23, 256))
d <- diagram(a, main = "Fe-S-O-H, after Majzlan et al., 2006")
water.lines(d, lwd = 2)
text(3, 22, describe.basis(2:3, digits = 2, oneline = TRUE))
text(3, 21, describe.property(c("T", "P"), c(25, 1), oneline = TRUE))

## Aqueous Al species, after Tagirov and Schott, 2001
basis(c("Al+3", "F-", "H+", "O2", "H2O"))
AlOH <- c("Al(OH)4-", "Al(OH)3", "Al(OH)2+", "AlOH+2")
Al <- "Al+3"
AlF <- c("AlF+2", "AlF2+", "AlF3", "AlF4-")
AlOHF <- c("Al(OH)2F2-", "Al(OH)2F", "AlOHF2")
species(c(AlOH, Al, AlF, AlOHF), "aq")
res <- 300
a <- affinity(pH = c(0.5, 6.5, res), `F-` = c(-2, -9, res), T = 200)
diagram(a, fill = "terrain")
dprop <- describe.property(c("T", "P"), c(200, "Psat"))
legend("topright", legend = dprop, bty = "n")
mtitle(c("Aqueous aluminum species",
         "After Tagirov and Schott, 2001 Fig. 4d"), cex = 0.95)

## Temperature-Pressure: kayanite-sillimanite-andalusite
# cf. Fig. 49 of Helgeson et al., 1978
# this is a system of one component (Al2SiO5), however:
# - number of basis species must be the same as of elements
# - avoid using H2O or other aqueous species because of
#     T/P limits of the water() calculations;
basis(c("corundum", "quartz", "oxygen"))
species(c("kyanite", "sillimanite", "andalusite"))
# Database has transition temperatures of kyanite and andalusite
# at 1 bar only, so we permit calculation at higher temperatures
a <- affinity(T = c(200, 900, 99), P = c(0, 9000, 101), exceed.Ttr = TRUE)
d <- diagram(a, fill = NULL)
slab <- syslab(c("Al2O3", "SiO2", "H2O"))
mtitle(c(as.expression(slab), "after Helgeson et al., 1978"))
# Find the approximate position of the triple point
tp <- find.tp(d$predominant)
Ttp <- a$vals[[1]][tp[1, 2]]
Ptp <- rev(a$vals[[2]])[tp[1, 1]]
points(Ttp, Ptp, pch = 10, cex = 5)

---

Regress equations-of-state parameters for aqueous species

Description

Fit experimental volumes and heat capacities using regression equations. Possible models include the Helgeson-Kirkham-Flowers (HKF) equations of state, or other equations defined using any combination of terms derived from the temperature, pressure and thermodynamic and electrostatic properties of water.

Usage

EOSregress(exptdata, var = "", T.max = 9999, ...)
  EOSvar(var, T, P, ...)
  EOScalc(coefficients, T, P, ...)
  EOSplot(exptdata, var = NULL, T.max = 9999, T.plot = NULL, 
    fun.legend = "topleft", coefficients = NULL, add = FALSE,
    lty = par("lty"), col=par("col"), ...)
  EOSlab(var, coeff = "")
  EOScoeffs(species, property, P=1)
  Cp_s_var(T = 298.15, P = 1, omega.PrTr = 0, Z = 0)
  V_s_var(T = 298.15, P = 1, omega.PrTr = 0, Z = 0)

Arguments

exptdata	dataframe, experimental data
var	character, name(s) of variables in the regression equations
T.max	numeric, maximum temperature for regression, in degrees Kelvin
T	numeric, temperature in Kelvin
P	numeric, pressure in bars
...	arguments specifying additional dependencies of the regression variables
T.plot	numeric, upper limit of temperature range to plot
fun.legend	character, where to place legend on plot
coefficients	dataframe, coefficients to use to make line on plot
add	logical, add lines to an existing plot?
lty	line style
col	color of lines
coeff	numeric, value of equation of state parameter for plot legend
species	character, name of aqueous species
property	character, ‘⁠Cp⁠’ or ‘⁠V⁠’
omega.PrTr	numeric, value of omega at reference T and P
Z	numeric, charge

Details

EOSregress uses a linear model (lm) to regress the experimental heat capacity or volume data in exptdata, which is a data frame with columns ‘⁠T⁠’ (temperature in degrees Kelvin), ‘⁠P⁠’ (pressure in bars), and ‘⁠Cp⁠’ or ‘⁠V⁠’ (heat capacity in cal/mol.K or volume in cm3/mol). The ‘⁠Cp⁠’ or ‘⁠V⁠’ data must be in the third column. Only data below the temperature of T.max are included in the regression. The regression formula is specified by a vector of names in var. The names of the variables can be any combination of the following (listed in the order of search): variables listed in the following table, any available property of water (e.g. ‘⁠V⁠’, ‘⁠alpha⁠’, ‘⁠QBorn⁠’), or the name of a function that can be found using get in the default environment. Examples of the latter are Cp_s_var, V_s_var, or functions defined by the user in the global environment; the arguments of these functions must include, but are not limited to, T and P.

T	$T$ (temperature)
P	$P$ (pressure)
TTheta	$(T-\Theta)$ ($\Theta$ = 228 K)
invTTheta	$1/(T-\Theta)$
TTheta2	$(T-\Theta)^2$
invTTheta2	$1/(T-\Theta)^2$
invPPsi	$1/(P+\Psi)$ ($\Psi$ = 2600 bar)
invPPsiTTheta	$1/((P+\Psi)(T-\Theta))$
TXBorn	$TX$ (temperature times $X$ Born function)
drho.dT	$d\rho/dT$ (temperature derivative of density of water)
V.kT	$V\kappa_T$ (volume times isothermal compressibility of water)

EOSvar calculates the value of the variable named var (defined as described above) at the specified T (temperature in degrees Kelvin) and P (pressure in bars). This function is used by EOSregress to get the values of the variables used in the regression.

EOScalc calculates the predicted heat capacities or volumes using coefficients provided by the result of EOSregress, at the temperatures and pressures specified by T and P.

EOSplot takes a table of data in exptdata, runs EOSregress and EOScalc and plots the results. The experimental data are plotted as points, and the calculated values as a smooth line. The point symbols are filled circles where the calculated value is within 10% of the experimental value; open circles otherwise.

EOSlab produces labels for the variables listed above that can be used as.expressions in plots. The value of coeff is prefixed to the name of the variable (using substitute, with a multiplication symbol). For the properties listed in the table above, and selected properties listed in water, the label is formatted using plotmath expressions (e.g., with italicized symbols and Greek letters). If var is a user-defined function, the function can be given a ‘⁠label⁠’ attribute to provide plotmath-style formatting; in this case the appropriate multiplication or division symbol should be specified (see example below).

EOScoeffs retrieves coefficients in the Helgeson-Kirkham-Flowers equations from the thermodynamic database (thermo$OBIGT) for the given aqueous species. If the property is ‘⁠Cp⁠’, the resulting data frame has column names of ‘⁠(Intercept)⁠’, ‘⁠invTTheta2⁠’ and ‘⁠TX⁠’, respectively holding the coefficients $c_1$, $c_2$ and $\omega$ in the equation $Cp^\circ = c_1 + c_2/(T-\Theta)^2 + {\omega}TX$. If the property is ‘⁠V⁠’, the data frame has column names of ‘⁠(Intercept)⁠’, ‘⁠invTTheta⁠’ and ‘⁠Q⁠’, respectively holding the coefficients $\sigma$, $\xi$ and $\omega$ in $V^\circ = \sigma + \xi/(T-\Theta) - {\omega}Q$. Here, $\sigma$ and $\xi$ are calculated from a1, a2, a3 and a4 in thermo()$OBIGT at the pressure indicated by P (default 1 bar).

The original motivation for writing these functions was to explore alternatives or possible modifications to the revised Helgeson-Kirkham-Flowers equations applied to aqueous nonelectrolytes. As pointed out by Schulte et al., 2001, the functional forms of the equations do not permit retrieving values of the solvation parameter ($\omega$) that closely represent the observed trends in both heat capacity and volume at high temperatures (above ca. 200 °C).

The examples below assume that the $\omega$ parameter in the HKF functions is a constant (does not depend on T and P), as is appropriate for nonelectrolytes. For charged species, the variables Cp_s_var and V_s_var can be used in the regressions. They correspond to the solvation contribution to heat capacity or volume, respectively, in the HKF EOS, divided by the value of $\omega$ at the reference temperature and pressure. Because these variables are themselves a function of omega.PrTr, an iterative procedure is needed to perform the regression.

Note that variables QBorn and V_s_var are both negated, so that the value of $\omega$ has its proper sign in the corresponding equations.

Value

For EOSregress, an object of class “lm”. EOSvar and EOScalc both return numeric values. EOScoeffs returns a data frame.

References

Hnědkovský, L. and Wood, R. H. (1997) Apparent molar heat capacities of aqueous solutions of CH₄, CO₂, H₂S, and NH₃ at temperatures from 304 K to 704 K at a pressure of 28 MPa. J. Chem. Thermodyn. 29, 731–747. doi:10.1006/jcht.1997.0192

Schulte, M. D., Shock, E. L. and Wood, R. H. (1995) The temperature dependence of the standard-state thermodynamic properties of aqueous nonelectrolytes. Geochim. Cosmochim. Acta 65, 3919–3930. doi:10.1016/S0016-7037(01)00717-7

See Also

The vignette eos-regress has more references and examples, including an iterative method to retrieve omega.PrTr.

Examples

## Fit experimental heat capacities of CH4
## using revised Helgeson-Kirkham-Flowers equations
# Read the data from Hnedkovsky and Wood, 1997
f <- system.file("extdata/cpetc/HW97_Cp.csv", package = "CHNOSZ")
d <- read.csv(f)
# Use data for CH4
d <- d[d$species == "CH4", ]
d <- d[, -1]
# Convert J to cal and MPa to bar
d$Cp <- convert(d$Cp, "cal")
d$P <- convert(d$P, "bar")
# Specify the terms in the HKF equations
var <- c("invTTheta2", "TXBorn")
# Perform regression, with a temperature limit
EOSlm <- EOSregress(d, var, T.max = 600)
# Calculate the Cp at some temperature and pressure
EOScalc(EOSlm$coefficients, 298.15, 1)
# Get the database values of c1, c2 and omega for CH4(aq)
CH4coeffs <- EOScoeffs("CH4", "Cp")
## Make plots comparing the regressions
## with the accepted EOS parameters of CH4
opar <- par(mfrow = c(2,2))
EOSplot(d, T.max = 600)
title("Cp of CH4(aq), fit to 600 K")
legend("bottomleft", pch = 1, legend = "Hnedkovsky and Wood, 1997")
EOSplot(d, coefficients = CH4coeffs)
title("Cp from EOS parameters in database")
EOSplot(d, T.max = 600, T.plot = 600)
title("Cp fit to 600 K, plot to 600 K")
EOSplot(d, coefficients = CH4coeffs, T.plot = 600)
title("Cp from EOS parameters in database")
par(opar)

# Continuing from above, with user-defined variables
Theta <- 228  # K 
invTTTheta3 <- function(T, P) (2*T) / (T-T*Theta) ^ 3
invTX <- function(T, P) 1 / T * water("XBorn", T = T, P = P)[,1]
# Print the calculated values of invTTTheta3
EOSvar("invTTTheta3", d$T, d$P)
# Use invTTTheta and invTX in a regression
var <- c("invTTTheta3", "invTX")
EOSregress(d, var)
# Give them a "label" attribute for use in the legend
attr(invTTTheta3, "label") <-
  quote(phantom()%*%2 * italic(T) / (italic(T) - italic(T) * Theta) ^ 3)
attr(invTX, "label") <- quote(phantom() / italic(T * X))
# Uncomment the following to make the plot
#EOSplot(d, var)

## Model experimental volumes of CH4
## using HKF equation and an exploratory one
f <- system.file("extdata/cpetc/HWM96_V.csv", package = "CHNOSZ")
d <- read.csv(f)
# Use data for CH4 near 280 bar
d <- d[d$species == "CH4", ]
d <- d[, -1]
d <- d[abs(d$P - 28) < 0.1, ]
d$P <- convert(d$P, "bar")
# The HKF equation
varHKF <- c("invTTheta", "QBorn")
# alpha is the expansivity coefficient of water
varal <- c("invTTheta", "alpha")
opar <- par(mfrow = c(2, 2))
# For both HKF and the expansivity equation,
# we'll fit up to a temperature limit
EOSplot(d, varHKF, T.max = 663, T.plot = 625)
legend("bottomright", pch = 1, legend = "Hnedkovsky et al., 1996")
title("V of CH4(aq), HKF equation")
EOSplot(d, varal, T.max = 663, T.plot = 625)
title("V of CH4(aq), expansivity equation")
EOSplot(d, varHKF, T.max = 663)
title("V of CH4(aq), HKF equation")
EOSplot(d, varal, T.max = 663)
title("V of CH4(aq), expansivity equation")
par(opar)
# Note that the volume regression using the HKF gives
# a result for omega (coefficient on Q) that is
# not consistent with the high-T heat capacities

---

Equilibrium chemical activities of species

Description

Calculate equilibrium chemical activities of species from the affinities of formation of the species at unit activity.

Usage

equilibrate(aout, balance = NULL, loga.balance = NULL,
    ispecies = !grepl("cr", aout$species$state), normalize = FALSE,
    as.residue = FALSE, method = c("boltzmann", "reaction"),
    tol = .Machine$double.eps^0.25)
  equil.boltzmann(Astar, n.balance, loga.balance)
  equil.reaction(Astar, n.balance, loga.balance, tol = .Machine$double.eps^0.25)
  moles(eout)

Arguments

aout	list, output from affinity

or mosaic

balance	character or numeric, how to balance the transformations
ispecies	numeric, which species to include
normalize	logical, normalize the molar formulas of species by the balancing coefficients?
as.residue	logical, report results for the normalized formulas?
Astar	numeric, affinities of formation reactions excluding species contribution
n.balance	numeric, number of moles of balancing component in the formation reactions of the species of interest
loga.balance	numeric (single value or vector), logarithm of total activity of balanced quantity
method	character, equilibration method to use
tol	numeric, convergence tolerance for uniroot
eout	list, output from equilibrate

Details

equilibrate calculates the chemical activities of species in metastable equilibrium, for constant temperature, pressure and chemical activities of basis species, using specified balancing constraints on reactions between species.

It takes as input aout, the output from affinity, giving the chemical affinities of formation reaction of each species, which may be calculated on a multidimensional grid of conditions. Alternatively, aout can be the output from mosaic, in which case the equilibrium activities of the formed species are calculated and combined with those of the changing basis species to make an object that can be plotted with diagram.

The equilibrium chemical activities of species are calculated using either the equil.reaction or equil.boltzmann functions, the latter only if the balance is on one mole of species.

equilibrate needs to be provided constraints on how to balance the reactions representing transformations between the species. balance indicates the balancing component, according to the following scheme:

• ‘⁠NULL⁠’: autoselect

• name of basis species: balance on this basis species

• ‘⁠length⁠’: balance on length of proteins

• ‘⁠1⁠’: balance on one mole of species

• numeric vector: user-defined constraints

The default value of NULL for balance indicates to use the coefficients on the basis species that is present (i.e. with non-zero coefficients) in all formation reactions, or if that fails, to set the balance to ‘⁠1⁠’. However, if all the species (as listed in code aout$species) are proteins (have an underscore character in their names), the default value of NULL for balance indicates to use ‘⁠length⁠’ as the balance.

NOTE: The summation of activities assumes an ideal system, so ‘molality’ is equivalent to ‘activity’ here. loga.balance gives the logarithm of the total activity of balance (which is total activity of species for ‘⁠1⁠’ or total activity of amino acid residue-equivalents for ‘⁠length⁠’). If loga.balance is missing, its value is taken from the activities of species listed in aout; this default is usually the desired operation. The supplied value of loga.balance may also be a vector of values, with length corresponding to the number of conditions in the calculations of affinity.

normalize if TRUE indicates to normalize the molar formulas of species by the balance coefficients. This operation is intended for systems of proteins, whose conventional formulas are much larger than the basis speices. The normalization also applies to the balancing coefficients, which as a result consist of ‘⁠1⁠’s. After normalization and equilibration, the equilibrium activities are then re-scaled (for the original formulas of the species), unless as.residue is TRUE.

equil.boltzmann is used to calculate the equilibrium activities if balance is ‘⁠1⁠’ (or when normalize or as.residue is TRUE), otherwise equil.reaction is called. The default behavior can be overriden by specifying either ‘⁠boltzmann⁠’ or ‘⁠reaction⁠’ in method. Using equil.reaction may be needed for systems with huge (negative or positive) affinities, where equil.boltzmann produces a NaN result.

ispecies can be supplied to identify a subset of the species to include in the equilibrium calculation. By default, this is all species except solids (species with ‘⁠cr⁠’ state). However, the stability regions of solids are still calculated (by a call to diagram without plotting). At all points outside of their stability region, the logarithms of activities of solids are set to -999. Likewise, where any solid species is calculated to be stable, the logarithms of activities of all aqueous species are set to -999.

moles simply calculates the total number of moles of elements corresponding to the activities of formed species in the output from equilibrate.

Value

equil.reaction and equil.boltzmann each return a list with dimensions and length equal to those of Astar, giving the log10 of the equilibrium activities of the species of interest. equilibrate returns a list, containing first the values in aout, to which are appended m.balance (the balancing coefficients if normalize is TRUE, a vector of ‘⁠1⁠’s otherwise), n.balance (the balancing coefficients if normalize is FALSE, a vector of ‘⁠1⁠’s otherwise), loga.balance, Astar, and loga.equil (the calculated equilibrium activities of the species).

Algorithms

The input values to equil.reaction and equil.boltzmann are in a list, Astar, all elements of the list having the same dimensions; they can be vectors, matrices, or higher-dimensionsal arrays. Each list element contains the chemical affinities of the formation reactions of one of the species of interest (in dimensionless base-10 units, i.e. A/2.303RT), calculated at unit activity of the species of interest. The equilibrium base-10 logarithm activities of the species of interest returned by either function satisfy the constraints that 1) the final chemical affinities of the formation reactions of the species are all equal and 2) the total activity of the balancing component is equal to (loga.balance). The first constraint does not impose a complete equilibrium, where the affinities of the formation reactions are all equal to zero, but allows for a metastable equilibrium, where the affinities of the formation reactions are equal to each other.

In equil.reaction (the algorithm described in Dick, 2008 and the only one available prior to CHNOSZ_0.8), the calculations of relative abundances of species are based on a solving a system of equations representing the two constraints stated above. The solution is found using uniroot with a flexible method for generating initial guesses.

In equil.boltzmann, the chemical activities of species are calculated using the Boltzmann distribution. This calculation is faster than the algorithm of equil.reaction, but is limited to systems where the transformations are all balanced on one mole of species. If equil.boltzmann is called with balance other than ‘⁠1⁠’, it stops with an error.

Warning

Despite its name, this function does not generally produce a complete equilibrium. It returns activities of species such that the affinities of formation reactions are equal to each other (and transformations between species have zero affinity); this is a type of metastable equilibrium. Although they are equal to each other, the affinities are not necessarily equal to zero. Use solubility to find complete equilibrium, where the affinities of the formation reactions become zero.

References

Dick, J. M. (2008) Calculation of the relative metastabilities of proteins using the CHNOSZ software package. Geochem. Trans. 9:10. doi:10.1186/1467-4866-9-10

See Also

diagram has examples of using equilibrate to make equilibrium activity diagrams. palply is used by both equil.reaction and equil.boltzmann to parallelize intensive parts of the calculations.

See the vignette multi-metal for an example of balancing on two elements (N in the basis species, C in the formed species).

Examples

## Equilibrium in a simple system:
## ionization of carbonic acid
basis("CHNOS+")
species(c("CO2", "HCO3-", "CO3-2"))
# Set unit activity of the species (0 = log10(1))
species(1:3, 0)
# Calculate Astar (for unit activity)
res <- 100
Astar <- affinity(pH = c(0, 14, res))$values
# The logarithms of activity for a total activity
# of the balancing component (CO2) equal to 0.001
loga.boltz <- equil.boltzmann(Astar, c(1, 1, 1), 0.001)
# Calculated another way
loga.react <- equil.reaction(Astar, c(1, 1, 1), rep(0.001, 100))
# They should be pretty close
stopifnot(all.equal(loga.boltz, loga.react))
# The first ionization constant (pKa)
ipKa <- which.min(abs(loga.boltz[[1]] - loga.boltz[[2]]))
pKa.equil <- seq(0, 14, length.out = res)[ipKa]
# Calculate logK directly
logK <- subcrt(c("CO2","H2O","HCO3-","H+"), c(-1, -1, 1, 1), T = 25)$out$logK
# We could decrease tolerance here by increasing res
stopifnot(all.equal(pKa.equil, -logK, tolerance = 1e-2))

---

Run examples from the documentation

Description

Run the examples contained in each of the documentation topics.

Usage

examples(save.png = FALSE)
  demos(which = c("sources", "protein.equil", "affinity", "NaCl",
    "density", "ORP", "ionize", "buffer", "protbuff",
    "glycinate", "mosaic", "copper", "arsenic", "solubility", "gold",
    "contour", "sphalerite", "minsol", "Shh", "saturation",
    "adenine", "DEW", "lambda", "potassium", "TCA", "aluminum", "AD",
    "comproportionation", "Pourbaix", "E_coli", "yttrium", "rank.affinity"),
    save.png = FALSE)

Arguments

save.png	logical, generate PNG image files for the plots?
which	character, which example to run

Details

examples runs all the examples in the help pages for the package. example is called for each topic with ask set to FALSE (so all of the figures are shown without prompting the user).

demos runs all the demos in the package. The demo(s) to run is/are specified by which; the default is to run them in the order of the list below.

sources

Cross-check the reference list with the thermodynamic database

protein.equil

Chemical activities of two proteins in metastable equilibrium (Dick and Shock, 2011)

affinity

Affinities of metabolic reactions and amino acid synthesis (Amend and Shock, 1998, 2001)

NaCl

Equilibrium constant for aqueous NaCl dissociation (Shock et al., 1992)

density

Density of H₂O, inverted from IAPWS-95 equations (rho.IAPWS95)

ORP

Temperature dependence of oxidation-reduction potential for redox standards

ionize

ionize.aa(): contour plots of net charge and ionization properties of LYSC_CHICK

buffer

Minerals and aqueous species as buffers of hydrogen fugacity (Schulte and Shock, 1995)

protbuff

Chemical activities buffered by thiol peroxidases or sigma factors

glycinate

Metal-glycinate complexes (Shock and Koretsky, 1995; Azadi et al., 2019)

mosaic

Eh-pH diagram with two sets of changing basis species (Garrels and Christ, 1965)

copper

Another example of mosaic: complexation of Cu with glycine (Aksu and Doyle, 2001)

arsenic

Another example of mosaic: Eh-pH diagram for the system As-O-H-S (Lu and Zhu, 2011)

solubility

Solubility of calcite (cf. Manning et al., 2013) and CO₂ (cf. Stumm and Morgan, 1996)

gold

Solubility of gold (Akinfiev and Zotov; 2001; Stefánsson and Seward, 2004; Williams-Jones et al., 2009)

contour

Gold solubility contours on a log fO2 - pH diagram (Williams-Jones et al., 2009)

sphalerite

Solubility of sphalerite (Akinfiev and Tagirov, 2014)

minsol

Solubilities of multiple minerals

dehydration

log K

of dehydration reactions; SVG file contains tooltips and links

Shh

Affinities of transcription factors relative to Sonic hedgehog (Dick, 2015)

saturation

Equilibrium activity diagram showing activity ratios and mineral saturation limits (Bowers et al., 1984)

adenine

HKF regression of heat capacity and volume of aqueous adenine (Lowe et al., 2017)

DEW

Deep Earth Water (DEW) model for high pressures (Sverjensky et al., 2014a and 2014b)

lambda

Effects of lambda transition on thermodynamic properties of quartz (Berman, 1988)

potassium

Comparison of thermodynamic datasets for predicting mineral stabilities (Sverjensky et al., 1991)

TCA

Standard Gibbs energies of the tricarboxylic (citric) acid cycle (Canovas and Shock, 2016)

aluminum

Reactions involving Al-bearing minerals (Zimmer et al., 2016; Tutolo et al., 2014)

carboxylase

Rank abundance distribution for RuBisCO and acetyl-CoA carboxylase

AD

Dissolved gases: Henry's constant, volume, and heat capacity (Akinfiev and Diamond, 2003)

comproportionation

Gibbs energy of sulfur comproportionation (Amend et al., 2020)

Pourbaix

Eh-pH diagram for Fe-O-H with equisolubility lines (Pourbaix, 1974)

E_coli

Gibbs energy of biomass synthesis in E. coli (LaRowe and Amend, 2016)

rank.affinity

Affinity ranking for proteins in yeast nutrient limitation (data from Tai et al., 2005)

yttrium

logB.to.OBIGT fits at 800 and 1000 bar and Y speciation in NaCl solution at varying pH (Guan et al., 2020)

For either function, if save.png is TRUE, the plots are saved in png files whose names begin with the names of the help topics or demos.

Two of the demos have external dependencies and are not automatically run by demos. ‘⁠dehydration⁠’ creates an interactive SVG file; this demo depends on RSVGTipsDevice, which is not available for Windows. ‘⁠carboxylase⁠’ creates an animated GIF; this demo requires that the ImageMagick convert commmand be available on the system (tested on Linux and Windows).

‘⁠carboxylase⁠’ animates diagrams showing rankings of calculated chemical activities along a combined T and loga_H2 gradient, or makes a single plot on the default device (without conversion to animated GIF) if a single temperature (T) is specified in the code. To run this demo, an empty directory named ‘⁠png⁠’ must be present (as a subdirectory of the R working directory). The proteins in the calculation are 24 carboxylases from a variety of organisms. There are 12 ribulose phosphate carboxylase and 12 acetyl-coenzyme A carboxylase; 6 of each type are from nominally mesophilic organisms and 6 from nominally thermophilic organisms, shown as blue and red symbols on the diagrams. The activities of hydrogen at each temperature are calculated using $\log a_{\mathrm{H_{2}}_{\left(aq\right)}}=-11+3/\left(40\times T\left(^{\circ}C\right)\right)$; this equation comes from a model of relative stabilities of proteins in a hot-spring environment (Dick and Shock, 2011).

Warning

The discontinuities apparent in the plot made by the NaCl demo illustrate limitations of the "g function" for charged species in the revised HKF model (the 355 °C boundary of region II in Figure 6 of Shock et al., 1992). Note that SUPCRT92 (Johnson et al., 1992) gives similar output at 500 bar. However, SUPCRT does not output thermodynamic properties above 350 °C at P_SAT; see Warning in subcrt.

References

Akinfiev, N. N. and Diamond, L. W. (2003) Thermodynamic description of aqueous nonelectrolytes at infinite dilution over a wide range of state parameters. Geochim. Cosmochim. Acta 67, 613–629. doi:10.1016/S0016-7037(02)01141-9

Akinfiev, N. N. and Tagirov, B. R. (2014) Zn in hydrothermal systems: Thermodynamic description of hydroxide, chloride, and hydrosulfide complexes. Geochem. Int. 52, 197–214. doi:10.1134/S0016702914030021

Akinfiev, N. N. and Zotov, A. V. (2001) Thermodynamic description of chloride, hydrosulfide, and hydroxo complexes of Ag(I), Cu(I), and Au(I) at temperatures of 25-500°C and pressures of 1-2000 bar. Geochem. Int. 39, 990–1006.

Aksu, S. and Doyle, F. M. (2001) Electrochemistry of copper in aqueous glycine solutions. J. Electrochem. Soc. 148, B51–B57.

Amend, J. P. and Shock, E. L. (1998) Energetics of amino acid synthesis in hydrothermal ecosystems. Science 281, 1659–1662. doi:10.1126/science.281.5383.1659

Amend, J. P. and Shock, E. L. (2001) Energetics of overall metabolic reactions of thermophilic and hyperthermophilic Archaea and Bacteria. FEMS Microbiol. Rev. 25, 175–243. doi:10.1016/S0168-6445(00)00062-0

Amend, J. P., Aronson, H. S., Macalady, J. and LaRowe, D. E. (2020) Another chemolithotrophic metabolism missing in nature: sulfur comproportionation. Environ. Microbiol. 22, 1971–1976. doi:10.1111/1462-2920.14982

Azadi, M. R., Karrech, A., Attar, M. and Elchalakani, M. (2019) Data analysis and estimation of thermodynamic properties of aqueous monovalent metal-glycinate complexes. Fluid Phase Equilib. 480, 25-40. doi:10.1016/j.fluid.2018.10.002

Berman, R. G. (1988) Internally-consistent thermodynamic data for minerals in the system Na₂O-K₂O-CaO-MgO-FeO-Fe₂O₃-Al₂O₃-SiO₂-TiO₂-H₂O-CO₂. J. Petrol. 29, 445-522. doi:10.1093/petrology/29.2.445

Bowers, T. S., Jackson, K. J. and Helgeson, H. C. (1984) Equilibrium Activity Diagrams for Coexisting Minerals and Aqueous Solutions at Pressures and Temperatures to 5 kb and 600°C, Springer-Verlag, Berlin, 397 p. https://www.worldcat.org/oclc/11133620

Canovas, P. A., III and Shock, E. L. (2016) Geobiochemistry of metabolism: Standard state thermodynamic properties of the citric acid cycle. Geochim. Cosmochim. Acta 195, 293–322. doi:10.1016/j.gca.2016.08.028

Dick, J. M. and Shock, E. L. (2011) Calculation of the relative chemical stabilities of proteins as a function of temperature and redox chemistry in a hot spring. PLOS One 6, e22782. doi:10.1371/journal.pone.0022782

Dick, J. M. (2015) Chemical integration of proteins in signaling and development. bioRxiv. doi:10.1101/015826

Garrels, R. M. and Christ, C. L. (1965) Solutions, Minerals, and Equilibria, Harper & Row, New York, 450 p. https://www.worldcat.org/oclc/517586

Guan, Q., Mei, Y., Etschmann, B., Testemale, D., Louvel, M. and Brugger, J. (2020) Yttrium complexation and hydration in chloride-rich hydrothermal fluids: A combined ab initio molecular dynamics and in situ X-ray absorption spectroscopy study. Geochim. Cosmochim. Acta 281, 168–189. doi:10.1016/j.gca.2020.04.015

Johnson, J. W., Oelkers, E. H. and Helgeson, H. C. (1992) SUPCRT92: A software package for calculating the standard molal thermodynamic properties of minerals, gases, aqueous species, and reactions from 1 to 5000 bar and 0 to 1000°C. Comp. Geosci. 18, 899–947. doi:10.1016/0098-3004(92)90029-Q

LaRowe, D. E. and Amend, J. P. (2016) The energetics of anabolism in natural settings. ISME J. 10, 1285–1295. doi:10.1038/ismej.2015.227

Lowe, A. R., Cox, J. S. and Tremaine, P. R. (2017) Thermodynamics of aqueous adenine: Standard partial molar volumes and heat capacities of adenine, adeninium chloride, and sodium adeninate from T = 278.15 K to 393.15 K. J. Chem. Thermodyn. 112, 129–145. doi:10.1016/j.jct.2017.04.005

Lu, P. and Zhu, C. (2011) Arsenic Eh–pH diagrams at 25°C and 1 bar. Environ. Earth Sci. 62, 1673–1683. doi:10.1007/s12665-010-0652-x

Manning, C. E., Shock, E. L. and Sverjensky, D. A. (2013) The chemistry of carbon in aqueous fluids at crustal and upper-mantle conditions: Experimental and theoretical constraints. Rev. Mineral. Geochem. 75, 109–148. doi:10.2138/rmg.2013.75.5

Pourbaix, M. (1974) Atlas of Electrochemical Equilibria in Aqueous Solutions, NACE, Houston, TX and CEBELCOR, Brussels. https://www.worldcat.org/oclc/563921897

Schulte, M. D. and Shock, E. L. (1995) Thermodynamics of Strecker synthesis in hydrothermal systems. Orig. Life Evol. Biosph. 25, 161–173. doi:10.1007/BF01581580

Shock, E. L. and Koretsky, C. M. (1995) Metal-organic complexes in geochemical processes: Estimation of standard partial molal thermodynamic properties of aqueous complexes between metal cations and monovalent organic acid ligands at high pressures and temperatures. Geochim. Cosmochim. Acta 59, 1497–1532. doi:10.1016/0016-7037(95)00058-8

Shock, E. L., Oelkers, E. H., Johnson, J. W., Sverjensky, D. A. and Helgeson, H. C. (1992) Calculation of the thermodynamic properties of aqueous species at high pressures and temperatures: Effective electrostatic radii, dissociation constants and standard partial molal properties to 1000 °C and 5 kbar. J. Chem. Soc. Faraday Trans. 88, 803–826. doi:10.1039/FT9928800803

Stefánsson, A. and Seward, T. M. (2004) Gold(I) complexing in aqueous sulphide solutions to 500°C at 500 bar. Geochim. Cosmochim. Acta 68, 4121–4143. doi:10.1016/j.gca.2004.04.006

Stumm, W. and Morgan, J. J. (1996) Aquatic Chemistry: Chemical Equilibria and Rates in Natural Waters, John Wiley & Sons, New York, 1040 p. https://www.worldcat.org/oclc/31754493

Sverjensky, D. A., Harrison, B. and Azzolini, D. (2014a) Water in the deep Earth: The dielectric constant and the solubilities of quartz and corundum to 60 kb and 1,200 °C. Geochim. Cosmochim. Acta 129, 125–145. doi:10.1016/j.gca.2013.12.019

Sverjensky, D. A., Hemley, J. J. and D'Angelo, W. M. (1991) Thermodynamic assessment of hydrothermal alkali feldspar-mica-aluminosilicate equilibria. Geochim. Cosmochim. Acta 55, 989-1004. doi:10.1016/0016-7037(91)90157-Z

Sverjensky, D. A., Stagno, V. and Huang, F. (2014b) Important role for organic carbon in subduction-zone fluids in the deep carbon cycle. Nat. Geosci. 7, 909–913. doi:10.1038/ngeo2291

Tai, S. L., Boer, V. M., Daran-Lapujade, P., Walsh, M. C., de Winde, J. H., Daran, J.-M. and Pronk, J. T. (2005) Two-dimensional transcriptome analysis in chemostat cultures: Combinatorial effects of oxygen availability and macronutrient limitation in Saccharomyces cerevisiae. J. Biol. Chem. 280, 437–447. doi:10.1074/jbc.M410573200

Tutolo, B. M., Kong, X.-Z., Seyfried, W. E., Jr. and Saar, M. O. (2014) Internal consistency in aqueous geochemical data revisited: Applications to the aluminum system. Geochim. Cosmochim. Acta 133, 216–234. doi:10.1016/j.gca.2014.02.036

Williams-Jones, A. E., Bowell, R. J. and Migdisov, A. A. (2009) Gold in solution. Elements 5, 281–287. doi:10.2113/gselements.5.5.281

Zimmer, K., Zhang, Y., Lu, P., Chen, Y., Zhang, G., Dalkilic, M. and Zhu, C. (2016) SUPCRTBL: A revised and extended thermodynamic dataset and software package of SUPCRT92. Comp. Geosci. 90, 97–111. doi:10.1016/j.cageo.2016.02.013

Examples

demos(c("ORP", "NaCl"))

---

Extra data

Description

The files in the subdirectories of extdata provide additional thermodynamic data and other data to support the examples in the package documentation and vignettes. See thermo for a description of the files in extdata/OBIGT, which are used to generate the thermodynamic database.

Details

Files in Berman contain thermodynamic data for minerals using the Berman formulation:

• Ber88_1988.csv contains thermodynamic data for minerals taken from Berman (1988).

• Other files with names like xxx_yyyy.csv contain thermodynamic data from other sources; xxx in the filename corresponds to the reference in thermo$OBIGT and yyyy gives the year of publication. Berman uses these data for the calculation of thermodynamic properties at specified P and T, which are then available for use in subcrt. If there are any duplicated mineral names in the files, only the most recent data are used, as determined by the year in the file name. Following conventions used SUPCRT92 (see Helgeson et al., 1978), the names of sanidine and microcline were changed to K-feldspar,high and K-feldspar,low (by using the same names in all data files, loading the optional SUPCRT92 data file updates these minerals rather than makes new ones).

• sympy.R is an R script that uses rSymPy to symbolically integrate Bermans's equations for heat capacity and volume to write experessions for enthalpy, entropy and Gibbs energy.

• The testing directory contains data files based on Berman and Aranovich (1996). These are used to demonstrate the addition of data from a user-supplied file (see Berman).

Files in cpetc contain experimental and calculated thermodynamic and environmental data:

• PM90.csv Heat capacities of four unfolded aqueous proteins taken from Privalov and Makhatadze, 1990. Temperature in °C is in the first column, and heat capacities of the proteins in J mol$^{-1}$ K$^{-1}$ in the remaining columns. See ionize.aa and the vignette anintro for examples that use this file.

• RH95.csv Heat capacity data for iron taken from Robie and Hemingway, 1995. Temperature in Kelvin is in the first column, heat capacity in J K$^{-1}$ mol$^{-1}$ in the second. See subcrt for an example that uses this file.

• SOJSH.csv Experimental equilibrium constants for the reaction NaCl(aq) = Na+ + Cl- as a function of temperature and pressure taken from Fig. 1 of Shock et al., 1992. See demo("NaCl") for an example that uses this file.

• HWM96_V.csv, HW97_Cp.csv Apparent molar volumes and heat capacities of CH₄, CO₂, H₂S, and NH₃ in dilute aqueous solutions reported by Hnědkovský et al., 1996 and Hnědkovský and Wood, 1997. Units are Kelvin, MPa, J/K/mol, and cm3/mol. See demo("AD"), EOSregress and the vignette eos-regress for examples that use these files.

• SC10_Rainbow.csv Values of temperature (°C), pH and logarithms of activity of CO₂, H₂, NH₄⁺, H₂S and CH₄ for mixing of seawater and hydrothermal fluid at Rainbow field (Mid-Atlantic Ridge), taken from Shock and Canovas, 2010. See the vignette anintro for an example that uses this file.

• SS98_Fig5a.csv, SS98_Fig5b.csv Values of logarithm of fugacity of O₂ and pH as a function of temperature for mixing of seawater and hydrothermal fluid, digitized from Figs. 5a and b of Shock and Schulte, 1998. See the vignette anintro for an example that uses this file.

• rubisco.csv UniProt IDs for Rubisco, ranges of optimal growth temperature of organisms, domain and name of organisms, and URL of reference for growth temperature, from Dick, 2014. See rank.affinity and the vignette anintro for examples that use this file.

• bluered.txt Blue - light grey - red color palette, computed using colorspace::diverge_hcl(1000, c = 100, l = c(50, 90), power = 1). This is used by ZC.col.

• AD03_Fig1?.csv Experimental data points digitized from Figure 1 of Akinfiev and Diamond, 2003, used in demo("AD").

• TKSS14_Fig2.csv Experimental data points digitized from Figure 2 of Tutolo et al., 2014, used in demo("aluminum").

• Mer75_Table4.csv Values of log(aK+/aH+) and log(aNa+/aH+) from Table 4 of Merino, 1975, used in demo("aluminum").

Files in protein contain protein sequences and amino acid compositions for proteins.

• rubisco.fasta Sequences of Rubisco obtained from UniProt (see Dick, 2014). See the vignette anintro for an example that uses this file.

• POLG.csv Amino acid compositions of a few proteins used for some tests and examples. These are various subunits of the Poliovirus type 1 polyprotein (POLG_POL1M in UniProt).

• TBD+05.csv lists genes with transcriptomic expression changes in carbon limitation stress response experiments in yeast (Tai et al., 2005).

• TBD+05_aa.csv has the amino acid compositions of proteins coded by those genes. The last two files are used in demo{"rank.affinity"}.

Files in taxonomy contain taxonomic data files:

• names.dmp and nodes.dmp are excerpts of NCBI taxonomy files (https://ftp.ncbi.nih.gov/pub/taxonomy/taxdump.tar.gz, accessed 2010-02-15). These files contain only the entries for Escherichia coli K-12, Saccharomyces cerevisiae, Homo sapiens, Pyrococcus furisosus and Methanocaldococcus jannaschii (taxids 83333, 4932, 9606, 186497, 243232) and the higher-ranking nodes (genus, family, etc.) in the respective lineages. See taxonomy for examples that use these files.

Files in adds contain additional thermodynamic data and group additivity definitions:

• BZA10.csv contains supplementary thermodynamic data taken from Bazarkina et al. (2010). The data can be added to the database in the current session using add.OBIGT. See add.OBIGT for an example that uses this file.

• OBIGT_check.csv contains the results of running check.OBIGT to check the internal consistency of entries in the default and optional datafiles.

• RH98_Table15.csv Group stoichiometries for high molecular weight crystalline and liquid organic compounds taken from Table 15 of Richard and Helgeson, 1998. The first three columns have the compound name, formula and physical state (‘⁠cr⁠’ or ‘⁠liq⁠’). The remaining columns have the numbers of each group in the compound; the names of the groups (columns) correspond to species in thermo$OBIGT. The compound named ‘⁠5a(H),14a(H)-cholestane⁠’ in the paper has been changed to ‘⁠5a(H),14b(H)-cholestane⁠’ here to match the group stoichiometry given in the table. See RH2OBIGT for a function that uses this file.

• SK95.csv contains thermodynamic data for alanate, glycinate, and their complexes with metals, taken from Amend and Helgeson (1997) and Shock and Koretsky (1995) as corrected in slop98.dat. These data are used in demo("copper") and demo("glycinate").

• LA19_test.csv contains thermodynamic data for dimethylamine and trimethylamine from LaRowe and Amend (2019) in energy units of both J and cal. This file is used in test-util.data.R) to check the messages produced by check.GHS and check.EOS.

References

Akinfiev, N. N. and Diamond, L. W. (2003) Thermodynamic description of aqueous nonelectrolytes at infinite dilution over a wide range of state parameters. Geochim. Cosmochim. Acta 67, 613–629. doi:10.1016/S0016-7037(02)01141-9

Amend, J. P. and Helgeson, H. C. (1997) Calculation of the standard molal thermodynamic properties of aqueous biomolecules at elevated temperatures and pressures. Part 1. L-α-amino acids. J. Chem. Soc., Faraday Trans. 93, 1927–1941. doi:10.1039/A608126F

Bazarkina, E. F., Zotov, A. V. and Akinfiev, N. N. (2010) Pressure-dependent stability of cadmium chloride complexes: Potentiometric measurements at 1–1000 bar and 25°C. Geol. Ore Deposits 52, 167–178. doi:10.1134/S1075701510020054

Berman, R. G. (1988) Internally-consistent thermodynamic data for minerals in the system Na₂O-K₂O-CaO-MgO-FeO-Fe₂O₃-Al₂O₃-SiO₂-TiO₂-H₂O-CO₂. J. Petrol. 29, 445-522. doi:10.1093/petrology/29.2.445

Berman, R. G. and Aranovich, L. Ya. (1996) Optimized standard state and solution properties of minerals. I. Model calibration for olivine, orthopyroxene, cordierite, garnet, and ilmenite in the system FeO-MgO-CaO-Al₂O₃-TiO₂-SiO₂. Contrib. Mineral. Petrol. 126, 1-24. doi:10.1007/s004100050233

Dick, J. M. (2014) Average oxidation state of carbon in proteins. J. R. Soc. Interface 11, 20131095. doi:10.1098/rsif.2013.1095

Gattiker, A., Michoud, K., Rivoire, C., Auchincloss, A. H., Coudert, E., Lima, T., Kersey, P., Pagni, M., Sigrist, C. J. A., Lachaize, C., Veuthey, A.-L., Gasteiger, E. and Bairoch, A. (2003) Automatic annotation of microbial proteomes in Swiss-Prot. Comput. Biol. Chem. 27, 49–58. doi:10.1016/S1476-9271(02)00094-4

Helgeson, H. C., Delany, J. M., Nesbitt, H. W. and Bird, D. K. (1978) Summary and critique of the thermodynamic properties of rock-forming minerals. Am. J. Sci. 278-A, 1–229. https://www.worldcat.org/oclc/13594862

Hnědkovský, L., Wood, R. H. and Majer, V. (1996) Volumes of aqueous solutions of CH₄, CO₂, H₂S, and NH₃ at temperatures from 298.15 K to 705 K and pressures to 35 MPa. J. Chem. Thermodyn. 28, 125–142. doi:10.1006/jcht.1996.0011

Hnědkovský, L. and Wood, R. H. (1997) Apparent molar heat capacities of aqueous solutions of CH₄, CO₂, H₂S, and NH₃ at temperatures from 304 K to 704 K at a pressure of 28 MPa. J. Chem. Thermodyn. 29, 731–747. doi:10.1006/jcht.1997.0192

Joint Genome Institute (2007) Bison Pool Environmental Genome. Protein sequence files downloaded from IMG/M (https://img.jgi.doe.gov/)

LaRowe, D. E. and Amend, J. P. (2019) The energetics of fermentation in natural settings. Geomicrobiol. J. 36, 492–505. doi:10.1080/01490451.2019.1573278

Merino, E. (1975) Diagenesis in teriary sandstones from Kettleman North Dome, California. II. Interstitial solutions: distribution of aqueous species at 100°C and chemical relation to diagenetic mineralogy. Geochim. Cosmochim. Acta 39, 1629–1645. doi:10.1016/0016-7037(75)90085-X

Privalov, P. L. and Makhatadze, G. I. (1990) Heat capacity of proteins. II. Partial molar heat capacity of the unfolded polypeptide chain of proteins: Protein unfolding effects. J. Mol. Biol. 213, 385–391. doi:10.1016/S0022-2836(05)80198-6

Richard, L. and Helgeson, H. C. (1998) Calculation of the thermodynamic properties at elevated temperatures and pressures of saturated and aromatic high molecular weight solid and liquid hydrocarbons in kerogen, bitumen, petroleum, and other organic matter of biogeochemical interest. Geochim. Cosmochim. Acta 62, 3591–3636. doi:10.1016/S0016-7037(97)00345-1

Robie, R. A. and Hemingway, B. S. (1995) Thermodynamic Properties of Minerals and Related Substances at 298.15 K and 1 Bar ($10^5$ Pascals) Pressure and at Higher Temperatures. U. S. Geol. Surv., Bull. 2131, 461 p. https://www.worldcat.org/oclc/32590140

Shock, E. and Canovas, P. (2010) The potential for abiotic organic synthesis and biosynthesis at seafloor hydrothermal systems. Geofluids 10, 161–192. doi:10.1111/j.1468-8123.2010.00277.x

Shock, E. L. and Koretsky, C. M. (1995) Metal-organic complexes in geochemical processes: Estimation of standard partial molal thermodynamic properties of aqueous complexes between metal cations and monovalent organic acid ligands at high pressures and temperatures. Geochim. Cosmochim. Acta 59, 1497–1532. doi:10.1016/0016-7037(95)00058-8

Shock, E. L., Oelkers, E. H., Johnson, J. W., Sverjensky, D. A. and Helgeson, H. C. (1992) Calculation of the thermodynamic properties of aqueous species at high pressures and temperatures: Effective electrostatic radii, dissociation constants and standard partial molal properties to 1000 °C and 5 kbar. J. Chem. Soc. Faraday Trans. 88, 803–826. doi:10.1039/FT9928800803

Shock, E. L. and Schulte, M. D. (1998) Organic synthesis during fluid mixing in hydrothermal systems. J. Geophys. Res. 103, 28513–28527. doi:10.1029/98JE02142

Tai, S. L., Boer, V. M., Daran-Lapujade, P., Walsh, M. C., de Winde, J. H., Daran, J.-M. and Pronk, J. T. (2005) Two-dimensional transcriptome analysis in chemostat cultures: Combinatorial effects of oxygen availability and macronutrient limitation in Saccharomyces cerevisiae. J. Biol. Chem. 280, 437–447. doi:10.1074/jbc.M410573200

Tutolo, B. M., Kong, X.-Z., Seyfried, W. E., Jr. and Saar, M. O. (2014) Internal consistency in aqueous geochemical data revisited: Applications to the aluminum system. Geochim. Cosmochim. Acta 133, 216–234. doi:10.1016/j.gca.2014.02.036

---

Properties of water from IAPWS-95

Description

Calculate thermodynamic properties of water following the IAPWS-95 formulation.

Usage

IAPWS95(property, T = 298.15, rho = 1000)

Arguments

property	character, name(s) of property(s) to calculate
T	numeric, temperature (K)
rho	numeric, density (kg m$^{-3}$)

Details

IAPWS95 provides an implementation of the IAPWS-95 formulation for properties (including pressure) calculated as a function of temperature and density.

The IAPWS95 function returns values of thermodynamic properties in specific units (per gram). The IAPWS-95 formulation follows the triple point convention used in engineering (values of internal energy and entropy are taken to be zero at the triple point).

For IAPWS95 the upper temperature limit of validity is 1000 °C, but extrapolation to much higher temperatures is possible (Wagner and Pruss, 2002). Valid pressures are from the greater of zero bar or the melting pressure at temperature to 10000 bar (with the provision for extrapolation to more extreme conditions). The function does not check these limits and will attempt calculations for any range of input parameters, but may return NA for properties that fail to be calculated at given temperatures and pressures and/or produce warnings or even errors when problems are encountered.

Value

A data frame the number of rows of which corresponds to the number of input temperature, pressure and/or density values.

References

Wagner, W. and Pruss, A. (2002) The IAPWS formulation 1995 for the thermodynamic properties of ordinary water substance for general and scientific use. J. Phys. Chem. Ref. Data 31, 387–535. doi:10.1063/1.1461829

See Also

util.water for properties along the saturation curve (WP02.auxiliary) and calculation of density from pressure and temperature (rho.IAPWS92). water.IAPWS95 is a wrapper around IAPWS95 and the utility functions, which converts the specific units to molar quantities, and is used in higher-level functions (water).

Examples

# Calculate pressure for given temperature and density
IAPWS95("P", T = 500, rho = 838.0235)

---

Search the thermodynamic database

Description

Search for species by name or formula, retrieve their thermodynamic properties and parameters, and add proteins to the thermodynamic database.

Usage

info(species = NULL, state = NULL, check.it=TRUE)

Arguments

species	character, names or formulas of species, or (for info only) numeric with same meaning as ispecies
state	character, physical states of the species
check.it	logical, check GHS and EOS parameters for self-consistency?

Details

info is the primary function used for querying the thermodynamic database (thermo()$OBIGT). It is often called recursively; first with a character value (or values) for species indicating the name(s) or formula(s) of the species of interest. The result of this call is a numeric value, which can be provided as an argument in a second call to info in order to retrieve a data frame of the thermodynamic properties of the species.

The text of species is searched in the names, chemical formulas, and abbreviations (in the ‘⁠abbrv⁠’ column) in the thermodynamic database. If the text of the species is matched, the index of that species is returned. If there are multiple matches for the species, and state is NULL, the index of first match is returned. The order of entries in the database is grouped by states in the order ‘⁠aq⁠’, ‘⁠cr⁠’, ‘⁠gas⁠’, ‘⁠liq⁠’. Therefore, for substances represented by both aqueous and gaseous species, the index of the aqueous species is returned, unless state is set to ‘⁠gas⁠’. Note that names (not formulas) of inorganic species, such as ‘⁠oxygen⁠’ and ‘⁠methane⁠’, are used only for the gas.

Names of species including an underscore character are indicative of proteins, e.g. ‘⁠LYSC_CHICK⁠’. If the name of a protein is provided to info and the composition of the protein can be found using pinfo, the thermodyamic properties and parameters of the nonionized protein (calculated using amino acid group additivity) are added to the thermodynamic database. Included in the return value, as for other species, is the index of the protein in the thermodynamic database or NA if the protein is not found. Names of proteins and other species can be mixed.

If no exact matches are found, info searches the database for similar names or formulas using agrep. If any of these are found, the results are summarized on the screen, but the function always returns NA in this case.

With a numeric argument, the rows of thermo()$OBIGT indicated by ispecies are returned, after removing any order-of-magnitude scaling factors (see thermo). If these species are all aqueous or are all not aqueous, the compounded column names used in thermo()$OBIGT are replaced with names appropriate for the corresponding equations of state. A missing value of one of the standard molal Gibbs energy (G) or enthalpy (H) of formation from the elements or entropy (S) is calculated from the other two, if available. If check.it is TRUE, several checks of self consistency among the thermodynamic properties and parameters are performed using check.GHS and check.EOS.

See Also

retrieve for searching species by element; check.OBIGT for checking self-consistency of each species.

Examples

## Summary of available data
info()

## Species information
# Search for something named (or whose formula is) "Fe"
si <- info("Fe")
# Use the number to get the full entry
info(si)
# Show data for the higher-temperature phases
info(si:(si+3))

## Dealing with states
# Order of precedence for names:
# aq > cr > gas > liq
info(c("ethanol", "adenosine"))  # aq, aq
# State argument overrides the default
info(c("ethanol", "adenosine"), state = c("gas", "cr"))
# Exceptions: gases have precedence for names of methane and inorganic gases
info(c("methane", "oxygen")) # gas, gas

# Formulas default to aqueous species, if available
i1 <- info(c("CH4", "CO2", "CS2", "MgO"))
info(i1)$state  # aq, aq, gas, cr
# State argument overrides the default
i2 <- info(c("CH4", "CO2", "MgO"), "gas") 
info(i2)$state  # gas, gas, NA

## Partial name or formula searches
info("ATP")
info("thiol")
info("MgC")
# Add an extra character to refine a search
# or to search using terms that have exact matches
info("MgC ")
info("acetate ")
info(" H2O")

---

Properties of ionization of proteins

Description

Calculate the charges of proteins and contributions of ionization to the thermodynamic properties of proteins.

Usage

ionize.aa(aa, property = "Z", T = 25, P = "Psat", pH = 7,
    ret.val = NULL, suppress.Cys = FALSE)

Arguments

aa	data frame, amino acid composition in the format of thermo()$protein
property	character, property to calculate
T	numeric, temperature in °C
P	numeric, pressure in bar, or ‘⁠Psat⁠’ for vapor pressure of H2O above 100 °C
pH	numeric, pH
ret.val	character, return the indicated value from intermediate calculations
suppress.Cys	logical, suppress (ignore) the ionization of the cysteine groups?

Details

The properties of ionization of proteins calculated by this function take account of the standard molal thermodynamic properties of ionizable amino acid sidechain groups and the terminal groups in proteins ([AABB]) and their equations of state parameters taken from Dick et al., 2006. The values of the ionization constants (pK) are calculated as a function of temperature, and the charges and the ionization contributions of other thermodynamic properties to the proteins are calculated additively, without consideration of electrostatic interactions, so they are best applied to the unfolded protein reference state.

For each amino acid composition in aa, the additive value of the property is calculated as a function of T, P and pH. property can be NULL to denote net charge, or if not NULL is one of the properties available in subcrt, or is ‘⁠A⁠’ to calculate the dimensionless chemical affinity (A/2.303RT) of the ionization reaction for the protein. If ret.val is one of ‘⁠pK⁠’, ‘⁠alpha⁠’, or ‘⁠aavals⁠’ it indicates to return the value of the ionization constant, degree of formation, or the values of the property for each ionizable group rather than taking their sums for the amino acid compositions in aa.

Value

The function returns a matrix (possibly with only one row or column) with number of rows corresponding to the longest of T, P or pH (values of any of these with shorter length are recycled) and a column for each of the amino acid compositions in aa.

References

Dick, J. M., LaRowe, D. E. and Helgeson, H. C. (2006) Temperature, pressure, and electrochemical constraints on protein speciation: Group additivity calculation of the standard molal thermodynamic properties of ionized unfolded proteins. Biogeosciences 3, 311–336. doi:10.5194/bg-3-311-2006

Makhatadze, G. I. and Privalov, P. L. (1990) Heat capacity of proteins. 1. Partial molar heat capacity of individual amino acid residues in aqueous solution: Hydration effect. J. Mol. Biol. 213, 375–384. doi:10.1016/S0022-2836(05)80197-4

Privalov, P. L. and Makhatadze, G. I. (1990) Heat capacity of proteins. II. Partial molar heat capacity of the unfolded polypeptide chain of proteins: Protein unfolding effects. J. Mol. Biol. 213, 385–391. doi:10.1016/S0022-2836(05)80198-6

See Also

pinfo, affinity

Examples

## Heat capacity of LYSC_CHICK as a function of T
pH <- c(5, 9, 3)
T <- seq(0, 100)
# Cp of non-ionized protein
Cp.nonion <- subcrt("LYSC_CHICK", T = T)$out[[1]]$Cp
plot(T, Cp.nonion, xlab = axis.label("T"), type = "l",
  ylab = axis.label("Cp"), ylim = c(20000, 35000))
# Cp of ionization and ionized protein
aa <- pinfo(pinfo("LYSC_CHICK"))
for(pH in c(5, 9, 3)) {
  Cp.ionized <- Cp.nonion + ionize.aa(aa, "Cp", T = T, pH = pH)[, 1]
  lines(T, Cp.ionized, lty = 2)
  text(80, Cp.ionized[70], paste("pH =", pH) )
}
# Makhatadze and Privalov's group contributions
T <- c(5, 25, 50, 75, 100, 125)
points(T, MP90.cp("LYSC_CHICK", T))
# Privalov and Makhatadze's experimental values
e <- read.csv(system.file("extdata/cpetc/PM90.csv", package = "CHNOSZ"))
points(e$T, e$LYSC_CHICK, pch = 16)
legend("bottomright", pch = c(16, 1, NA, NA), lty = c(NA, NA, 1, 2),
  legend = c("PM90 experiment", "MP90 groups", 
  "DLH06 groups no ion", "DLH06 groups ionized"))
title("Heat capacity of unfolded LYSC_CHICK")

---

Fit thermodynamic parameters to formation constants (log β)

Description

Fit thermodynamic parameters to experimental formation constants for an aqueous species and add the parameters to OBIGT.

Usage

logB.to.OBIGT(logB, species, coeffs, T, P, npar = 3,
    optimize.omega = FALSE, tolerance = 0.05)

Arguments

logB	Values of log β
species	Species in the formation reaction (the formed species must be last)
coeffs	Coefficients of the formation reaction
T	Temperature (units of T.units)
P	Temperature (‘⁠Psat⁠’ or numerical values with units of P.units)
npar	numeric, number of parameters to use in fitting
optimize.omega	logical, optimize the omega parameter (relevant for charged species)?
tolerance	Tolerance for checking equivalence of input and calculated log β values

Details

This function updates the OBIGT thermodynamic database with parameters fit to formation constants of aqueous species as a function of temperature. The logB argument should have decimal logarithm of formation constants for an aqueous complex (log β). The formation reaction is defined by species and coeffs. All species in the formation reaction must be present in OBIGT except for the last species, which is the newly formed species.

The function works as follows. First, the properties of the known species are combined with log β to calculate the standard Gibbs energy (G[T]) of the formed species at each value of T and P. Then, lm is used to fit one or more of the thermodynamic parameters G, S, c1, c2, and omega to the values of G[T]. The first two of these parameters are standard-state values at 25 °C and 1 bar, and the last three are parameters in the revised HKF equations (e.g. Eq. B25 of Shock et al., 1992). The fitted parameters for the formed species are then added to OBIGT. Finally, all.equal is used to test for approximate equivalence of the input values of log β and calculated equilibrium constants; if the mean absolute difference exceeds tolerance, an error occurs.

To avoid overfitting, only the first three parameters (G, S, and c1) are used by default. More parameters (up to 5) or fewer (down to 1) can be selected by changing npar. Volumetric parameters (a₁ to a₄) in the HKF equations currently aren't included, so the resulting fits are valid only at the input pressure values.

Independent of npar, the number of parameters used in the fit is not more than the number of data values (n). If n is less than 5, then the values of the unused parameters are set to 0 for addition to OBIGT. For instance, a single value of log β would be fit only with G, implying that computed values of G[T] have no temperature dependence.

The value of ω is a constant in the revised HKF equations for uncharged species, but for charged species, it is a function of T and P as described by the “g function” (Shock et al., 1992). An optimization step is available to refine the value of omega at 25 °C and 1 bar for charged species taking its temperature dependence into account for the fitting. However, representative datasets are not well represented by variable omega (see first example), so this step is skipped by default. Furthermore, logB.to.OBIGT by default also sets the z parameter in OBIGT to 0; this enforces a constant ω for charged species in calculations with subcrt (note that this is a parameter for the HKF equations and does not affect reaction balancing). Set optimize.omega to TRUE to override the defaults and activate variable ω for charged species; this takes effect only if npar == 5 and n > 5.

Value

The species index of the new species in OBIGT.

References

Migdisov, Art. A., Zezin, D. and Williams-Jones, A. E. (2011) An experimental study of Cobalt (II) complexation in Cl^- and H₂S-bearing hydrothermal solutions. Geochim. Cosmochim. Acta 75, 4065–4079. doi:10.1016/j.gca.2011.05.003

Mei, Y., Sherman, D. M., Liu, W., Etschmann, B., Testemale, D. and Brugger, J. (2015) Zinc complexation in chloride-rich hydrothermal fluids (25–600 °C): A thermodynamic model derived from ab initio molecular dynamics. Geochim. Cosmochim. Acta 150, 264–284. doi:10.1016/j.gca.2014.09.023

Shock, E. L., Oelkers, E. H., Johnson, J. W., Sverjensky, D. A. and Helgeson, H. C. (1992) Calculation of the thermodynamic properties of aqueous species at high pressures and temperatures: Effective electrostatic radii, dissociation constants and standard partial molal properties to 1000 °C and 5 kbar. J. Chem. Soc. Faraday Trans. 88, 803–826. doi:10.1039/FT9928800803

See Also

logB.to.OBIGT calls mod.OBIGT with zap = TRUE to clear the parameters of a formed species if it already exists in the OBIGT database. If preexisting parameters (e.g. volumetric HKF coefficients) weren't cleared, they would interfere with the fitting done here, which uses only selected parameters.

Examples

## CoHS+ from Migdisov et al. (2011)
logB <- c(6.24, 6.02, 5.84, 5.97, 6.52)
T <- c(120, 150, 200, 250, 300)
P <- "Psat"
species <- c("Co+2", "HS-", "CoHS+")
coeffs <- c(-1, -1, 1)
opar <- par(mfrow = c(2, 1))
for(o.o in c(TRUE, FALSE)) { 
  # Fit the parameters with or without variable omega
  inew <- logB.to.OBIGT(logB, species, coeffs, T, P, npar = 5, optimize.omega = o.o)
  # Print the new database entry
  info(inew)
  # Plot experimental logB
  plot(T, logB, "n", c(100, 320), c(5.8, 6.8),
    xlab = axis.label("T"), ylab = quote(log~beta))
  points(T, logB, pch = 19, cex = 2)
  # Plot calculated values
  Tfit <- seq(100, 320, 10)
  sres <- subcrt(species, coeffs, T = Tfit)
  lines(sres$out$T, sres$out$logK, col = 4)
  title(describe.reaction(sres$reaction))
  legend <- c("Migdisov et al. (2011)",paste0("logB.to.OBIGT(optimize.omega = ",o.o,")"))
  legend("top", legend, pch = c(19, NA), lty = c(0, 1), col = c(1, 4),
    pt.cex = 2, bg = "#FFFFFFB0")
}
par(opar)
# NB. Optimizing omega leads to unphysical oscillations in the logK (first plot)

## ZnCl+ from Mei et al. (2015)
# Values for 5000 bar calculated with the modified Ryzhenko-Bryzgalin (RB) model
logB <- c(-1.93,-1.16,-0.38,0.45,1.15,1.76,2.30,2.80,3.26,3.70,4.12,4.53,4.92)
species <- c("Zn+2", "Cl-", "ZnCl+")
coeffs <- c(-1, -1, 1)
T <- c(25, 60, 100, 150, 200, 250, 300, 350, 400, 450, 500, 550, 600)
P <- 5000
# Note: ZnCl+ is in the default OBIGT database;
# logB.to.OBIGT() "zaps" the previous parameters before adding the fitted ones
inew <- logB.to.OBIGT(logB, species, coeffs, T, P, npar = 5)
# Plot RB and logB.to.OBIGT values
plot(T, logB, xlab = axis.label("T"), ylab = axis.label("logB"), pch = 19, cex = 2)
Tfit <- seq(25, 600, 25)
sres <- subcrt(species, coeffs, T = Tfit, P = P)
lines(sres$out$T, sres$out$logK, col = 4)
title(describe.reaction(sres$reaction), line = 3)
title("5000 bar", font.main = 1, line = 1)
legend <- c("Mei et al. (2015)", "logB.to.OBIGT()")
legend("topleft", legend, pch = c(19, NA), lty = c(0, 1), col = c(1, 4), pt.cex = 2)

---

Parse chemical formulas

Description

Count the elements and charges in a chemical formula.

Usage

makeup(formula, multiplier = 1, sum = FALSE, count.zero = FALSE)
  count.elements(formula)

Arguments

formula	character, a chemical formula
multiplier	numeric, multiplier for the elemental counts in each formula
sum	logical, add together the elemental counts in all formulas?
count.zero	logical, include zero counts for elements?

Details

makeup parses a chemical formula expressed in string notation, returning the numbers of each element in the formula. The formula may carry a charge, indicated by a + or - sign, possibly followed by a magnitude, after the uncharged part of the formula. The formula may have multiple subformulas enclosed in parentheses (but the parentheses may not be nested), each one optionally followed by a numeric coefficient. The formula may have one suffixed subformula, separated by ‘⁠*⁠’ or ‘⁠:⁠’, optionally preceded by a numeric coefficient. All numbers may contain a decimal point.

Each subformula (or the entire formula without subformulas) should be a simple formula. A simple formula, processed by count.elements, must adhere to the following pattern: it starts with an elemental symbol; all elemental symbols start with an uppercase letter, and are followed by another elemental symbol, a number (possibly fractional, possibly signed), or nothing (the end of the formula). Any sequence of one uppercase letter followed by zero or more lowercase letters is recognized as an elemental symbol. makeup will issue a warning for elemental symbols that are not present in thermo$element.

makeup can handle numeric and length > 1 values for the formula argument. If the argument is numeric, it identifies row number(s) in thermo()$OBIGT from which to take the formulas of species. If formula has length > 1, the function returns a list containing the elemental counts in each of the formulas. If count.zero is TRUE, the elemental counts for each formula include zeros to indicate elements that are only present in any of the other formulas.

The multiplier argument must have either length = 1 or length equal to the number of formulas. The elemental count in each formula is multiplied by the respective value. If sum is true, the elemental counts in all formulas (after any multiplying) are summed together to yield a single bulk formula.

Value

A numeric vector with names refering to each of the elemental symbols in the formula. If more than one formula is provided, a list of numeric vectors is returned, unless sum is TRUE.

See Also

mass, entropy, basis, i2A

Examples

# Elemental composition of a simple compound
makeup("CO2")     # 1 carbon, 2 oxygen
# Formula of lawsonite, with a parenthetical part and a suffix
makeup("CaAl2Si2O7(OH)2*H2O")
# Fractional coefficients are OK
reddiv10 <- makeup("C10.6N1.6P0.1")
10*reddiv10 # 106, 16, 1 (Redfield ratio)

# The coefficient for charge is a number with a *preceding* sign
# e.g., ferric iron, with a charge of +3 is expressed as
makeup("Fe+3")
# Transcribing the formula the way it appears in many 
# publications produces a likely unintended result:
# 3 iron atoms and a charge of +1
makeup("Fe3+")

# These all represent a single negative charge, i.e., electron
makeup("-1")
makeup("Z-1+0")
makeup("Z0-1")   # the "old" formula for the electron in thermo()$OBIGT
makeup("(Z-1)")  # the current formula in thermo()$OBIGT

# Hypothetical compounds with negative numbers of elements
makeup("C-4(O-2)")   # -4 carbon, -2 oxygen
makeup("C-4O-2")     # -4 carbon,  1 oxygen, -2 charge
makeup("C-4O-2-2")   # -4 carbon, -2 oxygen, -2 charge

# The 'sum' argument can be used to check mass and charge
# balance in a chemical reaction
formula <- c("H2O", "H+", "(Z-1)", "O2")
mf <- makeup(formula, c(-1, 2, 2, 0.5), sum = TRUE)
all(mf == 0)  # TRUE

---