[CHNOSZ-commits] r22 - in pkg/CHNOSZ: . R inst inst/tests man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Sep 24 16:13:02 CEST 2012
Author: jedick
Date: 2012-09-24 16:13:01 +0200 (Mon, 24 Sep 2012)
New Revision: 22
Removed:
pkg/CHNOSZ/inst/FAQ
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/R/revisit.R
pkg/CHNOSZ/R/util.expression.R
pkg/CHNOSZ/R/util.misc.R
pkg/CHNOSZ/inst/NEWS
pkg/CHNOSZ/inst/tests/test-revisit.R
pkg/CHNOSZ/man/revisit.Rd
pkg/CHNOSZ/man/util.expression.Rd
pkg/CHNOSZ/man/util.misc.Rd
Log:
rename DGT() to DGxf() and add formatting of plotting symbols in expr.property()
expr.species() adds default subscripts for gases and charged aqueous species
remove inst/FAQ
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2012-09-23 15:11:26 UTC (rev 21)
+++ pkg/CHNOSZ/DESCRIPTION 2012-09-24 14:13:01 UTC (rev 22)
@@ -1,4 +1,4 @@
-Date: 2012-09-23
+Date: 2012-09-24
Package: CHNOSZ
Version: 0.9-7.98
Title: Chemical Thermodynamics and Activity Diagrams
Modified: pkg/CHNOSZ/R/revisit.R
===================================================================
--- pkg/CHNOSZ/R/revisit.R 2012-09-23 15:11:26 UTC (rev 21)
+++ pkg/CHNOSZ/R/revisit.R 2012-09-24 14:13:01 UTC (rev 22)
@@ -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", "dgt"))
+ if(tolower(target) %in% c("sd", "sd.log", "cv", "cv.log", "rmsd", "cvrmsd", "dgxf"))
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", "dgt"))
+ if(tolower(target) %in% c("sd", "sd.log", "cv", "cv.log", "rmsd", "cvrmsd", "dgxf"))
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, DGT.swap12=FALSE) {
+ pch=1, legend="", legend.x=NULL, lpch=NULL, main=NULL, lograt.ref=NULL, plot.ext=TRUE, DGxf.swap12=FALSE) {
# calculate and plot diversity indices of relative abundances
# 20090316 jmd
# d can be the output from diagram (enables plotting)
@@ -88,7 +88,7 @@
}
# check that all needed arguments are present
target.lower <- tolower(target)
- if(target.lower %in% c("richness", "cvrmsd", "spearman", "pearson", "logact", "dgt") & is.null(loga.ref))
+ if(target.lower %in% c("richness", "cvrmsd", "spearman", "pearson", "logact", "dgxf") & 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)
@@ -117,7 +117,7 @@
# "spearman" spearman correlation coefficient
# "pearson" pearson correlation coefficient
# "logact" maximize the activity of a species
- # "DGT" minimize the Gibbs energy of transformation
+ # "DGxf" minimize the Gibbs energy of transformation
# vectorize the entries in the logact list
for(i in 1:ns) {
@@ -241,7 +241,7 @@
# where the activity of a species is maximal
H <- logact[[loga.ref]]
- } else if(target.lower=="dgt") {
+ } else if(target.lower=="dgxf") {
# Gibbs energy of transformation to the observed assemblage
actarr <- list2array(logact)
# select species, vectorize, then put the Astar values into an array
@@ -254,8 +254,8 @@
# 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)
+ if(DGxf.swap12) out <- -DGxf(loga.ref, loga.equil, Astar)
+ else out <- DGxf(loga.equil, loga.ref, Astar)
return(out)
}
H <- as.numeric(palply(1:length(H), Gfun, actarr, Astararr))
@@ -283,12 +283,12 @@
# a 0-D plot
if(plot.it) {
plotted <- FALSE
- if(target=="qqr") {
+ if(target.lower=="qqr") {
actarr <- list2array(logact)
qqnorm(actarr,col=col,pch=pch,main="")
qqline(actarr)
plotted <- TRUE
- } else if(target %in% c("rmsd","cvrmsd","spearman","pearson")) {
+ } else if(target.lower %in% c("rmsd","cvrmsd","spearman","pearson")) {
# plot the points
ylab <- "loga.calc"
xlab <- "loga.ref"
@@ -323,7 +323,7 @@
# a 1-D plot
if(plot.it) {
if(is.null(ylim)) {
- if(target=="richness") {
+ if(target.lower=="richness") {
ylim <- 0
if(max(H) > ylim) ylim <- max(H) + 1
ylim <- c(0,ylim)
@@ -331,7 +331,10 @@
else ylim <- extendrange(H,f=0.075)
}
if(is.null(xlim)) xlim <- xrange
- if(!add) thermo.plot.new(xlim=xlim,ylim=ylim,xlab=axis.label(xname),ylab=target,yline=yline,
+ # format the target name if it's DGxf
+ if(target.lower=="dgxf") ylab <- expr.property("DGxf/2.303RT")
+ else ylab <- target
+ if(!add) thermo.plot.new(xlim=xlim,ylim=ylim,xlab=axis.label(xname),ylab=ylab,yline=yline,
cex=cex,lwd=lwd,mar=mar,side=side)
# plot the values
lines(xs,as.vector(H),col=col)
@@ -342,7 +345,7 @@
} else if(nd==2) {
# a 2-D plot
- iext <- where.extreme(H,target,do.sat=TRUE)
+ iext <- where.extreme(H,target.lower,do.sat=TRUE)
# what is the extreme value
ix <- iext$ix
iy <- iext$iy
@@ -366,7 +369,7 @@
# plot the location(s) of the extremum
points(xs[iext$ix],ys[iext$iy],pch=8,cex=2)
# show trajectories of the extrema
- iexts <- extremes(H,target)
+ iexts <- extremes(H,target.lower)
# take out large jumps
yext <- ys[iexts$y]
yext.1 <- c(yext[2:length(yext)],yext[length(yext)])
Modified: pkg/CHNOSZ/R/util.expression.R
===================================================================
--- pkg/CHNOSZ/R/util.expression.R 2012-09-23 15:11:26 UTC (rev 21)
+++ pkg/CHNOSZ/R/util.expression.R 2012-09-24 14:13:01 UTC (rev 22)
@@ -34,6 +34,9 @@
}
}
# write a designation of physical state
+ # use the state given in log if it's a gas or charged aqueous species
+ if(log %in% c("g", "gas")) state <- "g"
+ if("Z" %in% names(elements)) state <- log
if(state != "") {
# subscript it if we're not in a log expression
if(log != "") expr <- substitute(a*group('(',italic(b),')'),list(a=expr, b=state))
@@ -70,23 +73,27 @@
if(property=="IS") return(expression(IS))
if(property=="ZC") return(expression(bar(italic(Z))[C]))
# process each character in the property abbreviation
+ prevchar <- character()
for(i in 1:length(propchar)) {
+ if(i > 1) prevchar <- thischar
thischar <- propchar[i]
# D for greek Delta
# A for bold A (affinity)
# p for subscript italic P (in Cp)
- # 0 for degree sign
+ # 0 for degree sign (but not immediately following a number e.g. 2.303)
+ # f for subscript italic f (property of formation)
# r for subscript italic r (property of reaction)
- # f for subscript italic f (property of formation)
+ # x for subscript italic x (xf - property of transformation)
# all other letters are italicized
if(thischar %in% c(letters, LETTERS)) thisexpr <- substitute(italic(a), list(a=thischar))
else thisexpr <- substitute(a, list(a=thischar))
if(thischar=='D') thisexpr <- substitute(Delta)
if(thischar=='A') thisexpr <- substitute(bold(A))
if(thischar=='p') thisexpr <- substitute(a[italic(P)], list(a=""))
- if(thischar=='0') thisexpr <- substitute(degree)
+ if(thischar=='0' & !can.be.numeric(prevchar)) thisexpr <- substitute(degree)
if(thischar=='f') thisexpr <- substitute(a[italic(f)], list(a=""))
if(thischar=='r') thisexpr <- substitute(a[italic(r)], list(a=""))
+ if(thischar=='x') thisexpr <- substitute(a[italic(x)], list(a=""))
# put it together
expr <- substitute(a*b, list(a=expr, b=thisexpr))
}
Modified: pkg/CHNOSZ/R/util.misc.R
===================================================================
--- pkg/CHNOSZ/R/util.misc.R 2012-09-23 15:11:26 UTC (rev 21)
+++ pkg/CHNOSZ/R/util.misc.R 2012-09-24 14:13:01 UTC (rev 22)
@@ -141,7 +141,7 @@
# done!
}
-DGT <- function(loga1, loga2, Astar) {
+DGxf <- 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
@@ -156,7 +156,7 @@
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
+ DG <- G2 - G1
# return the sum
- return(sum(GT))
+ return(sum(DG))
}
Deleted: pkg/CHNOSZ/inst/FAQ
===================================================================
--- pkg/CHNOSZ/inst/FAQ 2012-09-23 15:11:26 UTC (rev 21)
+++ pkg/CHNOSZ/inst/FAQ 2012-09-24 14:13:01 UTC (rev 22)
@@ -1,24 +0,0 @@
-*************************
-start of a FAQ for CHNOSZ
-*************************
-
-- How do I load my own updated data files?
-
- Put the file in your working directory (use getwd() to find it).
- Then use add.obigt("mydatafile.csv").
- The file has to be in the same format as OBIGT.csv that comes with the package.
- Use system.file("data/OBIGT.csv", package="CHNOSZ") to find where that file is installed.
- But do not modify that file itself... it gets overwritten with package updates!
-
-- can I make a plot with log(aNa+/aH+) on an axis?
-
- Short answer: no, but plotting log(aNa+) would be numerically equivalent if you
- set aH+=1 (pH=0) (if your problem permits).
- Long answer: not easily; only logarithms of chemical activities of basis species
- (and temperature or pressure) can be the variables on the axes on plots made
- by diagram(). You could theoretically add a species with the composition
- of Na+ minus H+ to the database (it would have the formula "NaH-1+0" and the same
- properties as Na+) and use that as a basis species instead of Na+.
- Also, if you wanted an axis label like 'log(aNa+/aH+)' you'd have to add it
- manually, or modify the axis.label() function used by diagram().
-
Modified: pkg/CHNOSZ/inst/NEWS
===================================================================
--- pkg/CHNOSZ/inst/NEWS 2012-09-23 15:11:26 UTC (rev 21)
+++ pkg/CHNOSZ/inst/NEWS 2012-09-24 14:13:01 UTC (rev 22)
@@ -1,4 +1,4 @@
-CHANGES IN CHNOSZ 0.9-7.98 (2012-09-23)
+CHANGES IN CHNOSZ 0.9-7.98 (2012-09-24)
---------------------------------------
SIGNIFICANT USER-VISIBLE CHANGES:
Modified: pkg/CHNOSZ/inst/tests/test-revisit.R
===================================================================
--- pkg/CHNOSZ/inst/tests/test-revisit.R 2012-09-23 15:11:26 UTC (rev 21)
+++ pkg/CHNOSZ/inst/tests/test-revisit.R 2012-09-24 14:13:01 UTC (rev 22)
@@ -29,7 +29,7 @@
expect_equal(r.logact$ix, 71)
})
-test_that("DGT target gives zero at equilibrium and >0 not at equilibrium", {
+test_that("DGxf 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")
@@ -38,22 +38,22 @@
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
+ # calculate the DGxf/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)
+ r <- revisit(d, "DGxf", loga.ref=loga.ref, DGxf.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
+ # we can also call the DGxf 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()
+ DGxf.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))
+ DGxf.out <- c(DGxf.out, DGxf(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)
+ expect_equal(min(DGxf.out), 0)
+ expect_equal(a$xvals[which.min(DGxf.out)], -7.5)
})
Modified: pkg/CHNOSZ/man/revisit.Rd
===================================================================
--- pkg/CHNOSZ/man/revisit.Rd 2012-09-23 15:11:26 UTC (rev 21)
+++ pkg/CHNOSZ/man/revisit.Rd 2012-09-24 14:13:01 UTC (rev 22)
@@ -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, DGT.swap12 = FALSE)
+ lograt.ref = NULL, plot.ext = TRUE, DGxf.swap12 = FALSE)
extremes(z, target)
where.extreme(z, target, do.sat = FALSE)
}
@@ -44,7 +44,7 @@
\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{DGT.swap12}{logical, swap order of initial and final states?}
+ \item{DGxf.swap12}{logical, swap order of initial and final states?}
}
\details{
@@ -62,7 +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
+ \samp{DGxf} \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.
@@ -71,7 +71,7 @@
\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}}.
+ \code{DGxf} 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{DGxf}} 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{DGxf.swap12} reverses the order, so that the equilibrium state is the final state. Unlike the other targets, \code{DGxf} 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.
Modified: pkg/CHNOSZ/man/util.expression.Rd
===================================================================
--- pkg/CHNOSZ/man/util.expression.Rd 2012-09-23 15:11:26 UTC (rev 21)
+++ pkg/CHNOSZ/man/util.expression.Rd 2012-09-24 14:13:01 UTC (rev 22)
@@ -53,10 +53,14 @@
\samp{0} \tab degree sign (standard-state property) \cr
\samp{f} \tab subscript italic f (property of formation) \cr
\samp{r} \tab subscript italic r (property of reaction) \cr
+ \samp{x} \tab subscript italic x (property of transformation, xf) \cr
}
-Every other character that is one of the \code{\link{letters}} or \code{\link{LETTERS}} in the description of the property is italicized in the expression; other characters such as numerals or mathematical operators are shown without any special formatting. Special cases for the \code{property} argument (\samp{logK}, \samp{Eh}, \samp{pH}, \samp{pe}, \samp{IS} and \samp{ZC}) are interpreted as simple expressions, and are not parsed according to the above rules.
+A \samp{0} gets interpreted as a degree sign only if it does not immediately follow a number (so that e.g. \samp{2.303} can be included in an expression).
+Every other character that is one of the \code{\link{letters}} or \code{\link{LETTERS}} in the description of the property is italicized in the expression; other characters such as numerals or mathematical operators are shown without any special formatting.
+Special cases for the \code{property} argument (\samp{logK}, \samp{Eh}, \samp{pH}, \samp{pe}, \samp{IS} and \samp{ZC}) are interpreted as simple expressions, and are not parsed according to the above rules.
+
\code{expr.units} returns an expression for the units, based on one or more characters appearing in the \code{property}:
\tabular{ll}{
Modified: pkg/CHNOSZ/man/util.misc.Rd
===================================================================
--- pkg/CHNOSZ/man/util.misc.Rd 2012-09-23 15:11:26 UTC (rev 21)
+++ pkg/CHNOSZ/man/util.misc.Rd 2012-09-24 14:13:01 UTC (rev 22)
@@ -5,7 +5,7 @@
\alias{nonideal}
\alias{which.balance}
\alias{unitize}
-\alias{DGT}
+\alias{DGxf}
\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 Gibbs energy of transformation of a system.
@@ -17,7 +17,7 @@
nonideal(species, proptable, IS, T)
which.balance(species)
unitize(logact = NULL, length = NULL, logact.tot = 0)
- DGT(loga1, loga2, Astar)
+ DGxf(loga1, loga2, Astar)
}
\arguments{
@@ -50,7 +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.
+ \code{DGxf} 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