[CHNOSZ-commits] r43 - in pkg/CHNOSZ: . R inst inst/doc inst/tests man vignettes
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Feb 13 15:28:35 CET 2013
Author: jedick
Date: 2013-02-13 15:28:35 +0100 (Wed, 13 Feb 2013)
New Revision: 43
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/R/EOSregress.R
pkg/CHNOSZ/inst/NEWS
pkg/CHNOSZ/inst/doc/wjd.pdf
pkg/CHNOSZ/inst/tests/test-EOSregress.R
pkg/CHNOSZ/man/EOSregress.Rd
pkg/CHNOSZ/vignettes/wjd.Rnw
pkg/CHNOSZ/vignettes/wjd.lyx
Log:
EOSvar() finds user-defined variables
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2013-02-12 22:35:25 UTC (rev 42)
+++ pkg/CHNOSZ/DESCRIPTION 2013-02-13 14:28:35 UTC (rev 43)
@@ -1,6 +1,6 @@
Date: 2013-02-13
Package: CHNOSZ
-Version: 0.9-9.5
+Version: 0.9-9.6
Title: Chemical Thermodynamics and Activity Diagrams
Author: Jeffrey M. Dick
Maintainer: Jeffrey M. Dick <jmdick at asu.edu>
Modified: pkg/CHNOSZ/R/EOSregress.R
===================================================================
--- pkg/CHNOSZ/R/EOSregress.R 2013-02-12 22:35:25 UTC (rev 42)
+++ pkg/CHNOSZ/R/EOSregress.R 2013-02-13 14:28:35 UTC (rev 43)
@@ -19,9 +19,18 @@
"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],
- # the "default": get a variable that is a property of water
- (if(var %in% water.props()) water(var, T, P)[, 1]
- else stop(paste("can't find a variable named", var)))
+ # 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.props()) water(var, T, P)[, 1]
+ else if(exists(var)) {
+ if(is.function(get(var))) {
+ if(identical(names(formals(get(var))), c("T", "P"))) get(var)(T, P)
+ else stop(paste("the arguments of ", var, "() are not T, P", sep=""))
+ }
+ else stop(paste("an object named", var, "is not a function"))
+ }
+ else stop(paste("can't find a variable named", var))
+ )
)
return(out)
}
@@ -29,24 +38,27 @@
EOSlab <- function(var,coeff="") {
# make pretty labels for the variables
lab <- switch(EXPR = var,
- "(Intercept)" = substitute(YYY*" ",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)),
- "T" = substitute(YYY%*%italic(XXX),list(XXX=var,YYY=coeff)),
- "P" = substitute(YYY%*%italic(XXX),list(XXX=var,YYY=coeff)),
- "V" = substitute(YYY%*%italic(XXX),list(XXX=var,YYY=coeff)),
- "E" = substitute(YYY%*%italic(XXX),list(XXX=var,YYY=coeff)),
- "XBorn" = substitute(YYY%*%italic(XXX),list(XXX=var,YYY=coeff)),
- "QBorn" = substitute(YYY%*%italic(XXX),list(XXX=var,YYY=coeff)),
- "TXBorn" = substitute(YYY%*%italic(XXX),list(XXX=var,YYY=coeff)),
- "kT" = substitute(YYY%*%kappa[italic(T)],list(YYY=coeff)),
- "alpha" = substitute(YYY%*%alpha,list(YYY=coeff)),
- "beta" = substitute(YYY%*%beta,list(YYY=coeff)),
- "drho.dT" = substitute(YYY%*%(d~rho/dT),list(YYY=coeff)),
- "V.kT" = substitute(YYY%*%V~kappa[italic(T)],list(YYY=coeff)),
- NA
+ # 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)),
+ "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)),
+ # the rest are properties of water listed in water.Rd
+ "V" = substitute(YYY%*%italic(V), list(YYY=coeff)),
+ "E" = substitute(YYY%*%italic(E), list(YYY=coeff)),
+ "kT" = substitute(YYY%*%kappa[italic(T)], list(YYY=coeff)),
+ "alpha" = substitute(YYY%*%alpha, list(YYY=coeff)),
+ "beta" = substitute(YYY%*%beta, list(YYY=coeff)),
+ "XBorn" = substitute(YYY%*%italic(X), list(YYY=coeff)),
+ "QBorn" = substitute(YYY%*%italic(Q), list(YYY=coeff)),
+ # fallback, use the name of the variable (may be the name of a user-defined function)
+ substitute(YYY%*%italic(XXX), list(YYY=coeff, XXX=var))
)
return(lab)
}
@@ -84,7 +96,7 @@
}
EOSplot <- function(exptdata,var=NULL,T.max=9999,T.plot=NULL,
- P=NULL,fun.legend="topleft",coefficients=NULL) {
+ fun.legend="topleft",coefficients=NULL) {
# plot experimental and modelled volumes and heat capacities
# first figure out the property (Cp or V) from the exptdata
prop <- colnames(exptdata)[3]
@@ -93,71 +105,52 @@
if(prop=="Cp") var <- c("invTTheta2","TXBorn")
if(prop=="V") var <- c("invTTheta","QBorn")
}
- expt <- exptdata
# perform the regression, only using temperatures up to T.max
if(is.null(coefficients)) {
- EOSlm <- EOSregress(expt,var,T.max)
+ EOSlm <- EOSregress(exptdata, var, T.max)
coefficients <- EOSlm$coefficients
}
# only plot points below a certain temperature
- iexpt <- 1:nrow(expt)
- if(!is.null(T.plot)) iexpt <- which(expt$T < T.plot)
- iX <- match(prop,colnames(expt))
- ylim <- extendrange(expt[iexpt,iX],f=0.1)
- xlim <- extendrange(expt$T[iexpt],f=0.1)
+ 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)
+ xlim <- extendrange(exptdata$T[iexpt], f=0.1)
# 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)
- # we group the data by pressure ranges;
- # assume increasing temperatures are in the
- # same pressure range but a decrease in temperature
- # signals the next pressure range
- idrop <- c(1,which(diff(expt$T)<0)+1,length(expt$T)+1)
- Plab <- character()
- pch.open <- c(1,0,2)
- pch.filled <- c(16,15,17)
- for(i in 1:(length(idrop)-1)) {
- ip <- idrop[i]:(idrop[i+1]-1)
- # find the calculated values at these conditions
- myT <- expt$T[ip]
- myP <- expt$P[ip]
- calc.X <- EOScalc(coefficients,myT,myP)
- expt.X <- expt[ip,iX]
- # are we within 10% of the values
- in10 <- which(abs((calc.X-expt.X)/expt.X) < 0.1)
- pch <- rep(pch.open[i],length(myT))
- pch[in10] <- pch.filled[i]
- points(myT,expt[ip,iX],pch=pch)
- # if we calculate lines at a constant P, do that
- xs <- seq(xlim[1],xlim[2],length.out=200)
- if(!is.null(P)) {
- myT <- xs
- myP <- P
- calc.X <- EOScalc(coefficients,myT,myP)
- }
- # take out NAs and Infinite values
- iNA <- is.na(calc.X) | is.infinite(calc.X)
- xs <- xs[!iNA]
- calc.X <- calc.X[!iNA]
- myT <- myT[!iNA]
- # plot regression line
- lines(xs,splinefun(myT,calc.X,method="monoH.FC")(xs))
- Plim <- range(expt$P[ip])
- Plab <- c(Plab,paste(Plim[1],"-",Plim[2],"bar"))
- }
+ 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)
+ # 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)
# make legend
if(!is.null(fun.legend)) {
- coeffs <- as.character(round(as.numeric(coefficients),4))
+ 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="")
+ 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)
+ legend(fun.legend, legend=fun.lab, pt.cex=0.1)
}
- return(invisible(list(xlim=range(expt$T[iexpt]))))
+ return(invisible(list(xlim=range(exptdata$T[iexpt]))))
}
EOScoeffs <- function(species, property) {
Modified: pkg/CHNOSZ/inst/NEWS
===================================================================
--- pkg/CHNOSZ/inst/NEWS 2013-02-12 22:35:25 UTC (rev 42)
+++ pkg/CHNOSZ/inst/NEWS 2013-02-13 14:28:35 UTC (rev 43)
@@ -1,4 +1,4 @@
-CHANGES IN CHNOSZ 0.9-9.5 (2013-02-13)
+CHANGES IN CHNOSZ 0.9-9.6 (2013-02-13)
--------------------------------------
- Fix calculation of free energy derivative in wjd().
@@ -45,9 +45,16 @@
- Separate rho.IAPWS95() from water.IAPWS95().
-- Any property of water can be used as a variable in EOSvar().
+- In addition to the original regression variables, EOSvar() recognizes
+ names for available properties in water(), or can use the name to
+ get a user-defined function of temperature and pressure.
+- Simplify EOSplot() somewhat (don't group data by pressure ranges).
+- [temporary] Deactivate code using 'central' method in guess()
+ (wjd.Rnw) as limSolve package on not available on R-Forge for windows.
+
+
CHANGES IN CHNOSZ 0.9-9 (2013-01-01)
------------------------------------
Modified: pkg/CHNOSZ/inst/doc/wjd.pdf
===================================================================
(Binary files differ)
Modified: pkg/CHNOSZ/inst/tests/test-EOSregress.R
===================================================================
--- pkg/CHNOSZ/inst/tests/test-EOSregress.R 2013-02-12 22:35:25 UTC (rev 42)
+++ pkg/CHNOSZ/inst/tests/test-EOSregress.R 2013-02-13 14:28:35 UTC (rev 43)
@@ -2,5 +2,43 @@
test_that("EOSvar stops with unknown variables", {
expect_error(EOSvar("TX", T=25, P=1), "can't find a variable named TX")
+ # 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")
})
+test_that("regressions return known HKF parameters", {
+ # 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)
+ # convert=FALSE means that temperature has units of K
+ CH4.prop <- subcrt("CH4", T=T, P=P, grid="T", convert=FALSE)$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 calories (divide by 41.84)
+ V.coeff <- convert(V.lm$coefficients, "calories")
+ ## the tests: did we get the HKF parameters that are in the database?
+ CH4.par <- info(info("CH4"))
+ # c1 and c2
+ expect_equal(Cp.coeff[1], CH4.par$c1, check.attr=FALSE)
+ expect_equal(Cp.coeff[2], CH4.par$c2, check.attr=FALSE)
+ # omega (from Cp)
+ expect_equal(Cp.coeff[3], CH4.par$omega, check.attr=FALSE)
+ # a1, a2, a3 and a4
+ expect_equal(V.coeff[1], CH4.par$a1, check.attr=FALSE)
+ expect_equal(V.coeff[2], CH4.par$a2, check.attr=FALSE)
+ expect_equal(V.coeff[3], CH4.par$a3, check.attr=FALSE)
+ expect_equal(V.coeff[4], CH4.par$a4, check.attr=FALSE)
+ # omega (from V) - note negative sign
+ expect_equal(-V.coeff[5], CH4.par$omega, check.attr=FALSE)
+})
Modified: pkg/CHNOSZ/man/EOSregress.Rd
===================================================================
--- pkg/CHNOSZ/man/EOSregress.Rd 2013-02-12 22:35:25 UTC (rev 42)
+++ pkg/CHNOSZ/man/EOSregress.Rd 2013-02-13 14:28:35 UTC (rev 43)
@@ -8,7 +8,8 @@
\alias{EOScoeffs}
\title{Regress Equations-of-State Parameters for Aqueous Species}
\description{
- Functions for fitting experimental volume and heat capacities using regression equations. Possible models include the Helgeson-Kirkham-Flowers (HKF) equations of state and other equations defined using any combination of terms derived from the temperature, pressure and thermodynamic and electrostatic properties of water.
+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.
}
\usage{
@@ -16,7 +17,7 @@
EOSvar(var, T, P)
EOScalc(coefficients, T, P)
EOSplot(exptdata, var = NULL, T.max = 9999, T.plot = NULL,
- P = NULL, fun.legend = "topleft", coefficients = NULL)
+ fun.legend = "topleft", coefficients = NULL)
EOSlab(var, coeff = "")
EOScoeffs(species, property)
}
@@ -36,7 +37,10 @@
}
\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). 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}; these names correspond to variables identified below:
+\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).
+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).
\tabular{ll}{
\code{T} \tab \eqn{T}{T} (temperature) \cr
@@ -45,32 +49,33 @@
\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)} \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{V} \tab \eqn{V}{V} (volume of water) \cr
- \code{E} \tab \eqn{E}{E} (isobaric expansivity of water) \cr
- \code{kT} \tab \eqn{\kappa_T}{kT} (isothermal compressibility of water) \cr
- \code{alpha} \tab \eqn{\alpha}{alpha} (coefficient of isobaric expansivity of water) \cr
- \code{beta} \tab \eqn{\beta}{beta} (coefficients of isothermal compressibility of water) \cr
- \code{X} \tab \eqn{X}{X} (Born function \eqn{X}{X}) \cr
- \code{Q} \tab \eqn{Q}{Q} (Born function \eqn{Q}{Q}) \cr
- \code{TX} \tab \eqn{TX}{TX} (temperature times \eqn{X}{X}) \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} takes as input \code{var} (one of the names of the variables listed above), and \code{T} (temperature in degrees Kelvin), \code{P} (pressure in bars). It returns the value of the variable at the specified temperature-pressure condition(s). 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{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{EOSplot} takes a table of data in \code{exptdata}, runs \code{EOSregress} and \code{EOSpred} 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{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{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 (using \code{\link{substitute}}) to the name of the variable.
+\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{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 dataframe 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 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}.
+\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 (using \code{\link{substitute}}) to the name of the variable.
- The motivation for writing these functions is 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).
+\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}.
+The motivation for writing these functions is 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).
+
}
\value{
@@ -88,40 +93,7 @@
\examples{
\dontshow{data(thermo)}
-## regress calculated heat capacities and volumes of CH4(aq)
-## to test that regressions return known HKF parameters
-# calculate the properties of CH4 on a T-P grid
-T <- convert(seq(0, 350, 25), "K")
-P <- seq(200, 1000, 100)
-# convert=FALSE means that temperature has units of K
-CH4.prop <- subcrt("CH4", T=T, P=P, grid="T", convert=FALSE)$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 calories (divide by 41.84)
-V.coeff <- convert(V.lm$coefficients, "calories")
-## the tests: did we get the HKF parameters that are in the database?
-CH4.par <- info(info("CH4"))
-# c1 and c2
-stopifnot(all.equal(Cp.coeff[1], CH4.par$c1, check.attr=FALSE))
-stopifnot(all.equal(Cp.coeff[2], CH4.par$c2, check.attr=FALSE))
-# omega (from Cp)
-stopifnot(all.equal(Cp.coeff[3], CH4.par$omega, check.attr=FALSE))
-# a1, a2, a3 and a4
-stopifnot(all.equal(V.coeff[1], CH4.par$a1, check.attr=FALSE))
-stopifnot(all.equal(V.coeff[2], CH4.par$a2, check.attr=FALSE))
-stopifnot(all.equal(V.coeff[3], CH4.par$a3, check.attr=FALSE))
-stopifnot(all.equal(V.coeff[4], CH4.par$a4, check.attr=FALSE))
-# omega (from V) - note negative sign
-stopifnot(all.equal(-V.coeff[5], CH4.par$omega, check.attr=FALSE))
-
-## regress experimental heat capacities of CH4
+## 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/Cp.CH4.HW97.csv", package="CHNOSZ")
@@ -151,6 +123,16 @@
EOSplot(d, coefficients=CH4coeffs, T.plot=600)
title("Cp from EOS parameters in database")
+# following from above, with user-defined variables
+invTTTheta3 <- function(T, P) (2*T)/(T-T*thermo$opt$Theta)^3
+invTX <- function(T, P) 1/T*water("XBorn", T=T, P=P)[,1]
+EOSvar("invTTTheta3", d$T, d$P)
+var <- c("invTTTheta3", "invTX")
+EOSregress(d, var)
+# the plot is commented for this toy example
+#EOSplot(d, var)
+
+
## model experimental volumes of CH4
## using HKF equation and an exploratory one
f <- system.file("extdata/cpetc/V.CH4.HWM96.csv", package="CHNOSZ")
Modified: pkg/CHNOSZ/vignettes/wjd.Rnw
===================================================================
--- pkg/CHNOSZ/vignettes/wjd.Rnw 2013-02-12 22:35:25 UTC (rev 42)
+++ pkg/CHNOSZ/vignettes/wjd.Rnw 2013-02-13 14:28:35 UTC (rev 43)
@@ -265,8 +265,8 @@
as.list(args(wjd))$Y
Y1 <- guess(method="stoich")
Y1
-Y2 <- guess(method="central")
-Y2
+#Y2 <- guess(method="central")
+#Y2
@
@@ -285,30 +285,15 @@
were \Sexpr{niterY1-niter} more iterations than using the default
initial solution in \texttt{wjd()}.
-\setkeys{Gin}{width=0.6\textwidth}
-<<central_guess>>=
-wY2 <- wjd(Y=Y2)
-niterY2 <- length(wY2$lambda)
-niterY2
-is.near.equil(wY2, tol=0.0001)
-@
-\setkeys{Gin}{width=0.6\textwidth}
-The initial guess generated using the ``central'' method lead to
-\Sexpr{niterY2-niter} more iterations than encountered using the
-argument defaults in \texttt{wjd()}, and the algorithm converged to
-an even lower variability in chemical potentials of the elements than
-using the ``stoich'' method. The differences in convergence could
-be coincidental, and the ``stoich'' method might be preferable for
-general usage because it allows multiple guesses to be tested automatically
-(see below).
+
Do the different initial guesses actually give similar results?
\setkeys{Gin}{width=0.6\textwidth}
<<compare_guesses, fig=TRUE>>=
plot(1:10, w$X, xlab="species", ylab="mole fraction")
points(1:10, wY1$X, pch=0)
-points(1:10, wY2$X, pch=2)
+#points(1:10, wY2$X, pch=2)
@
\setkeys{Gin}{width=0.6\textwidth}
Modified: pkg/CHNOSZ/vignettes/wjd.lyx
===================================================================
--- pkg/CHNOSZ/vignettes/wjd.lyx 2013-02-12 22:35:25 UTC (rev 42)
+++ pkg/CHNOSZ/vignettes/wjd.lyx 2013-02-13 14:28:35 UTC (rev 43)
@@ -79,6 +79,11 @@
\filename_suffix 0
\color #000000
\end_branch
+\branch inactive
+\selected 0
+\filename_suffix 0
+\color #000000
+\end_branch
\index Index
\shortcut idx
\color #008000
@@ -1063,12 +1068,12 @@
\begin_layout Chunk
-Y2 <- guess(method="central")
+#Y2 <- guess(method="central")
\end_layout
\begin_layout Chunk
-Y2
+#Y2
\end_layout
\begin_layout Chunk
@@ -1169,6 +1174,10 @@
\end_layout
\begin_layout Standard
+\begin_inset Branch inactive
+status collapsed
+
+\begin_layout Standard
\begin_inset Branch short
status open
@@ -1269,13 +1278,18 @@
guesses to be tested automatically (see below).
\end_layout
+\end_inset
+
+
+\end_layout
+
\begin_layout Standard
Do the different initial guesses actually give similar results?
\end_layout
\begin_layout Standard
\begin_inset Branch short
-status open
+status collapsed
\begin_layout Chunk
@@ -1303,7 +1317,7 @@
\begin_layout Chunk
-points(1:10, wY2$X, pch=2)
+#points(1:10, wY2$X, pch=2)
\end_layout
\begin_layout Chunk
More information about the CHNOSZ-commits
mailing list