[CHNOSZ-commits] r114 - in pkg/CHNOSZ: . R inst man vignettes
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Nov 1 06:15:31 CET 2016
Author: jedick
Date: 2016-11-01 06:15:29 +0100 (Tue, 01 Nov 2016)
New Revision: 114
Added:
pkg/CHNOSZ/vignettes/EOSregress.Rmd
pkg/CHNOSZ/vignettes/EOSregress.bib
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/R/EOSregress.R
pkg/CHNOSZ/inst/NEWS
pkg/CHNOSZ/man/EOSregress.Rd
Log:
add EOSregress.Rmd
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2016-08-26 12:52:36 UTC (rev 113)
+++ pkg/CHNOSZ/DESCRIPTION 2016-11-01 05:15:29 UTC (rev 114)
@@ -1,11 +1,11 @@
-Date: 2016-08-26
+Date: 2016-11-01
Package: CHNOSZ
-Version: 1.0.8-1
+Version: 1.0.8-2
Title: Chemical Thermodynamics and Activity Diagrams
Author: Jeffrey Dick
Maintainer: Jeffrey Dick <j3ffdick at gmail.com>
Depends: R (>= 3.1.0)
-Suggests: limSolve, testthat, knitr
+Suggests: limSolve, testthat, knitr, rmarkdown, tufte
Imports: grDevices, graphics, stats, utils
Description: Functions and data sets to support chemical thermodynamic modeling in biochemistry
and low-temperature geochemistry. The features include calculation of the standard molal
Modified: pkg/CHNOSZ/R/EOSregress.R
===================================================================
--- pkg/CHNOSZ/R/EOSregress.R 2016-08-26 12:52:36 UTC (rev 113)
+++ pkg/CHNOSZ/R/EOSregress.R 2016-11-01 05:15:29 UTC (rev 114)
@@ -3,6 +3,19 @@
# 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", T=T, P=P, eos=data.frame(omega=omega.PrTr, Z=Z), contrib="s")
+ return(Cp_s[[1]][, 1]/omega.PrTr)
+}
+
+V_s_var <- function(T=298.15, P=1, omega.PrTr=0, Z=0) {
+ # solvation contribution to volume in the HKF EOS, divided by omega(Pr,Tr) (cm3.bar)
+ # [the negative sign on this term as written in the HKF EOS is accounted for by hkf()]
+ V_s <- hkf("V", T=T, P=P, eos=data.frame(omega=omega.PrTr, Z=Z), contrib="s")
+ 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)
@@ -102,20 +115,20 @@
return(EOSlm)
}
-EOScalc <- function(coefficients,T,P) {
+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)
+ 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) {
+EOSplot <- function(exptdata, var=NULL, T.max=9999, T.plot=NULL,
+ fun.legend="topleft", coefficients=NULL, add=FALSE, lty=1, ...) {
# plot experimental and modelled volumes and heat capacities
# first figure out the property (Cp or V) from the exptdata
prop <- colnames(exptdata)[3]
@@ -126,39 +139,43 @@
}
# perform the regression, only using temperatures up to T.max
if(is.null(coefficients)) {
- EOSlm <- EOSregress(exptdata, var, T.max)
+ 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)
- iX <- match(prop, colnames(exptdata))
- ylim <- extendrange(exptdata[iexpt, iX], f=0.1)
+ # 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
- 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[, iX]
- # 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[, iX], pch=pch)
- # take out NAs and Infinite values
- iNA <- is.na(calc.X) | is.infinite(calc.X)
+ 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)
msgout("EOSplot: plotting line for P=", P, " bar\n")
xs <- seq(xlim[1], xlim[2], length.out=200)
- calc.X <- EOScalc(coefficients, xs, P)
- lines(xs, calc.X)
+ calc.X <- EOScalc(coefficients, xs, P, ...)
+ lines(xs, calc.X, lty=lty)
# make legend
- if(!is.null(fun.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)
@@ -169,7 +186,7 @@
#fun.lab <- paste(names(coeffs),round(as.numeric(coeffs),4))
legend(fun.legend, legend=fun.lab, pt.cex=0.1)
}
- return(invisible(list(xlim=range(exptdata$T[iexpt]))))
+ return(invisible(list(xrange=range(exptdata$T[iexpt]), coefficients=coefficients)))
}
EOScoeffs <- function(species, property, P=1) {
Modified: pkg/CHNOSZ/inst/NEWS
===================================================================
--- pkg/CHNOSZ/inst/NEWS 2016-08-26 12:52:36 UTC (rev 113)
+++ pkg/CHNOSZ/inst/NEWS 2016-11-01 05:15:29 UTC (rev 114)
@@ -1,9 +1,10 @@
-CHANGES IN CHNOSZ 1.0.8-1 (2016-08-26)
+CHANGES IN CHNOSZ 1.0.8-2 (2016-11-01)
--------------------------------------
- Add "AA" as a keyword for preset species in basis() (cysteine,
glutamic acid, glutamine, H2O, oxygen).
+- Add EOSregress.Rmd vignette; update related functions.
CHANGES IN CHNOSZ 1.0.8 (2016-05-28)
------------------------------------
Modified: pkg/CHNOSZ/man/EOSregress.Rd
===================================================================
--- pkg/CHNOSZ/man/EOSregress.Rd 2016-08-26 12:52:36 UTC (rev 113)
+++ pkg/CHNOSZ/man/EOSregress.Rd 2016-11-01 05:15:29 UTC (rev 114)
@@ -6,42 +6,52 @@
\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 and/or user-defined functions of temperature and pressure.
+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)
+ EOScalc(coefficients, T, P, ...)
EOSplot(exptdata, var = NULL, T.max = 9999, T.plot = NULL,
- fun.legend = "topleft", coefficients = NULL)
+ fun.legend = "topleft", coefficients = NULL, add = FALSE,
+ lty = 1, ...)
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 degrees Kelvin}
+ \item{T}{numeric, temperature in Kelvin}
\item{P}{numeric, pressure in bars}
- \item{...}{arguments defining additional properties which variables are a function of}
+ \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}{numeric, line style}
\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 \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).
+\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).
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 (e.g. a function defined by the user in the global environment; the arguments of the function should be \code{T} and \code{P}; see example).
+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
@@ -74,12 +84,19 @@
\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}.
+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 motivation for writing these functions is to explore alternatives or possible modifications to the revised Helgeson-Kirkham-Flowers equations applied to aqueous nonelectrolytes.
+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 \eqn{^{\circ}}{°}C).
+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.
+See the \code{EOSregress} vignette for an example.
+
+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{
@@ -87,7 +104,7 @@
}
\seealso{
- See \code{\link{lm}} for the details of the regression calculations.
+The \code{EOSregress.Rmd} vignette has more references and examples.
}
\references{
@@ -115,11 +132,11 @@
dcoeffs <- EOSlm$coefficients - CH4coeffs
stopifnot(all(abs(dcoeffs/CH4coeffs) < 0.1))
## make plots comparing the regressions
-## here with the accepted EOS parameters of CH4
+## with the accepted EOS parameters of CH4
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")
+legend("bottomleft", pch=1, legend="Hnědkovský and Wood, 1997")
EOSplot(d, coefficients=CH4coeffs)
title("Cp from EOS parameters in database")
EOSplot(d, T.max=600, T.plot=600)
@@ -155,7 +172,7 @@
# for both HKF and the expansivity equation
# we'll fit up to a temperature limit
EOSplot(d, varHKF, T.max=663, T.plot=625)
-legend("bottomright", pch=1, legend="Hnedkovsky et al., 1996")
+legend("bottomright", pch=1, legend="Hnědkovský et al., 1996")
title("V of CH4(aq), HKF equation")
EOSplot(d, varal, T.max=663, T.plot=625)
title("V of CH4(aq), expansivity equation")
Added: pkg/CHNOSZ/vignettes/EOSregress.Rmd
===================================================================
--- pkg/CHNOSZ/vignettes/EOSregress.Rmd (rev 0)
+++ pkg/CHNOSZ/vignettes/EOSregress.Rmd 2016-11-01 05:15:29 UTC (rev 114)
@@ -0,0 +1,304 @@
+---
+title: "Regressing Thermodynamic Data"
+subtitle: "EOSregress in [CHNOSZ](http://www.chnosz.net)"
+author: "Jeffrey M. Dick"
+date: "`r Sys.Date()`"
+output:
+ tufte::tufte_html: default
+ tufte::tufte_handout:
+ citation_package: natbib
+ latex_engine: xelatex
+ tufte::tufte_book:
+ citation_package: natbib
+ latex_engine: xelatex
+vignette: >
+ %\VignetteIndexEntry{Regressing thermodynamic data}
+ %\VignetteEngine{knitr::rmarkdown}
+ \usepackage[utf8]{inputenc}
+bibliography: EOSregress.bib
+link-citations: yes
+---
+
+```{r setup, include=FALSE}
+# invalidate cache when the tufte version changes
+knitr::opts_chunk$set(tidy = FALSE, cache.extra = packageVersion('tufte'))
+options(htmltools.dir.version = FALSE)
+```
+
+
+This vignette demonstrates [`EOSregress`](http://www.chnosz.net/manual/EOSregress.html) and related functions for the regression of heat capacity and volumetric data to obtain "equations of state" coefficients.
+
+# A note on the equations
+
+The CHNOSZ thermodynamic database uses the revised Helgeson-Kirkham-Flowers equations of state (EOS) for aqueous species.
+Different terms in these equations give the "non-solvation" and "solvation" contributions to the standard molal heat capacity ($C_P^\circ$) and volume ($V^\circ$) as a function of temperature ($T$) and pressure ($P$).
+
+ at HKF81 published the original description of the equations; for $C_P^\circ$ the equation is
+$$C_P^\circ = c_1 + \frac{c_2}{(T-\Theta)^2} + {\omega}TX$$
+Here, $c_1$ and $c_2$ are the non-solvation parameters, and $\omega$ is the solvation parameter.
+$\Theta$ is equal to 228 K, and $X$ is one of the "Born functions" that relates the solvation process to the dielectric properties of water.
+For neutral species, all of the parameters, including the "effective" value of $\omega$, are constant.
+However, **for charged species, $\omega$ has derivatives, larger with higher temperature**, that depend on the charge and ionic radius.
+Revisions by @TH88 and @SOJSH92 refined the equations to improve their high-$T$ and $P$ behavior.
+
+# A note on the algorithms
+
+The regression functions are essentially a wrapper around the [`water`](http://www.chnosz.net/manual/water.html)-related functions in CHNOSZ and `lm` in R.
+The former provide values of the Born functions.
+Accordingly, numerical values for all of the *variables* in the equations (e.g. $\frac{1}{(T-\Theta)^2}$ and $TX$) can be calculated at each experimental $T$ and $P$.
+
+Applying a linear model (`lm`), the *coefficients* in the model can be obtained.
+In this case, they correspond to the HKF parameters, e.g. $c_1$ (the intercept), $c_2$ and $\omega$.
+The coefficients are constants, meaning that **the linear model is applicable only to neutral (uncharged) species**, for which the high-temperature derivatives of $\omega$ are 0.
+Further below, a procedure is described to iteratively solve the equations for charged species.
+
+# An example for neutral species
+
+```{r options, echo=FALSE}
+options(width = 90)
+```
+
+```{r chnosz, message=FALSE}
+library(CHNOSZ)
+data(thermo)
+```
+
+This is from the first example from `?EOSregress`.
+`r tufte::margin_note("The ? indicates a documentation topic in R. To view it, type ?EOSregress or help(EOSregress) at the command line.")`
+Here, we regress experimental heat capacities of aqueous methane ($\mathrm{CH_4}$) using the revised HKF equations.
+
+
+First, read the data from @HW97.
+We also convert J to cal and MPa to bar.
+
+```{r Cpdat}
+file <- system.file("extdata/cpetc/Cp.CH4.HW97.csv", package="CHNOSZ")
+Cpdat <- read.csv(file)
+Cpdat$Cp <- convert(Cpdat$Cp, "cal")
+Cpdat$P <- convert(Cpdat$P, "bar")
+```
+
+```{r EOSregress_hide, fig.margin = TRUE, fig.cap = "Heat Capacity of Aqueous Methane", fig.width=3.5, fig.height=3.5, cache=TRUE, results="hide", message=FALSE, echo=FALSE, dpi=50, out.width=672, out.height=336}
+var <- c("invTTheta2", "TXBorn")
+Cplm <- EOSregress(Cpdat, var, T.max=600)
+EOSplot(Cpdat, coefficients=round(Cplm$coefficients, 1))
+EOSplot(Cpdat, add=TRUE, lty=3)
+```
+
+Next, we specify the terms in the HKF equations and perform the regression, using data below 600 K.
+The terms in `Cplm` correspond to $c_1$, $c_2$ and $\omega$ in the HKF equations.
+We use `EOSplot` to plot the data and fitted line and show the coefficients in the legend.
+
+```{r EOSregress, eval=FALSE}
+var <- c("invTTheta2", "TXBorn")
+Cplm <- EOSregress(Cpdat, var, T.max=600)
+EOSplot(Cpdat, coefficients=round(Cplm$coefficients, 1))
+EOSplot(Cpdat, add=TRUE, lty=3)
+```
+
+`r tufte::margin_note("Be aware that the lines shown by EOSplot are calculated for a single pressure only, despite the temperature- and pressure-dependence of the data and regressions.")`
+
+Note the second call to `EOSplot`, where the coefficients are not provided.
+This causes the data to be re-regressed, this time without the $T$ < 600 K limitation.
+Although the overall fit, represented by the dotted line, is better, the fit to the low-temperature data is considerably worse.
+For nonelectroyltes, consider excluding extreme near-critical values from the regression (or modify the procedure to give them lower weight) in order to obtain more generally useful results.
+
+`EOScoeffs` is a small function to retrieve the HKF parameters from the database in CHNOSZ.
+The `stopifnot` expression tests that the regressed values are within 10% of the database values of $c_1$, $c_2$ and $\omega$ for $\mathrm{CH_4}$.
+
+```{r Cpcoeffs}
+Cpcoeffs_database <- EOScoeffs("CH4", "Cp")
+diff_coeffs <- Cplm$coefficients - Cpcoeffs_database
+stopifnot(all(abs(diff_coeffs/Cpcoeffs_database) < 0.1))
+Cpcoeffs_database
+```
+
+Compare the database values [from @SH90] with the regressed values in the legend of figure above.
+Some differences are expected as the values in the database are derived from earlier experimental data.
+
+## Setting the value of omega
+
+Given both high-temperature volumetric and calorimetric data for neutral species, the effective value of $\omega$ is most reliably regressed from the latter [@SSW01].
+Let's regress volumetric data using a value of omega taken from the heat capacity regression.
+First, read the data from @HWM96.
+
+```{r Vdat}
+file <- system.file("extdata/cpetc/V.CH4.HWM96.csv", package="CHNOSZ")
+Vdat <- read.csv(file)
+Vdat$P <- convert(Vdat$P, "bar")
+```
+
+Compressibilities of species (measured or estimated) are implied by the full set of HKF volumetric parameters ($a_1$, $a_2$, $a_3$, $a_4$).
+In this example we model volumes at nearly constant $P$.
+Therefore, we can use a simpler equation for $V^\circ$ written in terms of the "isobaric fit parameters" (Tanger and Helgeson 1988, p. 35) $\sigma$ and $\xi$, together with the solvation contribution that depends on the $Q$ Born function:
+$$V^\circ = \sigma + \frac{\xi}{T-\Theta} - {\omega}Q$$
+
+`r tufte::margin_note("EOSvar actually returns the negative of $Q$, so the omega symbol here carries no negative sign.")`
+
+Now we calculate the $Q$ Born function using `EOSvar` and multiply by $\omega$ (from the heat capacity regression) to get the solvation volume at each experimental temperature and pressure.
+Subtract the solvation volume from the experimental volumes and create a new data frame holding the calculated "non-solvation" volume.
+
+`r tufte::margin_note("Because we are dealing with volumes, the units of $ω$ are converted according to 1 cal = 41.84 cm$^3$∙bar.")`
+
+```{r Vdat_non}
+QBorn <- EOSvar("QBorn", T=Vdat$T, P=Vdat$P)
+Vomega <- convert(Cplm$coefficients[["TXBorn"]], "cm3bar")
+V_sol <- Vomega * QBorn
+V_non <- Vdat$V - V_sol
+Vdat_non <- data.frame(T=Vdat$T, P=Vdat$P, V=V_non)
+```
+
+Next, regress the non-solvation volume using the non-solvation terms in the HKF model.
+As with $C_P^\circ$, also get the values of the parameters from the database for comparison with the regression results.
+
+
+```{r Vdat_non_regress}
+var <- "invTTheta"
+Vnonlm <- EOSregress(Vdat_non, var, T.max=450)
+Vcoeffs <- round(c(Vnonlm$coefficients, QBorn=Vomega), 1)
+Vcoeffs_database <- EOScoeffs("CH4", "V")
+```
+
+```{r Vplot_hide, fig.margin=TRUE, results="hide", message=FALSE, echo=FALSE, fig.width=3.5, fig.height=7, fig.cap="Volume of Aqueous Methane", dpi=50, out.width=672, out.height=672}
+par(mfrow=c(2, 1))
+# plot 1
+EOSplot(Vdat, coefficients=Vcoeffs)
+EOSplot(Vdat, coefficients=Vcoeffs_database, add=TRUE, lty=2)
+# plot 2
+EOSplot(Vdat, coefficients=Vcoeffs_database, T.plot=600, lty=2)
+EOSplot(Vdat, coefficients=Vcoeffs, add=TRUE)
+```
+
+Finally, plot the data and regressions.
+The first plot shows all the data, and the second the low-temperature subset ($T$ < 600 K).
+The solid line is the two-term fit for $\sigma$ and $\xi$ (using $\omega$ from the heat capacity regression), and the dashed line shows the volumes calculated using the parameters in the database.
+The plot legends give the parameters from the two-term fit (first plot), or from the database (second plot).
+
+```{r Vplot, eval=FALSE}
+par(mfrow=c(2, 1))
+# plot 1
+EOSplot(Vdat, coefficients=Vcoeffs)
+EOSplot(Vdat, coefficients=Vcoeffs_database, add=TRUE, lty=2)
+# plot 2
+EOSplot(Vdat, coefficients=Vcoeffs_database, T.plot=600, lty=2)
+EOSplot(Vdat, coefficients=Vcoeffs, add=TRUE)
+```
+
+The equation for $V^\circ$ provides a reasonable approximation of the trend of lowest-temperature data ($T$ < 450 K).
+However, the equation does not closely reproduce the trend of higher-temperature $V^\circ$ data ($T$ < 600 K), nor behavior in the critical region.
+Because of these issues, some researchers are exploring alternatives to the HKF model for aqueous nonelectrolytes.
+(See also an example in `?EOSregress`.)
+
+
+# An example for charged species
+
+For this example, let's generate synthetic data for Na$^+$ using its parameters in the database.
+In the call to `subcrt` below, `convert=FALSE` means to take $T$ in units of K.
+
+```{r Nadat}
+T <- convert(seq(0, 600, 50), "K")
+P <- 1000
+prop.PT <- subcrt("Na+", T=T, P=P, grid="T", convert=FALSE)$out[[1]]
+Nadat <- prop.PT[, c("T", "P", "Cp")]
+```
+
+As noted above, $\omega$ for electrolytes is not a constant.
+What happens if we apply the linear model anyway, knowing it's not applicable (especially at high temperature)?
+
+```{r Nalm, fig.margin=TRUE, fig.width=3.5, fig.height=3.5, fig.cap="Heat Capacity of Na$^+$ (Inapplicable: Constant $ω$)", dpi=50, out.width=672, out.height=336}
+var <- c("invTTheta2", "TXBorn")
+Nalm <- EOSregress(Nadat, var, T.max=600)
+EOSplot(Nadat, coefficients=Nalm$coefficients, fun.legend=NULL)
+EOSplot(Nadat, add=TRUE, lty=3)
+```
+
+As before, the solid line is a fit to relatively low-temperature ($T$ < 600 K) data, and the dotted line a fit to the entire temperature range of the data.
+The fits using constant $\omega$ are clearly not acceptable.
+
+There is, however, a way out.
+A different variable, `Cp_s_var`, can be used to specify the calculation of the "solvation" heat capacity in the HKF equations **using the temperature- and pressure-dependent corrections for charged species**.
+To use this variable, the values of $\omega_{P_r,T_r}$ (omega at the reference temperature and pressure) and $Z$ (charge) must be given, in addition to $T$ and $P$.
+Of course, right now we **don't know** the value of $\omega_{P_r,T_r}$ -- it is the purpose of the regression to find it!
+But we can make a first guess using the value of $\omega$ found above.
+
+```{r Navars1}
+var1 <- c("invTTheta2", "Cp_s_var")
+omega.guess <- coef(Nalm)[3]
+```
+
+Then, we can use an iterative procedure that refines successive guesses of $\omega_{P_r,T_r}$.
+The convergence criterion is measured by the difference in sequential regressed values of $\omega$.
+
+```{r Nawhile, fig.margin=TRUE, fig.width=3.5, fig.height=3.5, fig.cap="Heat Capacity of Na$^+$ (Variable $ω$)", dpi=50, out.width=672, out.height=336}
+diff.omega <- 999
+while(abs(diff.omega) > 1) {
+ Nalm1 <- EOSregress(Nadat, var1, omega.PrTr=tail(omega.guess, 1), Z=1)
+ omega.guess <- c(omega.guess, coef(Nalm1)[3])
+ diff.omega <- tail(diff(omega.guess), 1)
+}
+EOSplot(Nadat, coefficients=signif(coef(Nalm1), 6),
+ omega.PrTr=tail(omega.guess, 1), Z=1)
+EOScoeffs("Na+", "Cp")
+```
+
+Alrighty! We managed to obtain EOS parameters from synthetic data for Na$^+$.
+The regressed values of the EOS parameters (shown in the plot legend) are very close to the database values (printed by the call to `EOScoeffs`) used to generate the synthetic data.
+
+## Doing it for volume
+
+Just like above, but using synthetic $V^\circ$ data.
+Note that the regressed value of $\omega$ has volumetric units (cm$^3$∙bar/mol), while `omega.PrTr` is in caloric units (cal/mol).
+Compared to $C_P^\circ$, the regression of $V^\circ$ is very finicky.
+Given a starting guess of $\omega_{P_r,T_r}$ of 1400000 cm$^3$∙bar/mol, the iteration converges on 1394890 instead of the "true" database value of 1383230 (represented by dashed line in the plot).
+
+```{r NaVolume_hide, fig.margin=TRUE, fig.width=3.5, fig.height=3.5, fig.cap="Volume of Na$^+$ (Variable $ω$)", results="hide", message=FALSE, echo=FALSE, dpi=50, out.width=672, out.height=336}
+T <- convert(seq(0, 600, 25), "K")
+P <- 1000
+prop.PT <- subcrt("Na+", T=T, P=P, grid="T", convert=FALSE)$out[[1]]
+NaVdat <- prop.PT[, c("T", "P", "V")]
+var1 <- c("invTTheta", "V_s_var")
+omega.guess <- 1400000
+diff.omega <- 999
+while(abs(diff.omega) > 1) {
+ NaVlm1 <- EOSregress(NaVdat, var1,
+ omega.PrTr=tail(convert(omega.guess, "calories"), 1), Z=1)
+ omega.guess <- c(omega.guess, coef(NaVlm1)[3])
+ diff.omega <- tail(diff(omega.guess), 1)
+}
+EOSplot(NaVdat, coefficients=signif(coef(NaVlm1), 6),
+ omega.PrTr=tail(convert(omega.guess, "calories"), 1), Z=1,
+ fun.legend="bottomleft")
+coefficients <- EOScoeffs("Na+", "V", P=1000)
+names(coefficients)[3] <- "V_s_var"
+EOSplot(NaVdat, coefficients=coefficients, Z=1, add=TRUE, lty=2,
+ omega.PrTr=convert(coefficients["V_s_var"], "calories"))
+```
+
+```{r NaVolume, eval=FALSE}
+T <- convert(seq(0, 600, 25), "K")
+P <- 1000
+prop.PT <- subcrt("Na+", T=T, P=P, grid="T", convert=FALSE)$out[[1]]
+NaVdat <- prop.PT[, c("T", "P", "V")]
+var1 <- c("invTTheta", "V_s_var")
+omega.guess <- 1400000
+diff.omega <- 999
+while(abs(diff.omega) > 1) {
+ NaVlm1 <- EOSregress(NaVdat, var1,
+ omega.PrTr=tail(convert(omega.guess, "calories"), 1), Z=1)
+ omega.guess <- c(omega.guess, coef(NaVlm1)[3])
+ diff.omega <- tail(diff(omega.guess), 1)
+}
+EOSplot(NaVdat, coefficients=signif(coef(NaVlm1), 6),
+ omega.PrTr=tail(convert(omega.guess, "calories"), 1), Z=1,
+ fun.legend="bottomleft")
+coefficients <- EOScoeffs("Na+", "V", P=1000)
+names(coefficients)[3] <- "V_s_var"
+EOSplot(NaVdat, coefficients=coefficients, Z=1, add=TRUE, lty=2,
+ omega.PrTr=convert(coefficients["V_s_var"], "calories"))
+```
+
+# Other possibilities
+
+These functions are limited to treatment of calorimetric and volumetric data.
+Other software tools have been described recently for deriving HKF parameters and other thermodynamic properties from different types of experimental data [@MKDW15; @Shv15].
+
Added: pkg/CHNOSZ/vignettes/EOSregress.bib
===================================================================
--- pkg/CHNOSZ/vignettes/EOSregress.bib (rev 0)
+++ pkg/CHNOSZ/vignettes/EOSregress.bib 2016-11-01 05:15:29 UTC (rev 114)
@@ -0,0 +1,108 @@
+% Encoding: UTF-8
+
+ at Article{HKF81,
+ author = {Helgeson, Harold C. and Kirkham, David H. and Flowers, George C.},
+ journal = {American Journal of Science},
+ title = {{T}heoretical prediction of the thermodynamic behavior of aqueous electrolytes at high pressures and temperatures: {IV}. {C}alculation of activity coefficients, osmotic coefficients, and apparent molal and standard and relative partial molal properties to 600°{C} and 5 {K}b},
+ year = {1981},
+ volume = {281},
+ number = {10},
+ pages = {1249--1516},
+ doi = {10.2475/ajs.281.10.1249},
+ size = {268 p.},
+}
+
+ at Article{SSW01,
+ author = {Schulte, Mitchell D. and Shock, Everett L. and Wood, Robert H.},
+ journal = {Geochimica et Cosmochimica Acta},
+ title = {{T}he temperature dependence of the standard-state thermodynamic properties of aqueous nonelectrolytes},
+ year = {2001},
+ volume = {65},
+ number = {21},
+ pages = {3919--3930},
+ doi = {10.1016/S0016-7037(01)00717-7},
+}
+
+ at Article{SH90,
+ author = {Shock, Everett L. and Helgeson, Harold C.},
+ journal = {Geochimica et Cosmochimica Acta},
+ title = {{C}alculation of the thermodynamic and transport properties of aqueous species at high pressures and temperatures: {S}tandard partial molal properties of organic species},
+ year = {1990},
+ volume = {54},
+ number = {4},
+ pages = {915 -- 945},
+ doi = {10.1016/0016-7037(90)90429-O},
+}
+
+ at Article{SOJSH92,
+ author = {Shock, Everett L. and Oelkers, Eric H. and Johnson, James W. and Sverjensky, Dimitri A. and Helgeson, Harold C.},
+ journal = {Journal of the Chemical Society, Faraday Transactions},
+ title = {{C}alculation of the thermodynamic properties of aqueous species at high pressures and temperatures: {E}ffective electrostatic radii, dissociation constants, and standard partial molal properties to 1000 °{C} and 5 kbar},
+ year = {1992},
+ volume = {88},
+ number = {6},
+ pages = {803 -- 826},
+ doi = {10.1039/FT9928800803},
+}
+
+ at Article{TH88,
+ author = {Tanger, IV, John C. and Helgeson, Harold C.},
+ journal = {American Journal of Science},
+ title = {{C}alculation of the thermodynamic and transport properties of aqueous species at high pressures and temperatures: {R}evised equations of state for the standard partial molal properties of ions and electrolytes},
+ year = {1988},
+ volume = {288},
+ number = {1},
+ pages = {19--98},
+ doi = {10.2475/ajs.288.1.19},
+}
+
+ at Article{HW97,
+ author = {Hn\v{e}dkovsk\'y, Lubom\'ir and Wood, Robert H.},
+ journal = {Journal of Chemical Thermodynamics},
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/chnosz -r 114
More information about the CHNOSZ-commits
mailing list