[CHNOSZ-commits] r823 - in pkg/CHNOSZ: . R vignettes
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Jan 11 13:10:41 CET 2024
Author: jedick
Date: 2024-01-11 13:10:40 +0100 (Thu, 11 Jan 2024)
New Revision: 823
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/R/affinity.R
pkg/CHNOSZ/R/buffer.R
pkg/CHNOSZ/vignettes/anintro.Rmd
Log:
Use more descriptive variable names in affinity()
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2024-01-09 02:23:52 UTC (rev 822)
+++ pkg/CHNOSZ/DESCRIPTION 2024-01-11 12:10:40 UTC (rev 823)
@@ -1,6 +1,6 @@
-Date: 2024-01-06
+Date: 2024-01-10
Package: CHNOSZ
-Version: 2.0.0-42
+Version: 2.0.0-43
Title: Thermodynamic Calculations and Diagrams for Geochemistry
Authors at R: c(
person("Jeffrey", "Dick", , "j3ffdick at gmail.com", role = c("aut", "cre"),
Modified: pkg/CHNOSZ/R/affinity.R
===================================================================
--- pkg/CHNOSZ/R/affinity.R 2024-01-09 02:23:52 UTC (rev 822)
+++ pkg/CHNOSZ/R/affinity.R 2024-01-11 12:10:40 UTC (rev 823)
@@ -54,18 +54,19 @@
args <- energy.args(args.orig, transect = transect)
args <- c(args, list(sout = sout, exceed.Ttr = exceed.Ttr, exceed.rhomin = exceed.rhomin))
- # The species we're given
+ # The user-defined species (including basis species, formed species, and proteins)
thermo <- get("thermo", CHNOSZ)
- mybasis <- thermo$basis
myspecies <- thermo$species
if(!is.null(property)) {
+
# The user just wants an energy property
buffer <- FALSE
args$what <- property
- out <- do.call("energy", args)
- a <- out$a
- sout <- out$sout
+ energy_result <- do.call("energy", args)
+ affinity_values <- energy_result$a
+ energy_sout <- energy_result$sout
+
} else {
# Affinity calculations
@@ -80,7 +81,7 @@
# Add protein residues to the species list
resnames <- c("H2O", aminoacids(3))
# Residue activities set to zero; account for protein activities later
- resprot <- paste(resnames,"RESIDUE", sep = "_")
+ resprot <- paste(resnames, "RESIDUE", sep = "_")
species(resprot, 0)
# Re-read thermo because the preceding command changed the species
thermo <- get("thermo", CHNOSZ)
@@ -89,11 +90,9 @@
# Buffer stuff
buffer <- FALSE
- hasbuffer <- !can.be.numeric(mybasis$logact)
- # However, variables specified in the arguments take precedence over the buffer 20171013
- isnotvar <- !rownames(mybasis) %in% args$vars
- ibufbasis <- which(hasbuffer & isnotvar)
- if(!is.null(mybasis) & length(ibufbasis) > 0) {
+ # The buffered basis species are those that have non-numeric logact and are not listed in the arguments
+ which.basis.is.buffered <- which(!can.be.numeric(thermo$basis$logact) & !rownames(thermo$basis) %in% args$vars)
+ if(!is.null(thermo$basis) & length(which.basis.is.buffered) > 0) {
buffer <- TRUE
message('affinity: loading buffer species')
if(!is.null(thermo$species)) is.species <- 1:nrow(thermo$species) else is.species <- numeric()
@@ -107,9 +106,9 @@
}
# Here we call 'energy'
- aa <- do.call("energy", args)
- a <- aa$a
- sout <- aa$sout
+ energy_result <- do.call("energy", args)
+ affinity_values <- energy_result$a
+ energy_sout <- energy_result$sout
if(return.sout) return(sout)
@@ -116,13 +115,13 @@
# More buffer stuff
if(buffer) {
args$what <- "logact.basis"
- args$sout <- sout
+ args$sout <- energy_sout
logact.basis.new <- logact.basis <- do.call("energy", args)$a
ibasis.new <- numeric()
for(k in 1:length(buffers)) {
- ibasis <- which(as.character(mybasis$logact) == buffers[k])
+ ibasis <- which(as.character(thermo$basis$logact) == buffers[k])
# Calculate the logKs from the affinities
- logK <- a
+ logK <- affinity_values
for(i in 1:length(logK)) {
logK[[i]] <- logK[[i]] + thermo$species$logact[i]
for(j in 1:length(logact.basis.new)) {
@@ -129,30 +128,26 @@
logK[[i]] <- logK[[i]] - logact.basis.new[[j]] * thermo$species[i, j]
}
}
- lbn <- buffer(logK = logK, ibasis = ibasis, logact.basis = logact.basis.new,
- is.buffer = is.buffer[[k]], balance = balance)
- for(j in 1:length(logact.basis.new)) if(j %in% ibasis) logact.basis.new[[j]] <- lbn[[2]][[j]]
+ buffer_result <- buffer(logK = logK, ibasis = ibasis, logact.basis = logact.basis.new, is.buffer = is.buffer[[k]], balance = balance)
+ for(j in 1:length(logact.basis.new)) if(j %in% ibasis) logact.basis.new[[j]] <- buffer_result[[2]][[j]]
# Calculation of the buffered activities' effect on chemical affinities
is.only.buffer.new <- is.only.buffer[is.only.buffer %in% is.buffer[[k]]]
- for(i in 1:length(a)) {
+ for(i in 1:length(affinity_values)) {
if(i %in% is.only.buffer.new) next
- for(j in 1:nrow(mybasis)) {
- # Let's only do this for the basis species specified by the user
- # even if others could be buffered
- if(!j %in% ibufbasis) next
+ for(j in 1:nrow(thermo$basis)) {
+ # Let's only do this for the basis species specified by the user even if others could be buffered
+ if(!j %in% which.basis.is.buffered) next
if(!j %in% ibasis) next
- aa <- a[[i]]
- a[[i]] <- aa + (logact.basis.new[[j]] - logact.basis[[j]]) * thermo$species[i, j]
- #if(!identical(a[[i]], aa)) print(paste(i, j))
+ affinity_values[[i]] <- affinity_values[[i]] + (logact.basis.new[[j]] - logact.basis[[j]]) * thermo$species[i, j]
}
}
if(k == length(buffers) & return.buffer) {
- logact.basis.new <- lbn[[2]]
- ibasis.new <- c(ibasis.new, lbn[[1]])
+ logact.basis.new <- buffer_result[[2]]
+ ibasis.new <- c(ibasis.new, buffer_result[[1]])
} else ibasis.new <- c(ibasis.new, ibasis)
}
species(is.only.buffer, delete = TRUE)
- if(length(is.only.buffer) > 0) a <- a[-is.only.buffer]
+ if(length(is.only.buffer) > 0) affinity_values <- affinity_values[-is.only.buffer]
# To return the activities of buffered basis species
tb <- logact.basis.new[unique(ibasis.new)]
if(!is.null(ncol(tb[[1]]))) {
@@ -177,7 +172,7 @@
loga.protein <- rep(loga.protein, length.out = length(iprotein))
protein.fun <- function(ip) {
tpext <- as.numeric(thermo$protein[iprotein[ip], 5:25])
- return(Reduce("+", pprod(a[ires], tpext)) - loga.protein[ip])
+ return(Reduce("+", pprod(affinity_values[ires], tpext)) - loga.protein[ip])
}
# Use another level of indexing to let the function report on its progress
jprotein <- 1:length(iprotein)
@@ -194,7 +189,7 @@
state <- resspecies$state[1]
name <- paste(thermo$protein$protein[iprotein], thermo$protein$organism[iprotein], sep = "_")
# The numbers of basis species in formation reactions of the proteins
- protbasis <- t(t((resspecies[ires, 1:nrow(mybasis)])) %*% t((thermo$protein[iprotein, 5:25])))
+ protbasis <- t(t((resspecies[ires, 1:nrow(thermo$basis)])) %*% t((thermo$protein[iprotein, 5:25])))
# Put them together
protspecies <- cbind(protbasis, data.frame(ispecies = ispecies, logact = loga.protein, state = state, name = name))
myspecies <- rbind(myspecies, protspecies)
@@ -201,8 +196,8 @@
rownames(myspecies) <- 1:nrow(myspecies)
## Update the affinity values
names(protein.affinity) <- ispecies
- a <- c(a, protein.affinity)
- a <- a[-ires]
+ affinity_values <- c(affinity_values, protein.affinity)
+ affinity_values <- affinity_values[-ires]
}
}
@@ -262,10 +257,10 @@
# Add IS value only if it given as an argument 20171101
# (even if its value is 0, the presence of IS will trigger diagram() to use "m" instead of "a" in axis labels)
iIS <- match("IS", names(args.orig))
- if(!is.na(iIS)) a <- list(fun = "affinity", args = allargs, sout = sout, property = property,
- basis = mybasis, species = myspecies, T = T, P = P, IS = args$IS, vars = vars, vals = vals, values = a)
- else a <- list(fun = "affinity", args = allargs, sout = sout, property = property,
- basis = mybasis, species = myspecies, T = T, P = P, vars = vars, vals = vals, values = a)
- if(buffer) a <- c(a, list(buffer = tb))
- return(a)
+ if(!is.na(iIS)) out <- list(fun = "affinity", args = allargs, sout = energy_sout, property = property,
+ basis = thermo$basis, species = myspecies, T = T, P = P, IS = args$IS, vars = vars, vals = vals, values = affinity_values)
+ else out <- list(fun = "affinity", args = allargs, sout = energy_sout, property = property,
+ basis = thermo$basis, species = myspecies, T = T, P = P, vars = vars, vals = vals, values = affinity_values)
+ if(buffer) out <- c(out, list(buffer = tb))
+ return(out)
}
Modified: pkg/CHNOSZ/R/buffer.R
===================================================================
--- pkg/CHNOSZ/R/buffer.R 2024-01-09 02:23:52 UTC (rev 822)
+++ pkg/CHNOSZ/R/buffer.R 2024-01-11 12:10:40 UTC (rev 823)
@@ -7,13 +7,13 @@
thermo <- get("thermo", CHNOSZ)
if(is.null(species)) {
iname <- which(name == thermo$buffer$name)
- if(length(iname)>0) species <- thermo$buffer$species[iname]
+ if(length(iname) > 0) species <- thermo$buffer$species[iname]
else species <- character()
}
ls <- length(species)
if(ls < length(name) | ls < length(state) | ls < length(logact))
stop('species must be at least as long as the other arguments')
- if(length(name) != ls) name <- rep(name,length.out = ls)
+ if(length(name) != ls) name <- rep(name, length.out = ls)
add <- TRUE
if(TRUE %in% (name %in% thermo$buffer$name)) {
add <- FALSE
@@ -20,10 +20,9 @@
imod <- which(thermo$buffer$name %in% name & thermo$buffer$species %in% species)
if(length(imod)>0) {
if(state[1] == '') {
- thermo$buffer <- thermo$buffer[-imod,]
+ thermo$buffer <- thermo$buffer[-imod, ]
assign("thermo", thermo, CHNOSZ)
- message(paste('mod.buffer: removed ',c2s(species),' in ',
- c2s(unique(name)),' buffer',sep = ''))
+ message(paste('mod.buffer: removed ', c2s(species), ' in ', c2s(unique(name)), ' buffer', sep = ''))
} else {
if(missing(state)) state <- thermo$buffer$state[imod]
if(missing(logact)) logact <- thermo$buffer$logact[imod]
@@ -34,12 +33,10 @@
thermo$buffer$state[imod] <- state
thermo$buffer$logact[imod] <- logact
assign("thermo", thermo, CHNOSZ)
- if(identical(state.old,state) & identical(logact.old,logact)) {
- message(paste('mod.buffer: nothing changed for ',
- c2s(species),' in ',c2s(unique(name)),' buffer',sep = ''))
+ if(identical(state.old, state) & identical(logact.old, logact)) {
+ message(paste('mod.buffer: nothing changed for ', c2s(species), ' in ', c2s(unique(name)), ' buffer', sep = ''))
} else {
- message(paste('mod.buffer: changed state and/or logact of ',
- c2s(species),' in ',c2s(unique(name)),' buffer',sep = ''))
+ message(paste('mod.buffer: changed state and/or logact of ', c2s(species), ' in ', c2s(unique(name)), ' buffer', sep = ''))
}
}
} else {
@@ -47,18 +44,18 @@
}
}
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$buffer <- rbind(thermo$buffer,t)
+ if(state[1] == '') state <- rep(thermo$opt$state, length.out = ls)
+ t <- data.frame(name = name, species = species, state = state, logact = logact)
+ thermo$buffer <- rbind(thermo$buffer, t)
assign("thermo", thermo, CHNOSZ)
message(paste('mod.buffer: added',c2s(unique(name))))
}
- return(invisible(thermo$buffer[thermo$buffer$name %in% name,]))
+ return(invisible(thermo$buffer[thermo$buffer$name %in% name, ]))
}
### Unexported functions ###
-buffer <- function(logK = NULL,ibasis = NULL,logact.basis = NULL,is.buffer = NULL,balance = 'PBB') {
+buffer <- function(logK = NULL, ibasis = NULL, logact.basis = NULL, is.buffer = NULL, balance = 'PBB') {
thermo <- get("thermo", CHNOSZ)
# If logK is NULL load the buffer species,
# otherwise perform buffer calculations
@@ -76,9 +73,9 @@
#ibuff <- info(species,state,quiet = TRUE)
ispecies <- c(ispecies, species(species, state, index.return = TRUE, add = TRUE))
}
- ispecies.new <- c(ispecies.new,list(ispecies))
+ ispecies.new <- c(ispecies.new, list(ispecies))
# Make sure to set the activities
- species(ispecies,thermo$buffer$logact[ib])
+ species(ispecies, thermo$buffer$logact[ib])
}
names(ispecies.new) <- buffers
return(ispecies.new)
@@ -91,17 +88,17 @@
bufbasis <- species.basis(thermo$species$ispecies[is.buffer])
bufname <- thermo$basis$logact[ibasis[1]]
basisnames <- rownames(thermo$basis)
- are.proteins <- grep('_',as.character(thermo$species$name[is.buffer]))
- if((length(are.proteins)>0 & balance == 'PBB') | balance == 1) {
+ are.proteins <- grep('_', as.character(thermo$species$name[is.buffer]))
+ if((length(are.proteins) > 0 & balance == 'PBB') | balance == 1) {
if(balance == 1) {
- basisnames <- c('product',basisnames)
- nb <- rep(1,nrow(bufbasis))
- bufbasis <- cbind(data.frame(product = nb),bufbasis)
+ basisnames <- c('product', basisnames)
+ nb <- rep(1, nrow(bufbasis))
+ bufbasis <- cbind(data.frame(product = nb), bufbasis)
} else {
- basisnames <- c('PBB',basisnames)
+ basisnames <- c('PBB', basisnames)
# Prepend a PBB column to bufbasis and increment ibasis by 1
nb <- as.numeric(protein.length(thermo$species$name[is.buffer]))
- bufbasis <- cbind(data.frame(PBB = nb),bufbasis)
+ bufbasis <- cbind(data.frame(PBB = nb), bufbasis)
}
ibasis <- ibasis + 1
# Make logact.basis long enough
@@ -108,7 +105,7 @@
tl <- length(logact.basis)
logact.basis[[tl+1]] <- logact.basis[[tl]]
# Rotate the entries so that the new one is first
- ilb <- c(tl+1,1:tl)
+ ilb <- c(tl+1, 1:tl)
logact.basis <- logact.basis[ilb]
}
# Say hello
@@ -119,7 +116,7 @@
if( (length(ibasis)+1) != length(is.buffer) & length(is.buffer) > 1) {
# Try to add buffered activities the user didn't specify
# (e.g. H2S in PPM buffer if only H2 was requested)
- for(i in 1:(length(is.buffer)-(length(ibasis)+1))) {
+ for(i in 1:(length(is.buffer) - (length(ibasis)+1))) {
newbasis <- NULL
# We want to avoid any basis species that might be used as the conservant
# look for additional activities to buffer ... do columns in reverse
@@ -131,10 +128,10 @@
}
}
if(!is.null(newbasis)) {
- ibasis <- c(ibasis,newbasis)
- ibasisadded <- c(ibasisadded,newbasis)
+ ibasis <- c(ibasis, newbasis)
+ ibasisadded <- c(ibasisadded, newbasis)
} else {
- stop('can not find enough buffered basis species for ',thermo$basis$logact[ibasis[1]],'.',sep = '')
+ stop('can not find enough buffered basis species for ', thermo$basis$logact[ibasis[1]], '.', sep = '')
}
}
}
@@ -146,26 +143,26 @@
# First try to get one that is present in all species
for(i in ncol(bufbasis):1) {
if(i %in% ibasis) next
- if(!TRUE %in% (bufbasis[,i] == 0)) newbasis <- i
+ if(!TRUE %in% (bufbasis[, i] == 0)) newbasis <- i
}
# Or look for one that is present at all
if(is.null(newbasis)) for(i in ncol(bufbasis):1) {
if(i %in% ibasis) next
- if(FALSE %in% (bufbasis[,i] == 0)) newbasis <- i
+ if(FALSE %in% (bufbasis[, i] == 0)) newbasis <- i
}
if(!is.null(newbasis)) {
ibasis <- c(ibasis,newbasis)
- #message(paste('buffer: the conserved activity is ',basisnames[newbasis],'.',sep = ''))
+ #message(paste('buffer: the conserved activity is ', basisnames[newbasis], '.', sep = ''))
#thermo$basis$logact[newbasis] <<- thermo$basis$logact[ibasis[1]]
}
else stop('no conserved activity found in your buffer (not enough basis species?)!')
}
- if(is.null(newbasis)) context <- '' else context <- paste(', ',basisnames[newbasis],' (conserved)',sep = '')
- reqtext <- paste(c2s(basisnames[ibasisrequested]),' (active)',sep = '')
- if(length(ibasisadded) == 0) addtext <- '' else addtext <- paste(', ',c2s(basisnames[ibasisadded]),sep = '')
+ if(is.null(newbasis)) context <- '' else context <- paste(', ', basisnames[newbasis], ' (conserved)', sep = '')
+ reqtext <- paste(c2s(basisnames[ibasisrequested]), ' (active)', sep = '')
+ if(length(ibasisadded) == 0) addtext <- '' else addtext <- paste(', ', c2s(basisnames[ibasisadded]), sep = '')
message(paste("buffer: '", bufname, "' for activity of ", reqtext, addtext, context, sep = ""))
# There could still be stuff here (over-defined system?)
- xx <- bufbasis[,-ibasis,drop = FALSE]
+ xx <- bufbasis[, -ibasis, drop = FALSE]
# For the case when all activities are buffered
if(ncol(xx) == 0) xx <- data.frame(xx = 0)
# Our stoichiometric matrix - should be square
@@ -177,12 +174,12 @@
b <- -logK[[is.buffer[i]]] + thermo$species$logact[is.buffer[i]]
if(ncol(xx) > 0) {
if(is.list(xx)) xxx <- xx[[1]] else xxx <- xx
- if(ncol(xx) == 1 & identical(as.numeric(xxx),0)) {
+ if(ncol(xx) == 1 & identical(as.numeric(xxx), 0)) {
# Do nothing
} else {
for(j in 1:ncol(xx)) {
- bs <- xx[i,j] * logact.basis[[match(colnames(xx)[j],basisnames)]]
- if(!is.matrix(bs)) bs <- matrix(bs,byrow = TRUE,nrow = nrow(as.data.frame(logact.basis[[1]])))
+ bs <- xx[i, j] * logact.basis[[match(colnames(xx)[j], basisnames)]]
+ if(!is.matrix(bs)) bs <- matrix(bs, byrow = TRUE, nrow = nrow(as.data.frame(logact.basis[[1]])))
if(ncol(bs) == 1) b <- matrix(b)
b <- b - bs
}
@@ -196,22 +193,21 @@
for(i in 1:nrow(B[[1]])) {
for(j in 1:ncol(B[[1]])) {
b <- numeric()
- for(k in 1:length(B)) b <- c(b,B[[k]][i,j])
+ for(k in 1:length(B)) b <- c(b, B[[k]][i, j])
AAA <- A
# solve for the activities in the buffer
- t <- solve(AAA,b)
+ t <- solve(AAA, b)
for(k in 1:length(ibasis))
- X[[k]][i,j] <- t[k]
+ X[[k]][i, j] <- t[k]
}
}
# Store results
for(i in 1:length(ibasis)) {
if(ncol(X[[i]]) == 1) X[[i]] <- as.numeric(X[[i]])
- else if(nrow(X[[i]]) == 1) X[[i]] <- as.matrix(X[[i]],nrow = 1)
+ else if(nrow(X[[i]]) == 1) X[[i]] <- as.matrix(X[[i]], nrow = 1)
logact.basis[[ibasis[i]]] <- X[[i]]
}
names(logact.basis) <- basisnames
- return(list(ibasis = ibasis,logact.basis = logact.basis))
+ return(list(ibasis = ibasis, logact.basis = logact.basis))
}
-
Modified: pkg/CHNOSZ/vignettes/anintro.Rmd
===================================================================
--- pkg/CHNOSZ/vignettes/anintro.Rmd 2024-01-09 02:23:52 UTC (rev 822)
+++ pkg/CHNOSZ/vignettes/anintro.Rmd 2024-01-11 12:10:40 UTC (rev 823)
@@ -801,19 +801,17 @@
```{r PPM_basis, results="hide", message=FALSE}
basis(c("FeS2", "H2S", "O2", "H2O"))
species(c("pyrite", "magnetite"))
-species("pyrrhotite", "cr2")
+species("pyrrhotite", "cr2", add = TRUE)
```
```{marginfigure}
The affinity of formation of pyrite happens to be zero because it is identical to one of the selected basis species.
```
-```{r PPM_affinity, message=FALSE, echo=1}
+```{r PPM_affinity, message=FALSE}
unlist(affinity(T = 300, P = 100)$values)
-## 2031 1999 2036
-## 0 0 0
```
We use <span style="color:red">`mod.buffer()`</span> to choose the `cr2` phase of pyrrhotite, which is stable at this temperature ([see above](#mosaic-diagrams) for how to get this information for minerals with polymorphic transitions).
-Then, we set up H<sub>2</sub>S and `r o2` to be buffered by PPM, and inspect their buffered activities:
+Then, we set up H<sub>2</sub>S and `r o2` to be buffered by PPM, and inspect their buffered activities (logarithmic values):
```{r PPM_setup, results="hide"}
mod.buffer("PPM", "pyrrhotite", "cr2")
basis(c("H2S", "O2"), c("PPM", "PPM"))
@@ -826,11 +824,15 @@
```{r demo_buffer_noecho, fig.margin=TRUE, fig.width=4, fig.height=4, out.width="100%", message=FALSE, echo=FALSE, cache=TRUE, fig.cap="Values of log<i>f</i><sub>H<sub>2</sub></sub> corresponding to mineral buffers or to given activities of aqueous species.", pngquant=pngquant, timeit=timeit}
demo(buffer, echo = FALSE)
```
-Et voilà! We have found log*a*<sub>H<sub>2</sub>S</sub> and `r logfO2` that are compatible with the coexistence of the three minerals.
+We have found log*a*<sub>H<sub>2</sub>S</sub> and `r logfO2` that are compatible with the coexistence of the three minerals.
Under these conditions, the affinities of formation reactions of the minerals in the buffer are all equal to zero:
-<!-- show static (commented) results because demo(buffer) changes the system -->
-```{r PPM_affinity, eval=FALSE}
+<!-- re-run the commands from above because demo(buffer) changes the system -->
+```{r PPM_basis, echo = FALSE, results = "hide"}
```
+```{r PPM_setup, echo = FALSE, results = "hide", message = FALSE}
+```
+```{r PPM_affinity, message = FALSE}
+```
Another example, based on Figure 6 of @SS95, is given in [<span style="color:blue">`demo(buffer)`</span>](../demo).
Here, values of log*f*<sub>H<sub>2</sub></sub> buffered by minerals or set by equilibrium with given activities of aqueous species are calculated using the two methods:
More information about the CHNOSZ-commits
mailing list