[CHNOSZ-commits] r48 - in pkg/CHNOSZ: . R data inst inst/doc inst/tests man vignettes
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Mar 15 02:23:24 CET 2013
Author: jedick
Date: 2013-03-15 02:23:23 +0100 (Fri, 15 Mar 2013)
New Revision: 48
Removed:
pkg/CHNOSZ/R/util.supcrt.R
pkg/CHNOSZ/data/xxx.R
pkg/CHNOSZ/data/yyy.R
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/R/EOSregress.R
pkg/CHNOSZ/R/affinity.R
pkg/CHNOSZ/R/anim.R
pkg/CHNOSZ/R/basis.R
pkg/CHNOSZ/R/buffer.R
pkg/CHNOSZ/R/cgl.R
pkg/CHNOSZ/R/examples.R
pkg/CHNOSZ/R/findit.R
pkg/CHNOSZ/R/hkf.R
pkg/CHNOSZ/R/info.R
pkg/CHNOSZ/R/iprotein.R
pkg/CHNOSZ/R/makeup.R
pkg/CHNOSZ/R/protein.info.R
pkg/CHNOSZ/R/species.R
pkg/CHNOSZ/R/subcrt.R
pkg/CHNOSZ/R/swap.basis.R
pkg/CHNOSZ/R/transfer.R
pkg/CHNOSZ/R/util.affinity.R
pkg/CHNOSZ/R/util.args.R
pkg/CHNOSZ/R/util.data.R
pkg/CHNOSZ/R/util.expression.R
pkg/CHNOSZ/R/util.formula.R
pkg/CHNOSZ/R/util.misc.R
pkg/CHNOSZ/R/util.plot.R
pkg/CHNOSZ/R/util.seq.R
pkg/CHNOSZ/R/util.units.R
pkg/CHNOSZ/R/water.R
pkg/CHNOSZ/R/wjd.R
pkg/CHNOSZ/R/zzz.R
pkg/CHNOSZ/data/thermo.R
pkg/CHNOSZ/inst/NEWS
pkg/CHNOSZ/inst/doc/anintro.pdf
pkg/CHNOSZ/inst/doc/equilibrium.pdf
pkg/CHNOSZ/inst/doc/hotspring.pdf
pkg/CHNOSZ/inst/doc/wjd.pdf
pkg/CHNOSZ/inst/tests/test-EOSregress.R
pkg/CHNOSZ/inst/tests/test-affinity.R
pkg/CHNOSZ/inst/tests/test-basis.R
pkg/CHNOSZ/inst/tests/test-iprotein.R
pkg/CHNOSZ/inst/tests/test-protein.info.R
pkg/CHNOSZ/inst/tests/test-species.R
pkg/CHNOSZ/inst/tests/test-swap.basis.R
pkg/CHNOSZ/inst/tests/test-thermo.R
pkg/CHNOSZ/man/IAPWS95.Rd
pkg/CHNOSZ/man/anim.Rd
pkg/CHNOSZ/man/buffer.Rd
pkg/CHNOSZ/man/eos.Rd
pkg/CHNOSZ/man/iprotein.Rd
pkg/CHNOSZ/man/objective.Rd
pkg/CHNOSZ/man/protein.Rd
pkg/CHNOSZ/man/read.expr.Rd
pkg/CHNOSZ/man/revisit.Rd
pkg/CHNOSZ/man/sideeffects.Rd
pkg/CHNOSZ/man/species.Rd
pkg/CHNOSZ/man/subcrt.Rd
pkg/CHNOSZ/man/swap.basis.Rd
pkg/CHNOSZ/man/util.affinity.Rd
pkg/CHNOSZ/man/util.args.Rd
pkg/CHNOSZ/man/util.array.Rd
pkg/CHNOSZ/man/util.blast.Rd
pkg/CHNOSZ/man/util.character.Rd
pkg/CHNOSZ/man/util.data.Rd
pkg/CHNOSZ/man/util.expression.Rd
pkg/CHNOSZ/man/util.formula.Rd
pkg/CHNOSZ/man/util.matrix.Rd
pkg/CHNOSZ/man/util.misc.Rd
pkg/CHNOSZ/man/util.seq.Rd
pkg/CHNOSZ/man/util.units.Rd
pkg/CHNOSZ/man/water.Rd
pkg/CHNOSZ/man/wjd.Rd
pkg/CHNOSZ/vignettes/anintro.Rnw
pkg/CHNOSZ/vignettes/anintro.lyx
pkg/CHNOSZ/vignettes/equilibrium.Rnw
pkg/CHNOSZ/vignettes/equilibrium.lyx
pkg/CHNOSZ/vignettes/hotspring.Rnw
pkg/CHNOSZ/vignettes/hotspring.lyx
pkg/CHNOSZ/vignettes/wjd.Rnw
pkg/CHNOSZ/vignettes/wjd.lyx
Log:
follow "good practice" in ?data
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2013-03-14 01:06:37 UTC (rev 47)
+++ pkg/CHNOSZ/DESCRIPTION 2013-03-15 01:23:23 UTC (rev 48)
@@ -1,6 +1,6 @@
-Date: 2013-03-14
+Date: 2013-03-15
Package: CHNOSZ
-Version: 0.9-9.8
+Version: 0.9-9.9
Title: Chemical Thermodynamics and Activity Diagrams
Author: Jeffrey M. Dick
Maintainer: Jeffrey M. Dick <jmdick at asu.edu>
Modified: pkg/CHNOSZ/R/EOSregress.R
===================================================================
--- pkg/CHNOSZ/R/EOSregress.R 2013-03-14 01:06:37 UTC (rev 47)
+++ pkg/CHNOSZ/R/EOSregress.R 2013-03-15 01:23:23 UTC (rev 48)
@@ -6,16 +6,17 @@
EOSvar <- function(var, T, P) {
# get the variables of a term in a regression equation
# T (K), P (bar)
+ opt <- get("thermo")$opt
out <- switch(EXPR = var,
"(Intercept)" = rep(1, length(T)),
"T" = T,
"P" = P,
- "TTheta" = T-thermo$opt$Theta, # T-Theta
- "invTTheta" = (T-thermo$opt$Theta)^-1, # 1/(T-Theta)
- "TTheta2" = (T-thermo$opt$Theta)^2, # (T-Theta)^2
- "invTTheta2" = (T-thermo$opt$Theta)^-2, # 1/(T-Theta)^2
- "invPPsi" = (P+thermo$opt$Psi)^-1, # 1/(P-Psi)
- "invPPsiTTheta" = (P+thermo$opt$Psi)^-1 * (T-thermo$opt$Theta)^-1, # 1/[(P-Psi)(T-Theta)]
+ "TTheta" = T-opt$Theta, # T-Theta
+ "invTTheta" = (T-opt$Theta)^-1, # 1/(T-Theta)
+ "TTheta2" = (T-opt$Theta)^2, # (T-Theta)^2
+ "invTTheta2" = (T-opt$Theta)^-2, # 1/(T-Theta)^2
+ "invPPsi" = (P+opt$Psi)^-1, # 1/(P-Psi)
+ "invPPsiTTheta" = (P+opt$Psi)^-1 * (T-opt$Theta)^-1, # 1/[(P-Psi)(T-Theta)]
"TXBorn" = T*water("XBorn", T=T, P=P)[, 1],
"drho.dT" = -water("rho", T=T, P=P)[, 1]*water("E", T=T, P=P)[, 1],
"V.kT" = water("V", T=T, P=P)[, 1]*water("kT", T=T, P=P)[, 1],
Modified: pkg/CHNOSZ/R/affinity.R
===================================================================
--- pkg/CHNOSZ/R/affinity.R 2013-03-14 01:06:37 UTC (rev 47)
+++ pkg/CHNOSZ/R/affinity.R 2013-03-15 01:23:23 UTC (rev 48)
@@ -22,6 +22,7 @@
args <- c(args,list(sout=sout,exceed.Ttr=exceed.Ttr))
# the species we're given
+ thermo <- get("thermo")
mybasis <- thermo$basis
myspecies <- thermo$species
@@ -48,6 +49,7 @@
# account for protein activities later
resprot <- paste(resnames,"RESIDUE",sep="_")
species(resprot, 0)
+ thermo <- get("thermo", "CHNOSZ")
ires <- match(resprot, thermo$species$name)
}
@@ -59,6 +61,7 @@
msgout('affinity: loading buffer species\n')
if(!is.null(thermo$species)) is.species <- 1:nrow(thermo$species) else is.species <- numeric()
is.buffer <- buffer(logK=NULL)
+ thermo <- get("thermo", "CHNOSZ")
is.buff <- numeric()
for(i in 1:length(is.buffer)) is.buff <- c(is.buff,as.numeric(is.buffer[[i]]))
is.only.buffer <- is.buff[!is.buff %in% is.species]
Modified: pkg/CHNOSZ/R/anim.R
===================================================================
--- pkg/CHNOSZ/R/anim.R 2013-03-14 01:06:37 UTC (rev 47)
+++ pkg/CHNOSZ/R/anim.R 2013-03-15 01:23:23 UTC (rev 48)
@@ -7,10 +7,8 @@
# we depend on an empty png directory
if(!"png" %in% dir()) stop("directory 'png' not present")
else if(length(dir("png")) > 0) stop("directory 'png' not empty")
- # initialize the system
# add supplementary data (from default location of data/OBIGT-2.csv)
- # which includes properties from
- data(thermo)
+ # which includes properties for the metabolites
add.obigt()
# expand default logfO2 range if we're at high temperature
if(high.T & missing(redox)) redox <- list(O2=c(-100,-40))
@@ -95,7 +93,6 @@
notna <- !is.na(pdata$name)
pname <- pdata$name[notna]
# set up the system; use O2 aq instead of gas
- data(thermo)
basis(c("CO2","NH3","H2S","H2","O2","H+"))
basis("O2","aq")
basis(c("CO2","NH3","H2S","H+"),c(-3,-3,-10,-7))
@@ -180,7 +177,6 @@
ido <- c(rep(1,15),1:res,rep(res,20))
}
# set up system
- data(thermo)
basis(c("CO2","H2O","NH3","H2","H2S","H+"),
c("aq","liq","aq","aq","aq","aq"),c(-3,0,-4,-6,-7,-7))
species(c(rubisco,accoaco))
Modified: pkg/CHNOSZ/R/basis.R
===================================================================
--- pkg/CHNOSZ/R/basis.R 2013-03-14 01:06:37 UTC (rev 47)
+++ pkg/CHNOSZ/R/basis.R 2013-03-15 01:23:23 UTC (rev 48)
@@ -5,6 +5,7 @@
# to add the basis to thermo$obigt
put.basis <- function(ispecies, logact = rep(NA, length(ispecies))) {
+ thermo <- get("thermo")
state <- thermo$obigt$state[ispecies]
# make the basis matrix, revised 20120114
# get the elemental makeup of each species,
@@ -30,8 +31,9 @@
# store the basis definition in thermo$basis, including
# both numeric and character data, so we need to use a data frame
comp <- cbind(as.data.frame(comp), ispecies, logact, state, stringsAsFactors=FALSE)
- # assign to the global thermo object
- thermo$basis <<- comp
+ # ready to assign to the global thermo object
+ thermo$basis <- comp
+ assign("thermo", thermo, "CHNOSZ")
# remove the species since there's no guarantee the
# new basis includes all their elements
species(delete=TRUE)
@@ -40,6 +42,7 @@
# modify the states or logact values in the existing basis definition
mod.basis <- function(species, state=NULL, logact=NULL) {
+ thermo <- get("thermo")
# the basis must be defined
if(is.null(thermo$basis)) stop("basis is not defined")
# loop over each species to modify
@@ -70,18 +73,20 @@
"' are not in the basis\n",sep=""))
}
}
- thermo$basis$logact[ib] <<- state[i]
+ thermo$basis$logact[ib] <- state[i]
} else {
# look for a species with this name in the requested state
ispecies <- suppressMessages(info(thermo$obigt$name[thermo$basis$ispecies[ib]], state[i], check.it=FALSE))
if(is.na(ispecies) | is.list(ispecies))
stop(paste("state or buffer '", state[i], "' not found for ", thermo$obigt$name[thermo$basis$ispecies[ib]], "\n", sep=""))
- thermo$basis$ispecies[ib] <<- ispecies
- thermo$basis$state[ib] <<- state[i]
+ thermo$basis$ispecies[ib] <- ispecies
+ thermo$basis$state[ib] <- state[i]
}
}
# then modify the logact
- if(!is.null(logact)) thermo$basis$logact[ib] <<- as.numeric(logact[i])
+ if(!is.null(logact)) thermo$basis$logact[ib] <- as.numeric(logact[i])
+ # assign the result to the CHNOSZ environment
+ assign("thermo", thermo, "CHNOSZ")
}
return(thermo$basis)
}
@@ -124,9 +129,13 @@
## the actual basis() function
## delete, retrieve, define or modify the basis species of a thermodynamic system
basis <- function(species=NULL, state=NULL, logact=NULL, delete=FALSE) {
+ thermo <- get("thermo")
## delete the basis species if requested
oldbasis <- thermo$basis
- if(delete) thermo$basis <<- NULL
+ if(delete) {
+ thermo$basis <- NULL
+ assign("thermo", thermo, "CHNOSZ")
+ }
## return the basis definition if requested
if(is.null(species)) return(oldbasis)
## from now on we need something to work with
Modified: pkg/CHNOSZ/R/buffer.R
===================================================================
--- pkg/CHNOSZ/R/buffer.R 2013-03-14 01:06:37 UTC (rev 47)
+++ pkg/CHNOSZ/R/buffer.R 2013-03-15 01:23:23 UTC (rev 48)
@@ -2,8 +2,9 @@
# Calculate chemical activities of buffered species
# 20061102 jmd
-mod.buffer <- function(name,species=NULL,state=thermo$opt$state,logact=-3) {
+mod.buffer <- function(name,species=NULL,state=get("thermo")$opt$state,logact=-3) {
# 20071102 add or change a buffer system
+ thermo <- get("thermo")
if(is.null(species)) {
iname <- which(name==thermo$buffers$name)
if(length(iname)>0) species <- thermo$buffers$species[iname]
@@ -19,7 +20,8 @@
imod <- which(thermo$buffers$name %in% name & thermo$buffers$species %in% species)
if(length(imod)>0) {
if(state[1]=='') {
- thermo$buffers <<- thermo$buffers[-imod,]
+ thermo$buffers <- thermo$buffers[-imod,]
+ assign("thermo", thermo, "CHNOSZ")
msgout(paste('mod.buffer: removed ',c2s(species),' in ',
c2s(unique(name)),' buffer.\n',sep=''))
} else {
@@ -29,8 +31,9 @@
if(length(logact)!=ls) logact <- rep(logact,length.out=ls)
state.old <- thermo$buffers$state[imod]
logact.old <- thermo$buffers$logact[imod]
- thermo$buffers$state[imod] <<- state
- thermo$buffers$logact[imod] <<- logact
+ thermo$buffers$state[imod] <- state
+ thermo$buffers$logact[imod] <- logact
+ assign("thermo", thermo, "CHNOSZ")
if(identical(state.old,state) & identical(logact.old,logact)) {
msgout(paste('mod.buffer: nothing changed for ',
c2s(species),' in ',c2s(unique(name)),' buffer.\n',sep=''))
@@ -46,13 +49,15 @@
if(add) {
if(state[1]=='') state <- rep(thermo$opt$state,length.out=ls)
t <- data.frame(name=name,species=species,state=state,logact=logact)
- thermo$buffers <<- rbind(thermo$buffers,t)
+ thermo$buffers <- rbind(thermo$buffers,t)
+ assign("thermo", thermo, "CHNOSZ")
msgout(paste('mod.buffer: added ',c2s(unique(name)),'.\n',sep=''))
}
return(invisible(thermo$buffers[thermo$buffers$name %in% name,]))
}
buffer <- function(logK=NULL,ibasis=NULL,logact.basis=NULL,is.buffer=NULL,balance='PBB') {
+ thermo <- get("thermo")
# if logK is NULL load the buffer species
# otherwise perform buffer calculations.
if(is.null(logK)) {
Modified: pkg/CHNOSZ/R/cgl.R
===================================================================
--- pkg/CHNOSZ/R/cgl.R 2013-03-14 01:06:37 UTC (rev 47)
+++ pkg/CHNOSZ/R/cgl.R 2013-03-15 01:23:23 UTC (rev 48)
@@ -5,7 +5,9 @@
cgl <- function(property=NULL,T=298.15,P=1,ghs=NULL,eos=NULL) {
# calculate properties of crystalline, liquid (except H2O) and gas species
# argument handling
- Tr <- thermo$opt$Tr; Pr <- thermo$opt$Pr
+ thermo <- get("thermo")
+ Tr <- thermo$opt$Tr
+ Pr <- thermo$opt$Pr
eargs <- eos.args('mk',property=property)
prop <- eargs$prop
props <- eargs$props
Modified: pkg/CHNOSZ/R/examples.R
===================================================================
--- pkg/CHNOSZ/R/examples.R 2013-03-14 01:06:37 UTC (rev 47)
+++ pkg/CHNOSZ/R/examples.R 2013-03-15 01:23:23 UTC (rev 48)
@@ -25,7 +25,7 @@
}
if(is.character(do.png)) dev.off()
# at the end we attempt to restore the old par() (active as of the first call of thermo.plot.new)
- par(thermo$opar)
+ par(get("thermo")$opar)
cat("Time elapsed: ", proc.time() - .ptime, "\n")
}
Modified: pkg/CHNOSZ/R/findit.R
===================================================================
--- pkg/CHNOSZ/R/findit.R 2013-03-14 01:06:37 UTC (rev 47)
+++ pkg/CHNOSZ/R/findit.R 2013-03-15 01:23:23 UTC (rev 48)
@@ -34,7 +34,7 @@
}
# the initial values of the guesses (if midpoint==FALSE)
- basis <- thermo$basis
+ basis <- get("thermo")$basis
# a hack so that we can use pH as a variable
if("pH" %in% names(lims)) {
@@ -224,7 +224,7 @@
for(i in which) {
niter <- length(x$value[[i]])
ylab <- names(x$value)[i]
- if(ylab %in% c(rownames(thermo$basis),"T","P","pH","Eh")) ylab <- axis.label(ylab)
+ if(ylab %in% c(rownames(get("thermo")$basis),"T","P","pH","Eh")) ylab <- axis.label(ylab)
# the values
plot(1:niter,x$value[[i]],xlab=xlab,ylab=ylab,...)
lines(1:niter,x$value[[i]])
Modified: pkg/CHNOSZ/R/hkf.R
===================================================================
--- pkg/CHNOSZ/R/hkf.R 2013-03-14 01:06:37 UTC (rev 47)
+++ pkg/CHNOSZ/R/hkf.R 2013-03-15 01:23:23 UTC (rev 48)
@@ -6,6 +6,7 @@
H2O.PT=NULL,H2O.PrTr=NULL,domega=TRUE) {
# calculate G, H, S, Cp, V, kT, and/or E using
# the revised HKF equations of state
+ thermo <- get("thermo")
# constants
Tr <- thermo$opt$Tr
Pr <- thermo$opt$Pr
Modified: pkg/CHNOSZ/R/info.R
===================================================================
--- pkg/CHNOSZ/R/info.R 2013-03-14 01:06:37 UTC (rev 47)
+++ pkg/CHNOSZ/R/info.R 2013-03-15 01:23:23 UTC (rev 48)
@@ -9,7 +9,7 @@
info.text <- function(ispecies) {
# a textual description of species name, formula, source, e.g.
# CO2 [CO2(aq)] (SSW01, SHS89, 11.Oct.07)
- this <- thermo$obigt[ispecies, ]
+ this <- get("thermo")$obigt[ispecies, ]
sourcetext <- this$ref1
ref2 <- this$ref2
if(!is.na(ref2)) sourcetext <- paste(sourcetext, ref2, sep=", ")
@@ -24,6 +24,7 @@
# thermo$obigt$[species|abbrv|formula] or NA otherwise
# a match to thermo$obigt$state is also required if 'state' is not NULL
# (first occurence of a match to species is returned otherwise)
+ thermo <- get("thermo")
# all matches for the species
matches.species <- thermo$obigt$name==species |
thermo$obigt$abbrv==species | thermo$obigt$formula==species
@@ -49,6 +50,7 @@
eos <- aa2eos(aa, state)
# the real assignment work
nrows <- suppressMessages(mod.obigt(eos))
+ thermo <- get("thermo", "CHNOSZ")
matches.species <- rep(FALSE, nrows)
matches.species[nrows] <- TRUE
} else return(NA)
@@ -59,7 +61,7 @@
# special treatment for H2O: aq retrieves the liq
if(species %in% c("H2O", "water") & state=="aq") state <- "liq"
# the matches for both species and state
- matches.state <- matches.species & grepl(state, thermo$obigt$state)
+ matches.state <- matches.species & grepl(state, get("thermo")$obigt$state)
if(!any(matches.state)) {
# the requested state is not available for this species
available.states <- thermo$obigt$state[matches.species]
@@ -95,6 +97,7 @@
info.numeric <- function(ispecies, check.it=TRUE) {
# from a numeric species index in 'ispecies' return the
# thermodynamic properties and equations-of-state parameters
+ thermo <- get("thermo")
# if we're called with NA, return an empty row
if(is.na(ispecies)) {
this <- thermo$obigt[1,]
@@ -148,6 +151,7 @@
# returns species indices that have an approximate match of 'species'
# to thermo$obigt$[name|abbrv|formula],
# possibly restricted to a given state
+ thermo <- get("thermo")
if(!is.null(state)) this <- thermo$obigt[thermo$obigt$state==state, ]
else this <- thermo$obigt
# only look for fairly close matches
@@ -179,6 +183,7 @@
info <- function(species=NULL, state=NULL, check.it=TRUE) {
## return information for one or more species in thermo$obigt
## if no species are requested, summarize the available data 20101129
+ thermo <- get("thermo")
if(is.null(species)) {
msgout("info: 'species' is NULL; summarizing information about thermodynamic data...\n")
msgout(paste("thermo$obigt has", nrow(thermo$obigt[thermo$obigt$state=="aq", ]), "aqueous,",
Modified: pkg/CHNOSZ/R/iprotein.R
===================================================================
--- pkg/CHNOSZ/R/iprotein.R 2013-03-14 01:06:37 UTC (rev 47)
+++ pkg/CHNOSZ/R/iprotein.R 2013-03-15 01:23:23 UTC (rev 48)
@@ -15,6 +15,7 @@
# 'protein' numeric (the rownumber itself)
# 'protein' character, e.g. LYSC_CHICK
# 'protein' and 'organism', e.g. 'LYSC', 'CHICK'
+ thermo <- get("thermo")
if(is.numeric(protein)) {
iproteins <- 1:nrow(thermo$protein)
protein[!protein %in% iproteins] <- NA
@@ -42,13 +43,13 @@
iprotein <- iprotein(protein, organism)
# drop NA matches
iprotein <- iprotein[!is.na(iprotein)]
- out <- thermo$protein[iprotein, ]
+ out <- get("thermo")$protein[iprotein, ]
# compute per-residue counts
if(residue) out[, 5:25] <- out[, 5:25]/rowSums(out[, 6:25])
return(out)
}
-aa2eos <- function(aa, state=thermo$opt$state) {
+aa2eos <- function(aa, state=get("thermo")$opt$state) {
# display and return the properties of
# proteins calculated from amino acid composition
# the names of the protein backbone groups depend on the state
@@ -60,10 +61,11 @@
groups <- paste("[", groups, "]", sep="")
# the rownumbers of the groups in thermo$obigt
groups_state <- paste(groups, state)
- obigt_state <- paste(thermo$obigt$name, thermo$obigt$state)
+ obigt <- get("thermo")$obigt
+ obigt_state <- paste(obigt$name, obigt$state)
igroup <- match(groups_state, obigt_state)
# the properties are in columns 8-20 of thermo$obigt
- groupprops <- thermo$obigt[igroup, 8:20]
+ groupprops <- obigt[igroup, 8:20]
# the elements in each of the groups
groupelements <- i2A(igroup)
# a function to work on a single row of aa
@@ -150,7 +152,7 @@
read.aa <- function(file="protein.csv") {
# 20090428 added colClasses here
aa <- read.csv(file,colClasses=c(rep("character",4),rep("numeric",21)))
- if(!identical(colnames(aa), colnames(thermo$protein)))
+ if(!identical(colnames(aa), colnames(get("thermo")$protein)))
stop(paste("format of", file, "is incompatible with thermo$protein"))
return(aa)
}
@@ -158,16 +160,20 @@
add.protein <- function(aa) {
# add a properly constructed data frame of
# amino acid counts to thermo$protein
- if(!identical(colnames(aa), colnames(thermo$protein)))
+ thermo <- get("thermo")
+ if(!identical(colnames(aa), colnames(thermo$protein)))
stop("the value of 'aa' is not a data frame with the same columns as thermo$protein")
# find any protein IDs that are duplicated
po <- paste(aa$protein, aa$organism, sep="_")
ip <- suppressMessages(iprotein(po))
ipdup <- !is.na(ip)
# now we're ready to go
- if(!all(ipdup)) thermo$protein <<- rbind(thermo$protein, aa[!ipdup, ])
- if(any(ipdup)) thermo$protein[ip[ipdup], ] <<- aa[ipdup, ]
- rownames(thermo$protein) <<- NULL
+ tp.new <- thermo$protein
+ if(!all(ipdup)) tp.new <- rbind(tp.new, aa[!ipdup, ])
+ if(any(ipdup)) tp.new[ip[ipdup], ] <- aa[ipdup, ]
+ rownames(tp.new) <- NULL
+ thermo$protein <- tp.new
+ assign("thermo", thermo, "CHNOSZ")
# return the new rownumbers
ip <- iprotein(po)
# make some noise
Modified: pkg/CHNOSZ/R/makeup.R
===================================================================
--- pkg/CHNOSZ/R/makeup.R 2013-03-14 01:06:37 UTC (rev 47)
+++ pkg/CHNOSZ/R/makeup.R 2013-03-15 01:23:23 UTC (rev 48)
@@ -222,6 +222,7 @@
}
# if the formula argument is numeric, get the formula
# of that species number in thermo$obigt
+ thermo <- get("thermo")
if(is.numeric(formula)) formula <- thermo$obigt$formula[formula]
# first deal with charge
cc <- count.charge(formula)
Modified: pkg/CHNOSZ/R/protein.info.R
===================================================================
--- pkg/CHNOSZ/R/protein.info.R 2013-03-14 01:06:37 UTC (rev 47)
+++ pkg/CHNOSZ/R/protein.info.R 2013-03-15 01:23:23 UTC (rev 48)
@@ -96,6 +96,7 @@
G.Z <- rep(NA, length(pname))
ZC <- ZC(pf)
# run ionization calculations if we have H+
+ thermo <- get("thermo")
if(!is.null(thermo$basis)) {
iHplus <- match("H+", rownames(thermo$basis))
if(!is.na(iHplus)) {
@@ -143,6 +144,7 @@
# what are the coefficients of the basis species in the formation reactions
sb <- species.basis(pf)
# calculate ionization states if H+ is a basis species
+ thermo <- get("thermo")
iHplus <- match("H+", rownames(thermo$basis))
if(!is.na(iHplus)) {
pH <- -thermo$basis$logact[iHplus]
@@ -169,6 +171,7 @@
pname <- paste(aa$protein, aa$organism, sep="_")
plength <- protein.length(aa)
# use thermo$basis to decide whether to ionize the proteins
+ thermo <- get("thermo")
ionize.it <- FALSE
iword <- "nonionized"
bmat <- basis.matrix()
Modified: pkg/CHNOSZ/R/species.R
===================================================================
--- pkg/CHNOSZ/R/species.R 2013-03-14 01:06:37 UTC (rev 47)
+++ pkg/CHNOSZ/R/species.R 2013-03-15 01:23:23 UTC (rev 48)
@@ -44,6 +44,7 @@
# 20080925 default quiet=TRUE 20101003 default quiet=FALSE
# 20120128 remove 'quiet' argument (messages can be hidden with suppressMessages())
# 20120523 return thermo$species instead of rownumbers therein, and remove message showing thermo$species
+ thermo <- get("thermo")
## argument processing
# we can't deal with NA species
if(identical(species, NA)) {
@@ -55,7 +56,8 @@
if(delete) {
# delete the entire definition if requested
if(is.null(species)) {
- thermo$species <<- NULL
+ thermo$species <- NULL
+ assign("thermo", thermo, "CHNOSZ")
return(thermo$species)
}
# from here we're trying to delete already defined species
@@ -73,9 +75,10 @@
isp <- isp[!ina]
# go on to delete this/these species
if(length(isp) > 0) {
- thermo$species <<- thermo$species[-isp,]
- if(nrow(thermo$species)==0) thermo$species <<- NULL
- else rownames(thermo$species) <<- 1:nrow(thermo$species)
+ thermo$species <- thermo$species[-isp,]
+ if(nrow(thermo$species)==0) thermo$species <- NULL
+ else rownames(thermo$species) <- 1:nrow(thermo$species)
+ assign("thermo", thermo, "CHNOSZ")
}
return(thermo$species)
}
@@ -109,6 +112,8 @@
if(!any(is.na(ispecies)) & !is.null(logact)) return(species(ispecies, state=logact, index.return=index.return))
# look for species in thermo$obigt
iobigt <- suppressMessages(info(species, state))
+ # since that could have updated thermo$obigt (with proteins), re-read thermo
+ thermo <- get("thermo", "CHNOSZ")
# check if we got all the species
ina <- is.na(iobigt)
if(any(ina)) stop(paste("species not available:", paste(species[ina], collapse=" ")))
@@ -134,19 +139,19 @@
}
# create the new species
newspecies <- data.frame(f, ispecies=iobigt, logact=logact, state=state, name=name, stringsAsFactors=FALSE)
- # nasty for R, but "H2PO4-" looks better than "H2PO4."
+ # "H2PO4-" looks better than "H2PO4."
colnames(newspecies)[1:nrow(thermo$basis)] <- rownames(thermo$basis)
# initialize or add to species data frame
if(is.null(thermo$species)) {
- thermo$species <<- newspecies
+ thermo$species <- newspecies
ispecies <- 1:nrow(thermo$species)
} else {
# don't add species that already exist
idup <- newspecies$ispecies %in% thermo$species$ispecies
- thermo$species <<- rbind(thermo$species, newspecies[!idup, ])
+ thermo$species <- rbind(thermo$species, newspecies[!idup, ])
ispecies <- match(newspecies$ispecies, thermo$species$ispecies)
}
- rownames(thermo$species) <<- seq(nrow(thermo$species))
+ rownames(thermo$species) <- seq(nrow(thermo$species))
} else {
# update activities or states of existing species
# first get the rownumbers in thermo$species
@@ -160,7 +165,7 @@
} else ispecies <- match(species, thermo$species$name)
# replace activities?
if(!is.null(logact)) {
- thermo$species$logact[ispecies] <<- logact
+ thermo$species$logact[ispecies] <- logact
} else {
# change states, checking for availability of the desired state
for(i in 1:length(ispecies)) {
@@ -182,13 +187,14 @@
warning(paste("can't update state of species", ispecies[i], "to", state[i], "\n"), call.=FALSE)
else {
ii <- match(state[i], thermo$obigt$state[iobigt])
- thermo$species$state[ispecies[i]] <<- state[i]
- thermo$species$name[ispecies[i]] <<- thermo$obigt$name[iobigt[ii]]
- thermo$species$ispecies[ispecies[i]] <<- as.numeric(rownames(thermo$obigt)[iobigt[ii]])
+ thermo$species$state[ispecies[i]] <- state[i]
+ thermo$species$name[ispecies[i]] <- thermo$obigt$name[iobigt[ii]]
+ thermo$species$ispecies[ispecies[i]] <- as.numeric(rownames(thermo$obigt)[iobigt[ii]])
}
}
}
}
+ assign("thermo", thermo, "CHNOSZ")
# return the new species definition or the index(es) of affected species
if(index.return) return(ispecies)
else return(thermo$species)
Modified: pkg/CHNOSZ/R/subcrt.R
===================================================================
--- pkg/CHNOSZ/R/subcrt.R 2013-03-14 01:06:37 UTC (rev 47)
+++ pkg/CHNOSZ/R/subcrt.R 2013-03-15 01:23:23 UTC (rev 48)
@@ -96,6 +96,7 @@
}
# get species information
+ thermo <- get("thermo")
# pre-20110808, we sent numeric species argument through info() to
# get species name and state(s)
# but why slow things down if we already have a species index?
@@ -115,6 +116,8 @@
mysearch <- species[i]
if(can.be.numeric(mysearch)) mysearch <- thermo$obigt$name[as.numeric(mysearch)]
si <- info.character(mysearch, state[i])
+ # that could have the side-effect of adding a protein; re-read thermo
+ thermo <- get("thermo", "CHNOSZ")
if(is.na(si[1])) stop('no info found for ',species[i],' ',state[i])
if(!is.null(state[i])) is.cr <- state[i]=='cr' else is.cr <- FALSE
if(thermo$obigt$state[si[1]]=='cr1' & (is.null(state[i]) | is.cr)) {
Modified: pkg/CHNOSZ/R/swap.basis.R
===================================================================
--- pkg/CHNOSZ/R/swap.basis.R 2013-03-14 01:06:37 UTC (rev 47)
+++ pkg/CHNOSZ/R/swap.basis.R 2013-03-15 01:23:23 UTC (rev 48)
@@ -3,17 +3,17 @@
# extracted from basis() 20120114 jmd
# return the current basis matrix
-basis.matrix <- function(basis = thermo$basis) {
+basis.matrix <- function(basis = get("thermo")$basis) {
if(is.null(basis)) stop("basis species are not defined")
return(as.matrix(basis[, 1:nrow(basis), drop=FALSE]))
}
# calculate chemical potentials of elements from logarithms of activity of basis species
-element.mu <- function(basis = thermo$basis, T = 25) {
+element.mu <- function(basis = get("thermo")$basis, T = 25) {
# matrix part of the basis definition
basis.mat <- basis.matrix(basis)
# the standard Gibbs energies of the basis species
- if(T==25) G <- thermo$obigt$G[basis$ispecies]
+ if(T==25) G <- get("thermo")$obigt$G[basis$ispecies]
else G <- unlist(subcrt(basis$ispecies, T=T, property="G")$out)
# chemical potentials of the basis species
species.mu <- G - convert(basis$logact, "G")
@@ -25,7 +25,7 @@
}
# calculate logarithms of activity of basis species from chemical potentials of elements
-basis.logact <- function(emu, basis = thermo$basis, T = 25) {
+basis.logact <- function(emu, basis = get("thermo")$basis, T = 25) {
# matrix part of the basis definition
basis.mat <- basis.matrix(basis)
# elements in emu can't be less than the number in the basis
@@ -35,7 +35,7 @@
# check that elements of basis.mat and emu are identical
if(any(is.na(ielem))) stop(paste("element(s)", paste(names(emu)[is.na(ielem)], collapse=" "), "not found in basis"))
# the standard Gibbs energies of the basis species
- if(T==25) G <- thermo$obigt$G[basis$ispecies]
+ if(T==25) G <- get("thermo")$obigt$G[basis$ispecies]
else G <- unlist(subcrt(basis$ispecies, T=T, property="G")$out)
# the chemical potentials of the basis species in equilibrium
# with the chemical potentials of the elements
@@ -50,7 +50,7 @@
# swap in one basis species for another
swap.basis <- function(species, species2) {
# before we do anything, remember the old basis definition
- oldbasis <- thermo$basis
+ oldbasis <- get("thermo")$basis
# and the species definition
ts <- species()
# the delete the species
@@ -85,7 +85,7 @@
# restore species if they were defined
if(!is.null(ts)) {
suppressMessages(species(ts$ispecies))
- suppressMessages(species(1:nrow(thermo$species), ts$logact))
+ suppressMessages(species(1:nrow(get("thermo")$species), ts$logact))
}
# all done, return the new basis definition
return(mb)
Modified: pkg/CHNOSZ/R/transfer.R
===================================================================
--- pkg/CHNOSZ/R/transfer.R 2013-03-14 01:06:37 UTC (rev 47)
+++ pkg/CHNOSZ/R/transfer.R 2013-03-15 01:23:23 UTC (rev 48)
@@ -32,6 +32,7 @@
logpresent <- -50
# the starting (basis0) and current (basis)
# species and basis conditions
+ thermo <- get("thermo")
basis <- basis0 <- thermo$basis
species <- species0 <- thermo$species
dmode0 <- dmode
@@ -103,7 +104,8 @@
# the slow way, short version
# we have to use basis species w/o the PBB
tempbasis <- basis1[rownames(basis1)!="PBB",]
- thermo$basis <<- tempbasis
+ thermo$basis <- tempbasis
+ assign("thermo", thermo, "CHNOSZ")
# the slow way, long-winded version
# (recreating affinity's call to buffer so
# we can store the intermediate results)
@@ -135,8 +137,9 @@
#ibuf <- buffargs$ibasis
oldbasis <- thermo$basis
oldspecies <- thermo$species
- thermo$basis <<- buffstuff$bufbasis
- thermo$species <<- buffstuff$bufspecies
+ thermo$basis <- buffstuff$bufbasis
+ thermo$species <- buffstuff$bufspecies
+ assign("thermo", thermo, "CHNOSZ")
is.buffer <- buffargs$is.buffer
for(i in 1:length(ibuf)) {
ib <- is.buffer[[i]]
@@ -147,8 +150,9 @@
if(i==1) br <- bresult else br <- c(br,bresult)
}
bresult <- br
- thermo$basis <<- oldbasis
- thermo$species <<- oldspecies
+ thermo$basis <- oldbasis
+ thermo$species <- oldspecies
+ assign("thermo", thermo, "CHNOSZ")
}
for(i in 1:length(ibuf)) {
# reference to the moles of species
@@ -287,9 +291,10 @@
# get the affinities for the first step
getaff <- function(mybl,sout=NULL) {
# do it for unit activities of minerals (and proteins?)
- thermo$species$logact <<- 0
+ thermo$species$logact <- 0
# prevent the PBB from getting in here
- thermo$basis$logact <<- mybl[1:nrow(basis0)]
+ thermo$basis$logact <- mybl[1:nrow(basis0)]
+ assign("thermo", thermo, "CHNOSZ")
if(is.null(sout)) {
# on the first step only, grab the intermediate results
# they are kept around for reasons of speed
@@ -740,8 +745,9 @@
# this is the second place to be careful of PBB
if('PBB' %in% rownames(basis)) basis$logact <- c(basis0$logact,0)
else basis$logact <- basis0$logact
- thermo$basis <<- basis0
- thermo$species <<- species0
+ thermo$basis <- basis0
+ thermo$species <- species0
+ assign("thermo", thermo, "CHNOSZ")
# report the success rate and total progress
aaa <- alphas
@@ -800,7 +806,7 @@
# or feldspar("open")
# setup conditions for feldspar reaction
#basis(c('Al+3','SiO2','K+','H2O','H+','O2'))
- thermo$basis <<- NULL
+ basis(delete=TRUE)
# SLS89 use H4SiO4 instead of SiO2
basis(c('Al+3','H4SiO4','K+','H2O','H+','O2'))
# some of SLS89's initial conditions
@@ -849,7 +855,7 @@
# apc("many")
# apc("buffer")
# assign basis species
- thermo$basis <<- NULL
+ basis(delete=TRUE)
if(basis=="CO2") basis(c("CO2","H2O","NH3","H2","H2S"),c(-10,0,-4,-10,-7))
else if(basis=="acetic") basis(c("acetic acid","H2O","NH3","H2","H2S"),c(-5.5,0,-4,-10,-7))
basis("H2","aq")
Modified: pkg/CHNOSZ/R/util.affinity.R
===================================================================
--- pkg/CHNOSZ/R/util.affinity.R 2013-03-14 01:06:37 UTC (rev 47)
+++ pkg/CHNOSZ/R/util.affinity.R 2013-03-15 01:23:23 UTC (rev 48)
@@ -1,7 +1,7 @@
# CHNOSZ/util-affinity.R
# helper functions for affinity()
-energy <- function(what,vars,vals,lims,T=thermo$opt$Tr,P="Psat",IS=0,sout=NULL,exceed.Ttr=FALSE,transect=FALSE) {
+energy <- function(what,vars,vals,lims,T=get("thermo")$opt$Tr,P="Psat",IS=0,sout=NULL,exceed.Ttr=FALSE,transect=FALSE) {
# 20090329 extracted from affinity() and made to
# deal with >2 dimensions (variables)
@@ -21,7 +21,7 @@
mybasis <- basis()
nbasis <- nrow(mybasis)
## species definition / number of species
- myspecies <- thermo$species
+ myspecies <- get("thermo")$species
if(is.character(what)) {
if(is.null(myspecies)) stop('species properties requested, but species have not been defined')
nspecies <- nrow(myspecies)
@@ -222,6 +222,7 @@
# over which to calculate logQ, logK and affinity
# the names should be T, P, IS and names of basis species
# (or pH, pe, Eh)
+ thermo <- get("thermo")
## inputs are like c(T1,T2,res)
# and outputs are like seq(T1,T2,length.out=res)
# unless transect: do the variables specify a transect? 20090627
@@ -373,7 +374,7 @@
return(a)
}
-A.ionization <- function(iprotein, vars, vals, T=thermo$opt$Tr, P="Psat", pH=7, transect=FALSE) {
+A.ionization <- function(iprotein, vars, vals, T=get("thermo")$opt$Tr, P="Psat", pH=7, transect=FALSE) {
# a function to build a list of values of A/2.303RT of protein ionization
# that can be used by energy(); 20120527 jmd
# some of the variables might not affect the values (e.g. logfO2)
Modified: pkg/CHNOSZ/R/util.args.R
===================================================================
--- pkg/CHNOSZ/R/util.args.R 2013-03-14 01:06:37 UTC (rev 47)
+++ pkg/CHNOSZ/R/util.args.R 2013-03-15 01:23:23 UTC (rev 48)
@@ -8,7 +8,7 @@
# things we also get with water
props <- c(props,'A','U','Cv','Psat','rho','Q','X','Y','epsilon','w')
# they keep on coming: things we also get with SUPCRT92
- if(length(agrep(tolower(thermo$opt$water),'supcrt9',max.distance=0.3))>0)
+ if(get("thermo")$opt$water == "SUPCRT92")
props <- c(props,'Z','visc','tcond','tdiff','Prndtl','visck','albe','daldT','alpha','beta')
else
props <- c(props,'P','N','UBorn','de.dT','de.dP')
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/chnosz -r 48
More information about the CHNOSZ-commits
mailing list