[CHNOSZ-commits] r560 - in pkg/CHNOSZ: . R inst tests/testthat
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Jul 15 03:04:39 CEST 2020
Author: jedick
Date: 2020-07-15 03:04:37 +0200 (Wed, 15 Jul 2020)
New Revision: 560
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/R/equilibrate.R
pkg/CHNOSZ/R/mosaic.R
pkg/CHNOSZ/R/species.R
pkg/CHNOSZ/inst/NEWS.Rd
pkg/CHNOSZ/tests/testthat/test-equilibrate.R
Log:
Make equilibrate() work for only solids
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2020-07-14 13:55:14 UTC (rev 559)
+++ pkg/CHNOSZ/DESCRIPTION 2020-07-15 01:04:37 UTC (rev 560)
@@ -1,6 +1,6 @@
-Date: 2020-07-14
+Date: 2020-07-15
Package: CHNOSZ
-Version: 1.3.6-33
+Version: 1.3.6-34
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/equilibrate.R
===================================================================
--- pkg/CHNOSZ/R/equilibrate.R 2020-07-14 13:55:14 UTC (rev 559)
+++ pkg/CHNOSZ/R/equilibrate.R 2020-07-15 01:04:37 UTC (rev 560)
@@ -36,78 +36,88 @@
n.balance.orig <- n.balance <- bout$n.balance
balance <- bout$balance
## if solids (cr) species are present, find them on a predominance diagram 20191111
- hascr <- any(grepl("cr", aout$species$state))
- if(hascr) dout <- diagram(aout, balance = balance, normalize = normalize, as.residue = as.residue, plot.it = FALSE, limit.water = FALSE)
- ## take selected species in 'ispecies'
- if(length(ispecies)==0) stop("the length of ispecies is zero")
- if(is.logical(ispecies)) ispecies <- which(ispecies)
- # take out species that have NA affinities
- ina <- sapply(aout$values, function(x) all(is.na(x)))
- ispecies <- ispecies[!ina[ispecies]]
- if(length(ispecies)==0) stop("all species have NA affinities")
- if(!identical(ispecies, 1:nspecies)) {
- message(paste("equilibrate: using", length(ispecies), "of", nspecies, "species"))
- aout$species <- aout$species[ispecies, ]
- aout$values <- aout$values[ispecies]
- n.balance <- n.balance[ispecies]
- }
- ## number of species that are left
- nspecies <- length(aout$values)
- ## say what the balancing coefficients are
- if(length(n.balance) < 100) message(paste("equilibrate: n.balance is", c2s(n.balance)))
- ## logarithm of total activity of the balancing basis species
- if(is.null(loga.balance)) {
- # sum up the activities, then take absolute value
- # in case n.balance is negative
- sumact <- abs(sum(10^aout$species$logact * n.balance))
- loga.balance <- log10(sumact)
- }
- # make loga.balance the same length as the values of affinity
- loga.balance <- unlist(loga.balance)
- nvalues <- length(unlist(aout$values[[1]]))
- if(length(loga.balance) == 1) {
- # we have a constant loga.balance
- message(paste0("equilibrate: loga.balance is ", loga.balance))
- loga.balance <- rep(loga.balance, nvalues)
+ iscr <- grepl("cr", aout$species$state)
+ ncr <- sum(iscr)
+ if(ncr > 0) dout <- diagram(aout, balance = balance, normalize = normalize, as.residue = as.residue, plot.it = FALSE, limit.water = FALSE)
+ if(ncr == nspecies) {
+ ## we get here if there are only solids 20200714
+ m.balance <- NULL
+ Astar <- NULL
+ loga.equil <- aout$values
+ for(i in 1:length(loga.equil)) loga.equil[[i]][] <- NA
} else {
- # we are using a variable loga.balance (supplied by the user)
- if(!identical(length(loga.balance), nvalues)) stop("length of loga.balance (", length(loga.balance), ") doesn't match the affinity values (", nvalues, ")")
- message(paste0("equilibrate: loga.balance has same length as affinity values (", length(loga.balance), ")"))
- }
- ## normalize -- normalize the molar formula by the balance coefficients
- m.balance <- n.balance
- isprotein <- grepl("_", as.character(aout$species$name))
- if(any(normalize) | as.residue) {
- if(any(n.balance < 0)) stop("one or more negative balancing coefficients prohibit using normalized molar formulas")
- n.balance[normalize|as.residue] <- 1
- if(as.residue) message(paste("equilibrate: using 'as.residue' for molar formulas"))
- else message(paste("equilibrate: using 'normalize' for molar formulas"))
- # set the formula divisor (m.balance) to 1 for species whose formulas are *not* normalized
- m.balance[!(normalize|as.residue)] <- 1
- } else m.balance[] <- 1
- ## Astar: the affinities/2.303RT of formation reactions with
- ## formed species in their standard-state activities
- Astar <- lapply(1:nspecies, function(i) {
- # 'starve' the affinity of the activity of the species,
- # and normalize the value by the nor-molar ratio
- (aout$values[[i]] + aout$species$logact[i]) / m.balance[i]
- })
- ## chose a method and compute the equilibrium activities of species
- if(missing(method)) {
- if(all(n.balance==1)) method <- method[1]
- else method <- method[2]
- }
- message(paste("equilibrate: using", method[1], "method"))
- if(method[1]=="boltzmann") loga.equil <- equil.boltzmann(Astar, n.balance, loga.balance)
- else if(method[1]=="reaction") loga.equil <- equil.reaction(Astar, n.balance, loga.balance, tol)
- ## if we normalized the formulas, get back to activities to species
- if(any(normalize) & !as.residue) {
- loga.equil <- lapply(1:nspecies, function(i) {
- loga.equil[[i]] - log10(m.balance[i])
+ ## we get here if there are any aqueous species 20200714
+ ## take selected species in 'ispecies'
+ if(length(ispecies)==0) stop("the length of ispecies is zero")
+ if(is.logical(ispecies)) ispecies <- which(ispecies)
+ # take out species that have NA affinities
+ ina <- sapply(aout$values, function(x) all(is.na(x)))
+ ispecies <- ispecies[!ina[ispecies]]
+ if(length(ispecies)==0) stop("all species have NA affinities")
+ if(!identical(ispecies, 1:nspecies)) {
+ message(paste("equilibrate: using", length(ispecies), "of", nspecies, "species"))
+ aout$species <- aout$species[ispecies, ]
+ aout$values <- aout$values[ispecies]
+ n.balance <- n.balance[ispecies]
+ }
+ ## number of species that are left
+ nspecies <- length(aout$values)
+ ## say what the balancing coefficients are
+ if(length(n.balance) < 100) message(paste("equilibrate: n.balance is", c2s(n.balance)))
+ ## logarithm of total activity of the balancing basis species
+ if(is.null(loga.balance)) {
+ # sum up the activities, then take absolute value
+ # in case n.balance is negative
+ sumact <- abs(sum(10^aout$species$logact * n.balance))
+ loga.balance <- log10(sumact)
+ }
+ # make loga.balance the same length as the values of affinity
+ loga.balance <- unlist(loga.balance)
+ nvalues <- length(unlist(aout$values[[1]]))
+ if(length(loga.balance) == 1) {
+ # we have a constant loga.balance
+ message(paste0("equilibrate: loga.balance is ", loga.balance))
+ loga.balance <- rep(loga.balance, nvalues)
+ } else {
+ # we are using a variable loga.balance (supplied by the user)
+ if(!identical(length(loga.balance), nvalues)) stop("length of loga.balance (", length(loga.balance), ") doesn't match the affinity values (", nvalues, ")")
+ message(paste0("equilibrate: loga.balance has same length as affinity values (", length(loga.balance), ")"))
+ }
+ ## normalize -- normalize the molar formula by the balance coefficients
+ m.balance <- n.balance
+ isprotein <- grepl("_", as.character(aout$species$name))
+ if(any(normalize) | as.residue) {
+ if(any(n.balance < 0)) stop("one or more negative balancing coefficients prohibit using normalized molar formulas")
+ n.balance[normalize|as.residue] <- 1
+ if(as.residue) message(paste("equilibrate: using 'as.residue' for molar formulas"))
+ else message(paste("equilibrate: using 'normalize' for molar formulas"))
+ # set the formula divisor (m.balance) to 1 for species whose formulas are *not* normalized
+ m.balance[!(normalize|as.residue)] <- 1
+ } else m.balance[] <- 1
+ ## Astar: the affinities/2.303RT of formation reactions with
+ ## formed species in their standard-state activities
+ Astar <- lapply(1:nspecies, function(i) {
+ # 'starve' the affinity of the activity of the species,
+ # and normalize the value by the nor-molar ratio
+ (aout$values[[i]] + aout$species$logact[i]) / m.balance[i]
})
+ ## chose a method and compute the equilibrium activities of species
+ if(missing(method)) {
+ if(all(n.balance==1)) method <- method[1]
+ else method <- method[2]
+ }
+ message(paste("equilibrate: using", method[1], "method"))
+ if(method[1]=="boltzmann") loga.equil <- equil.boltzmann(Astar, n.balance, loga.balance)
+ else if(method[1]=="reaction") loga.equil <- equil.reaction(Astar, n.balance, loga.balance, tol)
+ ## if we normalized the formulas, get back to activities to species
+ if(any(normalize) & !as.residue) {
+ loga.equil <- lapply(1:nspecies, function(i) {
+ loga.equil[[i]] - log10(m.balance[i])
+ })
+ }
}
## process cr species 20191111
- if(hascr) {
+ if(ncr > 0) {
# cr species were excluded from equilibrium calculation, so get values back to original lengths
norig <- length(dout$values)
n.balance <- n.balance.orig
Modified: pkg/CHNOSZ/R/mosaic.R
===================================================================
--- pkg/CHNOSZ/R/mosaic.R 2020-07-14 13:55:14 UTC (rev 559)
+++ pkg/CHNOSZ/R/mosaic.R 2020-07-15 01:04:37 UTC (rev 560)
@@ -78,7 +78,7 @@
mysp <- species(bases[[i]])
# 20191111 include only aq species in total activity
iaq <- mysp$state == "aq"
- species(which(iaq), basis0$logact[ibasis0[i]])
+ if(any(iaq)) species(which(iaq), basis0$logact[ibasis0[i]])
A.bases[[i]] <- suppressMessages(affinity(..., sout = sout))
}
Modified: pkg/CHNOSZ/R/species.R
===================================================================
--- pkg/CHNOSZ/R/species.R 2020-07-14 13:55:14 UTC (rev 559)
+++ pkg/CHNOSZ/R/species.R 2020-07-15 01:04:37 UTC (rev 560)
@@ -86,7 +86,8 @@
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
+ # if species is numeric and a small number it refers to the species definition,
+ # else to species index in thermo()$OBIGT
nspecies <- nrow(thermo$species)
if(is.null(thermo$species)) nspecies <- 0
if(max(species) > nspecies) iOBIGT <- species
Modified: pkg/CHNOSZ/inst/NEWS.Rd
===================================================================
--- pkg/CHNOSZ/inst/NEWS.Rd 2020-07-14 13:55:14 UTC (rev 559)
+++ pkg/CHNOSZ/inst/NEWS.Rd 2020-07-15 01:04:37 UTC (rev 560)
@@ -5,7 +5,7 @@
% macros
\newcommand{\H2O}{\ifelse{latex}{\eqn{\mathrm{H_{2}O}}}{\ifelse{html}{\out{H<sub>2</sub>O}}{H2O}}}
-\section{Changes in CHNOSZ version 1.3.6-32 (2020-07-14)}{
+\section{Changes in CHNOSZ version 1.3.6-34 (2020-07-15)}{
\subsection{MAJOR CHANGES}{
\itemize{
@@ -95,11 +95,13 @@
\subsection{CHANGES TO FUNCTIONS}{
\itemize{
+ \item Change default resolution in \code{affinity()} from 128 to 256.
+
\item \code{subcrt()}: change \strong{action.unbalanced} argument to
\strong{autobalance}, which now provides the ability to prevent
autobalancing.
- \item Changing the water model with \code{water()} now updates the
+ \item Setting the water model with \code{water()} now updates the
literature references in \code{thermo()$OBIGT}.
\item \code{thermo.refs()} now shows CHNOSZ version and date.
@@ -154,8 +156,6 @@
\item TODO: OBIGT.Rmd: change "CHNOSZ" references to "OBIGT".
- \item TODO: change default plot resolution from 128 to 256.
-
}
}
Modified: pkg/CHNOSZ/tests/testthat/test-equilibrate.R
===================================================================
--- pkg/CHNOSZ/tests/testthat/test-equilibrate.R 2020-07-14 13:55:14 UTC (rev 559)
+++ pkg/CHNOSZ/tests/testthat/test-equilibrate.R 2020-07-15 01:04:37 UTC (rev 560)
@@ -203,6 +203,11 @@
e <- equilibrate(a)
epredom <- diagram(e, plot.it = FALSE)$predominant
expect_equal(apredom, epredom)
+ # also test that equilibrate() works with *only* solids 20200715
+ species(Cu_cr)
+ acr <- affinity(a)
+ ecr <- equilibrate(acr)
+ expect_identical(e$values[8:9], ecr$values)
})
# references
More information about the CHNOSZ-commits
mailing list