[CHNOSZ-commits] r341 - in pkg/CHNOSZ: . R data demo inst man tests/testthat
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Nov 5 10:27:01 CET 2018
Author: jedick
Date: 2018-11-05 10:27:01 +0100 (Mon, 05 Nov 2018)
New Revision: 341
Removed:
pkg/CHNOSZ/demo/oldsolub.R
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/NAMESPACE
pkg/CHNOSZ/R/examples.R
pkg/CHNOSZ/R/nonideal.R
pkg/CHNOSZ/R/subcrt.R
pkg/CHNOSZ/data/opt.csv
pkg/CHNOSZ/demo/00Index
pkg/CHNOSZ/demo/DEW.R
pkg/CHNOSZ/demo/solubility.R
pkg/CHNOSZ/inst/NEWS
pkg/CHNOSZ/man/examples.Rd
pkg/CHNOSZ/man/nonideal.Rd
pkg/CHNOSZ/tests/testthat/test-logmolality.R
pkg/CHNOSZ/tests/testthat/test-nonideal.R
Log:
nonideal(): reorganize options for activity coefficient calculations
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2018-11-05 07:08:12 UTC (rev 340)
+++ pkg/CHNOSZ/DESCRIPTION 2018-11-05 09:27:01 UTC (rev 341)
@@ -1,6 +1,6 @@
Date: 2018-11-05
Package: CHNOSZ
-Version: 1.1.3-48
+Version: 1.1.3-49
Title: Thermodynamic Calculations and Diagrams for Geo(bio)chemistry
Authors at R: c(
person("Jeffrey", "Dick", , "j3ffdick at gmail.com", role = c("aut", "cre"),
Modified: pkg/CHNOSZ/NAMESPACE
===================================================================
--- pkg/CHNOSZ/NAMESPACE 2018-11-05 07:08:12 UTC (rev 340)
+++ pkg/CHNOSZ/NAMESPACE 2018-11-05 09:27:01 UTC (rev 341)
@@ -54,7 +54,7 @@
# added 20170301 or later
"GHS_Tr", "calculateDensity", "calculateGibbsOfWater",
"calculateEpsilon", "calculateQ", "water.DEW", "berman",
- "maxdiff", "expect_maxdiff", "Bdot",
+ "maxdiff", "expect_maxdiff", "bgamma",
# added 20171121 or later
"dumpdata", "thermo.axis", "solubility", "NaCl"
)
Modified: pkg/CHNOSZ/R/examples.R
===================================================================
--- pkg/CHNOSZ/R/examples.R 2018-11-05 07:08:12 UTC (rev 340)
+++ pkg/CHNOSZ/R/examples.R 2018-11-05 09:27:01 UTC (rev 341)
@@ -28,7 +28,7 @@
demos <- function(which=c("sources", "protein.equil", "affinity", "NaCl", "density",
"ORP", "revisit", "findit", "ionize", "buffer", "protbuff", "yeastgfp", "mosaic",
- "copper", "oldsolub", "solubility", "gold", "wjd", "bugstab", "Shh", "activity_ratios",
+ "copper", "solubility", "gold", "wjd", "bugstab", "Shh", "activity_ratios",
"adenine", "DEW", "lambda", "TCA", "go-IU", "bison"), save.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)) {
Modified: pkg/CHNOSZ/R/nonideal.R
===================================================================
--- pkg/CHNOSZ/R/nonideal.R 2018-11-05 07:08:12 UTC (rev 340)
+++ pkg/CHNOSZ/R/nonideal.R 2018-11-05 09:27:01 UTC (rev 341)
@@ -7,26 +7,32 @@
# generate nonideal contributions to thermodynamic properties
# 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/Helgeson0 methods only
+ # arguments A_DH and B_DH are needed for all methods other than "Alberty", and P is needed for "bgamma"
mettext <- function(method) {
- mettext <- paste(method, "method")
- if(method=="Helgeson0") mettext <- "Helgeson method with B-dot = 0"
+ mettext <- paste(method, "equation")
+ if(method=="Bdot0") mettext <- "B-dot equation (B-dot = 0)"
mettext
}
# we can use this function to change the nonideal method option
- if(identical(species, "Helgeson") | identical(species, "Helgeson0") | identical(species, "Alberty")) {
- thermo <- get("thermo")
- oldnon <- thermo$opt$nonideal
- thermo$opt$nonideal <- species
- assign("thermo", thermo, "CHNOSZ")
- message("nonideal: setting nonideal option to use ", mettext(species))
- return(invisible(oldnon))
+ if(missing(speciesprops)) {
+ if(species[1] %in% c("Bdot", "Bdot0", "bgamma", "bgamma0", "Alberty")) {
+ thermo <- get("thermo")
+ oldnon <- thermo$opt$nonideal
+ thermo$opt$nonideal <- species[1]
+ assign("thermo", thermo, "CHNOSZ")
+ message("nonideal: setting nonideal option to use ", mettext(species))
+ return(invisible(oldnon))
+ } else stop(species[1], " is not a valid nonideality setting (Bdot, Bdot0, bgamma, bgamma0, or Alberty)")
}
+ # check if we have a valid method setting
+ thermo <- get("thermo")
+ if(!thermo$opt$nonideal %in% c("Alberty", "Bdot", "Bdot0", "bgamma", "bgamma0")) stop("invalid setting (", thermo$opt$nonideal, ") in thermo$opt$nonideal")
+
# function to calculate extended Debye-Huckel equation and derivatives using Alberty's parameters
- Alberty <- function(Z, I, T, prop = "loggamma") {
+ Alberty <- function(prop = "loggamma", Z, I, T) {
# extended Debye-Huckel equation ("log")
# and its partial derivatives ("G","H","S","Cp")
# T in Kelvin
@@ -54,20 +60,14 @@
}
# function for Debye-Huckel equation with B-dot extended term parameter (Helgeson, 1969)
- Helgeson <- function(Z, I, T, P, A_DH, B_DH, prop = "loggamma", acirc) {
- ## "distance of closest approach" of ions in NaCl solutions (HKF81 Table 2)
- #acirc <- 3.72e-8 # cm
- if(method=="Helgeson") loggamma <- - A_DH * Z^2 * I^0.5 / (1 + acirc * B_DH * I^0.5) + Bdot * I
- else if(method=="Helgeson0") loggamma <- - A_DH * Z^2 * I^0.5 / (1 + acirc * B_DH * I^0.5)
+ Helgeson <- function(prop = "loggamma", Z, I, T, A_DH, B_DH, acirc, bgamma) {
+ loggamma <- - A_DH * Z^2 * I^0.5 / (1 + acirc * B_DH * I^0.5) + bgamma * I
R <- 1.9872 # gas constant, cal K^-1 mol^-1
if(prop=="loggamma") return(loggamma)
else if(prop=="G") return(R * T * log(10) * loggamma)
# note the log(10) (=2.303) ... use natural logarithm to calculate G
}
- # 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")
# loop over species #1: get the charge
@@ -82,7 +82,7 @@
Z[i] <- thisZ
}
# get species formulas to assign acirc 20181105
- if(grepl("Helgeson", method)) {
+ if(grepl("Bdot", method)) {
formula <- get("thermo")$obigt$formula[species]
# "ion size paramter" taken from UT_SIZES.REF of HCh package (Shvarov and Bastrakov, 1999),
# based on Table 2.7 of Garrels and Christ, 1965
@@ -102,9 +102,16 @@
nZ <- sum(Z!=0)
if(nZ > 1) message("nonideal: using ", paste(acirc[Z!=0], collapse=" "), " for ion size parameters of ", paste(formula[Z!=0], collapse=" "))
else if(nZ==1) message("nonideal: using ", acirc[Z!=0], " for ion size parameter of ", formula[Z!=0])
- # use correct units for ion size paramter
+ # use correct units (cm) for ion size parameter
acirc <- acirc * 10^-8
+ } else if(grepl("bgamma", method)) {
+ # "distance of closest approach" of ions in NaCl solutions (HKF81 Table 2)
+ acirc <- rep(3.72e-8, length(species))
}
+ # get b_gamma (or Bdot)
+ if(method=="bgamma") bgamma <- bgamma(convert(T, "C"), P)
+ else if(method=="Bdot") bgamma <- Bdot(convert(T, "C"))
+ else if(method %in% c("Bdot0", "bgamma0")) bgamma <- 0
# loop over species #2: activity coefficient calculations
iH <- info("H+")
ie <- info("e-")
@@ -120,20 +127,19 @@
didit <- FALSE
for(j in 1:ncol(myprops)) {
pname <- colnames(myprops)[j]
- if(method=="Alberty" & pname %in% c("G", "H", "S", "Cp")) {
- myprops[, j] <- myprops[, j] + Alberty(Z[i], IS, T, pname)
+ if(!pname %in% c("G", "H", "S", "Cp")) next
+ if(method=="Alberty") {
+ myprops[, j] <- myprops[, j] + Alberty(pname, Z[i], IS, T)
didit <- TRUE
- } else if(grepl("Helgeson", method) & pname %in% c("G", "H", "S", "Cp")) {
- myprops[, j] <- myprops[, j] + Helgeson(Z[i], IS, T, P, A_DH, B_DH, pname, acirc[i])
+ } else {
+ myprops[, j] <- myprops[, j] + Helgeson(pname, Z[i], IS, T, A_DH, B_DH, acirc[i], bgamma)
didit <- TRUE
}
}
# append a loggam column if we did any nonideal calculations of thermodynamic properties
if(didit) {
- if(method=="Alberty") myprops <- cbind(myprops, loggam = Alberty(Z[i], IS, T))
- else if(grepl("Helgeson", method)) {
- myprops <- cbind(myprops, loggam = Helgeson(Z[i], IS, T, P, A_DH, B_DH, "loggamma", acirc[i]))
- }
+ if(method=="Alberty") myprops <- cbind(myprops, loggam = Alberty("loggamma", Z[i], IS, T))
+ else myprops <- cbind(myprops, loggam = Helgeson("loggamma", Z[i], IS, T, A_DH, B_DH, acirc[i], bgamma))
}
speciesprops[[i]] <- myprops
if(didit) ndid <- ndid + 1
@@ -142,8 +148,8 @@
return(speciesprops)
}
-Bdot <- function(TC = 25, P = 1, showsplines = "") {
- # 20171012 calculate B-dot (bgamma) using P, T, points from:
+bgamma <- function(TC = 25, P = 1, showsplines = "") {
+ # 20171012 calculate b_gamma using P, T, points from:
# Helgeson, 1969 (doi:10.2475/ajs.267.7.729)
# Helgeson et al., 1981 (doi:10.2475/ajs.281.10.1249)
# Manning et al., 2013 (doi:10.2138/rmg.2013.75.5)
@@ -292,7 +298,7 @@
legend=c("Helgeson, 1969", "Helgeson et al., 1981", "Manning et al., 2013", "spline control point", "high-P extrapolation"))
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"))
- title(main=expression("Deybe-H\u00FCckel "*italic(b)[gamma]*" ('B-dot') parameter"))
+ title(main=expression("Deybe-H\u00FCckel extended term ("*italic(b)[gamma]*") parameter"))
} else if(showsplines=="P") {
thermo.plot.new(c(0, 5), c(-.2, .7), xlab=expression(log~italic(P)*"(bar)"), ylab=expression(italic(b)[gamma]))
# pressures that are used to make the isothermal splines (see below)
@@ -308,44 +314,44 @@
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))
- 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))
+ points(log10(P25), bgamma(25, P25))
+ points(log10(P100), bgamma(100, P100))
+ points(log10(P200), bgamma(200, P200))
+ points(log10(P300), bgamma(300, P300))
+ points(log10(P400), bgamma(400, P400))
+ points(log10(P500), bgamma(500, P500))
+ points(log10(P600), bgamma(600, P600))
+ points(log10(P700), bgamma(700, P700))
+ points(log10(P800), bgamma(800, P800))
+ points(log10(P900), bgamma(900, P900))
+ points(log10(P1000), bgamma(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), 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])
+ P <- c(1, seq(50, 5000, 50)); lines(log10(P), bgamma(25, P), col=col[1])
+ P <- c(1, seq(50, 20000, 50)); lines(log10(P), bgamma(100, P), col=col[2])
+ P <- c(1, seq(50, 40000, 50)); lines(log10(P), bgamma(200, P), col=col[3])
+ P <- c(1, seq(50, 60000, 50)); lines(log10(P), bgamma(300, P), col=col[4])
+ P <- seq(500, 60000, 50); lines(log10(P), bgamma(400, P), col=col[5])
+ P <- seq(1000, 60000, 50); lines(log10(P), bgamma(500, P), col=col[6])
+ P <- seq(2000, 60000, 50); lines(log10(P), bgamma(600, P), col=col[7])
+ P <- seq(10000, 60000, 50); lines(log10(P), bgamma(700, P), col=col[8])
+ P <- seq(10000, 60000, 50); lines(log10(P), bgamma(800, P), col=col[9])
+ P <- seq(10000, 60000, 50); lines(log10(P), bgamma(900, P), col=col[10])
+ P <- seq(10000, 60000, 50); lines(log10(P), bgamma(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")
- title(main=expression("Deybe-H\u00FCckel "*italic(b)[gamma]*" ('B-dot') parameter"))
+ title(main=expression("Deybe-H\u00FCckel extended term ("*italic(b)[gamma]*") parameter"))
} else {
# make T and P the same length
ncond <- max(length(T), length(P))
T <- rep(T, length.out=ncond)
P <- rep(P, length.out=ncond)
# loop over P, T conditions
- Bdot <- numeric()
+ bgamma <- numeric()
lastT <- NULL
for(i in 1:length(T)) {
# make it fast: skip splines at 25 degC and 1 bar
- if(T[i]==25 & P[i]==1) Bdot <- c(Bdot, 0.041)
+ if(T[i]==25 & P[i]==1) bgamma <- c(bgamma, 0.041)
else {
if(!identical(T[i], lastT)) {
# get the spline fits from particular pressures for each T
@@ -381,9 +387,17 @@
# 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]))
+ bgamma <- c(bgamma, ST(P[i]))
}
}
- return(Bdot)
+ return(bgamma)
}
}
+
+### unexported functions ###
+
+Bdot <- function(TC) {
+ Bdot <- splinefun(c(25, 50, 100, 150, 200, 250, 300), c(0.0418, 0.0439, 0.0468, 0.0479, 0.0456, 0.0348, 0))(TC)
+ Bdot[TC > 300] <- 0
+ return(Bdot)
+}
Modified: pkg/CHNOSZ/R/subcrt.R
===================================================================
--- pkg/CHNOSZ/R/subcrt.R 2018-11-05 07:08:12 UTC (rev 340)
+++ pkg/CHNOSZ/R/subcrt.R 2018-11-05 09:27:01 UTC (rev 341)
@@ -280,7 +280,7 @@
# always get density
H2O.props <- "rho"
# calculate A_DH and B_DH if we're using the B-dot (Helgeson) equation
- if(any(IS != 0) & grepl("Helgeson", thermo$opt$nonideal)) H2O.props <- c(H2O.props, "A_DH", "B_DH")
+ if(any(IS != 0) & thermo$opt$nonideal %in% c("Bdot", "Bdot0", "bgamma", "bgamma0")) H2O.props <- c(H2O.props, "A_DH", "B_DH")
# get other properties for H2O only if it's in the reaction
if(any(isH2O)) H2O.props <- c(H2O.props, eosprop)
hkfstuff <- hkf(eosprop, parameters = param, T = T, P = P, H2O.props=H2O.props)
@@ -298,7 +298,7 @@
}
# calculate activity coefficients if ionic strength is not zero
if(any(IS != 0)) {
- if(grepl("Helgeson", thermo$opt$nonideal)) p.aq <- nonideal(iphases[isaq], p.aq, newIS, T, P, H2O.PT$A_DH, H2O.PT$B_DH)
+ if(thermo$opt$nonideal %in% c("Bdot", "Bdot0", "bgamma", "bgamma0")) p.aq <- nonideal(iphases[isaq], p.aq, newIS, T, P, H2O.PT$A_DH, H2O.PT$B_DH)
else if(thermo$opt$nonideal=="Alberty") p.aq <- nonideal(iphases[isaq], p.aq, newIS, T)
}
outprops <- c(outprops, p.aq)
Modified: pkg/CHNOSZ/data/opt.csv
===================================================================
--- pkg/CHNOSZ/data/opt.csv 2018-11-05 07:08:12 UTC (rev 340)
+++ pkg/CHNOSZ/data/opt.csv 2018-11-05 09:27:01 UTC (rev 341)
@@ -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,nonideal,Berman,maxcores
-1E-10,cal,C,bar,aq,SUPCRT92,100,1,1,FALSE,liquid,1000,TRUE,TRUE,Helgeson,NA,2
+1E-10,cal,C,bar,aq,SUPCRT92,100,1,1,FALSE,liquid,1000,TRUE,TRUE,Bdot,NA,2
Modified: pkg/CHNOSZ/demo/00Index
===================================================================
--- pkg/CHNOSZ/demo/00Index 2018-11-05 07:08:12 UTC (rev 340)
+++ pkg/CHNOSZ/demo/00Index 2018-11-05 09:27:01 UTC (rev 341)
@@ -12,7 +12,6 @@
yeastgfp Subcellular locations: log fO2 - log aH2O and log a - log fO2 diagrams
mosaic Eh-pH diagram for iron oxides, sulfides and carbonate with two sets of changing basis species
copper Another example of mosaic(): complexation of copper with glycine species
-oldsolub Old solubility calculations using uniroot()
solubility Solubility of calcite and CO2(gas) as a function of pH
gold Solubility of gold
wjd Gibbs energy minimization: prebiological atmospheres and cell periphery of yeast
Modified: pkg/CHNOSZ/demo/DEW.R
===================================================================
--- pkg/CHNOSZ/demo/DEW.R 2018-11-05 07:08:12 UTC (rev 340)
+++ pkg/CHNOSZ/demo/DEW.R 2018-11-05 09:27:01 UTC (rev 341)
@@ -155,8 +155,8 @@
IS <- c(0.39, 0.57, 0.88, 1.45, 2.49)
pH <- c(3.80, 3.99, 4.14, 4.25, 4.33)
molC <- c(0.03, 0.2, 1, 4, 20)
-## use Debye-Huckel equation with B-dot set to zero
-nonideal("Helgeson0")
+## use extended Debye-Huckel equation with b_gamma set to zero
+nonideal("bgamma0")
## calculate affinities on the T-logfO2-pH-IS transect
a <- affinity(T = T, O2 = buf$O2 - 2, IS = IS, pH = pH, P = 50000)
## calculate metastable equilibrium activities using the total
Deleted: pkg/CHNOSZ/demo/oldsolub.R
===================================================================
--- pkg/CHNOSZ/demo/oldsolub.R 2018-11-05 07:08:12 UTC (rev 340)
+++ pkg/CHNOSZ/demo/oldsolub.R 2018-11-05 09:27:01 UTC (rev 341)
@@ -1,87 +0,0 @@
-# CHNOSZ/demo/oldsolub.R
-# 20150306 added to CHNOSZ as solubility.R
-# 20181030 renamed to oldsolub.R
-
-# demo showing how to calculate CO2(gas) or calcite solubility and aqueous carbonate speciation
-# the affinity() ... equilibrate() sequence in CHNOSZ gives *metastable equilibrium activities*
-# (activities for a total activity of the balanced component given in loga.balance) ...
-# here we are interested in finding the value of loga.balance itself
-# (the total activity of [CO2, HCO3-, CO3-2] species in the aqueous phase).
-# this total activity is the solubility of CO2(gas) or calcite if the affinities of the
-# aqueous species as formed from CO2(gas) or calcite are all equal to zero.
-# note that the affinities for species in metastable equilibrium are all equal.
-# Afun() calculates the metastable equilibrium affinities for a given loga.balance
-# and uniroot() finds the loga.balance where they are zero
-# additionally, if we are reacting calcite, the activity of Ca+2 should be set equal to loga.balance
-
-# for comparison with published calcite solubility plot, see Fig. 4A in
-# Manning et al., 2013, Reviews in Mineralogy & Geochemistry, v. 75, pp. 109-148
-# (doi: 10.2138/rmg.2013.75.5)
-
-# for comparison with published CO2 solubility plot, see Fig. 4.5 in
-# Stumm and Morgan, 1996, Aquatic Chemistry: Chemical Equilibria and Rates in Natural Waters
-# (New York: John Wiley & Sons), 3rd edition
-
-# set this to CO2 or calcite
-what <- "calcite"
-#what <- "CO2"
-
-# function to return the affinity of the metastable equilibrium species
-Afun <- function(loga.balance=-3, T=25) {
- if(what=="calcite") basis("Ca+2", loga.balance)
- a <- affinity(T=T)
- e <- equilibrate(a, loga.balance=loga.balance)
- # set metastable activities and re-calculate the affinity
- species(1:3, unlist(e$loga.equil))
- a <- affinity(T=T)
- # check they're actually equal
- stopifnot(all(abs(unlist(a$values) - as.vector(a$values[[1]])) < 1e-10))
- return(a$values[[1]])
-}
-
-# set up system
-if(what=="CO2") {
- basis("CHNOS+")
- basis("CO2", "gas")
- # ca. atmospheric PCO2
- basis("CO2", -3.5)
-} else if(what=="calcite") {
- basis(c("calcite", "Ca+2", "H2O", "O2", "H+"))
-}
-species(c("CO2", "HCO3-", "CO3-2"))
-T <- 25
-# decrease this for higher resolution
-pHstep <- 1
-
-# where we'll store the results
-loga.tot <- numeric()
-loga.CO2 <- loga.HCO3 <- loga.CO3 <- numeric()
-
-# loop over pH range
-pHs <- seq(0, 14, pHstep)
-for(pH in pHs) {
- print(paste("pH =", pH))
- basis("pH", pH)
- # this is for the solubility
- loga.balance <- suppressMessages(uniroot(Afun, c(-10, 10), T=T)$root)
- loga.tot <- c(loga.tot, loga.balance)
- # this is for the speciation
- if(what=="calcite") basis("Ca+2", loga.balance)
- a <- affinity(T=T)
- e <- equilibrate(a, loga.balance=loga.balance)
- loga <- unlist(e$loga.equil)
- loga.CO2 <- c(loga.CO2, loga[1])
- loga.HCO3 <- c(loga.HCO3, loga[2])
- loga.CO3 <- c(loga.CO3, loga[3])
-}
-
-# make plot
-ylim <- c(-10, 4)
-thermo.plot.new(xlim=range(pHs), ylim=ylim, xlab="pH", ylab="log a")
-lines(pHs, loga.tot, lwd=4, col="green2")
-lines(pHs, loga.CO2, lwd=2)
-lines(pHs, loga.HCO3, lty=2, lwd=2)
-lines(pHs, loga.CO3, lty=3, lwd=2)
-legend(ifelse(what=="calcite", "topright", "topleft"), lty=c(1, 1:3), lwd=c(4, 2, 2, 2), col=c("green2", rep("black", 3)),
- legend=as.expression(c("total", expr.species("CO2", state="aq"), expr.species("HCO3-"), expr.species("CO3-2"))))
-title(main=substitute("Solubility of"~what~"at"~T~degree*"C", list(what=what, T=T)))
Modified: pkg/CHNOSZ/demo/solubility.R
===================================================================
--- pkg/CHNOSZ/demo/solubility.R 2018-11-05 07:08:12 UTC (rev 340)
+++ pkg/CHNOSZ/demo/solubility.R 2018-11-05 09:27:01 UTC (rev 341)
@@ -1,6 +1,6 @@
-# CHNOSZ/demo/solubility.R: vectorized solubility without uniroot
-# 20181030 adapted from CHNOSZ/demo/oldsolub.R
-# 20181031 use new solubility(); add T-pH plots
+# CHNOSZ/demo/solubility.R: solubility of CO2 and calcite
+# 20150306 jmd first version; used uniroot() to find zero affinity
+# 20181031 use new vectorized, non-uniroot solubility(); add T-pH plots
# for comparison with published CO2 solubility plot, see Fig. 4.5 in
# Stumm and Morgan, 1996, Aquatic Chemistry: Chemical Equilibria and Rates in Natural Waters
Modified: pkg/CHNOSZ/inst/NEWS
===================================================================
--- pkg/CHNOSZ/inst/NEWS 2018-11-05 07:08:12 UTC (rev 340)
+++ pkg/CHNOSZ/inst/NEWS 2018-11-05 09:27:01 UTC (rev 341)
@@ -1,30 +1,35 @@
-CHANGES IN CHNOSZ 1.1.3-48 (2018-11-05)
+CHANGES IN CHNOSZ 1.1.3-49 (2018-11-05)
---------------------------------------
NEW FEATURES
-- Add solubility(). Run this after equilibrate() to calculate the
- solubility (loga.balance) of the balanced basis species.
+- Reorganize and expand options for activity coefficient calculations
+ (set in thermo$opt$nonideal: Bdot, Bdot0, bgamma, bgamma0, or Alberty).
+ The previous default, which corresponds to 'bgamma' (T- and
+ P-dependent extended term parameter with single ion-size parameter),
+ has been replaced by 'Bdot' (T-dependent extended term parameter and
+ species-dependent ion-size parameter; see below).
+- nonideal() with the 'Bdot' or 'Bdot0' equation uses specific
+ ion-size parameters for different ions, in accord with the HCh package
+ (Shvarov and Bastrakov, 1999). Parameters are from Table 2.7 of
+ Garrels and Christ, 1965.
+
- Add NaCl(), implementing a first-order calculation of the speciation
of NaCl in water, taking account of activity coefficients and the
reaction Na+ + Cl- = NaCl(aq).
-- nonideal() now uses specific ion-size parameters for different ions,
- in accord with the HCh package (Shvarov and Bastrakov, 1999).
- Parameters are from Table 2.7 of Garrels and Christ, 1965.
+- Add solubility(). Run this after equilibrate() to calculate the
+ solubility (loga.balance) of the balanced basis species.
+- Revise demo/solubility.R to show solubility calculations for CO2(gas)
+ and calcite as a function of T and pH.
+
- Add demo/gold.R for calculations of Au solubility in hydrothermal
chloride and sulfide solutions (based on diagrams from Akinfiev and
Zotov, 2001, Stefánsson and Seward, 2004, and Williams-Jones et al.,
2009).
-- Revise demo/solubility.R to show solubility calculations for CO2(gas)
- and calcite as a function of T and pH.
-
-- The old solubility demo, which uses uniroot() instead of the
- vectorized calculations in solubility(), has been renamed oldsolub.R.
-
THERMODYNAMIC DATA
- The Berman data (Berman, 1988 and later additions) have replaced the
Modified: pkg/CHNOSZ/man/examples.Rd
===================================================================
--- pkg/CHNOSZ/man/examples.Rd 2018-11-05 07:08:12 UTC (rev 340)
+++ pkg/CHNOSZ/man/examples.Rd 2018-11-05 09:27:01 UTC (rev 341)
@@ -15,9 +15,9 @@
examples(save.png = FALSE)
demos(which = c("sources", "protein.equil", "affinity", "NaCl",
"density", "ORP", "revisit", "findit", "ionize", "buffer",
- "protbuff", "yeastgfp", "mosaic", "copper", "oldsolub",
- "solubility", "gold", "wjd", "bugstab", "Shh", "activity_ratios",
- "adenine", "DEW", "lambda", "TCA", "go-IU", "bison"),
+ "protbuff", "yeastgfp", "mosaic", "copper", "solubility",
+ "gold", "wjd", "bugstab", "Shh", "activity_ratios", "adenine",
+ "DEW", "lambda", "TCA", "go-IU", "bison"),
save.png=FALSE)
}
@@ -45,7 +45,6 @@
\code{yeastgfp} \tab * Subcellular locations: \logfO2 - \logaH2O and \loga - \logfO2 diagrams (Dick, 2009) \cr
\code{mosaic} \tab * Eh-pH diagram with two sets of changing basis species (Garrels and Christ, 1965) \cr
\code{copper} \tab * Another example of \code{\link{mosaic}}: complexation of Cu with glycine (Aksu and Doyle, 2001) \cr
- \code{oldsolub} \tab Old solubility calculations using \code{\link{uniroot}} \cr
\code{solubility} \tab * Solubility of calcite (cf. Manning et al., 2013) and \CO2 (cf. Stumm and Morgan, 1996) \cr
\code{gold} \tab * Solubility of gold (Akinfiev and Zotov; 2001; Stef{\aacute}nsson and Seward, 2004; Williams-Jones et al., 2009) \cr
\code{wjd} \tab * \eqn{G}{G} minimization: prebiological atmospheres (Dayhoff et al., 1964) and cell periphery of yeast \cr
Modified: pkg/CHNOSZ/man/nonideal.Rd
===================================================================
--- pkg/CHNOSZ/man/nonideal.Rd 2018-11-05 07:08:12 UTC (rev 340)
+++ pkg/CHNOSZ/man/nonideal.Rd 2018-11-05 09:27:01 UTC (rev 341)
@@ -1,7 +1,7 @@
\encoding{UTF-8}
\name{nonideal}
\alias{nonideal}
-\alias{Bdot}
+\alias{bgamma}
\title{Activity coefficients of aqueous species}
\description{
Calculate activity coefficients and adjusted (non-ideal) molal properties of aqueous species.
@@ -10,7 +10,7 @@
\usage{
nonideal(species, speciesprops, IS, T, P, A_DH, B_DH,
method=get("thermo")$opt$nonideal)
- Bdot(TC, P, showsplines = "")
+ bgamma(TC, P, showsplines = "")
}
\arguments{
@@ -18,34 +18,39 @@
\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}, \samp{Helgeson}, or \samp{Helgeson0}}
+ \item{P}{numeric, pressure (bar); required for B-dot or b_gamma equation}
+ \item{A_DH}{numeric, A Debye-Huckel coefficient; required for B-dot or b_gamma equation}
+ \item{B_DH}{numeric, B Debye-Huckel coefficient; required for B-dot or b_gamma equation}
+ \item{method}{character, \samp{Alberty}, \samp{Bdot}, \samp{Bdot0}, or \samp{bgamma}}
\item{TC}{numeric, temperature (\degC)}
\item{showsplines}{character, show isobaric (\samp{T}) or isothermal (\samp{P}) splines}
}
\details{
\code{nonideal} takes a list of dataframes (in \code{speciesprops}) containing the standard molal properties of the identified \code{species}.
-The function calculates the *adjusted* properties for given ionic strength (\code{IS}); they are equal to the *standard* values at IS=0.
+The function calculates the *adjusted* properties for given ionic strength (\code{IS}); they are equal to the *standard* values only 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 (\Hplus) and electron (\eminus) are also bypassed by default; this makes sense if you are setting the pH, i.e. activity of \Hplus, to some value.
To apply the calculations to H+ and/or e-, change \code{\link{thermo}$opt$ideal.H} or \code{ideal.e} to FALSE.
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.
-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 adjusted molal properties at the specified ionic strength(s) and temperature(s).
+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; Hückel, 1925) and its derivatives, to calculate adjusted molal properties at the specified ionic strength(s) and temperature(s).
The adjusted 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 is used.
-This equation seems to have been originally proposed by Huckel, 1925; parameters were derived for use at high temperature and pressure by Helgeson, 1969; Helgeson et al., 1981; Manning, 2013.
+In addition to \code{IS} and \code{T}, the following two methods depend on values of \code{P}, \code{A_DH}, and \code{B_DH} given in the arguments.
+For these methods, \samp{G} is currently the only adjusted molal property that is calculated (but this can be used by \code{\link{subcrt}} to calculate adjusted equilibrium constants).
+
+If \code{method} is \samp{Bdot} (the default setting in CHNOSZ), the \dQuote{B-dot} form of the extended Debye-Hückel equation equation is used.
The distance of closest approach (the \dQuote{ion size parameter}) is taken from the UT_SIZES.REF file of the HCh package (Shvarov and Bastrakov, 1992), which is based on Table 2.7 of Garrels and Christ, 1965.
-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.
-For some uses, it is desirable to set the \dQuote{B-dot} parameter to zero; this can be done by setting the method to \code{Helgeson0}.
-Currently, \samp{G} is the only adjusted molal property that is calculated (but this can be used by \code{\link{subcrt}} to calculate adjusted equilibrium constants).
+The extended term parameter is expressed as \dQuote{B-dot}, which is a function only of temperature (Helgeson, 1969).
+For some uses, it is desirable to set the extended term parameter to zero; this can be done by setting the method to \samp{Bdot0}.
-\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.
+If \code{method} is \samp{bgamma}, the \dQuote{b_gamma} equation is used.
+The distance of closest approach is set to a constant (3.72e-8 cm) (e.g., Manning et al., 2013).
+The extended term parameter is calculated by calling the \code{bgamma} function.
+Alternatively, set the extended term parameter to zero with \samp{bgamma0}.
+
+\code{bgamma} calculates the extended term parameter (written as b_gamma; Helgeson et al., 1981) for activity coefficients in NaCl-dominated 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.
@@ -63,7 +68,7 @@
### Examples following Alberty, 2003
### (page numbers given below)
-## the default method setting is Helgeson;
+## the default method setting is Bdot;
## change it to Alberty
oldnon <- nonideal("Alberty")
@@ -132,10 +137,8 @@
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")
+### finished with Alberty equation, let's look at b_gamma
+nonideal("bgamma")
## activity coefficients for monovalent ions at 700 degC, 10 kbar
# after Manning, 2013, Fig. 7
@@ -155,12 +158,15 @@
gamma <- 10^nonidealprops[[1]]$loggam
# the error is larger at higher IS
stopifnot(maxdiff(gamma[1], Manning_gamma[1]) < 0.01)
-stopifnot(maxdiff(gamma, Manning_gamma) < 0.36)
+stopifnot(maxdiff(gamma, Manning_gamma) < 0.23)
-## data and splines used for calculating B-dot
+## data and splines used for calculating b_gamma
## (extended term parameter)
-Bdot(showsplines = "T")
-Bdot(showsplines = "P")
+bgamma(showsplines = "T")
+bgamma(showsplines = "P")
+
+# done with examples, restore default setting
+nonideal(oldnon) # == nonideal("Bdot")
}
\references{
Modified: pkg/CHNOSZ/tests/testthat/test-logmolality.R
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/chnosz -r 341
More information about the CHNOSZ-commits
mailing list