# 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=as.character(Sys.Date()))
# 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()
# }
Run the code above in your browser using DataLab