[CHNOSZ-commits] r235 - in pkg/CHNOSZ: . R man tests/testthat
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Oct 2 08:22:50 CEST 2017
Author: jedick
Date: 2017-10-02 08:22:50 +0200 (Mon, 02 Oct 2017)
New Revision: 235
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/R/berman.R
pkg/CHNOSZ/R/cgl.R
pkg/CHNOSZ/R/subcrt.R
pkg/CHNOSZ/man/berman.Rd
pkg/CHNOSZ/man/extdata.Rd
pkg/CHNOSZ/tests/testthat/test-berman.R
Log:
make berman() usable by subcrt()
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2017-10-02 04:45:01 UTC (rev 234)
+++ pkg/CHNOSZ/DESCRIPTION 2017-10-02 06:22:50 UTC (rev 235)
@@ -1,6 +1,6 @@
Date: 2017-10-02
Package: CHNOSZ
-Version: 1.1.0-33
+Version: 1.1.0-34
Title: Thermodynamic Calculations for Geobiochemistry
Author: Jeffrey Dick
Maintainer: Jeffrey Dick <j3ffdick at gmail.com>
Modified: pkg/CHNOSZ/R/berman.R
===================================================================
--- pkg/CHNOSZ/R/berman.R 2017-10-02 04:45:01 UTC (rev 234)
+++ pkg/CHNOSZ/R/berman.R 2017-10-02 06:22:50 UTC (rev 235)
@@ -1,7 +1,7 @@
# CHNOSZ/berman.R 20170930
# calculate thermodynamic properties of minerals using Berman formulation
-berman <- function(name, T = 298.15, P = 1) {
+berman <- function(name, T = 298.15, P = 1, thisinfo=NULL, check.G=FALSE) {
# reference temperature and pressure
Pr <- 1
Tr <- 298.15
@@ -23,15 +23,17 @@
d0 <- d1 <- d2 <- d3 <- d4 <- d5 <- k0 <- k1 <- k2 <- k3 <- v1 <- v2 <- v3 <- v4 <- NA
# assign values to the variables used below
for(i in 1:ncol(dat)) assign(colnames(dat)[i], dat[irow, i])
+ # get the entropy of the elements using the chemical formula in thermo$obigt
+ if(is.null(thisinfo)) thisinfo <- info(info(name, "cr_Berman", check.it=FALSE))
+ SPrTr_elements <- convert(entropy(thisinfo$formula), "J")
# check that G in data file is the G of formation from the elements --> Benson-Helgeson convention (DG = DH - T*DS)
- # we get the entropy of the elements using the chemical formula in thermo$obigt
- iname <- info(name, "cr_Berman", check.it=FALSE)
- SPrTr_elements <- convert(entropy(info(iname)$formula), "J")
- GfPrTr_calc <- HfPrTr - Tr * (SPrTr - SPrTr_elements)
- Gdiff <- GfPrTr_calc - GfPrTr
- if(abs(Gdiff) >= 1000) warning(paste0(name, ": GfPrTr(calc) - GfPrTr(table) is too big! == ",
- round(GfPrTr_calc - GfPrTr), " J/mol"), call.=FALSE)
- # (the tabulated GfPrTr is unused below)
+ if(check.G) {
+ GfPrTr_calc <- HfPrTr - Tr * (SPrTr - SPrTr_elements)
+ Gdiff <- GfPrTr_calc - GfPrTr
+ if(abs(Gdiff) >= 1000) warning(paste0(name, ": GfPrTr(calc) - GfPrTr(table) is too big! == ",
+ round(GfPrTr_calc - GfPrTr), " J/mol"), call.=FALSE)
+ # (the tabulated GfPrTr is unused below)
+ }
### thermodynamic properties ###
# calculate Cp and V (Berman, 1988 Eqs. 4 and 5)
Modified: pkg/CHNOSZ/R/cgl.R
===================================================================
--- pkg/CHNOSZ/R/cgl.R 2017-10-02 04:45:01 UTC (rev 234)
+++ pkg/CHNOSZ/R/cgl.R 2017-10-02 06:22:50 UTC (rev 235)
@@ -14,88 +14,95 @@
for(k in 1:nrow(parameters)) {
# the parameters for *this* species
PAR <- parameters[k, ]
- # start with NA values
- values <- data.frame(matrix(NA, ncol = length(property), nrow=ncond))
- colnames(values) <- property
- # additional calculations for quartz and coesite
- qtz <- quartz_coesite(PAR, T, P)
- isqtz <- !identical(qtz$V, 0)
- for(i in 1:length(property)) {
- PROP <- property[i]
- # a test for availability of the EoS parameters
- # here we assume that the parameters are in the same columns as in thermo$obigt
- # leave T transition (in 20th column) alone
- hasEOS <- any(!is.na(PAR[, 13:19]))
- # if at least one of the EoS parameters is available, zero out any NA's in the rest
- if(hasEOS) PAR[, 13:19][, is.na(PAR[, 13:19])] <- 0
- # equations for lambda adapted from HOK+98
- if(PROP == "Cp") {
- # use constant Cp if the EoS parameters are not available
- if(!hasEOS) p <- PAR$Cp
- else p <- PAR$a + PAR$b * T + PAR$c * T^-2 + PAR$d * T^-0.5 + PAR$e * T^2 + PAR$f * T^PAR$lambda
- }
- if(PROP == "V") {
- if(isqtz) p <- qtz$V
- else p <- rep(PAR$V, ncond)
- }
- if(PROP %in% c("E", "kT")) {
- p <- rep(NA, ncond)
- warning("cgl: E and/or kT of cr, gas and/or liq species are NA.")
- }
- if(PROP == "G") {
- # use constant Cp if the EoS parameters are not available
- if(!hasEOS) p <- PAR$Cp * (T - Tr - T * log(T/Tr)) else {
- # Gibbs energy integral: the value at Tref plus heat capacity terms
- p <- PAR$a * (T - Tr - T * log(T/Tr)) -
- PAR$b * (T - Tr)^2 / 2 - PAR$c * (1/T + T/Tr^2 - 2/Tr) / 2 -
- PAR$d * (T^0.5 - 0.5 * T * Tr^-0.5 - 0.5 * Tr^0.5) / -0.25 -
- PAR$e * (T^3 - 3 * T * Tr^2 + 2 * Tr^3) / 6
+ if(PAR$state=="cr_Berman") {
+ # use Berman equations (parameters not in thermo$obigt)
+ properties <- berman(PAR$name, T=T, P=P, thisinfo=PAR)
+ iprop <- match(property, colnames(properties))
+ values <- properties[, iprop, drop=FALSE]
+ } else {
+ # start with NA values
+ values <- data.frame(matrix(NA, ncol = length(property), nrow=ncond))
+ colnames(values) <- property
+ # additional calculations for quartz and coesite
+ qtz <- quartz_coesite(PAR, T, P)
+ isqtz <- !identical(qtz$V, 0)
+ for(i in 1:length(property)) {
+ PROP <- property[i]
+ # a test for availability of the EoS parameters
+ # here we assume that the parameters are in the same columns as in thermo$obigt
+ # leave T transition (in 20th column) alone
+ hasEOS <- any(!is.na(PAR[, 13:19]))
+ # if at least one of the EoS parameters is available, zero out any NA's in the rest
+ if(hasEOS) PAR[, 13:19][, is.na(PAR[, 13:19])] <- 0
+ # equations for lambda adapted from HOK+98
+ if(PROP == "Cp") {
+ # use constant Cp if the EoS parameters are not available
+ if(!hasEOS) p <- PAR$Cp
+ else p <- PAR$a + PAR$b * T + PAR$c * T^-2 + PAR$d * T^-0.5 + PAR$e * T^2 + PAR$f * T^PAR$lambda
}
- # use additional heat capacity term if it's defined
- if(!is.na(PAR$f) & !is.na(PAR$lambda)) if(PAR$f != 0) {
- if(PAR$lambda == -1) p <- p + PAR$f * (log(T/Tr) - T * (1/Tr - 1/T))
- else p <- p + PAR$f * ( T^(PAR$lambda + 1) - (PAR$lambda + 1) * T * Tr^PAR$lambda +
- PAR$lambda * Tr^(PAR$lambda + 1) ) / ( PAR$lambda * (PAR$lambda + 1) )
+ if(PROP == "V") {
+ if(isqtz) p <- qtz$V
+ else p <- rep(PAR$V, ncond)
}
- # entropy and volume terms
- if(!is.na(PAR$S)) p <- p - PAR$S * (T - Tr)
- if(isqtz) p <- p + qtz$G
- else if(!is.na(PAR$V)) p <- p + convert(PAR$V * (P - Pr), "calories")
- p <- PAR$G + p
- }
- if(PROP == "H") {
- # use constant Cp if the EoS parameters are not available
- if(!hasEOS) p <- PAR$Cp * (T - Tr) else {
- p <- PAR$a * (T - Tr) + PAR$b * (T^2 - Tr^2) / 2 +
- PAR$c * (1/T - 1/Tr) / -1 + PAR$d * (T^0.5 - Tr^0.5) / 0.5 +
- PAR$e * (T^3 - Tr^3) / 3
+ if(PROP %in% c("E", "kT")) {
+ p <- rep(NA, ncond)
+ warning("cgl: E and/or kT of cr, gas and/or liq species are NA.")
}
- if(!is.na(PAR$f) & !is.na(PAR$lambda)) if(PAR$f != 0) {
- if(PAR$lambda == -1) p <- p + PAR$f * log(T/Tr)
- else p <- p - PAR$f * ( T^(PAR$lambda + 1) - Tr^(PAR$lambda + 1) ) / (PAR$lambda + 1)
+ if(PROP == "G") {
+ # use constant Cp if the EoS parameters are not available
+ if(!hasEOS) p <- PAR$Cp * (T - Tr - T * log(T/Tr)) else {
+ # Gibbs energy integral: the value at Tref plus heat capacity terms
+ p <- PAR$a * (T - Tr - T * log(T/Tr)) -
+ PAR$b * (T - Tr)^2 / 2 - PAR$c * (1/T + T/Tr^2 - 2/Tr) / 2 -
+ PAR$d * (T^0.5 - 0.5 * T * Tr^-0.5 - 0.5 * Tr^0.5) / -0.25 -
+ PAR$e * (T^3 - 3 * T * Tr^2 + 2 * Tr^3) / 6
+ }
+ # use additional heat capacity term if it's defined
+ if(!is.na(PAR$f) & !is.na(PAR$lambda)) if(PAR$f != 0) {
+ if(PAR$lambda == -1) p <- p + PAR$f * (log(T/Tr) - T * (1/Tr - 1/T))
+ else p <- p + PAR$f * ( T^(PAR$lambda + 1) - (PAR$lambda + 1) * T * Tr^PAR$lambda +
+ PAR$lambda * Tr^(PAR$lambda + 1) ) / ( PAR$lambda * (PAR$lambda + 1) )
+ }
+ # entropy and volume terms
+ if(!is.na(PAR$S)) p <- p - PAR$S * (T - Tr)
+ if(isqtz) p <- p + qtz$G
+ else if(!is.na(PAR$V)) p <- p + convert(PAR$V * (P - Pr), "calories")
+ p <- PAR$G + p
}
- if(isqtz) p <- p + qtz$H
- ## SUPCRT seems to ignore this term? ... 20070802
- #else p <- p + convert(PAR$V*(P-Pr),'calories')
- p <- PAR$H + p
- }
- if(PROP=="S") {
- # use constant Cp if the EoS parameters are not available
- if(!hasEOS) p <- PAR$Cp * log(T/Tr) else {
- p <- PAR$a * log(T / Tr) + PAR$b * (T - Tr) +
- PAR$c * (T^-2 - Tr^-2) / -2 + PAR$e * (T^2 - Tr^2) / 2 +
- PAR$d * (T^-0.5 - Tr^-0.5) / -0.5
+ if(PROP == "H") {
+ # use constant Cp if the EoS parameters are not available
+ if(!hasEOS) p <- PAR$Cp * (T - Tr) else {
+ p <- PAR$a * (T - Tr) + PAR$b * (T^2 - Tr^2) / 2 +
+ PAR$c * (1/T - 1/Tr) / -1 + PAR$d * (T^0.5 - Tr^0.5) / 0.5 +
+ PAR$e * (T^3 - Tr^3) / 3
+ }
+ if(!is.na(PAR$f) & !is.na(PAR$lambda)) if(PAR$f != 0) {
+ if(PAR$lambda == -1) p <- p + PAR$f * log(T/Tr)
+ else p <- p - PAR$f * ( T^(PAR$lambda + 1) - Tr^(PAR$lambda + 1) ) / (PAR$lambda + 1)
+ }
+ if(isqtz) p <- p + qtz$H
+ ## SUPCRT seems to ignore this term? ... 20070802
+ #else p <- p + convert(PAR$V*(P-Pr),'calories')
+ p <- PAR$H + p
}
- if(!is.na(PAR$f) & !is.na(PAR$lambda)) if(PAR$f != 0) {
- p <- p + PAR$f * (T^PAR$lambda - Tr^PAR$lambda) / PAR$lambda
+ if(PROP=="S") {
+ # use constant Cp if the EoS parameters are not available
+ if(!hasEOS) p <- PAR$Cp * log(T/Tr) else {
+ p <- PAR$a * log(T / Tr) + PAR$b * (T - Tr) +
+ PAR$c * (T^-2 - Tr^-2) / -2 + PAR$e * (T^2 - Tr^2) / 2 +
+ PAR$d * (T^-0.5 - Tr^-0.5) / -0.5
+ }
+ if(!is.na(PAR$f) & !is.na(PAR$lambda)) if(PAR$f != 0) {
+ p <- p + PAR$f * (T^PAR$lambda - Tr^PAR$lambda) / PAR$lambda
+ }
+ p <- PAR$S + p + qtz$S
}
- p <- PAR$S + p + qtz$S
+ values[, i] <- p
}
- values[, i] <- p
- }
- out[[k]] <- values
- }
- return(out)
+ } # end calculations using parameters from thermo$obigt
+ out[[k]] <- values
+ } # end loop over species
+ return(out)
}
### unexported function ###
Modified: pkg/CHNOSZ/R/subcrt.R
===================================================================
--- pkg/CHNOSZ/R/subcrt.R 2017-10-02 04:45:01 UTC (rev 234)
+++ pkg/CHNOSZ/R/subcrt.R 2017-10-02 06:22:50 UTC (rev 235)
@@ -293,8 +293,8 @@
}
# crystalline, gas, liquid (except water) species
- iscgl <- reaction$state %in% c('liq','cr','gas','cr1','cr2','cr3',
- 'cr4','cr5','cr6','cr7','cr8','cr9') & reaction$name != 'water'
+ cglstates <- c("liq", "cr", "gas", "cr1", "cr2", "cr3", "cr4", "cr5", "cr6", "cr7", "cr8", "cr9", "cr_Berman")
+ iscgl <- reaction$state %in% cglstates & reaction$name != "water"
if(TRUE %in% iscgl) {
#si <- info(inpho[iscgl],quiet=TRUE)
@@ -314,9 +314,12 @@
# name and state
myname <- reaction$name[i]
mystate <- reaction$state[i]
+ # don't proceed if the state is cr_Berman
+ if(mystate=="cr_Berman") next
# check if we're below the transition temperature
if(!(reaction$state[i] %in% c('cr1','liq','cr','gas'))) {
Ttr <- Ttr(inpho[i]-1,P=P,dPdT=dPdTtr(inpho[i]-1))
+ if(all(is.na(Ttr))) next
if(any(T < Ttr)) {
status.Ttr <- "(extrapolating G)"
if(!exceed.Ttr) {
@@ -332,10 +335,9 @@
Ttr <- Ttr(inpho[i],P=P,dPdT=dPdTtr(inpho[i]))
else {
Ttr <- thermo$obigt$z.T[inpho[i]]
- if(is.na(Ttr)) next
}
- if(all(Ttr==0)) next
- if(any(T >= Ttr)) {
+ if(all(is.na(Ttr))) next
+ if(!all(Ttr==0) & any(T >= Ttr)) {
status.Ttr <- "(extrapolating G)"
if(!exceed.Ttr) {
p.cgl[[ncgl[i]]]$G[T>=Ttr] <- NA
@@ -357,14 +359,13 @@
# use variable-pressure standard Gibbs energy for gases
isgas <- reaction$state %in% "gas"
- if(any(isgas) & "g" %in% eosprop & 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% calcprop) {
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)))
colnames(out[[i]][ncol(out[[i]])]) <- 'logK'
}
Modified: pkg/CHNOSZ/man/berman.Rd
===================================================================
--- pkg/CHNOSZ/man/berman.Rd 2017-10-02 04:45:01 UTC (rev 234)
+++ pkg/CHNOSZ/man/berman.Rd 2017-10-02 06:22:50 UTC (rev 235)
@@ -7,13 +7,15 @@
}
\usage{
- berman(name, T = 298.15, P = 1)
+ berman(name, T = 298.15, P = 1, thisinfo = NULL, check.G = FALSE)
}
\arguments{
\item{name}{character, name of mineral}
\item{T}{numeric, temperature(s) at which to calculate properties (K)}
\item{P}{numeric, pressure(s) at which to calculate properties (bar)}
+ \item{thisinfo}{dataframe, row for mineral from \code{thermo$obigt}}
+ \item{check.G}{logical, check consistency of G, H, and S?}
}
\details{
@@ -23,6 +25,11 @@
These entropies are used to convert the apparent Gibbs energies from the Berman-Brown convention to the the Benson-Helgeson convention.
Becuase they use a different set of parameters than Helgeson et al., 1978 (see \code{\link{cgl}}), the standard state thermodynamic properties and parameters for the calculations are stored in files under \code{extdata/Berman}.
+
+If \code{check.G} is TRUE, the tabulated value of DGfTrPr is compared with one calculated from DHfPrTr - T*DSPrTr (DS is the difference between the summed entropies of the elements and the tabulated entropy for the mineral).
+A warning is produced if the absolute value of the difference between tabulated and calculated DGfTrPr is greater than 1000 J/mol.
+
+Providing \code{thisinfo} means that the mineral name is not searched in \code{thermo$obigt}, potentially saving some running time.
}
\examples{
@@ -34,13 +41,25 @@
berman("quartz")
# Gibbs energies of aQz and coesite at higher T and P
T <- seq(200, 1300, 100)
-P <- seq(23000, 32000, length.out=length(T))
+P <- seq(22870, 31900, length.out=length(T))
G_aQz <- berman("quartz", T=T, P=P)$G
G_Cs <- berman("coesite", T=T, P=P)$G
# that is close to the univariant curve (Ber88 Fig. 4),
# so the difference in G is close to 0
DGrxn <- G_Cs - G_aQz
stopifnot(all(abs(DGrxn) < 100))
+
+# calculate the properties of a "reaction" between
+# the Helgeson and Berman versions of quartz
+T <- 1000
+P <- c(10000, 20000)
+subcrt(rep("quartz", 2), c("cr", "cr_Berman"), c(-1, 1), T=T, P=P)
+
+# make a P-T diagram for SiO2 minerals (Ber88 Fig. 4)
+basis(c("SiO2", "O2"), c("cr_Berman", "gas"))
+species(c("quartz", "quartz,beta", "coesite"), "cr_Berman")
+a <- affinity(T=c(200, 1700, 200), P=c(0, 50000, 200))
+diagram(a)
}
\references{
Modified: pkg/CHNOSZ/man/extdata.Rd
===================================================================
--- pkg/CHNOSZ/man/extdata.Rd 2017-10-02 04:45:01 UTC (rev 234)
+++ pkg/CHNOSZ/man/extdata.Rd 2017-10-02 06:22:50 UTC (rev 235)
@@ -41,7 +41,7 @@
\item \code{SOJSH.csv} Experimental equilibrium constants for the reaction NaCl(aq) = Na+ + Cl- as a function of temperature and pressure taken from Fig. 1 of Shock et al., 1992. Data were extracted from the figure using g3data (\url{http://www.frantz.fi/software/g3data.php}). See \code{demo("NaCl")} for an example that uses this file.
\item \code{Cp.CH4.HW97.csv}, \code{V.CH4.HWM96.csv} Apparent molar heat capacities and volumes of CH4 in dilute aqueous solutions reported by Hnědkovský and Wood, 1997 and Hnědkovský et al., 1996. See \code{\link{EOSregress}} and the vignette \code{eos-regress.Rmd} for examples that use these files.
\item \code{SC10_Rainbow.csv} Values of temperature (\eqn{^{\circ}}{°}C), pH and logarithms of activity of \CO2, \H2, \NH4plus, \H2S and \CH4 for mixing of seawater and hydrothermal fluid at Rainbow field (Mid-Atlantic Ridge), taken from Shock and Canovas, 2010. See the vignette \code{anintro.Rmd} for an example that uses this file.
- \item \code{SS98_Fig5a.csv}, \code{SS98_Fig5b.csv} Values of logarithm of fugacity of \eqn{\mathrm{O_2}}{O2} and pH as a function of temperature for mixing of seawater and hydrothermal fluid, digitized from Figs. 5a and b of Shock and Schulte, 1998. See the vignette \code{anintro.Rmd} for an example that uses this file.
+ \item \code{SS98_Fig5a.csv}, \code{SS98_Fig5b.csv} Values of logarithm of fugacity of \O2 and pH as a function of temperature for mixing of seawater and hydrothermal fluid, digitized from Figs. 5a and b of Shock and Schulte, 1998. See the vignette \code{anintro.Rmd} for an example that uses this file.
\item \code{rubisco.csv} UniProt IDs for Rubisco, ranges of optimal growth temperature of organisms, domain and name of organisms, and URL of reference for growth temperature, from Dick, 2014. See the vignette \code{anintro.Rmd} for an example that uses this file.
}
Modified: pkg/CHNOSZ/tests/testthat/test-berman.R
===================================================================
--- pkg/CHNOSZ/tests/testthat/test-berman.R 2017-10-02 04:45:01 UTC (rev 234)
+++ pkg/CHNOSZ/tests/testthat/test-berman.R 2017-10-02 06:22:50 UTC (rev 235)
@@ -11,7 +11,7 @@
# running this without error means that:
# - formulas for the minerals are found in thermo$obigt
# - there are no warnings for minerals with GfPrTr(calc) >= 1000 J/cal different from GfPrTr(table)
- expect_silent(properties <- lapply(mineral, berman))
+ expect_silent(properties <- lapply(mineral, berman, check.G=TRUE))
# save the results so we can use them in the next tests
assign("prop_Berman", properties, inherits=TRUE)
More information about the CHNOSZ-commits
mailing list