[CHNOSZ-commits] r21 - in pkg/CHNOSZ: . R inst inst/tests man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sun Sep 23 17:11:26 CEST 2012
Author: jedick
Date: 2012-09-23 17:11:26 +0200 (Sun, 23 Sep 2012)
New Revision: 21
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/R/cgl.R
pkg/CHNOSZ/R/examples.R
pkg/CHNOSZ/R/hkf.R
pkg/CHNOSZ/R/util.data.R
pkg/CHNOSZ/R/water.R
pkg/CHNOSZ/inst/NEWS
pkg/CHNOSZ/inst/tests/test-thermo.R
pkg/CHNOSZ/man/eos.Rd
pkg/CHNOSZ/man/subcrt.Rd
pkg/CHNOSZ/man/util.data.Rd
Log:
water.SUPCRT92() and subcrt() return NA properties at Psat and temperature > critical temperature
mod.obigt() defaults to NA, not 0, for properties not provided for a new species
cgl() and hkf() revert to tabulated (constant) Cp and V if EoS parameters are NA
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2012-09-23 11:08:03 UTC (rev 20)
+++ pkg/CHNOSZ/DESCRIPTION 2012-09-23 15:11:26 UTC (rev 21)
@@ -1,4 +1,4 @@
-Date: 2012-09-20
+Date: 2012-09-23
Package: CHNOSZ
Version: 0.9-7.98
Title: Chemical Thermodynamics and Activity Diagrams
Modified: pkg/CHNOSZ/R/cgl.R
===================================================================
--- pkg/CHNOSZ/R/cgl.R 2012-09-23 11:08:03 UTC (rev 20)
+++ pkg/CHNOSZ/R/cgl.R 2012-09-23 15:11:26 UTC (rev 21)
@@ -20,29 +20,42 @@
w <- NULL
for(i in 1:length(prop)) {
property <- prop[i]
+ # a test for availability of the EoS parameters
+ # here we assume that the parameters are in the same position as in thermo$obigt
+ # leave T transition (in 20th column) alone
+ hasEOS <- any(!is.na(EOS[, 13:19]))
+ # if at least one of the EoS parameters is available, zero out any NA's in the rest
+ if(hasEOS) EOS[, 13:19][, is.na(EOS[, 13:19])] <- 0
# equations for lambda adapted from HOK+98
- if(property=='cp')
- p <- EOS$a + EOS$b*T + EOS$c*T^-2 + EOS$d*T^-0.5 + EOS$e*T^2 + EOS$f*T^EOS$lambda
- if(property=='v')
+ if(property=='cp') {
+ # use constant Cp if the EoS parameters are not available
+ if(!hasEOS) p <- EOS$Cp
+ else p <- EOS$a + EOS$b*T + EOS$c*T^-2 + EOS$d*T^-0.5 + EOS$e*T^2 + EOS$f*T^EOS$lambda
+ } else if(property=='v') {
p <- rep(EOS$V,length(T))
- if(property %in% c('e','kt')) {
+ } else if(property %in% c('e','kt')) {
p <- rep(NA,length(T))
warning('cgl: E and/or kT of cr, gas and/or liq species are NA.')
- }
- if(property=='g') {
+ } else if(property=='g') {
+ # use constant Cp if the EoS parameters are not available
+ if(!hasEOS) p <- GHS$G + EOS$Cp*(T-Tr-T*log(T/Tr)) else
+ # Gibbs energy integral: the value at Tref plus heat capacity terms
p <- GHS$G + EOS$a*(T-Tr-T*log(T/Tr)) -
EOS$b*(T-Tr)^2/2 - EOS$c*(1/T + T/Tr^2 - 2/Tr)/2 -
EOS$d*(T^0.5-0.5*T*Tr^-0.5-0.5*Tr^0.5)/-0.25 -
EOS$e*(T^3-3*T*Tr^2+2*Tr^3)/6
+ # entropy and volume terms
if(!is.na(GHS$S)) p <- p - GHS$S*(T-Tr)
if(!is.na(EOS$V)) p <- p + convert(EOS$V*(P-Pr),'calories')
+ # use additional heat capacity term if it's defined
if(!is.na(EOS$f) & !is.na(EOS$lambda)) if(EOS$f!=0) {
if(EOS$lambda== -1) p <- p + EOS$f*(log(T/Tr)-T*(1/Tr-1/T))
else p <- p + EOS$f * ( T^(EOS$lambda+1) - (EOS$lambda+1)*T*Tr^EOS$lambda +
EOS$lambda*Tr^(EOS$lambda+1) ) / ( EOS$lambda*(EOS$lambda+1) )
}
- }
- if(property=='h') {
+ } else if(property=='h') {
+ # use constant Cp if the EoS parameters are not available
+ if(!hasEOS) p <- EOS$Cp*(T-Tr) else
p <- EOS$a*(T-Tr) + EOS$b*(T^2-Tr^2)/2 +
EOS$c*(1/T-1/Tr)/-1 + EOS$d*(T^0.5-Tr^0.5)/0.5 +
EOS$e*(T^3-Tr^3)/3
@@ -53,8 +66,9 @@
if(EOS$lambda == -1) p <- p + EOS$f*log(T/Tr)
else p <- p - EOS$f * ( T^(EOS$lambda+1) - Tr^(EOS$lambda+1) ) / (EOS$lambda+1)
}
- }
- if(property=='s') {
+ } else if(property=='s') {
+ # use constant Cp if the EoS parameters are not available
+ if(!hasEOS) p <- EOS$Cp*log(T/Tr) else
p <- EOS$a*log(T/Tr) + EOS$b*(T-Tr) +
EOS$c*(T^(-2)-Tr^(-2))/(-2) + EOS$e*(T^2-Tr^2)/2 +
EOS$d*(T^-0.5-Tr^-0.5)/-0.5
Modified: pkg/CHNOSZ/R/examples.R
===================================================================
--- pkg/CHNOSZ/R/examples.R 2012-09-23 11:08:03 UTC (rev 20)
+++ pkg/CHNOSZ/R/examples.R 2012-09-23 15:11:26 UTC (rev 21)
@@ -98,6 +98,7 @@
P <- c(list("Psat"), as.list(seq(500, 4000, by=500)))
# for each of those what's the range of temperature
T <- list()
+ # T > 350 degC at Psat is possibly inappropriate; see "Warning" of subcrt.Rd
T[[1]] <- seq(0, 370, 5)
T[[2]] <- seq(265, 465, 5)
T[[3]] <- seq(285, 760, 5)
@@ -113,8 +114,6 @@
Texpt <- expt$T[iexpt]
logK <- c(logK, splinefun(s$out$T, s$out$logK)(Texpt))
}
- # FIXME: if there are zig-zags in the lines, those
- # are discontinuities in the EoS code
legend("bottomleft",pch=unique(expt$pch),
legend=c(unique(expt$source),tail(expt$source,1)))
title(main=paste("NaCl(aq) = Na+ + Cl-\n",
Modified: pkg/CHNOSZ/R/hkf.R
===================================================================
--- pkg/CHNOSZ/R/hkf.R 2012-09-23 11:08:03 UTC (rev 20)
+++ pkg/CHNOSZ/R/hkf.R 2012-09-23 15:11:26 UTC (rev 21)
@@ -45,6 +45,16 @@
# loop over each species
GHS <- ghs[k,]
EOS <- eos[k,]
+ # substitute Cp and V for missing EoS parameters
+ # here we assume that the parameters are in the same position as in thermo$obigt
+ # put the heat capacity in for c1 if both c1 and c2 are missing
+ if(all(is.na(EOS[, 17:18]))) EOS[, 17] <- EOS$Cp
+ # put the volume in for a1 if a1, a2, a3 and a4 are missing
+ if(all(is.na(EOS[, 13:16]))) EOS[, 13] <- convert(EOS$V, "calories")
+ # test for availability of the EoS parameters
+ hasEOS <- any(!is.na(EOS[, 13:20]))
+ # if at least one of the EoS parameters is available, zero out any NA's in the rest
+ if(hasEOS) EOS[, 13:20][, is.na(EOS[, 13:20])] <- 0
# compute values of omega(P,T) from those of omega(Pr,Tr)
# using g function etc. (Shock et al., 1992 and others)
omega <- EOS$omega # omega.PrTr
@@ -82,31 +92,33 @@
p.a <- EOS$a1*(P-Pr) + EOS$a2*log((Psi+P)/(Psi+Pr)) +
((2*T-Theta)/(T-Theta)^2)*(EOS$a3*(P-Pr)+EOS$a4*log((Psi+P)/(Psi+Pr)))
p <- p.c + p.a
- }
- if(prop=="s") {
+ } else if(prop=="s") {
p.c <- EOS$c1*log(T/Tr) -
(EOS$c2/Theta)*( 1/(T-Theta)-1/(Tr-Theta) +
log( (Tr*(T-Theta))/(T*(Tr-Theta)) )/Theta )
p.a <- (T-Theta)^(-2)*(EOS$a3*(P-Pr)+EOS$a4*log((Psi+P)/(Psi+Pr)))
p <- p.c + p.a
- }
- if(prop=="g") {
+ } else if(prop=="g") {
p.c <- -EOS$c1*(T*log(T/Tr)-T+Tr) -
EOS$c2*( (1/(T-Theta)-1/(Tr-Theta))*((Theta-T)/Theta) -
(T/Theta^2)*log((Tr*(T-Theta))/(T*(Tr-Theta))) )
p.a <- EOS$a1*(P-Pr) + EOS$a2*log((Psi+P)/(Psi+Pr)) +
(EOS$a3*(P-Pr) + EOS$a4*log((Psi+P)/(Psi+Pr)))/(T-Theta)
p <- p.c + p.a
+ # nonsolvation cp v kt e equations
+ } else if(prop=='cp') {
+ p <- EOS$c1 + EOS$c2 * ( T - Theta ) ^ (-2)
+ } else if(prop=='v') {
+ p <- convert(EOS$a1,'cm3bar') +
+ convert(EOS$a2,'cm3bar') / ( Psi + P) +
+ ( convert(EOS$a3,'cm3bar') + convert(EOS$a4,'cm3bar') / ( Psi + P ) ) / ( T - Theta)
+ } else if(prop=='kt') {
+ p <- ( convert(EOS$a2,'cm3bar') +
+ convert(EOS$a4,'cm3bar') / (T - Theta) ) * (Psi + P) ^ (-2)
+ } else if(prop=='e') {
+ p <- convert( - ( EOS$a3 + EOS$a4 / convert((Psi + P),'calories') ) *
+ (T - Theta) ^ (-2),'cm3bar')
}
- # nonsolvation cp v kt e equations
- if(prop=='cp') p <- EOS$c1 + EOS$c2 * ( T - Theta ) ^ (-2)
- if(prop=='v') p <- convert(EOS$a1,'cm3bar') +
- convert(EOS$a2,'cm3bar') / ( Psi + P) +
- ( convert(EOS$a3,'cm3bar') + convert(EOS$a4,'cm3bar') / ( Psi + P ) ) / ( T - Theta)
- if(prop=='kt') p <- ( convert(EOS$a2,'cm3bar') +
- convert(EOS$a4,'cm3bar') / (T - Theta) ) * (Psi + P) ^ (-2)
- if(prop=='e') p <- convert( - ( EOS$a3 + EOS$a4 / convert((Psi + P),'calories') ) *
- (T - Theta) ^ (-2),'cm3bar')
}
if( icontrib=="s") {
# solvation ghs equations
@@ -164,7 +176,7 @@
out0 <- numeric(length(rhohat))
out <- list(g=out0, dgdT=out0, d2gdT2=out0, dgdP=out0)
# only rhohat less than 1 will give results other than zero
- idoit <- rhohat < 1
+ idoit <- rhohat < 1 & !is.na(rhohat)
rhohat <- rhohat[idoit]
Tc <- Tc[idoit]
P <- P[idoit]
@@ -196,6 +208,8 @@
( af2 * (1000 - P) ^ 3 + af3 * (1000 - P) ^ 4 )
# limits of the f function (region II of Fig. 6)
ifg <- Tc > 155 & P < 1000 & Tc < 355
+ # in case any T or P are NA
+ ifg <- ifg & !is.na(ifg)
# Eq. 32
g[ifg] <- g[ifg] - f[ifg]
## now we have g at P, T
Modified: pkg/CHNOSZ/R/util.data.R
===================================================================
--- pkg/CHNOSZ/R/util.data.R 2012-09-23 11:08:03 UTC (rev 20)
+++ pkg/CHNOSZ/R/util.data.R 2012-09-23 15:11:26 UTC (rev 21)
@@ -52,9 +52,8 @@
if(length(inew) > 0) {
# the right number of blank rows of thermo$obigt
newrows <- thermo$obigt[1:length(inew), ]
- # missing values are NA for character, 0 for numeric
- newrows[, 1:6] <- NA
- newrows[, 7:20] <- 0
+ # if we don't know something it's NA
+ newrows[] <- NA
# put in a default state
newrows$state <- thermo$opt$state
# fill in the columns
Modified: pkg/CHNOSZ/R/water.R
===================================================================
--- pkg/CHNOSZ/R/water.R 2012-09-23 11:08:03 UTC (rev 20)
+++ pkg/CHNOSZ/R/water.R 2012-09-23 15:11:26 UTC (rev 21)
@@ -694,23 +694,22 @@
for(i in 1:length(T)) {
states[1] <- Tc[i]
states[2] <- P[i]
- if(NaN %in% c(Tc[i],P[i])) {
+ if(any(is.na(c(Tc[i],P[i])))) {
+ # if T or P is NA, all properties are NA
w <- matrix(rep(NA,23),nrow=1)
w.out[i,] <- w
p.out[i] <- NA
rho.out[i] <- NA
- if(isat) err <- 0 else err <- 1
- err.out[i] <- err
} else {
inc <- 0
- t <- .Fortran('H2O92',as.integer(specs),as.double(states),
+ h2o <- .Fortran('H2O92',as.integer(specs),as.double(states),
as.double(rep(0,46)),as.integer(0),PACKAGE='CHNOSZ')
# errors
- err <- t[[4]]
+ err <- h2o[[4]]
err.out[i] <- err
# density
- rho <- t[[2]][3]
- rho2 <- t[[2]][4]
+ rho <- h2o[[2]][3]
+ rho2 <- h2o[[2]][4]
if(rho2 > rho) {
# liquid is denser than vapor
rho <- rho2
@@ -719,15 +718,16 @@
}
rho.out[i] <- rho
# most of the properties we're interested in
- w <- t(t[[3]][iprop+inc])
+ w <- t(h2o[[3]][iprop+inc])
if(err==1) w[1,] <- NA
# update the ith row of the output matrix
w.out[i,] <- w
# Psat
if(isat | 'psat' %in% tolower(property)) {
- p <- t[[2]][2]
- #if(T[i] < 373.124) p <- 1
- if(p < 1) p <- 1
+ p <- h2o[[2]][2]
+ p[p==0] <- NA
+ # Psat specifies P=1 below 100 degC
+ p[p < 1] <- 1
p.out[i] <- p
} else {
p.out[i] <- P[i]
@@ -754,8 +754,8 @@
if(isat) msgout(paste("water.SUPCRT92: error",plural2," calculating ",
nerr," of ",length(T)," point",plural,"; for Psat we need T < 647.067 K\n",sep=""))
else msgout(paste("water.SUPCRT92: error",plural2," calculating ",nerr,
- " of ",length(T)," point",plural,"; T and/or P are NA, ",
- "or T < Tfusion at P, T > 2250 degC, or P > 30kb.\n",sep=""))
+ " of ",length(T)," point",plural,
+ "; T < Tfusion at P, T > 2250 degC, or P > 30kb.\n",sep=""))
# that last bit is taken from SUP92D.f in the SUPCRT92 distribution
}
} else {
Modified: pkg/CHNOSZ/inst/NEWS
===================================================================
--- pkg/CHNOSZ/inst/NEWS 2012-09-23 11:08:03 UTC (rev 20)
+++ pkg/CHNOSZ/inst/NEWS 2012-09-23 15:11:26 UTC (rev 21)
@@ -1,4 +1,4 @@
-CHANGES IN CHNOSZ 0.9-7.98 (2012-09-20)
+CHANGES IN CHNOSZ 0.9-7.98 (2012-09-23)
---------------------------------------
SIGNIFICANT USER-VISIBLE CHANGES:
@@ -114,6 +114,9 @@
unlist()ing the results easier (used in element.mu() and
basis.logact()).
+- subcrt() now outputs NA values for properties at temperatures above
+ the critical temperature of H2O, when Psat is being used.
+
- New function i2A() for generating a stoichiometric matrix from indices
of species in the thermodynamic database.
Modified: pkg/CHNOSZ/inst/tests/test-thermo.R
===================================================================
--- pkg/CHNOSZ/inst/tests/test-thermo.R 2012-09-23 11:08:03 UTC (rev 20)
+++ pkg/CHNOSZ/inst/tests/test-thermo.R 2012-09-23 15:11:26 UTC (rev 21)
@@ -11,12 +11,12 @@
expect_equal(c.Ala <- info(info("[Ala]", "cr"))$c, 0)
# when we make a protein, its G depends on temperature
expect_true(all(diff(subcrt("LYSC_CHICK", "cr")$out[[1]]$G) < 0))
- # turn the value of c for [Ala](cr) into NA
- mod.obigt(name="[Ala]", state="cr", c=NA)
- # now when we make a protein, its G is NA
+ # turn the values of G and S for [Ala](cr) into NA
+ mod.obigt(name="[Ala]", state="cr", G=NA, S=NA)
+ # now when we make a protein(cr), its G is NA
expect_true(all(is.na(subcrt("RNAS1_BOVIN", "cr")$out[[1]]$G)))
# also check propagation of NA for aqueous species
- mod.obigt(name="[Ala]", state="aq", c=NA)
+ mod.obigt(name="[Ala]", state="aq", G=NA, S=NA)
expect_true(all(is.na(subcrt("[Ala]", "aq")$out[[1]]$G)))
# be nice and restore the database
suppressPackageStartupMessages(data(thermo))
@@ -27,10 +27,10 @@
expect_error(mod.obigt(list(name="test")))
# the default state is aq
expect_message(itest <- mod.obigt(list(name="test", date=today())), "added test\\(aq\\)")
- # we should get zero values of G for a species with zero properties
- expect_true(all(subcrt(itest)$out[[1]]$G==0))
- # a value for Cp does not integrate to G, but a value for c1 does
- # (this is probably unexpected behaviour and could be changed ...)
- expect_true(all(subcrt(mod.obigt(list(name="test", Cp=100)))$out[[1]]$G==0))
- expect_false(all(subcrt(mod.obigt(list(name="test", c1=100)))$out[[1]]$G==0))
+ # we should get NA values of G for a species with NA properties
+ expect_true(all(is.na(subcrt(itest)$out[[1]]$G)))
+ # values for Cp and c1 integrate to the same values of G
+ G.Cp <- subcrt(mod.obigt(list(name="test", formula="C0", G=0, S=0, Cp=100)))$out[[1]]$G
+ G.c1 <- subcrt(mod.obigt(list(name="test", formula="C0", G=0, S=0, c1=100)))$out[[1]]$G
+ expect_equal(G.Cp, G.c1)
})
Modified: pkg/CHNOSZ/man/eos.Rd
===================================================================
--- pkg/CHNOSZ/man/eos.Rd 2012-09-23 11:08:03 UTC (rev 20)
+++ pkg/CHNOSZ/man/eos.Rd 2012-09-23 15:11:26 UTC (rev 21)
@@ -44,12 +44,20 @@
The parameters in the \code{cgl} equations of state for crystalline, gas and liquid species (except liquid water) include \code{V}, \code{a}, \code{b}, \code{c}, \code{d}, \code{e}, \code{f} and \code{lambda}. The terms denoted by \code{a}, \code{b} and \code{c} correspond to the Maier-Kelley equation for heat capacity (Maier and Kelley, 1932); the additional terms are useful for representing heat capacities of minerals (Robie and Hemingway, 1995) and gaseous or liquid organic species (Helgeson et al., 1998). The standard molal volumes (\samp{V}) of species in these calculations are taken to be independent of temperature and pressure.
- The temperature and pressure range of validity of the revised HKF equations of state for aqueous species corresponds to the stability region of liquid water or the supercritical fluid at conditions between 0 to 1000 \eqn{^{\circ}}{degrees }C and 1 to 5000 bar (Tanger and Helgeson, 1988; Shock and Helgeson, 1988). The \code{hkf} function does not check these limits and will compute properties as long as the requisite electrostatic properties of water are available. There are conceptually no temperature limits (other than 0 Kelvin) for the validity of the \code{cgl} equations of state. However, the actual working upper temperature limits correspond to the temperatures of phase transitions of minerals or to those temperatures beyond which extrapolations from experimental data become untenable. These temperature limits are stored in the thermodynamic database, but \code{cgl} ignores them (\code{\link{subcrt}} warns if they are exceeded).
+For both \code{hkf} and \code{cgl}, if at least one equations-of-state parameter for a species is provided, any NA values of the other parameters are reset to zero.
+If all equations-of-state parameters are NA, but values of \samp{Cp} and/or \samp{V} are available, those values are used in the integration of \samp{G}, \samp{H} and \samp{S} as a function of temperature.
- The \code{T} and \code{P} arguments (and \code{rhohat}, \code{Tc}, \code{alpha}, \code{daldT}, \code{beta} for \code{gfun}) should all be the same length; the functions do no argument validating. An error produced by \code{gfun} for multiple values of \code{T} but only one \code{P} was \dQuote{NAs are not allowed in subscripted assignments}.
+The \code{T} and \code{P} arguments (and \code{rhohat}, \code{Tc}, \code{alpha}, \code{daldT}, \code{beta} for \code{gfun}) should all be the same length; the functions perform no argument validating.
}
+\section{Warning}{
+The temperature and pressure range of validity of the revised HKF equations of state for aqueous species corresponds to the stability region of liquid water or the supercritical fluid at conditions between 0 to 1000 \eqn{^{\circ}}{degrees }C and 1 to 5000 bar (Tanger and Helgeson, 1988; Shock and Helgeson, 1988).
+The \code{hkf} function does not check these limits and will compute properties as long as the requisite electrostatic properties of water are available. There are conceptually no temperature limits (other than 0 Kelvin) for the validity of the \code{cgl} equations of state.
+However, the actual working upper temperature limits correspond to the temperatures of phase transitions of minerals or to those temperatures beyond which extrapolations from experimental data become highly uncertain.
+These temperature limits are stored in the thermodynamic database for some minerals, but \code{cgl} ignores them; however, \code{\link{subcrt}} warns if they are exceeded.
+}
+
\value{
A list of length equal to the number of species (i.e., number rows of supplied \code{ghs} and \code{eos} values). Each element of the list contains a dataframe, each column of which corresponds to one of the specified properties; the number of rows is equal to the number of pressure-temperature points.
}
Modified: pkg/CHNOSZ/man/subcrt.Rd
===================================================================
--- pkg/CHNOSZ/man/subcrt.Rd 2012-09-23 11:08:03 UTC (rev 20)
+++ pkg/CHNOSZ/man/subcrt.Rd 2012-09-23 15:11:26 UTC (rev 21)
@@ -65,7 +65,11 @@
If \code{IS} is set to a single value other than zero, \code{\link{nonideal}} is used to calculate the apparent properties (\code{G}, \code{H}, \code{S} and \code{Cp}) of charged aqueous species at the given ionic strength. To perform calculations at a single \code{P} and \code{T} and for multiple values of ionic strength, supply these values in \code{IS}. Calculations can also be performed on a \code{P}-\code{IS}, \code{T}-\code{IS} or \code{P,T}-\code{IS} grid. Values of \code{logK} of reactions calculated for \code{IS} not equal to zero are consistent with the apparent Gibbs energies of the charged aqueous species.
- \code{subcrt} is modeled after the functionality of the \acronym{SUPCRT92} package (Johnson et al., 1992). Certain features of \acronym{SUPCRT92} are not available here, for example, calculations as a function of density of \eqn{\mathrm{H_2O}}{H2O} instead of pressure, or calculations of temperatures of univariant curves (i.e. for which \code{logK} is zero), or calculations of the molar volumes of quartz and coesite as a function of temperature (taken here to be constant). The informative messages produced by \code{SUPCRT92} when temperature or pressure limits of the equations of state are exceeded generally are not reproduced here. However, \code{NA}s may be produced in the output of \code{subcrt} if the requisite thermodynamic or electrostatic properties of water can not be calculated at given conditions.
+\code{subcrt} is modeled after the functionality of the \acronym{SUPCRT92} package (Johnson et al., 1992).
+Certain features of \acronym{SUPCRT92} are not available here, for example, calculations as a function of density of \eqn{\mathrm{H_2O}}{H2O} instead of pressure, or calculations of temperatures of univariant curves (i.e. for which \code{logK} is zero), or calculations of the molar volumes of quartz and coesite as a function of temperature (taken here to be constant).
+The informative messages produced by \code{SUPCRT92} when temperature or pressure limits of the equations of state are exceeded generally are not reproduced here.
+However, \code{NA}s may be produced in the output of \code{subcrt} if the requisite thermodynamic or electrostatic properties of water can not be calculated at given conditions.
+Specifically, \code{NA}s are produced for calculations at \samp{Psat} when the temperature exceeds the critical temperature of H2O.
\code{read.supcrt} and \code{write.supcrt} are used to read and write thermodynamic data \code{file}s in the format of \acronym{SUPCRT92}. \code{read.supcrt} parses the thermodynamic data contained in a \acronym{SUPCRT92}-format file into a dataframe that is compatible with the format used in CHNOSZ (see \code{\link{thermo}$obigt}). \code{write.supcrt} creates a \acronym{SUPCRT92} \code{file} using the dataframe given in the \code{obigt} argument (if missing the entire species database is processed). An edited version of the \samp{slop98.dat} data file (Shock et al., 1998) that is compatible with \code{read.supcrt} is available at \url{http://chnosz.net/download/} (small formatting quirks in the original version cause glitches with this function).
@@ -73,7 +77,10 @@
\section{Warning}{
- Comparison of the output of \code{subcrt} with that of \acronym{SUPCRT92} indicates the two give similar values of properties for neutral aqueous species up to 1000 \eqn{^{\circ}}{degrees }C and 5000 bar and for charged aqueous species over the temperature range 0 to 300 \eqn{^{\circ}}{degrees }C. At higher temperatures, there are significant divergences for charged species. For example, there is less than a 2 percent difference in the value of \eqn{{\Delta}G^{\circ}}{Delta G0} for K+(aq) at \eqn{^{\circ}}{degrees }C and 5000 bar, but this deviation increases with decreasing pressure at the same temperature. Further testing is necessary in CHNOSZ in connection with the \eqn{g}{g}- and \eqn{f}{f}-functions (Shock et al., 1992) for the partial derivatives of the omega parameter which are incorporated into the \code{\link{hkf}} function. (Note in particular the NaCl dissociation example under \code{\link{water}}, the results of which are noticeably different from the calculations of Shock et al., 1992.)
+Comparison of the output of \code{subcrt} with that of \acronym{SUPCRT92} indicates the two give similar values of properties for neutral aqueous species up to 1000 \eqn{^{\circ}}{degrees }C and 5000 bar.
+Changes were made in CHNOSZ 0.9-8 to improve the calculation of the \eqn{g}{g}- and \eqn{f}{f}-functions (Shock et al., 1992) for the partial derivatives of the omega parameter which are used by the \code{\link{hkf}} function, so thermodynamic properties of charged aqueous species over the temperature range 0 to 1000 \eqn{^{\circ}}{degrees }C are now mostly consistent with \acronym{SUPCRT92}.
+Note, however, that while \acronym{SUPCRT92} limits calculations at Psat to 350 \eqn{^{\circ}}{degrees }C (\dQuote{beyond range of applicability of aqueous species equations}), there is no corresponding limit in place in \code{subcrt} (or \code{\link{hkf}}), so that inapplicable calculations will be performed up to the critical temperature (373.917 \eqn{^{\circ}}{degrees }C) without any warning.
+It is probably for this reason that there is a noticeable jump in the value of logK at Psat shown in the one of the examples (\code{\link{longex}("NaCl")}).
A known bug is misidentification of the stable polymorph of some minerals at high temperature; an example of this bug is shown in the \code{\link{stopifnot}} at the end of the skarn example below.
Modified: pkg/CHNOSZ/man/util.data.Rd
===================================================================
--- pkg/CHNOSZ/man/util.data.Rd 2012-09-23 11:08:03 UTC (rev 20)
+++ pkg/CHNOSZ/man/util.data.Rd 2012-09-23 15:11:26 UTC (rev 21)
@@ -48,7 +48,14 @@
\code{add.obigt} reads data from the specified \code{file} and adds it to \code{\link{thermo}$obigt}. Set \code{force} to TRUE to replace species that exist in the thermodynamic database (each unique combination of a name and a state in the database is considered to be a unique species). \code{force}, if not specified, reverts to TRUE if the \code{file} is left at its default. Given the default setting of \code{E.units}, the function does not perform any unit conversions. If \code{E.units} is set to \samp{J}, then the thermodynamic parameters are converted from units of Joules to calories, as used in the CHNOSZ database.
- \code{mod.obigt} changes one or more of the properties of one or more species or adds species to the thermodynamic database. These changes are lost if you reload the database by calling \code{\link{data}(thermo)} or if you quit the \R session without saving it. The name of the species to add or change must be supplied as the first argument of \code{...} or as a named argument (named \samp{name}). When adding new species, a \code{formula} should be included \code{state} should also be included along with the values of any of the thermodynamic properties. Additional arguments refer to the name of the property(s) to be updated and are matched to any part of compound column names in \code{\link{thermo}$obigt}, such as \samp{z.T}. The values provided should be in the units specifed in the documentation for the \code{thermo} data object, including any order-of-magnitude scaling factors. When adding species, properties that are not specified become NA.
+\code{mod.obigt} changes one or more of the properties of species or adds species to the thermodynamic database.
+These changes are lost if you reload the database by calling \code{\link{data}(thermo)} or if you quit the \R session without saving it.
+The name of the species to add or change must be supplied as the first argument of \code{...} or as a named argument (named \samp{name}).
+When adding new species, a \samp{formula} should be included along with the values of any of the thermodynamic properties.
+Additional arguments refer to the name of the property(s) to be updated and are matched to any part of compound column names in \code{\link{thermo}$obigt}, such as \samp{z} or \samp{T} in \samp{z.T}.
+Unless \samp{state} is specified as one of the properties, its value is taken from \code{thermo$opt$state}.
+When adding species, properties that are not specified become NA (except for \samp{state}).
+The values provided should be in the units specifed in the documentation for the \code{thermo} data object, including any order-of-magnitude scaling factors.
\code{today} returns the current date in the format adopted for \code{thermo$obigt} (inherited from SUPCRT-format data files) e.g. \samp{13.May.12} for May 13, 2012.
More information about the CHNOSZ-commits
mailing list