[CHNOSZ-commits] r267 - in pkg/CHNOSZ: . demo inst man tests/testthat vignettes
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Oct 25 17:50:57 CEST 2017
Author: jedick
Date: 2017-10-25 17:50:57 +0200 (Wed, 25 Oct 2017)
New Revision: 267
Added:
pkg/CHNOSZ/tests/testthat/test-logmolality.R
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/demo/DEW.R
pkg/CHNOSZ/inst/NEWS
pkg/CHNOSZ/man/CHNOSZ-package.Rd
pkg/CHNOSZ/vignettes/anintro.Rmd
pkg/CHNOSZ/vignettes/eos-regress.Rmd
pkg/CHNOSZ/vignettes/vig.bib
Log:
document activity-molality transformation in anintro.Rmd and test-logmolality.R
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2017-10-24 13:29:31 UTC (rev 266)
+++ pkg/CHNOSZ/DESCRIPTION 2017-10-25 15:50:57 UTC (rev 267)
@@ -1,6 +1,6 @@
-Date: 2017-10-24
+Date: 2017-10-25
Package: CHNOSZ
-Version: 1.1.0-65
+Version: 1.1.0-66
Title: Thermodynamic Calculations for Geobiochemistry
Author: Jeffrey Dick
Maintainer: Jeffrey Dick <j3ffdick at gmail.com>
@@ -11,12 +11,11 @@
compositional biology. Thermodynamic properties are taken from a database for minerals
and inorganic and organic aqueous species including biomolecules, or from amino acid
group additivity for proteins. High-temperature properties are calculated using the
- revised Helgeson-Kirkham-Flowers equations of state for aqueous species. Functions are
+ revised Helgeson-Kirkham-Flowers equations of state for aqueous species, and activity
+ coefficients can be calculated for specified ionic strength. Functions are
provided to define a system using basis species, automatically balance reactions,
- calculate the chemical affinities of reactions for selected species, and plot the results
- on potential diagrams or equilibrium activity diagrams. Experimental features are
- available to calculate activity coefficients for aqueous species or for multidimensional
- optimization of thermodynamic variables using an objective function.
+ calculate the chemical affinities of formation reactions for selected species, calculate
+ equilibrium activities, and plot the results on chemical activity diagrams.
License: GPL (>= 2)
BuildResaveData: no
VignetteBuilder: knitr
Modified: pkg/CHNOSZ/demo/DEW.R
===================================================================
--- pkg/CHNOSZ/demo/DEW.R 2017-10-24 13:29:31 UTC (rev 266)
+++ pkg/CHNOSZ/demo/DEW.R 2017-10-25 15:50:57 UTC (rev 267)
@@ -134,9 +134,8 @@
# pH set by jadeite + kyanite + coesite
# output from EQ3NR calculations (SSH14 Supporting Information)
# dissolved carbon: 0.03, 0.2, 1, 4, 20 molal
-# ionic strength: 0.39, 0.57, 0.88, 1.45, 2.49
+# true ionic strength: 0.39, 0.57, 0.88, 1.45, 2.49
# pH: 3.80, 3.99, 4.14, 4.25, 4.33
-T <- seq(600, 1000, 5)
## activate DEW model
data(thermo)
water("DEW")
@@ -151,33 +150,33 @@
mod.buffer("QFM_Berman", c("quartz", "fayalite", "magnetite"), "cr_Berman", 0)
## calculate logfO2 in QFM buffer
basis("O2", "QFM_Berman")
+T <- seq(600, 1000, 5); T100 <- seq(600, 1000, 100)
buf <- affinity(T=T, P=50000, return.buffer=TRUE)
## add species
species(c(inorganics, organics))
## generate spline functions from IS, pH, and molC values at every 100 degC
-IS <- splinefun(T[!T%%100], c(0.39, 0.57, 0.88, 1.45, 2.49))
-pH <- splinefun(T[!T%%100], c(3.80, 3.99, 4.14, 4.25, 4.33))
-molC <- splinefun(T[!T%%100], c(0.03, 0.2, 1, 4, 20))
+IS <- splinefun(T100, c(0.39, 0.57, 0.88, 1.45, 2.49))
+pH <- splinefun(T100, c(3.80, 3.99, 4.14, 4.25, 4.33))
+molC <- splinefun(T100, c(0.03, 0.2, 1, 4, 20))
## use Debye-Huckel equation with B-dot set to zero
nonideal("Helgeson0")
## calculate affinities on the T-logfO2-pH-IS transect
-a <- affinity(T=T, O2=buf$O2 - 2, IS=IS(T), pH=pH(T), P=50000)
+a <- affinity(T = T, O2 = buf$O2 - 2, IS = IS(T), pH = pH(T), P = 50000)
## calculate metastable equilibrium activities using the total
## carbon molality as an approximation of total activity
-e <- equilibrate(a, loga.balance=log10(molC(T)))
+e <- equilibrate(a, loga.balance = log10(molC(T)))
## 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", ylab="carbon fraction", names=names, col=col, ylim=c(0, 0.8))
+diagram(e, alpha = "balance", ylab = "carbon fraction", names = names, col = col, ylim = c(0, 0.8), ylab="carbon fraction")
## add legend and title
ltxt1 <- "P = 50000 bar"
-ltxt2 <- substitute(logfO2=="QFM-2", list(logfO2=axis.label("O2")))
-ltxt3 <- "pH = 4"
-pH <- seq(3.8, 4.3, length.out=length(T))
-legend("left", legend=as.expression(c(ltxt1, ltxt2, ltxt3)), bty="n")
+ltxt2 <- substitute(logfO2=="QFM-2", list(logfO2 = axis.label("O2")))
+pH <- seq(3.8, 4.3, length.out = length(T))
+legend("left", legend = as.expression(c(ltxt1, ltxt2)), bty = "n")
t1 <- "Aqueous carbon speciation"
t2 <- "after Sverjensky et al., 2014b"
mtitle(c(t1, t2))
@@ -192,7 +191,7 @@
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 <- -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)
@@ -201,7 +200,7 @@
# activity coefficients for monovalent species from EQ3NR output
loggamma <- c(-0.15, -0.18, -0.22, -0.26, -0.31)
# activity coefficients calculated in CHNOSZ
-sres <- subcrt("propanoate", T=seq(600, 1000, 100), P=50000, IS=c(0.39, 0.57, 0.88, 1.45, 2.49))
+sres <- subcrt("propanoate", T = seq(600, 1000, 100), P = 50000, IS = c(0.39, 0.57, 0.88, 1.45, 2.49))
stopifnot(maxdiff(sres$out[[1]]$loggam, loggamma) < 0.004)
###########
Modified: pkg/CHNOSZ/inst/NEWS
===================================================================
--- pkg/CHNOSZ/inst/NEWS 2017-10-24 13:29:31 UTC (rev 266)
+++ pkg/CHNOSZ/inst/NEWS 2017-10-25 15:50:57 UTC (rev 267)
@@ -1,4 +1,4 @@
-CHANGES IN CHNOSZ 1.1.0-64 (2017-10-23)
+CHANGES IN CHNOSZ 1.1.0-66 (2017-10-25)
---------------------------------------
MAJOR CHANGES:
@@ -143,6 +143,11 @@
- Add demo TCA.R for standard Gibbs energies of steps of the
citric acid cycle (Canovas and Shock, 2016).
+- Add test-logmolality.R to demonstrate transformation of variables from
+ activity to molality in the main workflow; refer to this test and
+ describe these transformations in anintro.Rmd; remove "experimental"
+ labeling of activity coefficient calculations in CHNOSZ-package.Rd.
+
CLEANUP:
- To save space, taxid_names.csv has been trimmed to hold only those
Modified: pkg/CHNOSZ/man/CHNOSZ-package.Rd
===================================================================
--- pkg/CHNOSZ/man/CHNOSZ-package.Rd 2017-10-24 13:29:31 UTC (rev 266)
+++ pkg/CHNOSZ/man/CHNOSZ-package.Rd 2017-10-25 15:50:57 UTC (rev 267)
@@ -60,7 +60,7 @@
\item Buffer calculations - compute activities of basis species that are determined by a buffer of one or more species (e.g., pyrite-pyrrhotite-magnetite; acetic acid-\CO2) (\code{\link{buffer}}).
- \item Activity coefficients (\bold{experimental}) - calculate activity coefficients of aqueous species using the extended Debye-Hückel equation (\code{\link{nonideal}}).
+ \item Activity coefficients - calculate activity coefficients of aqueous species using the extended Debye-Hückel equation (\code{\link{nonideal}}).
\item Activity statistics (\bold{experimental}) - calculate summary statistics for equilibrium activities of species (\code{\link{revisit}}).
Added: pkg/CHNOSZ/tests/testthat/test-logmolality.R
===================================================================
--- pkg/CHNOSZ/tests/testthat/test-logmolality.R (rev 0)
+++ pkg/CHNOSZ/tests/testthat/test-logmolality.R 2017-10-25 15:50:57 UTC (rev 267)
@@ -0,0 +1,124 @@
+context("logmolality")
+
+test_that("non-zero ionic strength transforms variables from activity to molality", {
+ # what happens with activity coefficients when using subcrt() to calculate affinity,
+ # and in the rest of the main workflow of CHNOSZ?
+ # 20171025
+
+ # first get the activity coefficients of H+ and HCO3-
+ # the long way...
+ wprop <- water(c("A_DH", "B_DH"), P=1)
+ nonid <- nonideal(c("H+", "HCO3-"), subcrt(c("H+", "HCO3-"), T=25)$out, IS=1, T=298.15, P=1, A_DH=wprop$A_DH, B_DH=wprop$B_DH)
+ # here we have a small difference due to rounding of the expected value:
+ expect_maxdiff(nonid[[2]]$loggam, -0.1890084, 1e-7)
+ # the short way...
+ loggam <- subcrt(c("H+", "HCO3-"), T=25, IS=1)$out[[2]]$loggam
+ # we get that loggam(H+)=0 and loggam(HCO3-)=-0.189
+ expect_equal(nonid[[2]]$loggam, loggam)
+
+ ## take-home lesson 0: with default settings,
+ ## the activity coefficient of H+ is always 1
+
+ # how do activity coefficient affect the value of G?
+ # let's step back and look at the *standard Gibbs energy* at IS = 0
+ out0 <- subcrt(c("H+", "HCO3-"), T=25)$out
+ # and at IS = 1
+ out1 <- subcrt(c("H+", "HCO3-"), T=25, IS=1)$out
+ # the apparent standard Gibbs energy is less than the standard Gibbs energy
+ # by an amount determined by the activity coefficient
+ expect_equal(out1[[2]]$G - out0[[2]]$G, -convert(loggam, "G"))
+
+ ## take-home lesson 1.5: setting IS in subcrt() gives apparent standard Gibbs energy
+
+ # now, what is the equilibrium constant for the reaction CO2 + H2O = H+ + HCO3-?
+ # (this is the standard state property at IS=0)
+ logK <- subcrt(c("CO2", "H2O", "H+", "HCO3-"), c(-1, -1, 1, 1), T=25)$out$logK
+ # we get logK = -6.344694 (rounded)
+ expect_maxdiff(logK, -6.344694, 1e-6)
+
+ # now, what is the affinity of the reaction at pH=7 and molalities of HCO3- and CO2 = 10^-3?
+
+ ## case 1: ionic strength = 0, so gamma = 0 and activity = molality
+ # first calculate it by hand from 2.303RTlog(K/Q)
+ # logQ = (logaH+ + logaHCO3-) - (logaH2O + logaCO2)
+ logQ0 <- (-7 + -3) - (0 + -3) # i.e. 7
+ # convert logs to energy according to G = -2.303RTlogK
+ # (we negate that to get affinity)
+ A0manual <- -convert(logK - logQ0, "G")
+ # now the subcrtmatic calculation with subcrt
+ A0subcrt <- subcrt(c("CO2", "H2O", "H+", "HCO3-"), c(-1, -1, 1, 1), T=25, logact=c(-3, 0, -7, -3))$out$A
+ # we get the same affinity!
+ expect_equal(A0subcrt, A0manual)
+
+ ## case 1: ionic strength = 0, so activity = molality * gamma
+ logaHCO3 = -3 + loggam
+ logQ1 <- (-7 + logaHCO3) - (0 + -3)
+ A1manual <- -convert(logK - logQ1, "G")
+ A1subcrt <- subcrt(c("CO2", "H2O", "H+", "HCO3-"), c(-1, -1, 1, 1), T=25, logact=c(-3, 0, -7, -3), IS=1)$out$A
+ expect_equal(A1subcrt, A1manual)
+
+ ## take-home lesson 1: using subcrt with IS not equal to zero, the "logact"
+ ## argument is logmolal in affinity calculations for charged aqueous species
+
+ # now, calculate the affinities using affinity()
+ basis("CHNOS+") # pH=7, logaCO2 = -3
+ species(c("CO2", "HCO3-")) # logactivities = -3
+
+ ## case 1: IS = 0
+ a0 <- affinity()
+ # that gives us values in log units; convert to energy
+ # (HCO3- is species #2)
+ A0affinity <- -convert(a0$values[[2]], "G")
+ expect_equal(A0affinity[[1]], A0subcrt)
+
+ ## case 2: IS = 1
+ a1 <- affinity(IS=1)
+ A1affinity <- -convert(a1$values[[2]], "G")
+ expect_equal(A1affinity[[1]], A1subcrt)
+
+ ## take-home lesson 2: using affinity() with IS not equal to zero, the "logact"
+ ## set by species() is logmolal in affinity calculations for charged aqueous species
+
+ # now, swap HCO3- for CO2 in the basis
+ swap.basis("CO2", "HCO3-")
+ basis("HCO3-", -3)
+ a0 <- affinity()
+ a1 <- affinity(IS=1)
+ # look at HCO3 formation affinity:
+ # they're both zero at IS = 0 or 1
+ expect_equal(a0$values[[2]][1], 0)
+ expect_equal(a1$values[[2]][1], 0)
+ # look at CO2 formation affinity:
+ ACO2_0affinity <- -convert(a0$values[[1]], "G")
+ ACO2_1affinity <- -convert(a1$values[[1]], "G")
+ # what's that?? we're looking at the reverse of the above
+ # reaction, i.e. H+ + HCO3- = CO2 + H2O
+ # so, logK = 6.345
+ logKrev <- -logK
+ logQrev0 <- -logQ0
+ logQrev1 <- -logQ1
+ ACO2_0manual <- -convert(logKrev - logQrev0, "G")
+ ACO2_1manual <- -convert(logKrev - logQrev1, "G")
+ expect_equal(ACO2_0manual, ACO2_0affinity[[1]])
+ expect_equal(ACO2_1manual, ACO2_1affinity[[1]])
+
+ ## take-home lesson 3: using affinity() with IS not equal to zero, the "logact"
+ ## set by basis() is logmolal in affinity calculations for charged aqueous species
+
+ # now look at equilibrate()
+ e0 <- equilibrate(a0)
+ e1 <- equilibrate(a1)
+ # using the equilibrated values, calculate affinity of the reaction CO2 + H2O = H+ + HCO3-
+ # case 1: IS = 0
+ logQeq0 <- (-7 + e0$loga.equil[[2]]) - (e0$loga.equil[[1]] + 0)
+ Aeq0 <- -convert(logK - logQeq0, "G") # zero!
+ expect_equal(Aeq0[[1]], 0)
+ # case 2: IS = 1
+ # here, loga.equil is the *molality*, so we must multiply by loggam
+ logQeq1 <- (-7 + e1$loga.equil[[2]] + loggam) - (e1$loga.equil[[1]] + 0)
+ Aeq1 <- -convert(logK - logQeq1, "G") # zero!
+ expect_equal(Aeq1[[1]], 0)
+
+ ## take-home lesson 4: using affinity() with IS not equal to zero, the "loga.equil"
+ ## returned by equilibrate() is logmolal for speciation calculations with charged aqueous species
+})
Modified: pkg/CHNOSZ/vignettes/anintro.Rmd
===================================================================
--- pkg/CHNOSZ/vignettes/anintro.Rmd 2017-10-24 13:29:31 UTC (rev 266)
+++ pkg/CHNOSZ/vignettes/anintro.Rmd 2017-10-25 15:50:57 UTC (rev 267)
@@ -796,8 +796,8 @@
CHNOSZ does not take account of all possible reactions in the speciation of a system.
Instead, it assumes that the total activity of species in the system is set by the activity of *one* basis species.
```{marginfigure}
-If activity coefficients are assumed to be zero, activities are equal to concentration and we can refer to "total activity".
-<span style="color:green">`nonideal()`</span> ([see below](#activity-coefficients)) is an experimental feature to deal with activity coefficients.
+When activity coefficients are assumed to be zero, activities are equal to concentration and we can refer to "total activity".
+If the ionic strength is specified, <span style="color:green">`nonideal()`</span> ([see below](#activity-coefficients)) can be used to calculate activity coefficients.
```
This balanced, or conserved, basis species must be present (with a positive or negative coefficient) in the formation reactions of all species considered.
@@ -981,6 +981,179 @@
For instance, the widespread rule of thumb for balancing mineral reactions on a chemical component is unrealistic for processes where volume is conserved [@MD98].
While choosing an inappropriate balance leads to infeasible models, consideration of the different possibilities might give insight into the conditions affecting the dynamics of some systems.
+# Activity coefficients
+
+<span style="color:green">`nonideal()`</span> uses the extended Debye--Hückel equation as formulated by @HKF81 for NaCl-dominated solutions, or optionally using parameters described in Chapter 3 of @Alb03, to calculate activity coefficients of charged species.
+The activity coefficients are calculated as a function of ionic strength (*I*), temperature, and charge of each species, without any other species-specific parameters.
+Using the default `Helgeson` method, the extended term parameter ("B-dot") is derived from data of @Hel69 and Helgeson et al. (1981), and extrapolations of @MSS13.
+
+## Transformation of variables
+
+Following the main workflow of CHNOSZ, <span style="color:green">`nonideal()`</span> normally does not need to be used directly.
+Intead, invoke the calculations by setting the `IS` argument in <span style="color:green">`subcrt()`</span> or <span style="color:green">`affinity()`</span>.
+There are a few things to remember when using activity coefficients:
+
+* H<sup>+</sup> is assumed to behave ideally, so its activity coefficient is 1 for any ionic strength. You can calculate activity coefficients of H<sup>+</sup> by setting `thermo$opt$ideal.H <<- FALSE`.
+
+* Using <span style="color:green">`subcrt()`</span> with `IS` not equal to zero, calculated values of `G` are the **apparent** standard Gibbs energy at specified ionic strength.
+
+* Using <span style="color:green">`subcrt()`</span> with `IS` not equal to zero, values in the `logact` argument stand for **log molality** of charged aqueous species in affinity calculations.
+
+* Using <span style="color:green">`affinity()`</span> with `IS` not equal to zero, the following values stand for **log molality** of charged aqueous species:
+ + values of `logact` set by <span style="color:green">`basis()`</span>;
+ + values of `logact` set by <span style="color:green">`species()`</span>;
+ + value of `loga.equil` returned by <span style="color:green">`equilibrate()`</span>
+
+In other words, the activation of activity coefficients transforms variables from activity to molality in the main workflow.
+A simple but detailed example demonstrating each of the above tranformations is in `tests/testthat/test-logmolality.R`.
+
+Because it is not possible to dynamically change the names of arguments in the functions, the user should be aware of the effects mentioned above.
+As another consequence, activities of basis species or equilibrium activities of species shown on diagrams should be relabeled as molalities when `IS` is used in the calculation of <span style="color:green">`affinity()`</span>.
+
+## Biochemical example
+
+For the following calculations, we change the nonideality method to `Alberty`; this is a simpler formulation with parameters that are suitable for biochemical species at relatively low temperatures:
+```{r Alberty}
+oldnon <- nonideal("Alberty")
+```
+
+Let's take a look at calculated activity coefficients at two temperatures and their effect on the standard Gibbs energies of formation (Δ*G*°<sub>*f*</sub>) of species with different charge:
+```{r subcrt_IS}
+subcrt(c("MgATP-2", "MgHATP-", "MgH2ATP"),
+ T = c(25, 100), IS = c(0, 0.25), property = "G")$out
+```
+
+The logarithms of the activity coefficients (`loggam`) are more negative for the higher-charged species, as well as at higher temperature, and have a stabilizing effect.
+That is, the apparent Gibbs energies at *I* > 0 are less than the standard Gibbs energies at *I* = 0.
+
+We can use these calculations to make some speciation plots, similar to Figures 1.2--1.5 in Alberty (2003).
+These figures show the distribution of differently charged species of adenosine triphosphate (ATP) as a function of pH, and the average number of H<sup>+</sup> and Mg<sup>+2</sup> bound to ATP in solution as a function of pH or pMg (-log*a*<sub>Mg<sup>+2</sup></sub>).
+
+Use <span style="color:green">`info()`</span> to see what ATP species are available.
+The sources of high-temperature thermodynamic data for these species are two papers by LaRowe and Helgeson [- at LH06a; - at LH06b].
+```{r info_ATP, results="hide"}
+info(" ATP")
+```
+
+The plots for this system in Alberty's book were made for *I* = 0.25 M and *T* = 25 °C.
+As a demonstration of CHNOSZ's capabilities, we can assign a temperature of 100 °C.
+```{r T_100}
+T <- 100
+```
+
+Use the following commands to set the basis species, add the variously protonated ATP species, calculate the affinities of the formation reactions, equilibrate the system, and make a degree of formation (α) or mole fraction diagram.
+This is similar to Figure 1.3 of Alberty (2003), but is calculated for *I* = 0 M and *T* = 100 °C:
+```{marginfigure}
+To make the code more readable, commands for plotting titles and legends are not shown.
+All of the commands are available in the source of this document.
+```
+
+```{r ATP, eval=FALSE, echo=2:6}
+par(mfrow = c(1, 4), mar = c(3.1, 3.6, 2.1, 1.6), mgp = c(1.8, 0.5, 0))
+basis("MgCHNOPS+")
+species(c("ATP-4", "HATP-3", "H2ATP-2", "H3ATP-", "H4ATP"))
+a <- affinity(pH = c(3, 9), T = T)
+e <- equilibrate(a)
+d <- diagram(e, alpha = TRUE, tplot = FALSE)
+title(main = describe.property("T", T))
+alphas <- do.call(rbind, d$plotvals)
+nH <- alphas * 0:4
+Hlab <- substitute(italic(N)[H^`+`])
+plot(a$vals[[1]], colSums(nH), type = "l", xlab = "pH", ylab=Hlab, lty=2, col=2)
+a <- affinity(pH = c(3, 9), IS = 0.25, T = T)
+e <- equilibrate(a)
+d <- diagram(e, alpha = TRUE, plot.it = FALSE)
+alphas <- do.call(rbind, d$plotvals)
+nH <- alphas * 0:4
+lines(a$vals[[1]], colSums(nH))
+legend("topright", legend = c("I = 0 M", "I = 0.25 M"), lty = 2:1, col = 2:1, cex = 0.8)
+ATP.H <- substitute("ATP and H"^`+`)
+title(main = ATP.H)
+species(c("MgATP-2", "MgHATP-", "MgH2ATP", "Mg2ATP"))
+Hplot <- function(pMg, IS = 0.25) {
+ basis("Mg+2", -pMg)
+ a <- affinity(pH = c(3, 9), IS = IS, T = T)
+ e <- equilibrate(a)
+ d <- diagram(e, alpha = TRUE, plot.it = FALSE)
+ alphas <- do.call(rbind, d$plotvals)
+ NH <- alphas * c(0:4, 0, 1, 2, 0)
+ lines(a$vals[[1]], colSums(NH), lty = 7 - pMg, col = 7 - pMg)
+}
+plot(c(3, 9), c(0, 2), type = "n", xlab = "pH", ylab = Hlab)
+lapply(2:6, Hplot)
+legend("topright", legend = paste("pMg = ", 2:6), lty = 5:1, col = 5:1, cex = 0.8)
+ATP.H.Mg <- substitute("ATP and H"^`+`~"and Mg"^`+2`)
+title(main = ATP.H.Mg)
+Mgplot <- function(pH, IS = 0.25) {
+ basis("pH", pH)
+ a <- affinity(`Mg+2` = c(-2, -7), IS = IS, T = T)
+ e <- equilibrate(a)
+ d <- diagram(e, alpha = TRUE, plot.it = FALSE)
+ alphas <- do.call(rbind, d$plotvals)
+ NMg <- alphas * species()$`Mg+`
+ lines(-a$vals[[1]], colSums(NMg), lty = 10 - pH, col = 10 - pH)
+}
+Mglab <- substitute(italic(N)[Mg^`+2`])
+plot(c(2, 7), c(0, 1.2), type = "n", xlab = "pMg", ylab = Mglab)
+lapply(3:9, Mgplot)
+legend("topright", legend = paste("pH = ", 3:9), lty = 7:1, col = 7:1, cex = 0.8)
+title(main = ATP.H.Mg)
+```
+
+Note that we have saved the numeric results of <span style="color:green">`diagram()`</span>, i.e. the degrees of formation of the species (α).
+With that, we can calculate and plot the average number of protons bound per ATP molecule.
+To do so, we use R's `rbind()` and `do.call()` to turn `alpha` into a matrix, then multiply by the number of protons bound to each species, and sum the columns to get the total (i.e. average proton number, *N*<sub>H<sup>+</sup></sub>):
+```{r ATP, eval=FALSE, echo=8:11}
+```
+
+Adding the `IS` argument to <span style="color:green">`affinity()`</span>, we can now plot *N*<sub>H<sup>+</sup></sub> at the given ionic strength.
+Here we set `plot.it = FALSE` in `diagram()` because we use the computed α to make our own plot.
+This is similar to Figure 1.3 of Alberty (2003), but at higher temperature:
+```{r ATP, eval=FALSE, echo=12:17}
+```
+
+Next, we add the Mg<sup>+2</sup>-complexed ATP species:
+```{r ATP, eval=FALSE, echo=21}
+```
+
+Here is a function to calculate and plot *N*<sub>H<sup>+</sup></sub> for a given pMg:
+```{r ATP, eval=FALSE, echo=22:30}
+```
+
+With that function in hand, we plot the lines corresponding to pMg = 2 to 6.
+This is similar to Figure 1.4 of Alberty (2003):
+```{r ATP, eval=FALSE, echo=31:32}
+```
+
+The next function calculates and plots the average number of Mg<sup>+2</sup> bound to ATP (*N*<sub>Mg<sup>+2</sup></sub>) for a given pH.
+Here we multiply `alpha` by the number of Mg<sup>+2</sup> in each species, and negate log*a*<sub>Mg<sup>+2</sup></sub> (the variable used in <span style="color:green">`affinity()`</span>) to get pMg:
+```{r ATP, eval=FALSE, echo=36:44}
+```
+
+Using that function, we plot the lines corresponding to pH = 3 to 9.
+This is similar to Figure 1.5 of Alberty (2003):
+```{r ATP, eval=FALSE, echo=45:47}
+```
+
+<p>
+```{r ATP, fig.fullwidth=TRUE, fig.width=10, fig.height=2.5, dpi=ifelse(dpi==50, 50, 100), out.width="100%", echo=FALSE, message=FALSE, results="hide", fig.cap="Binding of H<sup>+</sup> and Mg<sup>+2</sup> to ATP at 100 °C and *I* = 0 M (first plot) or *I* = 0.25 M (third and fourth plots).", cache=TRUE, pngquant=pngquant, timeit=timeit}
+```
+</p>
+
+We have calculated the distribution of ATP species and average binding number of H<sup>+</sup> and Mg<sup>+2</sup> for given pH, pMg, ionic strength, and temperature.
+Accounting for the distribution of chemical species lends itself to thermodynamic models for reactions between reactants that have multiple ionized and complexed states.
+In contrast, Alberty (2003) and others propose models for biochemical reactions where the ionized and complexed species are combined into a single representation.
+Those models invoke Legendre-transformed thermodynamic properties, such as transformed Gibbs energies that are tabulated for specified pH, pMg, and ionic strength.
+Although the conceptual pathways are different, the two approaches lead to equivalent results concerning the energetics of the overall reactions and the conditions for equilibrium [@SVI12].
+The example here shows how the required calculations can be performed at the species level using conventional standard Gibbs energies for species referenced to infinite dilution (zero ionic strength).
+The effects of ionic strength are modeled "on the fly" in CHNOSZ by setting the `IS` argument in <span style="color:green">`subcrt()`</span> or <span style="color:green">`affinity()`</span> to invoke the nonideality model on top of the standard Gibbs energies of species.
+
+Now that we're finished, we can reset the nonideality method to the default.
+(This really isn't needed here, because there aren't any nonideality calculations below):
+```{r oldnon}
+nonideal(oldnon)
+```
+
# Proteins
Proteins in CHNOSZ are handled a little bit differently from other species.
@@ -1505,161 +1678,10 @@
```{r bison_transferase, eval=FALSE, echo=12:15}
```
-# Experimental features
+# Optimization of chemical activities
-Some experimental features in CHNOSZ implement provisional algorithms and/or have not been extensively tested.
+**NOTE: This is an experimental feature.**
-## Activity coefficients
-
-<span style="color:green">`nonideal()`</span> uses the extended Debye--Hückel equation, described in Chapter 3 of @Alb03, to calculate activity coefficients of charged species.
-The calculations can be invoked by setting the `IS` argument in <span style="color:green">`subcrt()`</span> or <span style="color:green">`affinity()`</span>.
-The activity coefficients are calculated as a function of ionic strength (*I*), temperature, and charge of each species, without any other species-specific parameters.
-Due to these assumptions, as well as the limited applicability of the implemented equations to very high *I* and/or *T*, <span style="color:green">`nonideal()`</span> is presented as an experimental feature.
-
-For the following examples, we change the nonideality method to `Alberty`:
-```{r Alberty}
-oldnon <- nonideal("Alberty")
-```
-
-Let's take a look at calculated activity coefficients at two temperatures and their effect on the standard Gibbs energies of formation (Δ*G*°<sub>*f*</sub>) of species with different charge:
-```{r subcrt_IS}
-subcrt(c("MgATP-2", "MgHATP-", "MgH2ATP"),
- T = c(25, 100), IS = c(0, 0.25), property = "G")$out
-```
-
-The logarithms of the activity coefficients (`loggam`) are more negative for the higher-charged species, as well as at higher temperature, and have a stabilizing effect.
-That is, the apparent Gibbs energies at *I* > 0 are less than the standard Gibbs energies at *I* = 0.
-
-We can use these calculations to make some speciation plots, similar to Figures 1.2--1.5 in Alberty (2003).
-These figures show the distribution of differently charged species of adenosine triphosphate (ATP) as a function of pH, and the average number of H<sup>+</sup> and Mg<sup>+2</sup> bound to ATP in solution as a function of pH or pMg (-log*a*<sub>Mg<sup>+2</sup></sub>).
-
-Use <span style="color:green">`info()`</span> to see what ATP species are available.
-The sources of high-temperature thermodynamic data for these species are two papers by LaRowe and Helgeson [- at LH06a; - at LH06b].
-```{r info_ATP, results="hide"}
-info(" ATP")
-```
-
-The plots for this system in Alberty's book were made for *I* = 0.25 M and *T* = 25 °C.
-As a demonstration of CHNOSZ's capabilities, we can assign a temperature of 100 °C.
-```{r T_100}
-T <- 100
-```
-
-Use the following commands to set the basis species, add the variously protonated ATP species, calculate the affinities of the formation reactions, equilibrate the system, and make a degree of formation (α) or mole fraction diagram.
-This is similar to Figure 1.3 of Alberty (2003), but is calculated for *I* = 0 M and *T* = 100 °C:
-```{marginfigure}
-To make the code more readable, commands for plotting titles and legends are not shown.
-All of the commands are available in the source of this document.
-```
-
-```{r ATP, eval=FALSE, echo=2:6}
-par(mfrow = c(1, 4), mar = c(3.1, 3.6, 2.1, 1.6), mgp = c(1.8, 0.5, 0))
-basis("MgCHNOPS+")
-species(c("ATP-4", "HATP-3", "H2ATP-2", "H3ATP-", "H4ATP"))
-a <- affinity(pH = c(3, 9), T = T)
-e <- equilibrate(a)
-d <- diagram(e, alpha = TRUE, tplot = FALSE)
-title(main = describe.property("T", T))
-alphas <- do.call(rbind, d$plotvals)
-nH <- alphas * 0:4
-Hlab <- substitute(italic(N)[H^`+`])
-plot(a$vals[[1]], colSums(nH), type = "l", xlab = "pH", ylab=Hlab, lty=2, col=2)
-a <- affinity(pH = c(3, 9), IS = 0.25, T = T)
-e <- equilibrate(a)
-d <- diagram(e, alpha = TRUE, plot.it = FALSE)
-alphas <- do.call(rbind, d$plotvals)
-nH <- alphas * 0:4
-lines(a$vals[[1]], colSums(nH))
-legend("topright", legend = c("I = 0 M", "I = 0.25 M"), lty = 2:1, col = 2:1, cex = 0.8)
-ATP.H <- substitute("ATP and H"^`+`)
-title(main = ATP.H)
-species(c("MgATP-2", "MgHATP-", "MgH2ATP", "Mg2ATP"))
-Hplot <- function(pMg, IS = 0.25) {
- basis("Mg+2", -pMg)
- a <- affinity(pH = c(3, 9), IS = IS, T = T)
- e <- equilibrate(a)
- d <- diagram(e, alpha = TRUE, plot.it = FALSE)
- alphas <- do.call(rbind, d$plotvals)
- NH <- alphas * c(0:4, 0, 1, 2, 0)
- lines(a$vals[[1]], colSums(NH), lty = 7 - pMg, col = 7 - pMg)
-}
-plot(c(3, 9), c(0, 2), type = "n", xlab = "pH", ylab = Hlab)
-lapply(2:6, Hplot)
-legend("topright", legend = paste("pMg = ", 2:6), lty = 5:1, col = 5:1, cex = 0.8)
-ATP.H.Mg <- substitute("ATP and H"^`+`~"and Mg"^`+2`)
-title(main = ATP.H.Mg)
-Mgplot <- function(pH, IS = 0.25) {
- basis("pH", pH)
- a <- affinity(`Mg+2` = c(-2, -7), IS = IS, T = T)
- e <- equilibrate(a)
- d <- diagram(e, alpha = TRUE, plot.it = FALSE)
- alphas <- do.call(rbind, d$plotvals)
- NMg <- alphas * species()$`Mg+`
- lines(-a$vals[[1]], colSums(NMg), lty = 10 - pH, col = 10 - pH)
-}
-Mglab <- substitute(italic(N)[Mg^`+2`])
-plot(c(2, 7), c(0, 1.2), type = "n", xlab = "pMg", ylab = Mglab)
-lapply(3:9, Mgplot)
-legend("topright", legend = paste("pH = ", 3:9), lty = 7:1, col = 7:1, cex = 0.8)
-title(main = ATP.H.Mg)
-```
-
-Note that we have saved the numeric results of <span style="color:green">`diagram()`</span>, i.e. the degrees of formation of the species (α).
-With that, we can calculate and plot the average number of protons bound per ATP molecule.
-To do so, we use R's `rbind()` and `do.call()` to turn `alpha` into a matrix, then multiply by the number of protons bound to each species, and sum the columns to get the total (i.e. average proton number, *N*<sub>H<sup>+</sup></sub>):
-```{r ATP, eval=FALSE, echo=8:11}
-```
-
-Adding the `IS` argument to <span style="color:green">`affinity()`</span>, we can now plot *N*<sub>H<sup>+</sup></sub> at the given ionic strength.
-Here we set `plot.it = FALSE` in `diagram()` because we use the computed α to make our own plot.
-This is similar to Figure 1.3 of Alberty (2003), but at higher temperature:
-```{r ATP, eval=FALSE, echo=12:17}
-```
-
-Next, we add the Mg<sup>+2</sup>-complexed ATP species:
-```{r ATP, eval=FALSE, echo=21}
-```
-
-Here is a function to calculate and plot *N*<sub>H<sup>+</sup></sub> for a given pMg:
-```{r ATP, eval=FALSE, echo=22:30}
-```
-
-With that function in hand, we plot the lines corresponding to pMg = 2 to 6.
-This is similar to Figure 1.4 of Alberty (2003):
-```{r ATP, eval=FALSE, echo=31:32}
-```
-
-The next function calculates and plots the average number of Mg<sup>+2</sup> bound to ATP (*N*<sub>Mg<sup>+2</sup></sub>) for a given pH.
-Here we multiply `alpha` by the number of Mg<sup>+2</sup> in each species, and negate log*a*<sub>Mg<sup>+2</sup></sub> (the variable used in <span style="color:green">`affinity()`</span>) to get pMg:
-```{r ATP, eval=FALSE, echo=36:44}
-```
-
-Using that function, we plot the lines corresponding to pH = 3 to 9.
-This is similar to Figure 1.5 of Alberty (2003):
-```{r ATP, eval=FALSE, echo=45:47}
-```
-
-<p>
-```{r ATP, fig.fullwidth=TRUE, fig.width=10, fig.height=2.5, dpi=ifelse(dpi==50, 50, 100), out.width="100%", echo=FALSE, message=FALSE, results="hide", fig.cap="Binding of H<sup>+</sup> and Mg<sup>+2</sup> to ATP at 100 °C and *I* = 0 M (first plot) or *I* = 0.25 M (third and fourth plots).", cache=TRUE, pngquant=pngquant, timeit=timeit}
-```
-</p>
-
-We have calculated the distribution of ATP species and average binding number of H<sup>+</sup> and Mg<sup>+2</sup> for given pH, pMg, ionic strength, and temperature.
-Accounting for the distribution of chemical species lends itself to thermodynamic models for reactions between reactants that have multiple ionized and complexed states.
-In contrast, Alberty (2003) and others propose models for biochemical reactions where the ionized and complexed species are combined into a single representation.
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/chnosz -r 267
More information about the CHNOSZ-commits
mailing list