[CHNOSZ-commits] r212 - in pkg/CHNOSZ: . R inst man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Sep 25 08:49:54 CEST 2017
Author: jedick
Date: 2017-09-25 08:49:53 +0200 (Mon, 25 Sep 2017)
New Revision: 212
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/R/EOSregress.R
pkg/CHNOSZ/R/hkf.R
pkg/CHNOSZ/R/subcrt.R
pkg/CHNOSZ/inst/NEWS
pkg/CHNOSZ/man/eos.Rd
Log:
subcrt(): simplify calculations of water properties
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2017-09-25 04:12:03 UTC (rev 211)
+++ pkg/CHNOSZ/DESCRIPTION 2017-09-25 06:49:53 UTC (rev 212)
@@ -1,6 +1,6 @@
-Date: 2017-09-24
+Date: 2017-09-25
Package: CHNOSZ
-Version: 1.1.0-10
+Version: 1.1.0-11
Title: Thermodynamic Calculations for Geobiochemistry
Author: Jeffrey Dick
Maintainer: Jeffrey Dick <j3ffdick at gmail.com>
Modified: pkg/CHNOSZ/R/EOSregress.R
===================================================================
--- pkg/CHNOSZ/R/EOSregress.R 2017-09-25 04:12:03 UTC (rev 211)
+++ pkg/CHNOSZ/R/EOSregress.R 2017-09-25 06:49:53 UTC (rev 212)
@@ -5,14 +5,14 @@
Cp_s_var <- function(T=298.15, P=1, omega.PrTr=0, Z=0) {
# solvation contribution to heat capacity in the HKF EOS, divided by omega(Pr,Tr) (calories)
- Cp_s <- hkf("Cp", T=T, P=P, parameters=data.frame(omega=omega.PrTr, Z=Z), contrib="s")
+ Cp_s <- hkf("Cp", T=T, P=P, parameters=data.frame(omega=omega.PrTr, Z=Z), contrib="s")$aq
return(Cp_s[[1]][, 1]/omega.PrTr)
}
V_s_var <- function(T=298.15, P=1, omega.PrTr=0, Z=0) {
# solvation contribution to volume in the HKF EOS, divided by omega(Pr,Tr) (cm3.bar)
# [the negative sign on this term as written in the HKF EOS is accounted for by hkf()]
- V_s <- hkf("V", T=T, P=P, parameters=data.frame(omega=omega.PrTr, Z=Z), contrib="s")
+ V_s <- hkf("V", T=T, P=P, parameters=data.frame(omega=omega.PrTr, Z=Z), contrib="s")$aq
return(V_s[[1]][, 1]/convert(omega.PrTr, "cm3bar"))
}
Modified: pkg/CHNOSZ/R/hkf.R
===================================================================
--- pkg/CHNOSZ/R/hkf.R 2017-09-25 04:12:03 UTC (rev 211)
+++ pkg/CHNOSZ/R/hkf.R 2017-09-25 06:49:53 UTC (rev 212)
@@ -3,9 +3,10 @@
# 11/17/03 jmd
hkf <- function(property = NULL, T = 298.15, P = 1, parameters = NULL,
- contrib = c('n', 's', 'o'), H2O.PT = NULL) {
+ contrib = c('n', 's', 'o'), H2O.props="rho") {
# calculate G, H, S, Cp, V, kT, and/or E using
# the revised HKF equations of state
+ # H2O.props - H2O properties needed for subcrt() output
thermo <- get("thermo")
# constants
Tr <- thermo$opt$Tr
@@ -19,28 +20,28 @@
EOS.Props <- eargs$Prop
# nonsolvation, solvation, and origination contribution
notcontrib <- ! contrib %in% c('n', 's', 'o')
- if(TRUE %in% notcontrib)
- stop(paste("contrib must be in c('n', 's', 'o'); got", c2s(contrib[notcontrib])))
- # get water properties, if they weren't supplied in arguments (and we want solvation props)
- dosupcrt <- thermo$opt$water != "IAPWS95"
- if('s' %in% contrib) {
- H2O.props <- c("QBorn", "XBorn", "YBorn", "diel")
- if(dosupcrt) {
- # using H2O92D.f from SUPCRT92
- # (rho, alpha, daldT, beta - for partial derivatives of omega (g function))
- H2O.props <- c(H2O.props, "rho", "alpha", "daldT", "beta")
- } else {
- # using IAPWS-95
- # (NBorn, UBorn - for compressibility, expansibility)
- H2O.props <- c(H2O.props, 'NBorn', 'UBorn')
- }
- if(is.null(H2O.PT)) H2O.PT <- water(H2O.props, T = T, P = P)
- H2O.PrTr <- water(H2O.props, T = thermo$opt$Tr, P = thermo$opt$Pr)
- ZBorn <- -1 / H2O.PT$diel
- ZBorn.PrTr <- -1 / H2O.PrTr$diel
+ if(TRUE %in% notcontrib) stop(paste("contrib must be in c('n', 's', 'o'); got", c2s(contrib[notcontrib])))
+ # get water properties
+ # rho - for subcrt() output and derivatives of omega (if needed)
+ # Born functions and epsilon - for HKF calculations
+ H2O.props <- c(H2O.props, "QBorn", "XBorn", "YBorn", "diel")
+ if(grepl("SUPCRT", thermo$opt$water)) {
+ # using H2O92D.f from SUPCRT92
+ # (alpha, daldT, beta - for partial derivatives of omega (g function))
+ H2O.props <- c(H2O.props, "alpha", "daldT", "beta")
}
+ if(grepl("IAPWS", thermo$opt$water)) {
+ # using IAPWS-95
+ # (NBorn, UBorn - for compressibility, expansibility)
+ H2O.props <- c(H2O.props, 'NBorn', 'UBorn')
+ }
+ H2O <- water(H2O.props, T = c(thermo$opt$Tr, T), P = c(thermo$opt$Pr, P))
+ H2O.PrTr <- H2O[1, ]
+ H2O.PT <- H2O[-1, ]
+ ZBorn <- -1 / H2O.PT$diel
+ ZBorn.PrTr <- -1 / H2O.PrTr$diel
# a list to store the result
- x <- list()
+ aq.out <- list()
nspecies <- nrow(parameters)
for(k in 1:nspecies) {
# loop over each species
@@ -65,7 +66,7 @@
dwdP <- dwdT <- d2wdT2 <- numeric(length(T))
Z <- PAR$Z
omega.PT <- rep(PAR$omega, length(T))
- if(!identical(Z, 0) & !PAR$name=="H+" & dosupcrt) {
+ if(!identical(Z, 0) & !identical(PAR$name, "H+") & grepl("SUPCRT", thermo$opt$water)) {
# compute derivatives of omega: g and f functions (Shock et al., 1992; Johnson et al., 1992)
rhohat <- H2O.PT$rho/1000 # just converting kg/m3 to g/cm3
g <- gfun(rhohat, convert(T, "C"), P, H2O.PT$alpha, H2O.PT$daldT, H2O.PT$beta)
@@ -144,7 +145,7 @@
H2O.PT$QBorn + convert(dwdP,'cm3bar') * (-ZBorn - 1)
# TODO: the partial derivatives of omega are not included here here for kt and e
# (to do it, see p. 820 of SOJ+92 ... but kt requires d2wdP2 which we don't have yet)
- if(PROP == 'kt') p <- convert(omega,'cm3bar') * H2O.PT$N
+ if(PROP == 'kt') p <- convert(omega,'cm3bar') * H2O.PT$NBorn
if(PROP == 'e') p <- -convert(omega,'cm3bar') * H2O.PT$UBorn
}
if(icontrib == 'o') {
@@ -166,9 +167,9 @@
if(i > 1) w <- cbind(w, wnew) else w <- wnew
}
colnames(w) <- EOS.Props
- x[[k]] <- w
+ aq.out[[k]] <- w
}
- return(x)
+ return(list(aq=aq.out, H2O=H2O.PT))
}
### unexported functions ###
Modified: pkg/CHNOSZ/R/subcrt.R
===================================================================
--- pkg/CHNOSZ/R/subcrt.R 2017-09-25 04:12:03 UTC (rev 211)
+++ pkg/CHNOSZ/R/subcrt.R 2017-09-25 06:49:53 UTC (rev 212)
@@ -2,6 +2,13 @@
# calculate standard molal thermodynamic propertes
# 20060817 jmd
+## if we interactively source this file, the following are also needed to provide unexported functions:
+#source("util.args.R")
+#source("util.character.R")
+#source("info.R")
+#source("util.units.R")
+#source("util.data.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, check.Ttr=TRUE, exceed.Ttr=FALSE,
logact=NULL, action.unbalanced='warn', IS=0) {
@@ -42,8 +49,8 @@
# allowed properties
properties <- c('rho','logK','G','H','S','Cp','V','kT','E')
# property checking
- prop <- tolower(property)
- notproperty <- property[!prop %in% tolower(properties)]
+ outprop <- tolower(property)
+ notproperty <- property[!outprop %in% tolower(properties)]
if(length(notproperty) > 0) stop(paste(notproperty,
'are not valid properties\ntry rho, logK, G, H, S, V, Cp, kT, or E (or their lowercase equivalents)'))
# length checking
@@ -254,41 +261,35 @@
}
# calculate the properties
- # if we want affinities we must have logK
- if(!is.null(logact)) if(!'logk' %in% prop) prop <- c('logk',prop)
- # if logK but not g was requested, get g ...
- if('logk' %in% prop & ! 'g' %in% prop) eprop <- c(prop,'g') else eprop <- prop
- # don't request logk from the eos ...
- eosprop <- eprop[!eprop %in% c('logk','rho')]
+ # if we want affinities we must have logK; include it in the ouput
+ if(!is.null(logact)) if(!'logk' %in% outprop) outprop <- c('logk', outprop)
+ # if logK but not g was requested, we need to calculate g
+ eosprop <- outprop
+ if('logk' %in% outprop & ! 'g' %in% outprop) eosprop <- c(eosprop, 'g')
# also get g if we are dealing with mineral phases
- if(!'g' %in% eprop & length(inpho) > length(sinfo)) eosprop <- c(eosprop,'g')
- # the reaction result is in out
+ if(!'g' %in% eosprop & length(inpho) > length(sinfo)) eosprop <- c(eosprop, 'g')
+ # don't request logk or rho from the eos ...
+ eosprop <- eosprop[!eosprop %in% c('logk','rho')]
+ # the reaction result will go here
out <- list()
- # aqueous species
- if(TRUE %in% isaq | 'rho' %in% eprop) {
- # load the water properties
- wprop.PT <- character()
- wprop.PrTr <- 'rho'
- dosupcrt <- thermo$opt$water != "IAPWS95"
- if(TRUE %in% (prop %in% c('logk','g','h','s'))) wprop.PrTr <- c(wprop.PrTr,'YBorn')
- if(dosupcrt | TRUE %in% (prop %in% c('logk','g','h'))) wprop.PrTr <- c(wprop.PrTr,'diel')
- if(TRUE %in% (prop %in% c('cp'))) {wprop.PT <- c(wprop.PT,'XBorn','YBorn')}
- if(TRUE %in% (prop %in% c('v'))) {wprop.PT <- c(wprop.PT,'QBorn')}
- if(TRUE %in% (prop %in% c('kt'))) {wprop.PT <- c(wprop.PT,'NBorn')}
- if(TRUE %in% (prop %in% c('e'))) {wprop.PT <- c(wprop.PT,'UBorn')}
- # get additional properties required for omega derivatives
- if(dosupcrt) wprop.PT <- c(wprop.PT,'alpha','daldT','beta','diel')
- H2O.PT <- water(c(wprop.PrTr,wprop.PT),T=T,P=P)
- if(TRUE %in% isaq) {
- # now the species stuff
- # 20110808 if inpho are the species indices let's avoid
- # the overhead of info() and use new obigt2eos() instead
- #si <- info(inpho[isaq],quiet=TRUE)
- si <- obigt2eos(thermo$obigt[inpho[isaq],], "aq", fixGHS = TRUE)
- p.aq <- hkf(eosprop, T=T, P=P, parameters=si, H2O.PT=H2O.PT)
- if(any(IS!=0)) p.aq <- nonideal(inpho[isaq],p.aq,newIS,T)
- out <- c(out,p.aq)
- }
+ # aqueous species and H2O properties
+ if(TRUE %in% isaq) {
+ # 20110808 if inpho are the species indices let's avoid
+ # the overhead of info() and use new obigt2eos() instead
+ #si <- info(inpho[isaq],quiet=TRUE)
+ si <- obigt2eos(thermo$obigt[inpho[isaq],], "aq", fixGHS = TRUE)
+ # always get density
+ H2O.props <- "rho"
+ # get other properties for H2O only if it's in the reaction
+ if(any(isH2O)) H2O.props <- c(H2O.props, eosprop)
+ hkfstuff <- hkf(eosprop, T = T, P = P, parameters = si, H2O.props=H2O.props)
+ p.aq <- hkfstuff$aq
+ H2O.PT <- hkfstuff$H2O
+ if(any(IS != 0)) p.aq <- nonideal(inpho[isaq], p.aq, newIS, T)
+ out <- c(out, p.aq)
+ } else if(any(isH2O)) {
+ # we're not using the HKF, but still want water
+ H2O.PT <- water(c("rho", eosprop), T = T, P = P)
}
# crystalline, gas, liquid (except water) species
iscgl <- reaction$state %in% c('liq','cr','gas','cr1','cr2','cr3',
@@ -349,22 +350,20 @@
}
# water
- if(TRUE %in% isH2O) {
- if(!exists('H2O.PT',inherits=FALSE)) H2O.PT <- water('rho',T=T,P=P)
- if(length(eosprop)==0) eosprop <- 'rho'
- #message(paste('subcrt: water equation of state:',c2s(eosprop)))
- p.H2O <- list(tmp=water(eosprop,T=T,P=P))
- out <- c(out,rep(p.H2O,length(which(isH2O==TRUE))))
+ if(any(isH2O)) {
+ p.H2O <- H2O.PT[, match(eosprop, tolower(colnames(H2O.PT))), drop=FALSE]
+ p.H2O <- list(p.H2O)
+ out <- c(out, rep(p.H2O, sum(isH2O == TRUE)))
}
# use variable-pressure standard Gibbs energy for gases
isgas <- reaction$state %in% "gas"
- if(TRUE %in% isgas & "g" %in% eprop & thermo$opt$varP) {
+ if(any(isgas) & "g" %in% eosprop & thermo$opt$varP) {
for(i in which(isgas)) out[[i]]$G <- out[[i]]$G - convert(log10(P), "G", T=T)
}
# logK
- if('logk' %in% prop) {
+ if('logk' %in% outprop) {
for(i in 1:length(out)) {
# NOTE: the following depends on the water function renaming g to G
out[[i]] <- cbind(out[[i]],data.frame(logK=convert(out[[i]]$G,'logK',T=T)))
@@ -460,11 +459,10 @@
iscgl <- iscgl.new
isH2O <- isH2O.new
- newprop <- eprop[eprop!='rho']
# the order of the properties
- if(length(newprop)>1) for(i in 1:length(out)) {
+ if(length(outprop)>1) for(i in 1:length(out)) {
# keep state/loggam columns if they exists
- ipp <- match(newprop,tolower(colnames(out[[i]])))
+ ipp <- match(outprop,tolower(colnames(out[[i]])))
if('state' %in% colnames(out[[i]])) ipp <- c(ipp,match('state',colnames(out[[i]])))
if('loggam' %in% colnames(out[[i]])) ipp <- c(ipp,match('loggam',colnames(out[[i]])))
out[[i]] <- out[[i]][,ipp,drop=FALSE]
@@ -527,7 +525,7 @@
if(convert) P.out <- outvert(P,"bar") else P.out <- P
# try to stuff in a column of rho if we have aqueous species
# watch out! supcrt-ish densities are in g/cc not kg/m3
- if('rho' %in% prop | (missing(property) & any(c(isaq,isH2O))) & (names(out)[i])!='state')
+ if('rho' %in% outprop | (missing(property) & any(c(isaq,isH2O))) & (names(out)[i])!='state')
out[[i]] <- cbind(data.frame(T=T.out,P=P.out,rho=H2O.PT$rho/1000),out[[i]])
else
out[[i]] <- cbind(data.frame(T=T.out,P=P.out,out[[i]]))
Modified: pkg/CHNOSZ/inst/NEWS
===================================================================
--- pkg/CHNOSZ/inst/NEWS 2017-09-25 04:12:03 UTC (rev 211)
+++ pkg/CHNOSZ/inst/NEWS 2017-09-25 06:49:53 UTC (rev 212)
@@ -1,4 +1,4 @@
-CHANGES IN CHNOSZ 1.1.0-10 (2017-09-25)
+CHANGES IN CHNOSZ 1.1.0-11 (2017-09-25)
--------------------------------------
- water.lines() now works for diagrams of Eh, pe, logfO2, logaO2,
@@ -29,6 +29,10 @@
omega are now calculated within the function, thereby simplifying the
function call from subcrt().
+- Add 'H2O.props' argument to hkf(). This lists the properties of water
+ to calculate (at P, T) so that the calculations aren't duplicated by
+ subcrt().
+
CHANGES IN CHNOSZ 1.1.0 (2017-05-04)
------------------------------------
Modified: pkg/CHNOSZ/man/eos.Rd
===================================================================
--- pkg/CHNOSZ/man/eos.Rd 2017-09-25 04:12:03 UTC (rev 211)
+++ pkg/CHNOSZ/man/eos.Rd 2017-09-25 06:49:53 UTC (rev 212)
@@ -10,7 +10,7 @@
\usage{
cgl(property = NULL, T = 298.15, P = 1, parameters = NULL)
hkf(property = NULL, T = 298.15, P = 1, parameters = NULL,
- contrib = c("n", "s", "o"), H2O.PT = NULL)
+ contrib = c("n", "s", "o"), H2O.props = "rho")
}
\arguments{
@@ -19,7 +19,7 @@
\item{P}{numeric, pressure(s) at which to calculate properties (bar)}
\item{parameters}{dataframe, species parameters as one or more rows from \code{thermo$obigt}}
\item{contrib}{character, which contributions to consider in the revised HKF equations equations of state: (\code{n})onsolvation, (\code{s})olvation (the \eqn{\omega}{omega} terms), or (o)rigination contributions (i.e., the property itself at 25 \eqn{^{\circ}}{°}C and 1 bar). Default is \code{c("n","s","o")}, for all contributions}
- \item{H2O.PT}{dataframe, values of the electrostatic properties of water at the temperature(s) and pressure(s) of interest}
+ \item{H2O.props}{character, properties to calculate for water}
}
\details{
@@ -31,8 +31,7 @@
\code{hkf} implements the revised HKF equations of state (Helgeson et al., 1981; Tanger and Helgeson, 1988; Shock and Helgeson, 1988).
The equations-of-state parameters are \code{a1}, \code{a2}, \code{a3}, \code{a4}, \code{c1}, \code{c2}, \code{omega} and \code{Z}; the units of these parameters are as indicated for \code{\link{thermo}$obigt}, without the order of magnitude multipliers.
Note that the equation-of-state parameter \code{Z} (appearing in the \eqn{g}{g}-function for the temperature derivatives of the omega parameter; Shock et al., 1992) is taken from \code{thermo$obigt} and \emph{not} from the \code{\link{makeup}} of the species.
-\code{H2O.PT} is an optional argument that contains the electrostatic properties of \eqn{\mathrm{H_2O}}{H2O} required for the calculations.
-If either of these is \code{NULL} (the default), the required values are retrieved using \code{\link{water}}.
+\code{H2O.props} is an optional argument that lists the properties of water that should be returned; it is used by \code{\link{subcrt}} so that the time-consuming water calculations are only performed once.
The temperature and pressure derivatives of the \code{omega} parameter for charged species (where \code{Z != 0}, but not for the aqueous proton, H+) are calculated using the \eqn{g}{g}- and \eqn{f}{f}-functions described by Shock et al., 1992 and Johnson et al., 1992.
This option is turned off when the IAPWS-95 equations are activated (see \code{\link{water}}).
@@ -57,6 +56,7 @@
\value{
A list of length equal to the number of species (i.e., number rows of \code{parameters}).
Each element of the list contains a dataframe, each column of which corresponds to one of the specified properties; the number of rows is equal to the number of pressure-temperature points.
+Furthermore, in \code{hkf}, the output is a list consisting of the above-described object (named \samp{aq}) and a data frame of the calculated properties of water (named \samp{H2O}).
}
\seealso{
More information about the CHNOSZ-commits
mailing list