[CHNOSZ-commits] r821 - in pkg/CHNOSZ: . R inst/tinytest
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Dec 6 03:10:37 CET 2023
Author: jedick
Date: 2023-12-06 03:10:37 +0100 (Wed, 06 Dec 2023)
New Revision: 821
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/R/affinity.R
pkg/CHNOSZ/inst/tinytest/test-buffer.R
Log:
Minor code cleanup in affinity()
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2023-12-06 00:00:43 UTC (rev 820)
+++ pkg/CHNOSZ/DESCRIPTION 2023-12-06 02:10:37 UTC (rev 821)
@@ -1,6 +1,6 @@
-Date: 2023-12-05
+Date: 2023-12-06
Package: CHNOSZ
-Version: 2.0.0-40
+Version: 2.0.0-41
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 2023-12-06 00:00:43 UTC (rev 820)
+++ pkg/CHNOSZ/R/affinity.R 2023-12-06 02:10:37 UTC (rev 821)
@@ -11,12 +11,15 @@
#source("util.args.R")
#source("util.data.R")
#source("species.R")
+#source("info.R")
+#source("hkf.R")
+#source("cgl.R")
affinity <- function(..., property = NULL, sout = NULL, exceed.Ttr = FALSE, exceed.rhomin = FALSE,
return.buffer = FALSE, return.sout = FALSE, balance = "PBB", iprotein = NULL, loga.protein = 0, transect = NULL) {
# ...: variables over which to calculate
# property: what type of energy
- # (G.basis,G.species,logact.basis,logK,logQ,A)
+ # (G.basis, G.species, logact.basis, logK, logQ, A)
# return.buffer: return buffered activities
# balance: balance protein buffers on PBB
# exceed.Ttr: extrapolate Gibbs energies
@@ -60,7 +63,7 @@
# The user just wants an energy property
buffer <- FALSE
args$what <- property
- out <- do.call("energy",args)
+ out <- do.call("energy", args)
a <- out$a
sout <- out$sout
} else {
@@ -68,7 +71,7 @@
# Affinity calculations
property <- args$what
- # iprotein stuff
+ # Protein stuff
# Note that affinities of the residues are modified by ionization calculations in energy(), not here
if(!is.null(iprotein)) {
# Check all proteins are available
@@ -75,10 +78,11 @@
if(any(is.na(iprotein))) stop("`iprotein` has some NA values")
if(!all(iprotein %in% 1:nrow(thermo$protein))) stop("some value(s) of `iprotein` are not rownumbers of thermo()$protein")
# Add protein residues to the species list
- resnames <- c("H2O",aminoacids(3))
+ 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)
ires <- match(resprot, thermo$species$name)
}
@@ -93,18 +97,17 @@
buffer <- TRUE
message('affinity: loading buffer species')
if(!is.null(thermo$species)) is.species <- 1:nrow(thermo$species) else is.species <- numeric()
+ # Load the species in the buffer and get their species number(s)
is.buffer <- buffer(logK = NULL)
+ # Re-read thermo because the preceding command changed the species
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]
buffers <- names(is.buffer)
- # Reorder the buffers according to thermo$buffer
- buffers <- buffers[order(match(buffers,thermo$buffer$name))]
+ # Find species that are only in the buffer, not in the starting species list
+ is.only.buffer <- setdiff(unlist(is.buffer), is.species)
}
# Here we call 'energy'
- aa <- do.call("energy",args)
+ aa <- do.call("energy", args)
a <- aa$a
sout <- aa$sout
@@ -114,7 +117,7 @@
if(buffer) {
args$what <- "logact.basis"
args$sout <- sout
- logact.basis.new <- logact.basis <- do.call("energy",args)$a
+ 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])
@@ -123,11 +126,11 @@
for(i in 1:length(logK)) {
logK[[i]] <- logK[[i]] + thermo$species$logact[i]
for(j in 1:length(logact.basis.new)) {
- logK[[i]] <- logK[[i]] - logact.basis.new[[j]] * thermo$species[i,j]
+ 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 = as.numeric(is.buffer[[which(names(is.buffer) == buffers[k])]]),balance = balance)
+ 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]]
# Calculation of the buffered activities' effect on chemical affinities
is.only.buffer.new <- is.only.buffer[is.only.buffer %in% is.buffer[[k]]]
@@ -139,16 +142,16 @@
if(!j %in% ibufbasis) 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))
+ a[[i]] <- aa + (logact.basis.new[[j]] - logact.basis[[j]]) * thermo$species[i, j]
+ #if(!identical(a[[i]], aa)) print(paste(i, j))
}
}
if(k == length(buffers) & return.buffer) {
logact.basis.new <- lbn[[2]]
- ibasis.new <- c(ibasis.new,lbn[[1]])
- } else ibasis.new <- c(ibasis.new,ibasis)
+ ibasis.new <- c(ibasis.new, lbn[[1]])
+ } else ibasis.new <- c(ibasis.new, ibasis)
}
- species(is.only.buffer,delete = TRUE)
+ species(is.only.buffer, delete = TRUE)
if(length(is.only.buffer) > 0) a <- a[-is.only.buffer]
# To return the activities of buffered basis species
tb <- logact.basis.new[unique(ibasis.new)]
@@ -159,9 +162,9 @@
for(i in 1:length(tb)) {
#tb[[i]] <- as.data.frame(tb[[i]])
if(nd > 0) rownames(tb[[i]]) <-
- seq(args$lims[[1]][1],args$lims[[1]][2],length.out = args$lims[[1]][3])
+ seq(args$lims[[1]][1], args$lims[[1]][2], length.out = args$lims[[1]][3])
if(nd > 1) colnames(tb[[i]]) <-
- seq(args$lims[[2]][1],args$lims[[2]][2],length.out = args$lims[[2]][3])
+ seq(args$lims[[2]][1], args$lims[[2]][2], length.out = args$lims[[2]][3])
}
}
}
@@ -171,10 +174,10 @@
if(!is.null(iprotein)) {
# Fast protein calculations 20090331
# Function to calculate affinity of formation reactions from those of residues
- loga.protein <- rep(loga.protein,length.out = length(iprotein))
+ 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])
+ tpext <- as.numeric(thermo$protein[iprotein[ip], 5:25])
+ return(Reduce("+", pprod(a[ires], tpext)) - loga.protein[ip])
}
# Use another level of indexing to let the function report on its progress
jprotein <- 1:length(iprotein)
@@ -186,19 +189,19 @@
# The current species list, containing the residues
resspecies <- thermo$species
# Now we can delete the residues from the species list
- species(ires,delete = TRUE)
+ species(ires, delete = TRUE)
# State and protein names
state <- resspecies$state[1]
- name <- paste(thermo$protein$protein[iprotein],thermo$protein$organism[iprotein],sep = "_")
+ 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(mybasis)])) %*% 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)
+ protspecies <- cbind(protbasis, data.frame(ispecies = ispecies, logact = loga.protein, state = state, name = name))
+ myspecies <- rbind(myspecies, protspecies)
rownames(myspecies) <- 1:nrow(myspecies)
## Update the affinity values
names(protein.affinity) <- ispecies
- a <- c(a,protein.affinity)
+ a <- c(a, protein.affinity)
a <- a[-ires]
}
@@ -266,4 +269,3 @@
if(buffer) a <- c(a, list(buffer = tb))
return(a)
}
-
Modified: pkg/CHNOSZ/inst/tinytest/test-buffer.R
===================================================================
--- pkg/CHNOSZ/inst/tinytest/test-buffer.R 2023-12-06 00:00:43 UTC (rev 820)
+++ pkg/CHNOSZ/inst/tinytest/test-buffer.R 2023-12-06 02:10:37 UTC (rev 821)
@@ -33,6 +33,7 @@
expect_equivalent(unique(as.vector(A2$buffer[[1]])), logaSiO2[1], info = info)
# A test for 0 dimensions with a two-mineral buffer 20201103
+info <- "2-mineral buffer at a single point"
# (buffer errored trying to index columns of non-matrix prior to 1.3.6-85)
T <- 400
P <- 1000
More information about the CHNOSZ-commits
mailing list