[CHNOSZ-commits] r744 - in pkg/CHNOSZ: . R inst inst/extdata/OBIGT inst/tinytest man vignettes
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Sep 20 04:33:43 CEST 2022
Author: jedick
Date: 2022-09-20 04:33:43 +0200 (Tue, 20 Sep 2022)
New Revision: 744
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/R/add.OBIGT.R
pkg/CHNOSZ/R/cgl.R
pkg/CHNOSZ/R/hkf.R
pkg/CHNOSZ/R/subcrt.R
pkg/CHNOSZ/R/util.data.R
pkg/CHNOSZ/inst/NEWS.Rd
pkg/CHNOSZ/inst/extdata/OBIGT/AD.csv
pkg/CHNOSZ/inst/tinytest/test-AD.R
pkg/CHNOSZ/inst/tinytest/test-eos.R
pkg/CHNOSZ/inst/tinytest/test-info.R
pkg/CHNOSZ/man/add.OBIGT.Rd
pkg/CHNOSZ/man/subcrt.Rd
pkg/CHNOSZ/man/thermo.Rd
pkg/CHNOSZ/vignettes/multi-metal.Rmd
Log:
Start documenting 'model' and add it to reaction data frame in subcrt()
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2022-09-19 09:32:39 UTC (rev 743)
+++ pkg/CHNOSZ/DESCRIPTION 2022-09-20 02:33:43 UTC (rev 744)
@@ -1,6 +1,6 @@
-Date: 2022-09-19
+Date: 2022-09-20
Package: CHNOSZ
-Version: 1.9.9-35
+Version: 1.9.9-36
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/add.OBIGT.R
===================================================================
--- pkg/CHNOSZ/R/add.OBIGT.R 2022-09-19 09:32:39 UTC (rev 743)
+++ pkg/CHNOSZ/R/add.OBIGT.R 2022-09-20 02:33:43 UTC (rev 744)
@@ -58,6 +58,8 @@
newrows$E_units <- thermo$opt$E.units
# Fill in the columns
newrows[, icol] <- args[inew, ]
+ # Guess model from state 20220919
+ if(is.na(newrows$model)) newrows$model <- ifelse(newrows$state == "aq", "HKF", "CGL")
# Now check the formulas
e <- tryCatch(makeup(newrows$formula), error=function(e) e)
if(inherits(e, "error")) {
@@ -73,7 +75,7 @@
ntotal <- nrow(thermo$OBIGT)
ispecies[inew] <- (ntotal-length(inew)+1):ntotal
# Inform user
- message(paste("mod.OBIGT: added ", newrows$name, "(", newrows$state, ")", " with energy units of ", newrows$E_units, sep="", collapse="\n"))
+ message(paste("mod.OBIGT: added ", newrows$name, "(", newrows$state, ")", " with ", newrows$model, " model and energy units of ", newrows$E_units, sep="", collapse="\n"))
}
if(length(iold) > 0) {
# Loop over species
@@ -81,10 +83,12 @@
# The old values and the state
oldprop <- thermo$OBIGT[ispecies[iold[i]], icol]
state <- thermo$OBIGT$state[ispecies[iold[i]]]
- # Zap (clear) all preexisting values except for state 20220324
+ model <- thermo$OBIGT$model[ispecies[iold[i]]]
+ # Zap (clear) all preexisting values except for state and model 20220324
if(zap) {
thermo$OBIGT[ispecies[iold[i]], ] <- NA
thermo$OBIGT$state[ispecies[iold[i]]] <- state
+ thermo$OBIGT$model[ispecies[iold[i]]] <- model
}
# Tell user if they're the same, otherwise update the data entry
if(isTRUE(all.equal(oldprop, args[iold[i], ], check.attributes = FALSE)))
@@ -110,13 +114,6 @@
sysnosuffix <- sapply(strsplit(sysfiles, "\\."), "[", 1)
isys <- match(file, sysnosuffix)
if(!is.na(isys)) file <- system.file(paste0("extdata/OBIGT/", sysfiles[isys]), package="CHNOSZ")
-# else {
-# # We also match single system files with the state suffix removed
-# # (e.g. "DEW" for "DEW_aq", but not "organic" because we have "organic_aq", "organic_cr", etc.)
-# sysnostate <- sapply(strsplit(sysnosuffix, "_"), "[", 1)
-# isys <- which(file==sysnostate)
-# if(length(isys)==1) file <- system.file(paste0("extdata/OBIGT/", sysfiles[isys]), package="CHNOSZ")
-# }
# Read data from the file
to2 <- read.csv(file, as.is=TRUE)
Etxt <- paste(unique(to2$E_units), collapse = " and ")
Modified: pkg/CHNOSZ/R/cgl.R
===================================================================
--- pkg/CHNOSZ/R/cgl.R 2022-09-19 09:32:39 UTC (rev 743)
+++ pkg/CHNOSZ/R/cgl.R 2022-09-20 02:33:43 UTC (rev 744)
@@ -14,11 +14,11 @@
for(k in 1:nrow(parameters)) {
# The parameters for *this* species
PAR <- parameters[k, ]
- if(all(is.na(PAR[10:22]))) {
+ if(PAR$model == "Berman") {
# Use Berman equations (parameters not in thermo()$OBIGT)
properties <- Berman(PAR$name, T = T, P = P)
iprop <- match(property, colnames(properties))
- values <- properties[, iprop, drop=FALSE]
+ values <- properties[, iprop, drop = FALSE]
} else {
# In CHNOSZ, we have
# 1 cm^3 bar --> convert(1, "calories") == 0.02390057 cal
Modified: pkg/CHNOSZ/R/hkf.R
===================================================================
--- pkg/CHNOSZ/R/hkf.R 2022-09-19 09:32:39 UTC (rev 743)
+++ pkg/CHNOSZ/R/hkf.R 2022-09-20 02:33:43 UTC (rev 744)
@@ -38,6 +38,9 @@
if(grepl("DEW", thermo$opt$water)) {
# Using DEW model: get beta to calculate dgdP
H2O.props <- c(H2O.props, "beta")
+ } else {
+ # Don't allow DEW aqueous species if DEW water model isn't loaded 20220919
+ if(any(grepl("DEW", toupper(parameters$model)))) stop('The DEW water model is needed for one or more aqueous species. Run water("DEW") first.')
}
H2O <- water(H2O.props, T = c(Tr, T), P = c(Pr, P))
H2O.PrTr <- H2O[1, ]
Modified: pkg/CHNOSZ/R/subcrt.R
===================================================================
--- pkg/CHNOSZ/R/subcrt.R 2022-09-19 09:32:39 UTC (rev 743)
+++ pkg/CHNOSZ/R/subcrt.R 2022-09-20 02:33:43 UTC (rev 744)
@@ -11,6 +11,8 @@
#source("species.R")
#source("AD.R")
#source("nonideal.R")
+#source("hkf.R")
+#source("cgl.R")
subcrt <- function(species, coeff = 1, state = NULL, property = c("logK", "G", "H", "S", "V", "Cp"),
T = seq(273.15, 623.15, 25), P = "Psat", grid = NULL, convert = TRUE, exceed.Ttr = FALSE,
@@ -181,16 +183,23 @@
}
# Where we keep info about the species involved
+ # Add model information 20220919
+ model <- thermo$OBIGT$model[iphases]
+ # Label specific water model;
+ # this is also how we record a "wet" reaction
+ isH2O <- model == "H2O"
+ isH2O[is.na(isH2O)] <- FALSE
+ model[isH2O] <- paste0("water.", thermo$opt$water)
reaction <- data.frame(coeff = coeff.new, name = thermo$OBIGT$name[iphases],
formula = thermo$OBIGT$formula[iphases], state = thermo$OBIGT$state[iphases],
- ispecies = iphases, stringsAsFactors = FALSE)
+ ispecies = iphases, model = model, stringsAsFactors = FALSE)
# Make the rownames readable ... but they have to be unique
if(length(unique(iphases))==length(iphases)) rownames(reaction) <- as.character(iphases)
+ # Which are aqueous species
+ isaq <- reaction$state == "aq"
+ isaq.model <- toupper(reaction$model) %in% c("HKF", "AD", "DEW")
+ stopifnot(all(isaq == isaq.model))
- # Wetness etc.
- isH2O <- reaction$name=='water' & reaction$state=='liq'
- isaq <- reaction$state=='aq'
-
# Produce message about conditions
if(length(species)==1 & convert==FALSE) {
# No message produced here (internal calls from other functions)
@@ -282,12 +291,12 @@
# Aqueous species and H2O properties
if(TRUE %in% isaq) {
# 20110808 get species parameters using OBIGT2eos()
- # (this faster than using info() and is how we get everything in the same units)
+ # (this is faster than using info() and is how we get everything in the same units)
param <- OBIGT2eos(thermo$OBIGT[iphases[isaq],], "aq", fixGHS = TRUE, toJoules = TRUE)
- # Aqueous species with abbrv = "AD" use the AD model 20210407
- abbrv <- thermo$OBIGT$abbrv[iphases[isaq]]
- abbrv[is.na(abbrv)] <- ""
- isAD <- abbrv == "AD"
+ # Aqueous species with model = "AD" use the AD model 20210407
+ model <- thermo$OBIGT$model[iphases[isaq]]
+ model[is.na(model)] <- ""
+ isAD <- model == "AD"
# Always get density
H2O.props <- "rho"
# Calculate A_DH and B_DH if we're using the B-dot (Helgeson) equation
Modified: pkg/CHNOSZ/R/util.data.R
===================================================================
--- pkg/CHNOSZ/R/util.data.R 2022-09-19 09:32:39 UTC (rev 743)
+++ pkg/CHNOSZ/R/util.data.R 2022-09-20 02:33:43 UTC (rev 744)
@@ -404,17 +404,17 @@
# remove scaling factors from EOS parameters
# and apply column names depending on the EOS
if(identical(state, "aq")) {
- # Aqueous species with abbrv = "AD" use the AD model 20210407
- abbrv <- OBIGT$abbrv
- abbrv[is.na(abbrv)] <- ""
- isAD <- abbrv == "AD"
- # remove scaling factors for the HKF species, but not for the AD species
+ # Aqueous species with model = "AD" use the AD model 20210407
+ model <- OBIGT$model
+ model[is.na(model)] <- ""
+ isAD <- model == "AD"
+ # Remove scaling factors for the HKF species, but not for the AD species;
# protect this by an if statement to workaround error in subassignment to empty subset of data frame in R < 3.6.0
# (https://bugs.r-project.org/bugzilla/show_bug.cgi?id=17483) 20190302
if(any(!isAD)) OBIGT[!isAD, 15:22] <- t(t(OBIGT[!isAD, 15:22]) * 10^c(-1,2,0,4,0,4,5,0))
- # for AD species, set NA values in remaining columns (for display only)
+ # For AD species, set NA values in remaining columns (for display only)
if(any(isAD)) OBIGT[isAD, 18:21] <- NA
- # if all of the species are AD, change the variable names
+ # If all of the species are AD, change the variable names
if(all(isAD)) colnames(OBIGT)[15:22] <- c('a','b','xi','XX1','XX2','XX3','XX4','Z')
else colnames(OBIGT)[15:22] <- c('a1','a2','a3','a4','c1','c2','omega','Z')
} else {
@@ -434,14 +434,13 @@
}
}
if(fixGHS) {
- # fill in one of missing G, H, S
- # for use esp. by subcrt because NA for one of G, H or S
- # will preclude calculations at high T
- # which entries are missing just one
+ # Fill in one of missing G, H, S;
+ # for use esp. by subcrt because NA for one of G, H or S precludes calculations at high T
+ # Which entries are missing just one
imiss <- which(rowSums(is.na(OBIGT[, 10:12])) == 1)
if(length(imiss) > 0) {
for(i in 1:length(imiss)) {
- # calculate the missing value from the others
+ # Calculate the missing value from the others
ii <- imiss[i]
GHS <- as.numeric(GHS(as.character(OBIGT$formula[ii]), G = OBIGT[ii, 10], H = OBIGT[ii, 11], S = OBIGT[ii, 12],
E_units = ifelse(toJoules, "J", OBIGT$E_units[ii])))
Modified: pkg/CHNOSZ/inst/NEWS.Rd
===================================================================
--- pkg/CHNOSZ/inst/NEWS.Rd 2022-09-19 09:32:39 UTC (rev 743)
+++ pkg/CHNOSZ/inst/NEWS.Rd 2022-09-20 02:33:43 UTC (rev 744)
@@ -12,7 +12,7 @@
% links to vignettes 20220723
\newcommand{\viglink}{\ifelse{html}{\out{<a href="../CHNOSZ/doc/#1.html"><strong>#1.Rmd</strong></a>}}{\bold{#1.Rmd}}}
-\section{Changes in CHNOSZ version 1.9.9-35 (2022-09-19)}{
+\section{Changes in CHNOSZ version 1.9.9-36 (2022-09-20)}{
\subsection{MAJOR USER-VISIBLE CHANGES}{
\itemize{
@@ -19,12 +19,12 @@
\item A \code{model} column has been added to the OBIGT thermodynamic
database. This specifies the model for calculating standard thermodynamic
- properties and supersedes the heuristics for determining the model to use
- based on values in other columns. Currently used models are \samp{H2O}
- (water), \samp{HKF} (Helgeson-Kirkham-Flowers), \samp{CGL} (general heat
- capacity equation for crystalline, gas, and liquid species),
- \samp{Berman} (Berman equations), \samp{AD} (Akinfiev-Diamond), and
- \samp{DEW} (Deep-Earth Water model). All OBIGT data files that can be
+ properties and will eventually supersede the heuristics for determining
+ the model to use based on values in other columns. Currently used models
+ are \samp{H2O} (water), \samp{HKF} (Helgeson-Kirkham-Flowers), \samp{CGL}
+ (general heat capacity equation for crystalline, gas, and liquid
+ species), \samp{Berman} (Berman equations), \samp{AD} (Akinfiev-Diamond),
+ and \samp{DEW} (Deep Earth Water model). All OBIGT data files that can be
read using \code{add.OBIGT()} must have this column; there is no
back-compatibility support for data files in the older format.
@@ -49,7 +49,7 @@
\item Modify \code{OBIGT(no.organics = TRUE)} to prune the database of
organic compounds \emph{except} CH\s{4}. The ability to remove hundreds
of more complex organic compounds while keeping CH\s{4} is useful for
- some models of non-biological reactions in hydrothermal systems.
+ some calculations for abiotic reactions in hydrothermal systems.
\item Merge \file{biotic_aq.csv} into \file{organic_aq.csv}.
@@ -98,9 +98,9 @@
\subsection{NEW FEATURES: PROTEINS}{
\itemize{
- \item Add function and demo for \code{rank_affinity()} to calculate means
- of affinity rankings for specified groups of proteins, intended for
- evolutionary comparisons. The output of the new function can be used by
+ \item Add function \code{rank_affinity()} to calculate means of affinity
+ rankings for specified groups of proteins, intended for evolutionary
+ comparisons. The output of the new function can be used by
\code{diagram()}.
\item Add an \strong{as.residue} argument to \code{add.protein()} to
@@ -149,11 +149,11 @@
less call to \code{diagram()} in Mosaic Stacking 1 and 2 in
\strong{multi-metal.Rmd}, and now the entire chalcopyrite field in Mosaic
Stacking 2 is bounded by a thick orange line, instead of just the border
- shared with bornite.
+ with bornite.
\item Improve handling of non-integer coefficients in
- \code{expr.species()}. The result for FeS1.33 was previously equivalent
- to FeS[.33] but is now shown correctly as FeS[1.33].
+ \code{expr.species()}. The result for FeS1.33 was previously shown as
+ FeS[.33] but is now shown correctly as FeS[1.33].
}
}
@@ -161,10 +161,17 @@
\subsection{OTHER CHANGES}{
\itemize{
- \item As much as possible, functions in CHNOSZ have been modified to use
- Joules in internal variables. In particular, the value of the gas
- constant is based on Joules in \code{convert()} and other functions.
+ \item The \code{\var{reaction}} component of the output of
+ \code{subcrt()} contains the \code{model} for each species.
+ \item \code{subcrt()} produces an error if an aqueous species with the
+ \samp{DEW} \code{model} is requested but the DEW water model isn't
+ activated.
+
+ \item Most functions have been modified to use Joules in internal
+ variables. In particular, the value of the gas constant is expressed in
+ J K\S{-1} mol\S{-1} in \code{convert()} and other functions.
+
\item Add a \strong{zap} argument to \code{mod.OBIGT()} to clear
parameters of preexisting species (used by \code{logB_to_OBIGT()}).
@@ -171,8 +178,8 @@
\item Change default for \code{affinity()} to \code{loga.protein = 0}
(was -3).
- \item Add tests stack_mosaic.R and stack_solubility.R (these create PDF
- files for visual inspection of results).
+ \item Add tests \file{stack_mosaic.R} and \file{stack_solubility.R}
+ (these create PDF files for visual inspection of results).
\item Change license from GPL (>= 2) to GPL-3.
Modified: pkg/CHNOSZ/inst/extdata/OBIGT/AD.csv
===================================================================
--- pkg/CHNOSZ/inst/extdata/OBIGT/AD.csv 2022-09-19 09:32:39 UTC (rev 743)
+++ pkg/CHNOSZ/inst/extdata/OBIGT/AD.csv 2022-09-20 02:33:43 UTC (rev 744)
@@ -1,25 +1,25 @@
name,abbrv,formula,state,ref1,ref2,date,model,E_units,G,H,S,Cp,V,a1.a,a2.b,a3.c,a4.d,c1.e,c2.f,omega.lambda,z.T
-Ar,AD,Ar,aq,AD03.1,NA,2019-02-20,AD,J,NA,NA,NA,NA,NA,-8.5139,11.921,0.0733,NA,NA,NA,NA,NA
-H2S,AD,H2S,aq,AD03.1,NA,2019-02-20,AD,J,NA,NA,NA,NA,NA,-13.4046,13.8582,-0.2029,NA,NA,NA,NA,NA
-O2,AD,O2,aq,AD03.1,NA,2019-02-20,AD,J,NA,NA,NA,NA,NA,-9.754,12.9411,0.026,NA,NA,NA,NA,NA
-N2,AD,N2,aq,AD03.1,NA,2019-02-20,AD,J,NA,NA,NA,NA,NA,-11.538,14.6278,-0.032,NA,NA,NA,NA,NA
-NH3,AD,NH3,aq,AD03.1,NA,2019-02-20,AD,J,NA,NA,NA,NA,NA,-4.7245,4.9782,-0.0955,NA,NA,NA,NA,NA
-H2,AD,H2,aq,AD03.1,NA,2019-02-20,AD,J,NA,NA,NA,NA,NA,-8.4596,10.8301,0.309,NA,NA,NA,NA,NA
-CH4,AD,CH4,aq,AD03.1,NA,2019-02-20,AD,J,NA,NA,NA,NA,NA,-11.8462,14.8615,-0.1131,NA,NA,NA,NA,NA
-CO2,AD,CO2,aq,AD03.1,NA,2019-02-20,AD,J,NA,NA,NA,NA,NA,-8.8321,11.2684,-0.085,NA,NA,NA,NA,NA
-benzene,AD,C6H6,aq,AD03.1,NA,2019-02-20,AD,J,NA,NA,NA,NA,NA,-21.0084,22.934,-1.101,NA,NA,NA,NA,NA
-HCl,AD,HCl,aq,AD03.1,NA,2019-02-20,AD,J,NA,NA,NA,NA,NA,11.642,-7.4244,-0.28,NA,NA,NA,NA,NA
-Ne,AD,Ne,aq,AD03.2,NA,2019-02-21,AD,J,4565,NA,16.74,NA,20.4,1.0014,4.7976,0.5084,NA,NA,NA,NA,NA
-ethylene,AD,C2H4,aq,AD03.2,NA,2019-02-21,AD,J,19450,NA,28.7,NA,45.5,-16.8037,18.846,-0.4499,NA,NA,NA,NA,NA
-ethane,AD,C2H6,aq,AD03.2,NA,2019-02-21,AD,J,-4141,NA,26.75,NA,51.2,-16.3482,20.0628,-0.6091,NA,NA,NA,NA,NA
-propane,AD,C3H8,aq,AD03.2,NA,2019-02-21,AD,J,-2021,NA,33.49,NA,67,-25.3879,28.2616,-1.1471,NA,NA,NA,NA,NA
-butane,AD,C4H10,aq,AD03.2,NA,2019-02-21,AD,J,99,NA,39.66,NA,82.8,-33.8492,36.1457,-1.6849,NA,NA,NA,NA,NA
-benzene,AD,C6H6,aq,AD03.2,NA,2019-02-21,AD,J,32000,NA,35.62,NA,83.5,-39.109,37.5421,-1.9046,NA,NA,NA,NA,NA
-HF,AD,HF,aq,AD03.2,NA,2019-02-21,AD,J,-71662,NA,22.5,NA,12.5,3.0888,-3.5714,0.1008,NA,NA,NA,NA,NA
-SO2,AD,SO2,aq,AD03.2,NA,2019-02-21,AD,J,-71980,NA,38.7,NA,38.5,-14.5223,14.3512,-0.4295,NA,NA,NA,NA,NA
-B(OH)3,AD,B(OH)3,aq,AP14,NA,2019-02-22,AD,J,NA,NA,NA,NA,NA,-4.2561,4.0194,-1.057,NA,NA,NA,NA,NA
-Si(OH)4,AD,Si(OH)4,aq,AP14,NA,2019-02-22,AD,J,NA,NA,NA,NA,NA,0.9285,-0.9409,-1.8933,NA,NA,NA,NA,NA
-As(OH)3,AD,As(OH)3,aq,AP14,NA,2019-02-22,AD,J,NA,NA,NA,NA,NA,-9.903,7.6818,-1.23,NA,NA,NA,NA,NA
+Ar,NA,Ar,aq,AD03.1,NA,2019-02-20,AD,J,NA,NA,NA,NA,NA,-8.5139,11.921,0.0733,NA,NA,NA,NA,NA
+H2S,NA,H2S,aq,AD03.1,NA,2019-02-20,AD,J,NA,NA,NA,NA,NA,-13.4046,13.8582,-0.2029,NA,NA,NA,NA,NA
+O2,NA,O2,aq,AD03.1,NA,2019-02-20,AD,J,NA,NA,NA,NA,NA,-9.754,12.9411,0.026,NA,NA,NA,NA,NA
+N2,NA,N2,aq,AD03.1,NA,2019-02-20,AD,J,NA,NA,NA,NA,NA,-11.538,14.6278,-0.032,NA,NA,NA,NA,NA
+NH3,NA,NH3,aq,AD03.1,NA,2019-02-20,AD,J,NA,NA,NA,NA,NA,-4.7245,4.9782,-0.0955,NA,NA,NA,NA,NA
+H2,NA,H2,aq,AD03.1,NA,2019-02-20,AD,J,NA,NA,NA,NA,NA,-8.4596,10.8301,0.309,NA,NA,NA,NA,NA
+CH4,NA,CH4,aq,AD03.1,NA,2019-02-20,AD,J,NA,NA,NA,NA,NA,-11.8462,14.8615,-0.1131,NA,NA,NA,NA,NA
+CO2,NA,CO2,aq,AD03.1,NA,2019-02-20,AD,J,NA,NA,NA,NA,NA,-8.8321,11.2684,-0.085,NA,NA,NA,NA,NA
+benzene,NA,C6H6,aq,AD03.1,NA,2019-02-20,AD,J,NA,NA,NA,NA,NA,-21.0084,22.934,-1.101,NA,NA,NA,NA,NA
+HCl,NA,HCl,aq,AD03.1,NA,2019-02-20,AD,J,NA,NA,NA,NA,NA,11.642,-7.4244,-0.28,NA,NA,NA,NA,NA
+Ne,NA,Ne,aq,AD03.2,NA,2019-02-21,AD,J,4565,NA,16.74,NA,20.4,1.0014,4.7976,0.5084,NA,NA,NA,NA,NA
+ethylene,NA,C2H4,aq,AD03.2,NA,2019-02-21,AD,J,19450,NA,28.7,NA,45.5,-16.8037,18.846,-0.4499,NA,NA,NA,NA,NA
+ethane,NA,C2H6,aq,AD03.2,NA,2019-02-21,AD,J,-4141,NA,26.75,NA,51.2,-16.3482,20.0628,-0.6091,NA,NA,NA,NA,NA
+propane,NA,C3H8,aq,AD03.2,NA,2019-02-21,AD,J,-2021,NA,33.49,NA,67,-25.3879,28.2616,-1.1471,NA,NA,NA,NA,NA
+butane,NA,C4H10,aq,AD03.2,NA,2019-02-21,AD,J,99,NA,39.66,NA,82.8,-33.8492,36.1457,-1.6849,NA,NA,NA,NA,NA
+benzene,NA,C6H6,aq,AD03.2,NA,2019-02-21,AD,J,32000,NA,35.62,NA,83.5,-39.109,37.5421,-1.9046,NA,NA,NA,NA,NA
+HF,NA,HF,aq,AD03.2,NA,2019-02-21,AD,J,-71662,NA,22.5,NA,12.5,3.0888,-3.5714,0.1008,NA,NA,NA,NA,NA
+SO2,NA,SO2,aq,AD03.2,NA,2019-02-21,AD,J,-71980,NA,38.7,NA,38.5,-14.5223,14.3512,-0.4295,NA,NA,NA,NA,NA
+B(OH)3,NA,B(OH)3,aq,AP14,NA,2019-02-22,AD,J,NA,NA,NA,NA,NA,-4.2561,4.0194,-1.057,NA,NA,NA,NA,NA
+Si(OH)4,NA,Si(OH)4,aq,AP14,NA,2019-02-22,AD,J,NA,NA,NA,NA,NA,0.9285,-0.9409,-1.8933,NA,NA,NA,NA,NA
+As(OH)3,NA,As(OH)3,aq,AP14,NA,2019-02-22,AD,J,NA,NA,NA,NA,NA,-9.903,7.6818,-1.23,NA,NA,NA,NA,NA
B(OH)3,NA,B(OH)3,gas,AP14,NA,2019-02-22,CGL,cal,-222523.9,NA,66.589,NA,0,46.7734,-0.097,0.5932,-522.4665,0,0,0,1000
Si(OH)4,NA,Si(OH)4,gas,AP14,NA,2019-02-22,CGL,cal,-296285.9,NA,82.784,NA,0,22.9685,37.6434,-4.4407,0,-2.7782,7.923E-09,3,1000
As(OH)3,NA,As(OH)3,gas,AP14,NA,2019-02-22,CGL,cal,-143112.8,NA,73.599,NA,0,22.6338,-16.3623,-3.7213,0,-1.0194,2.7414E-09,3,1000
Modified: pkg/CHNOSZ/inst/tinytest/test-AD.R
===================================================================
--- pkg/CHNOSZ/inst/tinytest/test-AD.R 2022-09-19 09:32:39 UTC (rev 743)
+++ pkg/CHNOSZ/inst/tinytest/test-AD.R 2022-09-20 02:33:43 UTC (rev 744)
@@ -8,7 +8,7 @@
# Get the indices of the modified aqueous species
iaq <- iAD[info(iAD)$state == "aq"]
# Each aqueous species should be associated with the AD model
-expect_equal(unique(info(iaq)$abbrv), "AD", info = info)
+expect_equal(unique(info(iaq)$model), "AD", info = info)
# Each aqueous species should have a gaseous counterpart
igas <- info(info(iAD)$name, "gas")
expect_true(!any(is.na(igas)), info = info)
Modified: pkg/CHNOSZ/inst/tinytest/test-eos.R
===================================================================
--- pkg/CHNOSZ/inst/tinytest/test-eos.R 2022-09-19 09:32:39 UTC (rev 743)
+++ pkg/CHNOSZ/inst/tinytest/test-eos.R 2022-09-20 02:33:43 UTC (rev 744)
@@ -81,6 +81,8 @@
# in demo/DEW.R (b/c hkf() is no longer exported)
# Load SiO2 and Si2O4 data taken from DEW spreadsheet
iSi <- add.OBIGT("DEW", c("SiO2", "Si2O4"))
+# Override check for DEW water model 20220920
+mod.OBIGT(iSi, model = rep("HKF", 2))
Vn1 <- Vn2 <- numeric()
species <- c("CO3-2", "BO2-", "MgCl+", "SiO2", "HCO3-", "Si2O4")
# First method: use hkf() function
Modified: pkg/CHNOSZ/inst/tinytest/test-info.R
===================================================================
--- pkg/CHNOSZ/inst/tinytest/test-info.R 2022-09-19 09:32:39 UTC (rev 743)
+++ pkg/CHNOSZ/inst/tinytest/test-info.R 2022-09-20 02:33:43 UTC (rev 744)
@@ -33,7 +33,7 @@
info <- "info() gives correct column names for species using the AD model"
# add an aqueous species conforming to the AD model
-iCO2 <- mod.OBIGT("CO2", abbrv = "AD", a = -8.8321, b = 11.2684, c = -0.0850)
+iCO2 <- mod.OBIGT("CO2", model = "AD", a = -8.8321, b = 11.2684, c = -0.0850)
params <- info(iCO2)
expect_equal(params$a, -8.8321, info = info)
expect_equal(params$b, 11.2684, info = info)
Modified: pkg/CHNOSZ/man/add.OBIGT.Rd
===================================================================
--- pkg/CHNOSZ/man/add.OBIGT.Rd 2022-09-19 09:32:39 UTC (rev 743)
+++ pkg/CHNOSZ/man/add.OBIGT.Rd 2022-09-20 02:33:43 UTC (rev 744)
@@ -50,10 +50,13 @@
The formula is taken from the \samp{formula} argument, or if that is missing, is taken to be the same as the \samp{name} of the species.
An error occurs if the formula is not valid (i.e. can not be parsed by \code{\link{makeup}}).
For new species, properties that are not specified become NA, except for \samp{state} and \samp{E_units}, which take default values from \code{thermo()$opt}.
-These defaults can be overridden by giving a value for \samp{state} or \samp{E_units} in the new data.
+These defaults can be overridden by giving a value for \samp{state} or \samp{E_units} in the arguments.
+\samp{model}, if missing, is set to \samp{HKF} for \code{state == "aq"} or \samp{CGL} otherwise.
+When modifying some existing minerals in OBIGT, \code{model = "CGL"} should be explicitly given in order to override the Berman model.
+
When modifying species, the parameters indicated by the named arguments of \code{mod.OBIGT} are updated.
-Use \code{zap = TRUE} to replace all prexisting parameters (except for \code{state}) with NA values.
+Use \code{zap = TRUE} to replace all prexisting parameters (except for \code{state} and \code{model}) with NA values.
}
\value{
Modified: pkg/CHNOSZ/man/subcrt.Rd
===================================================================
--- pkg/CHNOSZ/man/subcrt.Rd 2022-09-19 09:32:39 UTC (rev 743)
+++ pkg/CHNOSZ/man/subcrt.Rd 2022-09-20 02:33:43 UTC (rev 744)
@@ -33,9 +33,15 @@
\details{
\code{subcrt} calculates the standard molal thermodynamic properties of species and reactions as a function of temperature and pressure.
For each of the \code{species} (a formula or name), optionally identified in a given \code{state}, the standard molal thermodynamic properties and equations-of-state parameters are retrieved via \code{\link{info}} (except for \H2O liquid).
-The standard molal properties of the species are computed using the revised Helgeson-Kirkham-Flowers (HKF) equations of state for aqueous species, a heat capacity equation for crystalline, gas, and liquid species (CGL), and dedicated functions for liquid and supercritical \H2O (\code{\link{water}}).
+The standard molal properties of the species are computed using the thermodynamic \code{model} given for each species (see \code{\link{thermo}}).
+This function also calculates the thermodynamic properties of \emph{reactions} by summing those of the respective species.
+This functionality is inspired the \acronym{SUPCRT92} package (Johnson et al., 1992).
-\code{T} and \code{P} denote the temperature and pressure conditions for the calculations and should generally be of the same length, unless \code{P} is \samp{Psat} (the default; see \code{\link{water}}).
+\code{T} and \code{P} denote the temperature and pressure for the calculations.
+The only valid non-numeric value is \samp{Psat} for \code{P}, which is the default (see \code{\link{water}}).
+For calculations below 273.16 K, \code{P} set to 1, as \Psat is not defined at subzero (\degC) temperatures.
+At temperatures above the critical point of water, \code{P} must be set to a numeric value; unless \code{exceed.rhomin} is TRUE, this should correspond to a fluid density \ge 0.35 g cm\S{-3}.
+
Argument \code{\link{grid}} if present can be one of \code{T} or \code{P} to perform the computation of a \code{T}\eqn{\times}{x}\code{P} or \code{P}\eqn{\times}{x}\code{T} grid.
The \code{property}s to be calculated can be one or more of those shown below:
@@ -51,10 +57,6 @@
\code{kT} \tab Isothermal compressibility \tab cm\eqn{^3} bar\eqn{^{-1}}{^-1} \cr
}
-Note that \code{E} and \code{kT} can only be calculated for aqueous species and only if the option (\code{\link{thermo}()$opt$water}) for calculations of properties using \code{water} is set to \code{IAPWS}.
-On the other hand, if the \code{water} option is \samp{SUPCRT} (the default), \code{E} and \code{kT} can be calculated for water but not for aqueous species.
-(This is not an inherent limitation in either formulation, but it is just a matter of implementation.)
-
If \code{convert} is \code{TRUE} (the default), the input values of \code{T} and \code{P} are interpreted to have the units given by \code{\link{T.units}} and \code{\link{P.units}} (default: \degC and bar), and the output values of \code{G}, \code{H}, \code{S} and \code{Cp} are based on the units given in \code{\link{E.units}} (default: Joules).
If \code{convert} is \code{FALSE}, the user units (\code{\link{T.units}}, \code{\link{P.units}}, and \code{\link{E.units}}) are ignored, and \code{T} and \code{P} are taken to be in Kelvin and bar, and the returned values of \code{G}, \code{H}, \code{S} and \code{Cp} are in Joules.
@@ -64,12 +66,6 @@
Alternatively, if \code{autobalance} is \code{TRUE}, the \code{\link{basis}} species of a system were previously defined, and all elements in the reaction are represented by the basis species, an unbalanced reaction given in the arguments to \code{subcrt} will be balanced automatically.
The auto balancing doesn't change the reaction coefficients of any species that are do not correspond to basis species.
-Minerals with polymorphic transitions (denoted by having states \samp{cr} (lowest T phase), \samp{cr2}, \samp{cr3} etc.) can be defined generically by \samp{cr} in the \code{state} argument with a character value for \code{species}.
-\code{subcrt} in this case simultaneously calculates the requested properties of all the phases of each such mineral (phase species) and, using the values of the transition temperatures calculated from those at P = 1 bar given in the thermodynamic database together with functions of the entropies and volumes of transitions (see \code{\link{dPdTtr}}), determines the stable phase of the mineral at any grid point and substitutes the properties of this phase at that point for further calculations.
-If individual phase species of minerals are specified (by \samp{cr}, \samp{cr2} etc. in \code{state}), and \code{exceed.Ttr} is \code{FALSE} (the default), the Gibbs energies of these minerals are assigned values of NA at conditions beyond their transition temperature, where another of the phases is stable.
-If you set \code{exceed.Ttr} to \code{TRUE} to calculate the properties of mineral polymorphs (i.e., using \samp{cr}) the function will identify the stable polymorph using the calculated Gibbs energies of the phase species instead of the tabulated transition temperatures.
-This is not generally advised, since the computed standard molal properties of a phase species lose their physical meaning beyond the transition temperatures of the phase.
-
If \code{logact} is provided, the chemical affinities of reactions are calculated.
\code{logact} indicates the logarithms of activities (fugacities for gases) of species in the reaction; if there are fewer values of \code{logact} than number of species those values are repeated as necessary.
If the reaction was unbalanced to start, the logarithms of activities of any basis species added to the reaction are taken from the current definition of the \code{\link{basis}} species.
@@ -82,11 +78,6 @@
Calculations can also be performed on a \code{P}-\code{IS}, \code{T}-\code{IS} or \code{P,T}-\code{IS} grid.
Values of \code{logK} of reactions calculated for \code{IS} not equal to zero are consistent with the adjusted Gibbs energies of the charged aqueous species.
-\code{subcrt} is modeled after the functionality of the \acronym{SUPCRT92} package (Johnson et al., 1992).
-Certain features of \acronym{SUPCRT92} are not available here, for example, calculations as a function of density of \H2O instead of pressure, or calculations of temperatures of univariant curves (i.e. for which \code{logK} is zero).
-
-For calculations below 273.16 K, the pressure should be set to 1, as \Psat is not defined in these conditions.
-
If \code{thermo()$opt$varP} is \code{TRUE}, standard Gibbs energies of gases will be converted from a standard state at 1 bar (as used in SUPCRT) to a variable pressure standard state (see chapter 12 in Anderson and Crerar, 1993).
This is useful for constructing e.g. boiling curves for organic compounds.
}
@@ -103,6 +94,21 @@
\code{NA}s are also output if the T, P conditions are otherwise beyond the capabilities of the water equations of state derived from SUPCRT92 (H2O92D.f), but the messages about this are produced by \code{\link{water.SUPCRT92}} rather than \code{subcrt}.
}
+\section{Limitations}{
+Note that \code{E} and \code{kT} can only be calculated for aqueous species and only if the option (\code{\link{thermo}()$opt$water}) for calculations of properties using \code{water} is set to \code{IAPWS}.
+On the other hand, if the \code{water} option is \samp{SUPCRT} (the default), \code{E} and \code{kT} can be calculated for water but not for aqueous species.
+(This is not an inherent limitation in either formulation, but it is just a matter of implementation.)
+}
+
+\section{Note on phase transitions}{
+Minerals with polymorphic transitions (denoted by having states \samp{cr} (lowest T phase), \samp{cr2}, \samp{cr3} etc.) can be defined generically by \samp{cr} in the \code{state} argument with a character value for \code{species}.
+\code{subcrt} in this case simultaneously calculates the requested properties of all the phases of each such mineral (phase species) and, using the values of the transition temperatures calculated from those at P = 1 bar given in the thermodynamic database together with functions of the entropies and volumes of transitions (see \code{\link{dPdTtr}}), determines the stable phase of the mineral at any grid point and substitutes the properties of this phase at that point for further calculations.
+If individual phase species of minerals are specified (by \samp{cr}, \samp{cr2} etc. in \code{state}), and \code{exceed.Ttr} is \code{FALSE} (the default), the Gibbs energies of these minerals are assigned values of NA at conditions beyond their transition temperature, where another of the phases is stable.
+If you set \code{exceed.Ttr} to \code{TRUE} to calculate the properties of mineral polymorphs (i.e., using \samp{cr}) the function will identify the stable polymorph using the calculated Gibbs energies of the phase species instead of the tabulated transition temperatures.
+This is not generally advised, since the computed standard molal properties of a phase species lose their physical meaning beyond the transition temperatures of the phase.
+
+}
+
\value{
For \code{subcrt}, a list of length two or three.
If the properties of a reaction were calculated, the first element of the list (named \samp{reaction}) contains a dataframe with the reaction parameters; the second element, named \samp{out}, is a dataframe containing the calculated properties.
Modified: pkg/CHNOSZ/man/thermo.Rd
===================================================================
--- pkg/CHNOSZ/man/thermo.Rd 2022-09-19 09:32:39 UTC (rev 743)
+++ pkg/CHNOSZ/man/thermo.Rd 2022-09-20 02:33:43 UTC (rev 744)
@@ -20,7 +20,8 @@
All the data are stored in the \code{thermo} data object in an environment named \code{CHNOSZ}.
\code{thermo()} is a convenience function to access or modify parts of this object, in particular some computational settings, for example, \code{thermo("opt$ideal.H" = FALSE)} (see \code{\link{nonideal}}).
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/chnosz -r 744
More information about the CHNOSZ-commits
mailing list