[CHNOSZ-commits] r457 - in pkg/CHNOSZ: . R inst man tests/testthat
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sun May 5 04:49:43 CEST 2019
Author: jedick
Date: 2019-05-05 04:49:37 +0200 (Sun, 05 May 2019)
New Revision: 457
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/R/mosaic.R
pkg/CHNOSZ/inst/NEWS
pkg/CHNOSZ/man/mosaic.Rd
pkg/CHNOSZ/tests/testthat/test-mosaic.R
Log:
mosaic(): bug fix: mole fractions of basis species affect their activities in affinity calculations
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2019-05-04 11:38:52 UTC (rev 456)
+++ pkg/CHNOSZ/DESCRIPTION 2019-05-05 02:49:37 UTC (rev 457)
@@ -1,6 +1,6 @@
-Date: 2019-05-04
+Date: 2019-05-05
Package: CHNOSZ
-Version: 1.3.2-4
+Version: 1.3.2-5
Title: Thermodynamic Calculations and Diagrams for Geochemistry
Authors at R: c(
person("Jeffrey", "Dick", , "j3ffdick at gmail.com", role = c("aut", "cre"),
Modified: pkg/CHNOSZ/R/mosaic.R
===================================================================
--- pkg/CHNOSZ/R/mosaic.R 2019-05-04 11:38:52 UTC (rev 456)
+++ pkg/CHNOSZ/R/mosaic.R 2019-05-05 02:49:37 UTC (rev 457)
@@ -3,6 +3,7 @@
# 20141220 jmd initial version
# 20190129 complete rewrite to use any number of groups of changing basis species
# and improve speed by pre-calculating subcrt values (sout)
+# 20190505 bug fix: adjust affinities of species formation reactions for mole fractions of basis species
## if this file is interactively sourced, the following are also needed to provide unexported functions:
#source("basis.R")
@@ -10,7 +11,7 @@
#source("util.args.R")
# function to calculate affinities with mosaic of basis species
-mosaic <- function(bases, bases2 = NULL, blend = TRUE, mixing = TRUE, ...) {
+mosaic <- function(bases, bases2 = NULL, blend = TRUE, ...) {
# argument recall 20190120
# if the first argument is the result from a previous mosaic() calculation,
@@ -41,7 +42,7 @@
hasbases2 <- TRUE
}
otherargs <- list(...)
- allargs <- c(list(bases = bases, blend = blend, mixing = mixing), otherargs)
+ allargs <- c(list(bases = bases, blend = blend), otherargs)
out <- do.call(mosaic, allargs)
# replace A.bases (affinity calculations for all groups of basis species) with backwards-compatbile A.bases and A.bases2
if(hasbases2) A.bases2 <- out$A.bases[[2]]
@@ -103,7 +104,7 @@
# calculate equilibrium mole fractions for each group of basis species
group.fraction <- list()
if(blend) {
- e.out <- list()
+ e.bases <- list()
for(i in 1:length(A.bases)) {
# this isn't needed (and doesn't work) if all the affinities are NA 20180925
if(any(!sapply(A.bases[[1]]$values, is.na))) {
@@ -113,7 +114,7 @@
a.tot <- Reduce("+", a.equil)
group.fraction[[i]] <- lapply(a.equil, function(x) x / a.tot)
# include the equilibrium activities in the output of this function 20190504
- e.out[[1]] <- e
+ e.bases[[1]] <- e
} else {
group.fraction[[i]] <- A.bases[[i]]$values
}
@@ -139,17 +140,28 @@
for(i in 1:length(ispecies)) ind.mat[[i]] <- 1:length(ispecies[[i]])
ind.mat <- as.matrix(expand.grid(ind.mat))
- # calculate mole fractions for each combination of basis species
- for(i in 1:nrow(ind.mat)) {
- # multiply fractions from each group
- for(j in 1:ncol(ind.mat)) {
- if(j==1) x <- group.fraction[[j]][[ind.mat[i, j]]]
- else x <- x * group.fraction[[j]][[ind.mat[i, j]]]
+ # loop over combinations of basis species
+ for(icomb in 1:nrow(ind.mat)) {
+ # loop over groups of changing basis species
+ for(igroup in 1:ncol(ind.mat)) {
+ # get mole fractions for this particular basis species
+ basisx <- group.fraction[[igroup]][[ind.mat[icomb, igroup]]]
+ # loop over species
+ for(jspecies in 1:length(aff.species[[icomb]]$values)) {
+ # get coefficient of this basis species in the formation reaction for this species
+ nbasis <- aff.species[[icomb]]$species[jspecies, ibasis0[igroup]]
+ # adjust affinity of species for mole fractions (i.e. lower activity) of basis species 20190505
+ aff.adjust <- nbasis * log10(basisx)
+ # avoid infinite values (from log10(0))
+ isfin <- is.finite(aff.adjust)
+ aff.species[[icomb]]$values[[jspecies]][isfin] <- aff.species[[icomb]]$values[[jspecies]][isfin] + aff.adjust[isfin]
+ }
+ # multiply fractions of basis species from each group to get overall fraction
+ if(igroup==1) groupx <- basisx
+ else groupx <- groupx * basisx
}
# multiply affinities by the mole fractions of basis species
- # include mixing term (-x*log10(x)) 20190121
- if(blend & mixing) aff.species[[i]]$values <- lapply(aff.species[[i]]$values, function(values) values * x - x * log10(x))
- else aff.species[[i]]$values <- lapply(aff.species[[i]]$values, function(values) values * x)
+ aff.species[[icomb]]$values <- lapply(aff.species[[icomb]]$values, function(values) values * groupx)
}
# get total affinities for the species
@@ -162,8 +174,8 @@
}
# for argument recall, include all arguments in output 20190120
- allargs <- c(list(bases = bases, blend = blend, mixing = mixing), list(...))
+ allargs <- c(list(bases = bases, blend = blend), list(...))
# return the affinities for the species and basis species
- if(blend) return(list(fun = "mosaic", args = allargs, A.species = A.species, A.bases = A.bases, e.out = e.out))
+ if(blend) return(list(fun = "mosaic", args = allargs, A.species = A.species, A.bases = A.bases, e.bases = e.bases))
else return(list(fun = "mosaic", args = allargs, A.species = A.species, A.bases = A.bases))
}
Modified: pkg/CHNOSZ/inst/NEWS
===================================================================
--- pkg/CHNOSZ/inst/NEWS 2019-05-04 11:38:52 UTC (rev 456)
+++ pkg/CHNOSZ/inst/NEWS 2019-05-05 02:49:37 UTC (rev 457)
@@ -1,14 +1,41 @@
-CHANGES IN CHNOSZ 1.3.2-4 (2019-05-04)
+CHANGES IN CHNOSZ 1.3.2-5 (2019-05-05)
--------------------------------------
-- demo/aluminum.R: add calculations using Si(OH)4 from Akinfiev-Diamond
- model (SiO2 in these reactions is replaced by Si(OH)4 - 2 H2O).
+BUG FIXES
+- In mosaic(), affinities of formation of species should not only be
+ weighted by the mole fraction of the changing basis species, but they
+ should be calculated using activities of basis species that are
+ adjusted (i.e. lowered) according to their mole fractions. Previously,
+ constant activities were used for each changing basies species,
+ leading to an artificial inflation of affinities of species,
+ especially near the transitions where the changing basis species are
+ nearly equal in abundance.
+
+- The 'mixing' argument of mosaic() has been removed, as it no longer
+ provides correct results and is not needed with the improved handling
+ of activities of basis species.
+
+OTHER CHANGES
+
+- Tests have been added to test-mosaic.R to check that activitives
+ produced by mosaic() - equilibrate() and mosaic() - solubility()
+ really are equilibrium activities, i.e. that the affinities of
+ reactions between species are equivalent to zero, especially near the
+ transitions of basis species.
+
- In mosaic(), equilibrate groups of changing basis species using total
activity taken from the corresponding basis species in the incoming
- basis() definition, and return the result of equilibrate() in the
- output.
+ basis() definition. This is not directly related to the bug above, but
+ provides a means to constrain the amount of an element in the basis
+ species.
+- In mosaic(), include the result of equilibrate() for each group of
+ basis species in the output ('e.bases').
+
+- demo/aluminum.R: add calculations using Si(OH)4 from Akinfiev-Diamond
+ model (SiO2 in these reactions is replaced by Si(OH)4 - 2 H2O).
+
CHANGES IN CHNOSZ 1.3.2 (2019-04-20)
------------------------------------
Modified: pkg/CHNOSZ/man/mosaic.Rd
===================================================================
--- pkg/CHNOSZ/man/mosaic.Rd 2019-05-04 11:38:52 UTC (rev 456)
+++ pkg/CHNOSZ/man/mosaic.Rd 2019-05-05 02:49:37 UTC (rev 457)
@@ -7,7 +7,7 @@
}
\usage{
- mosaic(bases, bases2 = NULL, blend = TRUE, mixing = TRUE, ...)
+ mosaic(bases, bases2 = NULL, blend = TRUE, ...)
}
\arguments{
@@ -14,7 +14,6 @@
\item{bases}{character, basis species to be changed in the calculation, or list, containing vectors for each group of changing basis species}
\item{bases2}{character, second set of changing basis species}
\item{blend}{logical, use relative abundances of basis species?}
- \item{mixing}{logical, include a term for the Gibbs energy of mixing?}
\item{...}{additional arguments to be passed to \code{\link{affinity}}}
}
@@ -32,9 +31,7 @@
If \code{blend} is TRUE (the default), the relative abundances of the basis species in each group are calculated using \code{\link{equilibrate}}, with the total activity taken from the corresponding basis species in the incoming \code{\link{basis}} definition.
Then, the function calculates overall affinities of the formation reactions of each species by combining reactions written using individual basis species in proportion to the relative abundances of the basis species.
-Additionally, if \code{mixing} is TRUE (the default), a term is included to account for the Gibbs energy of ideal mixing of the basis species.
-See the second example in \code{\link{solubility}} for a numerical test of the calculations using \code{blend} and \code{mixing}.
-If \code{blend} is FALSE, the function returns the affinities calculated using the single predominant basis species in \code{bases} at each condition (in this case, the \code{mixing} argument has no effect).
+If \code{blend} is FALSE, the function returns the affinities calculated using the single predominant basis species in \code{bases} at each condition.
A more flexible method of specifying multiple sets of basis species is now available.
Instead of using \code{bases} and \code{bases2}, supply a list for just the \code{bases} argument.
@@ -46,12 +43,13 @@
\value{
A list containing \code{A.species} (affinities of formation of the species with changing basis species) and \code{A.bases} (affinities of formation of the basis species in terms of the first basis species), each having same structure as the list returned by \code{\link{affinity}}.
-If \code{blend} is TRUE, the output also contains \code{e.out} (the output of \code{\link{equilibrate}} for each group of basis species)
+If \code{blend} is TRUE, the output also contains \code{e.bases} (the output of \code{\link{equilibrate}} for each group of basis species)
If \code{bases2} is provided, the list also contains \code{A.bases2} (affinities of formation of the second group of basis species).
}
\seealso{
\code{demo("mosaic")}, extending the example below by addition of carbonate species in \code{bases2}, and using thermodynamic data from Garrels and Christ, 1965.
+The help page of \code{\link{solubility}} has an example combining \code{mosaic} with solubility calculations.
}
\examples{
Modified: pkg/CHNOSZ/tests/testthat/test-mosaic.R
===================================================================
--- pkg/CHNOSZ/tests/testthat/test-mosaic.R 2019-05-04 11:38:52 UTC (rev 456)
+++ pkg/CHNOSZ/tests/testthat/test-mosaic.R 2019-05-05 02:49:37 UTC (rev 457)
@@ -35,13 +35,10 @@
# calculate affinities using the predominant basis species
pH <- c(0, 14, 29)
m1 <- mosaic(bases, pH = pH, blend = FALSE)
- # calculate affinities with smooth transitions between basis species, including a mixing energy
+ # calculate affinities with smooth transitions between basis species
m2 <- mosaic(bases, pH = pH)
# these species have no S so the results should be similar,
- # 20190121 except for a negative free energy of mixing (positive affinity)
- expect_true(all(m2$A.species$values[[1]] - m1$A.species$values[[1]] > 0))
- # the differences increase, then decrease
- expect_equal(unique(sign(diff(as.numeric(m2$A.species$values[[1]] - m1$A.species$values[[1]])))), c(1, -1))
+ expect_equivalent(m2$A.species$values[[1]], m1$A.species$values[[1]])
# now with S-bearing species ...
species(c("pyrrhotite", "pyrite"))
m3 <- mosaic(bases, pH = pH, blend = FALSE)
@@ -49,8 +46,71 @@
# the results are different ...
expect_equal(sapply(m3$A.species$values, "[", 13), sapply(m4$A.species$values, "[", 13), tol=1e-1)
# but more similar at extreme pH values
- expect_equal(sapply(m3$A.species$values, "[", 1), sapply(m4$A.species$values, "[", 1), tol=1e-6)
- expect_equal(sapply(m3$A.species$values, "[", 29), sapply(m4$A.species$values, "[", 29), tol=1e-10)
+ expect_equal(sapply(m3$A.species$values, "[", 1), sapply(m4$A.species$values, "[", 1), tol=1e-7)
+ expect_equal(sapply(m3$A.species$values, "[", 29), sapply(m4$A.species$values, "[", 29), tol=1e-13)
})
+test_that("mosaic() - equilibrate() produces equilibrium activities", {
+ # test added 20190505, based on an calculation sent by Kirt Robinson
+ basis(c("CO2", "NH3", "O2", "H2O", "H+"))
+ species(c("acetamide", "acetic acid", "acetate"))
+ m <- mosaic(c("NH3", "NH4+"), pH = c(0, 14))
+ e <- equilibrate(m$A.species)
+ # calculate logK for form acetamide from predominant species at low pH
+ s1 <- subcrt(c("acetic acid", "NH4+", "acetamide", "water", "H+"), c(-1, -1, 1, 1, 1), T = 25)
+ logK1 <- s1$out$logK
+ # values of activities
+ loga_acetic <- e$loga.equil[[2]]
+ loga_NH4 <- m$e.bases[[1]]$loga.equil[[2]]
+ loga_acetamide <- e$loga.equil[[1]]
+ loga_H2O <- m$e.bases[[1]]$basis$logact[[4]]
+ loga_Hplus <- - m$e.bases[[1]]$vals$pH
+ logQ1 <- - loga_acetic - loga_NH4 + loga_acetamide + loga_H2O + loga_Hplus
+ A1 <- logQ1 - logK1
+ ## in CHNOSZ versions before 1.3.2-5 (20190505), the affinity was zero at the pH extremes,
+ ## but peaked with a value of 0.3 (log10(2)) at pH 9.2 (equal activities of NH3 and NH4+)
+ #plot(m$e.bases[[1]]$vals$pH, A1, type = "l")
+ #title(main = describe.reaction(s1$reaction))
+ expect_equivalent(A1, rep(0, length(A1)))
+})
+
+test_that("mosaic() - solubility() produces equilibrium activities", {
+ # test added 20190505, simplified from demo/contour.R with varying pH at constant logfO2
+ # define temperature and pressure
+ T <- 250
+ P <- 500
+ # set up system
+ basis(c("Au", "Cl-", "H2S", "H2O", "oxygen", "H+"))
+ species(c("Au(HS)2-", "AuHS", "AuOH", "AuCl2-"))
+ # this get us close to total S = 0.01 m
+ basis("H2S", -2)
+ # set a low logfO2 to get into H2S - HS- fields
+ basis("O2", -40)
+ # calculate solution composition for 1 mol/kg NaCl
+ NaCl <- NaCl(T = T, P = P, m_tot = 1)
+ basis("Cl-", log10(NaCl$m_Cl))
+ # calculate affinity with changing basis species
+ bases <- c("H2S", "HS-", "HSO4-", "SO4-2")
+ m <- mosaic(bases, pH = c(2, 10), T = 250, P = 500, IS = NaCl$IS)
+ # calculate solubility
+ s <- solubility(m$A.species)
+ # calculate logK to form Au(HS)2- in H2S stability region
+ # include IS here to compute adjusted logK
+ s1 <- subcrt(c("Au", "H2S", "oxygen", "Au(HS)2-", "H2O", "H+"), c(-1, -2, -0.25, 1, 0.5, 1), T = T, P = P, IS = NaCl$IS)
+ logK1 <- s1$out$logK
+ # calculate logQ with the given or computed activities
+ loga_Au <- m$A.bases$basis$logact[[1]]
+ loga_H2S <- m$e.bases[[1]]$loga.equil[[1]]
+ logf_O2 <- m$A.bases$basis$logact[[5]]
+ loga_AuHS2minus <- s$loga.equil[[1]]
+ loga_H2O <- m$A.bases$basis$logact[[4]]
+ loga_Hplus <- - m$A.bases$vals$pH
+ logQ1 <- - 1 * loga_Au - 2 * loga_H2S - 0.25 * logf_O2 + 1 * loga_AuHS2minus + 0.5 * loga_H2O + 1 * loga_Hplus
+ # calculate affinity - should be zero!
+ A1 <- logQ1 - logK1
+ #plot(m$A.bases$vals$pH, A1, type = "l")
+ #title(main = describe.reaction(s1$reaction))
+ expect_equivalent(A1, rep(0, length(A1)))
+})
+
# TODO: test that basis specifications can be exchanged between bases and bases2 without altering output
More information about the CHNOSZ-commits
mailing list