[CHNOSZ-commits] r348 - in pkg/CHNOSZ: . R demo inst man tests/testthat
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Nov 8 14:59:57 CET 2018
Author: jedick
Date: 2018-11-08 14:59:57 +0100 (Thu, 08 Nov 2018)
New Revision: 348
Removed:
pkg/CHNOSZ/demo/QtzMsKfs.R
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/R/examples.R
pkg/CHNOSZ/R/nonideal.R
pkg/CHNOSZ/R/subcrt.R
pkg/CHNOSZ/R/util.affinity.R
pkg/CHNOSZ/demo/00Index
pkg/CHNOSZ/demo/gold.R
pkg/CHNOSZ/inst/NEWS
pkg/CHNOSZ/man/berman.Rd
pkg/CHNOSZ/man/examples.Rd
pkg/CHNOSZ/man/nonideal.Rd
pkg/CHNOSZ/man/subcrt.Rd
pkg/CHNOSZ/tests/testthat/test-logmolality.R
Log:
revert 'is.basis' changes of revision 346; use molality of Cl- in demo/gold.R
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2018-11-08 06:50:22 UTC (rev 347)
+++ pkg/CHNOSZ/DESCRIPTION 2018-11-08 13:59:57 UTC (rev 348)
@@ -1,6 +1,6 @@
Date: 2018-11-08
Package: CHNOSZ
-Version: 1.1.3-55
+Version: 1.1.3-56
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/examples.R
===================================================================
--- pkg/CHNOSZ/R/examples.R 2018-11-08 06:50:22 UTC (rev 347)
+++ pkg/CHNOSZ/R/examples.R 2018-11-08 13:59:57 UTC (rev 348)
@@ -28,7 +28,7 @@
demos <- function(which=c("sources", "protein.equil", "affinity", "NaCl", "density",
"ORP", "revisit", "findit", "ionize", "buffer", "protbuff", "yeastgfp", "mosaic",
- "copper", "solubility", "gold", "QtzMsKfs", "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-08 06:50:22 UTC (rev 347)
+++ pkg/CHNOSZ/R/nonideal.R 2018-11-08 13:59:57 UTC (rev 348)
@@ -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, is.basis=FALSE) {
+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
@@ -123,9 +123,6 @@
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+")
@@ -144,10 +141,10 @@
pname <- colnames(myprops)[j]
if(!pname %in% c("G", "H", "S", "Cp")) next
if(get("thermo")$opt$Setchenow == "bgamma") {
- myprops[, j] <- myprops[, j] + species.sign[i] * Setchenow(pname, IS, T, m_star, bgamma)
+ myprops[, j] <- myprops[, j] + Setchenow(pname, IS, T, m_star, bgamma)
didneutral <- TRUE
} else if(get("thermo")$opt$Setchenow == "bgamma0") {
- myprops[, j] <- myprops[, j] + species.sign[i] * Setchenow(pname, IS, T, m_star, bgamma = 0)
+ myprops[, j] <- myprops[, j] + Setchenow(pname, IS, T, m_star, bgamma = 0)
didneutral <- TRUE
}
}
@@ -156,10 +153,10 @@
pname <- colnames(myprops)[j]
if(!pname %in% c("G", "H", "S", "Cp")) next
if(method=="Alberty") {
- myprops[, j] <- myprops[, j] + species.sign[i] * Alberty(pname, Z[i], IS, T)
+ myprops[, j] <- myprops[, j] + Alberty(pname, Z[i], IS, T)
didcharged <- TRUE
} else {
- myprops[, j] <- myprops[, j] + species.sign[i] * Helgeson(pname, Z[i], IS, T, A_DH, B_DH, acirc[i], m_star, bgamma)
+ myprops[, j] <- myprops[, j] + 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-08 06:50:22 UTC (rev 347)
+++ pkg/CHNOSZ/R/subcrt.R 2018-11-08 13:59:57 UTC (rev 348)
@@ -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, is.basis = FALSE) {
+ exceed.rhomin = FALSE, logact = NULL, action.unbalanced = "warn", IS = 0) {
# revise the call if the states have
# come as the second argument
@@ -45,8 +45,6 @@
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")
@@ -300,11 +298,8 @@
}
# calculate activity coefficients if ionic strength is not zero
if(any(IS != 0)) {
- # 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)
+ 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)
} else if(any(isH2O)) {
Modified: pkg/CHNOSZ/R/util.affinity.R
===================================================================
--- pkg/CHNOSZ/R/util.affinity.R 2018-11-08 06:50:22 UTC (rev 347)
+++ pkg/CHNOSZ/R/util.affinity.R 2018-11-08 13:59:57 UTC (rev 348)
@@ -133,12 +133,10 @@
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, is.basis = is.basis)
+ 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)
return(do.call("subcrt",s.args))
}
}
Modified: pkg/CHNOSZ/demo/00Index
===================================================================
--- pkg/CHNOSZ/demo/00Index 2018-11-08 06:50:22 UTC (rev 347)
+++ pkg/CHNOSZ/demo/00Index 2018-11-08 13:59:57 UTC (rev 348)
@@ -14,7 +14,6 @@
copper Another example of mosaic(): complexation of copper with glycine species
solubility Solubility of calcite and CO2(gas) as a function of pH
gold Solubility of gold
-QtzMsKfs Helgeson and Berman minerals and calculations with molality
wjd Gibbs energy minimization: prebiological atmospheres and cell periphery of yeast
dehydration log K of dehydration reactions; SVG file contains tooltips and links
bugstab Formation potential of microbial proteins in colorectal cancer
Deleted: pkg/CHNOSZ/demo/QtzMsKfs.R
===================================================================
--- pkg/CHNOSZ/demo/QtzMsKfs.R 2018-11-08 06:50:22 UTC (rev 347)
+++ pkg/CHNOSZ/demo/QtzMsKfs.R 2018-11-08 13:59:57 UTC (rev 348)
@@ -1,69 +0,0 @@
-# CHNOSZ/demo/QtzMsKfs.R
-# T - log(K+/H+) diagram, after Sverjensky et al., 1991
-# (doi:10.1016/0016-7037(91)90157-Z)
-# 20171009 diagram added to berman.Rd
-# 20181108 moved to demo/QtzMsKfs.R; add molality calculations
-
-# this demo compares diagrams made using the Berman and Helgeson datasets,
-# and shows the use of nonideal calculations to set molalities in the basis species
-
-## set up the system: basis species
-basis(c("K+", "Al+3", "quartz", "H2O", "O2", "H+"))
-# use pH = 0 so that aK+ = aK+/aH+
-basis("pH", 0)
-# load the species
-species(c("K-feldspar", "muscovite", "kaolinite",
- "pyrophyllite", "andalusite"), "cr")
-## the "b_gamma" equation gets closer to the published diagram than "B-dot"
-thermo$opt$nonideal <<- "bgamma"
-
-## start with the data from Helgeson et al., 1978
-add.obigt("SUPCRT92")
-# calculate affinities in aK+ - temperature space
-# exceed.Tr: we go above stated temperature limit of pyrophyllite
-# (this is above its stability field on the diagram, so pyrophyllite doesn't appear in this region,
-# but its properties are needed needed to calculate relative stabilities of all minerals)
-res <- 400
-a <- affinity(`K+` = c(0, 5, res), T = c(200, 650, res), P = 1000, exceed.Ttr = TRUE)
-# make base plot with colors and no lines
-diagram(a, xlab = ratlab("K+", use.molality = TRUE), lty = 0, fill = "terrain")
-# add the lines, extending into the low-density region (exceed.rhomin = TRUE)
-a <- affinity(`K+` = c(0, 5, res), T = c(200, 650, res), P = 1000, exceed.Ttr = TRUE, exceed.rhomin = TRUE)
-diagram(a, add = TRUE, names = NULL, col = "red", lty = 2, lwd = 1.5)
-# calculate and plot the lines for 1 molal chloride
-a <- affinity(`K+` = c(0, 5, res), T = c(200, 650, res), P = 1000, exceed.Ttr = TRUE, exceed.rhomin = TRUE, IS = 1)
-diagram(a, add = TRUE, names = NULL, col = "red", lwd = 1.5)
-# the list of references:
-ref1 <- thermo.refs(species()$ispecies)$key
-
-## now use the (default) data from Berman, 1988
-# this resets the thermodynamic database
-# without affecting the basis and species settings
-data(OBIGT)
-# we can check that we have Berman's quartz
-# and not coesite or some other phase of SiO2
-iSiO2 <- rownames(basis()) == "SiO2"
-stopifnot(info(basis()$ispecies[iSiO2])$name == "quartz")
-# Berman's dataset doesn't have the upper temperature limits, so we don't need exceed.Ttr here
-a <- affinity(`K+` = c(0, 5, res), T = c(200, 650, res), P = 1000, exceed.rhomin = TRUE)
-diagram(a, add = TRUE, names = NULL, col = "blue", lty = 2, lwd = 1.5)
-a <- affinity(`K+` = c(0, 5, res), T = c(200, 650, res), P = 1000, exceed.rhomin = TRUE, IS = 1)
-diagram(a, add = TRUE, names = NULL, col = "blue", lwd = 1.5)
-# the list of references:
-ref2 <- thermo.refs(species()$ispecies)$key
-ref2 <- paste(ref2, collapse = ", ")
-
-# add experimental points for 1000 bar (Table 1 of Sverjensky et al., 1991)
-expt.T <- c(300, 400, 500, 550, # KFs-Ms-Qtz
- 400, 450, 500, 550, # Ms-And-Qtz
- 300, 350, # Ms-P-Qtz
- 300, 600) # Kaol-Ms-Qtz, KFs-And-Qtz
-expt.KH <- c(3.50, 2.75, 1.95, 1.40, 1.60, 1.57, 1.47, 1.38, 1.94, 1.80, 1.90, 0.63)
-points(expt.KH, expt.T, pch = 19, cex = 1.2)
-# add legend and title
-legend("top", "low-density region", text.font = 3, bty = "n")
-legend("top", describe.property(c(NA, NA, "P", "IS"), c(NA, NA, 1000, 1)), bty = "n")
-legend("left", c(ref1, ref2, "ion molality", "ion activity", "experiments"),
- lty = c(1, 1, 1, 2, 0), lwd = 1.5, col = c(2, 4, 1, 1, 1), pch = c(NA, NA, NA, NA, 19), bty = "n")
-title(main = syslab(c("K2O", "Al2O3", "SiO2", "H2O", "HCl")), line = 1.8)
-title(main = "Helgeson and Berman minerals, after Sverjensky et al., 1991", line = 0.3, font.main = 1)
Modified: pkg/CHNOSZ/demo/gold.R
===================================================================
--- pkg/CHNOSZ/demo/gold.R 2018-11-08 06:50:22 UTC (rev 347)
+++ pkg/CHNOSZ/demo/gold.R 2018-11-08 13:59:57 UTC (rev 348)
@@ -147,7 +147,7 @@
dK <- describe.basis(ibasis=5, use.molality=TRUE)
legend("topleft", c(dP, dNaCl, dK), bty = "n")
dbasis <- describe.basis(ibasis = c(9, 7, 10))
- legend("topright", dbasis, bty = "n")
+ legend(320, -4, dbasis, bty = "n")
title(main=("After Williams-Jones et al., 2009, Fig. 2B"), font.main = 1)
}
@@ -156,7 +156,7 @@
# (doi:10.1144/SP402.4)
Au_T2 <- function() {
species(c("Au(HS)2-", "Au(HS)", "AuOH", "AuCl2-"))
- # approximate activity of H2S for total S = 0.01 m
+ # total S = 0.01 m
basis("H2S", -2)
# apply HM buffer for fO2
basis("O2", "HM")
@@ -164,12 +164,12 @@
basis("H+", "QMK")
basis("K+", log10(0.5))
# calculate solution composition for 2 mol/kg NaCl
- NaCl <- NaCl(T = seq(150, 550, 10), P = 1000, m_tot=1)
+ NaCl <- NaCl(T = seq(150, 550, 10), P = 1000, m_tot=2)
# calculate affinity and solubility
a <- affinity(T = seq(150, 550, 10), `Cl-` = log10(NaCl$m_Cl), P = 1000, IS = NaCl$IS)
s <- solubility(a)
# make diagram and show total log molality
- diagram(s, ylim = c(-10, -2), col = col, lwd = 2, lty = 1)
+ diagram(s, ylim = c(-10, -3), col = col, lwd = 2, lty = 1)
diagram(s, add = TRUE, type = "loga.balance", lwd = 3, lty = 2)
# make legend and title
dP <- describe.property("P", 1000)
@@ -177,7 +177,7 @@
dK <- describe.basis(ibasis=5, use.molality=TRUE)
legend("topleft", c(dP, dNaCl, dK), bty = "n")
dbasis <- describe.basis(ibasis = c(9, 7, 10))
- legend("topright", dbasis, bty = "n")
+ legend(320, -3, dbasis, bty = "n")
title(main=("After Williams-Jones et al., 2009, Fig. 2A"), font.main = 1)
}
Modified: pkg/CHNOSZ/inst/NEWS
===================================================================
--- pkg/CHNOSZ/inst/NEWS 2018-11-08 06:50:22 UTC (rev 347)
+++ pkg/CHNOSZ/inst/NEWS 2018-11-08 13:59:57 UTC (rev 348)
@@ -1,23 +1,6 @@
-CHANGES IN CHNOSZ 1.1.3-55 (2018-11-08)
+CHANGES IN CHNOSZ 1.1.3-56 (2018-11-08)
---------------------------------------
-MAJOR BUG FIXED
-
-- Previously, with calculations using nonideal(), the Gibbs energies of
- basis species and formed species were transformed in the same
- direction, causing an incorrct conversion of activities to molalities.
- This has been fixed by adding the 'is.basis' argument to nonideal(),
- which is used by affinity() to transform the Gibbs energies in opposite
- directions for basis species and formed species.
-
-- 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.
- The latter demo also provides a comparison of the superseded Helgeson
- (SUPCRT92) and newer Berman datasets for minerals. See below for more
- information on demo/gold.R.
-
NEW FEATURES
- Add solubility(). Run this after affinity() to calculate the
Modified: pkg/CHNOSZ/man/berman.Rd
===================================================================
--- pkg/CHNOSZ/man/berman.Rd 2018-11-08 06:50:22 UTC (rev 347)
+++ pkg/CHNOSZ/man/berman.Rd 2018-11-08 13:59:57 UTC (rev 348)
@@ -94,6 +94,53 @@
DGrxn <- G_Cs - G_aQz
stopifnot(all(abs(DGrxn) < 100))
+### compare mineral stabilities in the Berman and Helgeson datasets
+### on a T - log(K+/H+) diagram, after Sverjensky et al., 1991
+### (doi:10.1016/0016-7037(91)90157-Z)
+## set up the system: basis species
+basis(c("K+", "Al+3", "quartz", "H2O", "O2", "H+"))
+# use pH = 0 so that aK+ = aK+/aH+
+basis("pH", 0)
+# load the species
+species(c("K-feldspar", "muscovite", "kaolinite",
+ "pyrophyllite", "andalusite"), "cr")
+## start with the data from Helgeson et al., 1978
+add.obigt("SUPCRT92")
+# calculate affinities in aK+ - temperature space
+# exceed.Tr: enable calculations above stated temperature limit of pyrophyllite
+res <- 400
+a <- affinity(`K+` = c(0, 5, res), T = c(200, 650, res), P = 1000, exceed.Ttr = TRUE)
+# make base plot with colors and no lines
+diagram(a, xlab = ratlab("K+", use.molality = TRUE), lty = 0, fill = "terrain")
+# add the lines, extending into the low-density region (exceed.rhomin = TRUE)
+a <- affinity(`K+` = c(0, 5, res), T = c(200, 650, res), P = 1000,
+ exceed.Ttr = TRUE, exceed.rhomin = TRUE)
+diagram(a, add = TRUE, names = NULL, col = "red", lwd = 1.5)
+# the list of references:
+ref1 <- thermo.refs(species()$ispecies)$key
+## now use the (default) data from Berman, 1988
+# this resets the thermodynamic database
+# without affecting the basis and species settings
+data(OBIGT)
+# we can check that we have Berman's quartz
+# and not coesite or some other phase of SiO2
+iSiO2 <- rownames(basis()) == "SiO2"
+stopifnot(info(basis()$ispecies[iSiO2])$name == "quartz")
+# Berman's dataset doesn't have the upper temperature limits, so we don't need exceed.Ttr here
+a <- affinity(`K+` = c(0, 5, res), T = c(200, 650, res), P = 1000, exceed.rhomin = TRUE)
+diagram(a, add = TRUE, names = NULL, col = "blue", lwd = 1.5)
+# the list of references:
+ref2 <- thermo.refs(species()$ispecies)$key
+ref2 <- paste(ref2, collapse = ", ")
+# add legend and title
+legend("top", "low-density region", text.font = 3, bty = "n")
+legend("topleft", describe.property(c("P", "IS"), c(1000, 1)), bty = "n")
+legend("left", c(ref1, ref2),
+ lty = c(1, 1), lwd = 1.5, col = c(2, 4), bty = "n")
+title(main = syslab(c("K2O", "Al2O3", "SiO2", "H2O", "HCl")), line = 1.8)
+title(main = "Helgeson and Berman minerals, after Sverjensky et al., 1991",
+ line = 0.3, font.main = 1)
+
# make a P-T diagram for SiO2 minerals (Ber88 Fig. 4)
basis(c("SiO2", "O2"), c("cr", "gas"))
species(c("quartz", "quartz,beta", "coesite"), "cr")
Modified: pkg/CHNOSZ/man/examples.Rd
===================================================================
--- pkg/CHNOSZ/man/examples.Rd 2018-11-08 06:50:22 UTC (rev 347)
+++ pkg/CHNOSZ/man/examples.Rd 2018-11-08 13:59:57 UTC (rev 348)
@@ -16,7 +16,7 @@
demos(which = c("sources", "protein.equil", "affinity", "NaCl",
"density", "ORP", "revisit", "findit", "ionize", "buffer",
"protbuff", "yeastgfp", "mosaic", "copper", "solubility",
- "gold", "QtzMsKfs", "wjd", "bugstab", "Shh", "activity_ratios",
+ "gold", "wjd", "bugstab", "Shh", "activity_ratios",
"adenine", "DEW", "lambda", "TCA", "go-IU", "bison"),
save.png=FALSE)
}
@@ -47,7 +47,6 @@
\code{copper} \tab * Another example of \code{\link{mosaic}}: complexation of Cu with glycine (Aksu and Doyle, 2001) \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{QtzMsKfs} \tab * Helgeson and Berman minerals and calculations with molality (Sverjensky et al., 1991) \cr
\code{wjd} \tab * \eqn{G}{G} minimization: prebiological atmospheres (Dayhoff et al., 1964) and cell periphery of yeast \cr
\code{dehydration} \tab * \logK of dehydration reactions; SVG file contains tooltips and links \cr
\code{bugstab} \tab * Formation potential of microbial proteins in colorectal cancer (Dick, 2016) \cr
Modified: pkg/CHNOSZ/man/nonideal.Rd
===================================================================
--- pkg/CHNOSZ/man/nonideal.Rd 2018-11-08 06:50:22 UTC (rev 347)
+++ pkg/CHNOSZ/man/nonideal.Rd 2018-11-08 13:59:57 UTC (rev 348)
@@ -9,7 +9,7 @@
\usage{
nonideal(species, speciesprops, IS, T, P, A_DH, B_DH,
- m_star=NULL, method=get("thermo")$opt$nonideal, is.basis=FALSE)
+ m_star=NULL, method=get("thermo")$opt$nonideal)
bgamma(TC, P, showsplines = "")
}
@@ -23,7 +23,6 @@
\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}
}
@@ -76,23 +75,6 @@
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-08 06:50:22 UTC (rev 347)
+++ pkg/CHNOSZ/man/subcrt.Rd 2018-11-08 13:59:57 UTC (rev 348)
@@ -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, is.basis = FALSE)
+ logact = NULL, action.unbalanced = "warn", IS = 0)
}
\arguments{
@@ -28,7 +28,6 @@
\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-08 06:50:22 UTC (rev 347)
+++ pkg/CHNOSZ/tests/testthat/test-logmolality.R 2018-11-08 13:59:57 UTC (rev 348)
@@ -5,7 +5,6 @@
# 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...
@@ -71,9 +70,7 @@
## case 2: IS = 1
a1 <- affinity(IS=1)
A1affinity <- -convert(a1$values[[2]], "G")
- # 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)
+ expect_equal(A1affinity[[1]], A1subcrt)
## 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
@@ -94,10 +91,7 @@
# so, logK = 6.345
logKrev <- -logK
logQrev0 <- -logQ0
- # 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)
+ logQrev1 <- -logQ1
ACO2_0manual <- -convert(logKrev - logQrev0, "G")
ACO2_1manual <- -convert(logKrev - logQrev1, "G")
expect_equal(ACO2_0manual, ACO2_0affinity[[1]])
@@ -118,9 +112,8 @@
# case 2: IS = 1
logact_HCO3 <- e1$loga.equil[[2]]
logact_CO2 <- e1$loga.equil[[1]]
- # 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)
+ # here, loga.equil is the *molality*, so we must multiply 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