[CHNOSZ-commits] r346 - in pkg/CHNOSZ: . R inst man man/macros tests/testthat
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Nov 7 16:50:34 CET 2018
Author: jedick
Date: 2018-11-07 16:50:34 +0100 (Wed, 07 Nov 2018)
New Revision: 346
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/R/nonideal.R
pkg/CHNOSZ/R/subcrt.R
pkg/CHNOSZ/R/util.affinity.R
pkg/CHNOSZ/inst/NEWS
pkg/CHNOSZ/man/macros/macros.Rd
pkg/CHNOSZ/man/nonideal.Rd
pkg/CHNOSZ/man/subcrt.Rd
pkg/CHNOSZ/tests/testthat/test-logmolality.R
Log:
nonideal(): add 'is.basis' argument to handle transformed Gibbs energies for basis species
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2018-11-07 13:14:21 UTC (rev 345)
+++ pkg/CHNOSZ/DESCRIPTION 2018-11-07 15:50:34 UTC (rev 346)
@@ -1,6 +1,6 @@
-Date: 2018-11-06
+Date: 2018-11-07
Package: CHNOSZ
-Version: 1.1.3-53
+Version: 1.1.3-54
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/nonideal.R
===================================================================
--- pkg/CHNOSZ/R/nonideal.R 2018-11-07 13:14:21 UTC (rev 345)
+++ pkg/CHNOSZ/R/nonideal.R 2018-11-07 15:50:34 UTC (rev 346)
@@ -3,7 +3,7 @@
# moved to nonideal.R from util.misc.R 20151107
# added Helgeson method 20171012
-nonideal <- function(species, speciesprops, IS, T, P, A_DH, B_DH, m_star=NULL, method=get("thermo")$opt$nonideal) {
+nonideal <- function(species, speciesprops, IS, T, P, A_DH, B_DH, m_star=NULL, method=get("thermo")$opt$nonideal, is.basis=FALSE) {
# generate nonideal contributions to thermodynamic properties
# number of species, same length as speciesprops list
# T in Kelvin, same length as nrows of speciespropss
@@ -123,6 +123,9 @@
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
+ # different signs for basis species and species of interest 20181107
+ species.sign <- ifelse(is.basis, -1, 1)
+ species.sign <- rep(species.sign, length.out=length(species))
# loop over species #2: activity coefficient calculations
if(is.null(m_star)) m_star <- IS
iH <- info("H+")
@@ -141,10 +144,10 @@
pname <- colnames(myprops)[j]
if(!pname %in% c("G", "H", "S", "Cp")) next
if(get("thermo")$opt$Setchenow == "bgamma") {
- myprops[, j] <- myprops[, j] + Setchenow(pname, IS, T, m_star, bgamma)
+ myprops[, j] <- myprops[, j] + species.sign[i] * Setchenow(pname, IS, T, m_star, bgamma)
didneutral <- TRUE
} else if(get("thermo")$opt$Setchenow == "bgamma0") {
- myprops[, j] <- myprops[, j] + Setchenow(pname, IS, T, m_star, bgamma = 0)
+ myprops[, j] <- myprops[, j] + species.sign[i] * Setchenow(pname, IS, T, m_star, bgamma = 0)
didneutral <- TRUE
}
}
@@ -153,10 +156,10 @@
pname <- colnames(myprops)[j]
if(!pname %in% c("G", "H", "S", "Cp")) next
if(method=="Alberty") {
- myprops[, j] <- myprops[, j] + Alberty(pname, Z[i], IS, T)
+ myprops[, j] <- myprops[, j] + species.sign[i] * Alberty(pname, Z[i], IS, T)
didcharged <- TRUE
} else {
- myprops[, j] <- myprops[, j] + Helgeson(pname, Z[i], IS, T, A_DH, B_DH, acirc[i], m_star, bgamma)
+ myprops[, j] <- myprops[, j] + species.sign[i] * Helgeson(pname, Z[i], IS, T, A_DH, B_DH, acirc[i], m_star, bgamma)
didcharged <- TRUE
}
}
Modified: pkg/CHNOSZ/R/subcrt.R
===================================================================
--- pkg/CHNOSZ/R/subcrt.R 2018-11-07 13:14:21 UTC (rev 345)
+++ pkg/CHNOSZ/R/subcrt.R 2018-11-07 15:50:34 UTC (rev 346)
@@ -12,7 +12,7 @@
subcrt <- function(species, coeff = 1, state = NULL, property = c("logK", "G", "H", "S", "V", "Cp"),
T = seq(273.15, 623.15, 25), P = "Psat", grid = NULL, convert = TRUE, exceed.Ttr = FALSE,
- exceed.rhomin = FALSE, logact = NULL, action.unbalanced = "warn", IS = 0) {
+ exceed.rhomin = FALSE, logact = NULL, action.unbalanced = "warn", IS = 0, is.basis = FALSE) {
# revise the call if the states have
# come as the second argument
@@ -45,6 +45,8 @@
if(length(species) > length(state)) state <- rep(state,length.out=length(species))
state <- state.args(state)
}
+ # make is.basis the same length as species
+ is.basis <- rep(is.basis, length.out=length(species))
# allowed properties
properties <- c("rho", "logK", "G", "H", "S", "Cp", "V", "kT", "E")
@@ -298,8 +300,11 @@
}
# calculate activity coefficients if ionic strength is not zero
if(any(IS != 0)) {
- 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)
+ # work out whether the species are basis species (from the is.basis argument) 20181107
+ isb <- is.basis[match(iphases[isaq], ispecies)]
+ 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, is.basis=isb)
+ else if(thermo$opt$nonideal=="Alberty") p.aq <- nonideal(iphases[isaq], p.aq, newIS, T, is.basis=isb)
}
outprops <- c(outprops, p.aq)
} else if(any(isH2O)) {
Modified: pkg/CHNOSZ/R/util.affinity.R
===================================================================
--- pkg/CHNOSZ/R/util.affinity.R 2018-11-07 13:14:21 UTC (rev 345)
+++ pkg/CHNOSZ/R/util.affinity.R 2018-11-07 15:50:34 UTC (rev 346)
@@ -133,10 +133,12 @@
if(!is.null(sout)) return(sout) else {
## subcrt arguments
species <- c(mybasis$ispecies,myspecies$ispecies)
+ is.basis <- c(rep(TRUE, length(mybasis$ispecies)), rep(FALSE, length(myspecies$ispecies)))
if("T" %in% vars) T <- vals[[which(vars=="T")]]
if("P" %in% vars) P <- vals[[which(vars=="P")]]
if("IS" %in% vars) IS <- vals[[which(vars=="IS")]]
- s.args <- list(species=species,property=property,T=T,P=P,IS=IS,grid=grid,convert=FALSE,exceed.Ttr=exceed.Ttr,exceed.rhomin=exceed.rhomin)
+ s.args <- list(species = species, property = property, T = T, P = P, IS = IS, grid = grid,
+ convert = FALSE, exceed.Ttr = exceed.Ttr, exceed.rhomin = exceed.rhomin, is.basis = is.basis)
return(do.call("subcrt",s.args))
}
}
Modified: pkg/CHNOSZ/inst/NEWS
===================================================================
--- pkg/CHNOSZ/inst/NEWS 2018-11-07 13:14:21 UTC (rev 345)
+++ pkg/CHNOSZ/inst/NEWS 2018-11-07 15:50:34 UTC (rev 346)
@@ -1,6 +1,20 @@
-CHANGES IN CHNOSZ 1.1.3-52 (2018-11-07)
+CHANGES IN CHNOSZ 1.1.3-54 (2018-11-07)
---------------------------------------
+MAJOR BUG FIXED
+
+- Previously, the virtual conversion of activities to molalities was
+ done incorrectly; the Gibbs energies of basis species and formed
+ species were transformed in the same direction. They should have
+ opposite transformations. This has been fixed by adding the 'is.basis'
+ argument to nonideal(), which is used by affinity() to transform the
+ Gibbs energies in the right direction.
+
+- For more information, see the new section of ?nonideal: 'is.basis and
+ the CHNOSZ workflow'.
+
+- Two new demos depend on the corrected behavior: gold.R and QtzMsKfs.R.
+
NEW FEATURES
- Add solubility(). Run this after affinity() to calculate the
Modified: pkg/CHNOSZ/man/macros/macros.Rd
===================================================================
--- pkg/CHNOSZ/man/macros/macros.Rd 2018-11-07 13:14:21 UTC (rev 345)
+++ pkg/CHNOSZ/man/macros/macros.Rd 2018-11-07 15:50:34 UTC (rev 346)
@@ -47,3 +47,4 @@
\newcommand{\Delta}{\ifelse{latex}{\eqn{\Delta}}{\ifelse{html}{\out{Δ}}{Δ}}}
\newcommand{\aacute}{\ifelse{latex}{\out{\'{a}}}{\ifelse{html}{\out{á}}{á}}}
\newcommand{\eacute}{\ifelse{latex}{\out{\'{e}}}{\ifelse{html}{\out{é}}{é}}}
+\newcommand{\gamma}{\ifelse{latex}{\eqn{\gamma}}{\ifelse{html}{\out{γ}}{γ}}}
Modified: pkg/CHNOSZ/man/nonideal.Rd
===================================================================
--- pkg/CHNOSZ/man/nonideal.Rd 2018-11-07 13:14:21 UTC (rev 345)
+++ pkg/CHNOSZ/man/nonideal.Rd 2018-11-07 15:50:34 UTC (rev 346)
@@ -9,7 +9,7 @@
\usage{
nonideal(species, speciesprops, IS, T, P, A_DH, B_DH,
- m_star=NULL, method=get("thermo")$opt$nonideal)
+ m_star=NULL, method=get("thermo")$opt$nonideal, is.basis=FALSE)
bgamma(TC, P, showsplines = "")
}
@@ -23,6 +23,7 @@
\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{is.basis}{logical, is (are) the specie basis species?}
\item{TC}{numeric, temperature (\degC)}
\item{showsplines}{character, show isobaric (\samp{T}) or isothermal (\samp{P}) splines}
}
@@ -75,6 +76,23 @@
This is a crude method of kriging the data, but produces fairly smooth interpolations without adding any external dependencies.
}
+\section{is.basis and the CHNOSZ workflow}{
+The main workflow in CHNOSZ (\code{\link{basis}} - \code{\link{species}} - \code{\link{affinity}} - ( \code{\link{solubility}} or \code{\link{equilibrate}} ) - \code{\link{diagram}}) is written in terms of chemical activities, not concentrations (i.e. molalities for aqueous species).
+To output molalities for the \emph{species of interest}, which are formed from the basis species, we would multiply CHNOSZ's activities by activity coefficients.
+But to obtain activities for the \emph{basis species}, we should divide the molalities that are desired in the input by activity coefficients.
+That is, to convert the entire workflow from activity to molality space requires opposite treatments for the basis species and the species being formed.
+To simplify the problem, CHNOSZ does not compute molalities by actually multiplying activities by activity coefficients (or vice versa) -- this would require complex calculations of activity coefficients at the input and output stages, considering the many possible dimensions of system variables -- a true mess!
+Instead, the same effect is obtained at the core of the workflow by using standard Gibbs energies adjusted for given ionic strength (i.e. transformed Gibbs energies).
+The transformation is very simple: by adding \emph{RT}\gamma (\gamma is the activity coefficient calculated at the appropriate \emph{T}, \emph{P}, and ionic strength) to the standard Gibbs energy, all expressions for activity of that species are converted to molality.
+That transformation is consistent with the requirements for the species of interest.
+The reverse transformation, subtracting \emph{RT}\gamma from the standard Gibbs energy, is needed for the basis species.
+
+The \code{is.basis} argument controls the direction of the transformation.
+In general, it should not be needed by the user, but is used by \code{affinity} to obtain the correctly transformed Gibbs energies.
+Thus, by activating nonideality calculations in \code{affinity} (with a non-zero \code{IS} argument), the activity variables, such as \code{logact} in the \code{basis} and \code{species} definitions, and \code{loga.equil} and \code{loga.balance} in the downstream calculations, are converted to molalities.
+Actually renaming the variables in the code is not possible, but \code{\link{diagram}} changes the plot labels to molalities if it is provided results from a calculation with non-zero \code{IS} set in \code{affinity}.
+}
+
\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 adjusted thermodynamic properties (at higher IS).
For all affected species, a column named \code{loggam} (common (base-10) logarithm of gamma, the activity coefficient) is appended to the output dataframe of species properties.
Modified: pkg/CHNOSZ/man/subcrt.Rd
===================================================================
--- pkg/CHNOSZ/man/subcrt.Rd 2018-11-07 13:14:21 UTC (rev 345)
+++ pkg/CHNOSZ/man/subcrt.Rd 2018-11-07 15:50:34 UTC (rev 346)
@@ -11,7 +11,7 @@
property = c("logK","G","H","S","V","Cp"),
T = seq(273.15,623.15,25), P = "Psat", grid = NULL,
convert = TRUE, exceed.Ttr = FALSE, exceed.rhomin = FALSE,
- logact = NULL, action.unbalanced = "warn", IS = 0)
+ logact = NULL, action.unbalanced = "warn", IS = 0, is.basis = FALSE)
}
\arguments{
@@ -28,6 +28,7 @@
\item{convert}{logical, are input and output units of T and P those of the user (\code{TRUE}) (see \code{\link{T.units}}), or are they Kelvin and bar (\code{FALSE})?}
\item{action.unbalanced}{character \samp{warn} or NULL, what action to take if unbalanced reaction is provided}
\item{IS}{numeric, ionic strength(s) at which to calculate adjusted molal properties, mol kg\eqn{^{-1}}{^-1}}
+ \item{is.basis}{logical, is (are) the species basis species? See \code{\link{nonideal}}}
}
\details{
Modified: pkg/CHNOSZ/tests/testthat/test-logmolality.R
===================================================================
--- pkg/CHNOSZ/tests/testthat/test-logmolality.R 2018-11-07 13:14:21 UTC (rev 345)
+++ pkg/CHNOSZ/tests/testthat/test-logmolality.R 2018-11-07 15:50:34 UTC (rev 346)
@@ -5,6 +5,7 @@
# and in the rest of the main workflow of CHNOSZ?
# 20171025 first version
# 20181106 include non-zero activity coefficient of CO2(aq)
+ # 20181107 include 'is.basis' and opposite transformations for basis species and formed species
### first get the activity coefficients of H+ and HCO3-
## the long way...
@@ -70,7 +71,9 @@
## case 2: IS = 1
a1 <- affinity(IS=1)
A1affinity <- -convert(a1$values[[2]], "G")
- expect_equal(A1affinity[[1]], A1subcrt)
+ # we had better use is.basis here, which indicates the direction of transformation of Gibbs energy 20181107
+ A1subcrt.trans <- subcrt(c("CO2", "H2O", "H+", "HCO3-"), c(-1, -1, 1, 1), T=25, logact=c(-3, 0, -7, -3), IS=1, is.basis=c(TRUE, TRUE, TRUE, FALSE))$out$A
+ expect_equal(A1affinity[[1]], A1subcrt.trans)
## take-home message 2: using affinity() with IS not equal to zero, the "logact"
## set by species() is logmolal in affinity calculations for charged aqueous species
@@ -91,7 +94,10 @@
# so, logK = 6.345
logKrev <- -logK
logQrev0 <- -logQ0
- logQrev1 <- -logQ1
+ # note the minus sign here, because HCO3 is now a basis species
+ # and has the opposite Gibbs energy transformation 20181107
+ logaHCO3 <- -3 - loggam_HCO3
+ logQrev1 <- (0 + logaCO2) - (-7 + logaHCO3)
ACO2_0manual <- -convert(logKrev - logQrev0, "G")
ACO2_1manual <- -convert(logKrev - logQrev1, "G")
expect_equal(ACO2_0manual, ACO2_0affinity[[1]])
@@ -112,8 +118,9 @@
# case 2: IS = 1
logact_HCO3 <- e1$loga.equil[[2]]
logact_CO2 <- e1$loga.equil[[1]]
- # here, loga.equil is the *molality*, so we must multiply by loggam
- logQeq1 <- (-7 + logact_HCO3 + loggam_HCO3) - (logact_CO2 + loggam_CO2 + 0)
+ # CO2 (formed species): convert log activity to log molality (multiply by loggam)
+ # HCO3- (basis species): convert log molality to log activity (divide by loggam)
+ logQeq1 <- (-7 + logact_HCO3 - loggam_HCO3) - (logact_CO2 + loggam_CO2 + 0)
Aeq1 <- -convert(logK - logQeq1, "G") # zero!
expect_equal(Aeq1[[1]], 0)
## take-home message 4: using affinity() with IS not equal to zero, the "loga.equil"
More information about the CHNOSZ-commits
mailing list