[CHNOSZ-commits] r398 - in pkg/CHNOSZ: . R demo inst man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Feb 18 12:59:48 CET 2019
Author: jedick
Date: 2019-02-18 12:59:47 +0100 (Mon, 18 Feb 2019)
New Revision: 398
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/R/util.expression.R
pkg/CHNOSZ/demo/NaCl.R
pkg/CHNOSZ/inst/NEWS
pkg/CHNOSZ/man/examples.Rd
pkg/CHNOSZ/man/util.expression.Rd
Log:
demo/NaCl.R: improve labels and documentation of discontinuities
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2019-02-17 21:43:50 UTC (rev 397)
+++ pkg/CHNOSZ/DESCRIPTION 2019-02-18 11:59:47 UTC (rev 398)
@@ -1,6 +1,6 @@
Date: 2019-02-18
Package: CHNOSZ
-Version: 1.2.0-4
+Version: 1.2.0-5
Title: Thermodynamic Calculations and Diagrams for Geochemistry
Authors at R: c(
person("Jeffrey", "Dick", , "j3ffdick at gmail.com", role = c("aut", "cre"),
Modified: pkg/CHNOSZ/R/util.expression.R
===================================================================
--- pkg/CHNOSZ/R/util.expression.R 2019-02-17 21:43:50 UTC (rev 397)
+++ pkg/CHNOSZ/R/util.expression.R 2019-02-18 11:59:47 UTC (rev 398)
@@ -255,7 +255,7 @@
if(i %in% iname) species <- reaction$name[i]
else {
# should the chemical formula have a state?
- if(identical(states,"all")) species <- expr.species(reaction$formula[i], state=reaction$state[i], use.state=TRUE)
+ if(identical(states,"all") | i %in% states) species <- expr.species(reaction$formula[i], state=reaction$state[i], use.state=TRUE)
else species <- expr.species(reaction$formula[i], state=reaction$state[i])
}
# get the absolute value of the reaction coefficient
@@ -278,7 +278,8 @@
}
}
# put an equals sign between reactants and products
- desc <- substitute(a==b, list(a=reactexpr, b=prodexpr))
+ # change this to unicode for the reaction double-arrow 20190218 \u21cc
+ desc <- substitute(a ~ "\u21cc" ~ b, list(a=reactexpr, b=prodexpr))
return(desc)
}
Modified: pkg/CHNOSZ/demo/NaCl.R
===================================================================
--- pkg/CHNOSZ/demo/NaCl.R 2019-02-17 21:43:50 UTC (rev 397)
+++ pkg/CHNOSZ/demo/NaCl.R 2019-02-18 11:59:47 UTC (rev 398)
@@ -4,40 +4,59 @@
## 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 degrees C and 5 kbar.
## J. Chem. Soc. Faraday Trans. 88, 803-826. https://doi.org/10.1039/FT9928800803 )
-species <- c("NaCl", "Na+", "Cl-")
-coeffs <- c(-1, 1, 1)
+
+## uncomment these lines to make the plot with the g-function disabled
+#mod.obigt("Cl-", z=0)
+#mod.obigt("Na+", z=0)
+
# start a new plot and show the experimental logK
thermo.plot.new(xlim=c(0, 1000), ylim=c(-5.5, 1),
xlab=axis.label("T"), ylab=axis.label("logK"))
expt <- read.csv(system.file("extdata/cpetc/SOJSH.csv",
package="CHNOSZ"), as.is=TRUE)
points(expt$T,expt$logK, pch=expt$pch)
+
# we'll be at 9 distinct pressure conditions, including Psat
-P <- c(list("Psat"), as.list(seq(500, 4000, by=500)))
+# Psat is repeated to show "not considered" region
+# (T >= 355 degC; Fig. 6 of Shock et al., 1992)
+P <- c(list("Psat", "Psat"), as.list(seq(500, 4000, by=500)))
# for each of those what's the range of temperature
T <- list()
-# T > 350 degC at Psat is possibly inappropriate; see "Warning" of subcrt.Rd
-T[[1]] <- seq(0, 370, 5)
-T[[2]] <- seq(265, 465, 5)
-T[[3]] <- seq(285, 760, 5)
-T[[4]] <- seq(395, 920, 5)
-T[[5]] <- T[[6]] <- T[[7]] <- T[[8]] <- T[[9]] <- seq(400, 1000, 5)
+T[[1]] <- seq(0, 354, 1)
+T[[2]] <- seq(354, 370, 1)
+T[[3]] <- seq(265, 465, 1)
+T[[4]] <- seq(285, 760, 1)
+T[[5]] <- seq(395, 920, 1)
+T[[6]] <- T[[7]] <- T[[8]] <- T[[9]] <- T[[10]] <- seq(400, 1000, 1)
+
# calculate and plot the logK
+species <- c("NaCl", "Na+", "Cl-")
+coeffs <- c(-1, 1, 1)
logK <- numeric()
for(i in 1:length(T)) {
s <- subcrt(species, coeffs, T=T[[i]], P=P[[i]])
- lines(s$out$T, s$out$logK)
- # keep the calculated values for each experimental condition
+ if(i==2) lty <- 3 else lty <- 1
+ lines(s$out$T, s$out$logK, lty=lty)
+ # keep the calculated values for each experimental condition (excluding Psat)
iexpt <- which(P[[i]]==expt$P)
Texpt <- expt$T[iexpt]
- logK <- c(logK, splinefun(s$out$T, s$out$logK)(Texpt))
+ if(i > 2) logK <- c(logK, splinefun(s$out$T, s$out$logK)(Texpt))
}
+
+# add title, labels, and legends
+title(describe.reaction(s$reaction, states = 1))
+text(150, -0.1, quote(italic(P)[sat]), cex=1.2)
+text(462, -4, "500 bar")
+text(620, -4.3, "1000 bar")
+text(796, -4.3, "1500 bar")
+text(813, -1.4, "4000 bar")
legend("bottomleft",pch=unique(expt$pch),
- legend=c(unique(expt$source),tail(expt$source,1)))
-mtitle(c(describe.reaction(s$reaction), expression(italic(P)[sat]~"or 500-4000 bar, after Shock et al., 1992")))
-# where do we diverge most from experiment?
-imaxdiff <- which.max(abs(logK - expt$logK))
-stopifnot(all.equal(c("Psat", 347.7),
- as.character(expt[imaxdiff,1:2])))
-# what's our average divergence?
+ legend=c(unique(expt$source),tail(expt$source,1)), bty="n")
+#mtitle(c(describe.reaction(s$reaction), expression(italic(P)[sat]~"and 500-4000 bar")))
+l1 <- quote("Revised HKF model with " * italic(g) * " function (Shock et al., 1992)")
+l2 <- "Non-recommended region (Shock et al., 1992, Fig. 6)"
+legend("topright", as.expression(c(l1, l2)), lty=c(1, 3), bty="n")
+
+# test for average divergence (excluding Psat)
+expt <- expt[!expt$P %in% "Psat", ]
stopifnot(mean(abs(logK - expt$logK)) < 0.09)
Modified: pkg/CHNOSZ/inst/NEWS
===================================================================
--- pkg/CHNOSZ/inst/NEWS 2019-02-17 21:43:50 UTC (rev 397)
+++ pkg/CHNOSZ/inst/NEWS 2019-02-18 11:59:47 UTC (rev 398)
@@ -1,4 +1,4 @@
-CHANGES IN CHNOSZ 1.2.0-4 (2019-02-18)
+CHANGES IN CHNOSZ 1.2.0-5 (2019-02-18)
--------------------------------------
CRAN COMPLIANCE
@@ -20,8 +20,24 @@
THERMODYNAMIC DATA
-- Move SiO2(aq) and H2AsO3- from SLOP98.csv to SUPCRT92.csv.
+- Revert to using SiO2(aq) from SUPCRT92 (i.e. Shock et al., 1989) in
+ the default database.
+- Move SiO2(aq) from Apps and Spycher, 2004 and recalculated HSiO3- to
+ new optional data file, OBIGT/AS04.csv.
+
+- Move H2AsO3- from SLOP98.csv to SUPCRT92.csv.
+
+DOCUMENTATION
+
+- In demo/NaCl.R, indicate region not considered by Shock et al., 1992
+ in developing the revised HKF model and presence of the associated
+ discontinuities in SUPCRT92.
+
+OTHER CHANGES
+
+- In describe.reaction(), change equals sign to reaction double arrow.
+
CHANGES IN CHNOSZ 1.2.0 (2019-02-09)
------------------------------------
Modified: pkg/CHNOSZ/man/examples.Rd
===================================================================
--- pkg/CHNOSZ/man/examples.Rd 2019-02-17 21:43:50 UTC (rev 397)
+++ pkg/CHNOSZ/man/examples.Rd 2019-02-18 11:59:47 UTC (rev 398)
@@ -73,12 +73,14 @@
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 \eqn{\log a_{\mathrm{H_{2}}_{\left(aq\right)}}=-11+3/\left(40\times T\left(^{\circ}C\right)\right)}{logaH2 = -11 + 3/40 * T(degC)}; this equation comes from a model of relative stabilities of proteins in a hot-spring environment (Dick and Shock, 2011).
-In the \samp{NaCl} demo, the \logK lines calculated at \Psat and P=500 bar show discontinuities at 355 \degC.
-Although not realistic, this behavior is consistent with the output of \acronym{SUPCRT92} (Johnson et al., 1992) at 500 bar.
-This is probably due to a transition between different regimes for the properties of water as coded in SUPCRT's \code{H2O92D.F}, which is used by CHNOSZ.
-(Note that SUPCRT does not output thermodynamic properties above 350 \degC at \Psat; see Warning in \code{\link{subcrt}}.)
}
+\section{Warning}{
+The discontinuities apparent in the plot made by the \code{NaCl} demo illustrate limitations of the "\emph{g} function" for charged species in the revised HKF model (the 355 \degC boundary of region II in Figure 6 of Shock et al., 1992).
+Note that \acronym{SUPCRT92} (Johnson et al., 1992) gives similar output at 500 bar.
+However, \acronym{SUPCRT} does not output thermodynamic properties above 350 \degC at \Psat; see Warning in \code{\link{subcrt}}.
+}
+
\references{
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\degC and pressures of 1-2000 bar. \emph{Geochem. Int.} \bold{39}, 990--1006. \url{http://pleiades.online/cgi-perl/search.pl/?type=abstract&name=geochem&number=10&year=1&page=990}
Modified: pkg/CHNOSZ/man/util.expression.Rd
===================================================================
--- pkg/CHNOSZ/man/util.expression.Rd 2019-02-17 21:43:50 UTC (rev 397)
+++ pkg/CHNOSZ/man/util.expression.Rd 2019-02-18 11:59:47 UTC (rev 398)
@@ -50,7 +50,7 @@
\item{ret.val}{logical, return only the value with the units?}
\item{reaction}{data frame, definition of reaction}
\item{iname}{numeric, show names instead of formulas for these species}
- \item{states}{character, if \samp{all}, show states for all species}
+ \item{states}{character, if \samp{all}, show states for all species; numeric, which species to show states for}
\item{system}{character, thermodynamic components}
\item{dash}{character to use for dash between components}
\item{ion}{character, an ion}
@@ -112,7 +112,12 @@
\code{describe.property} makes an expression summarizing the properties supplied in \code{property}, along with their \code{value}s. The expressions returned by both functions consist of a property, an equals sign, and a value (with units where appropriate); the expressions have a length equal to the number of property/value pairs. If \code{oneline} is TRUE, the property/value pairs are combined into a single line, separated by commas. The number of digits shown after the decimal point in the values is controlled by \code{digits}. If \code{ret.val} is TRUE, only the values and their units are returned; this is useful for labeling plots with values of temperature.
- \code{describe.reaction} makes an expression summarizing a chemical reaction. The \code{reaction} data frame can be generated using \code{\link{subcrt}}. Based on the sign of their reaction coefficients, species are placed on the reactant (left) or product (right) side of the reaction, where the species with their coefficients are separated by plus signs; the two sides of the reaction are separated by an equals sign. Coefficients equal to 1 are not shown. Chemical formulas of species include a designation of physical state if \code{states} is \samp{all}. Names of species (as provided in \code{reaction}) are shown instead of chemical formulas for the species identified by \code{iname}.
+\code{describe.reaction} makes an expression summarizing a chemical reaction.
+The \code{reaction} data frame can be generated using \code{\link{subcrt}}.
+Based on the sign of their reaction coefficients, species are placed on the reactant (left) or product (right) side of the reaction, where the species with their coefficients are separated by plus signs; the two sides of the reaction are separated by a reaction double arrow (Unicode U+21CC).
+Coefficients equal to 1 are not shown.
+Chemical formulas of species include the physical state if \code{states} is \samp{all}, or a numeric value indicating which species to label with the state.
+Names of species (as provided in \code{reaction}) are shown instead of chemical formulas for the species identified by \code{iname}.
\code{syslab} formats the given thermodynamic components (using \code{expr.species}) and adds intervening en dashes.
More information about the CHNOSZ-commits
mailing list