[CHNOSZ-commits] r514 - in pkg/CHNOSZ: . R inst man tests/testthat
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Nov 11 13:31:59 CET 2019
Author: jedick
Date: 2019-11-11 13:31:58 +0100 (Mon, 11 Nov 2019)
New Revision: 514
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/R/equilibrate.R
pkg/CHNOSZ/inst/NEWS
pkg/CHNOSZ/man/equilibrate.Rd
pkg/CHNOSZ/tests/testthat/test-equilibrate.R
Log:
equilibrate(): improve handling of solids
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2019-11-11 05:48:47 UTC (rev 513)
+++ pkg/CHNOSZ/DESCRIPTION 2019-11-11 12:31:58 UTC (rev 514)
@@ -1,6 +1,6 @@
Date: 2019-11-11
Package: CHNOSZ
-Version: 1.3.3-10
+Version: 1.3.3-11
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 2019-11-11 05:48:47 UTC (rev 513)
+++ pkg/CHNOSZ/R/equilibrate.R 2019-11-11 12:31:58 UTC (rev 514)
@@ -7,7 +7,7 @@
#source("util.character.R")
equilibrate <- function(aout, balance=NULL, loga.balance=NULL,
- ispecies=1:length(aout$values), normalize=FALSE, as.residue=FALSE,
+ ispecies = !grepl("cr", aout$species$state), normalize=FALSE, as.residue=FALSE,
method=c("boltzmann", "reaction"), tol=.Machine$double.eps^0.25) {
### calculate equilibrium activities of species from the affinities
### of their formation reactions from basis species at given activities
@@ -33,10 +33,14 @@
nspecies <- length(aout$values)
## get the balancing coefficients
bout <- balance(aout, balance)
- n.balance <- bout$n.balance
+ 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)
## 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]]
@@ -102,6 +106,39 @@
loga.equil[[i]] - log10(m.balance[i])
})
}
+ ## process cr species 20191111
+ if(hascr) {
+ # cr species were excluded from equilibrium calculation, so get values back to original lengths
+ norig <- length(dout$values)
+ n.balance <- n.balance.orig
+ imatch <- match(1:norig, ispecies)
+ m.balance <- m.balance[imatch]
+ Astar <- Astar[imatch]
+ loga.equil1 <- loga.equil[[1]]
+ loga.equil <- loga.equil[imatch]
+ # replace NULL loga.equil with input values (cr only)
+ ina <- which(is.na(imatch))
+ for(i in ina) {
+ loga.equil[[i]] <- loga.equil1
+ loga.equil[[i]][] <- dout$species$logact[[i]]
+ }
+ aout$species <- dout$species
+ aout$values <- dout$values
+ # find the grid points where any cr species is predominant
+ icr <- which(grepl("cr", dout$species$state))
+ iscr <- lapply(icr, function(x) dout$predominant == x)
+ iscr <- Reduce("|", iscr)
+ # at those grid points, make the aqueous species' activities practically zero
+ for(i in 1:norig) {
+ if(i %in% icr) next
+ loga.equil[[i]][iscr] <- -999
+ }
+ # at the other grid points, make the cr species' activities practically zero
+ for(i in icr) {
+ ispredom <- dout$predominant == i
+ loga.equil[[i]][!ispredom] <- -999
+ }
+ }
## put together the output
out <- c(aout, list(balance=balance, m.balance=m.balance, n.balance=n.balance,
loga.balance=loga.balance, Astar=Astar, loga.equil=loga.equil))
Modified: pkg/CHNOSZ/inst/NEWS
===================================================================
--- pkg/CHNOSZ/inst/NEWS 2019-11-11 05:48:47 UTC (rev 513)
+++ pkg/CHNOSZ/inst/NEWS 2019-11-11 12:31:58 UTC (rev 514)
@@ -1,4 +1,4 @@
-CHANGES IN CHNOSZ 1.3.3-10 (2019-11-11)
+CHANGES IN CHNOSZ 1.3.3-11 (2019-11-11)
---------------------------------------
- describe.reaction(): revert the change of using a double arrow. The
@@ -23,16 +23,22 @@
- Move some carbonates from SUPCRT92 back into default database:
artinite, azurite.
-- subcrt(): properties of minerals are now calculated including the
- upper T limit (or transition temperature, for e.g. cr1 -> cr2).
- Previously, properties were NA for temperature equal to the transition
- temperature.
+- subcrt(): properties of minerals are now returned at the upper T limit
+ (or transition temperature, for e.g. cr1 -> cr2). Previously,
+ properties were set to NA at (and not only above) the T limit.
-- mosaic(): set only the activities of aqueous basis species to the
- total activity (specified in the incoming basis() definition). This
+- mosaic(): set the activities of only aqueous basis species to the
+ total activity (taken from the incoming basis() definition). This
fixes a bug where activities of minerals (particularly, native sulfur)
- were changed from one.
+ were unexpectedly changed (their logact should be 0).
+- Improve the handling of solids in equilibrate(). They are now excluded
+ from the equilibrium calculation, but their stability fields are
+ calculated (using the maximum affinity method via diagram()). Where
+ solids are stable, logarithms of activities of aqueous species are set
+ to -999, and vice versa. Thanks to Feng Lai for the request and test
+ case.
+
CHANGES IN CHNOSZ 1.3.3 (2019-08-02)
------------------------------------
Modified: pkg/CHNOSZ/man/equilibrate.Rd
===================================================================
--- pkg/CHNOSZ/man/equilibrate.Rd 2019-11-11 05:48:47 UTC (rev 513)
+++ pkg/CHNOSZ/man/equilibrate.Rd 2019-11-11 12:31:58 UTC (rev 514)
@@ -10,11 +10,11 @@
}
\usage{
- equilibrate(aout, balance=NULL, loga.balance=NULL,
- ispecies=1:length(aout$values), normalize=FALSE, as.residue=FALSE,
- method=c("boltzmann", "reaction"), tol=.Machine$double.eps^0.25)
+ equilibrate(aout, balance = NULL, loga.balance = NULL,
+ ispecies = !grepl("cr", aout$species$state), normalize = FALSE, as.residue = FALSE,
+ method = c("boltzmann", "reaction"), tol = .Machine$double.eps^0.25)
equil.boltzmann(Astar, n.balance, loga.balance)
- equil.reaction(Astar, n.balance, loga.balance, tol=.Machine$double.eps^0.25)
+ equil.reaction(Astar, n.balance, loga.balance, tol = .Machine$double.eps^0.25)
moles(eout)
}
@@ -68,7 +68,11 @@
The default behavior can be overriden by specifying either \samp{boltzmann} or \samp{reaction} in \code{method}.
Using \code{equil.reaction} may be needed for systems with huge (negative or positive) affinities, where \code{equil.boltzmann} produces a NaN result.
-\code{ispecies} can be supplied to identify a subset of the species to include in the calculation.
+\code{ispecies} can be supplied to identify a subset of the species to include in the equilibrium calculation.
+By default, this is all species except solids (species with \samp{cr} state).
+However, the stability regions of solids are still calculated (by a call to \code{\link{diagram}} without plotting).
+At all points outside of their stability region, the logarithms of activities of solids are set to -999.
+Likewise, where any solid species is calculated to be stable, the logarithms of activities of all aqueous species are set to -999.
\code{moles} simply calculates the total number of moles of elements corresponding to the activities of formed species in the output from \code{equilibrate}.
}
Modified: pkg/CHNOSZ/tests/testthat/test-equilibrate.R
===================================================================
--- pkg/CHNOSZ/tests/testthat/test-equilibrate.R 2019-11-11 05:48:47 UTC (rev 513)
+++ pkg/CHNOSZ/tests/testthat/test-equilibrate.R 2019-11-11 12:31:58 UTC (rev 514)
@@ -190,6 +190,22 @@
expect_true(maxdiff(dloga_isoalkane_mix, dloga_isoalkane_norm) > maxdiff(dloga_isoalkane_mix, dloga_isoalkane_full))
})
+test_that("solids are not equilibrated, but their stability fields are calculated", {
+ # added 20191111; based on an example sent by Feng Lai on 20191020
+ Cu_aq <- c("CuCl", "CuCl2-", "CuCl3-2", "CuHS", "Cu(HS)2-", "CuOH", "Cu(OH)2-")
+ Cu_cr <- c("copper", "chalcocite")
+ basis(c("Cu+", "HS-", "Cl-", "H2O", "H+", "oxygen"))
+ basis("O2", -35)
+ basis("H+", -5)
+ species(Cu_aq, -3)
+ species(Cu_cr)
+ a <- affinity("Cl-" = c(-3, 0, 200), "HS-" = c(-10, 0, 200), T = 325, P = 500)
+ apredom <- diagram(a, plot.it = FALSE)$predominant
+ e <- equilibrate(a)
+ epredom <- diagram(e, plot.it = FALSE)$predominant
+ expect_equal(apredom, epredom)
+})
+
# references
# Seewald, J. S. (2001)
More information about the CHNOSZ-commits
mailing list