[CHNOSZ-commits] r38 - in pkg/CHNOSZ: . R inst inst/doc inst/tests man vignettes
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Jan 10 02:21:55 CET 2013
Author: jedick
Date: 2013-01-10 02:21:55 +0100 (Thu, 10 Jan 2013)
New Revision: 38
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/R/wjd.R
pkg/CHNOSZ/inst/NEWS
pkg/CHNOSZ/inst/doc/wjd.pdf
pkg/CHNOSZ/inst/tests/test-wjd.R
pkg/CHNOSZ/man/wjd.Rd
pkg/CHNOSZ/vignettes/wjd.Rnw
pkg/CHNOSZ/vignettes/wjd.lyx
Log:
fix calculation of free energy derivative in wjd()
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2013-01-03 04:48:59 UTC (rev 37)
+++ pkg/CHNOSZ/DESCRIPTION 2013-01-10 01:21:55 UTC (rev 38)
@@ -1,6 +1,6 @@
-Date: 2013-01-01
+Date: 2013-01-10
Package: CHNOSZ
-Version: 0.9-9
+Version: 0.9-9.1
Title: Chemical Thermodynamics and Activity Diagrams
Author: Jeffrey M. Dick
Maintainer: Jeffrey M. Dick <jmdick at asu.edu>
Modified: pkg/CHNOSZ/R/wjd.R
===================================================================
--- pkg/CHNOSZ/R/wjd.R 2013-01-03 04:48:59 UTC (rev 37)
+++ pkg/CHNOSZ/R/wjd.R 2013-01-10 01:21:55 UTC (rev 38)
@@ -163,7 +163,7 @@
# second constraint is that the derivative of free energy
# doesn't go positive ... Eq. 20
d.f.Y <- function(lambda,Y,D,C) {
- d.f.Y <- sum(D * C + log( (Y + lambda * D) / (sum(Y) + lambda * sum(D)) ))
+ d.f.Y <- sum(D * (C + log( (Y + lambda * D) / (sum(Y) + lambda * sum(D)) )))
return(d.f.Y)
}
# what are the free energy derivatives
Modified: pkg/CHNOSZ/inst/NEWS
===================================================================
--- pkg/CHNOSZ/inst/NEWS 2013-01-03 04:48:59 UTC (rev 37)
+++ pkg/CHNOSZ/inst/NEWS 2013-01-10 01:21:55 UTC (rev 38)
@@ -1,3 +1,8 @@
+CHANGES IN CHNOSZ 0.9-9.1 (2013-01-10)
+--------------------------------------
+
+- Fix calculation of free energy derivative in wjd().
+
CHANGES IN CHNOSZ 0.9-9 (2013-01-01)
------------------------------------
Modified: pkg/CHNOSZ/inst/doc/wjd.pdf
===================================================================
(Binary files differ)
Modified: pkg/CHNOSZ/inst/tests/test-wjd.R
===================================================================
--- pkg/CHNOSZ/inst/tests/test-wjd.R 2013-01-03 04:48:59 UTC (rev 37)
+++ pkg/CHNOSZ/inst/tests/test-wjd.R 2013-01-10 01:21:55 UTC (rev 38)
@@ -7,7 +7,7 @@
expect_equal(X, w$X, tol=1e-4)
})
-test_that("guess() operates on intermediate compositionsi but fails with endmembers", {
+test_that("guess() operates on intermediate compositions but fails with endmembers", {
alkanes <- c("n-hexane", "n-heptane", "n-octane", "n-nonane")
ialk <- info(alkanes, "liq")
expect_error(run.guess(ialk, "C6H14"), "there are only 0")
@@ -18,6 +18,48 @@
expect_error(run.guess(ialk, "C9H20"), "there are only 0")
})
+test_that("open-system equilibrium distributions reproduce the results of wjd()", {
+ ### set up system
+ # use proteins in the lipid particle (n=19)
+ y <- yeastgfp("lipid.particle")
+ # get the amino acid compositions of the proteins
+ aa <- more.aa(y$protein, "Sce")
+ # don't use those with NA abundance or sequence (leaves n=17)
+ ina <- is.na(y$abundance) | is.na(aa$chains)
+ aa <- aa[!ina, ]
+ # normalize the proteins to single residues; columns 6:25 are the amino acid counts
+ aa.625 <- aa[, 6:25]
+ aa[, 6:25] <- aa.625 / rowSums(aa.625)
+ # add proteins to thermo$protein
+ iprotein <- add.protein(aa)
+ # add proteins to thermo$obigt
+ iobigt <- info(paste(aa$protein, aa$organism, sep="_"))
+ ### closed system calculation (constant composition of elements)
+ # use equal initial abundances
+ Y <- rep(100, length(iobigt))
+ # run the Gibbs energy minimization (this did not iterate before 20130109,
+ # due to bug in calculation of free energy derivative)
+ w <- run.wjd(iobigt, Y=Y, Gfrac=1e-15, nlambda=1001)
+ # the molar abundances
+ X.closed <- w$X
+ # get the chemical potentials of the elements
+ ep <- equil.potentials(w)
+ # the corresponding logarithms of activities of the basis species
+ basis("CHNOS")
+ bl <- basis.logact(ep)
+ ### open system calculation (constant chemical potentials of basis species)
+ # set the activities of the basis species
+ basis(names(bl), bl)
+ # get the affinities of the formation reactions
+ a <- affinity(iprotein=iprotein)
+ # then the equilibrium abundances, with total moles of residues as used in wjd
+ e <- equilibrate(a, loga.balance=log10(sum(Y)))
+ X.open <- 10^unlist(e$loga.equil)
+ # the test: abundances calculated both ways are equal
+ expect_equal(X.closed, X.open, tol=0.019)
+ # seems that we could do better than that 1.9% mean difference!
+})
+
# see also test-swap.basis.R for an example using run.wjd() and
# equil.potentials() to generate chemical potentials of elements
Modified: pkg/CHNOSZ/man/wjd.Rd
===================================================================
--- pkg/CHNOSZ/man/wjd.Rd 2013-01-03 04:48:59 UTC (rev 37)
+++ pkg/CHNOSZ/man/wjd.Rd 2013-01-10 01:21:55 UTC (rev 38)
@@ -106,7 +106,9 @@
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.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 if 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.
Modified: pkg/CHNOSZ/vignettes/wjd.Rnw
===================================================================
--- pkg/CHNOSZ/vignettes/wjd.Rnw 2013-01-03 04:48:59 UTC (rev 37)
+++ pkg/CHNOSZ/vignettes/wjd.Rnw 2013-01-10 01:21:55 UTC (rev 38)
@@ -105,20 +105,21 @@
# dropping elements with zero stoichiometry
f1 <- as.chemical.formula(w$A[oX[1],w$A[oX[1],]!=0])
f2 <- as.chemical.formula(w$A[oX[2],w$A[oX[2],]!=0])
+f4 <- as.chemical.formula(w$A[oX[4],w$A[oX[4],]!=0])
@
We find that \Sexpr{f1} and \Sexpr{f2} are the most abundant species,
-after \Sexpr{niter} iterations.
+after \Sexpr{niter} iterations. The fourth most-abundant species,
+\Sexpr{f4}, is discussed below.
Let's track their mole fractions through the iterations. Write a function
-that returns the mole fractions of these two species, after a specified
-number of iterations.
+that returns the mole fractions of a species, after a specified number
+of iterations.
<<iterfun>>=
-iterfun <- function(imax) {
+iterfun <- function(imax, i) {
w <- wjd(imax=imax)
- x1 <- w$X[oX[1]]/sum(w$X)
- x2 <- w$X[oX[2]]/sum(w$X)
- return(list(x1=x1, x2=x2))
+ x <- w$X[i]/sum(w$X)
+ return(x)
}
@
Then apply the function over the different numbers of iterations,
@@ -126,28 +127,35 @@
\setkeys{Gin}{width=1.0\textwidth}
<<iterplot, fig=TRUE, height=3>>=
-sa <- sapply(0:niter, iterfun)
par(mfrow=c(1, 2))
-plot(0:niter, sa[1, ], xlab="iteration", ylab=paste("x", f1))
-plot(0:niter, sa[2, ], xlab="iteration", ylab=paste("x", f2))
+sa <- sapply(0:niter, iterfun, i=oX[1])
+plot(0:niter, sa, xlab="iteration", ylab=paste("x", f1))
+sa <- sapply(0:niter, iterfun, i=oX[4])
+plot(0:niter, sa, xlab="iteration", ylab=paste("x", f4))
@
\setkeys{Gin}{width=1.0\textwidth}
-A bit of winding: the mole fractions of \Sexpr{f1} and \Sexpr{f2}
-increase up to the third iteration, but at the fourth (and less so
-for the fifth and sixth), their mole fractions decrease. This behavior
-is consistent with a decrease in Gibbs energy at \emph{every} iteration;
-that can be verified by inspecting the values in the \Sexpr{niter}-iteration
+A bit of winding: the mole fraction of \Sexpr{f1} generally increases
+(there is a very slight decrease at the fifth iteration), but the
+mole fraction of \Sexpr{f4} increases in the first two iterations,
+then decreases in following ones. This behavior is consistent with
+a decrease in Gibbs energy at \emph{every} iteration; that can be
+verified by inspecting the values in the \Sexpr{niter}-iteration
result:
<<Gdown>>=
-diff(w$F.Y)
+all(diff(w$F.Y) < 0)
@
The decrease in Gibbs energy becomes smaller with every iteration,
-as an equilibrium distribution of species is approached.
+as an equilibrium distribution of species is approached:
+<<Gdownslower>>=
+identical(diff(w$F.Y), sort(diff(w$F.Y)))
+@
+
+
\section{Is it near the bottom? Testing for equilibrium}
@@ -215,13 +223,13 @@
species combinations. It returns TRUE if the maximum difference for
any element is below $\mu/RT=0.01$, but this limit can be changed
through the \texttt{tol} argument to the function. With the default
-setting, \texttt{is.near.equil()} does pass only 6-iteration result,
-but by reducing the tolerance, the tests can be made to fail:
+setting, \texttt{is.near.equil()} does pass the \Sexpr{niter}-iteration
+result, but by reducing the tolerance, the tests can be made to fail:
<<w_ep_plot>>=
-is.near.equil(w3)
-is.near.equil(w)
-is.near.equil(w, tol=0.0001)
+is.near.equil(w3) # 3 iterations
+is.near.equil(w) # 7 iterations
+is.near.equil(w, tol=0.00001)
@
Modified: pkg/CHNOSZ/vignettes/wjd.lyx
===================================================================
--- pkg/CHNOSZ/vignettes/wjd.lyx 2013-01-03 04:48:59 UTC (rev 37)
+++ pkg/CHNOSZ/vignettes/wjd.lyx 2013-01-10 01:21:55 UTC (rev 38)
@@ -397,6 +397,11 @@
\begin_layout Chunk
+f4 <- as.chemical.formula(w$A[oX[4],w$A[oX[4],]!=0])
+\end_layout
+
+\begin_layout Chunk
+
@
\end_layout
@@ -442,13 +447,26 @@
\end_inset
iterations.
-
+ The fourth most-abundant species,
+\begin_inset ERT
+status collapsed
+
+\begin_layout Plain Layout
+
+
+\backslash
+Sexpr{f4}
\end_layout
+\end_inset
+
+, is discussed below.
+\end_layout
+
\begin_layout Standard
Let's track their mole fractions through the iterations.
- Write a function that returns the mole fractions of these two species,
- after a specified number of iterations.
+ Write a function that returns the mole fractions of a species, after a
+ specified number of iterations.
\end_layout
\begin_layout Standard
@@ -462,7 +480,7 @@
\begin_layout Chunk
-iterfun <- function(imax) {
+iterfun <- function(imax, i) {
\end_layout
\begin_layout Chunk
@@ -472,21 +490,16 @@
\begin_layout Chunk
- x1 <- w$X[oX[1]]/sum(w$X)
+ x <- w$X[i]/sum(w$X)
\end_layout
\begin_layout Chunk
- x2 <- w$X[oX[2]]/sum(w$X)
+ return(x)
\end_layout
\begin_layout Chunk
- return(list(x1=x1, x2=x2))
-\end_layout
-
-\begin_layout Chunk
-
}
\end_layout
@@ -534,26 +547,31 @@
\begin_layout Chunk
-sa <- sapply(0:niter, iterfun)
+par(mfrow=c(1, 2))
\end_layout
\begin_layout Chunk
-par(mfrow=c(1, 2))
+sa <- sapply(0:niter, iterfun, i=oX[1])
\end_layout
\begin_layout Chunk
-plot(0:niter, sa[1, ], xlab="iteration", ylab=paste("x", f1))
+plot(0:niter, sa, xlab="iteration", ylab=paste("x", f1))
\end_layout
\begin_layout Chunk
-plot(0:niter, sa[2, ], xlab="iteration", ylab=paste("x", f2))
+sa <- sapply(0:niter, iterfun, i=oX[4])
\end_layout
\begin_layout Chunk
+plot(0:niter, sa, xlab="iteration", ylab=paste("x", f4))
+\end_layout
+
+\begin_layout Chunk
+
@
\end_layout
@@ -572,7 +590,7 @@
\end_layout
\begin_layout Standard
-A bit of winding: the mole fractions of
+A bit of winding: the mole fraction of
\begin_inset ERT
status collapsed
@@ -585,21 +603,21 @@
\end_inset
- and
+ generally increases (there is a very slight decrease at the fifth iteration),
+ but the mole fraction of
\begin_inset ERT
-status collapsed
+status open
\begin_layout Plain Layout
\backslash
-Sexpr{f2}
+Sexpr{f4}
\end_layout
\end_inset
- increase up to the third iteration, but at the fourth (and less so for
- the fifth and sixth), their mole fractions decrease.
+ increases in the first two iterations, then decreases in following ones.
This behavior is consistent with a decrease in Gibbs energy at
\emph on
every
@@ -631,7 +649,7 @@
\begin_layout Chunk
-diff(w$F.Y)
+all(diff(w$F.Y) < 0)
\end_layout
\begin_layout Chunk
@@ -642,14 +660,38 @@
\end_inset
The decrease in Gibbs energy becomes smaller with every iteration, as an
- equilibrium distribution of species is approached.
+ equilibrium distribution of species is approached:
\end_layout
+\begin_layout Standard
+\begin_inset Branch short
+status open
+
+\begin_layout Chunk
+
+<<Gdownslower>>=
+\end_layout
+
+\begin_layout Chunk
+
+identical(diff(w$F.Y), sort(diff(w$F.Y)))
+\end_layout
+
+\begin_layout Chunk
+
+@
+\end_layout
+
\end_inset
\end_layout
+\end_inset
+
+
+\end_layout
+
\begin_layout Section
Is it near the bottom? Testing for equilibrium
\end_layout
@@ -872,10 +914,23 @@
\family typewriter
is.near.equil()
\family default
- does pass only 6-iteration result, but by reducing the tolerance, the tests
- can be made to fail:
+ does pass the
+\begin_inset ERT
+status collapsed
+
+\begin_layout Plain Layout
+
+
+\backslash
+Sexpr{niter}
\end_layout
+\end_inset
+
+-iteration result, but by reducing the tolerance, the tests can be made
+ to fail:
+\end_layout
+
\begin_layout Standard
\begin_inset Branch short
status open
@@ -887,17 +942,17 @@
\begin_layout Chunk
-is.near.equil(w3)
+is.near.equil(w3) # 3 iterations
\end_layout
\begin_layout Chunk
-is.near.equil(w)
+is.near.equil(w) # 7 iterations
\end_layout
\begin_layout Chunk
-is.near.equil(w, tol=0.0001)
+is.near.equil(w, tol=0.00001)
\end_layout
\begin_layout Chunk
More information about the CHNOSZ-commits
mailing list