[CHNOSZ-commits] r96 - in pkg/CHNOSZ: . R data demo inst man tests/testthat
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sat Nov 7 09:53:26 CET 2015
Author: jedick
Date: 2015-11-07 09:53:26 +0100 (Sat, 07 Nov 2015)
New Revision: 96
Added:
pkg/CHNOSZ/R/nonideal.R
pkg/CHNOSZ/man/nonideal.Rd
Removed:
pkg/CHNOSZ/demo/CO2Ac.R
pkg/CHNOSZ/demo/diagram.R
pkg/CHNOSZ/demo/nonideal.R
pkg/CHNOSZ/demo/phosphate.R
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/R/diagram.R
pkg/CHNOSZ/R/equilibrate.R
pkg/CHNOSZ/R/examples.R
pkg/CHNOSZ/R/protein.info.R
pkg/CHNOSZ/R/species.R
pkg/CHNOSZ/R/swap.basis.R
pkg/CHNOSZ/R/util.expression.R
pkg/CHNOSZ/R/util.misc.R
pkg/CHNOSZ/R/util.plot.R
pkg/CHNOSZ/R/util.program.R
pkg/CHNOSZ/R/water.R
pkg/CHNOSZ/data/opt.csv
pkg/CHNOSZ/demo/00Index
pkg/CHNOSZ/demo/buffer.R
pkg/CHNOSZ/demo/copper.R
pkg/CHNOSZ/demo/nucleobase.R
pkg/CHNOSZ/demo/solubility.R
pkg/CHNOSZ/demo/wjd.R
pkg/CHNOSZ/inst/NEWS
pkg/CHNOSZ/man/CHNOSZ-package.Rd
pkg/CHNOSZ/man/buffer.Rd
pkg/CHNOSZ/man/data.Rd
pkg/CHNOSZ/man/examples.Rd
pkg/CHNOSZ/man/species.Rd
pkg/CHNOSZ/man/swap.basis.Rd
pkg/CHNOSZ/man/util.misc.Rd
pkg/CHNOSZ/man/util.plot.Rd
pkg/CHNOSZ/man/util.water.Rd
pkg/CHNOSZ/man/water.Rd
pkg/CHNOSZ/tests/testthat/test-makeup.R
Log:
create separate help page for nonideal()
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2015-11-05 05:18:37 UTC (rev 95)
+++ pkg/CHNOSZ/DESCRIPTION 2015-11-07 08:53:26 UTC (rev 96)
@@ -1,6 +1,6 @@
-Date: 2015-11-05
+Date: 2015-11-07
Package: CHNOSZ
-Version: 1.0.6-1
+Version: 1.0.6-2
Title: Chemical Thermodynamics and Activity Diagrams
Author: Jeffrey Dick
Maintainer: Jeffrey Dick <j3ffdick at gmail.com>
Modified: pkg/CHNOSZ/R/diagram.R
===================================================================
--- pkg/CHNOSZ/R/diagram.R 2015-11-05 05:18:37 UTC (rev 95)
+++ pkg/CHNOSZ/R/diagram.R 2015-11-07 08:53:26 UTC (rev 96)
@@ -222,7 +222,7 @@
myval <- sapply(plotvals, xfun)
ylim <- extendrange(myval)
}
- if(tplot) thermo.plot.new(xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab, cex=cex, mar=mar, yline=yline, side=side)
+ if(tplot) thermo.plot.new(xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab, cex=cex, mar=mar, yline=yline, side=side, ...)
else plot(0, 0, type="n", xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab, ...)
}
# draw the lines
Modified: pkg/CHNOSZ/R/equilibrate.R
===================================================================
--- pkg/CHNOSZ/R/equilibrate.R 2015-11-05 05:18:37 UTC (rev 95)
+++ pkg/CHNOSZ/R/equilibrate.R 2015-11-07 08:53:26 UTC (rev 96)
@@ -172,28 +172,24 @@
logadiff.max <- logadiff(Abar.range[2], i)
if(is.infinite(logadiff.max)) stop("FIXME: the second initial guess for Abar.max failed")
}
- # the change of logadiff with Abar
- # it's a weighted mean of the n.balance
- dlogadiff.dAbar <- (logadiff.max - logadiff.min) / diff(Abar.range)
- # change Abar to center logadiff (min/max) on zero
- logadiff.mean <- mean(c(logadiff.min, logadiff.max))
- Abar.range <- Abar.range - logadiff.mean / dlogadiff.dAbar
- # one iteration is enough for the examples in the package
- # but there might be a case where the range of logadiff doesn't cross zero
- # (e.g. for the carboxylic acid example in revisit.Rd)
- logadiff.min <- logadiff(Abar.range[1], i)
- logadiff.max <- logadiff(Abar.range[2], i)
- if(logadiff.min > 0 | logadiff.max < 0) {
- # do again what we did above
+ iter <- 0
+ while(logadiff.min > 0 | logadiff.max < 0) {
+ # the change of logadiff with Abar
+ # it's a weighted mean of the n.balance
dlogadiff.dAbar <- (logadiff.max - logadiff.min) / diff(Abar.range)
+ # change Abar to center logadiff (min/max) on zero
logadiff.mean <- mean(c(logadiff.min, logadiff.max))
Abar.range <- Abar.range - logadiff.mean / dlogadiff.dAbar
+ # one iteration is enough for the examples in the package
+ # but there might be a case where the range of logadiff doesn't cross zero
+ # (e.g. for the carboxylic acid example in revisit.Rd)
logadiff.min <- logadiff(Abar.range[1], i)
logadiff.max <- logadiff(Abar.range[2], i)
- if(logadiff.min > 0 | logadiff.max < 0) {
- stop("FIXME: make this function (Abarrange() in equil.reaction())
- iterate again to find a range of Abar such that the differences in
- logarithm of activity of the conserved component cross zero")
+ iter <- 1
+ if(iter > 5) {
+ stop("FIXME: we seem to be stuck! This function (Abarrange() in
+ equil.reaction()) can't find a range of Abar such that the differences
+ in logarithm of activity of the conserved component cross zero")
}
}
return(Abar.range)
Modified: pkg/CHNOSZ/R/examples.R
===================================================================
--- pkg/CHNOSZ/R/examples.R 2015-11-05 05:18:37 UTC (rev 95)
+++ pkg/CHNOSZ/R/examples.R 2015-11-07 08:53:26 UTC (rev 96)
@@ -30,8 +30,8 @@
}
demos <- function(which=c("sources", "NaCl", "density",
- "phosphate", "nucleobase", "ORP", "diagram", "revisit", "findit",
- "CO2Ac", "nonideal", "ionize", "buffer", "yeastgfp", "mosaic",
+ "nucleobase", "ORP", "revisit", "findit",
+ "ionize", "buffer", "yeastgfp", "mosaic",
"copper", "solubility", "wjd"), do.png=FALSE) {
# run one or more demos from CHNOSZ with ask=FALSE, and return the value of the last one
for(i in 1:length(which)) {
Added: pkg/CHNOSZ/R/nonideal.R
===================================================================
--- pkg/CHNOSZ/R/nonideal.R (rev 0)
+++ pkg/CHNOSZ/R/nonideal.R 2015-11-07 08:53:26 UTC (rev 96)
@@ -0,0 +1,68 @@
+# CHNOSZ/nonideal.R
+# first version of function: 20080308 jmd
+# moved to nonideal.R from util.misc.R 20151107
+
+nonideal <- function(species,proptable,IS,T) {
+ thermo <- get("thermo")
+ # generate nonideal contributions to thermodynamic properties
+ # number of species, same length as proptable list
+ # T in Kelvin, same length as nrows of proptables
+ # a function that does a lot of the work
+ loggamma2 <- function(Z,I,T,prop='log') {
+ # extended Debye-Huckel equation ('log')
+ # and its partial derivatives ('G','H','S','Cp')
+ # T in Kelvin
+ B <- 1.6
+ # equation for A from Clarke and Glew, 1980
+ #A <- expression(-16.39023 + 261.3371/T + 3.3689633*log(T)- 1.437167*(T/100) + 0.111995*(T/100)^2)
+ # equation for alpha from Alberty, 2003 p. 48
+ A <- alpha <- expression(1.10708 - 1.54508E-3 * T + 5.95584E-6 * T^2)
+ # from examples for deriv to take first and higher-order derivatives
+ DD <- function(expr,name, order = 1) {
+ if(order < 1) stop("'order' must be >= 1")
+ if(order == 1) D(expr,name)
+ else DD(D(expr, name), name, order - 1)
+ }
+ # Alberty, 2003 Eq. 3.6-1
+ loggamma <- function(a,Z,I,B) { - a * Z^2 * I^(1/2) / (1 + B * I^(1/2)) }
+ # XXX check the following equations 20080208 jmd
+ if(prop=='log') return(loggamma(eval(A),Z,I,B))
+ else if(prop=='G') return(thermo$opt$R * T * loggamma(eval(A),Z,I,B))
+ else if(prop=='H') return(thermo$opt$R * T^2 * loggamma(eval(DD(A,'T',1)),Z,I,B))
+ else if(prop=='S') return(- thermo$opt$R * T * loggamma(eval(DD(A,'T',1)),Z,I,B))
+ else if(prop=='Cp') return(thermo$opt$R * T^2 *loggamma(eval(DD(A,'T',2)),Z,I,B))
+ }
+ if(!is.numeric(species[[1]])) species <- info(species,'aq')
+ proptable <- as.list(proptable)
+ # which gamma function to use
+ #mlg <- get(paste('loggamma',which,sep=''))
+ ndid <- 0
+ for(i in 1:length(species)) {
+ myprops <- proptable[[i]]
+ # get the charge from the chemical formula
+ # force a charge count even if it's zero
+ mkp <- makeup(c("Z0", species[i]), sum=TRUE)
+ Z <- mkp[grep("Z", names(mkp))]
+ # don't do anything for neutral species
+ if(Z==0) next
+ # to keep unit activity coefficients of the proton and electron
+ if(species[i] == info("H+") & thermo$opt$ideal.H) next
+ if(species[i] == info("e-") & thermo$opt$ideal.e) next
+ didit <- FALSE
+ for(j in 1:ncol(myprops)) {
+ #if(colnames(myprops)[j]=='G') myprops[,j] <- myprops[,j] + thermo$opt$R * T * mlg(Z,IS,convert(T,'C'))
+ cname <- colnames(myprops)[j]
+ if(cname %in% c('G','H','S','Cp')) {
+ myprops[,j] <- myprops[,j] + loggamma2(Z,IS,T,cname)
+ didit <- TRUE
+ }
+ }
+ # append a loggam column if we did any nonideal calculations of thermodynamic properties
+ if(didit) myprops <- cbind(myprops,loggam=loggamma2(Z,IS,T))
+ proptable[[i]] <- myprops
+ if(didit) ndid <- ndid + 1
+ }
+ if(ndid > 0) msgout(paste('nonideal:',ndid,'species were nonideal\n'))
+ return(proptable)
+}
+
Modified: pkg/CHNOSZ/R/protein.info.R
===================================================================
--- pkg/CHNOSZ/R/protein.info.R 2015-11-05 05:18:37 UTC (rev 95)
+++ pkg/CHNOSZ/R/protein.info.R 2015-11-07 08:53:26 UTC (rev 96)
@@ -174,7 +174,7 @@
thermo <- get("thermo")
ionize.it <- FALSE
iword <- "nonionized"
- bmat <- basis.matrix()
+ bmat <- basis.elements()
if("H+" %in% rownames(bmat)) {
ionize.it <- TRUE
iword <- "ionized"
Modified: pkg/CHNOSZ/R/species.R
===================================================================
--- pkg/CHNOSZ/R/species.R 2015-11-05 05:18:37 UTC (rev 95)
+++ pkg/CHNOSZ/R/species.R 2015-11-07 08:53:26 UTC (rev 96)
@@ -2,9 +2,9 @@
# define species of interest
# to retrieve the coefficients of reactions to form the species from the basis species
-species.basis <- function(species) {
- # current basis matrix
- bmat <- basis.matrix()
+species.basis <- function(species=get("thermo")$species$ispecies) {
+ # current basis elements
+ bmat <- basis.elements()
tbmat <- t(bmat)
# what are the elements?
belem <- rownames(tbmat)
Modified: pkg/CHNOSZ/R/swap.basis.R
===================================================================
--- pkg/CHNOSZ/R/swap.basis.R 2015-11-05 05:18:37 UTC (rev 95)
+++ pkg/CHNOSZ/R/swap.basis.R 2015-11-07 08:53:26 UTC (rev 96)
@@ -2,8 +2,8 @@
# functions related to swapping basis species
# extracted from basis() 20120114 jmd
-# return the current basis matrix
-basis.matrix <- function(basis = get("thermo")$basis) {
+# return the current basis elements
+basis.elements <- function(basis = get("thermo")$basis) {
if(is.null(basis)) stop("basis species are not defined")
return(as.matrix(basis[, 1:nrow(basis), drop=FALSE]))
}
@@ -11,7 +11,7 @@
# calculate chemical potentials of elements from logarithms of activity of basis species
element.mu <- function(basis = get("thermo")$basis, T = 25) {
# matrix part of the basis definition
- basis.mat <- basis.matrix(basis)
+ basis.mat <- basis.elements(basis)
# the standard Gibbs energies of the basis species
if(T==25) G <- get("thermo")$obigt$G[basis$ispecies]
else G <- unlist(subcrt(basis$ispecies, T=T, property="G")$out)
@@ -27,7 +27,7 @@
# calculate logarithms of activity of basis species from chemical potentials of elements
basis.logact <- function(emu, basis = get("thermo")$basis, T = 25) {
# matrix part of the basis definition
- basis.mat <- basis.matrix(basis)
+ basis.mat <- basis.elements(basis)
# elements in emu can't be less than the number in the basis
if(length(emu) < ncol(basis.mat)) stop("number of elements in 'emu' is less than those in basis")
# sort names of emu in order of those in basis.mat
Modified: pkg/CHNOSZ/R/util.expression.R
===================================================================
--- pkg/CHNOSZ/R/util.expression.R 2015-11-05 05:18:37 UTC (rev 95)
+++ pkg/CHNOSZ/R/util.expression.R 2015-11-07 08:53:26 UTC (rev 96)
@@ -1,4 +1,4 @@
-# CHNOSZ/describe.R
+# CHNOSZ/util.expression.R
# write descriptions of chemical species, properties, reactions, conditions
# modified from describe(), axis.label() 20120121 jmd
@@ -7,8 +7,10 @@
# that include subscripts, superscripts (if charged)
# and optionally designations of states +/- loga or logf prefix
if(length(species) > 1) (stop("more than one species"))
- # the counts of elements in the species
- elements <- makeup(species)
+ # the counts of elements in the species:
+ # here we don't care too much if an "element" is a real element
+ # (listed in thermo$element), so we suppress warnings
+ elements <- suppressWarnings(makeup(species))
# where we'll put the expression
expr <- ""
# loop over elements
Modified: pkg/CHNOSZ/R/util.misc.R
===================================================================
--- pkg/CHNOSZ/R/util.misc.R 2015-11-05 05:18:37 UTC (rev 95)
+++ pkg/CHNOSZ/R/util.misc.R 2015-11-07 08:53:26 UTC (rev 96)
@@ -25,69 +25,6 @@
return(T + P / dPdT)
}
-nonideal <- function(species,proptable,IS,T) {
- thermo <- get("thermo")
- # generate nonideal contributions to thermodynamic properties
- # number of species, same length as proptable list
- # T in Kelvin, same length as nrows of proptables
- # a function that does a lot of the work
- loggamma2 <- function(Z,I,T,prop='log') {
- # extended Debye-Huckel equation ('log')
- # and its partial derivatives ('G','H','S','Cp')
- # T in Kelvin
- B <- 1.6
- # equation for A from Clarke and Glew, 1980
- #A <- expression(-16.39023 + 261.3371/T + 3.3689633*log(T)- 1.437167*(T/100) + 0.111995*(T/100)^2)
- # equation for alpha from Alberty, 2003 p. 48
- A <- alpha <- expression(1.10708 - 1.54508E-3 * T + 5.95584E-6 * T^2)
- # from examples for deriv to take first and higher-order derivatives
- DD <- function(expr,name, order = 1) {
- if(order < 1) stop("'order' must be >= 1")
- if(order == 1) D(expr,name)
- else DD(D(expr, name), name, order - 1)
- }
- # Alberty, 2003 Eq. 3.6-1
- loggamma <- function(a,Z,I,B) { - a * Z^2 * I^(1/2) / (1 + B * I^(1/2)) }
- # XXX check the following equations 20080208 jmd @@@
- if(prop=='log') return(loggamma(eval(A),Z,I,B))
- else if(prop=='G') return(thermo$opt$R * T * loggamma(eval(A),Z,I,B))
- else if(prop=='H') return(thermo$opt$R * T^2 * loggamma(eval(DD(A,'T',1)),Z,I,B))
- else if(prop=='S') return(- thermo$opt$R * T * loggamma(eval(DD(A,'T',1)),Z,I,B))
- else if(prop=='Cp') return(thermo$opt$R * T^2 *loggamma(eval(DD(A,'T',2)),Z,I,B))
- }
- if(!is.numeric(species[[1]])) species <- info(species,'aq')
- proptable <- as.list(proptable)
- # which gamma function to use
- #mlg <- get(paste('loggamma',which,sep=''))
- ndid <- 0
- for(i in 1:length(species)) {
- myprops <- proptable[[i]]
- # get the charge from the chemical formula
- # force a charge count even if it's zero
- mkp <- makeup(c("Z0", species[i]), sum=TRUE)
- Z <- mkp[grep("Z", names(mkp))]
- # don't do anything for neutral species
- if(Z==0) next
- # this keeps unit activity coefficient of the proton and electron
- if(species[i] %in% c(info("H+"), info("e-"))) next
- didit <- FALSE
- for(j in 1:ncol(myprops)) {
- #if(colnames(myprops)[j]=='G') myprops[,j] <- myprops[,j] + thermo$opt$R * T * mlg(Z,IS,convert(T,'C'))
- cname <- colnames(myprops)[j]
- if(cname %in% c('G','H','S','Cp')) {
- myprops[,j] <- myprops[,j] + loggamma2(Z,IS,T,cname)
- didit <- TRUE
- }
- }
- # append a loggam column if we did any nonideal calculations of thermodynamic properties
- if(didit) myprops <- cbind(myprops,loggam=loggamma2(Z,IS,T))
- proptable[[i]] <- myprops
- if(didit) ndid <- ndid + 1
- }
- if(ndid > 0) msgout(paste('nonideal:',ndid,'species were nonideal\n'))
- return(proptable)
-}
-
# which.balance is used by transfer(),
# keep it as a global function
which.balance <- function(species) {
Modified: pkg/CHNOSZ/R/util.plot.R
===================================================================
--- pkg/CHNOSZ/R/util.plot.R 2015-11-05 05:18:37 UTC (rev 95)
+++ pkg/CHNOSZ/R/util.plot.R 2015-11-07 08:53:26 UTC (rev 96)
@@ -200,11 +200,12 @@
}
}
-mtitle <- function(main,...) {
+mtitle <- function(main, line=0, ...) {
# make a possibly multi-line plot title
# useful for including expressions on multiple lines
+ # 'line' is the margin line of the last (bottom) line of the title
l <- length(main)
- for(i in 1:l) mtext(main[i],line=l-i,...)
+ for(i in 1:l) mtext(main[i], line=line+l-i, ...)
}
Modified: pkg/CHNOSZ/R/util.program.R
===================================================================
--- pkg/CHNOSZ/R/util.program.R 2015-11-05 05:18:37 UTC (rev 95)
+++ pkg/CHNOSZ/R/util.program.R 2015-11-07 08:53:26 UTC (rev 96)
@@ -27,9 +27,14 @@
# don't load methods package, to save startup time - ?makeCluster
cl <- parallel::makeCluster(nCores, methods=FALSE)
# export the variables and notify the user
- if(! "" %in% varlist) parallel::clusterExport(cl, varlist)
- msgout(paste("palply:", caller.name(4), "running", length(X), "calculations on",
- nCores, "cores with variable(s)", paste(varlist, collapse=", "), "\n"))
+ if(! "" %in% varlist) {
+ parallel::clusterExport(cl, varlist)
+ msgout(paste("palply:", caller.name(4), "running", length(X), "calculations on",
+ nCores, "cores with variable(s)", paste(varlist, collapse=", "), "\n"))
+ } else {
+ msgout(paste("palply:", caller.name(4), "running", length(X), "calculations on",
+ nCores, "cores\n"))
+ }
# run the calculations
out <- parallel::parLapply(cl, X, FUN, ...)
parallel::stopCluster(cl)
Modified: pkg/CHNOSZ/R/water.R
===================================================================
--- pkg/CHNOSZ/R/water.R 2015-11-05 05:18:37 UTC (rev 95)
+++ pkg/CHNOSZ/R/water.R 2015-11-07 08:53:26 UTC (rev 96)
@@ -285,7 +285,7 @@
# calculate values of P for Psat
if(identical(P, "Psat")) P <- psat()
msgout(" rho")
- my.rho <- rho.IAPWS95(T, P, get("thermo")$opt$IAPWS95.Psat.state)
+ my.rho <- rho.IAPWS95(T, P, get("thermo")$opt$IAPWS.sat)
rho <- function() my.rho
}
for(i in 1:length(property)) {
Modified: pkg/CHNOSZ/data/opt.csv
===================================================================
--- pkg/CHNOSZ/data/opt.csv 2015-11-05 05:18:37 UTC (rev 95)
+++ pkg/CHNOSZ/data/opt.csv 2015-11-07 08:53:26 UTC (rev 96)
@@ -1,2 +1,2 @@
-Tr,Pr,Theta,Psi,R,cutoff,E.units,T.units,P.units,state,water,G.tol,Cp.tol,V.tol,varP,IAPWS95.Psat.state,paramin
-298.15,1,228,2600,1.9872,1e-10,cal,C,bar,aq,SUPCRT92,100,1,1,FALSE,liquid,1000
+Tr,Pr,Theta,Psi,R,cutoff,E.units,T.units,P.units,state,water,G.tol,Cp.tol,V.tol,varP,IAPWS.sat,paramin,ideal.H,ideal.e
+298.15,1,228,2600,1.9872,1e-10,cal,C,bar,aq,SUPCRT92,100,1,1,FALSE,liquid,1000,TRUE,TRUE
Modified: pkg/CHNOSZ/demo/00Index
===================================================================
--- pkg/CHNOSZ/demo/00Index 2015-11-05 05:18:37 UTC (rev 95)
+++ pkg/CHNOSZ/demo/00Index 2015-11-07 08:53:26 UTC (rev 96)
@@ -1,14 +1,10 @@
sources cross-check the reference list with the thermodynamic database
NaCl equilibrium constant for aqueous NaCl dissociation
density density of H2O, inverted from IAPWS-95 equations
-phosphate phosphate speciation with pH, temperature and ionic strength
nucleobase relative stabilities of nucleobases and some amino acids
ORP oxidation-reduction potential of redox standards as a function of temperature
-diagram comparison of methods for calculating stability fields
revisit detailed example of usage of revisit()
findit detailed example of usage of findit()
-CO2Ac activity of CO2 buffered by acetic acid; comparing affinity(return.buffer=TRUE) with diagram(what="CO2")
-nonideal activity coefficient of charged species, using the IS argument of subcrt
ionize ionize.aa(): contour plots of net charge and ionization properties of LYSC_CHICK
buffer ionized proteins as a chemical activity buffer (1. thiol peroxidases 2. sigma factors)
yeastgfp logfO2-logaH2O diagrams for model proteins based on YeastGFP localizations
Deleted: pkg/CHNOSZ/demo/CO2Ac.R
===================================================================
--- pkg/CHNOSZ/demo/CO2Ac.R 2015-11-05 05:18:37 UTC (rev 95)
+++ pkg/CHNOSZ/demo/CO2Ac.R 2015-11-07 08:53:26 UTC (rev 96)
@@ -1,26 +0,0 @@
-# one can solve for the logarithm of activity of a
-# basis species using the 'what' argument of diagram
-basis("CHNOS")
-basis("CO2", 999)
-species("acetic acid", -3)
-a <- affinity(O2=c(-85, -70, 4), T=c(25, 100, 4))
-# hacking to write a title with formulas and subscripts
-lCO2 <- axis.label("CO2")
-main <- substitute(a~~b~~c,list(a=lCO2, b="buffered by",
- c="acetic acid"))
-d <- diagram(a, what="CO2", main=main)
-species(1, -10)
-a <- affinity(O2=c(-85, -70, 4), T=c(25, 100, 4))
-d <- diagram(a, what="CO2", add=TRUE, lty=2)
-# add a legend
-lAC <- expr.species("CH3COOH", log="aq")
-ltext <- c(as.expression(lAC), -3, -10)
-lty <- c(NA, 1, 2)
-legend("topright", legend=ltext, lty=lty, bg="white")
-# do return.buffer and diagram(what) give the same results?
-and <- as.numeric(d$plotvals[[1]])
-basis("CO2", "AC")
-mod.buffer("AC", logact=-10)
-a.buffer <- affinity(O2=c(-85, -70, 4), T=c(25, 100, 4), return.buffer=TRUE)
-ana <- as.numeric(unlist(a.buffer[[1]]))
-stopifnot(all.equal(ana, and))
Modified: pkg/CHNOSZ/demo/buffer.R
===================================================================
--- pkg/CHNOSZ/demo/buffer.R 2015-11-05 05:18:37 UTC (rev 95)
+++ pkg/CHNOSZ/demo/buffer.R 2015-11-07 08:53:26 UTC (rev 96)
@@ -11,11 +11,11 @@
basis(c("CO2", "H2O", "NH3", "O2"), "TPX")
a <- affinity(return.buffer=TRUE, T=50)
basis(c("CO2", "H2O", "NH3", "O2"), as.numeric(a[1:4]))
-a <- affinity(pH=c(0, 14, 200), T=c(25, 70, 200))
+a <- affinity(pH=c(4, 10, 300), T=c(40, 60, 300))
e <- equilibrate(a, normalize=TRUE)
diagram(e, fill=NULL)
title(main="Thiol peroxidases from bacteria")
-text(0.5, 66, describe.basis(thermo$basis[-6,], oneline=TRUE), adj=0)
+legend("topleft", describe.basis(thermo$basis[-6,]))
## Buffer + ionization: relative stabilities
## of E. coli sigma factors on a T-pH diagram
Modified: pkg/CHNOSZ/demo/copper.R
===================================================================
--- pkg/CHNOSZ/demo/copper.R 2015-11-05 05:18:37 UTC (rev 95)
+++ pkg/CHNOSZ/demo/copper.R 2015-11-07 08:53:26 UTC (rev 96)
@@ -61,8 +61,8 @@
# because the filling of fields masks it
water.lines()
box()
-title(main=paste("Aqueous Copper + Glycine, 25 deg C, 1 bar",
- "After Aksu and Doyle, 2001 Fig. 2b",sep="\n"))
+mtitle(expression("Copper-water-glycine at 25"~degree*"C and 1 bar",
+ "After Aksu and Doyle, 2001 (Fig. 2b)"), line=0.5)
# done!
data(thermo)
Deleted: pkg/CHNOSZ/demo/diagram.R
===================================================================
--- pkg/CHNOSZ/demo/diagram.R 2015-11-05 05:18:37 UTC (rev 95)
+++ pkg/CHNOSZ/demo/diagram.R 2015-11-07 08:53:26 UTC (rev 96)
@@ -1,17 +0,0 @@
-## a case where the maximum affinity method doesn't
-## reproduce an equilibrium predominance diagram
-basis("CHNOS+")
-# this adds data for some metabolites in the TCA cycle
-# from Dalla-Betta and Schulte, 2010
-add.obigt()
-species(c("oxaloacetate-2", "malate-2", "fumarate-2",
- "a-ketoglutarate-2", "citrate-3"))
-a <- affinity(O2=c(-80, -60), H2O=c(-5, 5))
-diagram(a, dotted=1, fill="heat")
-e <- equilibrate(a)
-diagram(e, add=TRUE, names=NULL, col="purple")
-e <- equilibrate(a, normalize=TRUE)
-diagram(e, add=TRUE, names=NULL)
-title(main=paste("maximum affinity method (fields)\n",
- "equilibrium calculations (lines)"))
-data(thermo)
Deleted: pkg/CHNOSZ/demo/nonideal.R
===================================================================
--- pkg/CHNOSZ/demo/nonideal.R 2015-11-05 05:18:37 UTC (rev 95)
+++ pkg/CHNOSZ/demo/nonideal.R 2015-11-07 08:53:26 UTC (rev 96)
@@ -1,40 +0,0 @@
-### non-ideality calculations -- activity coefficients of
-### aqueous species as a function of charge, temperature,
-### and ionic strength -- after Alberty, 2003
-## p. 16 Table 1.3 apparent pKa of acetic acid with
-## changing ionic strength
-subcrt(c("acetic acid","acetate","H+"),c(-1,1,1),
- IS=c(0,0.1,0.25),T=25,property="logK")
-# note that these *apparent* values of G and logK approach
-# their *standard* counterparts as IS goes to zero.
-## p. 95: basis and elemental stoichiometries of species
-## (a digression here from the nonideality calculations)
-# note coefficient of O2 and NH3 will be zero for these species
-basis(c("ATP-4","H+","H2O","HPO4-2","O2","NH3"))
-# cf Eq. 5.1-33: (basis composition)
-species(c("ATP-4","H+","H2O","HPO4-2","ADP-3","HATP-3","HADP-2","H2PO4-"))
-lb <- nrow(basis())
-# cf Eq. 5.1-32: (elemental composition)
-as.matrix(species()[,1:lb]) %*% as.matrix(basis()[,1:lb])
-## p. 273-275: activity coefficient (gamma)
-## as a function of ionic strength and temperature
-## (as of 20080304, these do look quantitatively different
-## from the plots in Alberty's book.)
-iplotfun <- function(T,col,add=TRUE) {
- IS <- seq(0,0.25,0.0025)
- s <- subcrt(c("H2PO4-","HADP-2","HATP-3","ATP-4"),IS=IS,grid="IS",T=T)
- if(!add) thermo.plot.new(xlim=range(IS),ylim=c(0,1),
- xlab=axis.label("IS"),ylab="gamma")
- for(i in 1:4) lines(IS,10^s$out[[i]]$loggam,col=col)
-}
-iplotfun(0,"blue",add=FALSE)
-text(0.1, 0.62, "Z = -1")
-iplotfun(25,"black")
-text(0.075, 0.18, "Z = -2")
-iplotfun(40,"red")
-text(0.05, 0.06, "Z = -3")
-title(main=paste("activity coefficient (gamma) of -1,-2,-3,-4",
- "charged species at 0, 25, 40 deg C, after Alberty, 2003",
- sep="\n"),cex.main=0.95)
-legend("topright", lty=c(NA, 1, 1, 1), col=c(NA, "blue", "black", "red"),
- legend=c(as.expression(axis.label("T")), 0, 25, 40))
Modified: pkg/CHNOSZ/demo/nucleobase.R
===================================================================
--- pkg/CHNOSZ/demo/nucleobase.R 2015-11-05 05:18:37 UTC (rev 95)
+++ pkg/CHNOSZ/demo/nucleobase.R 2015-11-07 08:53:26 UTC (rev 96)
@@ -1,21 +1,22 @@
## Nucleobase - Amino Acid Interaction Eh-H2O
-# for this example we try a unique basis definition
-basis(c("CO2","H2O","glutamine","e-","H+"),c(-3,0,-3,0,-7))
-species(c("uracil","cytosine","adenine","guanine",
- "phenylalanine","proline","lysine","glycine"),"aq")
+# for this example we try a custom basis set
+basis(c("CO2", "H2O", "glutamine", "e-", "H+"), c(-3, 0, -3, 0, -7))
+species(c("uracil", "cytosine", "adenine", "guanine",
+ "phenylalanine", "proline", "lysine", "glycine"), "aq")
# this loaded four nucleobases and four related amino acids
# (coded for by the homocodon triplets)
# check out the predominance diagrams
-a.1 <- affinity(H2O=c(-5,0),Eh=c(-0.5,0))
-diagram(a.1,fill=NULL)
+a.1 <- affinity(H2O=c(-5, 0, 300), Eh=c(-0.5, 0, 300))
+diagram(a.1, fill=NULL)
# overlay a different temperature
-a.2 <- affinity(H2O=c(-5,0),Eh=c(-0.5,0),T=100)
-diagram(a.2,col="red",add=TRUE,names=NULL)
-# start make a title for the plot
-tb <- thermo$basis # includes activities of basis species
+a.2 <- affinity(H2O=c(-5, 0, 300), Eh=c(-0.5, 0, 300), T=100)
+diagram(a.2, col="red", add=TRUE, names=NULL)
+# add title and legend
+title(main="Nucleobases and amino acids; P=Psat")
+# includes activities of basis species
+tb <- thermo$basis
# exclude those that are on the axes
-tb <- tb[!((rownames(tb) %in% c("e-","H2O"))),]
-title(main="Nucleobases and amino acids; P=Psat")
+tb <- tb[!((rownames(tb) %in% c("e-", "H2O"))),]
dp <- describe.property(c("T", "T"), c(25, 100))
db <- describe.basis(tb)
-legend("bottomleft", lty=c(1, 1, NA, NA, NA), col=c("black","red"), legend=c(dp, db))
+legend("bottomleft", lty=c(1, 1, NA, NA, NA), col=c("black", "red"), legend=c(dp, db))
Deleted: pkg/CHNOSZ/demo/phosphate.R
===================================================================
--- pkg/CHNOSZ/demo/phosphate.R 2015-11-05 05:18:37 UTC (rev 95)
+++ pkg/CHNOSZ/demo/phosphate.R 2015-11-07 08:53:26 UTC (rev 96)
@@ -1,32 +0,0 @@
-## speciation of phosphate as a function of ionic strength
-opar <- par(mfrow=c(2, 1))
-basis("CHNOPS+")
-T <- c(25, 100)
-species(c("PO4-3", "HPO4-2", "H2PO4-"))
-a25 <- affinity(IS=c(0, 0.14), T=T[1])
-e25 <- equilibrate(a25)
-d25 <- diagram(e25, ylim=c(-3.0, -2.6), legend.x=NULL)
-a100 <- affinity(IS=c(0, 0.14), T=T[2])
-e100 <- equilibrate(a100)
-d100 <- diagram(e100, ylim=c(-3.0, -2.6), add=TRUE, col="red")
-title(main="Non-ideality model for phosphate species")
-dp <- describe.property(c("pH", "T", "T"), c(7, T))
-legend("topright", lty=c(NA, 1, 1), col=c(NA, "black", "red"), legend=dp)
-text(0.07, -2.76, expr.species("HPO4-2"))
-text(0.07, -2.90, expr.species("H2PO4-"))
-# the crossing points of the logarithms of activity at the two temperatures
-# (it's higher IS at higher temperature)
-x25 <- which.min(abs(d25$loga.equil[[3]] - d25$loga.equil[[2]]))
-stopifnot(all.equal(x25, 27))
-x100 <- which.min(abs(d100$loga.equil[[3]] - d100$loga.equil[[2]]))
-stopifnot(all.equal(x100, 45))
-## phosphate predominance f(IS,pH)
-a <- affinity(IS=c(0, 0.14), pH=c(6, 13), T=T[1])
-d <- diagram(a, fill=NULL)
-a <- affinity(IS=c(0, 0.14), pH=c(6, 13), T=T[2])
-d <- diagram(a, add=TRUE, names=NULL, col="red")
-par(opar)
-# the most stable species are PO4-3, HPO4-2, H2PO4- with decreasing pH,
-# at any ionic strength (in the range used here)
-stopifnot(all.equal(unique(d$predominant[1, 1:128]), 3:1))
-stopifnot(all.equal(unique(d$predominant[128, 1:128]), 3:1))
Modified: pkg/CHNOSZ/demo/solubility.R
===================================================================
--- pkg/CHNOSZ/demo/solubility.R 2015-11-05 05:18:37 UTC (rev 95)
+++ pkg/CHNOSZ/demo/solubility.R 2015-11-07 08:53:26 UTC (rev 96)
@@ -74,10 +74,10 @@
# make plot
ylim <- c(-10, 4)
thermo.plot.new(xlim=range(pHs), ylim=ylim, xlab="pH", ylab="log a")
-lines(pHs, loga.tot, lwd=3)
+lines(pHs, loga.tot, lwd=4, col="green2")
lines(pHs, loga.CO2, lwd=2)
lines(pHs, loga.HCO3, lty=2, lwd=2)
lines(pHs, loga.CO3, lty=3, lwd=2)
-legend(ifelse(what=="calcite", "topright", "topleft"), lty=c(1, 1:3), lwd=c(3, 2, 2, 2),
+legend(ifelse(what=="calcite", "topright", "topleft"), lty=c(1, 1:3), lwd=c(4, 2, 2, 2), col=c("green2", rep("black", 3)),
legend=as.expression(c("total", expr.species("CO2", state="aq"), expr.species("HCO3-"), expr.species("CO3-2"))))
title(main=substitute(what~"solubility at"~T~degree*"C", list(what=what, T=T)))
Modified: pkg/CHNOSZ/demo/wjd.R
===================================================================
--- pkg/CHNOSZ/demo/wjd.R 2015-11-05 05:18:37 UTC (rev 95)
+++ pkg/CHNOSZ/demo/wjd.R 2015-11-07 08:53:26 UTC (rev 96)
@@ -19,13 +19,15 @@
# run the Gibbs energy minimization
w <- run.wjd(iobigt, Y=Y, imax=100)
# make a log-log plot
-plot(log10(y$abundance[!ina]), log10(w$X), xlim=c(1.5, 5), ylim=c(1.5, 5))
+plot(log10(y$abundance[!ina]), log10(w$X), xlim=c(1.5, 5), ylim=c(1.5, 5),
+ xlab="log10(abundance) reported in YeastGFP study",
+ ylab="log10(abundance) calculated using Gibbs energy minimization")
# get the element potentials (tolerating "close enough" to equilibrium)
emu <- equil.potentials(w, tol=1e7)
# then the logarithms of activities of the basis species
basis("CHNOS")
bl <- basis.logact(emu)
# make a title and legend
-title(main="calculated vs observed abundances: yeast cell periphery")
+title(main="relative abundances of proteins: yeast cell periphery")
basis(names(bl), bl)
legend("topleft", describe.basis(digits=2))
Modified: pkg/CHNOSZ/inst/NEWS
===================================================================
--- pkg/CHNOSZ/inst/NEWS 2015-11-05 05:18:37 UTC (rev 95)
+++ pkg/CHNOSZ/inst/NEWS 2015-11-07 08:53:26 UTC (rev 96)
@@ -1,10 +1,10 @@
-CHANGES IN CHNOSZ 1.0.6-1 (2015-11-05)
+CHANGES IN CHNOSZ 1.0.6-2 (2015-11-07)
--------------------------------------
- Add functions usrfig() (get figure limits in user coordinates) and
label.figure() (add label to figure outside of plot region).
-- Add copper.R demo (complexation of copper with glycine, uses
+- Add demo copper.R (complexation of copper with glycine, uses
mosaic()).
- Using new supporting function ibasis(), swap.basis() and mosaic()
@@ -17,6 +17,20 @@
concepts, organization of functions, and amino acid and protein
examples.
+- Remove demo diagram.R (concepts better shown in equilibrium.Rnw).
+
+- Rename basis.matrix() to basis.elements().
+
+- Argument of species.basis() defaults to current species indices.
+ ( species.basis() %*% basis.elements() returns elemental composition
+ of currently defined species.)
+
+- Move demos phosphate.R and nonideal.R to help page nonideal.Rd.
+
+- Add options to thermo$opt: 'ideal.H' and 'ideal.e'. Default TRUE tells
+ nonideal() to set activity coefficients of proton and electron to
+ zero.
+
CHANGES IN CHNOSZ 1.0.6 (2015-10-19)
------------------------------------
@@ -54,7 +68,7 @@
add specific importsFrom lines in NAMESPACE.
- Fix incorrect entry for entropy of aqueous methionine and [Met] in
- OBIGT.csv. Thanks to Apar Prasad for spotting this.
+ OBIGT.csv. Thanks to Apar Prasad for reporting this.
- Some fixes for compatibility with new version of testthat (increase
tolerance for one test in test-eos.R because of corrected alignment
Modified: pkg/CHNOSZ/man/CHNOSZ-package.Rd
===================================================================
--- pkg/CHNOSZ/man/CHNOSZ-package.Rd 2015-11-05 05:18:37 UTC (rev 95)
+++ pkg/CHNOSZ/man/CHNOSZ-package.Rd 2015-11-07 08:53:26 UTC (rev 96)
@@ -103,6 +103,8 @@
The values generated by \code{\link{buffer}} may not be applied correctly by \code{\link{affinity}} in calculating the affinities of the formation reactions. (The values returned by \code{affinity(..., return.buffer=TRUE)} do appear to be correct in the examples).
There is an unidentified inconsistency in \code{\link{transfer}} causing the reaction boundaries in one of the examples (\code{apc("closed")}) to be offset from the stability diagram. OTOH, \code{feldspar("closed")} appears to work correctly.
+
+ Values of activity coefficients may be affected by an unidentified bug in \code{\link{nonideal}}.
}
\examples{
Modified: pkg/CHNOSZ/man/buffer.Rd
===================================================================
--- pkg/CHNOSZ/man/buffer.Rd 2015-11-05 05:18:37 UTC (rev 95)
+++ pkg/CHNOSZ/man/buffer.Rd 2015-11-07 08:53:26 UTC (rev 96)
@@ -73,7 +73,7 @@
# change the activity of species in the buffer
mod.buffer("AC",logact=-10)
affinity(O2=c(-85,-70,4),T=c(25,100,4),return.buffer=TRUE)
-# see demo('CO2Ac') for a different strategy using the
+# see below for a different strategy using the
# 'what' argument of diagram
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/chnosz -r 96
More information about the CHNOSZ-commits
mailing list