[CHNOSZ-commits] r784 - in pkg/CHNOSZ: . R demo inst inst/tinytest man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon May 15 12:21:42 CEST 2023
Author: jedick
Date: 2023-05-15 12:21:42 +0200 (Mon, 15 May 2023)
New Revision: 784
Added:
pkg/CHNOSZ/R/EOSregress.R
pkg/CHNOSZ/demo/adenine.R
pkg/CHNOSZ/inst/tinytest/test-EOSregress.R
pkg/CHNOSZ/man/EOSregress.Rd
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/NAMESPACE
pkg/CHNOSZ/R/examples.R
pkg/CHNOSZ/demo/00Index
pkg/CHNOSZ/inst/NEWS.Rd
pkg/CHNOSZ/man/CHNOSZ-package.Rd
pkg/CHNOSZ/man/examples.Rd
pkg/CHNOSZ/man/extdata.Rd
Log:
Restore EOSregress.R and demo/adenine.R
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2023-05-15 09:04:46 UTC (rev 783)
+++ pkg/CHNOSZ/DESCRIPTION 2023-05-15 10:21:42 UTC (rev 784)
@@ -1,6 +1,6 @@
Date: 2023-05-15
Package: CHNOSZ
-Version: 2.0.0-3
+Version: 2.0.0-4
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/NAMESPACE
===================================================================
--- pkg/CHNOSZ/NAMESPACE 2023-05-15 09:04:46 UTC (rev 783)
+++ pkg/CHNOSZ/NAMESPACE 2023-05-15 10:21:42 UTC (rev 784)
@@ -27,6 +27,7 @@
"equil.boltzmann", "equil.reaction", "find.tp",
"ionize.aa", "MP90.cp", "aasum",
"ratlab",
+ "EOSregress", "EOScoeffs", "EOSplot", "EOSvar",
# demos
"protein.equil", "palply",
"label.plot",
@@ -39,6 +40,8 @@
# (no other functions are used in the tests)
# other exported functions that are not used above
"checkEOS", "checkGHS", "check.OBIGT",
+ "V_s_var", "Cp_s_var",
+ "EOSlab", "EOScalc",
"basis.elements", "element.mu", "ibasis",
"water.SUPCRT92",
"nonideal",
Copied: pkg/CHNOSZ/R/EOSregress.R (from rev 715, pkg/CHNOSZ/R/EOSregress.R)
===================================================================
--- pkg/CHNOSZ/R/EOSregress.R (rev 0)
+++ pkg/CHNOSZ/R/EOSregress.R 2023-05-15 10:21:42 UTC (rev 784)
@@ -0,0 +1,212 @@
+# CHNOSZ/EOSregress.R
+# Model volumes and heat capacities of aqueous species
+# 20091105 first version
+# 20110429 revise and merge with CHNOSZ package
+
+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", parameters = data.frame(omega = omega.PrTr, Z = Z), T = T, P = P, 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", parameters = data.frame(omega = omega.PrTr, Z = Z), T = T, P = P, contrib = "s")$aq
+ return(V_s[[1]][, 1]/convert(omega.PrTr, "cm3bar"))
+}
+
+EOSvar <- function(var, T, P, ...) {
+ # Get the variables of a term in a regression equation
+ # T (K), P (bar)
+ Theta <- 228 # K
+ Psi <- 2600 # bar
+ out <- switch(EXPR = var,
+ "(Intercept)" = rep(1, length(T)),
+ "T" = T,
+ "P" = P,
+ "TTheta" = T - Theta, # T-Theta
+ "invTTheta" = (T - Theta)^-1, # 1/(T-Theta)
+ "TTheta2" = (T - Theta)^2, # (T-Theta)^2
+ "invTTheta2" = (T - Theta)^-2, # 1/(T-Theta)^2
+ "invPPsi" = (P + Psi)^-1, # 1/(P+Psi)
+ "invPPsiTTheta" = (P + Psi)^-1 * (T - Theta)^-1, # 1/[(P+Psi)(T-Theta)]
+ "TXBorn" = T*water("XBorn", T = T, P = P)[, 1],
+ "drho.dT" = -water("rho", T = T, P = P)[, 1]*water("E", T = T, P = P)[, 1],
+ "V.kT" = water("V", T = T, P = P)[, 1]*water("kT", T = T, P = P)[, 1],
+ # Fallback: get a variable that is a property of water, or
+ # is any other function by name (possibly a user-defined function)
+ ( if(var %in% water.SUPCRT92()) water(var, T, P)[, 1]
+ else if(exists(var)) {
+ if(is.function(get(var))) {
+ if(all(c("T", "P") %in% names(formals(get(var))))) get(var)(T = T, P = P, ...)
+ else stop(paste("the arguments of ", var, "() do not contain T and P", sep = ""))
+ }
+ else stop(paste("an object named", var, "is not a function"))
+ }
+ else stop(paste("can't find a variable named", var))
+ )
+ )
+ # 20151126 Apply the negative sign in the HKF EOS for V to the variable
+ # (not to omega as previously assumed)
+ if(var == "QBorn") out <- -out
+ return(out)
+}
+
+EOSlab <- function(var, coeff = "") {
+ # Make pretty labels for the variables
+ lab <- switch(EXPR = var,
+ # These are regression variables listed in EOSregress.Rd
+ "(Intercept)" = substitute(YYY*" ", list(YYY = coeff)),
+ "T" = substitute(YYY%*%italic(T), list(YYY = coeff)),
+ "P" = substitute(YYY%*%italic(P), list(YYY = coeff)),
+ "TTheta" = substitute(YYY%*%(italic(T)-Theta), list(YYY = coeff)),
+ "invTTheta" = substitute(YYY/(italic(T)-Theta), list(YYY = coeff)),
+ "TTheta2" = substitute(YYY%*%(italic(T)-Theta)^2, list(YYY = coeff)),
+ "invTTheta2" = substitute(YYY/(italic(T)-Theta)^2, list(YYY = coeff)),
+ "invPPsi" = substitute(YYY/(italic(P)+Psi),list(YYY = coeff)),
+ "invPPsiTTheta" = substitute(YYY/((italic(P)+Psi)(italic(T)-Theta)), list(YYY = coeff)),
+ "TXBorn" = substitute(YYY%*%italic(TX), list(YYY = coeff)),
+ "drho.dT" = substitute(YYY%*%(d~rho/dT), list(YYY = coeff)),
+ "V.kT" = substitute(YYY%*%V~kappa[italic(T)], list(YYY = coeff)),
+ # These are non-single-letter properties of water as listed in water.Rd
+ "kT" = substitute(YYY%*%kappa[italic(T)], list(YYY = coeff)),
+ "alpha" = substitute(YYY%*%alpha, list(YYY = coeff)),
+ "beta" = substitute(YYY%*%beta, list(YYY = coeff)),
+ "epsilon" = substitute(YYY%*%epsilon, list(YYY = coeff)),
+ "rho" = substitute(YYY%*%rho, list(YYY = coeff)),
+ "NBorn" = substitute(YYY%*%italic(N), list(YYY = coeff)),
+ "QBorn" = substitute(YYY%*%italic(Q), list(YYY = coeff)),
+ "XBorn" = substitute(YYY%*%italic(X), list(YYY = coeff)),
+ "YBorn" = substitute(YYY%*%italic(Y), list(YYY = coeff)),
+ "ZBorn" = substitute(YYY%*%italic(Z), list(YYY = coeff)),
+ (
+ # If var is a function, does have an attribute named "label"?
+ if(exists(var)) {
+ if(is.function(get(var))) {
+ if(!is.null(attr(get(var), "label"))) {
+ return(substitute(YYY*XXX, list(YYY = coeff, XXX = attr(get(var), "label"))))
+ # Fallback, use the name of the variable
+ # (e.g. for a property of water such as A, G, S, U, H, or name of a user-defined function)
+ } else substitute(YYY%*%italic(XXX), list(YYY = coeff, XXX = var))
+ } else substitute(YYY%*%italic(XXX), list(YYY = coeff, XXX = var))
+ } else substitute(YYY%*%italic(XXX), list(YYY = coeff, XXX = var))
+ )
+ )
+ return(lab)
+}
+
+EOSregress <- function(exptdata, var = "", T.max = 9999, ...) {
+ # Regress exptdata using terms listed in fun
+ # Which values to use
+ iT <- which(exptdata$T <= T.max)
+ exptdata <- exptdata[iT, ]
+ # Temperature and pressure
+ T <- exptdata$T
+ P <- exptdata$P
+ # The third column is the property of interest: Cp or V
+ X <- exptdata[, 3]
+ # Now build a regression formula
+ if(length(var) == 0) stop("var is missing")
+ fmla <- as.formula(paste("X ~ ", paste(var, collapse = "+")))
+ # Retrieve the values of the variables
+ for(i in seq_along(var)) assign(var[i], EOSvar(var[i], T = T, P = P, ...))
+ # Now regress away!
+ EOSlm <- lm(fmla)
+ return(EOSlm)
+}
+
+EOScalc <- function(coefficients, T, P, ...) {
+ # Calculate values of volume or heat capacity from regression fit
+ X <- 0
+ for(i in 1:length(coefficients)) {
+ coeff.i <- coefficients[[i]]
+ fun.i <- EOSvar(names(coefficients)[i], T, P, ...)
+ X <- X + coeff.i * fun.i
+ }
+ return(X)
+}
+
+EOSplot <- function(exptdata, var = NULL, T.max = 9999, T.plot = NULL,
+ fun.legend = "topleft", coefficients = NULL, add = FALSE,
+ lty = par("lty"), col = par("col"), ...) {
+ # Plot experimental and modelled volumes and heat capacities
+ # First figure out the property (Cp or V) from the exptdata
+ prop <- colnames(exptdata)[3]
+ # If var is NULL use HKF equations
+ if(is.null(var)) {
+ if(prop == "Cp") var <- c("invTTheta2","TXBorn")
+ if(prop == "V") var <- c("invTTheta","QBorn")
+ }
+ # Perform the regression, only using temperatures up to T.max
+ if(is.null(coefficients)) {
+ EOSlm <- EOSregress(exptdata, var, T.max, ...)
+ coefficients <- EOSlm$coefficients
+ }
+ # Only plot points below a certain temperature
+ iexpt <- 1:nrow(exptdata)
+ if(!is.null(T.plot)) iexpt <- which(exptdata$T < T.plot)
+ # For a nicer plot, extend the ranges, but don't go below -20 degrees C
+ ylim <- extendrange(exptdata[iexpt, prop], f = 0.1)
+ xlim <- extendrange(exptdata$T[iexpt], f = 0.1)
+ xlim[xlim < 253.15] <- 253.15
+ # Start plot
+ if(!add) {
+ thermo.plot.new(xlim = xlim, ylim = ylim, xlab = axis.label("T", units = "K"),
+ ylab = axis.label(paste(prop, "0", sep = "")), yline = 2, mar = NULL)
+ # Different plot symbols to represent size of residuals
+ pch.open <- 1
+ pch.filled <- 16
+ # Find the calculated values at these conditions
+ calc.X <- EOScalc(coefficients, exptdata$T, exptdata$P, ...)
+ expt.X <- exptdata[, prop]
+ # Are we within 10% of the values?
+ in10 <- which(abs((calc.X-expt.X)/expt.X) < 0.1)
+ pch <- rep(pch.open, length(exptdata$T))
+ pch[in10] <- pch.filled
+ points(exptdata$T, exptdata[, prop], pch = pch)
+ }
+ # Plot regression line at a single P
+ P <- mean(exptdata$P)
+ message("EOSplot: plotting line for P = ", P, " bar")
+ xs <- seq(xlim[1], xlim[2], length.out = 200)
+ calc.X <- EOScalc(coefficients, xs, P, ...)
+ lines(xs, calc.X, lty = lty, col = col)
+ # Make legend
+ if(!is.null(fun.legend) & !add) {
+ # 20161101: negate QBorn and V_s_var
+ iQ <- names(coefficients) %in% c("QBorn", "V_s_var")
+ coefficients[iQ] <- -coefficients[iQ]
+ coeffs <- as.character(round(as.numeric(coefficients), 4))
+ # So that positive ones appear with a plus sign
+ ipos <- which(coeffs >= 0)
+ coeffs[ipos] <- paste("+", coeffs[ipos], sep = "")
+ # Make labels for the functions
+ fun.lab <- as.expression(lapply(1:length(coeffs),
+ function(x) {EOSlab(names(coefficients)[x],coeffs[x])} ))
+ #fun.lab <- paste(names(coeffs),round(as.numeric(coeffs),4))
+ legend(fun.legend, legend = fun.lab, pt.cex = 0.1)
+ }
+ return(invisible(list(xrange = range(exptdata$T[iexpt]), coefficients = coefficients)))
+}
+
+EOScoeffs <- function(species, property, P = 1) {
+ # Get the HKF coefficients for species in the database
+ iis <- info(info(species, "aq"))
+ if(property == "Cp") {
+ out <- as.numeric(iis[,c("c1", "c2", "omega")])
+ names(out) <- c("(Intercept)", "invTTheta2", "TXBorn")
+ } else if(property == "V") {
+ iis <- iis[,c("a1", "a2", "a3", "a4", "omega")]
+ # Calculate sigma and xi and convert to volumetric units: 1 J = 10 cm^3 bar
+ sigma <- convert( iis$a1 + iis$a2 / (2600 + P), "cm3bar" )
+ xi <- convert( iis$a3 + iis$a4 / (2600 + P), "cm3bar" )
+ omega <- convert( iis$omega, "cm3bar" )
+ # 20151126: We _don't_ put a negative sign on omega here;
+ # now, the negative sign in the HKF EOS is with the variable (QBorn or V_s_var)
+ out <- c(sigma, xi, omega)
+ names(out) <- c("(Intercept)", "invTTheta", "QBorn")
+ }
+ return(out)
+}
+
Modified: pkg/CHNOSZ/R/examples.R
===================================================================
--- pkg/CHNOSZ/R/examples.R 2023-05-15 09:04:46 UTC (rev 783)
+++ pkg/CHNOSZ/R/examples.R 2023-05-15 10:21:42 UTC (rev 784)
@@ -13,7 +13,7 @@
"solubility", "equilibrate",
"diagram", "mosaic", "mix",
"buffer", "nonideal", "NaCl",
- "add.protein", "ionize.aa", "rank.affinity",
+ "add.protein", "ionize.aa", "EOSregress", "rank.affinity",
"DEW", "logB.to.OBIGT", "stack_mosaic"
)
plot.it <- FALSE
@@ -33,7 +33,7 @@
demos <- function(which = c("sources", "protein.equil", "affinity", "NaCl", "density",
"ORP", "ionize", "buffer", "protbuff", "glycinate",
"mosaic", "copper", "arsenic", "solubility", "gold", "contour", "sphalerite", "minsol",
- "Shh", "saturation", "DEW", "lambda", "potassium", "TCA", "aluminum",
+ "Shh", "saturation", "adenine", "DEW", "lambda", "potassium", "TCA", "aluminum",
"AD", "comproportionation", "Pourbaix", "E_coli", "yttrium", "rank.affinity"), save.png = FALSE) {
# Run one or more demos from CHNOSZ with ask = FALSE, and return the value of the last one
out <- NULL
Modified: pkg/CHNOSZ/demo/00Index
===================================================================
--- pkg/CHNOSZ/demo/00Index 2023-05-15 09:04:46 UTC (rev 783)
+++ pkg/CHNOSZ/demo/00Index 2023-05-15 10:21:42 UTC (rev 784)
@@ -19,6 +19,7 @@
dehydration log K of dehydration reactions; SVG file contains tooltips and links
Shh Affinities of transcription factors relative to Sonic hedgehog
saturation Equilibrium activity diagram showing activity ratios and mineral saturation limits
+adenine HKF parameters regressed from heat capacity and volume of aqueous adenine
DEW Deep Earth Water (DEW) model for high pressures
lambda Thermodynamic properties of lambda transition in quartz
potassium Comparison of thermodynamic datasets for predicting mineral stabilities
Copied: pkg/CHNOSZ/demo/adenine.R (from rev 715, pkg/CHNOSZ/demo/adenine.R)
===================================================================
--- pkg/CHNOSZ/demo/adenine.R (rev 0)
+++ pkg/CHNOSZ/demo/adenine.R 2023-05-15 10:21:42 UTC (rev 784)
@@ -0,0 +1,122 @@
+# CHNOSZ/demos/adenine.R
+# Plot thermodynamic data and model fits for aqueous adenine 20170725
+library(CHNOSZ)
+
+# Notable functions in this demo:
+# EOSregress() - to regress HKF coefficients from Cp data
+# mod.OBIGT() - to modify the thermodynamic database for comparisons with an older set of parameters for adenine
+
+# LCT17 = Lowe et al., 2017 (J. Chem. Thermodyn., doi:10.1016/j.jct.2017.04.005)
+# LH06 = LaRowe and Helgeson, 2006 (Geochim. Cosmochim. Acta, doi:10.1016/j.gca.2006.04.010)
+# HKF = Helgeson-Kirkham-Flowers equations (e.g. Am. J. Sci., doi:10.2475/ajs.281.10.1249)
+
+# Cp0 and V0 of adenine from LCT17 Table 4
+AdH <- data.frame(
+ T = seq(288.15, 363.15, 5),
+ V = c(89.1, 89.9, 90.6, 91.3, 92, 92.7, 93.1, 93.6, 94.1, 94.9, 95.4, 95.9, 96.3, 96.9, 97.1, 97.8),
+ V_SD = c(1.1, 1.3, 1.1, 1, 1.1, 1, 0.9, 0.9, 0.8, 0.6, 0.7, 0.7, 0.7, 0.4, 1.1, 0.7),
+ Cp = c(207, 212, 216, 220, 224, 227, 230, 234, 237, 241, 243, 245, 248, 250, 252, 255),
+ Cp_SD = c(5, 7, 8, 7, 8, 7, 6, 7, 6, 5, 6, 6, 5, 4, 7, 5)
+)
+# Functions to calculate V and Cp using density model (LCT17 Eq. 28)
+Vfun <- function(v1, v2, k, T) {
+ # Gas constant (cm3 bar K-1 mol-1)
+ R <- 83.144598
+ # Isothermal compressibility (bar-1)
+ beta <- water("beta", TK)$beta
+ v1 + v2 / (T - 228) - k * R * beta
+}
+Cpfun <- function(c1, c2, k, T) {
+ # Gas constant (J K-1 mol-1)
+ R <- 8.3144598
+ # Isobaric temperature derivative of expansivity (K-2)
+ daldT <- water("daldT", TK)$daldT
+ c1 + c2 / (T - 228) ^ 2 - k * R * T * daldT
+}
+# Set up units (used for axis labels and HKF calculations)
+T.units("K")
+# Temperature and pressure points for calculations
+TK <- seq(275, 425)
+P <- water("Psat", TK)$Psat
+# Set up plots
+opar <- par(no.readonly = TRUE)
+layout(matrix(1:3), heights=c(1, 8, 8))
+# Title at top
+par(mar=c(0, 0, 0, 0), cex=1)
+plot.new()
+text(0.5, 0.5, "Heat capacity and volume of aqueous adenine\n(Lowe et al., 2017)", font=2)
+# Settings for plot areas
+par(mar = c(4, 4, 0.5, 1), mgp = c(2.4, 0.5, 0))
+# Location of x-axis tick marks
+xaxp <- c(275, 425, 6)
+
+### Cp plot (LCT17 Figures 4 and 12) ###
+plot(AdH$T, AdH$Cp, type = "p", xlim = range(TK), ylim = c(150, 350),
+ xlab = axis.label("T"), ylab = axis.label("Cp0"),
+ pch = 5, tcl = 0.3, xaxs = "i", yaxs = "i", las = 1, xaxp = xaxp
+)
+# Error bars (arrows trick from https://stackoverflow.com/questions/13032777/scatter-plot-with-error-bars)
+arrows(AdH$T, AdH$Cp - AdH$Cp_SD, AdH$T, AdH$Cp + AdH$Cp_SD, length = 0.05, angle = 90, code = 3)
+# Get LH06 predictions using HKF model;
+# this version of adenine parameters has been superseded by LCT17,
+# so we add them by hand
+mod.OBIGT("adenine-old", formula="C5H5N5", a1=21.5046, a2=8.501, a3=-2.6632, a4=-5.3561, c1=87.88, c2=-15.87, omega=0.065)
+LH06 <- subcrt("adenine-old", T = TK)$out$adenine
+lines(TK, LH06$Cp, lty = 3)
+# Density model (parameters from LCT17 Table 11)
+lines(TK, Cpfun(160.4, -653239, -7930.3, TK), lty = 2)
+
+# Regress HKF parameters
+# Specify the terms in the HKF equations
+var <- c("invTTheta2", "TXBorn")
+# Build a data frame with T, P, and Cp columns
+Cpdat <- data.frame(T = AdH[, "T"], P = 1, Cp = AdH[, "Cp"])
+# Convert Cp data from J to cal
+Cpdat$Cp <- convert(Cpdat$Cp, "cal")
+# Regress HKF parameters from Cp data
+HKFcoeffs <- EOSregress(Cpdat, var)$coefficients
+# Get predictions from the fitted model
+Cpfit <- EOScalc(HKFcoeffs, TK, P)
+# Plot the fitted model
+lines(TK, convert(Cpfit, "J"), lwd = 2, col = "green3")
+# Format coefficients for legend; use scientific notation for c2 and omega
+coeffs <- format(signif(HKFcoeffs, 4), scientific = TRUE)
+# Keep c1 out of scientific notation
+coeffs[1] <- signif(HKFcoeffs[[1]], 4)
+ipos <- which(coeffs >= 0)
+coeffs[ipos] <- paste("+", coeffs[ipos], sep = "")
+fun.lab <- as.expression(lapply(1:length(coeffs), function(x) {
+ EOSlab(names(HKFcoeffs)[x], coeffs[x])
+}))
+# Add legend: regressed HKF coefficients
+legend("topleft", legend = fun.lab, pt.cex = 0.1, box.col = "green3")
+# Add legend: lines
+legend("bottomright", lty = c(3, 2, 1), lwd = c(1, 1, 2), col = c("black", "black", "green3"), bty = "n",
+ legend = c("HKF model (LaRowe and Helgeson, 2006)",
+ "density model (Lowe et al., 2017)", "HKF model (fit using CHNOSZ)")
+)
+
+### V plot (LCT17 Figures 3 and 11) ###
+plot(AdH$T, AdH$V, type = "p", xlim = range(TK), ylim = c(85, 105),
+ xlab = axis.label("T"), ylab = axis.label("V0"),
+ pch = 5, tcl = 0.3, xaxs = "i", yaxs = "i", las = 1, xaxp = xaxp
+)
+axis(3, labels = FALSE, tcl = 0.3, xaxp = xaxp)
+axis(4, labels = FALSE, tcl = 0.3)
+arrows(AdH$T, AdH$V - AdH$V_SD, AdH$T, AdH$V + AdH$V_SD, length = 0.05, angle = 90, code = 3)
+# HKF model with coefficients from LH06
+lines(TK, LH06$V, lty = 3)
+# Density model with coefficients from LCT17
+lines(TK, Vfun(73.9, -917.0, -7930.3, TK), lty = 2)
+# HKF heat capacity coefficients from LCT17
+LCT17 <- subcrt("adenine", T = TK)$out$adenine
+lines(TK, LCT17$V, lwd = 2, col = "royalblue")
+legend("bottomright", lty = c(3, 2, 1), lwd = c(1, 1, 2), col = c("black", "black", "royalblue"), bty = "n",
+ legend=c("HKF model (LaRowe and Helgeson, 2006)",
+ "density model (Lowe et al., 2017)", "HKF model (fit by Lowe et al., 2017 using CHNOSZ)")
+)
+# Reset database and computational settings
+reset()
+
+layout(matrix(1))
+par(opar)
Modified: pkg/CHNOSZ/inst/NEWS.Rd
===================================================================
--- pkg/CHNOSZ/inst/NEWS.Rd 2023-05-15 09:04:46 UTC (rev 783)
+++ pkg/CHNOSZ/inst/NEWS.Rd 2023-05-15 10:21:42 UTC (rev 784)
@@ -12,6 +12,17 @@
% 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 2.0.0-4 (2023-05-15)}{
+
+ \itemize{
+
+ \item Restore EOSregress.R and demo/adenine.R.
+
+ }
+
+}
+
+
\section{Changes in CHNOSZ version 2.0.0 (2023-03-13)}{
\subsection{MAJOR USER-VISIBLE CHANGES}{
Copied: pkg/CHNOSZ/inst/tinytest/test-EOSregress.R (from rev 715, pkg/CHNOSZ/inst/tinytest/test-EOSregress.R)
===================================================================
--- pkg/CHNOSZ/inst/tinytest/test-EOSregress.R (rev 0)
+++ pkg/CHNOSZ/inst/tinytest/test-EOSregress.R 2023-05-15 10:21:42 UTC (rev 784)
@@ -0,0 +1,45 @@
+# Load default settings for CHNOSZ
+reset()
+
+info <- "EOSvar stops with unknown variables"
+expect_error(EOSvar("TX", T = 25, P = 1), "can't find a variable named TX", info = info)
+# Why can't the test find these?
+#TX <- 2
+#expect_error(EOSvar("TX", T = 25, P = 1), "an object named TX is not a function")
+#TX <- function(T) 2
+#expect_error(EOSvar("TX", T = 25, P = 1), "the arguments of TX\\(\\) are not T, P")
+
+info <- "Regressions return known HKF parameters (neutral species)"
+# Regress computed values of heat capacity and volume of CH4(aq)
+# calculated from HKF parameters on a T-P grid
+T <- convert(seq(0, 350, 25), "K")
+P <- seq(200, 1000, 100)
+# We use calories to compare with OBIGT data for CH4(aq) 20220325
+T.units("K")
+E.units("cal")
+CH4.prop <- subcrt("CH4", T = T, P = P, grid = "T")$out[[1]]
+# Terms in the HKF equations for Cp
+Cp.var <- c("invTTheta2", "TXBorn")
+# Get coefficients in Cp regression
+Cp.lm <- EOSregress(CH4.prop[, c("T", "P", "Cp")], Cp.var)
+Cp.coeff <- Cp.lm$coefficients
+# Terms in the HKF equations for V
+V.var <- c("invPPsi", "invTTheta", "invPPsiTTheta", "QBorn")
+# Get coefficients in V regression
+V.lm <- EOSregress(CH4.prop[, c("T", "P", "V")], V.var)
+# Use same units as HKF: convert from cm3.bar to joules (divide by 10) then to calories
+V.coeff <- convert(convert(V.lm$coefficients, "joules"), "cal")
+## The tests: did we get the HKF parameters that are in the database?
+CH4.par <- info(info("CH4"))
+# c1 and c2
+expect_equivalent(Cp.coeff[1], CH4.par$c1, info = info)
+expect_equivalent(Cp.coeff[2], CH4.par$c2, info = info)
+# omega (from Cp)
+expect_equivalent(Cp.coeff[3], CH4.par$omega, info = info)
+# a1, a2, a3 and a4
+expect_equivalent(V.coeff[1], CH4.par$a1, info = info)
+expect_equivalent(V.coeff[2], CH4.par$a2, info = info)
+expect_equivalent(V.coeff[3], CH4.par$a3, info = info)
+expect_equivalent(V.coeff[4], CH4.par$a4, info = info)
+# omega (from V)
+expect_equivalent(V.coeff[5], CH4.par$omega, info = info)
Modified: pkg/CHNOSZ/man/CHNOSZ-package.Rd
===================================================================
--- pkg/CHNOSZ/man/CHNOSZ-package.Rd 2023-05-15 09:04:46 UTC (rev 783)
+++ pkg/CHNOSZ/man/CHNOSZ-package.Rd 2023-05-15 10:21:42 UTC (rev 784)
@@ -20,7 +20,7 @@
Each help page (other than this one) has been given one of the following \dQuote{concept index entries}:
\itemize{
\item Main workflow: \code{\link{info}}, \code{\link{subcrt}}, \code{\link{basis}}, \code{\link{species}}, \code{\link{affinity}}, \code{\link{equilibrate}}, \code{\link{diagram}}
- \item Extended workflow: \code{\link{swap.basis}}, \code{\link{buffer}}, \code{\link{mosaic}}
+ \item Extended workflow: \code{\link{swap.basis}}, \code{\link{buffer}}, \code{\link{mosaic}}, \code{\link{EOSregress}}
\item Thermodynamic data: \code{\link{data}}, \code{\link{extdata}}, \code{\link{add.OBIGT}}, \code{\link{util.data}}
\item Thermodynamic calculations: \code{\link{util.formula}}, \code{\link{makeup}}, \code{\link{util.units}}, \code{\link{Berman}}, \code{\link{nonideal}}, \code{\link{util.misc}}
\item Water properties: \code{\link{water}}, \code{\link{util.water}}, \code{\link{DEW}}, \code{\link{IAPWS95}}
Copied: pkg/CHNOSZ/man/EOSregress.Rd (from rev 715, pkg/CHNOSZ/man/EOSregress.Rd)
===================================================================
--- pkg/CHNOSZ/man/EOSregress.Rd (rev 0)
+++ pkg/CHNOSZ/man/EOSregress.Rd 2023-05-15 10:21:42 UTC (rev 784)
@@ -0,0 +1,202 @@
+\encoding{UTF-8}
+\name{EOSregress}
+\alias{EOSregress}
+\alias{EOSvar}
+\alias{EOScalc}
+\alias{EOSplot}
+\alias{EOSlab}
+\alias{EOScoeffs}
+\alias{Cp_s_var}
+\alias{V_s_var}
+\title{Regress Equations-of-State Parameters for Aqueous Species}
+\description{
+Fit experimental volumes and heat capacities using regression equations.
+Possible models include the Helgeson-Kirkham-Flowers (HKF) equations of state, or other equations defined using any combination of terms derived from the temperature, pressure and thermodynamic and electrostatic properties of water.
+}
+
+\usage{
+ EOSregress(exptdata, var = "", T.max = 9999, ...)
+ EOSvar(var, T, P, ...)
+ EOScalc(coefficients, T, P, ...)
+ EOSplot(exptdata, var = NULL, T.max = 9999, T.plot = NULL,
+ fun.legend = "topleft", coefficients = NULL, add = FALSE,
+ lty = par("lty"), col=par("col"), ...)
+ EOSlab(var, coeff = "")
+ EOScoeffs(species, property, P=1)
+ Cp_s_var(T = 298.15, P = 1, omega.PrTr = 0, Z = 0)
+ V_s_var(T = 298.15, P = 1, omega.PrTr = 0, Z = 0)
+}
+
+\arguments{
+ \item{exptdata}{dataframe, experimental data}
+ \item{var}{character, name(s) of variables in the regression equations}
+ \item{T.max}{numeric, maximum temperature for regression, in degrees Kelvin}
+ \item{T}{numeric, temperature in Kelvin}
+ \item{P}{numeric, pressure in bars}
+ \item{...}{arguments specifying additional dependencies of the regression variables}
+ \item{T.plot}{numeric, upper limit of temperature range to plot}
+ \item{fun.legend}{character, where to place legend on plot}
+ \item{coefficients}{dataframe, coefficients to use to make line on plot}
+ \item{add}{logical, add lines to an existing plot?}
+ \item{lty}{line style}
+ \item{col}{color of lines}
+ \item{coeff}{numeric, value of equation of state parameter for plot legend}
+ \item{species}{character, name of aqueous species}
+ \item{property}{character, \samp{Cp} or \samp{V}}
+ \item{omega.PrTr}{numeric, value of omega at reference T and P}
+ \item{Z}{numeric, charge}
+}
+
+\details{
+\code{EOSregress} uses a linear model (\code{\link{lm}}) to regress the experimental heat capacity or volume data in \code{exptdata}, which is a data frame with columns \samp{T} (temperature in degrees Kelvin), \samp{P} (pressure in bars), and \samp{Cp} or \samp{V} (heat capacity in cal/mol.K or volume in cm3/mol).
+The \samp{Cp} or \samp{V} data must be in the third column.
+Only data below the temperature of \code{T.max} are included in the regression.
+The regression formula is specified by a vector of names in \code{var}.
+The names of the variables can be any combination of the following (listed in the order of search): variables listed in the following table, any available property of \code{\link{water}} (e.g. \samp{V}, \samp{alpha}, \samp{QBorn}), or the name of a function that can be found using \code{\link{get}} in the default environment.
+Examples of the latter are \code{Cp_s_var}, \code{V_s_var}, or functions defined by the user in the global environment; the arguments of these functions must include, but are not limited to, \code{T} and \code{P}.
+
+ \tabular{ll}{
+ \code{T} \tab \eqn{T}{T} (temperature) \cr
+ \code{P} \tab \eqn{P}{P} (pressure) \cr
+ \code{TTheta} \tab \eqn{(T-\Theta)}{(T-Theta)} (\eqn{\Theta}{Theta} = 228 K) \cr
+ \code{invTTheta} \tab \eqn{1/(T-\Theta)}{1/(T-Theta)} \cr
+ \code{TTheta2} \tab \eqn{(T-\Theta)^2}{(T-Theta)^2} \cr
+ \code{invTTheta2} \tab \eqn{1/(T-\Theta)^2}{1/(T-Theta)^2} \cr
+ \code{invPPsi} \tab \eqn{1/(P+\Psi)}{1/(P+Psi)} (\eqn{\Psi}{Psi} = 2600 bar) \cr
+ \code{invPPsiTTheta} \tab \eqn{1/((P+\Psi)(T-\Theta))}{1/((P+Psi)(T-Theta))} \cr
+ \code{TXBorn} \tab \eqn{TX}{TX} (temperature times \eqn{X}{X} Born function) \cr
+ \code{drho.dT} \tab \eqn{d\rho/dT}{drho/dT} (temperature derivative of density of water) \cr
+ \code{V.kT} \tab \eqn{V\kappa_T}{V.kT} (volume times isothermal compressibility of water)
+ }
+
+
+\code{EOSvar} calculates the value of the variable named \code{var} (defined as described above) at the specified \code{T} (temperature in degrees Kelvin) and \code{P} (pressure in bars).
+This function is used by \code{EOSregress} to get the values of the variables used in the regression.
+
+\code{EOScalc} calculates the predicted heat capacities or volumes using coefficients provided by the result of \code{EOSregress}, at the temperatures and pressures specified by \code{T} and \code{P}.
+
+\code{EOSplot} takes a table of data in \code{exptdata}, runs \code{EOSregress} and \code{EOScalc} and plots the results.
+The experimental data are plotted as points, and the calculated values as a smooth line.
+The point symbols are filled circles where the calculated value is within 10\% of the experimental value; open circles otherwise.
+
+\code{EOSlab} produces labels for the variables listed above that can be used \code{\link{as.expression}}s in plots.
+The value of \code{coeff} is prefixed to the name of the variable (using \code{\link{substitute}}, with a multiplication symbol).
+For the properties listed in the table above, and selected properties listed in \code{\link{water}}, the label is formatted using \code{\link{plotmath}} expressions (e.g., with italicized symbols and Greek letters).
+If \code{var} is a user-defined function, the function can be given a \samp{label} attribute to provide \code{\link{plotmath}}-style formatting; in this case the appropriate multiplication or division symbol should be specified (see example below).
+
+\code{EOScoeffs} retrieves coefficients in the Helgeson-Kirkham-Flowers equations from the thermodynamic database (\code{\link{thermo}$OBIGT}) for the given aqueous \code{species}.
+If the \code{property} is \samp{Cp}, the resulting data frame has column names of \samp{(Intercept)}, \samp{invTTheta2} and \samp{TX}, respectively holding the coefficients \eqn{c_1}{c1}, \eqn{c_2}{c2} and \eqn{\omega}{omega} in the equation \eqn{Cp^\circ = c_1 + c_2/(T-\Theta)^2 + {\omega}TX}{Cp = c1 + c2/(T-Theta)^2 + omega*TX}.
+If the \code{property} is \samp{V}, the data frame has column names of \samp{(Intercept)}, \samp{invTTheta} and \samp{Q}, respectively holding the coefficients \eqn{\sigma}{sigma}, \eqn{\xi}{xi} and \eqn{\omega}{omega} in \eqn{V^\circ = \sigma + \xi/(T-\Theta) - {\omega}Q}{V = sigma + xi/(T-Theta) - omega*Q}.
+Here, \eqn{\sigma}{sigma} and \eqn{\xi}{xi} are calculated from a1, a2, a3 and a4 in \code{thermo()$OBIGT} at the pressure indicated by \code{P} (default 1 bar).
+
+The original motivation for writing these functions was to explore alternatives or possible modifications to the revised Helgeson-Kirkham-Flowers equations applied to aqueous nonelectrolytes.
+As pointed out by Schulte et al., 2001, the functional forms of the equations do not permit retrieving values of the solvation parameter (\eqn{\omega}{omega}) that closely represent the observed trends in both heat capacity and volume at high temperatures (above ca. 200 \degC).
+
+The examples below assume that the \eqn{\omega}{omega} parameter in the HKF functions is a constant (does not depend on T and P), as is appropriate for nonelectrolytes.
+For charged species, the variables \code{Cp_s_var} and \code{V_s_var} can be used in the regressions.
+They correspond to the solvation contribution to heat capacity or volume, respectively, in the HKF EOS, divided by the value of \eqn{\omega}{omega} at the reference temperature and pressure.
+Because these variables are themselves a function of \code{omega.PrTr}, an iterative procedure is needed to perform the regression.
+
+Note that variables \code{QBorn} and \code{V_s_var} are both negated, so that the value of \eqn{\omega}{omega} has its proper sign in the corresponding equations.
+}
+
+\value{
+For \code{EOSregress}, an object of class \dQuote{lm}.
+\code{EOSvar} and \code{EOScalc} both return numeric values.
+\code{EOScoeffs} returns a data frame.
+}
+
+\seealso{
+The vignette \viglink{eos-regress} has more references and examples, including an iterative method to retrieve \code{omega.PrTr}.
+}
+
+\references{
+ Hnědkovský, L. and Wood, R. H. (1997) Apparent molar heat capacities of aqueous solutions of \CH4, \CO2, \H2S, and \NH3 at temperatures from 304 K to 704 K at a pressure of 28 MPa. \emph{J. Chem. Thermodyn.} \bold{29}, 731--747. \doi{10.1006/jcht.1997.0192}
+
+ Schulte, M. D., Shock, E. L. and Wood, R. H. (1995) The temperature dependence of the standard-state thermodynamic properties of aqueous nonelectrolytes. \emph{Geochim. Cosmochim. Acta} \bold{65}, 3919--3930. \doi{10.1016/S0016-7037(01)00717-7}
+}
+
+
+\examples{
+\dontshow{reset()}
+## Fit experimental heat capacities of CH4
+## using revised Helgeson-Kirkham-Flowers equations
+# Read the data from Hnedkovsky and Wood, 1997
+f <- system.file("extdata/cpetc/HW97_Cp.csv", package = "CHNOSZ")
+d <- read.csv(f)
+# Use data for CH4
+d <- d[d$species == "CH4", ]
+d <- d[, -1]
+# Convert J to cal and MPa to bar
+d$Cp <- convert(d$Cp, "cal")
+d$P <- convert(d$P, "bar")
+# Specify the terms in the HKF equations
+var <- c("invTTheta2", "TXBorn")
+# Perform regression, with a temperature limit
+EOSlm <- EOSregress(d, var, T.max = 600)
+# Calculate the Cp at some temperature and pressure
+EOScalc(EOSlm$coefficients, 298.15, 1)
+# Get the database values of c1, c2 and omega for CH4(aq)
+CH4coeffs <- EOScoeffs("CH4", "Cp")
+## Make plots comparing the regressions
+## with the accepted EOS parameters of CH4
+opar <- par(mfrow = c(2,2))
+EOSplot(d, T.max = 600)
+title("Cp of CH4(aq), fit to 600 K")
+legend("bottomleft", pch = 1, legend = "Hnedkovsky and Wood, 1997")
+EOSplot(d, coefficients = CH4coeffs)
+title("Cp from EOS parameters in database")
+EOSplot(d, T.max = 600, T.plot = 600)
+title("Cp fit to 600 K, plot to 600 K")
+EOSplot(d, coefficients = CH4coeffs, T.plot = 600)
+title("Cp from EOS parameters in database")
+par(opar)
+
+# Continuing from above, with user-defined variables
+Theta <- 228 # K
+invTTTheta3 <- function(T, P) (2*T) / (T-T*Theta) ^ 3
+invTX <- function(T, P) 1 / T * water("XBorn", T = T, P = P)[,1]
+# Print the calculated values of invTTTheta3
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/chnosz -r 784
More information about the CHNOSZ-commits
mailing list