[CHNOSZ-commits] r342 - in pkg/CHNOSZ: . R demo man tests/testthat
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Nov 5 12:32:04 CET 2018
Author: jedick
Date: 2018-11-05 12:32:04 +0100 (Mon, 05 Nov 2018)
New Revision: 342
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/R/NaCl.R
pkg/CHNOSZ/R/nonideal.R
pkg/CHNOSZ/demo/DEW.R
pkg/CHNOSZ/demo/gold.R
pkg/CHNOSZ/man/NaCl.Rd
pkg/CHNOSZ/man/nonideal.Rd
pkg/CHNOSZ/tests/testthat/test-logmolality.R
pkg/CHNOSZ/tests/testthat/test-nonideal.R
Log:
NaCl(): also consider activity coefficient of Na+
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2018-11-05 09:27:01 UTC (rev 341)
+++ pkg/CHNOSZ/DESCRIPTION 2018-11-05 11:32:04 UTC (rev 342)
@@ -1,6 +1,6 @@
Date: 2018-11-05
Package: CHNOSZ
-Version: 1.1.3-49
+Version: 1.1.3-50
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/R/NaCl.R
===================================================================
--- pkg/CHNOSZ/R/NaCl.R 2018-11-05 09:27:01 UTC (rev 341)
+++ pkg/CHNOSZ/R/NaCl.R 2018-11-05 11:32:04 UTC (rev 342)
@@ -3,17 +3,18 @@
# given a total molality of NaCl
# taking account of ion association: Na+ + Cl- = NaCl(aq)
# 20181102 jmd first version
+# 20181105 add activity coefficients of Na+
NaCl <- function(T=seq(100, 500, 100), P=1000, m_tot=2) {
# define a function for the reaction quotient
- logQ <- function(m_Cl, gamma) {
+ logQ <- function(m_Cl, gam_NaCl, gam_Na, gam_Cl) {
# starting with Q = a_NaCl / (a_Na+ * a_Cl-),
- # substitute gam_NaCl = 0, m_NaCl + m_Cl = m_tot, m_Cl = m_Na, gam_Cl = gam_Na = gamma
+ # substitute m_tot = m_NaCl + m_Cl and m_Cl = m_Na
# to write:
- log10( (m_tot - m_Cl) / (m_Cl * gamma) ^ 2 )
+ log10( (m_tot - m_Cl) * gam_NaCl / (m_Cl ^ 2 * gam_Na * gam_Cl) )
}
# define a function for affinity = log(K / Q)
- A <- function(m_Cl, gamma, logK) logK - logQ(m_Cl, gamma)
+ A <- function(m_Cl, gam_NaCl, gam_Na, gam_Cl, logK) logK - logQ(m_Cl, gam_NaCl, gam_Na, gam_Cl)
# calculate equilibrium constant at all temperatures (standard conditions: IS = 0)
logK <- subcrt(c("Na+", "Cl-", "NaCl"), c(-1, -1, 1), T = T, P = P)$out$logK
# calculate Debye-Huckel parameters at all temperatures
@@ -23,15 +24,18 @@
ISout <- a_Cl <- numeric(N)
# initial guess for m_Cl and ionic strength assuming complete dissociation of NaCl
IS <- m_Cl <- rep(m_tot, N)
- # the species index for Cl-
- iCl <- info("Cl-")
+ # the species indices for Cl- and Na+
+ ispecies <- info(c("Na+", "Cl-"))
# we start by doing calculations for all temperatures
doit <- !logical(N)
while(any(doit)) {
# calculate activity coefficient at given ionic strength
- gamma <- suppressMessages(10^nonideal(iCl, list(data.frame(G=numeric(N))), IS=IS, T=convert(T, "K"), P=P, A_DH=wout$A_DH, B_DH=wout$B_DH)[[1]]$loggam)
+ speciesprops <- rep(list(data.frame(G=numeric(N))), length(ispecies))
+ gammas <- suppressMessages(nonideal(ispecies, speciesprops, IS=IS, T=convert(T, "K"), P=P, A_DH=wout$A_DH, B_DH=wout$B_DH))
+ gam_Na <- 10^gammas[[1]]$loggam
+ gam_Cl <- 10^gammas[[2]]$loggam
# solve for m_Cl
- for(i in which(doit)) m_Cl[i] <- uniroot(A, c(0, m_tot), gamma=gamma[i], logK=logK[i])$root
+ for(i in which(doit)) m_Cl[i] <- uniroot(A, c(0, m_tot), gam_NaCl=1, gam_Na=gam_Na[i], gam_Cl=gam_Cl[i], logK=logK[i])$root
# calculate new ionic strength and deviation
ISnew <- m_Cl
dIS <- ISnew - IS
@@ -40,9 +44,10 @@
# keep going until the deviation in ionic strength at any temperature is less than 0.01
doit <- abs(dIS) > 0.01
}
- # assemble final gamma and activity of Cl-
- gamma <- suppressMessages(10^nonideal(iCl, list(data.frame(G=numeric(N))), IS=IS, T=convert(T, "K"), P=P, A_DH=wout$A_DH, B_DH=wout$B_DH)[[1]]$loggam)
- a_Cl <- m_Cl * gamma
+ # assemble final molality of Cl- and gammas
+ gammas <- suppressMessages(nonideal(ispecies, speciesprops, IS=IS, T=convert(T, "K"), P=P, A_DH=wout$A_DH, B_DH=wout$B_DH))
+ gam_Na <- 10^gammas[[1]]$loggam
+ gam_Cl <- 10^gammas[[2]]$loggam
# return the calculated values
- list(IS=IS, gamma=gamma, a_Cl=a_Cl)
+ list(IS=IS, m_Cl=m_Cl, gam_Na=gam_Na, gam_Cl=gam_Cl)
}
Modified: pkg/CHNOSZ/R/nonideal.R
===================================================================
--- pkg/CHNOSZ/R/nonideal.R 2018-11-05 09:27:01 UTC (rev 341)
+++ pkg/CHNOSZ/R/nonideal.R 2018-11-05 11:32:04 UTC (rev 342)
@@ -3,11 +3,12 @@
# moved to nonideal.R from util.misc.R 20151107
# added Helgeson method 20171012
-nonideal <- function(species, speciesprops, IS, T, P, A_DH, B_DH, method=get("thermo")$opt$nonideal) {
+nonideal <- function(species, speciesprops, IS, T, P, A_DH, B_DH, m_star=NULL, method=get("thermo")$opt$nonideal) {
# generate nonideal contributions to thermodynamic properties
# number of species, same length as speciesprops list
# T in Kelvin, same length as nrows of speciespropss
# arguments A_DH and B_DH are needed for all methods other than "Alberty", and P is needed for "bgamma"
+ # m_start is the total molality of all dissolved species; if not given, it is taken to be equal to ionic strength
mettext <- function(method) {
mettext <- paste(method, "equation")
@@ -28,8 +29,10 @@
}
# 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")
+ if(!method %in% c("Alberty", "Bdot", "Bdot0", "bgamma", "bgamma0")) {
+ if(missing(method)) stop("invalid setting (", thermo$opt$nonideal, ") in thermo$opt$nonideal")
+ else stop("invalid method (", thermo$opt$nonideal, ")")
+ }
# function to calculate extended Debye-Huckel equation and derivatives using Alberty's parameters
Alberty <- function(prop = "loggamma", Z, I, T) {
@@ -60,8 +63,8 @@
}
# function for Debye-Huckel equation with B-dot extended term parameter (Helgeson, 1969)
- 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
+ Helgeson <- function(prop = "loggamma", Z, I, T, A_DH, B_DH, acirc, m_star, bgamma) {
+ loggamma <- - A_DH * Z^2 * I^0.5 / (1 + acirc * B_DH * I^0.5) - log10(1 + 0.0180153 * m_star) + 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)
@@ -113,6 +116,7 @@
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
+ if(is.null(m_star)) m_star <- IS
iH <- info("H+")
ie <- info("e-")
speciesprops <- as.list(speciesprops)
@@ -132,14 +136,14 @@
myprops[, j] <- myprops[, j] + Alberty(pname, Z[i], IS, T)
didit <- TRUE
} else {
- myprops[, j] <- myprops[, j] + Helgeson(pname, Z[i], IS, T, A_DH, B_DH, acirc[i], bgamma)
+ myprops[, j] <- myprops[, j] + Helgeson(pname, Z[i], IS, T, A_DH, B_DH, acirc[i], m_star, 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("loggamma", Z[i], IS, T))
- else myprops <- cbind(myprops, loggam = Helgeson("loggamma", Z[i], IS, T, A_DH, B_DH, acirc[i], bgamma))
+ else myprops <- cbind(myprops, loggam = Helgeson("loggamma", Z[i], IS, T, A_DH, B_DH, acirc[i], m_star, bgamma))
}
speciesprops[[i]] <- myprops
if(didit) ndid <- ndid + 1
Modified: pkg/CHNOSZ/demo/DEW.R
===================================================================
--- pkg/CHNOSZ/demo/DEW.R 2018-11-05 09:27:01 UTC (rev 341)
+++ pkg/CHNOSZ/demo/DEW.R 2018-11-05 11:32:04 UTC (rev 342)
@@ -203,7 +203,9 @@
loggamma <- c(-0.15, -0.18, -0.22, -0.26, -0.31)
# activity coefficients calculated in CHNOSZ
sres <- subcrt("propanoate", T = seq(600, 1000, 100), P = 50000, IS = c(0.39, 0.57, 0.88, 1.45, 2.49))
-stopifnot(maxdiff(sres$out[[1]]$loggam, loggamma) < 0.004)
+stopifnot(maxdiff(sres$out[[1]]$loggam, loggamma) < 0.023)
+# if m_star in nonideal() was zero, we could decrease the tolerance here
+#stopifnot(maxdiff(sres$out[[1]]$loggam, loggamma) < 0.004)
###########
### all done!
Modified: pkg/CHNOSZ/demo/gold.R
===================================================================
--- pkg/CHNOSZ/demo/gold.R 2018-11-05 09:27:01 UTC (rev 341)
+++ pkg/CHNOSZ/demo/gold.R 2018-11-05 11:32:04 UTC (rev 342)
@@ -134,8 +134,9 @@
basis("H+", "QMK")
# calculate solution composition for 2 mol/kg NaCl
NaCl <- NaCl(T = seq(150, 550, 10), P = 1000, m_tot=2)
+ a_Cl <- NaCl$m_Cl * NaCl$gam_Cl
# calculate affinity, equilibrate, solubility
- a <- affinity(T = seq(150, 550, 10), `Cl-` = log10(NaCl$a_Cl), P = 1000, IS = NaCl$IS)
+ a <- affinity(T = seq(150, 550, 10), `Cl-` = log10(a_Cl), P = 1000, IS = NaCl$IS)
e <- equilibrate(a)
s <- solubility(e)
# make diagram and show total log molality
@@ -164,8 +165,9 @@
basis("H+", "QMK")
# calculate solution composition for 2 mol/kg NaCl
NaCl <- NaCl(T = seq(150, 550, 10), P = 1000, m_tot=2)
+ a_Cl <- NaCl$m_Cl * NaCl$gam_Cl
# calculate affinity, equilibrate, solubility
- a <- affinity(T = seq(150, 550, 10), `Cl-` = log10(NaCl$a_Cl), P = 1000, IS = NaCl$IS)
+ a <- affinity(T = seq(150, 550, 10), `Cl-` = log10(a_Cl), P = 1000, IS = NaCl$IS)
e <- equilibrate(a)
s <- solubility(e)
# make diagram and show total log molality
Modified: pkg/CHNOSZ/man/NaCl.Rd
===================================================================
--- pkg/CHNOSZ/man/NaCl.Rd 2018-11-05 09:27:01 UTC (rev 341)
+++ pkg/CHNOSZ/man/NaCl.Rd 2018-11-05 11:32:04 UTC (rev 342)
@@ -23,7 +23,7 @@
The algorithm starts by calculating the equilibrium constant (\emph{K}) of the reaction and assuming complete dissociation of NaCl(aq).
This also permits calculating the ionic strength from the molalities of Na\S{+} and Cl\S{-}.
Then, \code{\link{uniroot}} is used to find the equilibrium molality of Cl\S{-}; that is, where the affinity of the reaction (log(\emph{K}/\emph{Q})) becomes zero.
-The activity quotient (\emph{Q}) is evaluated taking account of activity coefficients calculated for the nominal ionic strength (see \code{\link{nonideal}}).
+The activity quotient (\emph{Q}) is evaluated taking account of activity coefficients of Na\S{+} and Cl\S{-} calculated for the nominal ionic strength (see \code{\link{nonideal}}).
The calculated molality of Cl\S{-} yields a new estimate of the ionic strength of the system.
The calculations are iterated until the deviation in ionic strength at all temperatures is less than 0.01.
}
@@ -34,7 +34,7 @@
}
\value{
-A list with components \samp{IS} (\dQuote{true} ionic strength from concentrations of unpaired ions), \samp{gamma} (activity coefficient of Cl\S{-}), \samp{a_Cl} (activity of Cl\S{-}).
+A list with components \samp{IS} (\dQuote{true} ionic strength from concentrations of unpaired ions), \samp{m_Cl} (molality of Cl\S{-}), \samp{gam_Na}, and \samp{gam_Cl} (activity coefficients of Na\S{+} and Cl\S{-}).
}
\seealso{
@@ -49,20 +49,20 @@
IS.HCh <- list(`100`=c(0.992, 1.969, 2.926, 3.858, 4.758, 5.619),
`300`=c(0.807, 1.499, 2.136, 2.739, 3.317, 3.875),
`500`=c(0.311, 0.590, 0.861, 1.125, 1.385, 1.642))
-gam.HCh <- list(`100`=c(0.565, 0.545, 0.551, 0.567, 0.589, 0.615),
- `300`=c(0.366, 0.307, 0.275, 0.254, 0.238, 0.224),
- `500`=c(0.19, 0.137, 0.111, 0.096, 0.085, 0.077))
+gam_Cl.HCh <- list(`100`=c(0.565, 0.545, 0.551, 0.567, 0.589, 0.615),
+ `300`=c(0.366, 0.307, 0.275, 0.254, 0.238, 0.224),
+ `500`=c(0.19, 0.137, 0.111, 0.096, 0.085, 0.077))
# total molality in the calculation with NaCl()
m_tot <- seq(1, 6, 0.5)
N <- length(m_tot)
# where we'll put the calculated values
IS.calc <- data.frame(`100`=numeric(N), `300`=numeric(N), `500`=numeric(N))
-gam.calc <- data.frame(`100`=numeric(N), `300`=numeric(N), `500`=numeric(N))
+gam_Cl.calc <- data.frame(`100`=numeric(N), `300`=numeric(N), `500`=numeric(N))
# NaCl() is *not* vectorized over m_tot, so we use a loop here
for(i in 1:length(m_tot)) {
NaCl.out <- NaCl(c(100, 300, 500), P=1000, m_tot=m_tot[i])
IS.calc[i, ] <- NaCl.out$IS
- gam.calc[i, ] <- NaCl.out$gamma
+ gam_Cl.calc[i, ] <- NaCl.out$gam_Cl
}
# plot ionic strength from HCh and NaCl() as points and lines
par(mfrow=c(2, 1))
@@ -78,11 +78,13 @@
legend("topleft", dprop, lty=1, pch=1, col=col)
title(main="H2O + NaCl; HCh (points) and 'NaCl()' (lines)")
# plot activity coefficient (gamma)
-plot(c(1,6), c(0,1), xlab="NaCl (mol/kg)", ylab="gamma", type="n")
+plot(c(1,6), c(0,1), xlab="NaCl (mol/kg)", ylab=expression(gamma[Cl^"-"]), type="n")
for(i in 1:3) {
- points(m.HCh, gam.HCh[[i]], col=col[i])
- lines(m_tot, gam.calc[, i], col=col[i])
+ points(m.HCh, gam_Cl.HCh[[i]], col=col[i])
+ lines(m_tot, gam_Cl.calc[, i], col=col[i])
}
+# we should be fairly close
+stopifnot(maxdiff(unlist(gam_Cl.calc[seq(1,11,2), ]), unlist(gam_Cl.HCh)) < 0.034)
}
\references{
Modified: pkg/CHNOSZ/man/nonideal.Rd
===================================================================
--- pkg/CHNOSZ/man/nonideal.Rd 2018-11-05 09:27:01 UTC (rev 341)
+++ pkg/CHNOSZ/man/nonideal.Rd 2018-11-05 11:32:04 UTC (rev 342)
@@ -9,7 +9,7 @@
\usage{
nonideal(species, speciesprops, IS, T, P, A_DH, B_DH,
- method=get("thermo")$opt$nonideal)
+ m_star=NULL, method=get("thermo")$opt$nonideal)
bgamma(TC, P, showsplines = "")
}
@@ -21,6 +21,7 @@
\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{m_star}{numeric, total molality of all dissolved species}
\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}
@@ -37,7 +38,9 @@
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.
-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.
+In addition to \code{IS} and \code{T}, the following two methods depend on values of \code{P}, \code{A_DH}, \code{B_DH}, and \code{m_star} given in the arguments.
+\code{m_star}, the total molality of all dissolved species, is used to compute the mole fraction to molality conversion factor.
+If \code{m_star} is NULL, it is taken to be equal to \code{IS}, which is an underestimate.
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.
Modified: pkg/CHNOSZ/tests/testthat/test-logmolality.R
===================================================================
--- pkg/CHNOSZ/tests/testthat/test-logmolality.R 2018-11-05 09:27:01 UTC (rev 341)
+++ pkg/CHNOSZ/tests/testthat/test-logmolality.R 2018-11-05 11:32:04 UTC (rev 342)
@@ -9,11 +9,10 @@
# the long way...
wprop <- water(c("A_DH", "B_DH"), P=1)
nonid <- nonideal(c("H+", "HCO3-"), subcrt(c("H+", "HCO3-"), T=25)$out, IS=1, T=298.15, P=1, A_DH=wprop$A_DH, B_DH=wprop$B_DH)
- # here we have a small difference due to rounding of the expected value:
- expect_maxdiff(nonid[[2]]$loggam, -0.1790625, 1e-7)
+ # compare with a precalculated value:
+ expect_maxdiff(nonid[[2]]$loggam, -0.1868168, 1e-7)
# the short way...
loggam <- subcrt(c("H+", "HCO3-"), T=25, IS=1)$out[[2]]$loggam
- # we get that loggam(H+)=0 and loggam(HCO3-)=-0.179
expect_equal(nonid[[2]]$loggam, loggam)
## take-home message -1: with default settings, the activity coefficient of H+ is always 1
Modified: pkg/CHNOSZ/tests/testthat/test-nonideal.R
===================================================================
--- pkg/CHNOSZ/tests/testthat/test-nonideal.R 2018-11-05 09:27:01 UTC (rev 341)
+++ pkg/CHNOSZ/tests/testthat/test-nonideal.R 2018-11-05 11:32:04 UTC (rev 342)
@@ -96,14 +96,14 @@
gamCl.300 <- 10^subcrt("Cl-", T=300, P=1000, IS=IS.HCh$`300`)$out$`Cl-`$loggam
gamCl.500 <- 10^subcrt("Cl-", T=500, P=1000, IS=IS.HCh$`500`)$out$`Cl-`$loggam
# TODO: get lower differences by adjusting the activity coefficient model in CHNOSZ
- expect_maxdiff(gamCl.100, gamCl.HCh$`100`, 0.13)
- expect_maxdiff(gamCl.300, gamCl.HCh$`300`, 0.05)
- expect_maxdiff(gamCl.500, gamCl.HCh$`500`, 0.02)
+ expect_maxdiff(gamCl.100, gamCl.HCh$`100`, 0.07)
+ expect_maxdiff(gamCl.300, gamCl.HCh$`300`, 0.03)
+ expect_maxdiff(gamCl.500, gamCl.HCh$`500`, 0.009)
# calculate activity coefficent of Cl- at each temperature
gamNa.100 <- 10^subcrt("Na+", T=100, P=1000, IS=IS.HCh$`100`)$out$`Na+`$loggam
gamNa.300 <- 10^subcrt("Na+", T=300, P=1000, IS=IS.HCh$`300`)$out$`Na+`$loggam
gamNa.500 <- 10^subcrt("Na+", T=500, P=1000, IS=IS.HCh$`500`)$out$`Na+`$loggam
- expect_maxdiff(gamNa.100, gamNa.HCh$`100`, 0.16)
- expect_maxdiff(gamNa.300, gamNa.HCh$`300`, 0.06)
- expect_maxdiff(gamNa.500, gamNa.HCh$`500`, 0.02)
+ expect_maxdiff(gamNa.100, gamNa.HCh$`100`, 0.08)
+ expect_maxdiff(gamNa.300, gamNa.HCh$`300`, 0.03)
+ expect_maxdiff(gamNa.500, gamNa.HCh$`500`, 0.013)
})
More information about the CHNOSZ-commits
mailing list