# NOT RUN {
# 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
stopifnot(all(abs(DGrxn) < 100))
### compare mineral stabilities in the Berman and Helgeson datasets
### on a T - log(K+/H+) diagram, after Sverjensky et al., 1991
### (doi:10.1016/0016-7037(91)90157-Z)
## set up the system: basis species
basis(c("K+", "Al+3", "quartz", "H2O", "O2", "H+"))
# use pH = 0 so that aK+ = aK+/aH+
basis("pH", 0)
# load the species
species(c("K-feldspar", "muscovite", "kaolinite",
"pyrophyllite", "andalusite"), "cr")
## start with the data from Helgeson et al., 1978
add.obigt("SUPCRT92")
# calculate affinities in aK+ - temperature space
# exceed.Tr: enable calculations above stated temperature limit of pyrophyllite
res <- 400
a <- affinity(`K+` = c(0, 5, res), T = c(200, 650, res), P = 1000, exceed.Ttr = TRUE)
# make base plot with colors and no lines
diagram(a, xlab = ratlab("K+", molality = TRUE), lty = 0, fill = "terrain")
# add the lines, extending into the low-density region (exceed.rhomin = TRUE)
a <- affinity(`K+` = c(0, 5, res), T = c(200, 650, res), P = 1000,
exceed.Ttr = TRUE, exceed.rhomin = TRUE)
diagram(a, add = TRUE, names = NULL, col = "red", lwd = 1.5)
# the list of references:
ref1 <- thermo.refs(species()$ispecies)$key
## now use the (default) data from Berman, 1988
# this resets the thermodynamic database
# without affecting the basis and species settings
obigt()
# we can check that we have Berman's quartz
# and not coesite or some other phase of SiO2
iSiO2 <- rownames(basis()) == "SiO2"
stopifnot(info(basis()$ispecies[iSiO2])$name == "quartz")
# Berman's dataset doesn't have the upper temperature limits,
# so we don't need exceed.Ttr here
a <- affinity(`K+` = c(0, 5, res), T = c(200, 650, res), P = 1000, exceed.rhomin = TRUE)
diagram(a, add = TRUE, names = NULL, col = "blue", lwd = 1.5)
# the list of references:
ref2 <- thermo.refs(species()$ispecies)$key
ref2 <- paste(ref2, collapse = ", ")
# add legend and title
legend("top", "low-density region", text.font = 3, bty = "n")
legend("topleft", describe.property(c("P", "IS"), c(1000, 1)), bty = "n")
legend("left", c(ref1, ref2),
lty = c(1, 1), lwd = 1.5, col = c(2, 4), bty = "n")
title(main = syslab(c("K2O", "Al2O3", "SiO2", "H2O", "HCl")), line = 1.8)
title(main = "Helgeson and Berman minerals, after Sverjensky et al., 1991",
line = 0.3, font.main = 1)
# cleanup for next example
reset()
# 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
E.units("J")
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
stopifnot((Gex_BA96 - Gex_Ber88)[1] > 0)
# the curves cross at about 725 deg C (BA96 Fig. 8)
# (actually, in our calculation they cross closer to 800 deg C)
stopifnot(T[which.min(abs(Gex_BA96 - Gex_Ber88))] == 800)
# reset the database (thermo()$opt$E.units, thermo()$obigt, and thermo()$opt$Berman)
reset()
# }
Run the code above in your browser using DataLab