[CHNOSZ-commits] r246 - in pkg/CHNOSZ: . R demo inst man tests/testthat vignettes
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sun Oct 8 16:27:07 CEST 2017
Author: jedick
Date: 2017-10-08 16:27:06 +0200 (Sun, 08 Oct 2017)
New Revision: 246
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/R/diagram.R
pkg/CHNOSZ/R/equilibrate.R
pkg/CHNOSZ/R/util.plot.R
pkg/CHNOSZ/R/water.R
pkg/CHNOSZ/demo/DEW.R
pkg/CHNOSZ/demo/berman.R
pkg/CHNOSZ/inst/NEWS
pkg/CHNOSZ/man/diagram.Rd
pkg/CHNOSZ/man/equilibrate.Rd
pkg/CHNOSZ/tests/testthat/test-equilibrate.R
pkg/CHNOSZ/vignettes/anintro.Rmd
Log:
allow varying values for loga.balance
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2017-10-07 06:08:37 UTC (rev 245)
+++ pkg/CHNOSZ/DESCRIPTION 2017-10-08 14:27:06 UTC (rev 246)
@@ -1,6 +1,6 @@
-Date: 2017-10-07
+Date: 2017-10-08
Package: CHNOSZ
-Version: 1.1.0-44
+Version: 1.1.0-45
Title: Thermodynamic Calculations for Geobiochemistry
Author: Jeffrey Dick
Maintainer: Jeffrey Dick <j3ffdick at gmail.com>
Modified: pkg/CHNOSZ/R/diagram.R
===================================================================
--- pkg/CHNOSZ/R/diagram.R 2017-10-07 06:08:37 UTC (rev 245)
+++ pkg/CHNOSZ/R/diagram.R 2017-10-08 14:27:06 UTC (rev 246)
@@ -48,7 +48,7 @@
} else if(what %in% rownames(eout$basis)) {
# to calculate the loga of basis species at equilibrium
if(!missing(groups)) stop("can't plot equilibrium activities of basis species for grouped species")
- if(alpha) stop("equilibrium activities of basis species not available with alpha=TRUE")
+ if(isTRUE(alpha) | is.character(alpha)) stop("equilibrium activities of basis species not available with alpha=TRUE")
plot.loga.basis <- TRUE
} else if(what=="loga.equil" & !"loga.equil" %in% names(eout)) stop("'eout' is not the output from equil()")
else if(what!="loga.equil") stop(what, " is not a basis species or 'loga.equil'")
@@ -127,10 +127,12 @@
}
## alpha: plot fractional degree of formation
- ## scale the activities to sum=1 ... 20091017
- if(alpha) {
+ # scale the activities to sum=1 ... 20091017
+ # allow scaling by balancing component 20171008
+ if(isTRUE(alpha) | is.character(alpha)) {
# remove the logarithms
act <- lapply(plotvals, function(x) 10^x)
+ if(identical(alpha, "balance")) for(i in 1:length(act)) act[[i]] <- act[[i]] * eout$n.balance[i]
# sum the activities
sumact <- Reduce("+", act)
# divide activities by the total
Modified: pkg/CHNOSZ/R/equilibrate.R
===================================================================
--- pkg/CHNOSZ/R/equilibrate.R 2017-10-07 06:08:37 UTC (rev 245)
+++ pkg/CHNOSZ/R/equilibrate.R 2017-10-08 14:27:06 UTC (rev 246)
@@ -2,6 +2,10 @@
# functions to calculation logarithm of activity
# of species in (metastable) equilibrium
+## if this file is interactively sourced, the following are also needed to provide unexported functions:
+#source("util.misc.R")
+#source("util.character.R")
+
equilibrate <- function(aout, balance=NULL, loga.balance=NULL,
ispecies=1:length(aout$values), normalize=FALSE, as.residue=FALSE,
method=c("boltzmann", "reaction"), tol=.Machine$double.eps^0.25) {
@@ -28,14 +32,25 @@
nspecies <- length(aout$values)
## say what the balancing coefficients are
if(length(n.balance) < 100) message(paste("equilibrate: n.balance is", c2s(n.balance)))
- ## logarithm of total activity of the balance
+ ## logarithm of total activity of the balancing basis species
if(is.null(loga.balance)) {
# sum up the activities, then take absolute value
# in case n.balance is negative
sumact <- abs(sum(10^aout$species$logact * n.balance))
loga.balance <- log10(sumact)
}
- message(paste0("equilibrate: loga.balance is ", loga.balance))
+ # make loga.balance the same length as the values of affinity
+ loga.balance <- unlist(loga.balance)
+ nvalues <- length(unlist(aout$values[[1]]))
+ if(length(loga.balance) == 1) {
+ # we have a constant loga.balance
+ message(paste0("equilibrate: loga.balance is ", loga.balance))
+ loga.balance <- rep(loga.balance, nvalues)
+ } else {
+ # we are using a variable loga.balance (supplied by the user)
+ if(!identical(length(loga.balance), nvalues)) stop("length of loga.balance (", length(loga.balance), ") doesn't match the affinity values (", nvalues, ")")
+ message(paste0("equilibrate: loga.balance has same length as affinity values (", length(loga.balance), ")"))
+ }
## normalize -- normalize the molar formula by the balance coefficients
m.balance <- n.balance
isprotein <- grepl("_", as.character(aout$species$name))
@@ -94,6 +109,7 @@
Anames <- names(Astar)
# first loop: make vectors
A <- palply("", 1:length(A), function(i) as.vector(A[[i]]))
+ loga.balance <- as.vector(loga.balance)
# second loop: get the exponentiated Astars (numerators)
# need to convert /2.303RT to /RT
#A[[i]] <- exp(log(10)*Astar[[i]]/n.balance[i])/n.balance[i]
@@ -103,7 +119,7 @@
At <- A[[1]]
At[] <- 0
for(i in 1:length(A)) At <- At + A[[i]]*n.balance[i]
- # fourth loop: calculate log abundances and replace the dimensions
+ # fourth loop: calculate log abundances
A <- palply("", 1:length(A), function(i) loga.balance + log10(A[[i]]/At))
# fifth loop: replace dimensions
for(i in 1:length(A)) dim(A[[i]]) <- Astardim
@@ -139,6 +155,7 @@
Anames <- names(Astar)
# make a matrix out of the list of Astar
Astar <- list2array(lapply(Astar, c))
+ if(length(loga.balance) != nrow(Astar)) stop("length of loga.balance must be equal to the number of conditions for affinity()")
# that produces the same result (other than colnames) and is much faster than
#Astar <- sapply(Astar, c)
# also, latter has NULL nrow for length(Astar[[x]])==1
@@ -148,7 +165,7 @@
# to calculate logact(thing) from Abar for the ith condition [2]
logactfun <- function(Abar, i) Astar[i,] - Abar * n.balance
# to calculate difference between logafun and loga.balance for the ith condition
- logadiff <- function(Abar, i) loga.balance - logafun(logactfun(Abar, i))
+ logadiff <- function(Abar, i) loga.balance[i] - logafun(logactfun(Abar, i))
# to calculate a range of Abar that gives negative and positive values of logadiff for the ith condition
Abarrange <- function(i) {
# starting guess of Abar (min/max) from range of Astar / n.balance
Modified: pkg/CHNOSZ/R/util.plot.R
===================================================================
--- pkg/CHNOSZ/R/util.plot.R 2017-10-07 06:08:37 UTC (rev 245)
+++ pkg/CHNOSZ/R/util.plot.R 2017-10-08 14:27:06 UTC (rev 246)
@@ -72,8 +72,14 @@
# (i.e. redox variable is on the y axis)
# get axes, T, P, and xpoints from output of affinity() or equilibrate()
if(missing(eout)) stop("'eout' (the output of affinity(), equilibrate(), or diagram()) is missing")
+ # number of variables used in affinity()
+ nvar1 <- length(eout$vars)
+ # if these were on a transect, the actual number of variables is less
+ dim <- dim(eout$loga.equil[[1]]) # for output from equilibrate()
+ if(is.null(dim)) dim <- dim(eout$values[[1]]) # for output from affinity()
+ nvar2 <- length(dim)
# we only work on diagrams with 2 variables
- if(length(eout$vars) != 2) return(NA)
+ if(nvar1 != 2 | nvar2 != 2) return(NA)
# if needed, swap axes so T or P is on x-axis
swapped <- FALSE
if(eout$vars[2] %in% c("T", "P")) {
Modified: pkg/CHNOSZ/R/water.R
===================================================================
--- pkg/CHNOSZ/R/water.R 2017-10-07 06:08:37 UTC (rev 245)
+++ pkg/CHNOSZ/R/water.R 2017-10-08 14:27:06 UTC (rev 246)
@@ -318,7 +318,7 @@
available_properties <- c("G", "epsilon", "QBorn", "V", "rho", "beta")
if(is.null(property)) return(available_properties)
# we can't use Psat here
- if(identical(P, "Psat")) stop("Psat isn't supported in this implementation of the DEW model")
+ if(identical(P, "Psat")) stop("Psat isn't supported in this implementation of the DEW model. Try setting P to at least 1000 bar.")
# use uppercase property names (including H, S, etc., so we get them from the SUPCRT properties)
wprop <- water.SUPCRT92()
iprop <- match(property, wprop)
Modified: pkg/CHNOSZ/demo/DEW.R
===================================================================
--- pkg/CHNOSZ/demo/DEW.R 2017-10-07 06:08:37 UTC (rev 245)
+++ pkg/CHNOSZ/demo/DEW.R 2017-10-08 14:27:06 UTC (rev 246)
@@ -7,7 +7,7 @@
oldwat <- water("DEW")
#### plot 1: quartz solubility at high pressure
-## after Figure 7D of Sverjensky et al., 2014
+## after Figure 7D of Sverjensky et al., 2014a [SHA14]
## (Geochim. Cosmochim. Acta, https://doi.org/10.1016/j.gca.2013.12.019)
# load SiO2 and Si2O4 data taken from DEW spreadsheet
@@ -41,14 +41,13 @@
text(lastT+25, tail(logm, 1), Pkb, adj=0)
}
t1 <- quote("Solubility of"~alpha*"-quartz")
-t2 <- "(after Sverjensky et al., 2014)"
+t2 <- "after Sverjensky et al., 2014a"
mtitle(as.expression(c(t1, t2)))
# TODO: lines are a little low at highest P and P ...
# does the Berman, 1988 quartz data increase high-PT solubilities?
#### plot 2: correlations between non-solvation volume and HKF a1 parameter
-## after Figures 12B and 12C of Sverjensky et al., 2014
-## (Geochim. Cosmochim. Acta, https://doi.org/10.1016/j.gca.2013.12.019)
+## after Figures 12B and 12C of Sverjensky et al., 2014a [SHA14]
# load the fitted parameters for species as used by SHA14
# TODO: also use their Ca+2??
# NOTE: don't load NaCl, NH4+, or HS- here because the DEW spreadsheet lists a1 from the correlation
@@ -86,8 +85,8 @@
t2 <- quote("volume and HKF "*italic(a)[1]*" parameter")
mtitle(as.expression(c(t1, t2)))
-#### plots 3 and 4: aqueous inorganic and organic carbon species at high pressure
-## after Figure 1b of Sverjensky et al., 2014
+#### plot 3: aqueous inorganic and organic carbon species at high pressure
+## after Figure 1b of Sverjensky et al., 2014b [SSH14]
## (Nature Geoscience, https://doi.org/10.1038/NGEO2291)
# define system with loga.species = 0
@@ -103,19 +102,100 @@
legend("bottomleft", legend=dp, bty="n")
}
-# first plot: CHNOSZ default database
data(OBIGT)
-dfun()
-t1 <- quote("CHNOSZ default database"[])
-t2 <- quote("(not recommended for high"~italic(P)*")")
-mtitle(as.expression(c(t1, t2)))
-# second plot: use CO2, HCO3-, CO3-2, and methane data from DEW spreadsheet
+## (not run) make diagram using CHNOSZ default database
+#dfun()
+#t1 <- quote("CHNOSZ default database"[])
+#t2 <- quote("(not recommended for high"~italic(P)*")")
+#mtitle(as.expression(c(t1, t2)))
+# make diagram using CO2, HCO3-, CO3-2, and methane data from DEW spreadsheet
add.obigt("DEW_aq", c("CO2", "HCO3-", "CO3-2", "methane"))
dfun()
CO2quote <- quote(list(CO[2], HCO[3]^"-", CO[3]^"-2"))
DEWexpr <- substitute("DEW data for"~x, list(x=CO2quote))
mtitle(as.expression(c(DEWexpr, "and methane")))
+#### plot 4: after SSH14 Fig. 3 (added 20171008)
+# conditions:
+# T = 600, 700, 800, 900, 1000 degC
+# P = 5.0GPa (50000 bar)
+# fO2 = QFM - 2
+# pH set by jadeite + kyanite + coesite
+# dissolved carbon 0.03, 0.2, 1, 4, 20 molal
+
+T <- seq(600, 1000, 5)
+
+## use DEW model
+data(thermo)
+water("DEW")
+# add species data for DEW
+inorganics <- c("methane", "CO2", "HCO3-", "CO3-2")
+organics <- c("formic acid", "formate", "acetic acid", "acetate", "propanoic acid", "propanoate")
+
+# skip updating acetate because the new data from the DEW spreadsheet give different logK
+add.obigt("DEW", c(inorganics, organics[-4]))
+## set basis species
+b.species <- c("Fe", "CO3-2", "H2O", "O2", "SiO2", "H+")
+b.state <- c("cr", "aq", "liq", "gas", "aq", "aq")
+basis(b.species, b.state)
+# for the time being we use a constant pH
+basis("H+", -4)
+
+## define a QFM buffer using Berman minerals
+mod.buffer("QFM_Berman", c("quartz", "fayalite", "magnetite"), "cr_Berman", 0)
+
+## calculate fO2 in QFM buffer minus 2
+basis("O2", "QFM_Berman")
+a <- affinity(T=T, P=50000, return.buffer=TRUE)
+QFM_2 <- a$O2 - 2
+
+## now that we have QFM-2, remove QFM buffer for calculations below
+basis("O2", 0)
+
+## add species
+species(c(inorganics, organics))
+
+## calculate affinities on the T-logfO2 transect
+a <- affinity(T=T, O2=QFM_2, P=50000)
+
+## use the total carbon molality as an approximation of total activity
+molC <- splinefun(seq(600, 1000, 100), c(0.03, 0.2, 1, 4, 20))(T)
+loga.C <- log10(molC)
+
+## calculate metastable equilibrium activities
+e <- equilibrate(a, loga.balance=loga.C)
+
+## make the diagram; don't plot names of low-abundance species
+names <- c(inorganics, organics)
+names[c(4, 5, 7, 9)] <- ""
+col <- rep("black", length(names))
+col[c(1, 3, 6, 8, 10)] <- c("red", "darkgreen", "purple", "orange", "navyblue")
+diagram(e, alpha="balance", names=names, col=col, ylim=c(0, 0.8))
+
+## add legend and title
+ltxt1 <- "P = 50000 bar"
+ltxt2 <- substitute(logfO2=="QFM-2", list(logfO2=axis.label("O2")))
+ltxt3 <- "pH = 4"
+legend("left", legend=as.expression(c(ltxt1, ltxt2, ltxt3)), bty="n")
+t1 <- "Aqueous carbon speciation"
+t2 <- "after Sverjensky et al., 2014b"
+mtitle(c(t1, t2))
+
+### additional checks
+# check that we're within 0.1 of the QFM-2 values used by SSH14
+stopifnot(maxdiff(QFM_2[T %% 100 == 0], c(-17.0, -14.5, -12.5, -10.8, -9.4)) < 0.1)
+# Here are the logKs of aqueous species dissociation reactions at 600 degC and 50000 bar,
+# taken from the Supporting Information of the paper (p. 103-109):
+inorganic.logK <- c(24.4765, -9.0784, -5.3468, 0)
+organic.logK <- c(1.7878, 2.5648, 15.3182, 16.9743, 30.4088, 28.9185)
+# calculate equilibrium constants of the reactions in CHNOSZ; use a negative sign to change from formation to dissociation
+logK.calc <- -unlist(affinity(T=600, P=50000, property="logK")$values)
+logK.calc - c(inorganic.logK, organic.logK)
+# check that we're within 0.021 of the logK values used by SSH14
+stopifnot(maxdiff(logK.calc, c(inorganic.logK, organic.logK)) < 0.021)
+
+#############
+### all done!
# reset the database and previous water computational option
data(OBIGT)
water(oldwat)
Modified: pkg/CHNOSZ/demo/berman.R
===================================================================
--- pkg/CHNOSZ/demo/berman.R 2017-10-07 06:08:37 UTC (rev 245)
+++ pkg/CHNOSZ/demo/berman.R 2017-10-08 14:27:06 UTC (rev 246)
@@ -1,6 +1,8 @@
# CHNOSZ/demo/berman.R 20171003
# make some mineral activity diagrams using Berman (1988) and related data
+res <- 200
+
# using the Helgeson data
# set up basis species
basis(c("K+", "Al+3", "quartz", "H2O", "O2", "H+"))
@@ -9,7 +11,7 @@
# load the species
species(c("K-feldspar", "muscovite", "kaolinite", "pyrophyllite", "andalusite"), "cr")
# calculate affinities in aK+ - temperature space
-a <- affinity(`K+`=c(0, 6), T=c(25, 650), P=1000)
+a <- affinity(`K+`=c(0, 5, res), T=c(200, 650, res), P=1000)
# note that we go just past the quartz transition, but it has no effect on the diagram
diagram(a)
@@ -22,7 +24,8 @@
# load the Berman minerals
species(c("K-feldspar", "muscovite", "kaolinite", "pyrophyllite", "andalusite"), "cr_Berman")
# calculate affinities in aK+ - temperature space
-a <- affinity(`K+`=c(0, 6), T=c(25, 650), P=1000)
+a <- affinity(`K+`=c(0, 5, res), T=c(200, 650, res), P=1000)
diagram(a, add=TRUE, names="", col="blue", lwd=2)
-legend("topleft", lty=c(1, 1), lwd=c(1, 2), col=c("black", "blue"), legend=c("Helgeson et al., 1978", "Berman, 1988"))
+legend("topleft", lty=c(1, 1), lwd=c(1, 2), col=c("black", "blue"),
+ legend=c("Helgeson et al., 1978 (unadjusted)", "Berman, 1988 (adjusted by Sverjensky et al., 1991)"))
Modified: pkg/CHNOSZ/inst/NEWS
===================================================================
--- pkg/CHNOSZ/inst/NEWS 2017-10-07 06:08:37 UTC (rev 245)
+++ pkg/CHNOSZ/inst/NEWS 2017-10-08 14:27:06 UTC (rev 246)
@@ -1,26 +1,40 @@
-CHANGES IN CHNOSZ 1.1.0-44 (2017-10-07)
+CHANGES IN CHNOSZ 1.1.0-45 (2017-10-08)
---------------------------------------
MAJOR CHANGES:
+- Add berman() function and extdata/Berman/*.csv files for calculating
+ thermodynamic properties of minerals using equations of Berman, 1988.
+
+- Calculations related to Berman's (1988) Figs. 1 and 2 for the lambda
+ transition of quartz are presented in the new demo lambda.R.
+
- Add functions implementing the Deep Earth Water (DEW) model
(Sverjensky et al., 2014): water.DEW() and its supporting functions
calculateDensity(), calculateGibbsofWater(), calculateEpsilon(),
calculateQ().
-- The computational setting for water (thermo$opt$par) can now be
- set using water("DEW"), water("IAPWS"), etc.
+- The computational setting for water (thermo$opt$par) can now be set
+ using water("DEW"), water("IAPWS"), etc.
-- Usage of the DEW model is demonstrated in the new demo DEW.R.
+- Usage of the DEW model is shown in the new demo DEW.R. This demo also
+ depends on the Berman equations (above) and, for the last diagram, the
+ following two NEWS items:
-- Add berman() function and extdata/Berman/*.csv files for calculating
- thermodynamic properties of mineral using equations of Berman, 1988.
+- In equilibrate(), it is now possible to combine affinity calculations
+ with changing activity of the balancing basis species (loga.balance).
+ For example, in the last plot of the DEW demo, the affinity transect
+ involves simultaneously varying temperature and logfO2 (given as
+ arguments to affinity()) as well as total concentration of carbon
+ (approximated by the loga.balance argument in equilibrate()).
-- Calculations related to Berman's (1988) Figs. 1 and 2 for the lambda
- transition of quartz are presented in the new demo lambda.R.
-
COMPUTATIONAL IMPROVEMENTS:
+- The 'alpha' argument of diagram() can be set to 'balance' to scale
+ the values by the balancing component. This is useful for making
+ "percent carbon" plots for systems where the species have different
+ carbon numbers.
+
- For minerals with phase transitions (states 'cr2' 'cr3' etc.) in
thermo$obigt (i.e. the Helgeson minerals), it is now possible to use
the minerals in basis(), species(), affinity() with proper accounting
Modified: pkg/CHNOSZ/man/diagram.Rd
===================================================================
--- pkg/CHNOSZ/man/diagram.Rd 2017-10-07 06:08:37 UTC (rev 245)
+++ pkg/CHNOSZ/man/diagram.Rd 2017-10-08 14:27:06 UTC (rev 246)
@@ -27,7 +27,7 @@
\arguments{
\item{eout}{list, object returned by \code{\link{equilibrate}} or \code{\link{affinity}}}
\item{what}{character, what property to calculate and plot}
- \item{alpha}{logical, for speciation diagrams, plot degree of formation instead of activities?}
+ \item{alpha}{logical or character (\samp{balance}), for speciation diagrams, plot degree of formation instead of activities?}
\item{normalize}{logical, divide chemical affinities by balance coefficients (rescale to whole formulas)?}
\item{as.residue}{logical, divide chemical affinities by balance coefficients (no rescaling)?}
\item{balance}{character, balancing constraint; see \code{\link{equilibrate}}}
@@ -80,6 +80,7 @@
The allowed variables are any that \code{\link{affinity}} accepts: temperature, pressure, or the chemical activities of the basis species
If \code{alpha} is TRUE, the fractional degrees of formation (ratios of activities to total activity) are plotted.
+Or, setting \code{alpha} to \samp{balance} allows the activities to be multiplied by the number of the balancing component; this is useful for making \dQuote{percent carbon} diagrams where the species differ in carbon number.
If \code{groups} is supplied, the activities of the species identified in each numeric element of this list are multiplied by the balance coefficients of the species, then summed together.
The names of the list are used to label the lines or fields for the summed activities of the resulting groups.
Modified: pkg/CHNOSZ/man/equilibrate.Rd
===================================================================
--- pkg/CHNOSZ/man/equilibrate.Rd 2017-10-07 06:08:37 UTC (rev 245)
+++ pkg/CHNOSZ/man/equilibrate.Rd 2017-10-08 14:27:06 UTC (rev 246)
@@ -23,8 +23,8 @@
\item{normalize}{logical, normalize the molar formulas of species by the balancing coefficients?}
\item{as.residue}{logical, report results for the normalized formulas?}
\item{Astar}{numeric, affinities of formation reactions excluding species contribution}
- \item{n.balance}{numeric, number of moles of conserved component in the formation reactions of the species of interest}
- \item{loga.balance}{numeric, logarithm of total activity of balanced quantity}
+ \item{n.balance}{numeric, number of moles of balancing component in the formation reactions of the species of interest}
+ \item{loga.balance}{numeric (single value or vector), logarithm of total activity of balanced quantity}
\item{method}{character, equilibration method to use}
\item{tol}{numeric, convergence tolerance for \code{\link{uniroot}}}
}
@@ -37,7 +37,7 @@
\code{aout} contains the chemical affinities of formation reactions of each species of interest.
\code{equilibrate} needs to be provided constraints on how to balance the reactions representing transformations between the species.
-\code{balance} indicates the balancing constraints, according to the following scheme:
+\code{balance} indicates the balancing component, according to the following scheme:
\itemize{
\item \samp{NULL}: autoselect
@@ -53,6 +53,7 @@
NOTE: The summation of activities assumes an ideal system, so \sQuote{molality} is equivalent to \sQuote{activity} here.
\code{loga.balance} gives the logarithm of the total activity of \code{balance} (which is total activity of species for \samp{1} or total activity of amino acid residue-equivalents for \samp{length}).
If \code{loga.balance} is missing, its value is taken from the activities of species listed in \code{aout}; this default is usually the desired operation.
+The supplied value of \code{loga.balance} may also be a vector of values, with length corresponding to the number of conditions in the calculations of \code{\link{affinity}}.
\code{normalize} if TRUE indicates to normalize the molar formulas of species by the balance coefficients.
This operation is intended for systems of proteins, whose conventional formulas are much larger than the basis speices.
@@ -70,15 +71,13 @@
The input values to \code{equil.reaction} and \code{equil.boltzmann} are in a list, \code{Astar}, all elements of the list having the same dimensions; they can be vectors, matrices, or higher-dimensionsal arrays.
Each list element contains the chemical affinities of the formation reactions of one of the species of interest (in dimensionless base-10 units, i.e. A/2.303RT), calculated at unit activity of the species of interest.
-The equilibrium activities (in base-10 log units) of the species of interest returned by either function satisfy the constraints that 1) the final chemical affinities of the formation reactions of the species are all equal and 2) the total activity of the conserved component is equal to (\code{loga.balance}).
+The equilibrium base-10 logarithm activities of the species of interest returned by either function satisfy the constraints that 1) the final chemical affinities of the formation reactions of the species are all equal and 2) the total activity of the balancing component is equal to (\code{loga.balance}).
The first constraint does \emph{not} impose a complete equilibrium, where the affinities of the formation reactions are all equal to zero, but allows for a metastable equilibrium, where the affinities of the formation reactions are equal to each other.
In \code{equil.reaction} (the algorithm described in Dick, 2008 and the only one available prior to CHNOSZ_0.8), the calculations of relative abundances of species are based on a solving a system of equations representing the two constraints stated above.
-A close approximation of the solution is found using \code{\link{uniroot}}.
-Prior to CHNOSZ_0.9-9, the values in the \code{Astar} were used to generate initial guesses of the logarithms of activities of species; values of \code{loga.balance} that were hugely different from these guesses could generate errors from \code{uniroot} such as "f() values at end points not of opposite sign".
-Now (from version 0.9-9), a more flexible method for generating guesses is in place.
+The solution is found using \code{\link{uniroot}} with a flexible method for generating initial guesses.
-In \code{equil.boltzmann} (algorithm available beginning with CHNOSZ_0.8), the chemical activities of species are calculated using the Boltzmann distribution.
+In \code{equil.boltzmann}, the chemical activities of species are calculated using the Boltzmann distribution.
This calculation is faster than the algorithm of \code{equil.reaction}, but is limited to systems where the transformations are all balanced on one mole of species.
If \code{equil.boltzmann} is called with \code{balance} other than \samp{1}, it stops with an error.
@@ -106,10 +105,10 @@
res <- 100
Astar <- affinity(pH=c(0, 14, res))$values
# the logarithms of activity for a total activity
-# of the balanced quantity (C or CO2) equal to 0.001
+# of the balancing component (CO2) equal to 0.001
loga.boltz <- equil.boltzmann(Astar, c(1, 1, 1), 0.001)
# calculated another way
-loga.react <- equil.reaction(Astar, c(1, 1, 1), 0.001)
+loga.react <- equil.reaction(Astar, c(1, 1, 1), rep(0.001, 100))
# probably close enough for most purposes
stopifnot(all.equal(loga.boltz, loga.react))
# the first ionization constant (pKa)
Modified: pkg/CHNOSZ/tests/testthat/test-equilibrate.R
===================================================================
--- pkg/CHNOSZ/tests/testthat/test-equilibrate.R 2017-10-07 06:08:37 UTC (rev 245)
+++ pkg/CHNOSZ/tests/testthat/test-equilibrate.R 2017-10-08 14:27:06 UTC (rev 246)
@@ -122,6 +122,38 @@
expect_error(equilibrate(a, balance=1.000001, method="boltzmann"), "won't run equil.boltzmann")
})
+test_that("equilibrate() can be used with a vector of loga.balance values", {
+ # system is balanced on CO2; species have different number of C, so loga.balance affects the equilibrium activities
+ basis("CHNOS")
+ species(c("formic acid", "acetic acid", "propanoic acid"))
+ # calculate reference values at logfO2 = -80, log(aCO2tot) = -6
+ basis("O2", -80)
+ a <- affinity()
+ e80.6 <- unlist(equilibrate(a, loga.balance = -6)$loga.equil)
+ # calculate reference values at logfO2 = -60, log(aCO2tot) = -8
+ basis("O2", -60)
+ a <- affinity()
+ e60.8 <- unlist(equilibrate(a, loga.balance = -8)$loga.equil)
+ # calculate affinity on a transect: logfO2 from -80 to -60
+ O2 <- seq(-80, -60)
+ aO2 <- affinity(O2=O2)
+ # values calculated at the ends of the transect should be the same as above
+ eO2.6 <- equilibrate(aO2, loga.balance = -6)$loga.equil
+ expect_equal(list2array(eO2.6)[1, ], e80.6)
+ eO2.8 <- equilibrate(aO2, loga.balance = -8)$loga.equil
+ expect_equal(list2array(eO2.8)[21, ], e60.8)
+ # make an vector of loga.balance to go with the logfO2 transect
+ # first we make a vector with non-matching length, and get an error
+ logaCO2.wronglen <- seq(-6, -8)
+ expect_error(equilibrate(aO2, loga.balance=logaCO2.wronglen), "length of loga.balance \\(3) doesn't match the affinity values \\(21)")
+ # now do it with the correct length
+ logaCO2 <- seq(-6, -8, length.out=length(O2))
+ eO2 <- equilibrate(aO2, loga.balance=logaCO2)
+ # now the first set of conditions is logfO2 = -80, log(aCO2tot) = -6
+ expect_equal(list2array(eO2$loga.equil)[1, ], e80.6)
+ # and the final set is logfO2 = -80, log(aCO2tot) = -6
+ expect_equal(list2array(eO2$loga.equil)[21, ], e60.8)
+})
# references
Modified: pkg/CHNOSZ/vignettes/anintro.Rmd
===================================================================
--- pkg/CHNOSZ/vignettes/anintro.Rmd 2017-10-07 06:08:37 UTC (rev 245)
+++ pkg/CHNOSZ/vignettes/anintro.Rmd 2017-10-08 14:27:06 UTC (rev 246)
@@ -930,7 +930,7 @@
That makes a diagram that is similar to Figure 6b of Shock and Schulte (1998).
```{marginfigure}
-Some differences from the original diagrams may be caused by the sensitivity of the calculations to log*f*<sub>O<sub>2</sub></sub>, for which the values we use here may carry artifacts introduced by digitization of their plot.
+Some differences from the original diagrams could be caused by the sensitivity of the calculations to log*f*<sub>O<sub>2</sub></sub>, for which the values we use here may carry artifacts introduced by digitization of their plot.
```
It is also possible to plot the distribution of species within individual groups, such as alcohols and ketones.
We do this by setting the names and line types for the *other* species to values that prevent them from being plotted:
More information about the CHNOSZ-commits
mailing list