[CHNOSZ-commits] r347 - in pkg/CHNOSZ: . R demo inst man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Nov 8 07:50:23 CET 2018
Author: jedick
Date: 2018-11-08 07:50:22 +0100 (Thu, 08 Nov 2018)
New Revision: 347
Added:
pkg/CHNOSZ/demo/QtzMsKfs.R
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/R/examples.R
pkg/CHNOSZ/R/util.expression.R
pkg/CHNOSZ/demo/00Index
pkg/CHNOSZ/demo/gold.R
pkg/CHNOSZ/inst/NEWS
pkg/CHNOSZ/man/berman.Rd
pkg/CHNOSZ/man/examples.Rd
pkg/CHNOSZ/man/util.expression.Rd
Log:
add demo/QtzMsKfs.R
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2018-11-07 15:50:34 UTC (rev 346)
+++ pkg/CHNOSZ/DESCRIPTION 2018-11-08 06:50:22 UTC (rev 347)
@@ -1,6 +1,6 @@
-Date: 2018-11-07
+Date: 2018-11-08
Package: CHNOSZ
-Version: 1.1.3-54
+Version: 1.1.3-55
Title: Thermodynamic Calculations and Diagrams for Geo(bio)chemistry
Authors at R: c(
person("Jeffrey", "Dick", , "j3ffdick at gmail.com", role = c("aut", "cre"),
Modified: pkg/CHNOSZ/R/examples.R
===================================================================
--- pkg/CHNOSZ/R/examples.R 2018-11-07 15:50:34 UTC (rev 346)
+++ pkg/CHNOSZ/R/examples.R 2018-11-08 06:50:22 UTC (rev 347)
@@ -28,7 +28,7 @@
demos <- function(which=c("sources", "protein.equil", "affinity", "NaCl", "density",
"ORP", "revisit", "findit", "ionize", "buffer", "protbuff", "yeastgfp", "mosaic",
- "copper", "solubility", "gold", "wjd", "bugstab", "Shh", "activity_ratios",
+ "copper", "solubility", "gold", "QtzMsKfs", "wjd", "bugstab", "Shh", "activity_ratios",
"adenine", "DEW", "lambda", "TCA", "go-IU", "bison"), save.png=FALSE) {
# run one or more demos from CHNOSZ with ask=FALSE, and return the value of the last one
for(i in 1:length(which)) {
Modified: pkg/CHNOSZ/R/util.expression.R
===================================================================
--- pkg/CHNOSZ/R/util.expression.R 2018-11-07 15:50:34 UTC (rev 346)
+++ pkg/CHNOSZ/R/util.expression.R 2018-11-08 06:50:22 UTC (rev 347)
@@ -82,6 +82,7 @@
propchar <- s2c(property)
expr <- ""
# some special cases
+ if(is.na(property)) return("")
if(property=="logK") return(quote(log~italic(K)))
# grepl here b/c diagram() uses "loga.equil" and "loga.basis"
if(grepl("loga", property)) {
@@ -223,7 +224,8 @@
propexpr <- valexpr <- character()
for(i in 1:length(property)) {
propexpr <- c(propexpr, expr.property(property[i]))
- if(value[i]=="Psat") thisvalexpr <- quote(italic(P)[sat])
+ if(is.na(value[i])) thisvalexpr <- ""
+ else if(value[i]=="Psat") thisvalexpr <- quote(italic(P)[sat])
else {
thisvalue <- format(round(as.numeric(value[i]), digits), nsmall=digits)
thisunits <- expr.units(property[i])
@@ -236,7 +238,8 @@
# write an equals sign between the property and value
desc <- character()
for(i in seq_along(propexpr)) {
- thisdesc <- substitute(a==b, list(a=propexpr[[i]], b=valexpr[[i]]))
+ if(is.na(value[i])) thisdesc <- propexpr[[i]]
+ else thisdesc <- substitute(a==b, list(a=propexpr[[i]], b=valexpr[[i]]))
if(oneline) {
# put all the property/value equations on one line, separated by commas
if(i==1) desc <- substitute(a, list(a=thisdesc))
@@ -284,7 +287,7 @@
}
# make formatted text for activity ratio 20170217
-ratlab <- function(ion="K+") {
+ratlab <- function(ion="K+", use.molality=FALSE) {
# the charge
Z <- makeup(ion)["Z"]
# the text for the exponent on aH+
@@ -292,9 +295,11 @@
# the expression for the ion and H+
expr.ion <- expr.species(ion)
expr.H <- expr.species("H+")
+ # with use.molality, change a to m
+ a <- ifelse(use.molality, "m", "a")
# the final expression
- if(exp.H=="1") substitute(log~(italic(a)[expr.ion] / italic(a)[expr.H]), list(expr.ion=expr.ion, expr.H=expr.H))
- else substitute(log~(italic(a)[expr.ion] / italic(a)[expr.H]^exp.H), list(expr.ion=expr.ion, expr.H=expr.H, exp.H=exp.H))
+ if(exp.H=="1") substitute(log~(italic(a)[expr.ion] / italic(a)[expr.H]), list(a=a, expr.ion=expr.ion, expr.H=expr.H))
+ else substitute(log~(italic(a)[expr.ion] / italic(a)[expr.H]^exp.H), list(a=a, expr.ion=expr.ion, expr.H=expr.H, exp.H=exp.H))
}
# make formatted text for thermodynamic system 20170217
Modified: pkg/CHNOSZ/demo/00Index
===================================================================
--- pkg/CHNOSZ/demo/00Index 2018-11-07 15:50:34 UTC (rev 346)
+++ pkg/CHNOSZ/demo/00Index 2018-11-08 06:50:22 UTC (rev 347)
@@ -14,6 +14,7 @@
copper Another example of mosaic(): complexation of copper with glycine species
solubility Solubility of calcite and CO2(gas) as a function of pH
gold Solubility of gold
+QtzMsKfs Helgeson and Berman minerals and calculations with molality
wjd Gibbs energy minimization: prebiological atmospheres and cell periphery of yeast
dehydration log K of dehydration reactions; SVG file contains tooltips and links
bugstab Formation potential of microbial proteins in colorectal cancer
Added: pkg/CHNOSZ/demo/QtzMsKfs.R
===================================================================
--- pkg/CHNOSZ/demo/QtzMsKfs.R (rev 0)
+++ pkg/CHNOSZ/demo/QtzMsKfs.R 2018-11-08 06:50:22 UTC (rev 347)
@@ -0,0 +1,69 @@
+# CHNOSZ/demo/QtzMsKfs.R
+# T - log(K+/H+) diagram, after Sverjensky et al., 1991
+# (doi:10.1016/0016-7037(91)90157-Z)
+# 20171009 diagram added to berman.Rd
+# 20181108 moved to demo/QtzMsKfs.R; add molality calculations
+
+# this demo compares diagrams made using the Berman and Helgeson datasets,
+# and shows the use of nonideal calculations to set molalities in the basis species
+
+## 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")
+## the "b_gamma" equation gets closer to the published diagram than "B-dot"
+thermo$opt$nonideal <<- "bgamma"
+
+## start with the data from Helgeson et al., 1978
+add.obigt("SUPCRT92")
+# calculate affinities in aK+ - temperature space
+# exceed.Tr: we go above stated temperature limit of pyrophyllite
+# (this is above its stability field on the diagram, so pyrophyllite doesn't appear in this region,
+# but its properties are needed needed to calculate relative stabilities of all minerals)
+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+", use.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", lty = 2, lwd = 1.5)
+# calculate and plot the lines for 1 molal chloride
+a <- affinity(`K+` = c(0, 5, res), T = c(200, 650, res), P = 1000, exceed.Ttr = TRUE, exceed.rhomin = TRUE, IS = 1)
+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
+data(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", lty = 2, lwd = 1.5)
+a <- affinity(`K+` = c(0, 5, res), T = c(200, 650, res), P = 1000, exceed.rhomin = TRUE, IS = 1)
+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 experimental points for 1000 bar (Table 1 of Sverjensky et al., 1991)
+expt.T <- c(300, 400, 500, 550, # KFs-Ms-Qtz
+ 400, 450, 500, 550, # Ms-And-Qtz
+ 300, 350, # Ms-P-Qtz
+ 300, 600) # Kaol-Ms-Qtz, KFs-And-Qtz
+expt.KH <- c(3.50, 2.75, 1.95, 1.40, 1.60, 1.57, 1.47, 1.38, 1.94, 1.80, 1.90, 0.63)
+points(expt.KH, expt.T, pch = 19, cex = 1.2)
+# add legend and title
+legend("top", "low-density region", text.font = 3, bty = "n")
+legend("top", describe.property(c(NA, NA, "P", "IS"), c(NA, NA, 1000, 1)), bty = "n")
+legend("left", c(ref1, ref2, "ion molality", "ion activity", "experiments"),
+ lty = c(1, 1, 1, 2, 0), lwd = 1.5, col = c(2, 4, 1, 1, 1), pch = c(NA, NA, NA, NA, 19), 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)
Modified: pkg/CHNOSZ/demo/gold.R
===================================================================
--- pkg/CHNOSZ/demo/gold.R 2018-11-07 15:50:34 UTC (rev 346)
+++ pkg/CHNOSZ/demo/gold.R 2018-11-08 06:50:22 UTC (rev 347)
@@ -132,15 +132,11 @@
basis("H2S", "PPM")
# apply QMK buffer for pH
basis("H+", "QMK")
+ basis("K+", log10(0.5))
# calculate solution composition for 2 mol/kg NaCl
NaCl <- NaCl(T = seq(150, 550, 10), P = 1000, m_tot=2)
- a_Cl <- NaCl$m_Cl * NaCl$gam_Cl
- # using this ionic strength, calculate the activity of K+
- # assuming complete dissociation of 0.5 mol/kg KCl
- gam_K <- 10^subcrt("K+", T = seq(150, 550, 10), P = 1000, IS=NaCl$IS)$out$`K+`$loggam
- a_K <- 0.5 * gam_K
# calculate affinity and solubility
- a <- affinity(T = seq(150, 550, 10), `Cl-` = log10(a_Cl), `K+` = log10(a_K), P = 1000, IS = NaCl$IS)
+ a <- affinity(T = seq(150, 550, 10), `Cl-` = log10(NaCl$m_Cl), P = 1000, IS = NaCl$IS)
s <- solubility(a)
# make diagram and show total log molality
diagram(s, ylim = c(-10, -4), col = col, lwd = 2, lty = 1)
@@ -166,15 +162,11 @@
basis("O2", "HM")
# apply QMK buffer for pH
basis("H+", "QMK")
+ basis("K+", log10(0.5))
# calculate solution composition for 2 mol/kg NaCl
- NaCl <- NaCl(T = seq(150, 550, 10), P = 1000, m_tot=2)
- a_Cl <- NaCl$m_Cl * NaCl$gam_Cl
- # using this ionic strength, calculate the activity of K+
- # assuming complete dissociation of 0.5 mol/kg KCl
- gam_K <- 10^subcrt("K+", T = seq(150, 550, 10), P = 1000, IS=NaCl$IS)$out$`K+`$loggam
- a_K <- 0.5 * gam_K
+ NaCl <- NaCl(T = seq(150, 550, 10), P = 1000, m_tot=1)
# calculate affinity and solubility
- a <- affinity(T = seq(150, 550, 10), `Cl-` = log10(a_Cl), `K+` = log10(a_K), P = 1000, IS = NaCl$IS)
+ a <- affinity(T = seq(150, 550, 10), `Cl-` = log10(NaCl$m_Cl), P = 1000, IS = NaCl$IS)
s <- solubility(a)
# make diagram and show total log molality
diagram(s, ylim = c(-10, -2), col = col, lwd = 2, lty = 1)
Modified: pkg/CHNOSZ/inst/NEWS
===================================================================
--- pkg/CHNOSZ/inst/NEWS 2018-11-07 15:50:34 UTC (rev 346)
+++ pkg/CHNOSZ/inst/NEWS 2018-11-08 06:50:22 UTC (rev 347)
@@ -1,19 +1,22 @@
-CHANGES IN CHNOSZ 1.1.3-54 (2018-11-07)
+CHANGES IN CHNOSZ 1.1.3-55 (2018-11-08)
---------------------------------------
MAJOR BUG FIXED
-- Previously, the virtual conversion of activities to molalities was
- done incorrectly; the Gibbs energies of basis species and formed
- species were transformed in the same direction. They should have
- opposite transformations. This has been fixed by adding the 'is.basis'
- argument to nonideal(), which is used by affinity() to transform the
- Gibbs energies in the right direction.
+- Previously, with calculations using nonideal(), the Gibbs energies of
+ basis species and formed species were transformed in the same
+ direction, causing an incorrct conversion of activities to molalities.
+ This has been fixed by adding the 'is.basis' argument to nonideal(),
+ which is used by affinity() to transform the Gibbs energies in opposite
+ directions for basis species and formed species.
- For more information, see the new section of ?nonideal: 'is.basis and
the CHNOSZ workflow'.
- Two new demos depend on the corrected behavior: gold.R and QtzMsKfs.R.
+ The latter demo also provides a comparison of the superseded Helgeson
+ (SUPCRT92) and newer Berman datasets for minerals. See below for more
+ information on demo/gold.R.
NEW FEATURES
Modified: pkg/CHNOSZ/man/berman.Rd
===================================================================
--- pkg/CHNOSZ/man/berman.Rd 2018-11-07 15:50:34 UTC (rev 346)
+++ pkg/CHNOSZ/man/berman.Rd 2018-11-08 06:50:22 UTC (rev 347)
@@ -100,36 +100,6 @@
a <- affinity(T=c(200, 1700, 200), P=c(0, 50000, 200))
diagram(a)
-## a longer example, comparing diagrams made using the
-## Berman and Helgeson datasets, after Sverjensky et al., 1991
-res <- 200
-# using the Helgeson data
-add.obigt("SUPCRT92")
-# set up 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")
-# calculate affinities in aK+ - temperature space
-a <- affinity(`K+`=c(0, 5, res), T=c(200, 650, res), P=1000)
-# note that we go just past the quartz transition,
-# but it has no effect on the diagram
-diagram(a, xlab=ratlab("K+"))
-# now using the Berman data: this resets the thermodynamic database
-# without affecting the basis and species settings
-data(OBIGT)
-# it might be good to check that we have Berman's quartz
-# and not coesite or some other SiO2 phase
-info(basis()$ispecies[3])
-a <- affinity(`K+`=c(0, 5, res), T=c(200, 650, res), P=1000)
-diagram(a, add=TRUE, names="", col="blue", lwd=2)
-legend("topleft", lty=c(1, 1, NA), lwd=c(1, 2, 0), col=c("black", "blue", ""),
- legend=c("Helgeson et al., 1978 (unadjusted)",
- "Berman, 1988", " (adjusted by Sverjensky et al., 1991)"), bty="n")
-title(main="Comparison of Helgeson and Berman datasets at 1000 bar")
-
## Getting data from a user-supplied file
## Ol-Opx exchange equilibrium, after Berman and Aranovich, 1996
E.units("J")
@@ -161,8 +131,6 @@
Berman, R. G. (2007) winTWQ (version 2.3): A software package for performing internally-consistent thermobarometric calculations. \emph{Open File} \bold{5462}, Geological Survey of Canada, 41 p. \url{https://doi.org/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. \emph{Am. J. Sci.} \bold{278-A}, 1--229. \url{http://www.worldcat.org/oclc/13594862}
-
-Sverjensky, D. A., Hemley, J. J. and D'Angelo, W. M. (1991) Thermodynamic assessment of hydrothermal alkali feldspar-mica-aluminosilicate equilibria. \emph{Geochim. Cosmochim. Acta} \bold{55}, 989-1004. \url{https://doi.org/10.1016/0016-7037(91)90157-Z}
}
\concept{Thermodynamic calculations}
Modified: pkg/CHNOSZ/man/examples.Rd
===================================================================
--- pkg/CHNOSZ/man/examples.Rd 2018-11-07 15:50:34 UTC (rev 346)
+++ pkg/CHNOSZ/man/examples.Rd 2018-11-08 06:50:22 UTC (rev 347)
@@ -16,8 +16,8 @@
demos(which = c("sources", "protein.equil", "affinity", "NaCl",
"density", "ORP", "revisit", "findit", "ionize", "buffer",
"protbuff", "yeastgfp", "mosaic", "copper", "solubility",
- "gold", "wjd", "bugstab", "Shh", "activity_ratios", "adenine",
- "DEW", "lambda", "TCA", "go-IU", "bison"),
+ "gold", "QtzMsKfs", "wjd", "bugstab", "Shh", "activity_ratios",
+ "adenine", "DEW", "lambda", "TCA", "go-IU", "bison"),
save.png=FALSE)
}
@@ -47,6 +47,7 @@
\code{copper} \tab * Another example of \code{\link{mosaic}}: complexation of Cu with glycine (Aksu and Doyle, 2001) \cr
\code{solubility} \tab * Solubility of calcite (cf. Manning et al., 2013) and \CO2 (cf. Stumm and Morgan, 1996) \cr
\code{gold} \tab * Solubility of gold (Akinfiev and Zotov; 2001; Stef{\aacute}nsson and Seward, 2004; Williams-Jones et al., 2009) \cr
+ \code{QtzMsKfs} \tab * Helgeson and Berman minerals and calculations with molality (Sverjensky et al., 1991) \cr
\code{wjd} \tab * \eqn{G}{G} minimization: prebiological atmospheres (Dayhoff et al., 1964) and cell periphery of yeast \cr
\code{dehydration} \tab * \logK of dehydration reactions; SVG file contains tooltips and links \cr
\code{bugstab} \tab * Formation potential of microbial proteins in colorectal cancer (Dick, 2016) \cr
@@ -122,6 +123,8 @@
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 \degC. \emph{Geochim. Cosmochim. Acta} \bold{129}, 125--145. \url{https://doi.org/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. \emph{Geochim. Cosmochim. Acta} \bold{55}, 989-1004. \url{https://doi.org/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. \emph{Nat. Geosci.} \bold{7}, 909--913. \url{https://doi.org/10.1038/ngeo2291}
Williams-Jones, A. E., Bowell, R. J. and Migdisov, A. A. (2009) Gold in solution. \emph{Elements} \bold{5}, 281--287. \url{https://doi.org/10.2113/gselements.5.5.281}
Modified: pkg/CHNOSZ/man/util.expression.Rd
===================================================================
--- pkg/CHNOSZ/man/util.expression.Rd 2018-11-07 15:50:34 UTC (rev 346)
+++ pkg/CHNOSZ/man/util.expression.Rd 2018-11-08 06:50:22 UTC (rev 347)
@@ -26,7 +26,7 @@
ret.val = FALSE)
describe.reaction(reaction, iname = numeric(), states = NULL)
syslab(system = c("K2O", "Al2O3", "SiO2", "H2O"), dash="\u2013")
- ratlab(ion = "K+")
+ ratlab(ion = "K+", use.molality = FALSE)
}
\arguments{
More information about the CHNOSZ-commits
mailing list