[CHNOSZ-commits] r72 - in pkg/CHNOSZ: . R inst man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sat Dec 20 08:36:08 CET 2014
Author: jedick
Date: 2014-12-20 08:36:08 +0100 (Sat, 20 Dec 2014)
New Revision: 72
Added:
pkg/CHNOSZ/R/mosaic.R
pkg/CHNOSZ/man/mosaic.Rd
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/R/affinity.R
pkg/CHNOSZ/inst/NEWS
pkg/CHNOSZ/man/affinity.Rd
Log:
add mosaic() for affinity calculations with changing basis species
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2014-12-15 08:45:44 UTC (rev 71)
+++ pkg/CHNOSZ/DESCRIPTION 2014-12-20 07:36:08 UTC (rev 72)
@@ -1,6 +1,6 @@
-Date: 2014-12-15
+Date: 2014-12-20
Package: CHNOSZ
-Version: 1.0.3-9
+Version: 1.0.3-10
Title: Chemical Thermodynamics and Activity Diagrams
Author: Jeffrey Dick
Maintainer: Jeffrey Dick <j3ffdick at gmail.com>
Modified: pkg/CHNOSZ/R/affinity.R
===================================================================
--- pkg/CHNOSZ/R/affinity.R 2014-12-15 08:45:44 UTC (rev 71)
+++ pkg/CHNOSZ/R/affinity.R 2014-12-20 07:36:08 UTC (rev 72)
@@ -26,6 +26,11 @@
mybasis <- thermo$basis
myspecies <- thermo$species
+ # stop if Eh or pe is requested but e- isn't in the basis
+ if(any(c("Eh", "pe") %in% names(args$lims))) {
+ if(!"e-" %in% rownames(mybasis)) stop("variable Eh or pe requested but e- isn't in the basis")
+ }
+
if(!is.null(property)) {
# the user just wants an energy property
buffer <- FALSE
Added: pkg/CHNOSZ/R/mosaic.R
===================================================================
--- pkg/CHNOSZ/R/mosaic.R (rev 0)
+++ pkg/CHNOSZ/R/mosaic.R 2014-12-20 07:36:08 UTC (rev 72)
@@ -0,0 +1,66 @@
+# CHNOSZ/mosaic.R
+# calculate affinities with changing basis species
+# 20141220 jmd
+
+# function to calculate affinities with mosaic of basis species
+mosaic <- function(bases, blend=FALSE, ...) {
+ # the arguments for affinity()
+ myargs <- list(...)
+ # are the swapped basis species on the plot?
+ # (the first one should be present in the starting basis set)
+ iswap <- match(bases[1], names(myargs))
+ # the log activity of the starting basis species
+ logact.swap <- basis()$logact[match(bases[1], row.names(basis()))]
+ # a list where we'll keep the affinity calculations
+ affs <- list()
+ for(i in seq_along(bases)) {
+ # set up argument list: name of swapped-in basis species
+ if(!is.na(iswap)) names(myargs)[iswap] <- bases[i]
+ # calculate affinities
+ affs[[i]] <- do.call(affinity, myargs)
+ # change the basis species; restore the original at the end of the loop
+ if(i < length(bases)) {
+ swap.basis(bases[i], bases[i+1])
+ basis(bases[i+1], logact.swap)
+ } else {
+ swap.basis(bases[i], bases[1])
+ basis(bases[1], logact.swap)
+ names <- row.names(basis())
+ }
+ }
+ # calculate affinities of formation of basis species
+ ispecies <- species()$ispecies
+ species.logact <- species()$logact
+ species(delete=TRUE)
+ species(bases)
+ A.basis <- do.call(affinity, myargs)
+ # restore original species with original activities
+ species(delete=TRUE)
+ species(ispecies, species.logact)
+ # affinities calculated using the first basis species
+ A.species <- affs[[1]]
+ if(blend) {
+ # calculate affinities using relative abundances of basis species
+ e <- equilibrate(A.basis)
+ for(j in seq_along(affs)) {
+ for(i in seq_along(A.species$values)) {
+ # start with zero affinity
+ if(j==1) A.species$values[[i]][] <- 0
+ # add affinity scaled by relative abundance of this basis species
+ A.species$values[[i]] <- A.species$values[[i]] + affs[[j]]$values[[i]] * 10^e$loga.equil[[j]]
+ }
+ }
+ } else {
+ # use affinities from the single predominant basis species
+ d <- diagram(A.basis, plot.it=FALSE)
+ # merge affinities using the second, third, ... basis species
+ for(j in tail(seq_along(affs), -1)) {
+ is.predominant <- d$predominant==j
+ for(i in seq_along(A.species$values)) {
+ A.species$values[[i]][is.predominant] <- affs[[j]]$values[[i]][is.predominant]
+ }
+ }
+ }
+ # return the affinities for the species and basis species
+ return(list(A.species=A.species, A.basis=A.basis))
+}
Modified: pkg/CHNOSZ/inst/NEWS
===================================================================
--- pkg/CHNOSZ/inst/NEWS 2014-12-15 08:45:44 UTC (rev 71)
+++ pkg/CHNOSZ/inst/NEWS 2014-12-20 07:36:08 UTC (rev 72)
@@ -1,5 +1,5 @@
-CHANGES IN CHNOSZ 1.0.3-9 (2014-12-15)
---------------------------------------
+CHANGES IN CHNOSZ 1.0.3-10 (2014-12-20)
+---------------------------------------
- Add files with average amino acid compositions of proteins from Bison
Pool grouped according to annotation keyword (DS11.csv) (moved here
@@ -24,6 +24,9 @@
- water.lines() gets 'O2state' argument for state of O2.
+- Add mosaic() function for affinity calculations with changing basis
+ species.
+
CHANGES IN CHNOSZ 1.0.3 (2014-01-12)
------------------------------------
Modified: pkg/CHNOSZ/man/affinity.Rd
===================================================================
--- pkg/CHNOSZ/man/affinity.Rd 2014-12-15 08:45:44 UTC (rev 71)
+++ pkg/CHNOSZ/man/affinity.Rd 2014-12-20 07:36:08 UTC (rev 72)
@@ -27,20 +27,31 @@
The calculation of chemical affinities relies on the current definitions of the \code{\link{basis}} species and \code{\link{species}} of interest.
It is possible to use the results of \code{affinity} to generate equilibrium activity diagrams using \code{\link{diagram}}.
- The equation used to calculate chemical affinity \emph{\bold{A}} is \emph{\bold{A}}=\eqn{RT\ln (K/Q)}{RT*ln(K/Q)} (Kondepudi and Prigogine, 1998), where \eqn{K} denotes the equilibrium constant of the reaction and \eqn{Q} stands for the activity product of the species in the reaction. (The equilibrium constant is related to standard Gibbs energy of reaction, \eqn{\Delta G^{\circ}_r}{DeltaG0r}, by \eqn{\Delta G^{\circ}_r = -2.303RT\log K}{DeltaG0r = -2.303*RT*logK}, where \eqn{R} and \eqn{T} stand for, respectively, the gas constant and temperature). With the approach of a given reaction to a state of equilibrium, the chemical affinity tends toward zero, or \eqn{K = Q}.
+The equation used to calculate chemical affinity \emph{\bold{A}} is \emph{\bold{A}}=\eqn{RT\ln (K/Q)}{RT*ln(K/Q)} (Kondepudi and Prigogine, 1998), where \eqn{K} denotes the equilibrium constant of the reaction and \eqn{Q} stands for the activity product of the species in the reaction.
+(The equilibrium constant is related to standard Gibbs energy of reaction, \eqn{\Delta G^{\circ}_r}{DeltaG0r}, by \eqn{\Delta G^{\circ}_r = -2.303RT\log K}{DeltaG0r = -2.303*RT*logK}, where \eqn{R} and \eqn{T} stand for, respectively, the gas constant and temperature).
+With the approach of a given reaction to a state of equilibrium, the chemical affinity tends toward zero, or \eqn{K = Q}.
- Valid properties are \samp{A} or NULL for chemical affinity, \samp{logK} or \samp{logQ} for logarithm of equilibrium constant and reaction activity product, or any of the properties available in \code{\link{subcrt}} except for \samp{rho}. The properties returned are those of the formation reactions of the species of interest from the basis species. It is also possible to calculate the properties of the species of interest themselves (not their formation reactions) by setting the \code{property} to \samp{G.species}, \samp{Cp.species}, etc. Except for \samp{A}, the properties of proteins or their reactions calculated in this manner are restricted to nonionized proteins.
+Valid properties are \samp{A} or NULL for chemical affinity, \samp{logK} or \samp{logQ} for logarithm of equilibrium constant and reaction activity product, or any of the properties available in \code{\link{subcrt}} except for \samp{rho}.
+The properties returned are those of the formation reactions of the species of interest from the basis species.
+It is also possible to calculate the properties of the species of interest themselves (not their formation reactions) by setting the \code{property} to \samp{G.species}, \samp{Cp.species}, etc.
+Except for \samp{A}, the properties of proteins or their reactions calculated in this manner are restricted to nonionized proteins.
- Zero, one, or more leading arguments to the function identify which of the chemical activities of basis species, temperature, pressure and/or ionic strength to vary. The names of each of these arguments may be the formula of any of the basis species of the system, or \samp{T}, \samp{P}, \samp{pe}, \samp{pH}, \samp{Eh}, or \samp{IS} (but names may not be repeated). To use the names of charged basis species such as \samp{K+} and \samp{SO4-2} as the arguments, they should be enclosed in quotes (see the example for aluminum speciation in \code{\link{diagram}}). The value of each argument is of the form \code{c(min,max)} or \code{c(min,max,res)} where \code{min} and \code{max} refer to the minimimum and maximum values of variable identified by the name of the argument, and \code{res} denotes the resolution, or number of points along which to do the calculations (this is assigned a default value of 128 if it is missing). For any arguments that refer to basis species, the numerical values are the logarithms of the activities of that basis species, or logarithms of fugacities if it is a gas. Unlike the \code{energy} function, the units of \samp{T} and \samp{P} in \code{affinity} are those the user has set using \code{\link{T.units}} and \code{\link{P.units} }(on program start-up these are \eqn{^{\circ}}{°}C and bar, respectively).
+Zero, one, or more leading arguments to the function identify which of the chemical activities of basis species, temperature, pressure and/or ionic strength to vary.
+The names of each of these arguments may be the formula of any of the basis species of the system, or \samp{T}, \samp{P}, \samp{pe}, \samp{pH}, \samp{Eh}, or \samp{IS} (but names may not be repeated).
+To use the names of charged basis species such as \samp{K+} and \samp{SO4-2} as the arguments, they should be enclosed in quotes (see the example for aluminum speciation in \code{\link{diagram}}).
+The value of each argument is of the form \code{c(min,max)} or \code{c(min,max,res)} where \code{min} and \code{max} refer to the minimimum and maximum values of variable identified by the name of the argument, and \code{res} denotes the resolution, or number of points along which to do the calculations (this is assigned a default value of 128 if it is missing).
+For any arguments that refer to basis species, the numerical values are the logarithms of the activities of that basis species, or logarithms of fugacities if it is a gas.
+Unlike the \code{energy} function, the units of \samp{T} and \samp{P} in \code{affinity} are those the user has set using \code{\link{T.units}} and \code{\link{P.units} }(on program start-up these are \eqn{^{\circ}}{°}C and bar, respectively).
- If one or more buffers are assigned to the definition of \code{\link{basis}} species, \code{affinity} calls \code{\link{buffer}} to calculate the logarithms of activities of these basis species from the buffer.
+If one or more buffers are assigned to the definition of \code{\link{basis}} species, \code{affinity} calls \code{\link{buffer}} to calculate the logarithms of activities of these basis species from the buffer.
The \code{iprotein} and \code{loga.protein} arguments can be used to compute the chemical affinities of formation reactions of proteins that are not in the current \code{\link{species}} definition.
This approach can be utilized in order to calculate the properties of many proteins in a fraction of the time it would take to calculate them individually.
The appropriate \code{\link{basis}} species still must be defined prior to calling \code{affinity}.
\code{iprotein} contains indices of desired proteins in \code{\link{thermo}$protein}; \code{affinity} adds to the species list the amino acid residues and and terminal H2O group (indicated by \dQuote{RESIDUE} in \code{thermo$protein}) then calculates the properties of the reactions for the residues and terminal group, including ionization effects, and adds them together to get those of the indicated proteins.
- In CHNOSZ version 0.9, \code{energy} gained a new argument \samp{transect} which is set to TRUE by \code{energy.args} when the length(s) of the variables is(are) greater than three. In this mode of operation, instead of performing the calculations on an \eqn{n}{n}-dimensional grid, the affinities are calculated on an \eqn{n}{n}-dimensional transect through chemical potential (possibly including T and/or P) space.
+In CHNOSZ version 0.9, \code{energy} gained a new argument \samp{transect} which is set to TRUE by \code{energy.args} when the length(s) of the variables is(are) greater than three.
+In this mode of operation, instead of performing the calculations on an \eqn{n}{n}-dimensional grid, the affinities are calculated on an \eqn{n}{n}-dimensional transect through chemical potential (possibly including T and/or P) space.
}
Added: pkg/CHNOSZ/man/mosaic.Rd
===================================================================
--- pkg/CHNOSZ/man/mosaic.Rd (rev 0)
+++ pkg/CHNOSZ/man/mosaic.Rd 2014-12-20 07:36:08 UTC (rev 72)
@@ -0,0 +1,71 @@
+\encoding{UTF-8}
+\name{mosaic}
+\alias{mosaic}
+\title{Chemical Affinities with Changing Basis Species}
+\description{
+Calculate chemical affinities of formation reactions of species using basis species that change with the conditions.
+}
+
+\usage{
+ mosaic(bases, blend=FALSE, ...)
+}
+
+\arguments{
+ \item{bases}{character, basis species to include in the calculation}
+ \item{blend}{logical, use relative abundances of basis species?}
+ \item{...}{additional arguments to be passed to \code{\link{affinity}}}
+}
+
+\details{
+
+\code{mosaic} can be used to calculate the reaction affinities when the basis species listed in \code{bases} change in relative abundance over the range of conditions, due to e.g. ionization, complexation or redox reactions.
+This is a way to \dQuote{speciate the basis species}.
+An example is consideration of the speciation of sulfur (\samp{SO4-2}, \samp{HSO4-}, \samp{HS-} and \samp{H2S}) as a function of oxygen fugacity (or Eh) and pH in calculating the relative stabilities of minerals and aqueous species in the system Fe-S-O-H.
+The first species listed in \code{bases} should correspond to one of the basis species currently defined, and the function only supports calculations using basis species that all derive from that species (i.e. all share the same element).
+The arguments in \code{...} are passed to \code{affinity} to specify the conditions.
+
+The function calculates the affinities using each basis species in turn, changing them via \code{\link{swap.basis}}.
+If \code{blend} is FALSE (the default), the function returns the affinities calculated using the single predominant basis species in \code{bases} at each condition.
+Chemical activity or predominance diagrams constructed using this method have been described as \dQuote{mosaic diagrams} in the literature.
+If \code{blend} is TRUE, the function combines the affinities of the formation reactions weighted by the relative abundances of the basis species at each condition.
+This tends to produce curved boundaries.
+
+}
+
+\value{
+A list containing \code{A.species} (affinities of formation of the species with changing basis speceis) and \code{A.basis} (affinities of formation of the basis species in terms of the first basis species), each having same structure as the list returned by \code{\link{affinity}}.
+}
+
+\examples{
+\dontshow{data(thermo)}# Fe-minerals and aqueous species in Fe-S-O-H system
+# speciate SO4-2, HSO4-, HS-, H2S as a function of Eh and pH
+# after Garrels and Christ, 1965 Figure 7.20
+pH <- c(0, 14, 200)
+Eh <- c(-1, 1, 200)
+T <- 25
+basis(c("FeO", "SO4-2", "H2O", "H+", "e-"))
+basis("SO4-2", -6)
+species(c("Fe+2", "Fe+3"), -6)
+species(c("pyrrhotite", "pyrite", "hematite", "magnetite"))
+# the basis species we'll swap through
+bases <- c("SO4-2", "HSO4-", "HS-", "H2S")
+# calculate affinities using the predominant basis species
+m1 <- mosaic(bases, pH=pH, Eh=Eh, T=T)
+# make a diagram and add water stability lines
+diagram(m1$A.species)
+water.lines("pH", "Eh", T=convert(T, "K"), col="seagreen", lwd=1.5)
+# show lines for Fe(aq) = 10^-4 M
+species(c("Fe+2", "Fe+3"), -4)
+m2 <- mosaic(bases, pH=pH, Eh=Eh, T=T)
+diagram(m2$A.species, add=TRUE, names=NULL, dotted=3)
+title(main=paste("Iron oxides and sulfides in water, log(total S) = -6",
+ "After Garrels and Christ, 1965", sep="\n"))
+# we could overlay the basis species predominance fields
+#diagram(m1$A.basis, add=TRUE, col="blue", col.names="blue", dotted=3)
+}
+
+\references{
+ Garrels, R. M. and Christ, C. L. (1965) \emph{Solutions, Minerals, and Equilibria}, Harper & Row, New York, 450 p. \url{http://www.worldcat.org/oclc/517586}
+}
+
+\keyword{secondary}
More information about the CHNOSZ-commits
mailing list