[CHNOSZ-commits] r255 - in pkg/CHNOSZ: . R demo inst man tests/testthat
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Oct 13 07:25:57 CEST 2017
Author: jedick
Date: 2017-10-13 07:25:56 +0200 (Fri, 13 Oct 2017)
New Revision: 255
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/R/nonideal.R
pkg/CHNOSZ/R/subcrt.R
pkg/CHNOSZ/demo/DEW.R
pkg/CHNOSZ/inst/NEWS
pkg/CHNOSZ/man/nonideal.Rd
pkg/CHNOSZ/tests/testthat/test-nonideal.R
Log:
subcrt(): recognize 'Helgeson'/'Helgeson0' methods for activity coefficients
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2017-10-12 14:56:08 UTC (rev 254)
+++ pkg/CHNOSZ/DESCRIPTION 2017-10-13 05:25:56 UTC (rev 255)
@@ -1,6 +1,6 @@
-Date: 2017-10-12
+Date: 2017-10-13
Package: CHNOSZ
-Version: 1.1.0-53
+Version: 1.1.0-54
Title: Thermodynamic Calculations for Geobiochemistry
Author: Jeffrey Dick
Maintainer: Jeffrey Dick <j3ffdick at gmail.com>
Modified: pkg/CHNOSZ/R/nonideal.R
===================================================================
--- pkg/CHNOSZ/R/nonideal.R 2017-10-12 14:56:08 UTC (rev 254)
+++ pkg/CHNOSZ/R/nonideal.R 2017-10-13 05:25:56 UTC (rev 255)
@@ -7,15 +7,21 @@
# generate nonideal contributions to thermodynamic properties
# number of species, same length as speciesprops list
# T in Kelvin, same length as nrows of speciespropss
- # arguments P, A_DH, B_DH are needed for Helgeson method only
+ # arguments P, A_DH, B_DH are needed for Helgeson/Helgeson0 methods only
+ mettext <- function(method) {
+ mettext <- paste(method, "method")
+ if(method=="Helgeson0") mettext <- "Helgeson method with B-dot = 0"
+ mettext
+ }
+
# we can use this function to change the nonideal method option
- if(identical(species, "Helgeson") | identical(species, "Alberty")) {
+ if(identical(species, "Helgeson") | identical(species, "Helgeson0") | identical(species, "Alberty")) {
thermo <- get("thermo")
oldnon <- thermo$opt$nonideal
thermo$opt$nonideal <- species
assign("thermo", thermo, "CHNOSZ")
- message("nonideal: setting method option to ", species)
+ message("nonideal: setting nonideal option to use ", mettext(species))
return(invisible(oldnon))
}
@@ -39,7 +45,8 @@
loggamma <- function(a, Z, I, B) - a * Z^2 * I^(1/2) / (1 + B * I^(1/2))
# TODO: check the following equations 20080208 jmd
R <- 1.9872 # gas constant, cal K^-1 mol^-1
- if(prop=="loggamma") return(loggamma(eval(A), Z, I, B))
+ # 20171013 convert loggamma to common logarithm
+ if(prop=="loggamma") return(loggamma(eval(A), Z, I, B) / log(10))
else if(prop=="G") return(R * T * loggamma(eval(A), Z, I, B))
else if(prop=="H") return(R * T^2 * loggamma(eval(DD(A, "T", 1)), Z, I, B))
else if(prop=="S") return(- R * T * loggamma(eval(DD(A, "T", 1)), Z, I, B))
@@ -50,10 +57,12 @@
Helgeson <- function(Z, I, T, P, A_DH, B_DH, prop = "loggamma") {
# "distance of closest approach" of ions in NaCl solutions (HKF81 Table 2)
acirc <- 3.72e-8 # cm
- loggamma <- - A_DH * Z^2 * I^0.5 / (1 + acirc * B_DH * I^0.5) + Bdot * I
+ if(method=="Helgeson") loggamma <- - A_DH * Z^2 * I^0.5 / (1 + acirc * B_DH * I^0.5) + Bdot * I
+ else if(method=="Helgeson0") loggamma <- - A_DH * Z^2 * I^0.5 / (1 + acirc * B_DH * I^0.5)
R <- 1.9872 # gas constant, cal K^-1 mol^-1
if(prop=="loggamma") return(loggamma)
- else if(prop=="G") return(R * T * loggamma)
+ else if(prop=="G") return(R * T * log(10) * loggamma)
+ # note the log(10) (=2.303) ... use natural logarithm to calculate G!!!
}
# get B-dot if we're using the Helgeson method
@@ -83,7 +92,7 @@
if(method=="Alberty" & pname %in% c("G", "H", "S", "Cp")) {
myprops[, j] <- myprops[, j] + Alberty(Z, IS, T, pname)
didit <- TRUE
- } else if(method=="Helgeson" & pname %in% c("G", "H", "S", "Cp")) {
+ } else if(grepl("Helgeson", method) & pname %in% c("G", "H", "S", "Cp")) {
myprops[, j] <- myprops[, j] + Helgeson(Z, IS, T, P, A_DH, B_DH, pname)
didit <- TRUE
}
@@ -91,14 +100,14 @@
# append a loggam column if we did any nonideal calculations of thermodynamic properties
if(didit) {
if(method=="Alberty") myprops <- cbind(myprops, loggam = Alberty(Z, IS, T))
- else if(method=="Helgeson") {
+ else if(grepl("Helgeson", method)) {
myprops <- cbind(myprops, loggam = Helgeson(Z, IS, T, P, A_DH, B_DH))
}
}
speciesprops[[i]] <- myprops
if(didit) ndid <- ndid + 1
}
- if(ndid > 0) message(paste("nonideal:", ndid, "species were nonideal"))
+ if(ndid > 0) message("nonideal: calculated activity coefficients for ", ndid, " species (", mettext(method), ")")
return(speciesprops)
}
Modified: pkg/CHNOSZ/R/subcrt.R
===================================================================
--- pkg/CHNOSZ/R/subcrt.R 2017-10-12 14:56:08 UTC (rev 254)
+++ pkg/CHNOSZ/R/subcrt.R 2017-10-13 05:25:56 UTC (rev 255)
@@ -56,10 +56,6 @@
# length checking
if(do.reaction & length(species)!=length(coeff))
stop('coeff must be same length as the number of species.')
- if(length(IS)>1) if(!identical(grid,'IS')) {
- if(is.null(grid)) grid <- 'IS'
- else stop('if you want length(IS) > 1, set grid=\'IS\'')
- }
if(!is.null(logact)) logact <- rep(logact,length.out=length(coeff))
# normalize temperature units
if(!missing(T)) {
@@ -100,6 +96,8 @@
# expansion of Psat and equivalence of argument lengths
tpargs <- TP.args(T=T,P=P)
T <- tpargs$T; P <- tpargs$P
+ if(length(newIS) > length(T)) T <- rep(T, length.out=length(newIS))
+ if(length(newIS) > length(P)) P <- rep(P, length.out=length(newIS))
}
# get species information
@@ -181,8 +179,6 @@
isH2O <- reaction$name=='water' & reaction$state=='liq'
isaq <- reaction$state=='aq'
- #if(length(T)==1) T.text <- paste(T,units('T')) else T.text <- paste(length(T),'values of T')
- #if(length(P)==1) P.text <- paste(P,units('P')) else P.text <- paste(length(P),'values of P')
ut <- T
if(identical(grid,'IS')) ut <- unique(ut)
if(length(ut)==1) T.text <- paste(ut,'K') else {
@@ -192,11 +188,10 @@
if(can.be.numeric(P)) P.text <- paste(round(as.numeric(P),2),'bar')
else P.text <- "P"
} else P.text <- 'P'
- #} else P.text <- paste(length(P),'values of P')
if(identical(P[[1]],'Psat')) P.text <- P
if(any(c(isH2O,isaq))) P.text <- paste(P.text,' (wet)',sep='')
if(length(species)==1 & convert==FALSE) {
- # we don't think we want messages here
+ # no message produced here
} else {
message(paste('subcrt:',length(species),'species at',T.text,'and',P.text))
}
@@ -280,12 +275,18 @@
si <- obigt2eos(thermo$obigt[inpho[isaq],], "aq", fixGHS = TRUE)
# always get density
H2O.props <- "rho"
+ # calculate A_DH and B_DH if we're using the B-dot (Helgeson) equation
+ if(any(IS != 0) & grepl("Helgeson", thermo$opt$nonideal)) H2O.props <- c(H2O.props, "A_DH", "B_DH")
# get other properties for H2O only if it's in the reaction
if(any(isH2O)) H2O.props <- c(H2O.props, eosprop)
hkfstuff <- hkf(eosprop, parameters = si, T = T, P = P, 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)
+ # calculate activity coefficients if ionic strength is not zero
+ if(any(IS != 0)) {
+ if(grepl("Helgeson", thermo$opt$nonideal)) p.aq <- nonideal(inpho[isaq], p.aq, newIS, T, P, H2O.PT$A_DH, H2O.PT$B_DH)
+ else if(thermo$opt$nonideal=="Alberty") 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
@@ -544,13 +545,9 @@
}
}
}
- # convert loggam to common logarithm and
# put ionic strength next to any loggam columns
for(i in 2:length(out)) {
- if('loggam' %in% colnames(out[[i]])) {
- out[[i]] <- cbind(out[[i]],IS=newIS)
- out[[i]][, "loggam"] <- out[[i]][, "loggam"]/log(10)
- }
+ if('loggam' %in% colnames(out[[i]])) out[[i]] <- cbind(out[[i]],IS=newIS)
}
# more fanagling for species
if(!do.reaction) {
Modified: pkg/CHNOSZ/demo/DEW.R
===================================================================
--- pkg/CHNOSZ/demo/DEW.R 2017-10-12 14:56:08 UTC (rev 254)
+++ pkg/CHNOSZ/demo/DEW.R 2017-10-13 05:25:56 UTC (rev 255)
@@ -139,7 +139,7 @@
T <- seq(600, 1000, 5)
-## use DEW model
+## activate DEW model
data(thermo)
water("DEW")
# add species data for DEW
@@ -167,12 +167,20 @@
## add species
species(c(inorganics, organics))
+## generate spline function for ionic strength
+IS <- splinefun(seq(600, 1000, 100), c(0.39, 0.57, 0.88, 1.45, 2.49))
+
+## use Debye-Huckel equation with B-dot set to zero
+nonideal("Helgeson0")
+## and calculate activity coefficient for the proton
+thermo$opt$ideal.H <<- FALSE
+
## calculate affinities on the T-logfO2 transect
-a <- affinity(T=T, O2=QFM_2, P=50000)
+a <- affinity(T=T, O2=QFM_2, P=50000, IS=IS(T))
## use the total carbon molality as an approximation of total activity
-molC <- splinefun(seq(600, 1000, 100), c(0.03, 0.2, 1, 4, 20))(T)
-loga.C <- log10(molC)
+molC <- splinefun(seq(600, 1000, 100), c(0.03, 0.2, 1, 4, 20))
+loga.C <- log10(molC(T))
## calculate metastable equilibrium activities
e <- equilibrate(a, loga.balance=loga.C)
@@ -193,19 +201,27 @@
t2 <- "after Sverjensky et al., 2014b"
mtitle(c(t1, t2))
-## additional checks
-# check that we're within 0.1 of the QFM-2 values used by SSH14
+### additional checks
+
+## check that we're within 0.1 of the QFM-2 values used by SSH14
stopifnot(maxdiff(QFM_2[T %% 100 == 0], c(-17.0, -14.5, -12.5, -10.8, -9.4)) < 0.1)
+
# Here are the logKs of aqueous species dissociation reactions at 600 degC and 50000 bar,
-# taken from the Supporting Information of the paper (p. 103-109):
+# values from EQ3NR output in Supporting Information of the paper (p. 103-109):
inorganic.logK <- c(24.4765, -9.0784, -5.3468, 0)
organic.logK <- c(1.7878, 2.5648, 15.3182, 16.9743, 30.4088, 28.9185)
# calculate equilibrium constants of the reactions in CHNOSZ; use a negative sign to change from formation to dissociation
logK.calc <- -unlist(affinity(T=600, P=50000, property="logK")$values)
logK.calc - c(inorganic.logK, organic.logK)
-# check that we're within 0.021 of the logK values used by SSH14
+## check that we're within 0.021 of the logK values used by SSH14
stopifnot(maxdiff(logK.calc, c(inorganic.logK, organic.logK)) < 0.021)
+## check that we get similar activity coefficients
+# values for monovalent species from EQ3NR output
+loggamma <- c(-0.15, -0.18, -0.22, -0.26, -0.31)
+sres <- subcrt("propanoate", T=seq(600, 1000, 100), P=50000, IS=c(0.39, 0.57, 0.88, 1.45, 2.49))
+stopifnot(maxdiff(sres$out[[1]]$loggam, loggamma) < 0.004)
+
###########
### all done!
# reset the database and previous water computational option
Modified: pkg/CHNOSZ/inst/NEWS
===================================================================
--- pkg/CHNOSZ/inst/NEWS 2017-10-12 14:56:08 UTC (rev 254)
+++ pkg/CHNOSZ/inst/NEWS 2017-10-13 05:25:56 UTC (rev 255)
@@ -1,4 +1,4 @@
-CHANGES IN CHNOSZ 1.1.0-53 (2017-10-12)
+CHANGES IN CHNOSZ 1.1.0-54 (2017-10-13)
---------------------------------------
MAJOR CHANGES:
@@ -18,8 +18,8 @@
using water("DEW"), water("IAPWS"), etc.
- Usage of the DEW model is shown in the new demo DEW.R. This demo also
- depends on the Berman equations (above) and, for the last diagram, the
- following *two* NEWS items:
+ depends on the Berman equations (above) and, for the last diagram in
+ the demo, the following *four* NEWS items:
- In equilibrate(), it is now possible to combine affinity calculations
with changing activity of the balancing basis species (loga.balance).
@@ -35,10 +35,11 @@
"percent carbon" plots for systems where the species have different
carbon numbers.
-- nonideal() now has two methods for calculating activity coefficients:
- 'Helgeson', which utilizes the "B-dot" equation of Helgeson, 1969,
- and 'Alberty' (the only method previously available). The 'Helgeson'
- method depends on the following *two* NEWS items:
+- nonideal() now has three methods for calculating activity coefficients:
+ 'Helgeson' and 'Helgeson0', which utilize the "B-dot" equation of
+ Helgeson, 1969, and 'Alberty' (the only method previously available).
+ The 'Helgeson' method depends on the following *two* NEWS items (but
+ the 'Helgeson0' method omits the B-dot term):
- All three water options (SUPCRT92, IAPWS95, DEW) can be used to
calculate 'A_DH' and 'B_DH', i.e. A and B coefficients in the extended
@@ -46,9 +47,9 @@
- Add Bdot() to calculate the "B-dot" (or b_gamma) extended term
parameter for activity coefficients in NaCl solutions at high
- temperature and pressure (Helgeson et al., 1981) including
- high-pressure extrapolations > 10 kbar based on data from Manning et
- al., 2013.
+ temperature and pressure (Helgeson et al., 1981), including
+ very high pressures (> 10 kbar) based on extrapolations from Manning
+ et al., 2013.
- For minerals with phase transitions (states 'cr2' 'cr3' etc.) in
thermo$obigt (i.e. the Helgeson minerals), it is now possible to use
@@ -60,7 +61,7 @@
coesite. Calculations for other minerals still assume constant
volume of each phase.
-- calculations of the g function are now enabled for DEW (with
+- Calculations of the g function are now enabled for DEW (with
pressure derivative) and IAPWS-95 (no derivatives included).
- water.lines() now works for diagrams of Eh, pe, logfO2, logaO2,
Modified: pkg/CHNOSZ/man/nonideal.Rd
===================================================================
--- pkg/CHNOSZ/man/nonideal.Rd 2017-10-12 14:56:08 UTC (rev 254)
+++ pkg/CHNOSZ/man/nonideal.Rd 2017-10-13 05:25:56 UTC (rev 255)
@@ -20,7 +20,7 @@
\item{P}{numeric, pressure (bar); required for Helgeson method}
\item{A_DH}{numeric, A Debye-Huckel coefficient; required for Helgeson method}
\item{B_DH}{numeric, B Debye-Huckel coefficient; required for Helgeson method}
- \item{method}{character, \samp{Alberty} or \samp{Helgeson}}
+ \item{method}{character, \samp{Alberty}, \samp{Helgeson}, or \samp{Helgeson0}}
\item{TC}{numeric, temperature (\degC)}
\item{showsplines}{character, show isobaric (\samp{T}) or isothermal (\samp{P}) splines}
}
@@ -35,9 +35,11 @@
If \code{method} is \samp{Alberty}, the values of \code{IS} are combined with Alberty's (2003) equation 3.6-1 (extended Debye-Hückel equation) and its derivatives, to calculate apparent molal properties at the specified ionic strength(s) and temperature(s).
The apparent molal properties that can be calculated include \samp{G}, \samp{H}, \samp{S} and \samp{Cp}; any columns in the dataframes of \code{speciesprops} with other names are left untouched.
-If \code{method} is \samp{Helgeson}, the \dQuote{B-dot} equation (Helgeson, 1969; Helgeson et al., 1981; Manning, 2013) is used.
+If \code{method} is \samp{Helgeson}, the \dQuote{B-dot} equation is used.
+This equation seems to have been originally proposed by Huckel, 1925; parameters were derived for use at high temperature and pressure by Helgeson, 1969; Helgeson et al., 1981; Manning, 2013.
In addition to \code{IS} and \code{T}, this method depends on values of \code{P}, \code{A_DH}, and \code{B_DH} given in the arguments.
The calculation of \dQuote{B-dot}, also used in the equations, is made within \code{nonideal} by calling the \code{Bdot} function.
+For some uses, it is desirable to set the \dQuote{B-dot} parameter to zero; this can be done by setting the method to \code{Helgeson0}.
Currently, \samp{G} is the only apparent molal property that is calculated (but this can be used by \code{\link{subcrt}} to calculate apparent equilibrium constants).
\code{Bdot} calculates the \dQuote{B-dot} deviation function (Helgeson, 1969) a.k.a. extended term parameter (written as b_gamma; Helgeson et al., 1981) for activity coefficients in NaCl solutions at high temperature and pressure.
@@ -49,15 +51,9 @@
This is a crude method of kriging the data, but produces fairly smooth interpolations without adding any external dependencies.
}
-\section{Warning}{
-The logarithms of activity coefficients (\code{loggam}) returned by \code{nonideal} use the natural logarithm (cf. Alberty, 2003 Eq. 3.6-1).
-To maintain consistency with the conventions used elsewhere in the package (i.e. for logarithms of equilibrium constants and of chemical activities), the values of \code{loggam} returned by \code{\link{subcrt}} are expressed using the common (base 10) logarithm.
-Note that the first example below uses \code{loggam} returned by \code{subcrt}, therefore requiring a base of 10 for calculating gamma.
-}
-
\value{
One (\samp{G}) or more (\samp{H}, \samp{S}, \samp{Cp}; currently only with the Alberty method) standard thermodynamic properties (at IS=0) in \code{speciesprops} are replaced by the corresponding apparent thermodynamic properties (at higher IS).
-For all affected species, a column named \code{loggam} (logarithm of gamma, the activity coefficient) is appended to the output dataframe of species properties.
+For all affected species, a column named \code{loggam} (common (base-10) logarithm of gamma, the activity coefficient) is appended to the output dataframe of species properties.
}
\examples{\dontshow{data(thermo)}
@@ -82,7 +78,7 @@
speciesprops <- subcrt(species, T = rep(T[j], length(IS)))$out
# use nonideal to calculate loggamma; this also adjusts G, H, S, Cp, but we don't use them here
nonidealprops <- nonideal(species, speciesprops, IS = IS, T = convert(T[j], "K"))
- for(i in 1:4) lines(IS, exp(nonidealprops[[i]]$loggam), lty=lty[j], col=col[i])
+ for(i in 1:4) lines(IS, 10^(nonidealprops[[i]]$loggam), lty=lty[j], col=col[i])
}
t1 <- "activity coefficient (gamma) of -1,-2,-3,-4 charged species"
t2 <- quote("at 0, 25, and 40 "*degree*"C, after Alberty, 2003")
@@ -170,6 +166,8 @@
Helgeson, H. C., Kirkham, D. H. and Flowers, G. C. (1981) Theoretical prediction of the thermodynamic behavior of aqueous electrolytes at high pressures and temperatures. IV. Calculation of activity coefficients, osmotic coefficients, and apparent molal and standard and relative partial molal properties to 600\degC and 5 Kb. \emph{Am. J. Sci.} \bold{281}, 1249--1516. \url{https://doi.org/10.2475/ajs.281.10.1249}
+Hückel, E. (1925). The theory of concentrated, aqueous solutions of strong electrolytes. \emph{Physikalische Zeitschrift} \bold{26}, 93--147.
+
Manning, C. E. (2013) Thermodynamic modeling of fluid-rock interaction at mid-crustal to upper-mantle conditions. \emph{Rev. Mineral. Geochem.} \bold{76}, 135--164. \url{https://doi.org/10.2138/rmg.2013.76.5}
Manning, C. E., Shock, E. L. and Sverjensky, D. A. (2013) The chemistry of carbon in aqueous fluids at crustal and upper-mantle conditions: Experimental and theoretical constraints. \emph{Rev. Mineral. Geochem.} \bold{75}, 109--148. \url{https://doi.org/10.2138/rmg.2013.75.5}
Modified: pkg/CHNOSZ/tests/testthat/test-nonideal.R
===================================================================
--- pkg/CHNOSZ/tests/testthat/test-nonideal.R 2017-10-12 14:56:08 UTC (rev 254)
+++ pkg/CHNOSZ/tests/testthat/test-nonideal.R 2017-10-13 05:25:56 UTC (rev 255)
@@ -49,3 +49,25 @@
expect_maxdiff(DEW30$A_DH, A30, 0.06)
expect_maxdiff(DEW30$B_DH / 1e8, B30, 0.024)
})
+
+#test_that("different methods give correct values of loggamma", {
+#
+#})
+
+test_that("affinity transect incorporates IS correctly", {
+ basis("CHNOS+")
+ species("acetate")
+ # calculations at single combinations of logfO2 and IS
+ basis("O2", -80); a80_0 <- affinity()
+ basis("O2", -60); a60_1 <- affinity(IS=1)
+ # calculations on a transect with those endpoints
+ a <- affinity(O2=seq(-80, -60, length.out=4), IS=seq(0, 1, length.out=4))
+ expect_equal(a$values[[1]][1], a80_0$values[[1]][1])
+ expect_equal(a$values[[1]][4], a60_1$values[[1]][1])
+ # 20171013: that was working fine, but how about a more complicated case involving T?
+ a25_0 <- affinity()
+ a50_1 <- affinity(T=50, IS=1)
+ a <- affinity(T=seq(25, 50, length.out=4), IS=seq(0, 1, length.out=4))
+ expect_equal(a$values[[1]][1], a25_0$values[[1]][1])
+ expect_equal(a$values[[1]][4], a50_1$values[[1]][1])
+})
More information about the CHNOSZ-commits
mailing list