[CHNOSZ-commits] r10 - in pkg: . R inst inst/doc inst/extdata/thermo man vignettes
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Sep 20 03:03:24 CEST 2012
Author: jedick
Date: 2012-09-20 03:03:23 +0200 (Thu, 20 Sep 2012)
New Revision: 10
Added:
pkg/inst/extdata/thermo/DLEN67.csv
Modified:
pkg/DESCRIPTION
pkg/R/wjd.R
pkg/inst/NEWS
pkg/inst/doc/wjd.pdf
pkg/man/extdata.Rd
pkg/man/wjd.Rd
pkg/vignettes/vig.bib
pkg/vignettes/wjd.Rnw
pkg/vignettes/wjd.lyx
Log:
wjd: add 'central' method to guess(), example from Dayhoff et al. to vignette
Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION 2012-09-18 15:06:42 UTC (rev 9)
+++ pkg/DESCRIPTION 2012-09-20 01:03:23 UTC (rev 10)
@@ -1,11 +1,11 @@
-Date: 2012-09-18
+Date: 2012-09-20
Package: CHNOSZ
Version: 0.9-7.98
Title: Chemical Thermodynamics and Activity Diagrams
Author: Jeffrey M. Dick
Maintainer: Jeffrey M. Dick <jmdick at asu.edu>
Depends: R (>= 2.10.0), utils
-Suggests: testthat, parallel
+Suggests: testthat, parallel, limSolve
Description: This package includes functions and data sets to support chemical thermodynamic
modeling in biochemistry and low-temperature geochemistry. The features include calculation
of the standard molal thermodynamic properties and chemical affinities of reactions involving
Modified: pkg/R/wjd.R
===================================================================
--- pkg/R/wjd.R 2012-09-18 15:06:42 UTC (rev 9)
+++ pkg/R/wjd.R 2012-09-20 01:03:23 UTC (rev 10)
@@ -225,7 +225,7 @@
return(ep)
}
-is.near.equil <- function(w, tol=0.01) {
+is.near.equil <- function(w, tol=0.01, quiet=FALSE) {
# given the output of wjd(), make a simple test for equilibrium
# that is, that the chemical potentials of the elements are nearly
# the same when calculated using different sets of species in the system
@@ -235,6 +235,12 @@
# equilibrium unless proven guilty
ine <- TRUE
for(i in 1:ncol(ep)) if(diff(range(ep[,i])) > tol) ine <- FALSE
+ if(!ine & !quiet) {
+ # talk about the differences in chemical potentials
+ epdiff <- abs(apply(apply(ep, 2, range), 2, diff))
+ imax <- which.max(epdiff)
+ msgout("is.near.equil: solution has variation of ", epdiff[imax], " in mu/RT of ", names(epdiff)[imax], "\n")
+ }
return(ine)
}
@@ -244,85 +250,108 @@
0,0,0,1,2,1,1,0,0,0,
0,0,1,0,0,0,1,1,2,1),ncol=3,
dimnames=list(NULL,c("H","N","O"))),
- B = c(2,1,1), iguess=1, ic=NULL
+ B = c(2,1,1), method=c("central", "stoich"), minX=0.001, iguess=1, ic=NULL
){
# given the elemental stoichiometries of a set of species (A)
# and the number of moles of elements (B)
# find moles of species that satisfy mass balance and are all positive
# generally this will be one of the solutions of an underdetermined system
- # 20111231 jmd
- # approach:
- # - select one of the (many) species combinations (ic) that
- # make a square, invertible stoichiometric matrix (the "variable" species)
- # - assign equal mole numbers to all the "other" species (Xother),
- # such that any element has at most max.frac fraction of the desired amount (B)
- # (max.frac is scanned from 0.01 to 0.99)
- # - calculate the mole numbers of the stoichiometry-setting species
- # that give the desired elemental composition; accept the provisional
- # solution if all numbers are positive
-
- # arguments:
- # A - the stoichiometric matrix
- # B - the moles of elements
- # iguess - which provisional guess to return (NULL for all)
- # ic - which specific combination of species to test (NULL for all)
-
# first of all, we can't do anything if all there are no elements
if(all(B==0)) stop("there are zero moles of all elements")
- # get the various combinations of species that are
- # stoichiometrically independent
- combs <- invertible.combs(A)
- # we will potentially process all of them unless a specific one is identified
- if(is.null(ic)) ic <- 1:nrow(combs)
- # a counter to keep track of the provisional guesses
- iprov <- 0
- # where to store the output if we want all guesses
- out <- list()
- for(i in ic) {
- # which species are the variable ones
- ivar <- combs[i,]
- # moles of elements for one mole of all of the other species
- Bother.1 <- colSums(A[-ivar, , drop=FALSE])
- # which element is proportionally most highly represented w.r.t. the desired composition
- imax <- which.max(Bother.1/B)
- # max.frac - the highest fraction of contribution to moles of elements by the "other" species
- for(max.frac in c(0.01, 0.05, 0.1, 0.25, 0.5, 0.75, 0.9, 0.95, 0.99)) {
- # the number of moles of all "other" species that give max.frac of any element
- Xother <- max.frac/(Bother.1/B)[imax]
- # moles of elements for this number of moles of all the other species
- Bother <- Bother.1 * Xother
- # moles of elements that are left for the variable species
- Bvar <- B - Bother
- # now solve for the number of moles of the variable species
- Xvar <- solve(t(A[ivar,]),Bvar)
- # stop the search if we found a positive solution
- if(all(Xvar > 0)) break
+
+ # if method="central" get central solution using limSolve package 20120919
+ if("central" %in% method) {
+ if(!"limSolve" %in% row.names(installed.packages())) {
+ msgout("guess: skipping 'central' method as limSolve package is not available\n")
+ } else {
+ require(limSolve)
+ # the inequality constraints for moles of species
+ G <- diag(nrow(A))
+ # minX is the minimum mole number we will accept
+ H <- rep(minX, nrow(A))
+ # get a solution
+ X <- xranges(E=t(A), F=B, G=G, H=H, central=TRUE, full=TRUE)[, "central"]
+ return(X)
}
- # put them together
- X <- numeric(nrow(A))
- X[-ivar] <- Xother
- X[ivar] <- Xvar
- # if all the moles are positive, this is a provisional
- # guess, otherwise make the result NA
- if(any(Xvar <= 0)) X <- NA
- else iprov <- iprov + 1
- # return the result if we're at the correct guess number
- if(is.null(iguess)) out <- c(out,list(X))
- else if(iprov==iguess) return(X)
}
- # if we're here, we should return all guesses, or
- # make an error (the requested guess number doesn't exist)
- if(is.null(iguess) & iprov > 0) return(out)
- else {
- if(is.null(iguess)) iguess <- "[ALL]"
- stop(paste("you asked for guess number ",iguess,
- " but there are only ",iprov,
- " that satisfy all stoichiometric constraints",sep=""))
+
+ if("stoich" %in% method) {
+ # if method="stoich" use a stoichiometric approach: 20111231 jmd
+ # - select one of the (many) species combinations (ic) that
+ # make a square, invertible stoichiometric matrix (the "variable" species)
+ # - assign equal mole numbers to all the "other" species (Xother),
+ # such that any element has at most max.frac fraction of the desired amount (B)
+ # (max.frac is scanned from 0.01 to 0.99)
+ # - calculate the mole numbers of the stoichiometry-setting species
+ # that give the desired elemental composition; accept the provisional
+ # solution if all numbers are positive
+
+ # arguments:
+ # A - the stoichiometric matrix
+ # B - the moles of elements
+ # iguess - which provisional guess to return (NULL for all)
+ # ic - which specific combination of species to test (NULL for all)
+
+ # get the various combinations of species that are
+ # stoichiometrically independent
+ combs <- invertible.combs(A)
+ # we will potentially process all of them unless a specific one is identified
+ if(is.null(ic)) ic <- 1:nrow(combs)
+ # a counter to keep track of the provisional guesses
+ iprov <- 0
+ # where to store the output if we want all guesses
+ out <- list()
+ for(i in ic) {
+ # which species are the variable ones
+ ivar <- combs[i,]
+ # moles of elements for one mole of all of the other species
+ Bother.1 <- colSums(A[-ivar, , drop=FALSE])
+ # which element is proportionally most highly represented w.r.t. the desired composition
+ imax <- which.max(Bother.1/B)
+ # max.frac - the highest fraction of contribution to moles of elements by the "other" species
+ for(max.frac in c(0.01, 0.05, 0.1, 0.25, 0.5, 0.75, 0.9, 0.95, 0.99)) {
+ # the number of moles of all "other" species that give max.frac of any element
+ Xother <- max.frac/(Bother.1/B)[imax]
+ # moles of elements for this number of moles of all the other species
+ Bother <- Bother.1 * Xother
+ # moles of elements that are left for the variable species
+ Bvar <- B - Bother
+ # now solve for the number of moles of the variable species
+ Xvar <- solve(t(A[ivar,]),Bvar)
+ # stop the search if we found a positive solution
+ if(all(Xvar > 0)) break
+ }
+ # put them together
+ X <- numeric(nrow(A))
+ X[-ivar] <- Xother
+ X[ivar] <- Xvar
+ # add names
+ names(X) <- rownames(A)
+ # if all the moles are positive, this is a provisional
+ # guess, otherwise make the result NA
+ if(any(Xvar <= 0)) X <- NA
+ else iprov <- iprov + 1
+ # return the result if we're at the correct guess number
+ if(is.null(iguess)) out <- c(out,list(X))
+ else if(iprov==iguess) return(X)
+ }
+ # if we're here, we should return all guesses, or
+ # make an error (the requested guess number doesn't exist)
+ if(is.null(iguess) & iprov > 0) return(out)
+ else {
+ if(is.null(iguess)) iguess <- "[ALL]"
+ stop(paste("you asked for guess number ",iguess,
+ " but there are only ",iprov,
+ " that satisfy all stoichiometric constraints",sep=""))
+ }
}
+
+ # if we're here we didn't find a guessing method
+ stop("no method found")
}
-run.wjd <- function(ispecies, B=NULL, Y=run.guess(ispecies, B), P=1, T=25, imax=10, Gfrac=1e-7, tol=0.01) {
+run.wjd <- function(ispecies, B=NULL, method="stoich", Y=run.guess(ispecies, B, method), P=1, T=25, imax=10, Gfrac=1e-7, tol=0.01) {
### set up a Gibbs energy minimization
### using compositions and standard Gibbs energies of species
### from database in CHNOSZ 20120101 jmd
@@ -341,22 +370,30 @@
# a single guess
w <- wjd(A, G0.RT, Y, P=P, imax=imax, Gfrac=Gfrac)
} else {
- # loop over all the guesses created by run.guess
- Y <- Y[!is.na(Y)]
- for(i in 1:length(Y)) {
- w <- wjd(A, G0.RT, Y[[i]], P=P, imax=imax, Gfrac=Gfrac)
- if(is.near.equil(w, tol=tol)) {
- msgout("run.wjd: got apparently near equilibrium on initial solution ", i, " of ", length(Y), "\n")
- break
- } else if(i==length(Y)) {
- stop(paste("couldn't get apparently near equilibrium after", length(Y), "initial solutions"))
+ # if we're using method "central" there is only one guess
+ if(method=="central") {
+ w <- wjd(A, G0.RT, Y, P=P, imax=imax, Gfrac=Gfrac)
+ } else {
+ # for method "stoich" loop over all the guesses created by run.guess
+ Y <- Y[!is.na(Y)]
+ for(i in 1:length(Y)) {
+ w <- wjd(A, G0.RT, Y[[i]], P=P, imax=imax, Gfrac=Gfrac)
+ if(is.near.equil(w, tol=tol)) {
+ msgout("run.wjd: got within tolerance on initial solution ", i, " of ", length(Y), "\n")
+ break
+ }
+ if(i==length(Y)) msgout("run.wjd: tested ", length(Y), " initial solutions\n")
}
}
+ # only return a near equilibrium solution
+ if(!is.near.equil(w, tol=tol)) {
+ stop(paste("couldn't find a solution within mu/RT tolerance of", tol))
+ }
}
return(w)
}
-run.guess <- function(ispecies, B=NULL, iguess=NULL) {
+run.guess <- function(ispecies, B=NULL, method="stoich", iguess=NULL) {
## run guess() using species from database 20120612
# get the stoichiometric matrix for the species
A <- i2A(ispecies)
@@ -372,7 +409,7 @@
B <- as.numeric(unlist(mB))
} else B <- as.vector(B)
# setup initial guess
- Y <- guess(A, B, iguess=iguess)
+ Y <- guess(A, B, method=method, iguess=iguess)
# take away NA guesses
Y <- Y[!is.na(Y)]
return(Y)
Modified: pkg/inst/NEWS
===================================================================
--- pkg/inst/NEWS 2012-09-18 15:06:42 UTC (rev 9)
+++ pkg/inst/NEWS 2012-09-20 01:03:23 UTC (rev 10)
@@ -1,11 +1,11 @@
-CHANGES IN CHNOSZ 0.9-7.98 (2012-09-18)
+CHANGES IN CHNOSZ 0.9-7.98 (2012-09-20)
---------------------------------------
SIGNIFICANT USER-VISIBLE CHANGES:
-- Add a Suggests dependency on 'testthat' and tests in 'inst/tests',
- particularly for functions that have been modified during this
- development cycle.
+- Create set of tests in 'inst/tests', particularly for functions that
+ have been modified during this development cycle, and add a Suggests
+ dependency on 'testthat'.
- The 'thermo' object, which holds the thermodynamic database,
and system definitions (made by the user), is now placed in an
@@ -19,6 +19,10 @@
linearly independent combinations of rows of a matrix,
is.near.equil(), and run.wjd().
+- Add guess() as another supporting function for wjd(), to produce
+ initial guesses of moles of species satisfying elemental bulk
+ composition, and a Suggests dependency on 'limSolve'.
+
MAJOR CHANGES TO FUNCTIONS:
- Refactor the primary functionality of makeup(), parsing of chemical
Modified: pkg/inst/doc/wjd.pdf
===================================================================
(Binary files differ)
Added: pkg/inst/extdata/thermo/DLEN67.csv
===================================================================
--- pkg/inst/extdata/thermo/DLEN67.csv (rev 0)
+++ pkg/inst/extdata/thermo/DLEN67.csv 2012-09-20 01:03:23 UTC (rev 10)
@@ -0,0 +1,19 @@
+compound,formula,G300,G500,G700,G1000
+water,H2O,-54.61,-52.36,-49.92,-46.04
+"carbon dioxide",CO2,-94.26,-94.39,-94.5,-94.61
+nitrogen,N2,0,0,0,0
+methane,CH4,-12.11,-7.84,-3.05,4.61
+ammonia,NH3,-3.94,1.11,6.51,14.93
+"carbon monoxide",CO,-32.85,-37.18,-41.53,-47.94
+ethane,C2H6,-7.79,1.17,10.9,26.13
+"formic acid",CH2O2,-83.85,-79.16,-74.2,-66.53
+propane,C3H8,-5.54,8.23,22.93,45.68
+"hydrogen cyanide",HCN,28.69,27.04,25.43,23.07
+methylamine,CH5N,6.65,16.59,27.27,44.66
+butane,C4H10,-3.94,14.54,34.19,64.5
+glycine,C2H5NO2,-76.76,-60.66,-43.84,-17.26
+benzene,C6H6,31.06,39.24,48.21,62.27
+hexane,C6H14,0.18,28.3,57.98,103.57
+naphthalene,C10H8,56.27,68.72,81.17,99.85
+anthracene,C14H10,82.35,100.23,118.49,146.58
+octane,C8H18,4.29,42.02,81.74,142.62
Modified: pkg/man/extdata.Rd
===================================================================
--- pkg/man/extdata.Rd 2012-09-18 15:06:42 UTC (rev 9)
+++ pkg/man/extdata.Rd 2012-09-20 01:03:23 UTC (rev 10)
@@ -74,6 +74,7 @@
\item \code{groups_big.csv} Group contribution matrix: five structural groups on the columns ([-CH3],[-CH2-],[-CH2OH],[-CO-],[-COOH]) and 24 compounds on the rows (alkanes, alcohols, ketones, acids, multiply substituted compounds).
\item \code{groups_small.csv} Group contribution matrix: twelve bond-specific groups on the columns, and 25 compounds on the rows (as above, plus isocitrate). Group identity and naming conventions adapted from Benson and Buss (1958) and Domalski and Hearing (1993). See the \sQuote{xadditivity} vignette for examples that use this file and \code{groups_big.csv}.
\item \code{RH98_Table15.csv} Group stoichiometries for high molecular weight crystalline and liquid organic compounds taken from Table 15 of Richard and Helgeson, 1998. The first three columns have the \code{compound} name, \code{formula} and physical \code{state} (\samp{cr} or \samp{liq}). The remaining columns have the numbers of each group in the compound; the names of the groups (columns) correspond to species in \code{\link{thermo}$obigt}. The compound named \samp{5a(H),14a(H)-cholestane} in the paper has been changed to \samp{5a(H),14b(H)-cholestane} here to match the group stoichiometry given in the table. See \code{\link{RH2obigt}} for a function that uses this file.
+ \item \code{DLEN67.csv} Standard Gibbs energies of formation, in kcal/mol, from Dayhoff et al., 1967, for nitrogen (N2) plus 17 compounds shown in Fig. 2 of Dayhoff et al., 1964, at 300, 500, 700 and 1000 K.
}
@@ -89,6 +90,10 @@
Benson, S. W. and Buss, J. H. (1958) Additivity rules for the estimation of molecular properties. Thermodynamic properties. \emph{J. Chem. Phys.} \bold{29}, 546--572. \url{http://dx.doi.org/10.1063/1.1744539}
+ Dayhoff, M. O. and Lippincott, E. R. and Eck, R. V. (1964) Thermodynamic Equilibria In Prebiological Atmospheres. \emph{Science} \bold{146}, 1461--1464. \url{http://dx.doi.org/10.1126/science.146.3650.1461}
+
+ Dayhoff, M. O. and Lippincott, E. R., Eck, R. V. and Nagarajan (1967) Thermodynamic Equilibrium In Prebiological Atmospheres of C, H, O, N, P, S, and Cl. Report SP-3040, National Aeronautics and Space Administration. \url{http://ntrs.nasa.gov/search.jsp?R=19670017966}
+
Domalski, E. S. and Hearing, E. D. (1993) Estimation of the thermodynamic properties of C-H-N-O-S-Halogen compounds at 298.15 K \emph{J. Phys. Chem. Ref. Data} \bold{22}, 805--1159. \url{http://dx.doi.org/10.1063/1.555927}
Gattiker, A., Michoud, K., Rivoire, C., Auchincloss, A. H., Coudert, E., Lima, T., Kersey, P., Pagni, M., Sigrist, C. J. A., Lachaize, C., Veuthey, A.-L., Gasteiger, E. and Bairoch, A. (2003) Automatic annotation of microbial proteomes in Swiss-Prot. \emph{Comput. Biol. Chem.} \bold{27}, 49--58. \url{http://dx.doi.org/10.1016/S1476-9271(02)00094-4}
Modified: pkg/man/wjd.Rd
===================================================================
--- pkg/man/wjd.Rd 2012-09-18 15:06:42 UTC (rev 9)
+++ pkg/man/wjd.Rd 2012-09-20 01:03:23 UTC (rev 10)
@@ -28,17 +28,18 @@
Gfrac = 1e-7
)
element.potentials(w, plot.it=FALSE, iplot=1:ncol(w$A))
- is.near.equil(w, tol=0.01)
+ is.near.equil(w, tol=0.01, quiet=FALSE)
guess(
A = matrix(c(
1,2,2,0,0,1,0,0,0,1,
0,0,0,1,2,1,1,0,0,0,
0,0,1,0,0,0,1,1,2,1),ncol=3,
dimnames=list(NULL,c("H","N","O"))),
- B = c(2,1,1), iguess = 1, ic = NULL
+ B = c(2,1,1), method=c("central", "stoich"), minX = 0.001, iguess = 1, ic = NULL
)
- run.wjd(ispecies, B = NULL, Y = run.guess(ispecies, B), P=1, T=25, imax = 10, Gfrac = 1e-7, tol = 0.01)
- run.guess(ispecies, B = NULL, iguess = NULL)
+ run.wjd(ispecies, B = NULL, method = "stoich", Y = run.guess(ispecies, B, method),
+ P=1, T=25, imax = 10, Gfrac = 1e-7, tol = 0.01)
+ run.guess(ispecies, B = NULL, method = "stoich", iguess = NULL)
equil.potentials(w)
}
@@ -54,7 +55,10 @@
\item{plot.it}{logical, make a plot?}
\item{iplot}{numeric, which elements for which to make plots}
\item{tol}{numeric, maximum difference in chemical potentials that counts as equilibrium}
+ \item{quiet}{logical, don't output messages?}
\item{B}{numeric, numbers of moles of the elements}
+ \item{method}{character, method used for generating an initial solution}
+ \item{minX}{numeric, minimum mole number for 'central' method}
\item{iguess}{numeric, which guess to return}
\item{ic}{numeric, which combination(s) of variable species to use (NULL for all)}
\item{ispecies}{numeric, species indices (rownumbers of \code{thermo}$obigt)}
@@ -63,18 +67,44 @@
\details{
- \code{wjd} implements the steepest descent algorithm for Gibbs energy minimization in a closed system described by White et al., 1958. Where those authors use the term \dQuote{free energy} we apply the contemporary usage of \dQuote{Gibbs energy}. The code in this function is independent of other functions or datasets in CHNOSZ. The function's default values of \code{A}, \code{G0.RT}, \code{Y} and \code{P} correspond to the example problem described by White et al., 1958 for gases in the H, N, O system at 3500 K. Note that for this, and for any other equilibrium problem that can be simulated using this function, the mole quantities in \code{Y} must all be positive numbers. Operationally, this vector defines not only the \dQuote{initial solution} but also the bulk composition of the system; it is not possible to define the bulk composition using mole numbers of elements alone. The \code{dimnames} attribute in the default value for \code{A} gives the names of the elements; this attribute is not necessary for the function to operate, but is used in the examples below to help label the plots.
+\code{wjd} implements the steepest descent algorithm for Gibbs energy minimization in a closed system described by White et al., 1958.
+\dQuote{Gibbs energy} (G) referred to here is the same as the \dQuote{free energy} (F) used by White et al., 1958.
+\code{wjd} itself is independent of other functions or datasets in CHNOSZ, but \code{run.wjd} and \code{run.guess} are provided to make access to the thermodynamic database in CHNOSZ easier.
+The default values of \code{A}, \code{G0.RT}, \code{Y} and \code{P} correspond to the example problem described by White et al., 1958 for gases in the H, N, O system at 3500 K.
+Note that for this, and for any other equilibrium problem that can be simulated using this function, the mole quantities in \code{Y} must all be positive numbers.
+Operationally, this vector defines not only the \dQuote{initial solution} but also the bulk composition of the system; it is not possible to define the bulk composition using mole numbers of elements alone.
+The \code{dimnames} attribute in the default value for \code{A} gives the names of the elements; this attribute is not necessary for the function to operate, but is used in the examples below to help label the plots.
+
White et al., 1958 describe in detail the computation of the direction of steepest descent by means of Lagrange multipliers. They propose an iterative solution to the energy minimization problem, where at each step the mole numbers of species are recomputed and a new steepest descent direction calculated from there. However, the authors only give general guidelines for computing the value of \eqn{\lambda}{lambda}, that is, the fraction of the total distance the system actually moves in the direction of steepest descent from the current point at each iteration. The two constraints given for determining the value of \eqn{\lambda}{lambda} are that all mole numbers of species are positive, and the Gibbs energy of the system actually decreases (the minimum point is not passed). In the code described here, at each iteration the minimum value of lambda, not exceeding unity, that violates the first condition is computed directly (a value of one is assigned if the mole numbers remain positive through the entire range). With the default setting of \code{nlambda}, 101 values of \eqn{\lambda}{lambda} at even intervals from 0 to this maximum permissible value are tested for the second condition at each iteration, and the highest conforming value is selected. If a value of 0 occurs, it means that the algorithm has reached an endpoint independently of the iteration and convergence tests (\code{rho} and \code{Gfrac}; see below). If this occurs, the value of \code{nlambda} might have to be increased depending on the user's needs.
The number of iterations is controlled by \code{imax} and \code{Gfrac}. The maximum number of iterations is set by \code{imax}; it can even be zero, though such a setting might only be useful in combination with \code{element.potentials} to characterize the initial state of a minimization problem. Within the limit of \code{imax}, the iterations continue until the magnitutde of the change in total Gibbs energy of the system, as a fraction of the system's energy in the previous iteration, is lower than the value specified in \code{Gfrac}. For the first example below, the default setting of \code{Gfrac} causes the algorithm to stop after six iterations.
Using the output of \code{wjd}, provided in the argument \code{w}, \code{element.potentials} calculates the chemical potentials of the elements in the system. It does so by combining the values of \code{G0.RT} of species with the inverses of stoichiometric matrices of combinations of species whose elemental compositions are linearly independent from each other. These possible combinations are constructed using the function \code{\link{invertible.combs}}. The value returned by \code{element.potentials} is a matrix, with each column corresponding to a different element and each row to a different combination of species. The entries in the matrix are the chemical potentials of the elements divided by \eqn{RT}{RT}. If \code{plot.it} is set to TRUE, the chemical potentials of the elements are plotted as a function of species combination number, with as many plots as elements, unless \code{iplot} is set to another value (e.g. \samp{c(1,3)} for only elements 1 and 3). In the first example below, the number of unique combinations of species is 120, but only 76 of these combinations provide stoichiometric independence.
- There is no guarantee that the algorithm will converge on a global (or even be close to a local) minimum. However, some tests are available to help assess the likelihood that a solution is close to equilibrium. A necessary condition of equilibrium is that the chemical potentials of the elements be independent of the particular combination of species used to compute them. \code{is.near.equil} compares the chemical poetntials for each element, computed using \code{element.potentials}, with the value of \code{tol}. If, for each element, the range of potentials/RT (difference between minimum and maximum) is less than \code{tol}, the result is TRUE, otherwise the function returns FALSE. The default value of \code{tol} corresponds to an energy of 0.01 * 1.9872 * 298.15 = ca. 6 cal/mol at 25 degrees C.
+There is no guarantee that the algorithm will converge on a global (or even be close to a local) minimum.
+However, some tests are available to help assess the likelihood that a solution is close to equilibrium.
+A necessary condition of equilibrium is that the chemical potentials of the elements be independent of the particular combination of species used to compute them.
+\code{is.near.equil} compares the chemical potentials for each element, computed using \code{element.potentials}, with the value of \code{tol}.
+If, for each element, the range of potentials/RT (difference between minimum and maximum) is less than \code{tol}, the result is TRUE, otherwise the function returns FALSE, and prints a message unless \code{quiet} is TRUE.
+The default value of \code{tol} corresponds to an energy of 0.01 * 1.9872 * 298.15 = ca. 6 cal/mol at 25 degrees C.
- One of the limitations of the algorithm coded in \code{wjd} is that the initial solution, and all iterations, require positive (non-zero) numbers of moles of every species. Often, when investigating an equilibrium problem, the stoichiometric constraints are expressed most readily in terms of bulk composition -- numbers of moles of each element. \code{guess} is a function to make initial guesses about the numbers of moles of all species in the system. Its system-specific arguments are the stoichiometric matrix \code{A} (as defined above for \code{wjd}) and the bulk composition vector \code{B}, giving the number of moles of each element, in the same order that the elements appear in \code{A}. The method implemented in this function is to test successive combinations of species whose compositions are linearly independent (all those from \code{\link{invertible.combs}} if \code{ic} is NULL, or only those identified by \code{ic}). Each combination is referred to as the \sQuote{variable} species. Then the function sets the moles of all \sQuote{other} species to a single value. This value is determined by the constraint that the greatest proportion, relative to the bulk composition in \code{B}, of any element contributed by all the \sQuote{other} species is equal to \code{max.frac} (see code). The function tests nine hard-coded values of \code{max.frac} from 0.01 to 0.99 until an accepted result is found: The function solves for the moles of the \sQuote{variable} species that make up the difference in numbers of moles of elements. If the numbers of moles of all the \sQuote{variable} species is positive, this guess is accepted. The first such guess is returned if \code{iguess} is 1 (the default), and higher values of this variable indicate which successive guess to return. If \code{iguess} is NULL, all results are returned in a list, with non-successful guesses indicated by NA. In the first example below, of the 76 combinations of species whose elemental compositions are linearly independent, 32 yield guesses following this method.
+One of the constraints of the algorithm coded in \code{wjd} is that the initial solution, and all iterations, require positive (non-zero) numbers of moles of every species.
+Often, when investigating an equilibrium problem, the stoichiometric constraints are expressed most readily in terms of bulk composition -- numbers of moles of each element.
+\code{guess} is a function to make initial guesses about the numbers of moles of all species in the system subject to the positivity constraints.
+Its system-specific arguments are the stoichiometric matrix \code{A} (as defined above for \code{wjd}) and the bulk composition vector \code{B}, giving the number of moles of each element, in the same order that the elements appear in \code{A}.
+The first \code{method} available in \code{guess} generates the \sQuote{central} solution of the system of linear equations using the \code{\link[limSolve]{xranges}} function from \pkg{limSolve}. The central solution is the mean of ranges of unknowns. The inequality constraint, or minimum number of moles of any species, is given by \code{minX}.
+The second method for \code{guess} \sQuote{stoich} (and the default for \code{run.guess} and \code{run.wjd}) is to test successive combinations of species whose elemental compositions are linearly independent.
+The linearly independent combinations tested are all those from \code{\link{invertible.combs}} if \code{ic} is NULL, or only those identified by \code{ic}.
+Each combination is referred to as the \sQuote{variable} species; the moles of all \sQuote{other} species are set to a single value.
+This value is determined by the constraint that the greatest proportion, relative to the bulk composition in \code{B}, of any element contributed by all the \sQuote{other} species is equal to a value in \code{max.frac} (see code).
+The function tests nine hard-coded values of \code{max.frac} from 0.01 to 0.99, at each one solving for the moles of the \sQuote{variable} species that make up the difference in numbers of moles of elements.
+If the numbers of moles of all the \sQuote{variable} species is positive, the guess is accepted.
+The first accepted guess is returned if \code{iguess} is 1 (the default); other values of \code{iguess} indicate which guess to return.
+If \code{iguess} is NULL, all results are returned in a list, with non-successful guesses indicated by NA.
+In the first example below, of the 76 combinations of species whose elemental compositions are linearly independent, 32 yield guesses following this method.
+
\code{run.wjd} is a wrapper function to run \code{wjd}, provided the \code{ispecies} in the thermodynamic database (\code{\link{thermo}$obigt}), the chemical composition of the system in \code{B}, and the pressure \code{P} and temperature \code{T}; the latter are passed to \code{\link{subcrt}} (with \code{exceed.Ttr = TRUE}) to generate the values of \code{G0.RT} for \code{wjd}. Alternatively to \code{B}, the initial guess of numbers of moles of species can be provided in \code{Y}; otherwise as many combinations of \code{Y} as returned from \code{run.guess} are tested until one is found that \code{is.near.equil}. The function gives an error in none of the guesses in \code{Y} is near equilibrium, within the tolerance set by \code{tol}.
\code{run.guess} is a wrapper function to call \code{guess} using the stoichiometric matrix \code{A} built from the \code{ispecies} indices in the thermodynamic database.
@@ -118,10 +148,10 @@
is.near.equil(w) # FALSE
## test all of the guesses of inititial mole quantities
-## provided by guess() using default setting (H2NO)
-# 2 of them give FALSE for is.near.equil the tolerance lowered to 0.001
+## provided by guess() using default bulk composition (H2NO)
+# 9 of them are not is.near.equil with the tolerance lowered to 0.0001
sapply( 1:32, function(i)
- is.near.equil(wjd(Y=guess(iguess=i)), tol=0.001) )
+ is.near.equil(wjd(Y=guess(method = "stoich", iguess=i)), tol=0.0001) )
\dontshow{par(opar)}
## using run.wjd(): 20 crystalline amino acids
Modified: pkg/vignettes/vig.bib
===================================================================
--- pkg/vignettes/vig.bib 2012-09-18 15:06:42 UTC (rev 9)
+++ pkg/vignettes/vig.bib 2012-09-20 01:03:23 UTC (rev 10)
@@ -1,4 +1,4 @@
-% This file was created with JabRef 2.8.
+% This file was created with JabRef 2.8.1.
% Encoding: ISO8859_1
@ARTICLE{ABR+06,
@@ -69,7 +69,7 @@
author = {Anderson, N. L. and Anderson, N. G.},
title = {{T}he human plasma proteome: {H}istory, character, and diagnostic
prospects (vol 1, pg 845, 2002)},
- journal = {Molecular \& Cellular Proteomics},
+ journal = {Mol. Cell. Proteomics},
year = {2003},
volume = {2},
pages = {50--50},
@@ -83,7 +83,7 @@
author = {Anderson, N. L. and Anderson, N. G.},
title = {{T}he human plasma proteome - {H}istory, character, and diagnostic
prospects},
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/chnosz -r 10
More information about the CHNOSZ-commits
mailing list