[CHNOSZ-commits] r31 - in pkg/CHNOSZ: . R demo inst inst/tests man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sun Dec 9 15:09:44 CET 2012
Author: jedick
Date: 2012-12-09 15:09:44 +0100 (Sun, 09 Dec 2012)
New Revision: 31
Removed:
pkg/CHNOSZ/demo/copper.R
pkg/CHNOSZ/demo/pie.R
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/R/examples.R
pkg/CHNOSZ/R/info.R
pkg/CHNOSZ/R/species.R
pkg/CHNOSZ/R/swap.basis.R
pkg/CHNOSZ/R/util.affinity.R
pkg/CHNOSZ/R/zzz.R
pkg/CHNOSZ/demo/00Index
pkg/CHNOSZ/inst/NEWS
pkg/CHNOSZ/inst/tests/test-affinity.R
pkg/CHNOSZ/inst/tests/test-species.R
pkg/CHNOSZ/inst/tests/test-util.affinity.R
pkg/CHNOSZ/man/examples.Rd
Log:
cleanup species(), and fix array permutation in A.ionization()
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2012-11-27 14:04:22 UTC (rev 30)
+++ pkg/CHNOSZ/DESCRIPTION 2012-12-09 14:09:44 UTC (rev 31)
@@ -1,6 +1,6 @@
-Date: 2012-11-27
+Date: 2012-12-09
Package: CHNOSZ
-Version: 0.9.8-5
+Version: 0.9.8-6
Title: Chemical Thermodynamics and Activity Diagrams
Author: Jeffrey M. Dick
Maintainer: Jeffrey M. Dick <jmdick at asu.edu>
Modified: pkg/CHNOSZ/R/examples.R
===================================================================
--- pkg/CHNOSZ/R/examples.R 2012-11-27 14:04:22 UTC (rev 30)
+++ pkg/CHNOSZ/R/examples.R 2012-12-09 14:09:44 UTC (rev 31)
@@ -29,9 +29,9 @@
cat("Time elapsed: ", proc.time() - .ptime, "\n")
}
-demos <- function(which=c("sources", "NaCl", "copper", "cordierite",
- "phosphate", "nucleobase", "pie", "orp",
- "findit", "CO2Ac", "nonideal", "TPX")) {
+demos <- function(which=c("sources", "NaCl", "cordierite",
+ "phosphate", "nucleobase", "orp",
+ "findit", "CO2Ac", "nonideal")) {
# run one or more demos from CHNOSZ with ask=FALSE, and return the value of the last one
for(i in 1:length(which)) {
# say something so the user sees where we are
Modified: pkg/CHNOSZ/R/info.R
===================================================================
--- pkg/CHNOSZ/R/info.R 2012-11-27 14:04:22 UTC (rev 30)
+++ pkg/CHNOSZ/R/info.R 2012-12-09 14:09:44 UTC (rev 31)
@@ -218,7 +218,7 @@
out <- sapply(seq_along(species), function(i) {
# first look for exact match
ispecies <- info.character(species[i], state[i])
- # if no exact match and it's not a protein, show but do not use approximate matches
+ # if no exact match and it's not a protein, show approximate matches (side effect of info.approx)
if(identical(ispecies, NA) & !grepl("_", species[i])) ispecies.notused <- info.approx(species[i], state[i])
# do not accept multiple matches
if(length(ispecies) > 1) ispecies <- NA
Modified: pkg/CHNOSZ/R/species.R
===================================================================
--- pkg/CHNOSZ/R/species.R 2012-11-27 14:04:22 UTC (rev 30)
+++ pkg/CHNOSZ/R/species.R 2012-12-09 14:09:44 UTC (rev 31)
@@ -39,36 +39,30 @@
return(out)
}
+# to add to or alter the species definition
species <- function(species=NULL, state=NULL, delete=FALSE, index.return=FALSE) {
-# 20080925 changed default to quiet=TRUE
-# 20101003 changed default to 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
- missingstate <- missing(state)
- state <- state.args(state)
-
+ # 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
# we can't deal with NA species
if(identical(species, NA)) {
cn <- caller.name()
if(length(cn) > 0) ctext <- paste("(calling function was ", cn, ")", sep="") else ctext <- ""
stop(paste("'species' is NA", ctext))
}
-
# delete the entire species definition or only selected species
if(delete) {
- # remember the old species definition
- oldspecies <- thermo$species
# delete the entire definition if requested
if(is.null(species)) {
thermo$species <<- NULL
- return(oldspecies)
+ return(thermo$species)
}
# from here we're trying to delete already defined species
- if(is.null(oldspecies)) stop("nonexistent species definition")
+ if(is.null(thermo$species)) stop("nonexistent species definition")
# match species to delete by name and species number
isp <- rep(NA, length(species))
ispname <- match(species, thermo$species$name)
- ispnum <- match(species, 1:nrow(oldspecies))
+ ispnum <- match(species, 1:nrow(thermo$species))
isp[!is.na(ispname)] <- ispname[!is.na(ispname)]
isp[!is.na(ispnum)] <- ispnum[!is.na(ispnum)]
# filter out non-matching species
@@ -83,141 +77,107 @@
else rownames(thermo$species) <<- 1:nrow(thermo$species)
}
return(thermo$species)
- }
-
+ }
+ # parse state argument
+ state <- state.args(state)
# if no species or states are given, just return the species list
if(is.null(species) & is.null(state)) return(thermo$species)
# if no species are given use all of them if available
if(is.null(species) & !is.null(thermo$species)) species <- 1:nrow(thermo$species)
-
+ # make species and state arguments the same length
if(length(species) > length(state) & !is.null(state)) state <- rep(state,length.out=length(species)) else
if(length(state) > length(species) & !is.null(species)) species <- rep(species,length.out=length(state))
-
+ # rename state as logact if it's numeric, or get default values of logact
+ if(is.numeric(state[1])) {
+ logact <- state
+ state <- NULL
+ } else logact <- NULL
# if they don't look like states (aq,gas,cr) or activities (numeric),
- # use them as a suffix for species name (e.g., a protein-organism)
- if( length(which(state %in% unique(as.character(thermo$obigt$state)))) <
- length(state) & !can.be.numeric(state[1]) & !can.be.numeric(species[1]) ) {
- for(i in 1:length(state)) species[i] <- paste(species[i],'_',state[i],sep='')
- state <- rep(thermo$opt$state,length.out=length(state))
+ # use them as a suffix for species name (i.e., a protein_organism)
+ allstates <- unique(thermo$obigt$state)
+ if( sum(state %in% allstates) < length(state) & !can.be.numeric(state[1]) & !can.be.numeric(species[1]) ) {
+ species <- paste(species, state, sep="_")
+ state <- rep(thermo$opt$state, length.out=length(state))
}
-
- # append/change species entries
- if(is.null(thermo$basis)) stop('basis species are not defined')
+ # character first argument, look for species in thermo$obigt
+ iobigt <- NULL
if(is.character(species[1])) {
- # character first argument, species in thermo$obigt
- # but only give states if they are numeric
- is <- NULL
- if(!can.be.numeric(state[[1]])) is <- state
- ispecies <- suppressMessages(info(species, is))
+ iobigt <- suppressMessages(info(species, state))
# check if we got all the species
- ina <- is.na(ispecies)
- # info() returns a list if any of the species had multiple approximate matches
- # we don't accept any of those
- if(is.list(ispecies)) ina <- ina | sapply(ispecies,length) > 1
- if(any(ina)) stop(paste("species not available:",paste(species[ina],collapse=" ")))
- if(length(ispecies)==0) return(species())
- was.character <- TRUE
- } else if(is.numeric(species[1])) {
- ispecies <- species
- ispecies <- ispecies[!is.na(ispecies)]
- if(length(ispecies)==0) return(species())
- was.character <- FALSE
+ ina <- is.na(iobigt)
+ if(any(ina)) stop(paste("species not available:", paste(species[ina], collapse=" ")))
+ } else {
+ # if species is numeric and low number it refers to the index of existing species, else to thermo$obigt
+ nspecies <- nrow(thermo$species)
+ if(is.null(thermo$species)) nspecies <- 0
+ if(max(species) > nspecies) iobigt <- species
}
- jspecies <- ispecies
- myspecies <- NULL
- if(!is.null(thermo$species))
- if(TRUE %in% (ispecies %in% thermo$species$ispecies)) {
- myspecies <- match(ispecies,thermo$species$ispecies)
- myspecies <- thermo$species$ispecies[myspecies[!is.na(myspecies)]]
- ispecies <- ispecies[!ispecies%in%myspecies]
- }
- # only add to an existing dataframe if the indices can't
- # all possibly refer to the rows of the dataframe
- doit <- TRUE
- ## might (not) work well when you want to add e.g. H2O after some others
- if(!is.null(thermo$species)) if(all(jspecies %in% 1:nrow(thermo$species)))
- if(!was.character) doit <- FALSE
- if(length(ispecies) > 0 & !(is.numeric(ispecies[1]) & is.numeric(state[1]) ) & doit) {
+ # create or add to species definition
+ if(!is.null(iobigt)) {
+ if(is.null(thermo$basis)) stop("basis species are not defined")
# the coefficients in reactions to form the species from basis species
- f <- (species.basis(ispecies))
- # the default states and activities
- state <- as.character(thermo$obigt$state[ispecies])
- logact <- numeric()
- for(i in 1:length(state)) {
- if(length(agrep('missing',rownames(f)[i]))>0) la <- NA
- else { if(state[i]=='aq') la <- -3 else la <- 0 }
- logact <- c(logact,la)
+ f <- (species.basis(iobigt))
+ # the states and species names
+ state <- as.character(thermo$obigt$state[iobigt])
+ name <- as.character(thermo$obigt$name[iobigt])
+ # get default values of logact
+ if(is.null(logact)) {
+ logact <- rep(0, length(species))
+ logact[state=="aq"] <- -3
}
- # yes, the species names too
- name <- as.character(thermo$obigt$name[ispecies])
- # add the ispecies values
- t <- data.frame(f,ispecies=ispecies,logact=logact,state=state,name=name,stringsAsFactors=FALSE)
+ # 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."
- colnames(t)[1:nrow(thermo$basis)] <- rownames(thermo$basis)
- if(is.null(thermo$species)) thermo$species <<- t else
- thermo$species <<- rbind(thermo$species,t)
- rownames(thermo$species) <<- seq(1:nrow(thermo$species))
- }
- # update activities or states
- if(!is.null(state)) {
- state <- rep(state,length.out=length(jspecies))
- # if is looks like species aren't set yet, try to do so
- if(is.null(thermo$species)) {
- species(jspecies)
+ colnames(newspecies)[1:nrow(thermo$basis)] <- rownames(thermo$basis)
+ # initialize or add to species data frame
+ if(is.null(thermo$species)) {
+ thermo$species <<- newspecies
+ ispecies <- 1:nrow(thermo$species)
} else {
- mj <- jspecies[!jspecies %in% thermo$species$ispecies]
- if(!can.be.numeric(species[1])) species(mj)
+ # don't add species that already exist
+ idup <- newspecies$ispecies %in% thermo$species$ispecies
+ thermo$species <<- rbind(thermo$species, newspecies[!idup, ])
+ ispecies <- match(newspecies$ispecies, thermo$species$ispecies)
}
- # we bet that the number of rows is smaller
- # than the indices of whatever species we have
- if(can.be.numeric(species[1]) & max(jspecies) <= nrow(thermo$species))
- jspecies <- thermo$species$ispecies[jspecies]
- mj <- match(jspecies,thermo$species$ispecies)
- if(can.be.numeric(state[1])) {
- if(NA %in% mj[1]) warning(paste('can\'t update activity of species',
- c2s(which(is.na(mj))),' requested'),call.=FALSE)
- thermo$species$logact[mj] <<- state
+ rownames(thermo$species) <<- seq(nrow(thermo$species))
+ } else {
+ # update activities or states of existing species
+ # first get the rownumbers in thermo$species
+ if(is.numeric(species[1])) ispecies <- species
+ else ispecies <- match(species, thermo$species$name)
+ # replace activities?
+ if(!is.null(logact)) {
+ thermo$species$logact[ispecies] <<- logact
} else {
- mj <- match(jspecies,thermo$species$ispecies)
- state <- rep(state,length.out=length(mj))
- name <- thermo$species$name[mj]
- # try to check that the states actually exist
- for(k in 1:length(mj)) {
- doit <- TRUE
- if(NA %in% mj[k]) doit <- FALSE
- myform <- thermo$obigt$formula[thermo$species$ispecies[mj[k]]]
- #iobigt <- which(thermo$obigt$name==thermo$species$name[mj[k]] | thermo$obigt$formula==myform)
+ # change states, checking for availability of the desired state
+ for(i in 1:length(ispecies)) {
+ myform <- thermo$obigt$formula[thermo$species$ispecies[ispecies[i]]]
+ #iobigt <- which(thermo$obigt$name==thermo$species$name[ispecies[k]] | thermo$obigt$formula==myform)
# 20080925 don't match formula -- two proteins might have the
# same formula (e.g. YLR367W and YJL190C)
- #iobigt <- which(thermo$obigt$name==thermo$species$name[mj[k]])
+ #iobigt <- which(thermo$obigt$name==thermo$species$name[ispecies[k]])
# 20091112 do match formula if it's not a protein -- be able to
# change "carbon dioxide(g)" to "CO2(aq)"
- if(length(grep("_",thermo$species$name[mj[k]])) > 0)
- iobigt <- which(thermo$obigt$name==thermo$species$name[mj[k]])
+ if(length(grep("_",thermo$species$name[ispecies[i]])) > 0)
+ iobigt <- which(thermo$obigt$name==thermo$species$name[ispecies[i]])
else {
- iobigt <- which(thermo$obigt$name==thermo$species$name[mj[k]] & thermo$obigt$state==state[k])
+ iobigt <- which(thermo$obigt$name==thermo$species$name[ispecies[i]] & thermo$obigt$state==state[i])
if(length(iobigt)==0)
- iobigt <- which(thermo$obigt$name==thermo$species$name[mj[k]] | thermo$obigt$formula==myform)
+ iobigt <- which(thermo$obigt$name==thermo$species$name[ispecies[i]] | thermo$obigt$formula==myform)
}
- if(!state[k] %in% thermo$obigt$state[iobigt])
- doit <- FALSE
- if(!doit) warning(paste('can\'t update state of species ',
- mj[k],' to ',state[k],'.\n',sep=''),call.=FALSE)
+ if(!state[i] %in% thermo$obigt$state[iobigt])
+ warning(paste("can't update state of species", ispecies[i], "to", state[i], "\n"), call.=FALSE)
else {
- ii <- match(state[k],thermo$obigt$state[iobigt])
- thermo$species$state[mj[k]] <<- state[k]
- thermo$species$name[mj[k]] <<- thermo$obigt$name[iobigt[ii]]
- thermo$species$ispecies[mj[k]] <<- as.numeric(rownames(thermo$obigt)[iobigt[ii]])
+ 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]])
}
}
}
- } else {
- # this message turns out to be kinda distracting. what should be here?
- #if(!is.null(myspecies))
- #cat(paste('species: keeping ',c2s(thermo$obigt$name[myspecies],sep=', '),'.\n',sep=''))
}
- # return the new species definition or index(es) of affected species
- if(index.return) return(match(jspecies, thermo$species$ispecies))
+ # 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/swap.basis.R
===================================================================
--- pkg/CHNOSZ/R/swap.basis.R 2012-11-27 14:04:22 UTC (rev 30)
+++ pkg/CHNOSZ/R/swap.basis.R 2012-12-09 14:09:44 UTC (rev 31)
@@ -52,7 +52,9 @@
# before we do anything, remember the old basis definition
oldbasis <- thermo$basis
# and the species definition
- ts <- species(delete=TRUE)
+ ts <- species()
+ # the delete the species
+ species(delete=TRUE)
if(is.null(oldbasis))
stop("swapping basis species requires an existing basis definition")
# both arguments must have length 1
Modified: pkg/CHNOSZ/R/util.affinity.R
===================================================================
--- pkg/CHNOSZ/R/util.affinity.R 2012-11-27 14:04:22 UTC (rev 30)
+++ pkg/CHNOSZ/R/util.affinity.R 2012-12-09 14:09:44 UTC (rev 31)
@@ -386,7 +386,11 @@
if(transect) TPpH <- list(T=T, P=P, pH=pH)
else {
# make a grid of all combinations
- TPpH <- expand.grid(T=T, P=P, pH=pH, stringsAsFactors=FALSE)
+ # put T, P, pH in same order as they show up vars
+ egargs <- list(T=T, P=P, pH=pH)
+ TPpHorder <- order(c(iT, iP, iHplus))
+ egargs <- c(egargs[TPpHorder], list(stringsAsFactors=FALSE))
+ TPpH <- do.call(expand.grid, egargs)
# figure out the dimensions of T-P-pH (making sure to drop any that aren't in vars)
TPpHdim <- numeric(3)
TPpHdim[iT] <- length(T)
Modified: pkg/CHNOSZ/R/zzz.R
===================================================================
--- pkg/CHNOSZ/R/zzz.R 2012-11-27 14:04:22 UTC (rev 30)
+++ pkg/CHNOSZ/R/zzz.R 2012-12-09 14:09:44 UTC (rev 31)
@@ -14,22 +14,6 @@
date <- um[nchar(um)>0][2]
# identify the program and version
packageStartupMessage(paste("CHNOSZ version ", version, " (", date, ")", sep=""))
- # message if 'parallel' package is available but not loaded
- if(!"parallel" %in% sessionInfo()$basePkgs)
- if(length(find.package("parallel", quiet=TRUE)) > 0)
- packageStartupMessage("Suggested package 'parallel' available but not loaded")
# load the 'thermo' data object
data(thermo)
-
- ## load data files in user's directory
- #if(file.exists('protein.csv')) {
- # tr <- try(rbind(read.csv('protein.csv',as.is=TRUE),thermo$protein),silent=TRUE)
- # if(identical(class(tr),'try-error')) cat("thermo: protein.csv in current directory is not compatible with thermo$protein data table.\n")
- # else add.protein("protein.csv")
- #}
- #if(file.exists('obigt.csv')) {
- # tr <- try(rbind(read.csv('obigt.csv',as.is=TRUE),thermo$obigt),silent=TRUE)
- # if(identical(class(tr),'try-error')) cat("thermo: obigt.csv in current directory is not compatible with thermo$obigt data table.\n")
- # else add.obigt("obigt.csv")
- #}
}
Modified: pkg/CHNOSZ/demo/00Index
===================================================================
--- pkg/CHNOSZ/demo/00Index 2012-11-27 14:04:22 UTC (rev 30)
+++ pkg/CHNOSZ/demo/00Index 2012-12-09 14:09:44 UTC (rev 31)
@@ -1,10 +1,8 @@
sources cross-check the reference list with the thermodynamic database
NaCl equilibrium constant for aqueous NaCl dissociation
-copper an Eh-pH diagram for the copper-water-glycine system
cordierite equilibrium constant of hydrous cordierite dehydration
phosphate phosphate speciation with pH, temperature and ionic strength
nucleobase relative stabilities of nucleobases and some amino acids
-pie pie charts comparing relative abundances of organisms and model proteins in metastable equilibrium
orp oxidation-reduction potential of redox standards as a function of temperature
findit example of findit() using log-normal distribution as an objective
CO2Ac activity of CO2 buffered by acetic acid; comparing affinity(return.buffer=TRUE) with diagram(what="CO2")
Deleted: pkg/CHNOSZ/demo/copper.R
===================================================================
--- pkg/CHNOSZ/demo/copper.R 2012-11-27 14:04:22 UTC (rev 30)
+++ pkg/CHNOSZ/demo/copper.R 2012-12-09 14:09:44 UTC (rev 31)
@@ -1,71 +0,0 @@
-## Eh-pH diagrams for copper-water-glycine
-## After Fig. 2 of Aksu and Doyle, 2001
-## (Aksu, S. and Doyle, F. M., 2001. Electrochemistry of copper in aqueous glycine
-## solutions. J. Electrochem. Soc., 148, B51-B57. doi:10.1149/1.1344532)
-## We need to add some species and change some Gibbs energies.
-# update rows of the database
-i <- info(c("Cu(Gly)+","glycinate","e-","H+"))
-# this was written before mod.obigt() was around... could use that instead
-n <- nrow(thermo$obigt <<- rbind(thermo$obigt,thermo$obigt[rep(i[1],2),]))
-thermo$obigt$name[n-1] <<- "Cu(Gly)2-"
-thermo$obigt$name[n] <<- "HCu(Gly)+2"
-thermo$obigt$formula[n-1] <<- as.chemical.formula(makeup(c(i[1],i[2],i[3]), sum=TRUE))
-thermo$obigt$formula[n] <<- as.chemical.formula(makeup(c(i[1],i[4]), sum=TRUE))
-# In Fig 2b, total log activities of Cu (Cu_T)
-# and glycine (L_T) are -4 and -1
-basis(c("Cu+2","H2O","H+","e-","glycine","CO2"),c(999,0,999,999,-1,999))
-# solid species
-species(c("copper","cuprite","tenorite"))
-# aqueous species
-species(c("glycinium","glycine","glycinate","Cu+","Cu+2","CuO2-2","HCuO2-",
- "Cu(Gly)+","Cu(Gly)2","Cu(Gly)2-","HCu(Gly)+2"),-4)
-ispecies <- species()$ispecies
-# update the Gibbs energies using A&D's Table 1 and Table II
-logK <- c(convert(convert(c(0,-146,-129.7,-384.061,-370.647,-314.833,
- 49.98,65.49,-183.6,-258.5,-298.2)*1000,"cal"),"logK"),15.64,10.1,2.92)
-# do it in order so later species take account of prev. species' values
-for(i in 1:length(logK)) {
- G <- convert(logK[i],"G")
- if(i==12) G <- G + thermo$obigt$G[ispecies[8]] +
- 2*thermo$obigt$G[ispecies[6]]
- if(i==13) G <- G + thermo$obigt$G[ispecies[7]] +
- 2*thermo$obigt$G[ispecies[6]]
- if(i==14) G <- G + thermo$obigt$G[ispecies[11]]
- thermo$obigt$G[ispecies[i]] <- G
-} # done with changing Gibbs free energies!
-# we have to get some leftovers out of there or diagram() gets confused
-species(c("glycinium","glycine","glycinate"),delete=TRUE)
-# make a plot to see if it's working
-ispecies <- ispecies[-(1:6)]
-afun <- function(cu,gly) {
- # from above: our fifth basis species is glycine(-ate,-ium)
- basis(rownames(basis())[5],gly)
- t <- match(ispecies,species()$ispecies)
- species(t,cu)
- affinity(pH=c(0,16),Eh=c(-0.6,1.0))
-}
-diagram(afun(-4,-1))
-title(main=paste("Aqueous Copper + Glycine, 25 deg C, 1 bar",
- "After Aksu and Doyle, 2001 Fig. 2b",sep="\n"))
-# What's missing? Try glycinate not glycine in reactions at ph > ~9.8
-basis(c("Cu+2","H2O","H+","e-","glycinate","CO2"),
- c(999,0,999,999,-2,999))
-species(c("copper","cuprite","tenorite","Cu+","Cu+2","CuO2-2","HCuO2-",
- "Cu(Gly)+","Cu(Gly)2","Cu(Gly)2-","HCu(Gly)+2"))
-loga_Cu <- -4
-loga_Gly <- -1
-diagram(afun(loga_Cu,loga_Gly),fill=NULL,col="blue",
- names=species()$name,col.names="blue",add=TRUE)
-water.lines()
-# the glycine ionization constants could be calculated using
-# subcrt, here they are taken from A&D Table II
-abline(v=c(2.35,9.778),lty=3)
-# now put glycinium (low pH) in the basis
-basis(c("Cu+2","H2O","H+","e-","glycinium","CO2"),c(999,0,999,999,-2,999))
-species(c("copper","cuprite","tenorite","Cu+","Cu+2","CuO2-2","HCuO2-",
- "Cu(Gly)+","Cu(Gly)2","Cu(Gly)2-","HCu(Gly)+2"))
-diagram(afun(loga_Cu,loga_Gly),fill=NULL,col="green",
- names=NULL,col.names="green",add=TRUE)
-# let's forget our changes to 'thermo' so that examples
-# below that use glycine will work as expected
-# data(thermo)
Deleted: pkg/CHNOSZ/demo/pie.R
===================================================================
--- pkg/CHNOSZ/demo/pie.R 2012-11-27 14:04:22 UTC (rev 30)
+++ pkg/CHNOSZ/demo/pie.R 2012-12-09 14:09:44 UTC (rev 31)
@@ -1,144 +0,0 @@
-# This is an attempt to predict some characteristics of the relative
-# abundances of organisms that are depicted in the pie charts shown by
-# Spear et al., 2005 (Spear, J. R., Walker, J. J., McCollom, T. M. and
-# Pace, N. R. Hydrogen and bioenergetics in the Yellowstone geothermal
-# ecosystem. Proc. Natl. Acad. Sci. U. S. A., 102, 2555-2560, 2005.).
-# For each type of organism present, we use a single model protein. We
-# take the reported relative abundances of the four most abundant
-# organisms to generate an assemblage of proteins that buffers the
-# activies of CO2, H2O and NH3. Then we calculate relative abundances of
-# all proteins at different values of oxygen fugacity and H2S activity
-# and display them on pie charts that can be compared with organismal
-# abundances. The results show that it is possible to predict conditions
-# that foster high (many slices in the pie diagrams) or low (few slices)
-# diversity, even if the make up of the community is not perfectly
-# replicated. jmd 2008-2009
-### Setup
-opar <- par(no.readonly=TRUE)
-# Bacterial species and abundances from Spear et al., 2005 (Fig. 4):
-O1 <- c('Aquificales','Thermotogales','Bacillus','Thermodesulfobacteria')
-O2 <- c('Thermus/Deinococcus','Hydrogenobacter', 'Hydrogenobaculum','Bacteroidetes')
-O3 <- c('a-proteobacterium','d-proteobacterium','g-proteobacterium','OD-1')
-ORGANISMS <- c(O1,O2,O3)
-# Which species are most abundant in low and high sulfur environments;
-# the first three of each are aquificales, thermotogales, thermodesulfobacteria:
-lowS.which <- c(1,2,4,3,5)
-lowS.values <- c(75,12,4,7,2)
-highS.which <- c(1,2,4,11,6,9,7,8,10,13)
-highS.values <- c(63,2,10,2,5,2,9,2,3,1)
-# What chemical species we use to represent each organism:
-P1 <- c('ACCA_AQUAE','A8F7V4_THELT','RBL1_THIDA')
-P2 <- c('Q93EV7_9BACT','ACCA_DEIRA','Q05KD0_HYDTH','A7WGI1_9AQUI')
-P3 <- c('Q64VW6_BACFR','ACCA_CAUCR','A1VC70_DESVV','ACCA_PSEAE')
-PROTEINS <- c(P1,P2,P3)
-### Function definitions
-# To make pie charts for calculated abundances:
-plot.pie.calc <- function(which="high",T=25,main='') {
- # first clean up the buffer definition in case we have
- # been run already
- thermo$buffers <<- thermo$buffers[thermo$buffers$name!='PROTEINS',]
- # we take four predominant proteins from SWM+05
- myprot <- PROTEINS[get(paste(which,"S.which",sep=""))][1:4]
- mypercent <- get(paste(which,"S.values",sep=""))[1:4]
- # use these four proteins to create a buffer
- mybufprot <- paste(myprot,'RESIDUE',sep='.')
- mod.buffer('PROTEINS',mybufprot,'aq',log10(mypercent/100))
- # our species are the residues of all proteins
- species(delete=TRUE)
- # 20120520 the next 3 lines take the place of previous protein.residue call
- aa <- ip2aa(PROTEINS, residue=TRUE)
- aa$organism <- paste(aa$organism, "RESIDUE", sep=".")
- add.protein(aa)
- species(paste(aa$protein, aa$organism, sep="_"))
- # assign the buffer to three basis species
- basis(c('CO2','H2O','NH3'),'PROTEINS')
- # calculate the buffered activities
- a <- affinity(return.buffer=TRUE,balance=1,T=T)
- # make the titles
- sub <- c2s(paste(names(a)[1:3],round(as.numeric(a)[1:3])))
- main <- paste('\n',main,'\n',sub)
- # set the total species activities to those in the buffer
- species(1:nrow(species()),-99)
- species(mybufprot,log10(mypercent/100))
- # get the activities back to numeric
- basis(names(a)[1:3],as.numeric(a)[1:3])
- thermo$basis$logact <<- as.numeric(thermo$basis$logact)
- # colors
- col <- rep('white',99)
- col[match(myprot,PROTEINS)] <- heat.colors(4)
- # calculate the distribution of species
- mylogaH2O <- thermo$basis$logact[rownames(thermo$basis)=='H2O']
- a <- affinity(H2O=c(mylogaH2O,mylogaH2O-1,2),T=T)
- logacts <- equilibrate(a, normalize=TRUE, balance=1)$loga.equil
- # assemble the names and logarithms of activities
- # of species that are above a minimum value
- names <- character()
- values <- numeric()
- cols <- character()
- logactmin <- -2
- for(i in 1:length(logacts)) {
- myvalue <- logacts[[i]][1]
- if(myvalue > logactmin) {
- names <- c(names,ORGANISMS[i])
- values <- c(values,myvalue)
- cols <- c(cols,col[i])
- }
- }
- # remove the logarithms
- values <- 10^values
- # sort them by abundance
- isort <- sort(values,index.return=TRUE,decreasing=TRUE)$ix
- names <- names[isort]
- values <- values[isort]
- cols <- cols[isort]
- # make a pie chart
- pie(values,names,clockwise=FALSE,main=main,col=cols,radius=0.7)
-}
-# To plot pie charts for observed abundances of organisms:
-plot.pie.obs <- function(which="low") {
- # the values from SWM+05
- names <- ORGANISMS[get(paste(which,"S.which",sep=""))]
- values <- get(paste(which,"S.values",sep=""))
- main <- paste("observed at",which,"H2S")
- # colors for the four dominant species
- mycol <- heat.colors(4)
- # colors for the rest
- mycol <- c(mycol,rep('white',length(names)-length(mycol)))
- # sort the values
- isort <- sort(values,index.return=TRUE,decreasing=TRUE)$ix
- values <- values[isort]
- names <- names[isort]
- mycol <- mycol[isort]
- pie(values,names,clockwise=FALSE,main=main,col=mycol,radius=0.7)
-}
-# To plot both types of pie diagrams (showing calculated protein
-# activities and observed abundances of organisms) at a given temperature:
-plot.pie <- function(T=80) {
- # first deal with the layout
- layout(matrix(1:4,byrow=TRUE,nrow=2))
- opar <- par(mar=c(1,1,1,1))
- # basis definition
- basis(c('CO2','H2O','NH3','H2S','hydrogen','H+'))
- basis('pH',7)
- val <- function(text)
- round(as.numeric(thermo$basis$logact[rownames(thermo$basis)==text]))
- # now to plotting
- # low sulfur and relatively oxidizing
- basis(c('H2','H2S'),c(-9,-9))
- plot.pie.calc("low",T=T,main=paste("calculated for H2",val("H2"),"H2S",val("H2S")))
- label.plot('a')
- plot.pie.obs("low")
- label.plot('b')
- # high sulfur and relatively reducing
- basis(c('H2','H2S'),c(-5,-3))
- plot.pie.calc("high",T=T,main=paste("calculated for H2",val("H2"),"H2S",val("H2S")))
- label.plot('c')
- plot.pie.obs("high")
- label.plot('d')
- par(opar)
-}
-### Now run it!
-# at 80 degrees C:
-plot.pie(80)
-# reset plot device
-par(opar)
Modified: pkg/CHNOSZ/inst/NEWS
===================================================================
--- pkg/CHNOSZ/inst/NEWS 2012-11-27 14:04:22 UTC (rev 30)
+++ pkg/CHNOSZ/inst/NEWS 2012-12-09 14:09:44 UTC (rev 31)
@@ -1,4 +1,4 @@
-CHANGES IN CHNOSZ 0.9.8-5 (2012-11-27)
+CHANGES IN CHNOSZ 0.9.8-6 (2012-12-09)
--------------------------------------
SIGNIFICANT USER-VISIBLE CHANGES:
Modified: pkg/CHNOSZ/inst/tests/test-affinity.R
===================================================================
--- pkg/CHNOSZ/inst/tests/test-affinity.R 2012-11-27 14:04:22 UTC (rev 30)
+++ pkg/CHNOSZ/inst/tests/test-affinity.R 2012-12-09 14:09:44 UTC (rev 31)
@@ -93,6 +93,8 @@
pH <- c(7.350, 7.678, 7.933, 7.995, 8.257)
# Eq. 24 of the paper
H2 <- -11+T*3/40
+ # remove "RESIDUE" entries in thermo$obigt (clutter from first test)
+ data(thermo)
basis(c("HCO3-", "H2O", "NH3", "HS-", "H2", "H+"),
"aq", c(-3, 0, -4, -7, 999, 999))
sites <- c("N", "S", "R", "Q", "P")
@@ -133,3 +135,12 @@
species("CSG_HALJP")
expect_equal(affinity()$values[[1]][1], A.2303RT.ionized, 1e-6)
})
+
+test_that("affinity() for proteins keeps track of pH on 2-D calculations", {
+ # (relates to the "thisperm" construction in A.ionization() )
+ basis("CHNOS+")
+ species("LYSC_CHICK")
+ a1 <- affinity(pH=c(6, 8, 3))
+ a2 <- affinity(pH=c(6, 8, 3), T=c(0, 75, 4))
+ expect_equal(as.numeric(a1$values[[1]]), a2$values[[1]][, 2])
+})
Modified: pkg/CHNOSZ/inst/tests/test-species.R
===================================================================
--- pkg/CHNOSZ/inst/tests/test-species.R 2012-11-27 14:04:22 UTC (rev 30)
+++ pkg/CHNOSZ/inst/tests/test-species.R 2012-12-09 14:09:44 UTC (rev 31)
@@ -26,4 +26,33 @@
species("H2O")
expect_warning(species("CO2", delete=TRUE), "not present, so can not be deleted")
expect_is(species("water", delete=TRUE), "NULL")
+ # we should also get NULL if *all* species are deleted
+ species("H2O")
+ expect_is(species(delete=TRUE), "NULL")
})
+
+test_that("non-available species cause error, and species can be added or modified by numeric index", {
+ basis("CHNOS")
+ expect_error(species("wate"), "species not available")
+ # add CO2, aq
+ species("CO2")
+ # we can't add the same species twice
+ expect_equal(nrow(species("CO2")), 1)
+ # change it to gas
+ expect_equal(species(1, "gas")$state, "gas")
+ # change its log fugacity to -5
+ expect_equal(species(1, -5)$logact, -5)
+ # add CO2, aq
+ expect_equal(nrow(species("CO2")), 2)
+ # add alanine by index in thermo$obigt
+ expect_equal(nrow(species(info("alanine"))), 3)
+})
+
+test_that("index_return provides indices for touched species", {
+ basis("CHNOS")
+ expect_equal(species("CO2", index.return=TRUE), 1)
+ # here it's "touched" (but not added or modified)
+ expect_equal(species("CO2", index.return=TRUE), 1)
+ expect_equal(species(c("H2O", "NH3"), index.return=TRUE), c(2, 3))
+ expect_equal(species(2, "gas", index.return=TRUE), 2)
+})
Modified: pkg/CHNOSZ/inst/tests/test-util.affinity.R
===================================================================
--- pkg/CHNOSZ/inst/tests/test-util.affinity.R 2012-11-27 14:04:22 UTC (rev 30)
+++ pkg/CHNOSZ/inst/tests/test-util.affinity.R 2012-12-09 14:09:44 UTC (rev 31)
@@ -11,7 +11,7 @@
expect_equal(dim(A[[1]]), c(3, 2, 4))
})
-test_that("A.ionization() handles complex arguments generated by energy.args()", {
+test_that("A.ionization() handles arguments generated by energy.args()", {
# specifically, when the user asks for Eh, energy.args() produces
# values of pe that are specified in all dimensions
# (because the conversion from Eh is a function of temperature)
@@ -19,9 +19,9 @@
# we need a basis() for the function to poke around in
# (no e- or H+ needed at this time though)
basis("CHNOS")
- ea <- energy.args(list(CO2=c(-5, 0, 6), pH=c(0, 14, 2), Eh=c(-1, 0, 3), T=c(25, 26, 2)))
+ ea <- energy.args(list(CO2=c(-5, 0, 6), pH=c(0, 14, 2), Eh=c(-1, 0, 3), T=c(25, 28, 4)))
A <- A.ionization(1, ea$vars, ea$vals)
- expect_equal(dim(A[[1]]), c(6, 2, 3, 2))
+ expect_equal(dim(A[[1]]), c(6, 2, 3, 4))
# a simpler set of arguments, for testing dimensions with identical extents
# (the reason for the "for(dim in c(TPpHdim, otherdim))" in A.ionization())
ea <- energy.args(list(pH=c(0, 14, 2), Eh=c(-1, 0, 2)))
Modified: pkg/CHNOSZ/man/examples.Rd
===================================================================
--- pkg/CHNOSZ/man/examples.Rd 2012-11-27 14:04:22 UTC (rev 30)
+++ pkg/CHNOSZ/man/examples.Rd 2012-12-09 14:09:44 UTC (rev 31)
@@ -13,9 +13,9 @@
\usage{
examples(do.png = FALSE)
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/chnosz -r 31
More information about the CHNOSZ-commits
mailing list