CHNOSZ (version 1.3.2)

add.obigt: Functions to Work with the Thermodynamic Database

Description

Add or modify species in the thermodynamic database.

Usage

add.obigt(file, species = NULL, force = TRUE, E.units = "cal")
  mod.obigt(...)
  today()

Arguments

file

character, path to a file

species

character, names of species to load from file

force

logical, force replacement of already existing species?

E.units

character, units of energy, cal or J

...

character or numeric, properties of species to modify in the thermodynamic database

Value

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

Details

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.

file is first matched against the names of files in the extdata/OBIGT directory packaged with CHNOSZ. In this case, the filename suffixes are removed, so DEW_aq, organic_aq, and organic_cr are valid names. The function also accepts single matches with the state suffix dropped (DEW but not organic). 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).

Given the default setting of E.units, the function does not perform any unit conversions. If E.units is set to J, then the thermodynamic parameters are converted from units of Joules to calories, as used in the CHNOSZ database.

When adding 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. These changes are lost if you reload the database by calling data(thermo) or if you quit the R session without saving it. The name of the species to add or change must be supplied as the first argument of ... or as a named argument (named name). 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 results if the formula is not valid (i.e. can not be parsed bymakeup). Additional arguments refer to the name of the property(s) to be updated and are matched to any part of compound column names in thermo$obigt, such as z or T in z.T. Unless state is specified as one of the properties, its value is taken from thermo()$opt$state. When adding species, properties that are not specified become NA, except for state, which takes a default value from thermo()$opt$state, and z.T, which for aqueous species is set to the charge calculated from the chemical formula (otherwise, NA charge for newly added species would trigger the AkDi model). The values provided should be in the units specifed in the documentation for the thermo data object, including any order-of-magnitude scaling factors.

today returns the current date in the format adopted for thermo()$obigt (inherited from SUPCRT-format data files) e.g. 13.May.12 for May 13, 2012.

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<U+00B0>C. Geology of Ore Deposits 52, 167--178. https://doi.org/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. https://doi.org/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. https://doi.org/10.1016/0016-7037(89)90341-4

Stefnsson, A. (2001) Dissolution of primary minerals of basalt in natural waters. I. Calculation of mineral solubilities from 0C to 350C. Chem. Geol. 172, 225--250. https://doi.org/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 <U+00B0>C and 5 kbar. Geochim. Cosmochim. Acta 61, 1359--1412. https://doi.org/10.1016/S0016-7037(97)00009-4

See Also

thermo, util.data, mod.buffer

Examples

# NOT RUN {
## modify an existing species (example only)
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 a species
iCl2O <- mod.obigt("Cl2O", G=20970)
info(iCl2O)
# add a species with a name that is different from the formula
mod.obigt("buckminsterfullerene", formula="C60", state="cr", date=today())
# retrieve the species data (thermodynamic properties in this toy example are NA)
info(info("C60"))
# 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
})
# reset default database
obigt()
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))
}))

# another use of add.obigt()
# compare Delta G of AABB = UPBB + H2O
# (Figure 9 of Kitadai, 2014)
E.units("J")
# default database has values from Kitadai, 2014
Kit14 <- subcrt(c("[AABB]", "[UPBB]", "H2O"), c(-1, 1, 1), T = seq(0, 300, 10))
# optional file OldAA has superseded values of [UPBB] from Dick et al., 2006
add.obigt("OldAA")
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")
# reset database *and* settings (units)
reset()

# Another use of add.obigt(): calculate Delta G of
# H4SiO4 = SiO2 + 2H2O using different 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(-100, 600))
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", legend = as.expression(stxt), bty = "n")
legend("topright", c("Stef\u00e1nsson, 2001", "Shock et al., 1989",
  "Apps and Spycher, 2004"), lty = c(0, 1, 2), bty = "n")
reset()
# }