[CHNOSZ-commits] r357 - in pkg/CHNOSZ: . R demo inst man tests/testthat
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Jan 17 16:40:18 CET 2019
Author: jedick
Date: 2019-01-17 16:40:17 +0100 (Thu, 17 Jan 2019)
New Revision: 357
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/R/solubility.R
pkg/CHNOSZ/demo/solubility.R
pkg/CHNOSZ/inst/NEWS
pkg/CHNOSZ/man/solubility.Rd
pkg/CHNOSZ/tests/testthat/test-solubility.R
Log:
solubility(): detect dissolution reactions and preliminary implementation of find.IS
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2019-01-17 07:32:56 UTC (rev 356)
+++ pkg/CHNOSZ/DESCRIPTION 2019-01-17 15:40:17 UTC (rev 357)
@@ -1,6 +1,6 @@
Date: 2019-01-17
Package: CHNOSZ
-Version: 1.1.3-64
+Version: 1.1.3-65
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/solubility.R
===================================================================
--- pkg/CHNOSZ/R/solubility.R 2019-01-17 07:32:56 UTC (rev 356)
+++ pkg/CHNOSZ/R/solubility.R 2019-01-17 15:40:17 UTC (rev 357)
@@ -1,51 +1,97 @@
# solubility.R: vectorized solubility calculations without uniroot
# 20181031 jmd
# 20181106 work on the output from affinity(); no "equilibrate()" needed!
+# 20190117 add find.IS and test for dissociation reaction
## if this file is interactively sourced, the following are also needed to provide unexported functions:
#source("equilibrate.R")
#source("util.misc.R")
-solubility <- function(aout, balance=NULL, split=FALSE) {
+solubility <- function(aout, dissociation=NULL, find.IS=FALSE) {
## concept: the logarithms of activities of species at equilibrium are equal to
## Astar, the affinities calculated for unit activities of species
-
- ## however, the values in aout can be calculated for other than
- ## unit activities of species, so we have to take away the activites
- Astar <- function(i) aout$values[[i]] + aout$species$logact[i]
- loga.equil <- lapply(1:length(aout$values), Astar)
- ## for a dissociation (split) on a *per reaction* (not system) basis,
- ## apply the divisor here and skip the if(split){} part below
- ## (can be used to reproduce Fig. 4 of Manning et al., 2013)
- if(is.numeric(split)) loga.equil <- lapply(loga.equil, "/", split)
+ ## does the system involve a dissociation reaction?
+ if(is.null(dissociation)) {
+ # assume FALSE unless determined otherwise
+ dissociation <- FALSE
+ # react the first basis species to form the first species
+ sres <- suppressMessages(subcrt(c(rownames(aout$basis)[1], aout$species$name[1]), c(-1, 1)))
+ # note that the reaction is auto-balanced using all the basis species
+ # if the reaction involves the second basis species, we consider it to be a dissociation reaction
+ if(rownames(aout$basis)[2] %in% sres$reaction$formula) {
+ # however, if the second basis species is H2O, H+, e-, O2 (or others?), we don't have enough information, so stop
+ if(rownames(aout$basis)[2] %in% c("H2O", "H+", "e-", "O2")) {
+ print(sres$reaction)
+ stop("unsure whether this is a dissociation reaction. if it is, redefine the basis to put a product ion second")
+ }
+ dissociation <- TRUE
+ }
+ message("solubility: test for dissociation reaction returns ", dissociation)
+ } else message("solubility: argument for dissociation reaction is ", dissociation)
+ # get starting ionic strength (probably zero, but could be anything set by user)
+ IS <- aout$IS
+ if(is.null(IS)) IS <- 0
+
## to output loga.balance we need the balancing coefficients
- bout <- balance(aout, balance)
+ bout <- balance(aout)
n.balance <- bout$n.balance
balance <- bout$balance
# get logarithm of total activity of the balancing basis species
logabfun <- function(loga.equil, n.balance) {
# exponentiate, multiply by n.balance, sum, logarithm
a.equil <- mapply("^", 10, loga.equil, SIMPLIFY = FALSE)
- a.balance <- mapply("*", a.equil, n.balance, SIMPLIFY=FALSE)
+ a.balance <- mapply("*", a.equil, n.balance, SIMPLIFY = FALSE)
a.balance <- Reduce("+", a.balance)
log10(a.balance)
}
- loga.balance <- logabfun(loga.equil, bout$n.balance)
+
+ # for find.IS=TRUE, iterate to converge on ionic strength
+ niter <- 1
+ while(TRUE) {
+ ## the values in aout can be calculated for other than
+ ## unit activities of species, so we have to take away the activites
+ Astar <- function(i) aout$values[[i]] + aout$species$logact[i]
+ loga.equil <- lapply(1:length(aout$values), Astar)
- # recalculate things for a 1:1 split species (like CaCO3 = Ca+2 + CO3+2)
- if(isTRUE(split)) {
- # the multiplicity becomes the exponent in the reaction quotient
- loga.split <- loga.balance / 2
- # the contribution to affinity
- Asplit <- lapply(n.balance, "*", loga.split)
- # adjust the affinity and get new equilibrium activities
- aout$values <- mapply("-", aout$values, Asplit, SIMPLIFY=FALSE)
- loga.equil <- lapply(1:length(aout$values), Astar)
- # check that the new loga.balance == loga.split; this might not work for non-1:1 species
- loga.balance <- logabfun(loga.equil, n.balance)
- stopifnot(all.equal(loga.balance, loga.split))
+ ## for a dissociation on a *per reaction* (not system) basis,
+ ## apply the divisor here and skip the if(dissociation){} part below
+ ## (can be used to reproduce Fig. 4 of Manning et al., 2013)
+ if(is.numeric(dissociation)) loga.equil <- lapply(loga.equil, "/", dissociation)
+
+ loga.balance <- logabfun(loga.equil, bout$n.balance)
+
+ # recalculate things for a dissociation reaction (like CaCO3 = Ca+2 + CO3+2)
+ if(isTRUE(dissociation)) {
+ # the multiplicity becomes the exponent in the reaction quotient
+ loga.split <- loga.balance / 2
+ # the contribution to affinity
+ Asplit <- lapply(n.balance, "*", loga.split)
+ # adjust the affinity and get new equilibrium activities
+ aout$values <- mapply("-", aout$values, Asplit, SIMPLIFY=FALSE)
+ loga.equil <- lapply(1:length(aout$values), Astar)
+ # check that the new loga.balance == loga.split
+ # TODO: does this work for non-1:1 species?
+ loga.balance <- logabfun(loga.equil, n.balance)
+ stopifnot(all.equal(loga.balance, loga.split))
+ }
+
+ # get the old IS
+ IS.old <- rep(IS, length.out = length(aout$values[[1]]))
+ # calculate the ionic strength (assuming no ion pairing)
+ # i.e. for Sr+2 and SO4-2
+ # TODO: determine the correct values for the actual species in the system
+ mol.cation <- mol.anion <- 10^loga.balance
+ Z <- 2
+ IS <- (mol.cation * Z^2 + mol.anion * Z^2) / 2
+ # report the current ionic strength
+ if(find.IS) message("solubility: (iteration ", niter, ") ionic strength range is ", paste(round(range(IS), 4), collapse=" "))
+ # stop iterating if we reached the tolerance (or find.IS=FALSE)
+ if(!find.IS | all(IS - IS.old < 1e-4)) break
+ # recalculate the affinity using the new IS
+ aout <- suppressMessages(do.call(aout$fun, list(aout, IS = IS)))
+ niter <- niter + 1
}
# make the output
Modified: pkg/CHNOSZ/demo/solubility.R
===================================================================
--- pkg/CHNOSZ/demo/solubility.R 2019-01-17 07:32:56 UTC (rev 356)
+++ pkg/CHNOSZ/demo/solubility.R 2019-01-17 15:40:17 UTC (rev 357)
@@ -51,7 +51,7 @@
basis(c("calcite", "Ca+2", "H2O", "O2", "H+"))
species(c("CO2", "HCO3-", "CO3-2"))
a <- affinity(pH = c(pH, res), T = T1, IS = IS)
-s <- solubility(a, split = TRUE)
+s <- solubility(a)
diagram(s, ylim = c(-10, 4), type = "loga.balance", lwd = 4, col = "green2")
diagram(s, add = TRUE, dy = 1)
legend("topright", lty = c(1, 1:3), lwd = c(4, 2, 2, 2),
@@ -64,6 +64,6 @@
# calcite T-pH plot
a <- affinity(pH = c(pH, res), T = c(T, res), IS = IS)
-s <- solubility(a, split = TRUE)
+s <- solubility(a)
diagram(s, type = "loga.balance")
title(main = "Solubility of calcite", font.main = 1)
Modified: pkg/CHNOSZ/inst/NEWS
===================================================================
--- pkg/CHNOSZ/inst/NEWS 2019-01-17 07:32:56 UTC (rev 356)
+++ pkg/CHNOSZ/inst/NEWS 2019-01-17 15:40:17 UTC (rev 357)
@@ -1,10 +1,14 @@
-CHANGES IN CHNOSZ 1.1.3-64 (2019-01-17)
+CHANGES IN CHNOSZ 1.1.3-65 (2019-01-17)
---------------------------------------
NEW FEATURES
- Add solubility(). Run this after affinity() to calculate the
- solubility of the conserved basis species.
+ solubility of a solid or gas defined as the conserved basis species,
+ which is involved in the formation of one or more dissolved species.
+ Features in development include automatic detection of dissociation
+ reactions and finding the final ionic strength for dissolution of a
+ mineral into pure water.
- Revise demo/solubility.R to show solubility calculations for CO2(gas)
and calcite as a function of T and pH.
Modified: pkg/CHNOSZ/man/solubility.Rd
===================================================================
--- pkg/CHNOSZ/man/solubility.Rd 2019-01-17 07:32:56 UTC (rev 356)
+++ pkg/CHNOSZ/man/solubility.Rd 2019-01-17 15:40:17 UTC (rev 357)
@@ -7,29 +7,48 @@
}
\usage{
- solubility(aout, balance = NULL, split = FALSE)
+ solubility(aout, dissociation = NULL, find.IS = FALSE)
}
\arguments{
\item{aout}{list, output from \code{\link{affinity}}}
- \item{balance}{character, basis species to conserve in reactions}
- \item{split}{logical, does the mineral undergo a dissociation reaction?}
+ \item{dissociation}{logical, does the mineral undergo a dissociation reaction?}
+ \item{find.IS}{logical, find the equilibrium ionic strength by iteration?}
}
\details{
This function performs a simple task: from the values of \code{\link{affinity}} of formation reactions of species at given activity, it works backward to find the activities of species that make the affinities zero.
This corresponds to complete equilibrium with all of the basis species.
-For solubility calculations, the basis species should be set up so that formation reactions of species are balanced on the thing being dissolved (a mineral such as CaCO\s3 or gas such as CO\s2).
-This is usually identified as the first basis species, but if there is any ambiguity, use \code{balance} to indicate the conserved basis species.
+Usually, the basis species should be set up so that the first basis species represents the substance being dissolved (a mineral such as CaCO\s3 or gas such as CO\s2).
+Internally, this is treated as the conserved basis species, so it must be present in all of the formation reactions of the species.
+It is also possible to set the conserved basis species as other than the first one (see \code{demo(gold)}), but this implies that dissociation reactions are not occurring (see below).
-\code{split} should be set when calculating the solubility of something that dissociates (not just dissolves).
-For example, to calculate the solubility of calcite (CaCO\s3), each form of carbonate is represented by a different reaction, but all of these reactions also release calcium ions.
-The equilibrium calculation must take account of the \emph{total} activity of the shared ion (Ca\S{+2}); naturally, that value was unknown for the calculation of \code{\link{affinity}}.
-The solution is accomplished by setting \code{split} to TRUE to recalculate the affinities (working backward, as if the split didn't occur).
-Then, the resulting activities correspond to equilibrium considering the system-wide activity of Ca\S{+2}.
-A \emph{not recommended} alternative is to set \code{split} to a numeric value (probably 2) to calculate activities on a per-reaction basis, where each reaction has its own activity of Ca\S{+2}.
-That does not give a complete equilibrium, but may be required to reproduce some published diagrams.
+The \code{\link{species}} should be defined to represent one set of ions (anions or cations or their complexes) formed in solution, all involving the conserved basis species.
+For a dissociation reaction, the second basis species should be used to represent the counterion (cation or anion).
+The function perfoms some additional steps to calculate the solubility of something that dissociates (not just dissolves).
+For example, the dissolution of calcite (CaCO\s3), involves the release of both calcium ions and different forms of carbonate in solution, depending on the pH.
+The equilibrium calculation must take account of the \emph{total} activity of the shared ion (Ca\S{+2}), which is unknown at the start of the calculation.
+The solution is accomplished by recalculating the affinities, essentially working backward from the assumption that the dissociation didn't occur.
+The resulting activities correspond to equilibrium considering the system-wide activity of Ca\S{+2}.
+
+The function attempts to automatically detect whether dissociation reactions are involved by looking at the first species.
+If the formation reaction of that species includes \emph{both} the first and second basis species, the \code{dissociation} flag is set to TRUE.
+An example reaction of this type can be found in demo(solubility): CaCO3 (first basis species) = Ca+2 (second basis species) + CO3-2 (first species).
+Note that if the conserved basis species is not the first basis species, then the automatic detection of \code{dissociation} will always return FALSE.
+Therefore, a reaction corresponding to Au (fourth basis species) + ... = ... gives \code{dissociation} = FALSE (see \code{demo(gold)}).
+This algorithm for determining whether dissociation occurs is prone to error, so \code{dissociation} can be explicitly set in the arguments.
+A \emph{not recommended} alternative is to set \code{dissociation} to a numeric value corresponding to the stoichiometry of released species (i.e. 2 for a 1:1 electrolyte).
+This setting indicates to calculate activities on a per-reaction basis, where each reaction has its own (independent) activity of Ca\S{+2}.
+That does not give a complete equilibrium in the system, but may be required to reproduce some published diagrams.
+
+Note that other variables (pH, ionic strength, activities of other basis species) should be defined in the preceding call to \code{\link{affinity}}.
+However, for dissolving a substance in pure water, \code{find.IS} can be set to TRUE to determine the final ionic strength.
+This works by calculating the ionic strength from the equilibrium solubility calculation, then re-running \code{affinity} with those values.
+Note that for dissociation reactions, the ionic strength is calculated from both the ions present in the species definition and the counter ion, which should be the second basis species.
+TODO: Speciation of counterions (e.g. ionized forms of carbonate or sulfate) can be accomplished by using the \code{\link{mosaic}} function instead of \code{affinity}.
+The calculation is iterated until the ionic strength deviation at every point is lower than a preset tolerance (1e-4).
+
The output of \code{solubility} has the same format as that of \code{equilibrate}, and can be used by \code{\link{diagram}} with \code{type = "loga.balance"} to plot the solubilities, or with \code{type = NULL} to plot the activities of species.
}
@@ -54,7 +73,7 @@
T <- 25
IS <- 0
-# start with CO2
+# start with CO2 (not a dissociation reaction)
basis(c("carbon dioxide", "H2O", "O2", "H+"))
# ca. atmospheric PCO2
basis("CO2", -3.5)
@@ -74,11 +93,11 @@
list(what = expr.species("CO2"), T = T)), line = 1.5)
mtext("cf. Fig. 4.5 of Stumm and Morgan, 1996")
-# now do calcite
+# now do calcite (a dissociation reaction)
basis(c("calcite", "Ca+2", "H2O", "O2", "H+"))
species(c("CO2", "HCO3-", "CO3-2"))
a <- affinity(pH = c(pH, res), T = T, IS = IS)
-s <- solubility(a, split = TRUE)
+s <- solubility(a)
diagram(s, ylim = c(-10, 4), type = "loga.balance", lwd = 4, col = "green2")
diagram(s, add = TRUE, dy = 1)
legend("topright", lty = c(1, 1:3), lwd = c(4, 2, 2, 2),
Modified: pkg/CHNOSZ/tests/testthat/test-solubility.R
===================================================================
--- pkg/CHNOSZ/tests/testthat/test-solubility.R 2019-01-17 07:32:56 UTC (rev 356)
+++ pkg/CHNOSZ/tests/testthat/test-solubility.R 2019-01-17 15:40:17 UTC (rev 357)
@@ -30,7 +30,7 @@
basis(c("calcite", "Ca+2", "H2O", "O2", "H+"))
species(c("CO2", "HCO3-", "CO3-2"))
a <- affinity(pH = c(pH, res), T = T, IS = IS)
- s <- solubility(a, split = TRUE)
+ s <- solubility(a)
# here we need to also set the activity of Ca+2
checkfun <- function(i) {
logact <- sapply(s$loga.equil, "[", i)
@@ -42,3 +42,11 @@
expect_equal(max(abs(unlist(checkfun(33)$values))), 0)
expect_equal(max(abs(unlist(checkfun(99)$values))), 0)
})
+
+test_that("solubility() catches some error conditions", {
+ # demo(solubility) has basis(c("calcite", "Ca+2", "H2O", "O2", "H+")), but what if the user puts H2O second?
+ basis(c("calcite", "H2O", "Ca+2", "O2", "H+"))
+ species(c("CO2", "HCO3-", "CO3-2"))
+ a <- affinity()
+ expect_error(solubility(a), "unsure whether this is a dissociation reaction")
+})
More information about the CHNOSZ-commits
mailing list