[CHNOSZ-commits] r9 - in pkg: . R inst inst/tests man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Sep 18 17:06:43 CEST 2012
Author: jedick
Date: 2012-09-18 17:06:42 +0200 (Tue, 18 Sep 2012)
New Revision: 9
Modified:
pkg/DESCRIPTION
pkg/R/diagram.R
pkg/R/revisit.R
pkg/R/util.misc.R
pkg/inst/NEWS
pkg/inst/tests/test-revisit.R
pkg/man/revisit.Rd
pkg/man/util.misc.Rd
Log:
add DGT() function and target in revisit()
fix drawing of 'dotted' lines in diagram()
Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION 2012-09-04 23:20:56 UTC (rev 8)
+++ pkg/DESCRIPTION 2012-09-18 15:06:42 UTC (rev 9)
@@ -1,4 +1,4 @@
-Date: 2012-09-05
+Date: 2012-09-18
Package: CHNOSZ
Version: 0.9-7.98
Title: Chemical Thermodynamics and Activity Diagrams
Modified: pkg/R/diagram.R
===================================================================
--- pkg/R/diagram.R 2012-09-04 23:20:56 UTC (rev 8)
+++ pkg/R/diagram.R 2012-09-18 15:06:42 UTC (rev 9)
@@ -24,10 +24,9 @@
if(is.null(ispecies)) {
ispecies <- 1:nspecies
# take out species that have NA affinities
- if(nd > 0) {
- for(i in 1:nspecies) if(all(is.na(aval[[i]]))) ispecies[i] <- NA
- ispecies <- ispecies[!is.na(ispecies)]
- }
+ for(i in 1:nspecies) if(all(is.na(aval[[i]]))) ispecies[i] <- NA
+ ispecies <- ispecies[!is.na(ispecies)]
+ if(length(ispecies)==0) stop("all species have NA affinities")
}
if(!identical(ispecies,1:nspecies)) {
if(length(ispecies)==0) {
@@ -239,26 +238,30 @@
vline <- function(out,ix) {
ny <- nrow(out)
xs <- rep(ix,ny*2+1)
- if(0 %in% (ix%%dotted))
- return(list(xs=xs,ys=rep(NA,length(xs))))
ys <- c(rep(ny:1,each=2),0)
y1 <- out[,ix]
y2 <- out[,ix+1]
+ # no line segment inside a stability field
iy <- which(y1==y2)
ys[iy*2] <- NA
+ # no line segment at a dotted position
+ iyd <- which(ys%%dotted==0)
+ ys[iyd] <- NA
return(list(xs=xs,ys=ys))
}
hline <- function(out,iy) {
nx <- nrow(out)
ny <- ncol(out)
ys <- rep(ny-iy,nx*2+1)
- if(0 %in% (iy%%dotted))
- return(list(xs=rep(NA,length(ys)),ys=ys))
xs <- c(0,rep(1:nx,each=2))
x1 <- out[iy,]
x2 <- out[iy+1,]
+ # no line segment inside a stability field
ix <- which(x1==x2)
xs[ix*2] <- NA
+ # no line segment at a dotted position
+ ixd <- which(xs%%dotted==0)
+ xs[ixd] <- NA
return(list(xs=xs,ys=ys))
}
clipfun <- function(z,zlim) {
@@ -441,10 +444,12 @@
else title(main=as.expression(axis.title(what,balance.title(balance))))
}
}
- if(what=="logact" | do.loga.equil) return(invisible(list(basis=affinity$basis,species=affinity$species,
+ out <- list(basis=affinity$basis,species=affinity$species,
T=affinity$T,P=affinity$P,xname=affinity$xname,xlim=affinity$xlim,yname=affinity$yname,
- ylim=affinity$ylim,logact=A)))
- else return(invisible(list(aval)))
+ ylim=affinity$ylim,aval=aval)
+ if(what=="logact" | do.loga.equil) out <- c(out, list(logact=A))
+ if(what=="logact") out <- c(out, list(Astar=Astar))
+ return(invisible(out))
}
### 1-D (property or speciation) diagram
if(nd==1) {
@@ -525,11 +530,11 @@
}
if(alpha) for(i in 1:length(A)) A[[i]] <- 10^A[[i]]
# 20090324 return list with single element 'logact'
- out <- list(basis=affinity$basis,species=affinity$species,
+ out <- list(basis=affinity$basis,species=affinity$species,aval=aval,
T=affinity$T,P=affinity$P,xname=affinity$xname,xlim=affinity$xlim,yname=affinity$yname,
- ylim=affinity$ylim)
+ ylim=affinity$ylim,aval=aval)
if(what=="logact" | do.loga.equil) out <- c(out,list(logact=A))
- else out <- c(out,list(aval=aval))
+ if(what=="logact") out <- c(out, list(Astar=Astar))
return(invisible(out))
}
@@ -596,13 +601,13 @@
if(mam) return(invisible(list(basis=affinity$basis,species=species,T=affinity$T,P=affinity$P,
xname=affinity$xname,xlim=affinity$xlim,yname=affinity$yname,ylim=affinity$ylim,out=out,aval=aval)))
else return(invisible(list(basis=affinity$basis,species=species,T=affinity$T,P=affinity$P,xname=affinity$xname,
- xlim=affinity$xlim,yname=affinity$yname,ylim=affinity$ylim,out=out,aval=aval,logact=A)))
+ xlim=affinity$xlim,yname=affinity$yname,ylim=affinity$ylim,out=out,aval=aval,Astar=Astar,logact=A)))
} else {
#msgout(paste('diagram: 2-D plot of',property,'not available\n'))
if(mam) return(invisible(list(basis=affinity$basis,species=species,T=affinity$T,P=affinity$P,
xname=affinity$xname,xlim=affinity$xlim,yname=affinity$yname,ylim=affinity$ylim,aval=aval)))
else return(invisible(list(basis=affinity$basis,species=species,T=affinity$T,P=affinity$P,xname=affinity$xname,
- xlim=affinity$xlim,yname=affinity$yname,ylim=affinity$ylim,aval=aval,logact=A)))
+ xlim=affinity$xlim,yname=affinity$yname,ylim=affinity$ylim,aval=aval,Astar=Astar,logact=A)))
}
} else {
# if we're not calculating predominances we can only make a contour plot for properties
Modified: pkg/R/revisit.R
===================================================================
--- pkg/R/revisit.R 2012-09-04 23:20:56 UTC (rev 8)
+++ pkg/R/revisit.R 2012-09-18 15:06:42 UTC (rev 9)
@@ -5,7 +5,7 @@
where.extreme <- function(z, target, do.sat=FALSE) {
if(missing(target)) stop("no target specified")
# are we interested in a maximum or minimum?
- if(tolower(target) %in% c("sd", "sd.log", "cv", "cv.log", "rmsd", "cvrmsd"))
+ if(tolower(target) %in% c("sd", "sd.log", "cv", "cv.log", "rmsd", "cvrmsd", "dgt"))
myext <- "minimum" else myext <- "maximum"
# do we care about the sign of the index?
if(tolower(target) %in% c("sd", "sd.log", "cv", "cv.log", "rmsd", "cvrmsd"))
@@ -40,7 +40,7 @@
extremes <- function(z, target) {
if(missing(target)) stop("no target specified")
# are we interested in a maximum or minimum?
- if(tolower(target) %in% c("sd", "sd.log", "cv", "cv.log", "rmsd", "cvrmsd"))
+ if(tolower(target) %in% c("sd", "sd.log", "cv", "cv.log", "rmsd", "cvrmsd", "dgt"))
myext <- "minimum" else myext <- "maximum"
# do we care about the sign of the index?
if(tolower(target) %in% c("sd", "sd.log", "cv", "cv.log", "rmsd", "cvrmsd"))
@@ -64,7 +64,7 @@
revisit <- function(d, target="cv", loga.ref=NULL,
plot.it=NULL, col=par("fg"), yline=2, ylim=NULL, ispecies=NULL, add=FALSE,
cex=par("cex"), lwd=par("lwd"), mar=NULL, side=1:4, xlim=NULL, labcex=0.6,
- pch=1, legend="", legend.x=NULL, lpch=NULL, main=NULL, lograt.ref=NULL, plot.ext=TRUE) {
+ pch=1, legend="", legend.x=NULL, lpch=NULL, main=NULL, lograt.ref=NULL, plot.ext=TRUE, DGT.swap12=FALSE) {
# calculate and plot diversity indices of relative abundances
# 20090316 jmd
# d can be the output from diagram (enables plotting)
@@ -87,8 +87,9 @@
logact <- d$logact
}
# check that all needed arguments are present
- if(target %in% c("richness", "cvrmsd", "spearman", "pearson", "logact") & is.null(loga.ref))
- stop(paste("for '",target,"' target, loga.ref must be supplied", sep=""))
+ target.lower <- tolower(target)
+ if(target.lower %in% c("richness", "cvrmsd", "spearman", "pearson", "logact", "dgt") & is.null(loga.ref))
+ stop(paste("for '", target, "' target, loga.ref must be supplied", sep=""))
# take a subset (or all) of the species
if(is.null(ispecies)) ispecies <- 1:length(logact)
logact <- logact[ispecies]
@@ -104,18 +105,19 @@
# given a list of logarithms of activities of species
# (as vectors or matrices or higher dimensional arrays)
# calculate a diversity index of the same dimensions
- # available targets:
+ # these targets only depend on the logarithms of activities from diagram() (logact):
# "shannon" shannon entropy
# "sd" standard deviation
# "cv" coefficient of variation
# "sd.log", "cv.log" SD/CV for the logarithms of activity
# "qqr" correlation coefficient on q-q plot (i.e., normality test)
-
- # "richness" species richness (loga.ref)
- # "cvrmsd" coefficient of variation of rmsd (loga.ref)
- # "spearman" spearman correlation coefficient (loga.ref)
- # "pearson" pearson correlation coefficient (log.ref)
- # "logact" maximize the activity of a species (loga.ref)
+ # these targets also depend on reference logarithms of activities (loga.ref):
+ # "richness" species richness
+ # "cvrmsd" coefficient of variation of rmsd
+ # "spearman" spearman correlation coefficient
+ # "pearson" pearson correlation coefficient
+ # "logact" maximize the activity of a species
+ # "DGT" minimize the Gibbs energy of transformation
# vectorize the entries in the logact list
for(i in 1:ns) {
@@ -128,7 +130,7 @@
H[] <- 0
# ratio-specific calculations
if(!is.null(lograt.ref)) {
- if(!target %in% c("rmsd", "cvrmsd", "spearman", "pearson"))
+ if(!target.lower %in% c("rmsd", "cvrmsd", "spearman", "pearson"))
stop(paste("target",target,"not available when comparing activity ratios"))
else(msgout("revisit: calculating activity ratios\n"))
# instead of logact we use calculated log activity ratio
@@ -136,7 +138,6 @@
loga.ref <- lograt.ref
}
# target-specific calculations
- target.lower <- tolower(target)
if(target.lower %in% c("sd", "cv")) {
# build a matrix; rows are species
@@ -239,6 +240,22 @@
} else if(target.lower=="logact") {
# where the activity of a species is maximal
H <- logact[[loga.ref]]
+
+ } else if(target.lower=="dgt") {
+ # Gibbs energy of transformation to the observed assemblage
+ actarr <- list2array(logact)
+ Astararr <- list2array(d$Astar)
+ Gfun <- function(i, actarr, Astararr) {
+ loga.equil <- actarr[i, ]
+ Astar <- Astararr[i, ]
+ # direction of transformation:
+ # swap12=FALSE: loga.equil(Astar) --> loga.ref
+ # swap12=TRUE: loga.ref(Astar) --> loga.equil
+ if(DGT.swap12) out <- -DGT(loga.ref, loga.equil, Astar)
+ else out <- DGT(loga.equil, loga.ref, Astar)
+ return(out)
+ }
+ H <- as.numeric(palply(1:length(H), Gfun, actarr, Astararr))
} else stop(paste("specified target '", target, "' not available", sep=""))
# replace dims
Modified: pkg/R/util.misc.R
===================================================================
--- pkg/R/util.misc.R 2012-09-04 23:20:56 UTC (rev 8)
+++ pkg/R/util.misc.R 2012-09-18 15:06:42 UTC (rev 9)
@@ -141,3 +141,22 @@
# done!
}
+DGT <- function(loga1, loga2, Astar) {
+ # calculate the Gibbs energy/2.303RT of transformation
+ # from an initial to a final assemblage at constant T, P and
+ # chemical potentials of basis species 20120917
+ # loga1 - logarithms of activity of species in the initial assemblage
+ # loga2 - logarithms of activity of species in the final assemblage
+ # Astar - starved (of activity of the species of interest) values of chemical affinity/2.303RT
+ # first remove the logarithms
+ a1 <- 10^loga1
+ a2 <- 10^loga2
+ # the molal Gibbs energy in the initial and final states
+ # (derived from integrating A = Astar - RTln(a) from a1 to a2)
+ G1 <- -a1*Astar + a1*loga1 - a1/log(10)
+ G2 <- -a2*Astar + a2*loga2 - a2/log(10)
+ # calculate the change in molal Gibbs energy for each species
+ GT <- G2 - G1
+ # return the sum
+ return(sum(GT))
+}
Modified: pkg/inst/NEWS
===================================================================
--- pkg/inst/NEWS 2012-09-04 23:20:56 UTC (rev 8)
+++ pkg/inst/NEWS 2012-09-18 15:06:42 UTC (rev 9)
@@ -1,4 +1,4 @@
-CHANGES IN CHNOSZ 0.9-7.98 (2012-09-05)
+CHANGES IN CHNOSZ 0.9-7.98 (2012-09-18)
---------------------------------------
SIGNIFICANT USER-VISIBLE CHANGES:
@@ -148,6 +148,10 @@
- Add protein.equil() for step-by-step calculation of chemical
activities of proteins in metastable equilibrium.
+- Add a function DGT() and a corresponding target in revisit() for
+ calculating the Gibbs energy of transformation of a system at constant
+ temperature, pressure and chemical activities of basis species.
+
DATA UPDATES:
- Add 148 liquid and 148 crystalline acyclic isoprenoids, polycyclic
Modified: pkg/inst/tests/test-revisit.R
===================================================================
--- pkg/inst/tests/test-revisit.R 2012-09-04 23:20:56 UTC (rev 8)
+++ pkg/inst/tests/test-revisit.R 2012-09-18 15:06:42 UTC (rev 9)
@@ -13,7 +13,7 @@
expect_error(revisit(d, "logact"), "for 'logact' target, loga.ref must be supplied")
})
-test_that("target calculations are working", {
+test_that("target calculations give expected results", {
r.cv <- revisit(d, "CV", plot.it=FALSE)
r.sd <- revisit(d, "sd", plot.it=FALSE)
r.shannon <- revisit(d, "shannon", plot.it=FALSE)
@@ -28,3 +28,32 @@
r.logact <- revisit(d, "logact", 3, plot.it=FALSE)
expect_equal(r.logact$ix, 71)
})
+
+test_that("DGT target gives zero at equilibrium and >0 not at equilibrium", {
+ # let's use n-alkanes
+ basis(c("CH4", "H2"), c("gas", "gas"))
+ species(c("methane", "ethane", "propane", "n-butane"), "liq")
+ # calculate equilibrium distribution over a range of logaH2
+ a <- affinity(H2=c(-10, -5, 101), exceed.Ttr=TRUE)
+ d <- diagram(a, plot.it=FALSE)
+ # the reference equilibrium distribution at logfH2 = -7.5
+ loga.ref <- list2array(d$logact)[51, ]
+ # calculate the DGT/RT relative to the reference distribution
+ # (we use swap12 so that the reference distribution is the initial one)
+ r <- revisit(d, "DGT", loga.ref=loga.ref, DGT.swap12=TRUE, plot.it=FALSE)
+ # we should find a minimum of zero at logfH2 = -7.5
+ expect_equal(min(r$H), 0)
+ expect_equal(a$xvals[which.min(r$H)], -7.5)
+ # we can also call the DGT function directly
+ # again with reference distribution as the initial one,
+ # but this time the Astar is kept constant (for logfH2 = -7.5)
+ Astar <- list2array(d$Astar)[51, ]
+ DGT.out <- numeric()
+ for(i in 1:length(a$xvals)) {
+ loga.equil <- list2array(d$logact)[i, ]
+ DGT.out <- c(DGT.out, DGT(loga.ref, loga.equil, Astar))
+ }
+ # we should find a minimum of zero at logfH2 = -7.5
+ expect_equal(min(DGT.out), 0)
+ expect_equal(a$xvals[which.min(DGT.out)], -7.5)
+})
Modified: pkg/man/revisit.Rd
===================================================================
--- pkg/man/revisit.Rd 2012-09-04 23:20:56 UTC (rev 8)
+++ pkg/man/revisit.Rd 2012-09-18 15:06:42 UTC (rev 9)
@@ -14,7 +14,7 @@
ispecies = NULL, add = FALSE, cex = par("cex"), lwd = par("lwd"),
mar = NULL, side = 1:4, xlim = NULL, labcex = 0.6, pch = 1,
legend = "", legend.x = NULL, lpch = NULL, main = NULL,
- lograt.ref = NULL, plot.ext = TRUE)
+ lograt.ref = NULL, plot.ext = TRUE, DGT.swap12 = FALSE)
extremes(z, target)
where.extreme(z, target, do.sat = FALSE)
}
@@ -43,7 +43,8 @@
\item{lograt.ref}{numeric, log10 of reference abundance ratios}
\item{plot.ext}{logical, show the location of the extreme value(s)?}
\item{z}{numeric, matrix of values}
- \item{do.sat}{logical, identify multiple extreme values}
+ \item{do.sat}{logical, identify multiple extreme values? }
+ \item{DGT.swap12}{logical, swap order of initial and final states?}
}
\details{
@@ -61,6 +62,7 @@
\samp{pearson} \tab Pearson correlation coefficient \tab max \tab \code{loga.ref} \cr
\samp{count} \tab count of species above a certain activity \tab max \tab \code{loga.ref} \cr
\samp{logact} \tab logarithm of activity of a single species \tab max \tab \code{loga.ref} \cr
+ \samp{DGT} \tab Gibbs energy of transformation of the system \tab min \tab {loga.ref} \cr
}
\code{sd}, \code{cv}, \code{shannon} and \code{qqr} all operate on the logarithms of activities of species.
@@ -69,6 +71,8 @@
\code{richness} counts the numbers of species whose logarithms of activities are above the single value of \code{loga.ref}. \code{logact} simply operates on the logarithm of activity of a single species (as supplied in \code{d}); \code{loga.ref} in this case is an integer value referring to the number of the species in the list in \code{d}.
+ \code{DGT} is the change in Gibbs energy of a system undergoing transformation among species of interest at constant T, P and chemical potentials of basis species. See \code{\link{DGT}} for the equation that is used. The transformation is from an initial state (the equilibrium values in \code{d}) to a final state (the values given in \code{loga.ref}); setting \code{DGT.swap12} reverses the order, so that the equilibrium state is the final state. Unlike the other targets, \code{DGT} also depends on the values of the starved affinity (\code{Astar}) that are returned by \code{\link{diagram}}.
+
Toether with the results of the calculation of the \code{target}, the function returns the position (in as many dimensions as the input values) and value of the extremum of the target value -- either the minimum or maximum, as indicated in the table.
If \code{plot.it} is TRUE, and \code{d} is the output from \code{diagram}, and the number of variables is 1 or 2, the results are plotted --- a line diagram in 1 dimension or a contour plot in 2 dimensions. No plot is made if the logarithms of activity have no dimensions (a single value for each species), except for the case of \code{qqr}, showing a quantile-quantile plot (\code{\link{qqnorm}}). On the plots the location of the extreme value is indicated by a dashed vertical line on a 1-D plot or a point marked by an asterisk on a 2-D plot. On 2-D plots the valleys (or ridges) leading to the location of the extremum are plotted. The ridges or valleys are plotted as dashed lines and are colored green for the \eqn{x}{x} values returned by \code{extremes} and blue for the \eqn{y}{y} values returned by \code{extremes}.
Modified: pkg/man/util.misc.Rd
===================================================================
--- pkg/man/util.misc.Rd 2012-09-04 23:20:56 UTC (rev 8)
+++ pkg/man/util.misc.Rd 2012-09-18 15:06:42 UTC (rev 9)
@@ -5,9 +5,10 @@
\alias{nonideal}
\alias{which.balance}
\alias{unitize}
+\alias{DGT}
\title{Functions for Miscellaneous Tasks}
\description{
- Calculate \eqn{dP/dT}{dP/dT} and temperature of phase transitions, calculate heat capacities of unfolded proteins using an equation from the literature, calculate non-ideal contributions to apparent standard molal properties, identify a conserved basis species, scale logarithms of activity to a desired total activity.
+ Calculate \eqn{dP/dT}{dP/dT} and temperature of phase transitions, calculate heat capacities of unfolded proteins using an equation from the literature, calculate non-ideal contributions to apparent standard molal properties, identify a conserved basis species, scale logarithms of activity to a desired total activity, calculate Gibbs energy of transformation of a system.
}
\usage{
@@ -16,6 +17,7 @@
nonideal(species, proptable, IS, T)
which.balance(species)
unitize(logact = NULL, length = NULL, logact.tot = 0)
+ DGT(loga1, loga2, Astar)
}
\arguments{
@@ -29,6 +31,9 @@
\item{logact}{numeric, logarithms of activity}
\item{length}{numeric, numbers of residues}
\item{logact.tot}{numeric, logarithm of total activity}
+ \item{loga1}{numeric, initial logarithms of activities of the species of interest}
+ \item{loga2}{numeric, final logarithms of activities of the species of interest}
+ \item{Astar}{numeric, the starved affinity}
}
\details{
@@ -45,6 +50,7 @@
\code{unitize} scales the logarithms of activities given in \code{logact} so that the logarithm of total activity of residues is equal to zero (i.e. total activity of residues is one), or to some other value set in \code{logact.tot}. \code{length} indicates the number of residues in each species. If \code{logact} is NULL, the function takes the logarithms of activities from the current species definition. If any of those species are proteins, the function gets their lengths using \code{protein.length}.
+ \code{DGT} calculates the change in Gibbs energy/2.303RT of a system in which species with initial logarithms of activitiy (\code{loga1}) are transformed to the same species with different final logarithms of activity (\code{loga2}) at constant temperature, pressure and chemical potentials of basis species. It is calculated as the sum over species of (G2-G1) where G1/RT = -a1*Astar + a1*loga1 - a1 + a constant (where a1 is 10^loga1), likewise for G2, and where \code{Astar} is the starved affinity, that is the affinity of the reaction to form one mole of the species at unit activity from the basis species in their defined activities. The equation used arises from integrating dG = -A/dxi = -A/dn where xi is the reaction progress variable, dn/dxi = 1 is the reaction coefficient on the species, and A = Astar - 2.303RTloga is the chemical affinity.
}
\value{
More information about the CHNOSZ-commits
mailing list