[CHNOSZ-commits] r716 - in pkg/CHNOSZ: . R demo inst inst/tinytest man vignettes
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sun Mar 27 05:25:03 CEST 2022
Author: jedick
Date: 2022-03-27 05:25:02 +0200 (Sun, 27 Mar 2022)
New Revision: 716
Removed:
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/R/util.data.R
pkg/CHNOSZ/demo/00Index
pkg/CHNOSZ/demo/copper.R
pkg/CHNOSZ/inst/NEWS.Rd
pkg/CHNOSZ/inst/tinytest/test-eos.R
pkg/CHNOSZ/man/CHNOSZ-package.Rd
pkg/CHNOSZ/man/examples.Rd
pkg/CHNOSZ/man/extdata.Rd
pkg/CHNOSZ/man/subcrt.Rd
pkg/CHNOSZ/man/util.formula.Rd
pkg/CHNOSZ/vignettes/anintro.Rmd
pkg/CHNOSZ/vignettes/customizing.Rmd
pkg/CHNOSZ/vignettes/mklinks.sh
pkg/CHNOSZ/vignettes/vig.bib
Log:
Remove EOSregress()
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2022-03-26 08:11:30 UTC (rev 715)
+++ pkg/CHNOSZ/DESCRIPTION 2022-03-27 03:25:02 UTC (rev 716)
@@ -1,6 +1,6 @@
-Date: 2022-03-26
+Date: 2022-03-27
Package: CHNOSZ
-Version: 1.4.3-8
+Version: 1.9.9-9
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 2022-03-26 08:11:30 UTC (rev 715)
+++ pkg/CHNOSZ/NAMESPACE 2022-03-27 03:25:02 UTC (rev 716)
@@ -28,7 +28,6 @@
"ionize.aa", "MP90.cp", "aasum",
"qqr", "RMSD", "CVRMSD", "spearman", "DGmix", "DDGmix", "DGtr",
"ratlab",
- "EOSregress", "EOScoeffs", "EOSplot", "EOSvar",
# demos
"protein.equil", "palply",
"label.plot",
@@ -41,9 +40,7 @@
# (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",
"DGinf", "SD", "pearson", "shannon", "CV", "logact",
- "EOSlab", "EOScalc",
"basis.elements", "element.mu", "ibasis",
"water.SUPCRT92", "plot_findit",
"nonideal", "uniprot.aa",
Deleted: pkg/CHNOSZ/R/EOSregress.R
===================================================================
--- pkg/CHNOSZ/R/EOSregress.R 2022-03-26 08:11:30 UTC (rev 715)
+++ pkg/CHNOSZ/R/EOSregress.R 2022-03-27 03:25:02 UTC (rev 716)
@@ -1,213 +0,0 @@
-# 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 2022-03-26 08:11:30 UTC (rev 715)
+++ pkg/CHNOSZ/R/examples.R 2022-03-27 03:25:02 UTC (rev 716)
@@ -15,7 +15,7 @@
"diagram", "mosaic", "mix",
"buffer", "nonideal", "NaCl",
"add.protein", "ionize.aa",
- "objective", "revisit", "EOSregress")
+ "objective", "revisit")
plot.it <- FALSE
if(is.character(save.png))
png(paste(save.png,"%d.png",sep=""),width=500,height=500,pointsize=12)
@@ -33,7 +33,7 @@
demos <- function(which=c("sources", "protein.equil", "affinity", "NaCl", "density",
"ORP", "findit", "ionize", "buffer", "protbuff", "glycinate",
"mosaic", "copper", "arsenic", "solubility", "gold", "contour", "sphalerite", "zinc",
- "Shh", "saturation", "adenine", "DEW", "lambda", "potassium", "TCA", "aluminum",
+ "Shh", "saturation", "DEW", "lambda", "potassium", "TCA", "aluminum",
"AD", "comproportionation", "Pourbaix", "E_coli"), save.png=FALSE) {
# run one or more demos from CHNOSZ with ask=FALSE, and return the value of the last one
for(i in 1:length(which)) {
Modified: pkg/CHNOSZ/R/util.data.R
===================================================================
--- pkg/CHNOSZ/R/util.data.R 2022-03-26 08:11:30 UTC (rev 715)
+++ pkg/CHNOSZ/R/util.data.R 2022-03-27 03:25:02 UTC (rev 716)
@@ -396,9 +396,9 @@
# Take a data frame in the format of thermo()$OBIGT of one or more rows,
# remove scaling factors from equations-of-state parameters,
# and apply new column names depending on the state.
-# And convert energy units from J to cal (used by subcrt()) 20190530
# If fixGHS is TRUE a missing one of G, H or S for any species is calculated
# from the other two and the chemical formula of the species.
+# If toJoules is TRUE, convert parameters to Joules 20220325
# This function is used by both info and subcrt when retrieving entries from the thermodynamic database.
OBIGT2eos <- function(OBIGT, state, fixGHS = FALSE, toJoules = FALSE) {
# remove scaling factors from EOS parameters
@@ -429,6 +429,8 @@
# We only convert column 20 for aqueous species (omega), not for cgl species (lambda) 20190903
if(identical(state, "aq")) OBIGT[ical, c(9:12, 14:20)] <- convert(OBIGT[ical, c(9:12, 14:20)], "J")
else OBIGT[ical, c(9:12, 14:19)] <- convert(OBIGT[ical, c(9:12, 14:19)], "J")
+ # Also update the E_units column 20220325
+ OBIGT$E_units[ical] <- "J"
}
}
if(fixGHS) {
Modified: pkg/CHNOSZ/demo/00Index
===================================================================
--- pkg/CHNOSZ/demo/00Index 2022-03-26 08:11:30 UTC (rev 715)
+++ pkg/CHNOSZ/demo/00Index 2022-03-27 03:25:02 UTC (rev 716)
@@ -20,7 +20,6 @@
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
Deleted: pkg/CHNOSZ/demo/adenine.R
===================================================================
--- pkg/CHNOSZ/demo/adenine.R 2022-03-26 08:11:30 UTC (rev 715)
+++ pkg/CHNOSZ/demo/adenine.R 2022-03-27 03:25:02 UTC (rev 716)
@@ -1,122 +0,0 @@
-# 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
-}
-# Aet 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/demo/copper.R
===================================================================
--- pkg/CHNOSZ/demo/copper.R 2022-03-26 08:11:30 UTC (rev 715)
+++ pkg/CHNOSZ/demo/copper.R 2022-03-27 03:25:02 UTC (rev 716)
@@ -33,7 +33,7 @@
if(i == 12) myG <- myG + getG("Cu+2") + 2*getG("glycinate")
if(i == 13) myG <- myG + getG("Cu+") + 2*getG("glycinate")
if(i == 14) myG <- myG + getG("Cu(Gly)+")
- # Energies are in Joules, so we have to change units of species in default OBIGT 20220326
+ # Energies are in Joules, so we have to change units of species in default OBIGT 20220325
mod.OBIGT(names[i], G = myG, E_units = "J")
}
Modified: pkg/CHNOSZ/inst/NEWS.Rd
===================================================================
--- pkg/CHNOSZ/inst/NEWS.Rd 2022-03-26 08:11:30 UTC (rev 715)
+++ pkg/CHNOSZ/inst/NEWS.Rd 2022-03-27 03:25:02 UTC (rev 716)
@@ -10,7 +10,7 @@
\newcommand{\s}{\ifelse{latex}{\eqn{_{#1}}}{\ifelse{html}{\out{<sub>#1</sub>}}{#1}}}
\newcommand{\S}{\ifelse{latex}{\eqn{^{#1}}}{\ifelse{html}{\out{<sup>#1</sup>}}{^#1}}}
-\section{Changes in CHNOSZ version 1.4.3-7 (2022-03-25)}{
+\section{Changes in CHNOSZ version 1.9.9-9 (2022-03-27)}{
\subsection{MAJOR CHANGE}{
\itemize{
@@ -22,11 +22,9 @@
to be modified to produce expected output.
\item As much as possible, functions in CHNOSZ have been modified to use
- Joules in internal variables, and the gas constant has been changed to
- Joules in \code{convert()} and other functions.
+ Joules in internal variables. In particular, the value of the gas
+ constant is based on Joules in \code{convert()} and other functions.
- \item The documentation is still catching up ...
-
}
}
@@ -33,15 +31,17 @@
\subsection{NEW FEATURES}{
\itemize{
- \item Add logB_to_OBIGT() to fit three thermodynamic parameters (G, S,
- and Cp) to formation constants of aqueous species as a function of
- temperature.
+ \item Add \code{logB_to_OBIGT()} to fit three thermodynamic parameters
+ (\samp{G}, \samp{S}, and \samp{Cp}) to formation constants of aqueous
+ species as a function of temperature.
- \item Add \code{zap} argument to mod.OBIGT() to clear parameters of
- preexisting species.
+ \item Add \code{zap} argument to \code{mod.OBIGT()} to clear parameters
+ of preexisting species.
- \item Add vignette customizing.Rmd about customizing the thermodynamic
- database.
+ \item Add vignette \strong{customizing.Rmd} with description of database
+ format, data-entry conventions, and examples of customizing the
+ thermodynamic database using \code{add.OBIGT()}, \code{mod.OBIGT()}, and
+ \code{logB_to_OBIGT()}.
}
}
@@ -49,13 +49,13 @@
\subsection{REMOVED FEATURES}{
\itemize{
- \item The eos-regress.Rmd vignette has been removed to debug problems
- associated with the switch to Joules.
-
- \item Remove exports of \code{cgl}, \code{hkf}, and \code{AD}.
- \code{subcrt} should be used for all calculations of thermodynamic
+ \item Remove exports of \code{cgl()}, \code{hkf()}, and \code{AD()}.
+ \code{subcrt()} should be used for all calculations of thermodynamic
properties.
+ \item \code{EOSregress()} and the associated demo and vignette have been
+ removed.
+
}
}
Deleted: pkg/CHNOSZ/inst/tinytest/test-EOSregress.R
===================================================================
--- pkg/CHNOSZ/inst/tinytest/test-EOSregress.R 2022-03-26 08:11:30 UTC (rev 715)
+++ pkg/CHNOSZ/inst/tinytest/test-EOSregress.R 2022-03-27 03:25:02 UTC (rev 716)
@@ -1,45 +0,0 @@
-# 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/inst/tinytest/test-eos.R
===================================================================
--- pkg/CHNOSZ/inst/tinytest/test-eos.R 2022-03-26 08:11:30 UTC (rev 715)
+++ pkg/CHNOSZ/inst/tinytest/test-eos.R 2022-03-27 03:25:02 UTC (rev 716)
@@ -76,7 +76,10 @@
expect_equal(gfun.1000$dgdP * 1e6, dgdP.1000.ref, tolerance = 1e-1, info = info)
expect_equal(gfun.4000$dgdP * 1e6, dgdP.4000.ref, tolerance = 1e-3, info = info)
-# load SiO2 and Si2O4 data taken from DEW spreadsheet
+info <- "hkf() and subcrt() give consistent values for non-solvation volume"
+# Added on 20220326 to check usage of subcrt() with omega = 0
+# 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"))
Vn1 <- Vn2 <- numeric()
species <- c("CO3-2", "BO2-", "MgCl+", "SiO2", "HCO3-", "Si2O4")
@@ -87,7 +90,7 @@
# Get the nonsolvation volume from hkf()
Vn1 <- c(Vn1, CHNOSZ:::hkf("V", par, contrib="n")$aq[[1]]$V)
}
-# In version 2.0.0, hkf() assumes the parameters are Joules,
+# In CHNOSZ 2.0.0, hkf() assumes the parameters are Joules,
# but we gave it calorie-based parameters, so we need to convert to Joules
Vn1 <- convert(Vn1, "J")
# Second method: subcrt() with omega = 0
@@ -97,8 +100,6 @@
}
expect_equal(Vn1, Vn2)
-
-
# Reference
# Shock, E. L., Oelkers, E. H., Johnson, J. W., Sverjensky, D. A. and Helgeson, H. C. (1992)
Modified: pkg/CHNOSZ/man/CHNOSZ-package.Rd
===================================================================
--- pkg/CHNOSZ/man/CHNOSZ-package.Rd 2022-03-26 08:11:30 UTC (rev 715)
+++ pkg/CHNOSZ/man/CHNOSZ-package.Rd 2022-03-27 03:25:02 UTC (rev 716)
@@ -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}}, \code{\link{objective}}, \code{\link{revisit}}, \code{\link{findit}}, \code{\link{EOSregress}}
+ \item Extended workflow: \code{\link{swap.basis}}, \code{\link{buffer}}, \code{\link{mosaic}}, \code{\link{objective}}, \code{\link{revisit}}, \code{\link{findit}}
\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}}
Deleted: pkg/CHNOSZ/man/EOSregress.Rd
===================================================================
--- pkg/CHNOSZ/man/EOSregress.Rd 2022-03-26 08:11:30 UTC (rev 715)
+++ pkg/CHNOSZ/man/EOSregress.Rd 2022-03-27 03:25:02 UTC (rev 716)
@@ -1,202 +0,0 @@
-\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.
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/chnosz -r 716
More information about the CHNOSZ-commits
mailing list