[CHNOSZ-commits] r344 - in pkg/CHNOSZ: . R demo inst man tests/testthat
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Nov 7 03:18:54 CET 2018
Author: jedick
Date: 2018-11-07 03:18:53 +0100 (Wed, 07 Nov 2018)
New Revision: 344
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/R/examples.R
pkg/CHNOSZ/R/solubility.R
pkg/CHNOSZ/demo/gold.R
pkg/CHNOSZ/demo/solubility.R
pkg/CHNOSZ/inst/NEWS
pkg/CHNOSZ/man/NaCl.Rd
pkg/CHNOSZ/man/equilibrate.Rd
pkg/CHNOSZ/man/nonideal.Rd
pkg/CHNOSZ/man/solubility.Rd
pkg/CHNOSZ/tests/testthat/test-solubility.R
Log:
solubility(): calculate activity directly from affinity (no equilibrate step)
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2018-11-06 05:24:26 UTC (rev 343)
+++ pkg/CHNOSZ/DESCRIPTION 2018-11-07 02:18:53 UTC (rev 344)
@@ -1,6 +1,6 @@
Date: 2018-11-06
Package: CHNOSZ
-Version: 1.1.3-51
+Version: 1.1.3-52
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-06 05:24:26 UTC (rev 343)
+++ pkg/CHNOSZ/R/examples.R 2018-11-07 02:18:53 UTC (rev 344)
@@ -9,8 +9,8 @@
"util.array", "util.blast", "util.data", "util.expression",
"util.fasta", "util.formula", "util.matrix", "util.misc", "util.seq", "util.units",
"util.water", "taxonomy", "info", "protein.info", "hkf", "water", "IAPWS95", "subcrt",
- "makeup", "basis", "swap.basis", "species", "affinity", "equil.boltzmann",
- "diagram", "buffer", "nonideal", "add.protein", "protein", "ionize.aa", "yeast.aa",
+ "makeup", "basis", "swap.basis", "species", "affinity", "solubility", "equilibrate",
+ "diagram", "buffer", "nonideal", "NaCl", "add.protein", "protein", "ionize.aa", "yeast.aa",
"objective", "revisit", "EOSregress", "wjd")
plot.it <- FALSE
if(is.character(save.png))
Modified: pkg/CHNOSZ/R/solubility.R
===================================================================
--- pkg/CHNOSZ/R/solubility.R 2018-11-06 05:24:26 UTC (rev 343)
+++ pkg/CHNOSZ/R/solubility.R 2018-11-07 02:18:53 UTC (rev 344)
@@ -1,30 +1,55 @@
# solubility.R: vectorized solubility calculations without uniroot
# 20181031 jmd
+# 20181106 work on the output from affinity(); no "equilibrate()" needed!
-solubility <- function(eout, exp = 1) {
- # exp = 1: e.g. dissolution of CO2
- # exp = 2: e.g. dissolution (and dissociation) of CaCO3
+## if this file is interactively sourced, the following are also needed to provide unexported functions:
+#source("equilibrate.R")
+#source("util.misc.R")
- # bookkeeping: track any single species
- itrack <- 1
- # the log activity used to calculate the affinity
- loga.species.track <- eout$species$logact[itrack]
- # the affinities at the starting loga.balance
- A.track <- eout$values[[itrack]]
- # the loga.equil at the starting loga.balance
- loga.equil.track <- eout$loga.equil[[itrack]]
+solubility <- function(aout, balance=NULL, split=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)
- # subjunctive: what would the affinities be if the
- # activity of the tracked species was set to loga.equil?
- A.whatif <- loga.species.track + A.track - loga.equil.track
+ # get the balancing coefficients
+ bout <- balance(aout, balance)
+ n.balance <- bout$n.balance
+ balance <- bout$balance
- # predictive: assuming the species distribution doesn't change,
- # what is the total loga that gives zero affinity?
- # TODO: modify this according to stoichiometry (species with > 1 of the balanced basis species)
- loga.total <- (eout$loga.balance + A.whatif) / exp
- message("solubility: calculated logarithm of total activity of ", eout$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 <- Reduce("+", a.balance)
+ log10(a.balance)
+ }
+ loga.balance <- logabfun(loga.equil, bout$n.balance)
- # use the predicted loga.total to re-calculate activities of species
- aout <- eout[1:which(names(eout)=="values")]
- equilibrate(aout, loga.balance = loga.total)
+ # 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))
+ }
+
+ # make the output
+ # (we don't deal with normalized formulas yet, so for now m.balance==n.balance)
+ c(aout, list(balance=bout$balance, m.balance=bout$n.balance, n.balance=bout$n.balance,
+ loga.balance=loga.balance, Astar=loga.equil, loga.equil=loga.equil))
}
Modified: pkg/CHNOSZ/demo/gold.R
===================================================================
--- pkg/CHNOSZ/demo/gold.R 2018-11-06 05:24:26 UTC (rev 343)
+++ pkg/CHNOSZ/demo/gold.R 2018-11-07 02:18:53 UTC (rev 344)
@@ -80,11 +80,10 @@
# apply PPM buffer for fO2 and aH2S
basis("O2", "PPM")
basis("H2S", "PPM")
- # calculate affinity, equilibrate, solubility
+ # calculate affinity and solubility
# (set IS = 0 for diagram to show "log m" instead of "log a")
a <- affinity(pH = c(3, 8), T = 300, P = 1000, IS = 0)
- e <- equilibrate(a)
- s <- solubility(e)
+ s <- solubility(a)
# make diagram and show total log molality
diagram(s, ylim = c(-10, -5), col = col, lwd = 2, lty = 1)
diagram(s, add = TRUE, type = "loga.balance", lwd = 3, lty = 2)
@@ -106,11 +105,10 @@
# apply PPM buffer for fO2 and aH2S
basis("O2", "PPM")
basis("H2S", "PPM")
- # calculate affinity, equilibrate, solubility
+ # calculate affinity and solubility
# (set IS = 0 for diagram to show "log m" instead of "log a")
a <- affinity(pH = c(3, 8), T = 450, P = 1000, IS = 0)
- e <- equilibrate(a)
- s <- solubility(e)
+ s <- solubility(a)
# make diagram and show total log molality
diagram(s, ylim = c(-8, -3), col = col, lwd = 2, lty = 1)
diagram(s, add = TRUE, type = "loga.balance", lwd = 3, lty = 2)
@@ -141,10 +139,9 @@
# assuming complete dissociation of 0.5 mol/kg KCl
gam_K <- 10^subcrt("K+", T = seq(150, 550, 10), P = 1000, IS=NaCl$IS)$out$`K+`$loggam
a_K <- 0.5 * gam_K
- # calculate affinity, equilibrate, solubility
+ # calculate affinity and solubility
a <- affinity(T = seq(150, 550, 10), `Cl-` = log10(a_Cl), `K+` = log10(a_K), P = 1000, IS = NaCl$IS)
- e <- equilibrate(a)
- s <- solubility(e)
+ s <- solubility(a)
# make diagram and show total log molality
diagram(s, ylim = c(-10, -4), col = col, lwd = 2, lty = 1)
diagram(s, add = TRUE, type = "loga.balance", lwd = 3, lty = 2)
@@ -176,10 +173,9 @@
# assuming complete dissociation of 0.5 mol/kg KCl
gam_K <- 10^subcrt("K+", T = seq(150, 550, 10), P = 1000, IS=NaCl$IS)$out$`K+`$loggam
a_K <- 0.5 * gam_K
- # calculate affinity, equilibrate, solubility
+ # calculate affinity and solubility
a <- affinity(T = seq(150, 550, 10), `Cl-` = log10(a_Cl), `K+` = log10(a_K), P = 1000, IS = NaCl$IS)
- e <- equilibrate(a)
- s <- solubility(e)
+ s <- solubility(a)
# make diagram and show total log molality
diagram(s, ylim = c(-10, -2), col = col, lwd = 2, lty = 1)
diagram(s, add = TRUE, type = "loga.balance", lwd = 3, lty = 2)
Modified: pkg/CHNOSZ/demo/solubility.R
===================================================================
--- pkg/CHNOSZ/demo/solubility.R 2018-11-06 05:24:26 UTC (rev 343)
+++ pkg/CHNOSZ/demo/solubility.R 2018-11-07 02:18:53 UTC (rev 344)
@@ -25,8 +25,7 @@
basis("CO2", -3.5)
species(c("CO2", "HCO3-", "CO3-2"))
a <- affinity(pH = c(pH, res), T = T1, IS = IS)
-e <- equilibrate(a)
-s <- solubility(e)
+s <- solubility(a)
# first plot total activity line
diagram(s, ylim = c(-10, 4), type = "loga.balance", lwd = 4, col = "green2")
# add activities of species
@@ -39,11 +38,12 @@
title(main = substitute("Solubility of"~what~"at"~T~degree*"C",
list(what = expr.species("CO2"), T = T1)), line = 1.6)
mtext("cf. Fig. 4.5 of Stumm and Morgan, 1996")
+# check the endpoints
+stopifnot(round(s$loga.balance[c(1, res)])==c(-5, 6))
# CO2 T-pH plot
a <- affinity(pH = c(pH, res), T = c(T, res), IS = IS)
-e <- equilibrate(a)
-s <- solubility(e)
+s <- solubility(a)
diagram(s, type = "loga.balance")
title(main = substitute("Solubility of"~what, list(what = expr.species("CO2"))))
@@ -51,8 +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)
-e <- equilibrate(a)
-s <- solubility(e, exp = 2)
+s <- solubility(a, split = TRUE)
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),
@@ -60,10 +59,11 @@
title(main = substitute("Solubility of"~what~"at"~T~degree*"C",
list(what = "calcite", T = T1)), line = 1.6)
mtext("cf. Fig. 4A of Manning et al., 2013")
+# check the endpoints
+stopifnot(round(s$loga.balance[c(1, res)])==c(4, -4))
# calcite T-pH plot
a <- affinity(pH = c(pH, res), T = c(T, res), IS = IS)
-e <- equilibrate(a)
-s <- solubility(e, exp = 2)
+s <- solubility(a, split = TRUE)
diagram(s, type = "loga.balance")
title(main = "Solubility of calcite", font.main = 1)
Modified: pkg/CHNOSZ/inst/NEWS
===================================================================
--- pkg/CHNOSZ/inst/NEWS 2018-11-06 05:24:26 UTC (rev 343)
+++ pkg/CHNOSZ/inst/NEWS 2018-11-07 02:18:53 UTC (rev 344)
@@ -1,10 +1,10 @@
-CHANGES IN CHNOSZ 1.1.3-51 (2018-11-06)
+CHANGES IN CHNOSZ 1.1.3-52 (2018-11-07)
---------------------------------------
NEW FEATURES
-- Add solubility(). Run this after equilibrate() to calculate the
- solubility (loga.balance) of the balanced basis species.
+- Add solubility(). Run this after affinity() to calculate the
+ solubility of the conserved basis species.
- Revise demo/solubility.R to show solubility calculations for CO2(gas)
and calcite as a function of T and pH.
Modified: pkg/CHNOSZ/man/NaCl.Rd
===================================================================
--- pkg/CHNOSZ/man/NaCl.Rd 2018-11-06 05:24:26 UTC (rev 343)
+++ pkg/CHNOSZ/man/NaCl.Rd 2018-11-07 02:18:53 UTC (rev 344)
@@ -65,7 +65,7 @@
gam_Cl.calc[i, ] <- NaCl.out$gam_Cl
}
# plot ionic strength from HCh and NaCl() as points and lines
-par(mfrow=c(2, 1))
+opar <- par(mfrow=c(2, 1))
col <- c("black", "red", "orange")
plot(c(1,6), c(0,6), xlab="NaCl (mol/kg)", ylab=axis.label("IS"), type="n")
for(i in 1:3) {
@@ -88,10 +88,11 @@
}
# we should be fairly close
stopifnot(maxdiff(unlist(gam_Cl.calc[seq(1,11,2), ]), unlist(gam_Cl.HCh)) < 0.033)
+par(opar)
}
\references{
-Shvarov, Y. and Bastrakov, E. (1999) HCh: A software package for geochemical equilibrium modelling. User's Guide. \emph{Australian Geological Survey Organisation} \bold{1999/25}.
+Shvarov, Y. and Bastrakov, E. (1999) HCh: A software package for geochemical equilibrium modelling. User's Guide. \emph{Australian Geological Survey Organisation} \bold{1999/25}. \url{http://pid.geoscience.gov.au/dataset/ga/25473}
}
\concept{Extended workflow}
Modified: pkg/CHNOSZ/man/equilibrate.Rd
===================================================================
--- pkg/CHNOSZ/man/equilibrate.Rd 2018-11-06 05:24:26 UTC (rev 343)
+++ pkg/CHNOSZ/man/equilibrate.Rd 2018-11-07 02:18:53 UTC (rev 344)
@@ -83,6 +83,13 @@
}
+\section{Warning}{
+Despite its name, this function does not generally produce a complete equilibrium.
+It returns activities of species such that the affinities of formation reactions are equal to each other (and transformations between species have zero affinity); this is a type of metastable equilibrium.
+Although they are equal to each other, the affinities are not necessarily equal to zero.
+Use \code{solubility} to find complete equilibrium, where the affinities of the formation reactions become zero.
+}
+
\seealso{
\code{\link{diagram}} has examples of using \code{equilibrate} to make equilibrium activity diagrams. \code{\link{revisit}} can be used to perform further analysis of the equilibrium activities.
\code{\link{palply}} is used by both \code{equil.reaction} and \code{equil.boltzmann} to parallelize intensive parts of the calculations.
Modified: pkg/CHNOSZ/man/nonideal.Rd
===================================================================
--- pkg/CHNOSZ/man/nonideal.Rd 2018-11-06 05:24:26 UTC (rev 343)
+++ pkg/CHNOSZ/man/nonideal.Rd 2018-11-07 02:18:53 UTC (rev 344)
@@ -2,7 +2,7 @@
\name{nonideal}
\alias{nonideal}
\alias{bgamma}
-\title{Activity coefficients of aqueous species}
+\title{Activity Coefficients of Aqueous Species}
\description{
Calculate activity coefficients and adjusted (transformed) molal properties of aqueous species.
}
@@ -39,7 +39,7 @@
Calculations for the proton (\Hplus) and electron (\eminus) are skipped by default; this makes sense if you are setting the pH, i.e. activity of \Hplus, to some value.
To apply the calculations to H+ and/or e-, change \code{\link{thermo}$opt$ideal.H} or \code{ideal.e} to FALSE.
-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).
+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) 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}, \code{B_DH}, and \code{m_star} given in the arguments.
@@ -194,13 +194,11 @@
Helgeson, H. C., Kirkham, D. H. and Flowers, G. C. (1981) Theoretical prediction of the thermodynamic behavior of aqueous electrolytes at high pressures and temperatures. IV. Calculation of activity coefficients, osmotic coefficients, and apparent molal and standard and relative partial molal properties to 600\degC and 5 Kb. \emph{Am. J. Sci.} \bold{281}, 1249--1516. \url{https://doi.org/10.2475/ajs.281.10.1249}
-Hückel, E. (1925). The theory of concentrated, aqueous solutions of strong electrolytes. \emph{Physikalische Zeitschrift} \bold{26}, 93--147.
-
Manning, C. E. (2013) Thermodynamic modeling of fluid-rock interaction at mid-crustal to upper-mantle conditions. \emph{Rev. Mineral. Geochem.} \bold{76}, 135--164. \url{https://doi.org/10.2138/rmg.2013.76.5}
Manning, C. E., Shock, E. L. and Sverjensky, D. A. (2013) The chemistry of carbon in aqueous fluids at crustal and upper-mantle conditions: Experimental and theoretical constraints. \emph{Rev. Mineral. Geochem.} \bold{75}, 109--148. \url{https://doi.org/10.2138/rmg.2013.75.5}
-Shvarov, Y. and Bastrakov, E. (1999) HCh: A software package for geochemical equilibrium modelling. User's Guide. \emph{Australian Geological Survey Organisation} \bold{1999/25}.
+Shvarov, Y. and Bastrakov, E. (1999) HCh: A software package for geochemical equilibrium modelling. User's Guide. \emph{Australian Geological Survey Organisation} \bold{1999/25}. \url{http://pid.geoscience.gov.au/dataset/ga/25473}
}
\concept{Thermodynamic calculations}
Modified: pkg/CHNOSZ/man/solubility.Rd
===================================================================
--- pkg/CHNOSZ/man/solubility.Rd 2018-11-06 05:24:26 UTC (rev 343)
+++ pkg/CHNOSZ/man/solubility.Rd 2018-11-07 02:18:53 UTC (rev 344)
@@ -7,44 +7,46 @@
}
\usage{
- solubility(eout, exp = 1)
+ solubility(aout, balance = NULL, split = FALSE)
}
\arguments{
- \item{eout}{list, output from \code{\link{equilibrate}}}
- \item{exp}{numeric, exponent characterizing the stoichiometry of the dissociation reaction}
+ \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?}
}
\details{
-Use this function to calculate solubilities of minerals (such as CaCO\s3) or gases (such as CO\s2).
-Start by using \code{\link{equilibrate}} to calculate equilibrium chemical activities of species given a constant value of \code{loga.balance} (the logarithm of total activity of the balanced basis species).
-Note that this produces affinities of formation reactions of species that are equal to each other, but are generally not equal to zero.
-\code{solubility} adjusts \code{loga.balance} such that the affinities of the formation reaction become zero.
-This corresponds to \dQuote{true} equilibrium for a solution in contact with the balanced basis species - i.e. the solubility of that species.
+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.
-Normally, the balance is automatically identified as the first basis species that is present in all of the species.
-If that is not adequate, it can be explicitly set via the \samp{balance} setting in \code{equilibrate}.
+\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 value of \code{exp} should be changed when calculating solubility of species that dissociate (not just dissolve).
-For example, set \code{exp} to 2 for calculating the solubility of calcite (CaCO\s3).
-
-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"}.
+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.
}
\section{Warning}{
This function has not been tested for systems that may form dimers or higher-order complexes (such as Au\s{2}S\s{2}\S{2-}).
-The lower figures in \code{demo("gold")} are incomplete, as they do not account for other possible reactions not involving Au, particularly the decrease of Cl\S{-} concentration owing to the rising stability of the NaCl\s{(aq)} complex at high temperature.
Except for relatively simple systems, even after careful refinement, the results from CHNOSZ, which considers chemical activities as the independent variables, will not match the results from speciation-solubility (or Gibbs energy minimization) codes, where the system is defined by its bulk composition.
}
\seealso{
\code{demo("solubility")} adds \T-pH diagrams to the CO\s{2} and calcite example here.
\code{demo("gold")} shows solubility calculations for Au in aqueous solutions with hydroxide, chloride, and hydrosulfide complexes.
+\code{\link{equilibrate}} calculates equilibrium chemical activities of species given a constant value of \code{loga.balance} (the logarithm of total activity of the conserved basis species).
}
\examples{\dontshow{data(thermo)}
# solubility of CO2 and calcite as a function of pH
-par(mfrow = c(1, 2))
+opar <- par(mfrow = c(1, 2))
# set pH range and resolution, constant temperature and ionic strength
pH <- c(0, 14)
@@ -58,8 +60,7 @@
basis("CO2", -3.5)
species(c("CO2", "HCO3-", "CO3-2"))
a <- affinity(pH = c(pH, res), T = T, IS = IS)
-e <- equilibrate(a)
-s <- solubility(e)
+s <- solubility(a)
# first plot total activity line
diagram(s, ylim = c(-10, 4), type = "loga.balance", lwd = 4, col = "green2")
# add activities of species
@@ -77,8 +78,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)
-e <- equilibrate(a)
-s <- solubility(e, exp = 2)
+s <- solubility(a, split = TRUE)
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),
@@ -86,6 +86,8 @@
title(main = substitute("Solubility of"~what~"at"~T~degree*"C",
list(what = "calcite", T = T)))
mtext("cf. Fig. 4A of Manning et al., 2013")
+
+par(opar)
}
\references{
Modified: pkg/CHNOSZ/tests/testthat/test-solubility.R
===================================================================
--- pkg/CHNOSZ/tests/testthat/test-solubility.R 2018-11-06 05:24:26 UTC (rev 343)
+++ pkg/CHNOSZ/tests/testthat/test-solubility.R 2018-11-07 02:18:53 UTC (rev 344)
@@ -7,29 +7,38 @@
T <- 25
IS <- 0
- # start with CO2
+ ## start with CO2
basis(c("carbon dioxide", "H2O", "O2", "H+"))
# ca. atmospheric PCO2
basis("CO2", -3.5)
species(c("CO2", "HCO3-", "CO3-2"))
a <- affinity(pH = c(pH, res), T = T, IS = IS)
- e <- equilibrate(a)
- s <- solubility(e)
+ s <- solubility(a)
+ # a function to check for stable conditions (affinity = 0)
+ # do this by setting activities in species() then calculating the affintiy
+ checkfun <- function(i) {
+ logact <- sapply(s$loga.equil, "[", i)
+ species(1:3, logact)
+ basis("pH", s$vals[[1]][i])
+ affinity(T = T, IS = IS)
+ }
+ # check any 'i' here - let's just take two
+ expect_equal(max(abs(unlist(checkfun(33)$values))), 0)
+ expect_equal(max(abs(unlist(checkfun(99)$values))), 0)
- # check for stable conditions (affinity = 0)
- species(1:3, 0)
- atest <- affinity(pH = s$vals[[1]], T = T, IS = IS)
- expect_true(all(sapply(unlist(atest$values) - unlist(s$loga.equil), all.equal, 0)))
-
# now do calcite
basis(c("calcite", "Ca+2", "H2O", "O2", "H+"))
species(c("CO2", "HCO3-", "CO3-2"))
a <- affinity(pH = c(pH, res), T = T, IS = IS)
- e <- equilibrate(a)
- s <- solubility(e, exp = 2)
-
- # check for stable conditions (affinity = 0)
- species(1:3, 0)
- atest <- affinity(pH = s$vals[[1]], `Ca+2` = s$loga.balance, T = T, IS = IS)
- expect_true(all(sapply(unlist(atest$values) - unlist(s$loga.equil), all.equal, 0)))
+ s <- solubility(a, split = TRUE)
+ # here we need to also set the activity of Ca+2
+ checkfun <- function(i) {
+ logact <- sapply(s$loga.equil, "[", i)
+ species(1:3, logact)
+ basis("pH", s$vals[[1]][i])
+ basis("Ca+2", s$loga.balance[i])
+ affinity(T = T, IS = IS)
+ }
+ expect_equal(max(abs(unlist(checkfun(33)$values))), 0)
+ expect_equal(max(abs(unlist(checkfun(99)$values))), 0)
})
More information about the CHNOSZ-commits
mailing list