[CHNOSZ-commits] r88 - in pkg/CHNOSZ: . R data demo inst inst/tests man tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Jun 9 12:21:09 CEST 2015
Author: jedick
Date: 2015-06-09 12:21:09 +0200 (Tue, 09 Jun 2015)
New Revision: 88
Added:
pkg/CHNOSZ/R/util.water.R
pkg/CHNOSZ/demo/density.R
pkg/CHNOSZ/inst/tests/test-util.program.R
pkg/CHNOSZ/man/util.water.Rd
Removed:
pkg/CHNOSZ/demo/cordierite.R
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/R/IAPWS95.R
pkg/CHNOSZ/R/affinity.R
pkg/CHNOSZ/R/equilibrate.R
pkg/CHNOSZ/R/examples.R
pkg/CHNOSZ/R/objective.R
pkg/CHNOSZ/R/util.fasta.R
pkg/CHNOSZ/R/util.program.R
pkg/CHNOSZ/R/water.R
pkg/CHNOSZ/data/opt.csv
pkg/CHNOSZ/demo/00Index
pkg/CHNOSZ/inst/NEWS
pkg/CHNOSZ/man/CHNOSZ-package.Rd
pkg/CHNOSZ/man/IAPWS95.Rd
pkg/CHNOSZ/man/data.Rd
pkg/CHNOSZ/man/examples.Rd
pkg/CHNOSZ/man/util.fasta.Rd
pkg/CHNOSZ/man/util.program.Rd
pkg/CHNOSZ/man/water.Rd
pkg/CHNOSZ/tests/test-all.R
Log:
improve density calculations for IAPWS-95
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2015-05-19 14:38:46 UTC (rev 87)
+++ pkg/CHNOSZ/DESCRIPTION 2015-06-09 10:21:09 UTC (rev 88)
@@ -1,11 +1,11 @@
-Date: 2015-05-19
+Date: 2015-06-09
Package: CHNOSZ
-Version: 1.0.5
+Version: 1.0.5-1
Title: Chemical Thermodynamics and Activity Diagrams
Author: Jeffrey Dick
Maintainer: Jeffrey Dick <j3ffdick at gmail.com>
-Depends: R (>= 3.0.0)
-Suggests: limSolve, parallel, testthat, knitr
+Depends: R (>= 3.1.0)
+Suggests: limSolve, testthat, knitr
Description: Functions and data sets to support chemical thermodynamic modeling in biochemistry
and low-temperature geochemistry. The features include calculation of the standard molal
thermodynamic properties and chemical affinities of reactions involving minerals and/or
Modified: pkg/CHNOSZ/R/IAPWS95.R
===================================================================
--- pkg/CHNOSZ/R/IAPWS95.R 2015-05-19 14:38:46 UTC (rev 87)
+++ pkg/CHNOSZ/R/IAPWS95.R 2015-06-09 10:21:09 UTC (rev 88)
@@ -246,54 +246,3 @@
return(ww)
}
-WP02.auxiliary <- function(property='rho.liquid',T=298.15) {
- # auxiliary equations for liquid-vapor phase boundary
- # from Wagner and Pruss, 2002
- # critical point
- T.critical <- 647.096 # K
- P.critical <- 22.064 # MPa
- rho.critical <- 322 # kg m-3
-
- if(property %in% c("P.sigma","dP.sigma.dT")) {
- # vapor pressure
- V <- 1 - T / T.critical # theta (dimensionless)
- a1 <- -7.85951783
- a2 <- 1.84408259
- a3 <- -11.7866497
- a4 <- 22.6807411
- a5 <- -15.9618719
- a6 <- 1.80122502
- ln.P.sigma.P.critical <- T.critical / T *
- ( a1*V + a2*V^1.5 + a3*V^3 + a4*V^3.5 + a5*V^4 + a6*V^7.5 )
- P.sigma <- P.critical * exp(ln.P.sigma.P.critical)
- if(property=="dP.sigma.dT") out <- - P.sigma / T * ( ln.P.sigma.P.critical +
- a1 + 1.5*a2*V^0.5 + 3*a3*V^2 + 3.5*a4*V^2.5 + 4*a5*V^3 + 7.5*a6*V^6.5 )
- else out <- P.sigma
- } else if(property=="rho.liquid") {
- # saturated liquid density
- V <- 1 - T / T.critical
- b1 <- 1.99274064
- b2 <- 1.09965342
- b3 <- -0.510839303
- b4 <- -1.75493479
- b5 <- -45.5170352
- b6 <- -6.74694450E5
- rho.liquid <- rho.critical * (
- 1 + b1*V^(1/3) + b2*V^(2/3) + b3*V^(5/3) + b4*V^(16/3) + b5*V^(43/3) + b6*V^(110/3) )
- out <- rho.liquid
- } else if(property=="rho.vapor") {
- # saturated vapor density
- V <- 1 - T / T.critical
- c1 <- -2.03150240
- c2 <- -2.68302940
- c3 <- -5.38626492
- c4 <- -17.2991605
- c5 <- -44.7586581
- c6 <- -63.9201063
- rho.vapor <- rho.critical * exp (
- c1*V^(2/6) + c2*V^(4/6) + c3*V^(8/6) + c4*V^(18/6) + c5*V^(37/6) + c6*V^(71/6) )
- out <- rho.vapor
- } else stop(paste('i can not calculate',property))
- return(out)
-}
-
Modified: pkg/CHNOSZ/R/affinity.R
===================================================================
--- pkg/CHNOSZ/R/affinity.R 2015-05-19 14:38:46 UTC (rev 87)
+++ pkg/CHNOSZ/R/affinity.R 2015-06-09 10:21:09 UTC (rev 88)
@@ -152,15 +152,13 @@
# from those of residues
loga.protein <- rep(loga.protein,length.out=length(iprotein))
protein.fun <- function(ip) {
- if(ip %% 50 == 0) msgout(paste(ip,"..",sep=""))
tpext <- as.numeric(thermo$protein[iprotein[ip],5:25])
- return(Reduce("+", pprod(a[ires],tpext)) - loga.protein[ip])
+ return(Reduce("+", CHNOSZ::pprod(a[ires],tpext)) - loga.protein[ip])
}
# use another level of indexing to let the function
# report on its progress
jprotein <- 1:length(iprotein)
- protein.affinity <- palply(jprotein,protein.fun)
- if(length(iprotein) > 49) msgout("\n")
+ protein.affinity <- palply("", jprotein, protein.fun)
## update the species list
# we use negative values for ispecies to denote that
# they index thermo$protein and not thermo$species
Modified: pkg/CHNOSZ/R/equilibrate.R
===================================================================
--- pkg/CHNOSZ/R/equilibrate.R 2015-05-19 14:38:46 UTC (rev 87)
+++ pkg/CHNOSZ/R/equilibrate.R 2015-06-09 10:21:09 UTC (rev 88)
@@ -89,18 +89,18 @@
Astardim <- dim(Astar[[1]])
Anames <- names(Astar)
# first loop: make vectors
- A <- palply(1:length(A), function(i) as.vector(A[[i]]))
+ A <- palply("", 1:length(A), function(i) as.vector(A[[i]]))
# second loop: get the exponentiated Astars (numerators)
# need to convert /2.303RT to /RT
#A[[i]] <- exp(log(10)*Astar[[i]]/n.balance[i])/n.balance[i]
- A <- palply(1:length(A), function(i) exp(log(10)*Astar[[i]]/n.balance[i]))
+ A <- palply("", 1:length(A), function(i) exp(log(10)*Astar[[i]]/n.balance[i]))
# third loop: accumulate the denominator
# initialize variable to hold the sum
At <- A[[1]]
At[] <- 0
for(i in 1:length(A)) At <- At + A[[i]]*n.balance[i]
# fourth loop: calculate log abundances and replace the dimensions
- A <- palply(1:length(A), function(i) loga.balance + log10(A[[i]]/At))
+ A <- palply("", 1:length(A), function(i) loga.balance + log10(A[[i]]/At))
# fifth loop: replace dimensions
for(i in 1:length(A)) dim(A[[i]]) <- Astardim
# add names and we're done!
@@ -206,7 +206,7 @@
return(Abar)
}
# calculate the logact(thing) for each condition
- logact <- palply(1:nrow(Astar), function(i) {
+ logact <- palply("", 1:nrow(Astar), function(i) {
# get the equilibrium Abar for each condition
Abar <- Abarfun(i)
return(logactfun(Abar, i))
Modified: pkg/CHNOSZ/R/examples.R
===================================================================
--- pkg/CHNOSZ/R/examples.R 2015-05-19 14:38:46 UTC (rev 87)
+++ pkg/CHNOSZ/R/examples.R 2015-06-09 10:21:09 UTC (rev 88)
@@ -15,10 +15,10 @@
"objective", "revisit", "transfer", "anim", "EOSregress", "wjd")
plot.it <- FALSE
if(is.character(do.png))
- png(paste(do.png,"%d.png",sep=""),width=700,height=700,pointsize=18)
+ png(paste(do.png,"%d.png",sep=""),width=500,height=500,pointsize=12)
else if(do.png) plot.it <- TRUE
for(i in 1:length(topics)) {
- if(plot.it) png(paste(topics[i],"%d.png",sep=""),width=700,height=700,pointsize=18)
+ if(plot.it) png(paste(topics[i],"%d.png",sep=""),width=500,height=500,pointsize=12)
myargs <- list(topic=topics[i],ask=FALSE)
do.call(example,myargs)
if(plot.it) dev.off()
@@ -29,16 +29,18 @@
cat("Time elapsed: ", proc.time() - .ptime, "\n")
}
-demos <- function(which=c("sources", "NaCl", "cordierite",
+demos <- function(which=c("sources", "NaCl", "density",
"phosphate", "nucleobase", "orp", "diagram", "revisit", "findit",
"CO2Ac", "nonideal", "ionize", "buffer", "yeastgfp", "mosaic",
- "solubility", "wjd")) {
+ "solubility", "wjd"), do.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)) {
# say something so the user sees where we are
msgout("------------\n")
msgout(paste("demos: running '", which[i], "'\n", sep=""))
+ if(do.png) png(paste(which[i],"%d.png",sep=""),width=500,height=500,pointsize=12)
out <- demo(which[i], package="CHNOSZ", character.only=TRUE, echo=FALSE, ask=FALSE)
+ if(do.png) dev.off()
}
return(invisible(out))
}
Modified: pkg/CHNOSZ/R/objective.R
===================================================================
--- pkg/CHNOSZ/R/objective.R 2015-05-19 14:38:46 UTC (rev 87)
+++ pkg/CHNOSZ/R/objective.R 2015-06-09 10:21:09 UTC (rev 88)
@@ -186,7 +186,7 @@
return(sum(DG))
}
# we need to index both loga1 and Astar
- DGtr <- unlist(palply(seq(nrow(loga1)), function(i) {
+ DGtr <- unlist(lapply(seq(nrow(loga1)), function(i) {
dgtr(loga1[i, ], loga2, Astar[i, ])
}))
return(DGtr)
Modified: pkg/CHNOSZ/R/util.fasta.R
===================================================================
--- pkg/CHNOSZ/R/util.fasta.R 2015-05-19 14:38:46 UTC (rev 87)
+++ pkg/CHNOSZ/R/util.fasta.R 2015-06-09 10:21:09 UTC (rev 88)
@@ -33,11 +33,11 @@
# use the system grep
if(is.null(startswith)) startswith <- "" else startswith <- paste("^",startswith,".*",sep="")
if(ignore.case) ic <- "-i" else ic <- ""
- out <- palply(1:length(pattern),sysgrep)
+ out <- lapply(1:length(pattern), sysgrep)
} else {
# use R grep
if(is.null(lines)) lines <- readLines(file)
- out <- palply(1:length(pattern),Rgrep)
+ out <- lapply(1:length(pattern), Rgrep)
}
# make numeric (NA for ones that aren't matched)
out <- as.numeric(sapply(out,as.numeric))
@@ -108,7 +108,7 @@
organism <- bnf
# protein/gene name is from header line for entry
# (strip the ">" and go to the first space)
- if(is.null(id)) id <- as.character(palply(1:length(i), function(j) {
+ if(is.null(id)) id <- as.character(palply("", 1:length(i), function(j) {
# get the text of the line
f1 <- linefun(i[j],i[j])
# stop if the first character is not ">"
@@ -206,7 +206,7 @@
return(count)
}
# counts for each sequence
- a <- palply(seq, countfun, start, stop)
+ a <- palply("", seq, countfun, start, stop)
a <- t(as.data.frame(a, optional=TRUE))
# clean up row/column names
colnames(a) <- letts
Modified: pkg/CHNOSZ/R/util.program.R
===================================================================
--- pkg/CHNOSZ/R/util.program.R 2015-05-19 14:38:46 UTC (rev 87)
+++ pkg/CHNOSZ/R/util.program.R 2015-06-09 10:21:09 UTC (rev 88)
@@ -17,17 +17,20 @@
return(name)
}
-palply <- function(X, FUN, ...) {
- # a wrapper function to run parLapply if length(X) > 10000
+palply <- function(varlist, X, FUN, ...) {
+ # a wrapper function to run parLapply if length(X) >= thermo$opt$paramin
# and package 'parallel' is available, otherwise run lapply
- if(length(X) > 10000 & "parallel" %in% (.packages())) {
- # the calculations - modified from ?parLapply
- ## Use option mc.cores to choose an appropriate cluster size.
- # or detectCores if that is NULL, and set max at 2 for now
- # (to be nice to CRAN etc.)
- nCores <- max(getOption("mc.cores", parallel::detectCores()), 2)
+ if(length(X) >= get("thermo")$opt$paramin) {
+ # Use option mc.cores to choose an appropriate cluster size.
+ # and set max at 2 for now (per CRAN policies)
+ nCores <- min(getOption("mc.cores"), 2)
# don't load methods package, to save startup time - ?makeCluster
cl <- parallel::makeCluster(nCores, methods=FALSE)
+ # export the variables and notify the user
+ if(! "" %in% varlist) parallel::clusterExport(cl, varlist)
+ msgout(paste("palply:", caller.name(4), "running", length(X), "calculations on",
+ nCores, "cores with variable(s)", paste(varlist, collapse=", "), "\n"))
+ # run the calculations
out <- parallel::parLapply(cl, X, FUN, ...)
parallel::stopCluster(cl)
} else out <- lapply(X, FUN, ...)
Added: pkg/CHNOSZ/R/util.water.R
===================================================================
--- pkg/CHNOSZ/R/util.water.R (rev 0)
+++ pkg/CHNOSZ/R/util.water.R 2015-06-09 10:21:09 UTC (rev 88)
@@ -0,0 +1,181 @@
+# CHNOSZ/util.water.R
+
+WP02.auxiliary <- function(property='rho.liquid',T=298.15) {
+ # auxiliary equations for liquid-vapor phase boundary
+ # from Wagner and Pruss, 2002
+ # critical point
+ T.critical <- 647.096 # K
+ P.critical <- 22.064 # MPa
+ rho.critical <- 322 # kg m-3
+
+ if(property %in% c("P.sigma","dP.sigma.dT")) {
+ # vapor pressure
+ V <- 1 - T / T.critical # theta (dimensionless)
+ a1 <- -7.85951783
+ a2 <- 1.84408259
+ a3 <- -11.7866497
+ a4 <- 22.6807411
+ a5 <- -15.9618719
+ a6 <- 1.80122502
+ ln.P.sigma.P.critical <- T.critical / T *
+ ( a1*V + a2*V^1.5 + a3*V^3 + a4*V^3.5 + a5*V^4 + a6*V^7.5 )
+ P.sigma <- P.critical * exp(ln.P.sigma.P.critical)
+ if(property=="dP.sigma.dT") out <- - P.sigma / T * ( ln.P.sigma.P.critical +
+ a1 + 1.5*a2*V^0.5 + 3*a3*V^2 + 3.5*a4*V^2.5 + 4*a5*V^3 + 7.5*a6*V^6.5 )
+ else out <- P.sigma
+ } else if(property=="rho.liquid") {
+ # saturated liquid density
+ V <- 1 - T / T.critical
+ b1 <- 1.99274064
+ b2 <- 1.09965342
+ b3 <- -0.510839303
+ b4 <- -1.75493479
+ b5 <- -45.5170352
+ b6 <- -6.74694450E5
+ rho.liquid <- rho.critical * (
+ 1 + b1*V^(1/3) + b2*V^(2/3) + b3*V^(5/3) + b4*V^(16/3) + b5*V^(43/3) + b6*V^(110/3) )
+ out <- rho.liquid
+ } else if(property=="rho.vapor") {
+ # saturated vapor density
+ V <- 1 - T / T.critical
+ c1 <- -2.03150240
+ c2 <- -2.68302940
+ c3 <- -5.38626492
+ c4 <- -17.2991605
+ c5 <- -44.7586581
+ c6 <- -63.9201063
+ rho.vapor <- rho.critical * exp (
+ c1*V^(2/6) + c2*V^(4/6) + c3*V^(8/6) + c4*V^(18/6) + c5*V^(37/6) + c6*V^(71/6) )
+ out <- rho.vapor
+ } else stop(paste('i can not calculate',property))
+ return(out)
+}
+
+# return a density in kg m-3
+# corresponding to the given pressure (bar) and temperature (K)
+rho.IAPWS95 <- function(T=298.15, P=1, state="", trace=0) {
+ # function for which to find a zero
+ dP <- function(rho, T, P.MPa) IAPWS95("P", rho=rho, T=T)[, 1] - P.MPa
+ # convert bar to MPa
+ P.MPa <- convert(P, "MPa")
+ rho <- numeric()
+ T.critical <- 647.096 # K
+ P.critical <- 22.064 # MPa
+ for(i in 1:length(T)) {
+ Psat <- WP02.auxiliary("P.sigma", T[i])
+ if(T[i] > T.critical) {
+ # above critical temperature
+ interval <- c(0.1, 1)
+ extendInt <- "upX"
+ if(trace > 0) msgout("supercritical (T) ")
+ } else if(P.MPa[i] > P.critical) {
+ # above critical pressure
+ rho.sat <- WP02.auxiliary("rho.liquid", T=T[i])
+ interval <- c(rho.sat, rho.sat + 1)
+ extendInt <- "upX"
+ if(trace > 0) msgout("supercritical (P) ")
+ } else if(P.MPa[i] <= 0.9999*Psat) {
+ # steam
+ rho.sat <- WP02.auxiliary("rho.vapor", T=T[i])
+ interval <- c(rho.sat*0.1, rho.sat)
+ extendInt <- "upX"
+ if(trace > 0) msgout("steam ")
+ } else if(P.MPa[i] >= 1.00005*Psat) {
+ # water
+ rho.sat <- WP02.auxiliary("rho.liquid", T=T[i])
+ interval <- c(rho.sat, rho.sat + 1)
+ extendInt <- "upX"
+ if(trace > 0) msgout("water ")
+ } else if(!state %in% c("liquid", "vapor")) {
+ # we're close to the saturation curve;
+ # calculate rho and G for liquid and vapor and return rho for stable phase
+ if(trace > 0) msgout("close to saturation; trying liquid and vapor\n")
+ rho.liquid <- rho.IAPWS95(T[i], P[i], state="liquid", trace=trace)
+ rho.vapor <- rho.IAPWS95(T[i], P[i], state="vapor", trace=trace)
+ G.liquid <- IAPWS95("G", rho=rho.liquid, T=T[i])
+ G.vapor <- IAPWS95("G", rho=rho.vapor, T=T[i])
+ if(G.liquid < G.vapor) {
+ this.rho <- rho.liquid
+ if(trace > 0) msgout(paste0("G.liquid(", G.liquid, ") < G.vapor(", G.vapor, ")\n"))
+ } else {
+ this.rho <- rho.vapor
+ if(trace > 0) msgout(paste0("G.vapor(", G.vapor, ") < G.liquid (", G.liquid, ")\n"))
+ }
+ rho <- c(rho, this.rho)
+ next
+ } else {
+ # we are looking at a specific state
+ if(trace > 0) msgout(paste("specified state:", state, " "))
+ if(state=="vapor") rho0 <- WP02.auxiliary("rho.vapor", T[i])
+ else if(state=="liquid") rho0 <- WP02.auxiliary("rho.liquid", T[i])
+ # a too-big range may cause problems e.g.
+ # interval <- c(rho0*0.9, rho0*1.1) fails for T=253.15, P=1
+ interval <- c(rho0*0.95, rho0*1.05)
+ # if P on the initial interval are both higher or lower than target P,
+ # set the direction of interval extension
+ P.init <- IAPWS95("P", rho=interval, T=c(T[i], T[i]))[, 1]
+ if(all(P.init < P.MPa[i])) extendInt <- "downX"
+ else if(all(P.init > P.MPa[i])) extendInt <- "upX"
+ else extendInt <- "yes"
+ }
+ if(trace > 0) msgout(paste0("T=", T[i], " P=", P[i], " rho=[", interval[1], ",", interval[2], "]\n"))
+ this.rho <- try(uniroot(dP, interval, extendInt=extendInt, trace=trace, T=T[i], P.MPa=P.MPa[i])$root, TRUE)
+ if(!is.numeric(this.rho)) {
+ warning("rho.IAPWS95: problems finding density at ", T[i], " K and ", P[i], " bar", call.=FALSE)
+ this.rho <- NA
+ }
+ rho <- c(rho, this.rho)
+ }
+ return(rho)
+}
+
+water.AW90 <- function(T=298.15,rho=1000,P=0.1) {
+ # Equations for the dielectric constant of water
+ # from Archer and Wang, 1990
+ # T in K
+ # rho in kg m-3
+ # p in MPa
+
+ # Table 2
+ b <- c(-4.044525E-2, 103.6180 , 75.32165 ,
+ -23.23778 ,-3.548184 ,-1246.311 ,
+ 263307.7 ,-6.928953E-1,-204.4473)
+ alpha <- 18.1458392E-30 # m^3
+ #alpha <- 14.7E-30
+ mu <- 6.1375776E-30 # C m
+ N.A <- 6.0221367E23 # mol-1
+ k <- 1.380658E-23 # Boltzmann constant, J K-1
+ M <- 0.0180153 # kg mol-1
+ rho.0 <- 1000 # kg m-3
+ # Equation 1
+ epsilon.0 <- 8.8541878E-12 # permittivity of vacuum, C^2 J-1 m-1
+ #epsfun.lhs <- function(e) (e-1)*(2*e+1)/(9*e)
+ epsfun.rhs <- function(T,V.m) N.A*(alpha+mufun()/(3*epsilon.0*k*T))/(3*V.m)
+ #epsfun <- function(e,T,V.m) epsfun.lhs(e) - epsfun.rhs(T,V.m)
+ mufun <- function() gfun()*mu^2
+ gfun <- function() rhofun()*rho/rho.0 + 1
+ # Equation 3
+ rhofun <- function() b[1]*P*T^-1 + b[2]*T^-0.5 + b[3]*(T-215)^-1 +
+ b[4]*(T-215)^-0.5 + b[5]*(T-215)^-0.25 +
+ exp(b[6]*T^-1 + b[7]*T^-2 + b[8]*P*T^-1 + b[9]*P*T^-2)
+ epsilon <- function(T,rho) {
+ #tu <- try(uniroot(epsfun,c(1E-1,1E3),T=T,V.m=M/rho)$root,TRUE)
+ epspoly <- function() epsfun.rhs(T,V.m=M/rho)
+ tu <- (9*epspoly() + 1 + ((9*epspoly()+1)*(9*epspoly()+1) + 8)^0.5) / 4 #Marc Neveu added 4/24/2013
+ if(!is.numeric(tu)) {
+ warning('water.AW90: no root for density at ',T,' K and ',rho,' kg m-3.',call.=FALSE,immediate.=TRUE)
+ tu <- NA
+ }
+ return(tu)
+ }
+ # get things the right length
+ our.T <- T; our.rho <- rho; our.P <- P
+ t <- numeric()
+ for(i in 1:length(our.T)) {
+ T <- our.T[i]
+ rho <- our.rho[i]
+ P <- our.P[i]
+ t <- c(t,epsilon(T,rho))
+ }
+ return(t)
+}
Modified: pkg/CHNOSZ/R/water.R
===================================================================
--- pkg/CHNOSZ/R/water.R 2015-05-19 14:38:46 UTC (rev 87)
+++ pkg/CHNOSZ/R/water.R 2015-06-09 10:21:09 UTC (rev 88)
@@ -2,57 +2,6 @@
# calculate thermodynamic and electrostatic properties of H2O
# 20061016 jmd
-water.AW90 <- function(T=298.15,rho=1000,P=0.1) {
- # Equations for the dielectric constant of water
- # from Archer and Wang, 1990
- # T in K
- # rho in kg m-3
- # p in MPa
-
- # Table 2
- b <- c(-4.044525E-2, 103.6180 , 75.32165 ,
- -23.23778 ,-3.548184 ,-1246.311 ,
- 263307.7 ,-6.928953E-1,-204.4473)
- alpha <- 18.1458392E-30 # m^3
- #alpha <- 14.7E-30
- mu <- 6.1375776E-30 # C m
- N.A <- 6.0221367E23 # mol-1
- k <- 1.380658E-23 # Boltzmann constant, J K-1
- M <- 0.0180153 # kg mol-1
- rho.0 <- 1000 # kg m-3
- # Equation 1
- epsilon.0 <- 8.8541878E-12 # permittivity of vacuum, C^2 J-1 m-1
- #epsfun.lhs <- function(e) (e-1)*(2*e+1)/(9*e)
- epsfun.rhs <- function(T,V.m) N.A*(alpha+mufun()/(3*epsilon.0*k*T))/(3*V.m)
- #epsfun <- function(e,T,V.m) epsfun.lhs(e) - epsfun.rhs(T,V.m)
- mufun <- function() gfun()*mu^2
- gfun <- function() rhofun()*rho/rho.0 + 1
- # Equation 3
- rhofun <- function() b[1]*P*T^-1 + b[2]*T^-0.5 + b[3]*(T-215)^-1 +
- b[4]*(T-215)^-0.5 + b[5]*(T-215)^-0.25 +
- exp(b[6]*T^-1 + b[7]*T^-2 + b[8]*P*T^-1 + b[9]*P*T^-2)
- epsilon <- function(T,rho) {
- #tu <- try(uniroot(epsfun,c(1E-1,1E3),T=T,V.m=M/rho)$root,TRUE)
- epspoly <- function() epsfun.rhs(T,V.m=M/rho)
- tu <- (9*epspoly() + 1 + ((9*epspoly()+1)*(9*epspoly()+1) + 8)^0.5) / 4 #Marc Neveu added 4/24/2013
- if(!is.numeric(tu)) {
- warning('water.AW90: no root for density at ',T,' K and ',rho,' kg m-3.',call.=FALSE,immediate.=TRUE)
- tu <- NA
- }
- return(tu)
- }
- # get things the right length
- our.T <- T; our.rho <- rho; our.P <- P
- t <- numeric()
- for(i in 1:length(our.T)) {
- T <- our.T[i]
- rho <- our.rho[i]
- P <- our.P[i]
- t <- c(t,epsilon(T,rho))
- }
- return(t)
-}
-
water <- function(property = NULL, T = get("thermo")$opt$Tr, P = "Psat") {
# calculate the properties of liquid H2O as a function of T and P
# T in Kelvin, P in bar
@@ -192,43 +141,6 @@
return(w.out[, iprop, drop=FALSE])
}
-rho.IAPWS95 <- function(T=298.15, P=1) {
- # return a density in kg m-3
- # corresponding to the given pressure (bar) and temperature (K)
- if(identical(P, "Psat")) stop("this function doesn't take P='Psat'")
- dP <- function(rho, T, P) {
- dP <- IAPWS95("P", rho=rho, T=T)[, 1] - convert(P, "MPa")
- return(dP)
- }
- rho <- numeric()
- for(i in 1:length(T)) {
- if(T[i] < 647.096) {
- rho.lower <- WP02.auxiliary('rho.liquid',T=T[i])-2
- rho.upper <- rho.lower + 400
- if(P[i] < 5000) rho.upper <- rho.lower + 300
- if(P[i] < 1000) rho.upper <- rho.lower + 200
- if(P[i] < 300) {
- rho.upper <- rho.lower + 30
- if(T[i] < 250) { #Marc Neveu added 4/23/2013
- rho.lower <- rho.lower - 10
- rho.upper <- rho.lower + 40
- }
- }
- } else {
- rho.lower <- 0.01
- rho.upper <- 1200
- }
- this.rho <- try(uniroot(dP, c(rho.lower, rho.upper), T=T[i], P=P[i])$root, TRUE)
- if(!is.numeric(this.rho)) {
- warning("rho.IAPWS95: no root for density between ", round(rho.lower, 1),
- " and ", round(rho.upper, 1), " kg m-3 at ", T[i], " K and ", P[i], " bar", call.=FALSE)
- this.rho <- NA
- }
- rho <- c(rho, this.rho)
- }
- return(rho)
-}
-
water.IAPWS95 <- function(property, T=298.15, P=1) {
# to get the properties of water via IAPWS-95
msgout(paste("water.IAPWS95: calculating", length(T), "values for"))
@@ -373,7 +285,7 @@
# calculate values of P for Psat
if(identical(P, "Psat")) P <- psat()
msgout(" rho")
- my.rho <- rho.IAPWS95(T, P)
+ my.rho <- rho.IAPWS95(T, P, get("thermo")$opt$IAPWS95.Psat.state)
rho <- function() my.rho
}
for(i in 1:length(property)) {
Modified: pkg/CHNOSZ/data/opt.csv
===================================================================
--- pkg/CHNOSZ/data/opt.csv 2015-05-19 14:38:46 UTC (rev 87)
+++ pkg/CHNOSZ/data/opt.csv 2015-06-09 10:21:09 UTC (rev 88)
@@ -1,2 +1,2 @@
-Tr,Pr,Theta,Psi,R,cutoff,E.units,T.units,P.units,state,water,G.tol,Cp.tol,V.tol,varP
-298.15,1,228,2600,1.9872,1e-10,cal,C,bar,aq,SUPCRT92,100,1,1,FALSE
+Tr,Pr,Theta,Psi,R,cutoff,E.units,T.units,P.units,state,water,G.tol,Cp.tol,V.tol,varP,IAPWS95.Psat.state,paramin
+298.15,1,228,2600,1.9872,1e-10,cal,C,bar,aq,SUPCRT92,100,1,1,FALSE,liquid,1000
Modified: pkg/CHNOSZ/demo/00Index
===================================================================
--- pkg/CHNOSZ/demo/00Index 2015-05-19 14:38:46 UTC (rev 87)
+++ pkg/CHNOSZ/demo/00Index 2015-06-09 10:21:09 UTC (rev 88)
@@ -1,6 +1,6 @@
sources cross-check the reference list with the thermodynamic database
NaCl equilibrium constant for aqueous NaCl dissociation
-cordierite equilibrium constant of hydrous cordierite dehydration
+density density of H2O, inverted from IAPWS-95 equations
phosphate phosphate speciation with pH, temperature and ionic strength
nucleobase relative stabilities of nucleobases and some amino acids
orp oxidation-reduction potential of redox standards as a function of temperature
Deleted: pkg/CHNOSZ/demo/cordierite.R
===================================================================
--- pkg/CHNOSZ/demo/cordierite.R 2015-05-19 14:38:46 UTC (rev 87)
+++ pkg/CHNOSZ/demo/cordierite.R 2015-06-09 10:21:09 UTC (rev 88)
@@ -1,19 +0,0 @@
-### 1-D property plot
-## for hydrous cordierite = cordierite + H2O
-## after Helgeson et al., 1978
-## (Summary and critique of the thermodynamic properties of
-## rock-forming minerals. Am. J. Sci., 278-A, 1-229)
-basis(c("cordierite,hydrous","Mg+2","SiO2","H2O","O2","H+"))
-species("cordierite")
-# water.SUPCRT92 can only get us up to 5000 bar
-# (lines for 7000 and 10000 bar are in the original diagram)
-P <- c(1,2,3,5)*1000
-col <- rainbow(length(P))
-for(i in 1:length(P)) {
- a <- affinity(property="logK",T=c(20,800),P=P[i])
- diagram(a,add=(i!=1),ylim=c(-4,2),legend.x=NULL,
- col=col[i],main="")
-}
-legend("topright",lty=1,col=col,legend=paste(P,"bar"))
-title(main=paste("hydrous cordierite = cordierite + H2O",
- "After Helgeson et al., 1978",sep="\n"),cex.main=0.9)
Added: pkg/CHNOSZ/demo/density.R
===================================================================
--- pkg/CHNOSZ/demo/density.R (rev 0)
+++ pkg/CHNOSZ/demo/density.R 2015-06-09 10:21:09 UTC (rev 88)
@@ -0,0 +1,90 @@
+# make T-P diagram for H2O, colored according to density
+
+# IAPWS95 or SUPCRT92
+method <- "IAPWS95"
+# low or high T,P range
+TPrange <- "low"
+
+blue <- "blue"
+if(TPrange=="low") {
+ T <- seq(300, 700, 5)
+ P <- seq(1, 351, 5)
+ bias <- 1.68
+} else {
+ # upper T,P limit for SUPCRT92: 2250 degC, 30000 bar
+ T <- seq(273, 2523, 50)
+ P <- seq(1, 30000, length.out=50)
+ # to attempt to match the colors using the different methods
+ # (ranges are different because IAPWS95 reports higher density in
+ # the high-P, low-T region, where SUPCRT92 doesn't give output)
+ if(method=="IAPWS95") bias <- 2.2
+ else if(method=="SUPCRT92") {
+ bias <- 2.1
+ blue <- "#0d0dff"
+ }
+}
+TP <- expand.grid(T=T, P=P)
+if(method=="IAPWS95") {
+ # the following should trigger parallel calculations
+ # if nrow(TP) (5751 for TPrange="low") is >= thermo$opt$paramin (default 1000)
+ rho <- palply("TP", 1:nrow(TP), function(i){CHNOSZ::rho.IAPWS95(TP$T[i], TP$P[i])})
+} else if(method=="SUPCRT92") {
+ rho <- water.SUPCRT92("rho", TP$T, TP$P)
+ # water.SUPCRT92 returns 0 when the density can't be calculated
+ rho[rho==0] <- NA
+}
+rho.num <- unlist(rho)
+rho.mat <- matrix(rho.num, nrow=length(T), ncol=length(P))
+# blueest for most dense, reddest for least dense
+# bias is adjusted to white for the critical density
+ncol <- 100
+col <- colorRampPalette(c("red", "white", blue), bias=bias)(ncol)
+# first make a background image (for debugging -
+# will be visible only if some density calculations fail)
+fill.mat <- matrix(0, nrow=length(T), ncol=length(P))
+image(T, P, fill.mat, col="black", xlab=axis.label("T", "K"), ylab=axis.label("P"))
+# now plot densities
+image(T, P, rho.mat, col=col, add=TRUE)
+# add a title and calculate saturation line
+if(method=="IAPWS95") {
+ title(main=expression("Density of"~H[2]*O~"inverted from IAPWS-95 equations"))
+ title(main=expression("Line calculated using auxiliary equations for saturation"), line=0.8)
+ Psat <- convert(WP02.auxiliary("P.sigma", T), "bar")
+} else if(method=="SUPCRT92") {
+ title(main=expression("Density of"~H[2]*O~"calculated using SUPCRT92"))
+ Psat <- water.SUPCRT92("Psat", T, "Psat")[,1]
+}
+# plot saturation line
+lines(T, Psat, lwd=6)
+lines(T, Psat, lwd=3, col="gold")
+# add a color key
+if(TPrange=="low") {
+ x <- c(333, 366, 370)
+ y <- c(100, 300)
+} else if(TPrange=="high") {
+ x <- c(600, 780, 800)
+ y <- c(10000, 25000)
+}
+ykey <- seq(y[1], y[2], length.out=ncol+1)
+for(i in 1:ncol) rect(x[1], ykey[i], x[2], ykey[i+1], col=col[i], border=NA)
+rect(x[1], ykey[1], x[2], rev(ykey)[1])
+# label the extrema
+rrange <- range(rho.num, na.rm=TRUE)
+text(x[3], ykey[1], as.expression(substitute(x~kg/m^3, list(x=round(rrange[1], 4)))), adj=0)
+text(x[3], rev(ykey)[1], as.expression(substitute(x~kg/m^3, list(x=round(rrange[2], 4)))), adj=0)
+# label the critical density
+rlevels <- seq(rrange[1], rrange[2], length.out=ncol+1)
+rho.critical <- 322
+icrit <- which.min(abs(rlevels-rho.critical))
+text(x[3], ykey[icrit], as.expression(substitute(x~kg/m^3~group("(", rho[c], ")"), list(x=rho.critical))), adj=0)
+
+#if(method=="IAPWS95") {
+# # the saturation line is very accurate but not quite perfect;
+# # we can show whether it is on the liquid or vapor side
+# ina <- is.na(P.sigma)
+# rho.sigma <- rho.IAPWS95(T[!ina], P.sigma[!ina])
+# col <- rep("blue", length(rho.sigma))
+# col[rho.sigma < 322] <- "red"
+# points(T[!ina], P.sigma[!ina], col=col, pch=20)
+#}
+
Modified: pkg/CHNOSZ/inst/NEWS
===================================================================
--- pkg/CHNOSZ/inst/NEWS 2015-05-19 14:38:46 UTC (rev 87)
+++ pkg/CHNOSZ/inst/NEWS 2015-06-09 10:21:09 UTC (rev 88)
@@ -1,3 +1,28 @@
+CHANGES IN CHNOSZ 1.0.5-1 (2015-06-09)
+--------------------------------------
+
+- Rewrite rho.IAPWS95() to be able to invert density from IAPWS-95
+ equations for a more extensive range of T,P values.
+
+- Update R dependency to R-3.1.0, needed for 'extendInt' argument
+ of uniroot() (used by rho.IAPWS95()).
+
+- Add demo/density.R to show density of H2O calculated using
+ rho.IAPWS95() (optionally using water.SUPCRT92()).
+
+- Remove the relatively trivial demo/cordierite.R.
+
+- Update wrapper function for parallel calculations palply() (export of
+ variables, used in demo/density.R) and add tests in
+ test-util.program.R.
+
+- Tests that initiate calls to palply() (and therefore makeCluster())
+ failed with 'cannot open file 'startup.Rs': No such file or directory'
+ Fixed by adding Sys.setenv("R_TESTS" = "") to test-all.R.
+ Issue discussed here:
+ https://github.com/hadley/testthat/issues/129
+ https://github.com/hadley/testthat/issues/144
+
CHANGES IN CHNOSZ 1.0.5 (2015-05-19)
------------------------------------
Added: pkg/CHNOSZ/inst/tests/test-util.program.R
===================================================================
--- pkg/CHNOSZ/inst/tests/test-util.program.R (rev 0)
+++ pkg/CHNOSZ/inst/tests/test-util.program.R 2015-06-09 10:21:09 UTC (rev 88)
@@ -0,0 +1,28 @@
+context("util.program")
+
+# these tests are inefficient uses of parallelization
+# (overhead is greater than the savings from multiple cores)
+# just here to test that the functions are working
+
+test_that("palply() launches calculations on multiple cores", {
+ if(min(getOption("mc.cores"), 2) > 1 & parallel::detectCores() > 1) {
+ x <- 1:1001
+ # for this we don't have to export any variables so varlist=""
+ expect_message(palply("", 1:length(x), function(i) i^2), "running 1001 calculations")
+ }
+})
+
+test_that("other functions are calling palply() properly", {
+ if(min(getOption("mc.cores"), 2) > 1 & parallel::detectCores() > 1) {
+ ff <- system.file("extdata/fasta/HTCC1062.faa.xz", package="CHNOSZ")
+ expect_message(aa <- read.fasta(ff), "read.fasta running 1354 calculations")
+ # ^^^ also messaged: count.aa running 1354 calculations
+ ip <- add.protein(aa)
+ basis("CHNOS")
+ expect_message(a <- affinity(O2=c(-90, -60, 1000), iprotein=ip), "affinity running 1354 calculations")
+ expect_message(e <- equilibrate(a), "equil.reaction running 1000 calculations")
+ expect_message(e <- equilibrate(a, normalize=TRUE), "equil.boltzmann running 1354 calculations")
+ # ^^^ above message repeated 2x
+ }
+})
+
Modified: pkg/CHNOSZ/man/CHNOSZ-package.Rd
===================================================================
--- pkg/CHNOSZ/man/CHNOSZ-package.Rd 2015-05-19 14:38:46 UTC (rev 87)
+++ pkg/CHNOSZ/man/CHNOSZ-package.Rd 2015-06-09 10:21:09 UTC (rev 88)
@@ -82,7 +82,8 @@
\section{Compatibility}{
- As of version 1.0.4, the package depends on \R version 3.0.0 or greater (previous versions use \code{\link{Stangle}} to extract R code from vignettes when installing the source package, leading to failure processing hotspring.Rnw, which now uses knitr instead of Sweave).
+ Starting with version 1.0.5-1, the package depends on \R version 3.1.0 or greater (for \code{extendInt} argument of \code{\link{uniroot}}, used in code{\link{rho.IAPWS95}}).
+ As of version 1.0.4 (release 1.0.5 on CRAN), the package depends on \R version 3.0.0 or greater (previous versions use \code{\link{Stangle}} to extract R code from vignettes when installing the source package, leading to failure processing hotspring.Rnw, which now uses knitr instead of Sweave).
Before version 1.0.4, the recommended version of \R was 2.14.0 or greater (to find vignettes in the \code{vignettes} directory, and for availability of \pkg{parallel} in the standard library).
As of version 0.9-9, the package depends on \R version 2.12.0 or greater (so useDynLib in the NAMESPACE can find the shared library on Windows).
Starting with version 0.9-6 of the package, the dependency was given as \R version 2.10.0 or greater (to read compressed data files).
Modified: pkg/CHNOSZ/man/IAPWS95.Rd
===================================================================
--- pkg/CHNOSZ/man/IAPWS95.Rd 2015-05-19 14:38:46 UTC (rev 87)
+++ pkg/CHNOSZ/man/IAPWS95.Rd 2015-06-09 10:21:09 UTC (rev 88)
@@ -3,7 +3,6 @@
\alias{IAPWS95}
\alias{IAPWS95.idealgas}
\alias{IAPWS95.residual}
-\alias{WP02.auxiliary}
\title{Properties of Water from IAPWS-95}
\description{
Calculate thermodynamic properties of water following the IAPWS-95 formulation.
@@ -13,7 +12,6 @@
IAPWS95(property, T = 298.15, rho = 1000)
IAPWS95.idealgas(p, delta, tau)
IAPWS95.residual(p, delta, tau)
- WP02.auxiliary(property, T = 298.15)
}
\arguments{
@@ -32,8 +30,6 @@
The \code{IAPWS95} function returns values of thermodynamic properties in specific units (per gram).
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/chnosz -r 88
More information about the CHNOSZ-commits
mailing list