[CHNOSZ-commits] r254 - in pkg/CHNOSZ: . R data demo inst man tests/testthat vignettes
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Oct 12 16:56:09 CEST 2017
Author: jedick
Date: 2017-10-12 16:56:08 +0200 (Thu, 12 Oct 2017)
New Revision: 254
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/R/nonideal.R
pkg/CHNOSZ/R/subcrt.R
pkg/CHNOSZ/data/opt.csv
pkg/CHNOSZ/demo/ORP.R
pkg/CHNOSZ/inst/NEWS
pkg/CHNOSZ/man/nonideal.Rd
pkg/CHNOSZ/man/water.Rd
pkg/CHNOSZ/tests/testthat/test-affinity.R
pkg/CHNOSZ/tests/testthat/test-nonideal.R
pkg/CHNOSZ/vignettes/anintro.Rmd
Log:
nonideal(): add 'Helgeson' method
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2017-10-12 09:10:59 UTC (rev 253)
+++ pkg/CHNOSZ/DESCRIPTION 2017-10-12 14:56:08 UTC (rev 254)
@@ -1,6 +1,6 @@
Date: 2017-10-12
Package: CHNOSZ
-Version: 1.1.0-52
+Version: 1.1.0-53
Title: Thermodynamic Calculations for Geobiochemistry
Author: Jeffrey Dick
Maintainer: Jeffrey Dick <j3ffdick at gmail.com>
Modified: pkg/CHNOSZ/R/nonideal.R
===================================================================
--- pkg/CHNOSZ/R/nonideal.R 2017-10-12 09:10:59 UTC (rev 253)
+++ pkg/CHNOSZ/R/nonideal.R 2017-10-12 14:56:08 UTC (rev 254)
@@ -3,14 +3,24 @@
# moved to nonideal.R from util.misc.R 20151107
# added Helgeson method 20171012
-nonideal <- function(species, proptable, IS, T, P) {
- thermo <- get("thermo")
+nonideal <- function(species, speciesprops, IS, T, P, A_DH, B_DH, method=get("thermo")$opt$nonideal) {
# generate nonideal contributions to thermodynamic properties
- # number of species, same length as proptable list
- # T in Kelvin, same length as nrows of proptables
+ # number of species, same length as speciesprops list
+ # T in Kelvin, same length as nrows of speciespropss
+ # arguments P, A_DH, B_DH are needed for Helgeson method only
+ # we can use this function to change the nonideal method option
+ if(identical(species, "Helgeson") | identical(species, "Alberty")) {
+ thermo <- get("thermo")
+ oldnon <- thermo$opt$nonideal
+ thermo$opt$nonideal <- species
+ assign("thermo", thermo, "CHNOSZ")
+ message("nonideal: setting method option to ", species)
+ return(invisible(oldnon))
+ }
+
# function to calculate extended Debye-Huckel equation and derivatives using Alberty's parameters
- Alberty <- function(Z, I, T, prop = "log") {
+ Alberty <- function(Z, I, T, prop = "loggamma") {
# extended Debye-Huckel equation ("log")
# and its partial derivatives ("G","H","S","Cp")
# T in Kelvin
@@ -29,7 +39,7 @@
loggamma <- function(a, Z, I, B) - a * Z^2 * I^(1/2) / (1 + B * I^(1/2))
# TODO: check the following equations 20080208 jmd
R <- 1.9872 # gas constant, cal K^-1 mol^-1
- if(prop=="log") return(loggamma(eval(A), Z, I, B))
+ if(prop=="loggamma") return(loggamma(eval(A), Z, I, B))
else if(prop=="G") return(R * T * loggamma(eval(A), Z, I, B))
else if(prop=="H") return(R * T^2 * loggamma(eval(DD(A, "T", 1)), Z, I, B))
else if(prop=="S") return(- R * T * loggamma(eval(DD(A, "T", 1)), Z, I, B))
@@ -37,21 +47,27 @@
}
# function for Debye-Huckel equation with B-dot extended term parameter (Helgeson, 1969)
- Helgeson <- function() {
- # "distance of closest approach" of ions in NaCl solutions
- # HKF81 Table 2
- acirc <- 3.72 # Angstrom
+ Helgeson <- function(Z, I, T, P, A_DH, B_DH, prop = "loggamma") {
+ # "distance of closest approach" of ions in NaCl solutions (HKF81 Table 2)
+ acirc <- 3.72e-8 # cm
+ loggamma <- - A_DH * Z^2 * I^0.5 / (1 + acirc * B_DH * I^0.5) + Bdot * I
+ R <- 1.9872 # gas constant, cal K^-1 mol^-1
+ if(prop=="loggamma") return(loggamma)
+ else if(prop=="G") return(R * T * loggamma)
}
+ # get B-dot if we're using the Helgeson method
+ if(method=="Helgeson") Bdot <- Bdot(convert(T, "C"), P)
+
# get species indices
if(!is.numeric(species[[1]])) species <- info(species, "aq")
iH <- info("H+")
ie <- info("e-")
- proptable <- as.list(proptable)
+ speciesprops <- as.list(speciesprops)
ndid <- 0
# loop over species
for(i in 1:length(species)) {
- myprops <- proptable[[i]]
+ myprops <- speciesprops[[i]]
# get the charge from the chemical formula
# force a charge count even if it's zero
mkp <- makeup(c("Z0", species[i]), sum=TRUE)
@@ -64,18 +80,26 @@
didit <- FALSE
for(j in 1:ncol(myprops)) {
pname <- colnames(myprops)[j]
- if(pname %in% c("G", "H", "S", "Cp")) {
- myprops[,j] <- myprops[,j] + Alberty(Z, IS, T, pname)
+ if(method=="Alberty" & pname %in% c("G", "H", "S", "Cp")) {
+ myprops[, j] <- myprops[, j] + Alberty(Z, IS, T, pname)
didit <- TRUE
+ } else if(method=="Helgeson" & pname %in% c("G", "H", "S", "Cp")) {
+ myprops[, j] <- myprops[, j] + Helgeson(Z, IS, T, P, A_DH, B_DH, pname)
+ didit <- TRUE
}
}
# append a loggam column if we did any nonideal calculations of thermodynamic properties
- if(didit) myprops <- cbind(myprops, loggam = Alberty(Z, IS, T))
- proptable[[i]] <- myprops
+ if(didit) {
+ if(method=="Alberty") myprops <- cbind(myprops, loggam = Alberty(Z, IS, T))
+ else if(method=="Helgeson") {
+ myprops <- cbind(myprops, loggam = Helgeson(Z, IS, T, P, A_DH, B_DH))
+ }
+ }
+ speciesprops[[i]] <- myprops
if(didit) ndid <- ndid + 1
}
if(ndid > 0) message(paste("nonideal:", ndid, "species were nonideal"))
- return(proptable)
+ return(speciesprops)
}
Bdot <- function(TC = 25, P = 1, showsplines = "") {
@@ -104,7 +128,7 @@
# values for Bdot x 100 from Helgeson (1969), Figure (P = Psat)
if(!isoP | showsplines != "") {
T0 <- c(23.8, 49.4, 98.9, 147.6, 172.6, 197.1, 222.7, 248.1, 268.7)
- B0 <- c(4.07, 4.27, 4.30, 4.62, 4.86, 4.73, 4.09, 3.61, 1.56)
+ B0 <- c(4.07, 4.27, 4.30, 4.62, 4.86, 4.73, 4.09, 3.61, 1.56) / 100
# we could use the values from Hel69 Table 2 but extrapolation of the
# their fitted spline function turns sharply upward above 300 degC
#T0a <- c(25, 50, 100, 150, 200, 250, 270, 300)
@@ -114,39 +138,39 @@
# values for bgamma x 100 from Helgeson et al., 1981 Table 27
if(is500 | !isoP | showsplines != "") {
T0.5 <- seq(0, 400, 25)
- B0.5 <- c(5.6, 7.1, 7.8, 8.0, 7.8, 7.5, 7.0, 6.4, 5.7, 4.8, 3.8, 2.6, 1.0, -1.2, -4.1, -8.4, -15.2)
+ B0.5 <- c(5.6, 7.1, 7.8, 8.0, 7.8, 7.5, 7.0, 6.4, 5.7, 4.8, 3.8, 2.6, 1.0, -1.2, -4.1, -8.4, -15.2) / 100
S0.5 <- splinefun(T0.5, B0.5)
if(is500) return(S0.5(T))
}
if(is1000 | !isoP | showsplines != "") {
T1 <- seq(0, 500, 25)
- B1 <- c(6.6, 7.7, 8.7, 8.3, 8.2, 7.9, 7.5, 7.0, 6.5, 5.9, 5.2, 4.4, 3.5, 2.5, 1.1, -0.6, -2.8, -5.7, -9.3, -13.7, -19.2)
+ B1 <- c(6.6, 7.7, 8.7, 8.3, 8.2, 7.9, 7.5, 7.0, 6.5, 5.9, 5.2, 4.4, 3.5, 2.5, 1.1, -0.6, -2.8, -5.7, -9.3, -13.7, -19.2) / 100
S1 <- splinefun(T1, B1)
if(is1000) return(S1(T))
}
if(is2000 | !isoP | showsplines != "") {
# 550 and 600 degC points from Manning et al., 2013 Fig. 11
T2 <- c(seq(0, 500, 25), 550, 600)
- B2 <- c(7.4, 8.3, 8.8, 8.9, 8.9, 8.7, 8.5, 8.1, 7.8, 7.4, 7.0, 6.6, 6.2, 5.8, 5.2, 4.6, 3.8, 2.9, 1.8, 0.5, -1.0, -3.93, -4.87)
+ B2 <- c(7.4, 8.3, 8.8, 8.9, 8.9, 8.7, 8.5, 8.1, 7.8, 7.4, 7.0, 6.6, 6.2, 5.8, 5.2, 4.6, 3.8, 2.9, 1.8, 0.5, -1.0, -3.93, -4.87) / 100
S2 <- splinefun(T2, B2)
if(is2000) return(S2(T))
}
if(is3000 | !isoP | showsplines != "") {
T3 <- seq(0, 500, 25)
- B3 <- c(6.5, 8.3, 9.2, 9.6, 9.7, 9.6, 9.4, 9.3, 9.2, 9.0, 8.8, 8.6, 8.3, 8.1, 7.8, 7.5, 7.1, 6.6, 6.0, 5.4, 4.8)
+ B3 <- c(6.5, 8.3, 9.2, 9.6, 9.7, 9.6, 9.4, 9.3, 9.2, 9.0, 8.8, 8.6, 8.3, 8.1, 7.8, 7.5, 7.1, 6.6, 6.0, 5.4, 4.8) / 100
S3 <- splinefun(T3, B3)
if(is3000) return(S3(T))
}
if(is4000 | !isoP | showsplines != "") {
T4 <- seq(0, 500, 25)
- B4 <- c(4.0, 7.7, 9.5, 10.3, 10.7, 10.8, 10.8, 10.8, 10.7, 10.6, 10.5, 10.4, 10.3, 10.2, 10.0, 9.8, 9.6, 9.3, 8.9, 8.5, 8.2)
+ B4 <- c(4.0, 7.7, 9.5, 10.3, 10.7, 10.8, 10.8, 10.8, 10.7, 10.6, 10.5, 10.4, 10.3, 10.2, 10.0, 9.8, 9.6, 9.3, 8.9, 8.5, 8.2) / 100
S4 <- splinefun(T4, B4)
if(is4000) return(S4(T))
}
if(is5000 | !isoP | showsplines != "") {
# 550 and 600 degC points from Manning et al., 2013 Fig. 11
T5 <- c(seq(0, 500, 25), 550, 600)
- B5 <- c(0.1, 6.7, 9.6, 11.1, 11.8, 12.2, 12.4, 12.4, 12.4, 12.4, 12.4, 12.3, 12.3, 12.2, 12.1, 11.9, 11.8, 11.5, 11.3, 11.0, 10.8, 11.2, 12.52)
+ B5 <- c(0.1, 6.7, 9.6, 11.1, 11.8, 12.2, 12.4, 12.4, 12.4, 12.4, 12.4, 12.3, 12.3, 12.2, 12.1, 11.9, 11.8, 11.5, 11.3, 11.0, 10.8, 11.2, 12.52) / 100
S5 <- splinefun(T5, B5)
if(is5000) return(S5(T))
}
@@ -154,44 +178,44 @@
# here, one control point at 10 degC is added to make the splines curve down at low T
if(is10000 | !isoP | showsplines != "") {
T10 <- c(25, seq(300, 1000, 50))
- B10 <- c(12, 17.6, 17.8, 18, 18.2, 18.9, 21, 23.3, 26.5, 28.8, 31.4, 34.1, 36.5, 39.2, 41.6, 44.1)
+ B10 <- c(12, 17.6, 17.8, 18, 18.2, 18.9, 21, 23.3, 26.5, 28.8, 31.4, 34.1, 36.5, 39.2, 41.6, 44.1) / 100
S10 <- splinefun(T10, B10)
if(is10000) return(S10(T))
}
if(is20000 | !isoP | showsplines != "") {
T20 <- c(25, seq(300, 1000, 50))
- B20 <- c(16, 21.2, 21.4, 22, 22.4, 23.5, 26.5, 29.2, 32.6, 35.2, 38.2, 41.4, 44.7, 47.7, 50.5, 53.7)
+ B20 <- c(16, 21.2, 21.4, 22, 22.4, 23.5, 26.5, 29.2, 32.6, 35.2, 38.2, 41.4, 44.7, 47.7, 50.5, 53.7) / 100
S20 <- splinefun(T20, B20)
if(is20000) return(S20(T))
}
if(is30000 | !isoP | showsplines != "") {
T30 <- c(25, seq(300, 1000, 50))
- B30 <- c(19, 23.9, 24.1, 24.6, 25.2, 26.7, 30.3, 32.9, 36.5, 39.9, 43, 46.4, 49.8, 53.2, 56.8, 60)
+ B30 <- c(19, 23.9, 24.1, 24.6, 25.2, 26.7, 30.3, 32.9, 36.5, 39.9, 43, 46.4, 49.8, 53.2, 56.8, 60) / 100
S30 <- splinefun(T30, B30)
if(is30000) return(S30(T))
}
# 40-60 kb points extrapolated from 10-30 kb points of Manning et al., 2013
if(is40000 | !isoP | showsplines != "") {
T40 <- c(seq(300, 1000, 50))
- B40 <- c(25.8, 26, 26.4, 27.2, 28.9, 33, 35.5, 39.2, 43.2, 46.4, 49.9, 53.4, 57.1, 61.2, 64.4)
+ B40 <- c(25.8, 26, 26.4, 27.2, 28.9, 33, 35.5, 39.2, 43.2, 46.4, 49.9, 53.4, 57.1, 61.2, 64.4) / 100
S40 <- splinefun(T40, B40)
if(is40000) return(S40(T))
}
if(is50000 | !isoP | showsplines != "") {
T50 <- c(seq(300, 1000, 50))
- B50 <- c(27.1, 27.3, 27.7, 28.5, 30.5, 34.8, 37.3, 41.1, 45.5, 48.7, 52.4, 55.9, 59.8, 64.3, 67.5)
+ B50 <- c(27.1, 27.3, 27.7, 28.5, 30.5, 34.8, 37.3, 41.1, 45.5, 48.7, 52.4, 55.9, 59.8, 64.3, 67.5) / 100
S50 <- splinefun(T50, B50)
if(is50000) return(S50(T))
}
if(is60000 | !isoP | showsplines != "") {
T60 <- c(seq(300, 1000, 50))
- B60 <- c(28, 28.2, 28.6, 29.5, 31.6, 36.1, 38.6, 42.5, 47.1, 50.4, 54.1, 57.6, 61.6, 66.5, 69.7)
+ B60 <- c(28, 28.2, 28.6, 29.5, 31.6, 36.1, 38.6, 42.5, 47.1, 50.4, 54.1, 57.6, 61.6, 66.5, 69.7) / 100
S60 <- splinefun(T60, B60)
if(is60000) return(S60(T))
}
# show points and spline(T) curves
if(showsplines == "T") {
- thermo.plot.new(c(0, 1000), c(-20, 70), xlab=axis.label("T"), ylab=expression(b[gamma]%*%100))
+ thermo.plot.new(c(0, 1000), c(-.2, .7), xlab=axis.label("T"), ylab=expression(b[gamma]))
points(T0, B0, pch=0)
points(T0.5, B0.5, pch=1)
points(T1, B1, pch=1)
@@ -229,12 +253,11 @@
legend("bottomright", col=c(NA, rev(col)), lty=1,
legend=c("kbar", "60", "50", "40", "30", "20", "10", "5", "4", "3", "2", "1", "0.5", "Psat"))
} else if(showsplines=="P") {
- #thermo.plot.new(c(0, 60000), c(-20, 70), xlab=axis.label("P"), ylab=expression(b[gamma]%*%100))
- thermo.plot.new(c(0, 5), c(-20, 70), xlab="log P", ylab=expression(b[gamma]%*%100))
+ thermo.plot.new(c(0, 5), c(-.2, .7), xlab="log P(bar)", ylab=expression(b[gamma]))
# pressures that are used to make the isothermal splines (see below)
P25 <- c(1, 500, 1000, 2000, 3000, 4000, 5000)
P100 <- c(1, 500, 1000, 2000, 3000, 4000, 5000, 10000, 20000)
- P200 <- c(15, 500, 1000, 2000, 3000, 4000, 5000, 10000, 20000, 30000, 40000)
+ P200 <- c(16, 500, 1000, 2000, 3000, 4000, 5000, 10000, 20000, 30000, 40000)
P300 <- c(86, 500, 1000, 2000, 3000, 4000, 5000, 10000, 20000, 30000, 40000, 50000, 60000)
P400 <- c(500, 1000, 2000, 3000, 4000, 5000, 10000, 20000, 30000, 40000, 50000, 60000)
P500 <- c(1000, 2000, 3000, 4000, 5000, 10000, 20000, 30000, 40000, 50000, 60000)
@@ -244,30 +267,30 @@
P900 <- c(10000, 20000, 30000, 40000, 50000, 60000)
P1000 <- c(10000, 20000, 30000, 40000, 50000, 60000)
# plot the pressure and B-dot values used to make the isothermal splines
- points(log10(P25), Bdot(25, P25) * 100)
- points(log10(P100), Bdot(100, P100) * 100)
- points(log10(P200), Bdot(200, P200) * 100)
- points(log10(P300), Bdot(300, P300) * 100)
- points(log10(P400), Bdot(400, P400) * 100)
- points(log10(P500), Bdot(500, P500) * 100)
- points(log10(P600), Bdot(600, P600) * 100)
- points(log10(P700), Bdot(700, P700) * 100)
- points(log10(P800), Bdot(800, P800) * 100)
- points(log10(P900), Bdot(900, P900) * 100)
- points(log10(P1000), Bdot(1000, P1000) * 100)
+ points(log10(P25), Bdot(25, P25))
+ points(log10(P100), Bdot(100, P100))
+ points(log10(P200), Bdot(200, P200))
+ points(log10(P300), Bdot(300, P300))
+ points(log10(P400), Bdot(400, P400))
+ points(log10(P500), Bdot(500, P500))
+ points(log10(P600), Bdot(600, P600))
+ points(log10(P700), Bdot(700, P700))
+ points(log10(P800), Bdot(800, P800))
+ points(log10(P900), Bdot(900, P900))
+ points(log10(P1000), Bdot(1000, P1000))
# plot the isothermal spline functions
col <- tail(rev(rainbow(12)), -1)
- P <- c(1, seq(50, 5000, 50)); lines(log10(P), Bdot(25, P) * 100, col=col[1])
- P <- c(1, seq(50, 20000, 50)); lines(log10(P), Bdot(100, P) * 100, col=col[2])
- P <- c(1, seq(50, 40000, 50)); lines(log10(P), Bdot(200, P) * 100, col=col[3])
- P <- c(1, seq(50, 60000, 50)); lines(log10(P), Bdot(300, P) * 100, col=col[4])
- P <- seq(500, 60000, 50); lines(log10(P), Bdot(400, P) * 100, col=col[5])
- P <- seq(1000, 60000, 50); lines(log10(P), Bdot(500, P) * 100, col=col[6])
- P <- seq(2000, 60000, 50); lines(log10(P), Bdot(600, P) * 100, col=col[7])
- P <- seq(10000, 60000, 50); lines(log10(P), Bdot(700, P) * 100, col=col[8])
- P <- seq(10000, 60000, 50); lines(log10(P), Bdot(800, P) * 100, col=col[9])
- P <- seq(10000, 60000, 50); lines(log10(P), Bdot(900, P) * 100, col=col[10])
- P <- seq(10000, 60000, 50); lines(log10(P), Bdot(1000, P) * 100, col=col[11])
+ P <- c(1, seq(50, 5000, 50)); lines(log10(P), Bdot(25, P), col=col[1])
+ P <- c(1, seq(50, 20000, 50)); lines(log10(P), Bdot(100, P), col=col[2])
+ P <- c(1, seq(50, 40000, 50)); lines(log10(P), Bdot(200, P), col=col[3])
+ P <- c(1, seq(50, 60000, 50)); lines(log10(P), Bdot(300, P), col=col[4])
+ P <- seq(500, 60000, 50); lines(log10(P), Bdot(400, P), col=col[5])
+ P <- seq(1000, 60000, 50); lines(log10(P), Bdot(500, P), col=col[6])
+ P <- seq(2000, 60000, 50); lines(log10(P), Bdot(600, P), col=col[7])
+ P <- seq(10000, 60000, 50); lines(log10(P), Bdot(700, P), col=col[8])
+ P <- seq(10000, 60000, 50); lines(log10(P), Bdot(800, P), col=col[9])
+ P <- seq(10000, 60000, 50); lines(log10(P), Bdot(900, P), col=col[10])
+ P <- seq(10000, 60000, 50); lines(log10(P), Bdot(1000, P), col=col[11])
legend("topleft", col=c(NA, col), lty=1, legend=c("degrees C", 25, 100, 200, 300, 400, 500, 600, 700, 800, 900, 1000))
legend("bottomright", pch=1, legend="points from iso-P splines")
} else {
@@ -316,7 +339,7 @@
# remember this T; if it's the same as the next one, we won't re-make the spline
lastT <- T[i]
}
- Bdot <- c(Bdot, ST(P[i]) / 100)
+ Bdot <- c(Bdot, ST(P[i]))
}
}
return(Bdot)
Modified: pkg/CHNOSZ/R/subcrt.R
===================================================================
--- pkg/CHNOSZ/R/subcrt.R 2017-10-12 09:10:59 UTC (rev 253)
+++ pkg/CHNOSZ/R/subcrt.R 2017-10-12 14:56:08 UTC (rev 254)
@@ -466,10 +466,11 @@
iscgl <- iscgl.new
isH2O <- isH2O.new
- # the order of the properties
- if(length(calcprop)>1) for(i in 1:length(out)) {
- # keep state/loggam columns if they exists
+ # adjust the output order of the properties
+ for(i in 1:length(out)) {
+ # the calculated properties are first
ipp <- match(calcprop, colnames(out[[i]]))
+ # move polymorph/loggam columns to end
if('polymorph' %in% colnames(out[[i]])) ipp <- c(ipp,match('polymorph',colnames(out[[i]])))
if('loggam' %in% colnames(out[[i]])) ipp <- c(ipp,match('loggam',colnames(out[[i]])))
out[[i]] <- out[[i]][,ipp,drop=FALSE]
@@ -492,9 +493,9 @@
# then to the user's units (outvert) from cal
A <- outvert(convert(-A,'G',T=T),'cal')
}
- # the addition of properties
+ # loop over reaction coefficients
for(i in 1:length(coeff)) {
- # assemble state columns if they exist
+ # assemble polymorph columns separately
if('polymorph' %in% colnames(out[[i]])) {
sc <- as.data.frame(out[[i]]$polymorph)
out[[i]] <- out[[i]][,-match('polymorph',colnames(out[[i]]))]
@@ -502,13 +503,13 @@
if(is.null(morphcols)) morphcols <- sc
else morphcols <- cbind(morphcols,sc)
}
- # include a zero loggam column if we need it
- # for those species that are ideal
+ # include a zero loggam column if needed (for those species that are ideal)
o.i <- out[[i]]
if('loggam' %in% colnames(o.i)) if(!'loggam' %in% colnames(o))
o <- cbind(o,loggam=0)
if('loggam' %in% colnames(o)) if(!'loggam' %in% colnames(o.i))
o.i <- cbind(o.i,loggam=0)
+ # the real addition of properties
o <- o + o.i * coeff[i]
}
# output for reaction (stack on polymorph columns if exist)
Modified: pkg/CHNOSZ/data/opt.csv
===================================================================
--- pkg/CHNOSZ/data/opt.csv 2017-10-12 09:10:59 UTC (rev 253)
+++ pkg/CHNOSZ/data/opt.csv 2017-10-12 14:56:08 UTC (rev 254)
@@ -1,2 +1,2 @@
-cutoff,E.units,T.units,P.units,state,water,G.tol,Cp.tol,V.tol,varP,IAPWS.sat,paramin,ideal.H,ideal.e
-1E-10,cal,C,bar,aq,SUPCRT92,100,1,1,FALSE,liquid,1000,TRUE,TRUE
+cutoff,E.units,T.units,P.units,state,water,G.tol,Cp.tol,V.tol,varP,IAPWS.sat,paramin,ideal.H,ideal.e,nonideal
+1E-10,cal,C,bar,aq,SUPCRT92,100,1,1,FALSE,liquid,1000,TRUE,TRUE,Helgeson
Modified: pkg/CHNOSZ/demo/ORP.R
===================================================================
--- pkg/CHNOSZ/demo/ORP.R 2017-10-12 09:10:59 UTC (rev 253)
+++ pkg/CHNOSZ/demo/ORP.R 2017-10-12 14:56:08 UTC (rev 254)
@@ -1,10 +1,12 @@
-# yell2010/orp.R 20100715 jmd
+# CHNOSZ/demo/ORP.R
+# first version 20100715 jmd
+
# calculate the temperature dependence of
# potentials vs. standard hydrogen electrode (SHE) of various electrodes (Ag/AgCl)
# and ORP standards (ZoBell, Light's, (tri)iodide)
-# CHNOSZ provides functions subcrt() and convert()
-# used in this example
+# use the extended Debye-Huckel equation with Alberty's parameters
+oldnon <- nonideal("Alberty")
# Bard et al.'s fit to the potential
# (Bard, Parson, Jordan, Standard Potentials In Aqueous Solution, 1985)
@@ -167,3 +169,6 @@
# finally, make the plot
figure()
+
+# reset the nonideality method to default
+nonideal(oldnon)
Modified: pkg/CHNOSZ/inst/NEWS
===================================================================
--- pkg/CHNOSZ/inst/NEWS 2017-10-12 09:10:59 UTC (rev 253)
+++ pkg/CHNOSZ/inst/NEWS 2017-10-12 14:56:08 UTC (rev 254)
@@ -1,4 +1,4 @@
-CHANGES IN CHNOSZ 1.1.0-51 (2017-10-12)
+CHANGES IN CHNOSZ 1.1.0-53 (2017-10-12)
---------------------------------------
MAJOR CHANGES:
@@ -19,7 +19,7 @@
- Usage of the DEW model is shown in the new demo DEW.R. This demo also
depends on the Berman equations (above) and, for the last diagram, the
- following two NEWS items:
+ following *two* NEWS items:
- In equilibrate(), it is now possible to combine affinity calculations
with changing activity of the balancing basis species (loga.balance).
@@ -35,6 +35,11 @@
"percent carbon" plots for systems where the species have different
carbon numbers.
+- nonideal() now has two methods for calculating activity coefficients:
+ 'Helgeson', which utilizes the "B-dot" equation of Helgeson, 1969,
+ and 'Alberty' (the only method previously available). The 'Helgeson'
+ method depends on the following *two* NEWS items:
+
- All three water options (SUPCRT92, IAPWS95, DEW) can be used to
calculate 'A_DH' and 'B_DH', i.e. A and B coefficients in the extended
Debye-Huckel equation of Helgeson, 1969.
Modified: pkg/CHNOSZ/man/nonideal.Rd
===================================================================
--- pkg/CHNOSZ/man/nonideal.Rd 2017-10-12 09:10:59 UTC (rev 253)
+++ pkg/CHNOSZ/man/nonideal.Rd 2017-10-12 14:56:08 UTC (rev 254)
@@ -4,34 +4,43 @@
\alias{Bdot}
\title{Activity coefficients of aqueous species}
\description{
-Calculate activity coefficients and non-ideal contributions to apparent standard molal properties of aqueous species.
+Calculate activity coefficients and apparent (non-ideal) molal properties of aqueous species.
}
\usage{
- nonideal(species, proptable, IS, T)
+ nonideal(species, speciesprops, IS, T, P, A_DH, B_DH, method=get("thermo")$opt$nonideal)
Bdot(TC, P, showsplines = "")
}
\arguments{
\item{species}{names or indices of species for which to calculate nonideal properties}
- \item{proptable}{list of dataframes of species properties}
+ \item{speciesprops}{list of dataframes of species properties}
\item{IS}{numeric, ionic strength(s) used in nonideal calculations, mol kg\eqn{^{-1}}{^-1}}
\item{T}{numeric, temperature (K)}
+ \item{P}{numeric, pressure (bar); required for Helgeson method}
+ \item{A_DH}{numeric, A Debye-Huckel coefficient; required for Helgeson method}
+ \item{B_DH}{numeric, B Debye-Huckel coefficient; required for Helgeson method}
+ \item{method}{character, \samp{Alberty} or \samp{Helgeson}}
\item{TC}{numeric, temperature (\degC)}
- \item{P}{numeric, pressure (bar)}
\item{showsplines}{character, show isobaric (\samp{T}) or isothermal (\samp{P}) splines}
}
\details{
-\code{nonideal} takes a list of dataframes (in \code{proptable}) containing the standard molal properties of the identified \code{species}.
+\code{nonideal} takes a list of dataframes (in \code{speciesprops}) containing the standard molal properties of the identified \code{species}.
+The function calculates the *apparent* properties for given ionic strength (\code{IS}); they are equal to the *standard* values at IS=0.
The function bypasses (leaves unchanged) properties of all species whose charge (determined by the number of Z in their \code{\link{makeup}}) is equal to zero.
The proton (H+) and electron (e-) are also bypassed by default; to apply the calculations to H+ and/or e-, change \code{\link{thermo}$opt$ideal.H} or \code{ideal.e} to FALSE.
-The values of \code{IS} are combined with Alberty's (2003) equation 3.6-1 (extended Debye-Hückel equation) and its derivatives, to calculate apparent molal properties at the specified ionic strength(s) and temperature(s).
-The lengths of \code{IS} and \code{T} supplied in the arguments should be equal to the number of rows of each dataframe in \code{proptable}, or one to use single values throughout.
-The apparent molal properties that can be calculated include \code{G}, \code{H}, \code{S} and \code{Cp}; any columns in the dataframes of \code{proptable} with other names are left untouched.
-A column named \code{loggam} (logarithm of gamma, the activity coefficient) is appended to the output dataframe of species properties.
+The lengths of \code{IS} and \code{T} supplied in the arguments should be equal to the number of rows of each dataframe in \code{speciesprops}, or length one to use single values throughout.
-\code{Bdot} calculates the \dQuote{B-dot} deviation function (Helgeson, 1969) or extended term parameter (written as b_gamma; Helgeson et al., 1981) for activity coefficients in NaCl solutions at high temperature and pressure.
+If \code{method} is \samp{Alberty}, the values of \code{IS} are combined with Alberty's (2003) equation 3.6-1 (extended Debye-Hückel equation) and its derivatives, to calculate apparent molal properties at the specified ionic strength(s) and temperature(s).
+The apparent molal properties that can be calculated include \samp{G}, \samp{H}, \samp{S} and \samp{Cp}; any columns in the dataframes of \code{speciesprops} with other names are left untouched.
+
+If \code{method} is \samp{Helgeson}, the \dQuote{B-dot} equation (Helgeson, 1969; Helgeson et al., 1981; Manning, 2013) is used.
+In addition to \code{IS} and \code{T}, this method depends on values of \code{P}, \code{A_DH}, and \code{B_DH} given in the arguments.
+The calculation of \dQuote{B-dot}, also used in the equations, is made within \code{nonideal} by calling the \code{Bdot} function.
+Currently, \samp{G} is the only apparent molal property that is calculated (but this can be used by \code{\link{subcrt}} to calculate apparent equilibrium constants).
+
+\code{Bdot} calculates the \dQuote{B-dot} deviation function (Helgeson, 1969) a.k.a. extended term parameter (written as b_gamma; Helgeson et al., 1981) for activity coefficients in NaCl solutions at high temperature and pressure.
Data at Psat and 0.5 to 5 kb are taken from Helgeson (1969, Table 2 and Figure 3) and Helgeson et al. (1981, Table 27) and extrapolated values at 10 to 30 kb from Manning et al. (2013, Figure 11).
Furthermore, the 10 to 30 kb data were used to generate super-extrapolated values at 40, 50, and 60 kb, which may be encountered using the \code{\link{water.DEW}} calculations.
If all \code{P} correspond to one of the isobaric conditions, the values of \code{Bdot} at \code{T} are calculated by spline fits to the isobaric data.
@@ -46,51 +55,59 @@
Note that the first example below uses \code{loggam} returned by \code{subcrt}, therefore requiring a base of 10 for calculating gamma.
}
+\value{
+One (\samp{G}) or more (\samp{H}, \samp{S}, \samp{Cp}; currently only with the Alberty method) standard thermodynamic properties (at IS=0) in \code{speciesprops} are replaced by the corresponding apparent thermodynamic properties (at higher IS).
+For all affected species, a column named \code{loggam} (logarithm of gamma, the activity coefficient) is appended to the output dataframe of species properties.
+}
+
\examples{\dontshow{data(thermo)}
-### Examples following Alberty, 2003
+### Examples following Alberty, 2003
+### (page numbers given below)
-## p. 273-276: activity coefficient (gamma)
-## as a function of ionic strength and temperature
+## the default method setting is Helgeson;
+## change it to Alberty
+oldnon <- nonideal("Alberty")
+
+## using nonideal() directly
+# p. 273-276: activity coefficient (gamma)
+# as a function of ionic strength and temperature
+IS <- seq(0, 0.25, 0.005)
T <- c(0, 25, 40)
-col <- c("blue", "black", "red")
-IS <- seq(0, 0.25, 0.0025)
-thermo.plot.new(xlim=range(IS), ylim=c(0, 1), xlab=axis.label("IS"), ylab="gamma")
+lty <- 1:3
+species <- c("H2PO4-", "HADP-2", "HATP-3", "ATP-4")
+col <- rainbow(4)
+thermo.plot.new(xlim = range(IS), ylim = c(0, 1), xlab = axis.label("IS"), ylab = "gamma")
for(j in 1:3) {
- s <- subcrt(c("H2PO4-", "HADP-2", "HATP-3", "ATP-4"), IS=IS, grid="IS", T=T[j])
- for(i in 1:4) lines(IS, 10^s$out[[i]]$loggam, col=col[j])
+ # use subcrt to generate speciesprops
+ speciesprops <- subcrt(species, T = rep(T[j], length(IS)))$out
+ # use nonideal to calculate loggamma; this also adjusts G, H, S, Cp, but we don't use them here
+ nonidealprops <- nonideal(species, speciesprops, IS = IS, T = convert(T[j], "K"))
+ for(i in 1:4) lines(IS, exp(nonidealprops[[i]]$loggam), lty=lty[j], col=col[i])
}
-text(0.125, 0.8, "Z = -1")
-text(0.1, 0.42, "Z = -2")
-text(0.075, 0.18, "Z = -3")
-text(0.05, 0.08, "Z = -4")
-title(main=paste("activity coefficient (gamma) of -1,-2,-3,-4",
- "charged species at 0, 25, 40 deg C, after Alberty, 2003",
- sep="\n"), cex.main=0.95)
-legend("topright", lty=c(NA, 1, 1, 1), col=c(NA, "blue", "black", "red"),
+t1 <- "activity coefficient (gamma) of -1,-2,-3,-4 charged species"
+t2 <- quote("at 0, 25, and 40 "*degree*"C, after Alberty, 2003")
+mtitle(as.expression(c(t1, t2)))
+legend("topright", lty=c(NA, 1:3), bty="n",
legend=c(as.expression(axis.label("T")), 0, 25, 40))
+legend("top", lty=1, col=col, bty="n",
+ legend = as.expression(lapply(species, expr.species)))
-## p. 16 Table 1.3: apparent pKa of acetic acid with
-## changing ionic strength
-# we set this option to FALSE so that nonideal() will calculate activity
-# coefficients for the proton (makes for better replication of the values
-# in Alberty's book)
+## more often, the 'IS' argument of subcrt() is used to compute
+## apparent properties at given ionic strength
+# p. 16 Table 1.3: apparent pKa of acetic acid
+# set ideal.H to FALSE to calculate activity coefficients for the proton
+# (makes for better replication of the values in Alberty's book)
thermo$opt$ideal.H <<- FALSE
-subcrt(c("acetic acid", "acetate", "H+"), c(-1, 1, 1),
- IS=c(0, 0.1, 0.25), T=25, property="logK")
-# note that *apparent* values equal *standard* values at IS=0
+(sres <- subcrt(c("acetate", "H+", "acetic acid"), c(-1, -1, 1),
+ IS=c(0, 0.1, 0.25), T=25, property="logK"))
+# we're within 0.01 log of Alberty's pK values
+Alberty_logK <- c(4.75, 4.54, 4.47)
+stopifnot(maxdiff(sres$out$logK, Alberty_logK) < 0.01)
# reset option to default
thermo$opt$ideal.H <<- TRUE
-## p. 95: basis and elemental stoichiometries of species
-# (this example doesn't use activity coefficients)
-basis(c("ATP-4", "H+", "H2O", "HPO4-2", "O2", "NH3"))
-# cf Eq. 5.1-33: basis composition
-species(c("ATP-4", "H+", "H2O", "HPO4-2", "ADP-3", "HATP-3", "HADP-2",
- "H2PO4-"))
-
-### A different example
-
-# speciation of phosphate as a function of ionic strength
+### An example using IS with affinity():
+## speciation of phosphate as a function of ionic strength
opar <- par(mfrow=c(2, 1))
basis("CHNOPS+")
Ts <- c(25, 100)
@@ -99,21 +116,47 @@
a <- affinity(IS=c(0, 0.14), T=T)
e <- equilibrate(a)
if(T==25) diagram(e, ylim=c(-3.0, -2.6), legend.x=NULL)
- else d <- diagram(e, ylim=c(-3.0, -2.6), add=TRUE, col="red")
+ else diagram(e, ylim=c(-3.0, -2.6), add=TRUE, col="red")
}
title(main="Non-ideality model for phosphate species")
dp <- describe.property(c("pH", "T", "T"), c(7, Ts))
legend("topright", lty=c(NA, 1, 1), col=c(NA, "black", "red"), legend=dp)
text(0.07, -2.76, expr.species("HPO4-2"))
text(0.07, -2.90, expr.species("H2PO4-"))
-#
-# phosphate predominance f(IS,pH)
+## phosphate predominance f(IS,pH)
a <- affinity(IS=c(0, 0.14), pH=c(6, 13), T=Ts[1])
d <- diagram(a, fill=NULL)
a <- affinity(IS=c(0, 0.14), pH=c(6, 13), T=Ts[2])
d <- diagram(a, add=TRUE, names=NULL, col="red")
par(opar)
+
+### finished with Alberty equation, let's look at Helgeson
+# this is the default setting, but is needed here because
+# we set "Alberty" above
+nonideal(oldnon) # same as nonideal("Helgeson")
+
+## activity coefficients for monovalent ions at 700 degC, 10 kbar
+# after Manning, 2010, Fig. 7
+IS <- c(0.001, 0.01, 0.1, 1, 2, 2.79)
+# we're above 5000 bar, so need to use IAPWS-95 or DEW
+oldwat <- water("DEW")
+wprop <- water(c("A_DH", "B_DH"), T=convert(700, "K"), P=10000)
+water(oldwat)
+# just use an empty table for a single species
+speciesprops <- list(data.frame(G=rep(0, length(IS))))
+# choose any monovalent species
+(nonidealprops <- nonideal("Na+", speciesprops, IS=IS,
+ T=convert(700, "K"), P=10000, A_DH=wprop$A_DH, B_DH=wprop$B_DH))
+# we get the nonideal Gibbs energy contribution and
+# the activity coefficient; check values of the latter
+Manning_gamma <- c(0.93, 0.82, 0.65, 0.76, 1.28, 2)
+gamma <- 10^nonidealprops[[1]]$loggam
+# we're getting progressively further from his values with
+# higher IS; not sure why
+stopifnot(maxdiff(gamma[1], Manning_gamma[1]) < 0.01)
+stopifnot(maxdiff(gamma, Manning_gamma) < 0.23)
+
## data and splines used for calculating B-dot
## (extended term parameter)
Bdot(showsplines = "T")
@@ -127,6 +170,8 @@
Helgeson, H. C., Kirkham, D. H. and Flowers, G. C. (1981) Theoretical prediction of the thermodynamic behavior of aqueous electrolytes at high pressures and temperatures. IV. Calculation of activity coefficients, osmotic coefficients, and apparent molal and standard and relative partial molal properties to 600\degC and 5 Kb. \emph{Am. J. Sci.} \bold{281}, 1249--1516. \url{https://doi.org/10.2475/ajs.281.10.1249}
+Manning, C. E. (2013) Thermodynamic modeling of fluid-rock interaction at mid-crustal to upper-mantle conditions. \emph{Rev. Mineral. Geochem.} \bold{76}, 135--164. \url{https://doi.org/10.2138/rmg.2013.76.5}
+
Manning, C. E., Shock, E. L. and Sverjensky, D. A. (2013) The chemistry of carbon in aqueous fluids at crustal and upper-mantle conditions: Experimental and theoretical constraints. \emph{Rev. Mineral. Geochem.} \bold{75}, 109--148. \url{https://doi.org/10.2138/rmg.2013.75.5}
}
Modified: pkg/CHNOSZ/man/water.Rd
===================================================================
--- pkg/CHNOSZ/man/water.Rd 2017-10-12 09:10:59 UTC (rev 253)
+++ pkg/CHNOSZ/man/water.Rd 2017-10-12 14:56:08 UTC (rev 254)
@@ -170,7 +170,7 @@
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/chnosz -r 254
More information about the CHNOSZ-commits
mailing list