[CHNOSZ-commits] r286 - in pkg/CHNOSZ: . R data inst man tests/testthat
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Nov 10 12:43:46 CET 2017
Author: jedick
Date: 2017-11-10 12:43:44 +0100 (Fri, 10 Nov 2017)
New Revision: 286
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/R/info.R
pkg/CHNOSZ/R/subcrt.R
pkg/CHNOSZ/data/opt.csv
pkg/CHNOSZ/inst/NEWS
pkg/CHNOSZ/man/data.Rd
pkg/CHNOSZ/man/info.Rd
pkg/CHNOSZ/man/subcrt.Rd
pkg/CHNOSZ/tests/testthat/test-subcrt.R
Log:
add thermo$opt$Berman and some checks to subcrt()
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2017-11-10 04:02:50 UTC (rev 285)
+++ pkg/CHNOSZ/DESCRIPTION 2017-11-10 11:43:44 UTC (rev 286)
@@ -1,6 +1,6 @@
Date: 2017-11-10
Package: CHNOSZ
-Version: 1.1.0-83
+Version: 1.1.0-84
Title: Thermodynamic Calculations for Geobiochemistry
Author: Jeffrey Dick
Maintainer: Jeffrey Dick <j3ffdick at gmail.com>
Modified: pkg/CHNOSZ/R/info.R
===================================================================
--- pkg/CHNOSZ/R/info.R 2017-11-10 04:02:50 UTC (rev 285)
+++ pkg/CHNOSZ/R/info.R 2017-11-10 11:43:44 UTC (rev 286)
@@ -131,7 +131,11 @@
# if a single name matches, use that one (useful for distinguishing pseudo-H4SiO4 and H4SiO4) 20171020
matches.name <- matches.species & thermo$obigt$name==species
if(sum(matches.name)==1) ispecies.out <- which(matches.name)
- else ispecies.out <- ispecies[1] # otherwise, return only the first species that matches
+ else {
+ # prefer the Berman minerals?
+ if(thermo$opt$Berman & "cr_Berman" %in% thermo$obigt$state[ispecies]) ispecies.out <- ispecies[thermo$obigt$state[ispecies]=="cr_Berman"]
+ else ispecies.out <- ispecies[1] # otherwise, return only the first species that matches
+ }
# let user know if there is more than one state for this species
mystate <- thermo$obigt$state[ispecies.out]
ispecies.other <- ispecies[!ispecies %in% ispecies.out]
Modified: pkg/CHNOSZ/R/subcrt.R
===================================================================
--- pkg/CHNOSZ/R/subcrt.R 2017-11-10 04:02:50 UTC (rev 285)
+++ pkg/CHNOSZ/R/subcrt.R 2017-11-10 11:43:44 UTC (rev 286)
@@ -37,7 +37,6 @@
}
do.reaction <- FALSE
- #if(!missing(coeff) & coeff!=1) do.reaction <- TRUE
if(!missing(coeff)) do.reaction <- TRUE
# species and states are made the same length
@@ -67,6 +66,9 @@
if(convert) P <- envert(P,'bar')
}
+ # warn for too high temperatures for Psat 20171110
+ if(identical(P, "Psat") & any(T > 647.067)) warning("attempting calculation at P = 'Psat' for some T > Tcritical; set P = 1 (or higher)")
+
# gridding?
do.grid <- FALSE
if(!is.null(grid)) if(!is.logical(grid)) do.grid <- TRUE
@@ -145,6 +147,15 @@
ton <- thermo$obigt$name
tos <- thermo$obigt$state
+ # warn if we're running a reaction with both Berman and Helgeson minerals 20171110
+ if(do.reaction) {
+ ref1 <- thermo$obigt$ref1
+ ref2 <- thermo$obigt$ref2
+ hasHelgeson <- any(grepl("HDNB78", ref1[sinfo])) | any(grepl("HDNB78", ref2[sinfo]))
+ hasBerman <- any(tos[sinfo]=="cr_Berman")
+ if(hasHelgeson & hasBerman) warning("the reaction has minerals from both the Helgeson and Berman datasets; data may not be internally consistent")
+ }
+
# stop if species not found
noname <- is.na(sinfo)
if(TRUE %in% noname)
@@ -174,9 +185,9 @@
}
# where we keep info about the species involved
- reaction <- data.frame( coeff=coeff.new,name=ton[inpho],
- formula = thermo$obigt$formula[inpho],state=tos[inpho],
- ispecies=inpho, stringsAsFactors=FALSE)
+ reaction <- data.frame(coeff = coeff.new, name = ton[inpho],
+ formula = thermo$obigt$formula[inpho], state = tos[inpho],
+ ispecies = inpho, stringsAsFactors = FALSE)
# make the rownames readable ... but they have to be unique
if(length(unique(inpho))==length(inpho)) rownames(reaction) <- as.character(inpho)
@@ -347,7 +358,7 @@
# we don't warn here about the transition
warn.above <- FALSE
}
- if(all(is.na(Ttr))) next
+ if(any(is.na(Ttr))) next
if(!all(Ttr==0) & any(T >= Ttr)) {
status.Ttr <- "(extrapolating G)"
if(!exceed.Ttr) {
@@ -425,9 +436,11 @@
for(j in 1:nrow(G)) {
ps <- which.min(as.numeric(G[j,]))
if(length(ps)==0) {
- # minimum not found: NAs have crept in (like something wrong with Psat?)
- # (or no non-NA value of G to begin with, e.g. aegerine)
- ps <- 1
+ # minimum not found (we have NAs)
+ # - no non-NA value of G to begin with, e.g. aegerine) --> probably should use lowest-T phase
+ #ps <- 1
+ # - above temperature limit for the highest-T phase (subcrt.Rd skarn example) --> use highest-T phase 20171110
+ ps <- ncol(G)
if(exceed.Ttr) warning('subcrt: stable phase for ',reaction$name[iphases[ps]],' at T-P point ',j,
' undetermined (using ',reaction$state[iphases[ps]],')',call.=FALSE)
}
Modified: pkg/CHNOSZ/data/opt.csv
===================================================================
--- pkg/CHNOSZ/data/opt.csv 2017-11-10 04:02:50 UTC (rev 285)
+++ pkg/CHNOSZ/data/opt.csv 2017-11-10 11:43:44 UTC (rev 286)
@@ -1,2 +1,2 @@
-cutoff,E.units,T.units,P.units,state,water,G.tol,Cp.tol,V.tol,varP,IAPWS.sat,paramin,ideal.H,ideal.e,nonideal
-1E-10,cal,C,bar,aq,SUPCRT92,100,1,1,FALSE,liquid,1000,TRUE,TRUE,Helgeson
+cutoff,E.units,T.units,P.units,state,water,G.tol,Cp.tol,V.tol,varP,IAPWS.sat,paramin,ideal.H,ideal.e,nonideal,Berman
+1E-10,cal,C,bar,aq,SUPCRT92,100,1,1,FALSE,liquid,1000,TRUE,TRUE,Helgeson,FALSE
Modified: pkg/CHNOSZ/inst/NEWS
===================================================================
--- pkg/CHNOSZ/inst/NEWS 2017-11-10 04:02:50 UTC (rev 285)
+++ pkg/CHNOSZ/inst/NEWS 2017-11-10 11:43:44 UTC (rev 286)
@@ -1,4 +1,4 @@
-CHANGES IN CHNOSZ 1.1.0-80 (2017-11-08)
+CHANGES IN CHNOSZ 1.1.0-84 (2017-11-10)
---------------------------------------
MAJOR CHANGES:
@@ -146,6 +146,9 @@
- Add 'tol' argument to equil.reaction() (convergence tolerance for
uniroot()).
+- Add thermo$opt$Berman to signal info(), subcrt(), etc. to prefer
+ minerals in Berman dataset over Helgeson dataset (default: FALSE).
+
REFACTORING AND CLEANUP:
- In hkf() and cgl(), combine 'ghs' and 'eos' arguments into single
Modified: pkg/CHNOSZ/man/data.Rd
===================================================================
--- pkg/CHNOSZ/man/data.Rd 2017-11-10 04:02:50 UTC (rev 285)
+++ pkg/CHNOSZ/man/data.Rd 2017-11-10 11:43:44 UTC (rev 286)
@@ -10,21 +10,23 @@
\alias{buffers}
\title{Thermodynamic Database and System Settings}
\description{
-Run \code{data(thermo)} to initialize all of the data used in CHNOSZ.
+Run \code{data(thermo)} to initialize or reset all of the data used in CHNOSZ.
This includes the computational settings, thermodynamic database, and system settings (chemical species).
The system settings are changed using \code{\link{basis}} and \code{\link{species}}.
-To restore the default system settings (no species loaded), use \code{basis("")}.
+To restore the default system settings (no species loaded), run \code{basis("")}.
The thermodynamic database is changed using \code{\link{add.obigt}} and \code{\link{mod.obigt}}.
-To restore the default database, use \code{data(OBIGT)}.
+To restore the default database, run \code{data(OBIGT)}.
-The computational settings are changed using \code{\link{water}}, \code{\link{P.units}}, \code{\link{T.units}}, \code{\link{E.units}} or by direct manipulation of \code{thermo$opt}, as well as some other commands (e.g. \code{\link{mod.buffer}}).
-To restore the default computational settings, thermodynamic database, and system settings, use \code{data(thermo)}.
+The computational settings are changed using \code{\link{water}}, \code{\link{P.units}}, \code{\link{T.units}}, \code{\link{E.units}}, and some other commands (e.g. \code{\link{mod.buffer}}).
-In an interactive session, use \code{data(thermo)} to restore the thermodynamic database and computational and system settings to their initial state.
-Or, use \code{basis("")} to clear the system settings (basis and species), or \code{data(OBIGT)} to restore the default database.
+Some settings can only be changed by direct manipulation of \code{thermo$opt}.
+In an interactive session, this should be done using the super-assignment operator (e.g. \code{thermo$opt$Berman <<- TRUE}) so that the \code{thermo} object is not copied to the global environment.
+(Doing so would cause problems, as many functions are designed to access the \code{thermo} object in the \code{CHNOSZ} environment.)
+To restore the default computational settings, thermodynamic database, and system settings, run \code{data(thermo)}.
+
The data files provided with CHNOSZ are in the \code{data} and \code{extdata/OBIGT} directories of the package.
The \code{*.csv} files in these directories are used to build the \code{thermo} data object in an environment named \code{CHNOSZ}.
The structure of the \code{thermo} object is described below.
@@ -37,8 +39,6 @@
\format{
- The items in the \code{thermo} object are documented below.
-
\itemize{
\item \code{thermo$opt}
@@ -57,7 +57,9 @@
\code{IAPWS.sat} \tab character \tab State of water for saturation properties [\samp{liquid}] (see \code{\link{util.water}})\cr
\code{paramin} \tab integer \tab Minimum number of calculations to launch parallel processes [1000] (see \code{\link{palply}}) \cr
\code{ideal.H} \tab logical \tab Should \code{\link{nonideal}} ignore the proton? [\code{TRUE}] \cr
- \code{ideal.e} \tab logical \tab Should \code{\link{nonideal}} ignore the electron? [\code{TRUE}]
+ \code{ideal.e} \tab logical \tab Should \code{\link{nonideal}} ignore the electron? [\code{TRUE}] \cr
+ \code{nonideal} \tab character \tab Method for \code{\link{nonideal}} [\code{Helgeson}] \cr
+ \code{Berman} \tab logical \tab Should \code{\link{info}} preferentially return Berman minerals? [\code{FALSE}]
}
\item \code{thermo$element}
@@ -217,7 +219,7 @@
} % end of format
\seealso{
-Other data files, supporting the examples and vignettes, are documented separately at \code{\link{extdata}}.
+Other data files, including those supporting the examples and vignettes, are documented separately at \code{\link{extdata}}.
}
\examples{
Modified: pkg/CHNOSZ/man/info.Rd
===================================================================
--- pkg/CHNOSZ/man/info.Rd 2017-11-10 04:02:50 UTC (rev 285)
+++ pkg/CHNOSZ/man/info.Rd 2017-11-10 11:43:44 UTC (rev 286)
@@ -26,6 +26,9 @@
If there are multiple matches for the \code{species}, and \code{state} is NULL, the index of first match is returned.
The order of entries in \code{thermo$obigt} is grouped by states in the order \samp{aq}, \samp{cr}, \samp{gas}, \samp{liq}, so for species in both aqueous and gaseous states the index of the aqueous species is returned, unless \code{state} is set to \samp{gas}.
+Note that entries for minerals using the \code{\link{berman}} equations are placed after the dataset derived from Helgeson et al., 1978, so the latter data are used by default in case of duplicates.
+This behavior can be changed by setting \code{\link{thermo}$opt$Berman} to TRUE; alternatively, individual minerals in the Berman dataset can be identifed by setting \code{state} to \samp{cr_Berman}.
+
Names of species including an underscore character are indicative of proteins, e.g. \samp{LYSC_CHICK}.
If the name of a protein is provided to \code{info} and the composition of the protein can be found using \code{\link{protein}}, the thermodyamic properties and parameters of the nonionized protein (calculated using amino acid group additivity) are added to the thermodynamic database.
Included in the return value, as for other species, is the index of the protein in the thermodynamic database or \code{NA} if the protein is not found. Names of proteins and other species can be mixed.
Modified: pkg/CHNOSZ/man/subcrt.Rd
===================================================================
--- pkg/CHNOSZ/man/subcrt.Rd 2017-11-10 04:02:50 UTC (rev 285)
+++ pkg/CHNOSZ/man/subcrt.Rd 2017-11-10 11:43:44 UTC (rev 286)
@@ -98,8 +98,6 @@
Changes were made in CHNOSZ 0.9-9 to improve the calculation of the \eqn{g}{g}- and \eqn{f}{f}-functions (Shock et al., 1992) for the partial derivatives of the omega parameter which are used by the \code{\link{hkf}} function, so thermodynamic properties of charged aqueous species over the temperature range 0 to 1000 \degC are now mostly consistent with \acronym{SUPCRT92}.
Note, however, that while \acronym{SUPCRT92} limits calculations at Psat to 350 \degC (\dQuote{beyond range of applicability of aqueous species equations}), there is no corresponding limit in place in \code{subcrt} (or \code{\link{hkf}}), so that inapplicable calculations will be performed up to the critical temperature (373.917 \degC) without any warning.
It is probably for this reason that there is a noticeable jump in the value of \logK at Psat shown in the one of the examples (\code{demo("NaCl")}).
-
-A known bug is misidentification of the stable polymorph of some minerals at high temperature; an example of this bug is shown in the \code{\link{stopifnot}} at the end of the skarn example below.
}
\value{
@@ -167,10 +165,13 @@
subcrt(c("pyrite", "pyrrhotite", "H2O", "H2S", "O2"), c(-1, 1, -1, 1, 0.5))
## these produce messages about problems with the calculation
+# above the T, P limits for the H2O equations of state
+subcrt("alanine", T=c(2250, 2251), P=c(30000, 30001), grid="T")
# Psat is not defined above the critical point
+\dontrun{
+## (also gives a warning)
subcrt("alanine", T=seq(0, 5000, by=1000))
-# above the T, P limits for the H2O equations of state
-subcrt("alanine", T=c(2250, 2251), P=c(30000, 30001), grid="T")
+}
## minerals with phase transitions
# compare calculated values of heat capacity of iron with
@@ -206,9 +207,9 @@
P=P, T=T, grid="P")
# The results are not identical to SUPCRT92, as CHNOSZ has updated
# parameters for species e.g. Cu+ from Shock et al., 1997.
-# FIXME: the last 1 should be a 3
+# Check the calculated phase transitions for chalcopyrite
stopifnot(all.equal(s$polymorphs$chalcopyrite[1:11],
- c(1, 1, 1, 1, 1, 1, 2, 3, 3, 3, 1)))
+ c(1, 1, 1, 1, 1, 1, 2, 3, 3, 3, 3)))
## Standard Gibbs energy of reactions with HCN and
## formaldehyde, after Schulte and Shock, 1995 Fig. 1
Modified: pkg/CHNOSZ/tests/testthat/test-subcrt.R
===================================================================
--- pkg/CHNOSZ/tests/testthat/test-subcrt.R 2017-11-10 04:02:50 UTC (rev 285)
+++ pkg/CHNOSZ/tests/testthat/test-subcrt.R 2017-11-10 11:43:44 UTC (rev 286)
@@ -156,6 +156,10 @@
#expect_equal(as.numeric(a$values[[1]]), c(0, 0))
})
+test_that("warning is produced for reaction with Helgeson and Berman minerals", {
+ expect_warning(subcrt(c("quartz", "quartz"), c(-1, 1), c("cr", "cr_Berman")), "data may not be internally consistent")
+})
+
# references
# Amend, J. P. and Shock, E. L. (2001)
More information about the CHNOSZ-commits
mailing list