[CHNOSZ-commits] r97 - in pkg/CHNOSZ: . R inst man vignettes
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sun Nov 8 11:54:56 CET 2015
Author: jedick
Date: 2015-11-08 11:54:55 +0100 (Sun, 08 Nov 2015)
New Revision: 97
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/R/protein.info.R
pkg/CHNOSZ/R/util.plot.R
pkg/CHNOSZ/inst/NEWS
pkg/CHNOSZ/man/protein.info.Rd
pkg/CHNOSZ/vignettes/equilibrium.Rnw
pkg/CHNOSZ/vignettes/equilibrium.lyx
pkg/CHNOSZ/vignettes/vig.bib
Log:
update vignette equilibrium.Rnw
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2015-11-07 08:53:26 UTC (rev 96)
+++ pkg/CHNOSZ/DESCRIPTION 2015-11-08 10:54:55 UTC (rev 97)
@@ -1,6 +1,6 @@
-Date: 2015-11-07
+Date: 2015-11-08
Package: CHNOSZ
-Version: 1.0.6-2
+Version: 1.0.6-3
Title: Chemical Thermodynamics and Activity Diagrams
Author: Jeffrey Dick
Maintainer: Jeffrey Dick <j3ffdick at gmail.com>
Modified: pkg/CHNOSZ/R/protein.info.R
===================================================================
--- pkg/CHNOSZ/R/protein.info.R 2015-11-07 08:53:26 UTC (rev 96)
+++ pkg/CHNOSZ/R/protein.info.R 2015-11-08 10:54:55 UTC (rev 97)
@@ -121,7 +121,7 @@
if(round.it) pf <- round(pf, 3)
}
# convert each protein formula to a single line
- formula <- as.chemical.formula(pf)
+ formula <- as.chemical.formula(round(pf, 3))
if(round.it) {
length <- round(length, 1)
G <- round(G, 3)
@@ -161,7 +161,7 @@
return(sb)
}
-protein.equil <- function(protein, T=25, loga.protein=0) {
+protein.equil <- function(protein, T=25, loga.protein=0, digits=4) {
# show the individual steps in calculating metastable equilibrium among proteins
msgout("protein.equil: temperature from argument is ", T, " degrees C\n")
TK <- convert(T, "K")
@@ -196,30 +196,30 @@
G0prot <- unlist(suppressMessages(subcrt(pname, T=T, property="G")$out))
# standard Gibbs energy of formation reaction of nonionized protein, cal/mol
G0protform <- G0prot - G0basissum
- msgout("protein.equil [1]: reaction to form nonionized protein from basis species has G0(cal/mol) of ", G0protform[1], "\n")
+ msgout("protein.equil [1]: reaction to form nonionized protein from basis species has G0(cal/mol) of ", signif(G0protform[1], digits), "\n")
if(ionize.it) {
# standard Gibbs energy of ionization of protein, cal/mol
G0ionization <- suppressMessages(ionize.aa(aa, property="G", T=T, pH=pH))[1, ]
- msgout("protein.equil [1]: ionization reaction of protein has G0(cal/mol) of ", G0ionization[1], "\n")
+ msgout("protein.equil [1]: ionization reaction of protein has G0(cal/mol) of ", signif(G0ionization[1], digits), "\n")
# standard Gibbs energy of formation reaction of ionized protein, cal/mol
G0protform <- G0protform + G0ionization
}
# standard Gibbs energy of formation reaction of non/ionized residue equivalents, dimensionless
G0res.RT <- G0protform/thermo$opt$R/TK/plength
- msgout("protein.equil [1]: per residue, reaction to form ", iword, " protein from basis species has G0/RT of ", G0res.RT[1], "\n")
+ msgout("protein.equil [1]: per residue, reaction to form ", iword, " protein from basis species has G0/RT of ", signif(G0res.RT[1], digits), "\n")
# coefficients of basis species in formation reactions of residues
resbasis <- suppressMessages(protein.basis(aa, T=T, normalize=TRUE))
# logQstar and Astar/RT
logQstar <- colSums(t(resbasis) * - thermo$basis$logact)
- msgout("protein.equil [1]: per residue, logQstar is ", logQstar[1], "\n")
+ msgout("protein.equil [1]: per residue, logQstar is ", signif(logQstar[1], digits), "\n")
Astar.RT <- -G0res.RT - log(10)*logQstar
- msgout("protein.equil [1]: per residue, Astar/RT = -G0/RT - 2.303logQstar is ", Astar.RT[1], "\n")
+ msgout("protein.equil [1]: per residue, Astar/RT = -G0/RT - 2.303logQstar is ", signif(Astar.RT[1], digits), "\n")
if(!is.numeric(protein)) msgout("protein.equil [1]: not comparing calculations with affinity() because 'protein' is not numeric\n")
else {
# for **Astar** we have to set the activities of the proteins to zero, not loga.protein!
a <- suppressMessages(affinity(iprotein=protein, T=T, loga.protein=0))
aAstar.RT <- log(10) * as.numeric(a$values) / plength
- msgout("check it! per residue, Astar/RT calculated using affinity() is ", aAstar.RT[1], "\n")
+ msgout("check it! per residue, Astar/RT calculated using affinity() is ", signif(aAstar.RT[1], digits), "\n")
if(!isTRUE(all.equal(Astar.RT, aAstar.RT, check.attributes=FALSE)))
stop("Bug alert! The same value for Astar/RT cannot be calculated manually as by using affinity()")
}
@@ -227,36 +227,36 @@
else {
## next set of output: equilibrium calculations
msgout("protein.equil [all]: lengths of all proteins are ", paste(plength, collapse=" "), "\n")
- msgout("protein.equil [all]: Astar/RT of all residue equivalents are ", paste(Astar.RT, collapse=" "), "\n")
+ msgout("protein.equil [all]: Astar/RT of all residue equivalents are ", paste(signif(Astar.RT, digits), collapse=" "), "\n")
expAstar.RT <- exp(Astar.RT)
sumexpAstar.RT <- sum(expAstar.RT)
- msgout("protein.equil [all]: sum of exp(Astar/RT) of all residue equivalents is ", sumexpAstar.RT, "\n")
+ msgout("protein.equil [all]: sum of exp(Astar/RT) of all residue equivalents is ", signif(sumexpAstar.RT, digits), "\n")
# boltzmann distribution
alpha <- expAstar.RT / sumexpAstar.RT
- msgout("protein.equil [all]: equilibrium degrees of formation (alphas) of residue equivalents are ", paste(alpha, collapse=" "), "\n")
+ msgout("protein.equil [all]: equilibrium degrees of formation (alphas) of residue equivalents are ", paste(signif(alpha, digits), collapse=" "), "\n")
# check with equilibrate()
if(is.numeric(protein)) {
loga.equil.protein <- unlist(suppressMessages(equilibrate(a, normalize=TRUE))$loga.equil)
# here we do have to convert from logarithms of activities of proteins to degrees of formation of residue equivalents
a.equil.residue <- plength*10^loga.equil.protein
ealpha <- a.equil.residue/sum(a.equil.residue)
- msgout("check it! alphas of residue equivalents from equilibrate() are ", paste(ealpha, collapse=" "), "\n")
+ msgout("check it! alphas of residue equivalents from equilibrate() are ", paste(signif(ealpha, digits), collapse=" "), "\n")
if(!isTRUE(all.equal(alpha, ealpha, check.attributes=FALSE)))
stop("Bug alert! The same value for alpha cannot be calculated manually as by using equilibrate()")
}
# total activity of residues
loga.residue <- log10(sum(plength * 10^loga.protein))
- msgout("protein.equil [all]: for activity of proteins equal to 10^", loga.protein, ", total activity of residues is 10^", loga.residue, "\n")
+ msgout("protein.equil [all]: for activity of proteins equal to 10^", signif(loga.protein, digits), ", total activity of residues is 10^", signif(loga.residue, digits), "\n")
# equilibrium activities of residues
loga.residue.equil <- log10(alpha*10^loga.residue)
- msgout("protein.equil [all]: log10 equilibrium activities of residue equivalents are ", paste(loga.residue.equil, collapse=" "), "\n")
+ msgout("protein.equil [all]: log10 equilibrium activities of residue equivalents are ", paste(signif(loga.residue.equil, digits), collapse=" "), "\n")
# equilibrium activities of proteins
loga.protein.equil <- log10(10^loga.residue.equil/plength)
- msgout("protein.equil [all]: log10 equilibrium activities of proteins are ", paste(loga.protein.equil, collapse=" "), "\n")
+ msgout("protein.equil [all]: log10 equilibrium activities of proteins are ", paste(signif(loga.protein.equil, digits), collapse=" "), "\n")
# check with equilibrate()
if(is.numeric(protein)) {
eloga.protein.equil <- unlist(suppressMessages(equilibrate(a, loga.balance=loga.residue, normalize=TRUE))$loga.equil)
- msgout("check it! log10 eq'm activities of proteins from equilibrate() are ", paste(eloga.protein.equil, collapse=" "), "\n")
+ msgout("check it! log10 eq'm activities of proteins from equilibrate() are ", paste(signif(eloga.protein.equil, digits), collapse=" "), "\n")
if(!isTRUE(all.equal(loga.protein.equil, eloga.protein.equil, check.attributes=FALSE)))
stop("Bug alert! The same value for log10 equilibrium activities of proteins cannot be calculated manually as by using equilibrate()")
}
Modified: pkg/CHNOSZ/R/util.plot.R
===================================================================
--- pkg/CHNOSZ/R/util.plot.R 2015-11-07 08:53:26 UTC (rev 96)
+++ pkg/CHNOSZ/R/util.plot.R 2015-11-08 10:54:55 UTC (rev 97)
@@ -152,13 +152,14 @@
ylim <- pu[3:4]
# exact lines
# warning: Eh calculations are reliable only at a single T
- if(xaxis=='pH' & (yaxis=='Eh' | yaxis=='O2' | yaxis=="pe")) {
+ if(xaxis=="O2" | (xaxis=='pH' & (yaxis=='Eh' | yaxis=='O2' | yaxis=="pe"))) {
if('reduction' %in% which) {
logfH2 <- 0
logK <- subcrt(c("H2O", "O2", "H2"), c(-1, 0.5, 1), c("liq", O2state, "gas"), T=T, P=P, convert=FALSE)$out$logK
# this is logfO2 if O2state=="gas", or logaO2 if O2state=="aq"
logfO2 <- 2 * logK - logfH2 + 2 * logaH2O
- if(yaxis=='O2') abline(h=logfO2,lty=lty,lwd=lwd,col=col)
+ if(xaxis=='O2') abline(v=logfO2,lty=lty,lwd=lwd,col=col)
+ else if(yaxis=='O2') abline(h=logfO2,lty=lty,lwd=lwd,col=col)
else if(yaxis=="Eh") lines(xlim,convert(logfO2,'E0',T=T,P=P,pH=xlim),lty=lty,lwd=lwd,col=col)
else if(yaxis=="pe") lines(xlim,convert(convert(logfO2,'E0',T=T,P=P,pH=xlim),"pe",T=T),lty=lty,lwd=lwd,col=col)
}
@@ -167,6 +168,7 @@
logK <- subcrt(c("O2", "O2"), c(-1, 1), c("gas", O2state), T=T, P=P, convert=FALSE)$out$logK
# this is logfO2 if O2state=="gas", or logaO2 if O2state=="aq"
logfO2 <- logfO2 + logK
+ if(xaxis=='O2') abline(v=logfO2,lty=lty,lwd=lwd,col=col)
if(yaxis=='O2') abline(h=logfO2,lty=lty,lwd=lwd,col=col)
else if(yaxis=="Eh") lines(xlim,convert(logfO2,'E0',T=T,P=P,pH=xlim),lty=lty,lwd=lwd,col=col)
else if(yaxis=="pe") lines(xlim,convert(convert(logfO2,'E0',T=T,P=P,pH=xlim),"pe",T=T),lty=lty,lwd=lwd,col=col)
Modified: pkg/CHNOSZ/inst/NEWS
===================================================================
--- pkg/CHNOSZ/inst/NEWS 2015-11-07 08:53:26 UTC (rev 96)
+++ pkg/CHNOSZ/inst/NEWS 2015-11-08 10:54:55 UTC (rev 97)
@@ -1,4 +1,4 @@
-CHANGES IN CHNOSZ 1.0.6-2 (2015-11-07)
+CHANGES IN CHNOSZ 1.0.6-3 (2015-11-08)
--------------------------------------
- Add functions usrfig() (get figure limits in user coordinates) and
@@ -13,9 +13,9 @@
- diagram() now has '...' argument to pass additional options to
plot() (useful with diagram(tplot=FALSE, ...)).
-- TODO: Update vignette 'equilibrium.Rnw' with better definitions of
- concepts, organization of functions, and amino acid and protein
- examples.
+- Update vignette 'equilibrium.Rnw' with better definitions of concepts,
+ organization of functions, and examples and applications. Now uses
+ knitr.
- Remove demo diagram.R (concepts better shown in equilibrium.Rnw).
Modified: pkg/CHNOSZ/man/protein.info.Rd
===================================================================
--- pkg/CHNOSZ/man/protein.info.Rd 2015-11-07 08:53:26 UTC (rev 96)
+++ pkg/CHNOSZ/man/protein.info.Rd 2015-11-08 10:54:55 UTC (rev 97)
@@ -18,7 +18,7 @@
protein.length(protein, organism = NULL)
protein.info(protein, T = 25, residue = FALSE, round.it = FALSE)
protein.basis(protein, T = 25, normalize = FALSE)
- protein.equil(protein, T=25, loga.protein = 0)
+ protein.equil(protein, T=25, loga.protein = 0, digits = 4)
MP90.cp(protein, T)
group.formulas()
}
@@ -31,6 +31,7 @@
\item{T}{numeric, temperature in \eqn{^{\circ}}{°}C}
\item{round.it}{logical, round the values in the output?}
\item{loga.protein}{numeric, decimal logarithms of reference activities of proteins}
+ \item{digits}{integer, number of significant digits (see \code{\link{signif}})}
}
\details{
Modified: pkg/CHNOSZ/vignettes/equilibrium.Rnw
===================================================================
--- pkg/CHNOSZ/vignettes/equilibrium.Rnw 2015-11-07 08:53:26 UTC (rev 96)
+++ pkg/CHNOSZ/vignettes/equilibrium.Rnw 2015-11-08 10:54:55 UTC (rev 97)
@@ -1,32 +1,30 @@
-%% LyX 2.1.3 created this file. For more info, see http://www.lyx.org/.
+%% LyX 2.1.4 created this file. For more info, see http://www.lyx.org/.
%% Do not edit unless you really know what you are doing.
-\documentclass[english,noae]{article}
+\documentclass[english,round]{article}
\usepackage{mathpazo}
+\renewcommand{\sfdefault}{lmss}
+\renewcommand{\ttdefault}{lmtt}
\usepackage[T1]{fontenc}
\usepackage[latin9]{inputenc}
-\usepackage[letterpaper]{geometry}
+\usepackage{geometry}
\geometry{verbose,tmargin=2.5cm,bmargin=2.5cm,lmargin=2.5cm,rmargin=2.5cm}
\usepackage{color}
\usepackage{babel}
\usepackage{amsbsy}
\usepackage{amssymb}
-\usepackage[numbers]{natbib}
+\usepackage[authoryear]{natbib}
\usepackage[unicode=true,pdfusetitle,
bookmarks=true,bookmarksnumbered=false,bookmarksopen=false,
- breaklinks=false,pdfborder={0 0 1},backref=false,colorlinks=true]
+ breaklinks=false,pdfborder={0 0 0},backref=false,colorlinks=true]
{hyperref}
\hypersetup{
- citecolor=black, linkcolor=black}
+ citecolor=magenta, urlcolor=blue}
\usepackage{breakurl}
\makeatletter
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Textclass specific LaTeX commands.
-<<echo=F>>=
- if(exists(".orig.enc")) options(encoding = .orig.enc)
-@
-
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% User specified LaTeX commands.
%\VignetteIndexEntry{Equilibrium in CHNOSZ}
+%\VignetteEngine{knitr::knitr}
% so DOIs in bibliography show up as hyperlinks
\newcommand*{\doi}[1]{\href{http://dx.doi.org/#1}{doi: #1}}
@@ -34,62 +32,683 @@
\makeatother
\begin{document}
+<<setup, include=FALSE, cache=FALSE>>=
+library(knitr)
+## set global chunk options
+opts_chunk$set(fig.align='center')
+opts_chunk$set(tidy=TRUE)
+## set code/output width
+options(width=85)
+## set number of digits
+options(digits=3)
+@
+
\title{Equilibrium in CHNOSZ}
\author{Jeffrey M. Dick}
\maketitle
+This document defines the concepts, explains the organization of functions,
+and provides examples of calculating equilibrium in CHNOSZ. It also
+highlights some applications of the methods (i.e. to reproduce published
+diagrams) and includes an Appendix on details of the equilibration
+calculations.
-\section{Background; terminology}
-\emph{Species of interest} (or simply \emph{species}) are chosen by
-the user to represent the possible species in a chemical system. \emph{Basis
-species} combine linearly to make up the species of interest. Basis
-species are analogous to thermodynamic components in that they are
-the minimum number to describe the compositional variation, but unlike
-components basis species can be charged.
+\section{Concepts}
+\begin{description}
+\item [{Species~of~interest}] Chemical species for which you want to
+calculate relative stabilities.
+\item [{Basis~species}] Species in terms of which you want to write all
+formation reactions of species of interest.
+\item [{Formation~reactions}] Stoichiometric chemical reactions showing
+the mass requirements for formation of 1 mole of each species of interest
+from the basis species.
+\item [{Chemical~affinity}] Change in Gibbs energy per increment of reaction
+progress; the negative of the overall Gibbs energy of reaction; $\boldsymbol{A}=2.303RT\log(K/Q)$,
+where $K$ is the equilibrium constant and $Q$ is the reaction activity
+quotient ($\log$ here is used for base-10 logarithms).
+\item [{(1)~Reference~activity}] User-defined (usually equal) activities
+of species of interest.
+\item [{(1)~Reference~affinity}] ($\boldsymbol{A}_{{\it ref}}$) Chemical
+affinity of formation reaction with a reference activity of the species
+of interest.
+\item [{(1)~Maximum~affinity~method}] Comparison of reference affinities
+for given balance coefficients in order to calculate stability regions
+on a predominance diagram.
+\item [{Balance~coefficients}] ($n_{balance}$) The number of moles of
+a basis species present in the formation reaction of each of the species
+of interest. Reactions between any two species of interest then are
+``balanced'' on this basis species. Can be a quantity other than
+basis species (e.g., balance = 1, or length of amino acid sequence
+of protein).
+\item [{Predominance~diagram}] Diagram showing fields of maximal stability
+(i.e. greatest activity at equilibrium) for species of interest as
+a function of two variables (aka equal activity diagram).
+\item [{(2)~Starred~affinity}] ($\boldsymbol{A}^{*}$) Chemical affinity
+of formation reaction with unit activity of the species of interest
+(aka ``starved'' affinity because the activity of the species of
+interest drops out of $Q$).
+\item [{(2)~Total~balance~activity}] The sum of activities of this basis
+species contributed by each of the species of interest. (In Appendix:
+activity of the immobile or conserved component; $a_{\mathrm{ic}}$.)
+\item [{(2)~Equilibration~method}] Comparison of starred affinities in
+order to calculate activities of species of interest for given balance
+coefficients and total balance activity.
+\item [{Speciation~diagram}] Diagram showing the activities of species
+of interest, usually as a function of 1 variable (aka activity diagram).
+\item [{Boltzmann~distribution}] Algorithm used for the equilibration
+method when the balance coefficients are 1.
+\item [{Reaction~matrix}] Algorithm used for the equilibration method
+when the balance coefficients are not all 1.
+\item [{Normalization}] Algorithm used for large molecules such as proteins;
+chemical formulas and affinities are scaled to an equal size (e.g.
+a single residue; ``residue equivalent'' in Appendix), activities
+are calculated using balance = 1, and formulas and activities are
+rescaled to the original size of the molecule.
+\item [{Mosaic}] Calculations of chemical affinities for making diagrams
+where the speciation of basis species depends on the variables.
+\end{description}
+The numbered groups above are connected with two distinct approaches
+to generating diagrams:
+\begin{enumerate}
+\item With the \textbf{maximum affinity method} for creating predominance
+diagrams, the user sets the reference activities of the species of
+interest; the program compares the reference affinities at these conditions
+to determine the most stable species (highest activity, i.e. predominant
+at equilibrium).
+\item With the \textbf{equilibration method} for creating predominance or
+activity diagrams, the user explicitly sets the total balance activity
+or the program takes it from the reference activities of the species.
+The starred affinities are used to calculate equilibrium activities
+using one of two techniques (Boltzmann distribution for balance =
+1, reaction matrix for balance $\ne$ 1).
+\end{enumerate}
+The affinities used in these calculations can be calculated using
+\texttt{affinity()}, which works with a single basis set, or with
+\texttt{mosaic()}, which uses multiple basis sets to account for basis
+species that themselves may change as a function of the variables
+of interest (e.g. ionization of carbonic acid as a function of pH).
+This document focuses primarily on the \texttt{affinity()} function;
+for more information on mosaic diagrams see the help page (type \texttt{?mosaic}
+at the R command line).
-The calculations of chemical equilibrium in CHNOSZ are formulated
-for a system that is open with respect to the basis species. Therefore,
-the natural variables are temperature, pressure and chemical potentials
-of the basis species.
+Step-by-step examples of some of the calculations, particularly the
+reaction matrix algorithm, are provided in the Appendix. For further
+description of the equilibration method applied to proteins see \citet{DS13}
+(also with a derivation of energetic distance from equilibrium using
+the \textbf{starred affinity}).
-To calculate the equilibrium distribution of species in a given system,
-a single linear balancing constraint must be specified. Therefore,
-we say things like ``the reactions are balanced on $\mathrm{CO_{2}}$''
-or ``the reactions are balanced on protein length''. The \emph{balancing
-coefficients} describe these constraints. At different times in the
-documentation of CHNOSZ, the balancing constraint has been associated
-with the conserved component, immobile component, or conserved basis
-species, all with the same meaning.
-The ``reactions'' in the preceding statements refer to \emph{transformations}
-between species. The actual calculations, however, start with the
-definitions of the formation reactions. The \emph{formation reaction}
-for any species has one mole of the species as a product, and the
-mass balance is made up of the basis species.
+\section{Organization}
-By definition, the formation reaction is written to form one mole
-of a species. For many systems, the extent of the molar formula of
-the species is not a matter of concern. However, for systems made
-of polymers, such as proteins, it is often desirable to normalize
-the molar formula by the balancing coefficients. This normalization
-has been referred to previously as using the residue equivalents of
-proteins \citep{Dic08}; here the terminology of \emph{normalize}
-will be used preferentially\footnote{The older style of function call using \texttt{diagram(..., residue=TRUE)}
-has been replaced by \texttt{equilibrate(..., normalize=TRUE)} or
-\texttt{diagram(..., normalize=TRUE)} starting with version 0.9-9
-of CHNOSZ.}.
+The function sequences below assume you have already defined the basis
+species and species of interest using \texttt{basis(...)} and \texttt{species(...)}
+(ellipses here and below indicate system-specific input).
+Note that if \texttt{equilibrate()} or \texttt{diagram()} is called
+without an explicit \texttt{balance} argument, the balance coefficients
+will be taken from the first basis species (in the current basis definition)
+that is present in all of the species. Depending on the system, this
+may coincide either with balance = 1 or with balance $\ne$ 1. In
+the case of \texttt{normalize = TRUE} or \texttt{as.residue = TRUE},
+the balance coefficients (for the purposes of the equilibration step)
+are temporarily set to 1.
+\begin{enumerate}
+\item Maximum affinity method, balance = 1
+
+\begin{enumerate}
+\item Typical use: simple mineral/aqueous species stability comparisons
+\item Function sequence:\\
+\texttt{a <- affinity(...)}\\
+\texttt{diagram(a, balance = 1)}
+\item Algorithm: $\max\left\{ \boldsymbol{A}_{{\it ref}}\right\} $
+\end{enumerate}
+\item Equilibration method, balance = 1
+
+\begin{enumerate}
+\item Typical use: simple aqueous species activity comparisons
+\item Function sequence:\\
+\texttt{a <- affinity(...)}\\
+\texttt{e <- equilibrate(e, balance = 1)}\\
+\texttt{diagram(e)}
+\item Algorithm: Boltzmann distribution
+\end{enumerate}
+\item Maximum affinity method, balance $\ne$ 1
+
+\begin{enumerate}
+\item Typical use: mineral/aqueous species stability comparisons
+\item Function sequence:\\
+\texttt{a <- affinity(...)}~\\
+\texttt{diagram(a, balance = ...)}
+\item Algorithm: $\max\left\{ \boldsymbol{A}_{{\it ref}}/n_{balance}\right\} $
+\end{enumerate}
+\item Equilibration method, balance $\ne$ 1
+
+\begin{enumerate}
+\item Typical use: aqueous species activity comparisons
+\item Function sequence:\\
+\texttt{a <- affinity(...)}~\\
+\texttt{e <- equilibrate(a, balance = ...)}~\\
+\texttt{diagram(e)}
+\item Algorithm: Reaction matrix
+\end{enumerate}
+\item Maximum affinity method, normalize = TRUE
+
+\begin{enumerate}
+\item Typical use: protein/polymer stability comparisons
+\item Function sequence:\\
+\texttt{a <- affinity(...)}~\\
+\texttt{diagram(a, normalize = TRUE)}
+\item Algorithm: $\max\left\{ \boldsymbol{A}^{*}/n_{balance}-\log n_{balance}\right\} $
+\end{enumerate}
+\item Equilibration method, normalize = TRUE
+
+\begin{enumerate}
+\item Typical use: protein/polymer activity comparisons
+\item Function sequence:\\
+\texttt{a <- affinity(...)}~\\
+\texttt{e <- equilibrate(a, normalize = TRUE)}~\\
+\texttt{diagram(e)}
+\item Algorithm: Scale formulas and affinities to residues; Boltzmann distribution
+(balance = 1); Scale activities to proteins
+\end{enumerate}
+\end{enumerate}
+
+\section{Examples}
+
+
+\subsection{Amino acids}
+
+Basis species: $\mathrm{CO_{2}}$, $\mathrm{H_{2}O}$, $\mathrm{NH_{3}}$,
+$\mathrm{H_{2}S}$, $\mathrm{O_{2}}$. Species of interest: 20 amino
+acids. (Only the first few lines of the data frame of amino acid species
+are shown.)
+
+<<AAsetup>>=
+library(CHNOSZ)
+data(thermo)
+basis("CHNOS")
+species(aminoacids(""))[1:5, ]
+@
+
+Code for figures. Function names are keyed to the subfigures.
+
+<<AAfigures>>=
+res <- 200
+aa <- aminoacids()
+
+aaA <- function() {
+ a <- affinity(O2=c(-90, -70, res), H2O=c(-20, 10, res))
+ diagram(a, balance=1, names=aa)
+}
+
+aaB <- function() {
+ a <- affinity(O2=c(-90, -70, res), H2O=c(-20, 10, res))
+ e <- equilibrate(a, balance=1)
+ diagram(e, names=aa)
+}
+
+aaC <- function() {
+ a <- affinity(O2=c(-71, -66, res), H2O=c(-8, 4, res))
+ diagram(a, balance="CO2", names=aa)
+}
+
+aaD <- function() {
+ a <- affinity(O2=c(-71, -66), H2O=c(-8, 4))
+ e <- equilibrate(a, balance="CO2")
+ diagram(e, names=aa)
+}
+
+aaE <- function() {
+ basis("O2", -66)
+ a <- affinity(H2O=c(-8, 4))
+ e <- equilibrate(a, balance="CO2")
+ diagram(e, ylim=c(-5, -1), names=aa)
+}
+
+aaF <- function() {
+ species(1:20, -4)
+ a <- affinity(H2O=c(-8, 4))
+ e <- equilibrate(a, balance="CO2")
+ diagram(e, ylim=c(-5, -1), names=aa)
+}
+@
+
+Note that for the plot we use the 1-letter abbreviations of the amino
+acids, unlike the full species names (\texttt{aminoacids()} is a function
+in CHNOSZ that returns their names or abbreviations).
+
+<<AAnames>>=
+AA <- aminoacids("")
+names(AA) <- aa
+AA
+@
+
+The annotated figure is shown next. The actual code used to set up
+the plots, add labels, etc. is in the source of this vignette (not
+shown in the PDF).
+
+<<showtime, include=FALSE>>=
+showtime <- function(st) {
+ # plot time in lower-right of figure region
+ f <- usrfig()
+ par(xpd=TRUE)
+ if(st[3] > 2) col <- "red" else col <- "black"
+ text(f$x[2], f$y[1], paste(round(st[3], 1), "s\n"), adj=1, col=col)
+ par(xpd=FALSE)
+}
+@
+
+<<AAplot, echo=FALSE, message=FALSE, results="hide", fig.width=13/2, fig.height=9/2, cache=TRUE>>=
+layout(t(matrix(c(1:7, 11, 8:10, 12), nrow=4)), widths=c(1, 4, 4, 4), heights=c(1, 4, 4))
+
+## row 0 (column titles)
+opar <- par(mar=c(0, 0, 0, 0))
+plot.new()
+plot.new()
+text(0.5, 0.5, "maximum affinity", cex=1.4)
+plot.new()
+text(0.5, 0.5, "equilibration", cex=1.4)
+plot.new()
+text(0.5, 0.5, "equilibration", cex=1.4)
+par(opar)
+
+## row 1 (balance = 1)
+opar <- par(mar=c(0, 0, 0, 0))
+plot.new()
+text(0.5, 0.5, "balance = 1", srt=90, cex=1.4)
+par(opar)
+# figure A
+st <- system.time(dA <- aaA())
+showtime(st)
+title(main="loga(species) = -3", cex.main=1)
+label.figure("A", col="blue", yfrac=0.9, xfrac=0.1)
+# figure B
+st <- system.time(dB <- aaB())
+showtime(st)
+title(main=paste("loga(total species) =", round(dB$loga.balance, 2)), cex.main=1)
+label.figure("B", col="blue", yfrac=0.9, xfrac=0.1)
+
+## row 2 (balance = nCO2)
+opar <- par(mar=c(0, 0, 0, 0))
+plot.new()
+text(0.5, 0.5, 'balance = "CO2"', srt=90, cex=1.4)
+par(opar)
+# figure C
+st <- system.time(dC <- aaC())
+showtime(st)
+title(main="loga(species) = -3", cex.main=1)
+label.figure("C", col="blue", yfrac=0.9, xfrac=0.1)
+# figure D
+st <- system.time(dD <- aaD())
+showtime(st)
+title(main=paste("loga(total CO2) =", round(dD$loga.balance, 2)), cex.main=1)
+label.figure("D", col="blue", yfrac=0.9, xfrac=0.1)
+
+## right (speciation at different total activity of CO2)
+par(xpd=NA)
+lines(c(-66, -64.5), c(4, 9), lty=2)
+lines(c(-66, -64.5), c(-8, -8.5), lty=2)
+par(xpd=FALSE)
+# figure E
+st <- system.time(dE <- aaE())
+showtime(st)
+title(main=paste("loga(total CO2) =", round(dE$loga.balance, 2)), cex.main=1)
+label.figure("E", col="blue", yfrac=0.9, xfrac=0.1)
+# figure F
+st <- system.time(dF <- aaF())
+showtime(st)
+title(main=paste("loga(total CO2) =", round(dF$loga.balance, 2)), cex.main=1)
+label.figure("F", col="blue", yfrac=0.9, xfrac=0.1)
+
+@
+
+Comments on the plots:
+\begin{itemize}
+\item The equal-activity lines in Figures \textcolor{blue}{A} and \textcolor{blue}{B}
+are \emph{identical}. For balance = 1, the maximum affinity method
+and the equilibration method should produce the same predominance
+diagrams. (More precisely, because balance = 1, the conditions of
+equal activity of any species of interest \textbf{are independent
+of} the actual value of that activity.)
+\item Figures \textcolor{blue}{C} and \textcolor{blue}{D} are \emph{different}.
+For balance $\ne$ 1, the maximum affinity method and the equilibration
+method will generally produce difference predominance diagrams. (Because
+balance $\ne$ 1, the conditions of equal activity of any species
+of interest \textbf{depend on} the actual value of that activity.)
+\item Both Figures \textcolor{blue}{E} and \textcolor{blue}{F} are constructed
+using the equilibration method, to calculate activities of species
+as a function of $\log a_{\mathrm{H_{2}O}}$ at $\log f_{\mathrm{O_{2}}}=-66$.
+Figure \textcolor{blue}{E} shows the results for the default settings
+($a_{\mathrm{CO_{2}}}$ is the sum of activities present in all species,
+taken from initial species activity of $10^{-3}$) and the crossing
+lines indicating equal activities \emph{are identical to the positions
+in Figure }\textcolor{blue}{\emph{D}} at $\log f_{\mathrm{O_{2}}}=-66$.
+\item Figure \textcolor{blue}{F} shows the results for a lower total activity
+of $\mathrm{CO_{2}}$. Consequently, the activities of the predominant
+species decrease from ca. $10^{-2}$ in Figure \textcolor{blue}{E}
+to ca. $10^{-3}$ in Figure \textcolor{blue}{F}. Also, the stability
+region of the smaller glycine has grown at the expense of the neighboring
+bigger amino acids, so that the crossing lines indicating equal activities
+in Figure \textcolor{blue}{F} \emph{are closer to those shown in Figure
+}\textcolor{blue}{\emph{C}} at $\log f_{\mathrm{O_{2}}}=-66$.
+\item In other words, a lower equal-activity value causes the stability
+region of the species with the smaller balance coefficient to invade
+that of the species with the larger balance coefficient.
+\item Figures \textcolor{blue}{A}, \textcolor{blue}{B}, \textcolor{blue}{C},
+and \textcolor{blue}{D} are all equal activity diagrams, but have
+different constraints on the activities:
+
+\begin{itemize}
+\item Maximum affinity method (Figures \textcolor{blue}{A}, \textcolor{blue}{C}):
+Equal activities of species set to a constant value.
+\item Equilibration method (Figures \textcolor{blue}{B}, \textcolor{blue}{D}):
+Equal activities of species determined by overall speciation of the
+system.
+\end{itemize}
+\end{itemize}
+
+\subsection{Proteins}
+
+Basis species: $\mathrm{CO_{2}}$, $\mathrm{H_{2}O}$, $\mathrm{NH_{3}}$,
+$\mathrm{H_{2}S}$, $\mathrm{O_{2}}$, $\mathrm{H^{+}}$. Species
+of interest: 6 proteins from the set of archaeal and bacterial surface
+layer proteins considered by \citet{Dic08}.
+
+<<PRsetup, results="hide">>=
+basis("CHNOS+")
+organisms <- c("METJA", "HALJP", "METVO", "ACEKI", "GEOSE", "BACLI")
+proteins <- c(rep("CSG", 3), rep("SLAP", 3))
+species(proteins, organisms)
+@
+
+Code for the figures.
+
+<<PRfigures>>=
+prA <- function() {
+ a <- affinity(O2=c(-90, -70, res), H2O=c(-20, 0, res))
+ e <- equilibrate(a, balance="length", loga.balance=0)
+ diagram(e, names=organisms)
+}
+
+prB <- function() {
+ a <- affinity(O2=c(-90, -70))
+ e <- equilibrate(a, balance="length", loga.balance=0)
+ diagram(e, names=organisms, ylim=c(-5, -1), legend.x=NA)
+}
+
+prC <- function() {
+ a <- affinity(O2=c(-90, -70, res), H2O=c(-20, 0, res))
+ e <- equilibrate(a, normalize=TRUE, loga.balance=0)
+ diagram(e, names=organisms)
+}
+
+prD <- function() {
+ a <- affinity(O2=c(-90, -70))
+ e <- equilibrate(a, normalize=TRUE, loga.balance=0)
+ diagram(e, names=organisms, ylim=c(-5, -1), legend.x=NA)
+}
+
+prE <- function() {
+ a <- affinity(O2=c(-90, -70, res), H2O=c(-20, 0, res))
+ e <- equilibrate(a, as.residue=TRUE, loga.balance=0)
+ diagram(e, names=organisms)
+}
+
+prF <- function() {
+ a <- affinity(O2=c(-90, -70))
+ e <- equilibrate(a, as.residue=TRUE, loga.balance=0)
+ diagram(e, names=organisms, ylim=c(-3, 1), legend.x=NA)
+}
+@
+
+The plots follow. As before, the code used to layout the figure and
+label the plots is not shown in the PDF.
+
+<<PRplot, message=FALSE, echo=FALSE, results="hide", fig.width=13/2, fig.height=9/2, cache=TRUE>>=
+layout(t(matrix(1:12, nrow=4)), widths=c(1, 4, 4, 4), heights=c(1, 4, 4))
+
+## row 0 (column titles)
+opar <- par(mar=c(0, 0, 0, 0))
+plot.new()
+plot.new()
+text(0.5, 0.5, 'balance = "length"', cex=1.4)
+plot.new()
+text(0.5, 0.5, "normalize = TRUE\n(balance = 1)", cex=1.4)
+plot.new()
+text(0.5, 0.5, "as.residue = TRUE\n(balance = 1)", cex=1.4)
+par(opar)
+
+## row 1 (equilibrate 2D)
+opar <- par(mar=c(0, 0, 0, 0))
+plot.new()
+text(0.5, 0.5, "equilibration", srt=90, cex=1.4)
+par(opar)
+# figure A (balance = "length")
+st <- system.time(dA <- prA())
+showtime(st)
+label.figure("A", col="blue", yfrac=0.9, xfrac=0.1)
+# figure C (normalize = TRUE)
+st <- system.time(dC <- prC())
+showtime(st)
+label.figure("C", col="blue", yfrac=0.9, xfrac=0.1)
+# figure E (as.residue = TRUE)
+st <- system.time(dE <- prE())
+showtime(st)
+label.figure("E", col="blue", yfrac=0.9, xfrac=0.1)
+
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/chnosz -r 97
More information about the CHNOSZ-commits
mailing list