[CHNOSZ-commits] r790 - in pkg/CHNOSZ: . R inst inst/tinytest man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Jun 20 13:47:58 CEST 2023
Author: jedick
Date: 2023-06-20 13:47:58 +0200 (Tue, 20 Jun 2023)
New Revision: 790
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/R/subcrt.R
pkg/CHNOSZ/inst/NEWS.Rd
pkg/CHNOSZ/inst/tinytest/test-subcrt.R
pkg/CHNOSZ/man/subcrt.Rd
Log:
Add 'use.polymorphs' argument to subcrt()
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2023-06-20 06:44:10 UTC (rev 789)
+++ pkg/CHNOSZ/DESCRIPTION 2023-06-20 11:47:58 UTC (rev 790)
@@ -1,6 +1,6 @@
Date: 2023-06-20
Package: CHNOSZ
-Version: 2.0.0-9
+Version: 2.0.0-10
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/subcrt.R
===================================================================
--- pkg/CHNOSZ/R/subcrt.R 2023-06-20 06:44:10 UTC (rev 789)
+++ pkg/CHNOSZ/R/subcrt.R 2023-06-20 11:47:58 UTC (rev 790)
@@ -16,7 +16,7 @@
subcrt <- function(species, coeff = 1, state = NULL, property = c("logK", "G", "H", "S", "V", "Cp"),
T = seq(273.15, 623.15, 25), P = "Psat", grid = NULL, convert = TRUE, exceed.Ttr = FALSE,
- exceed.rhomin = FALSE, logact = NULL, autobalance = TRUE, IS = 0) {
+ exceed.rhomin = FALSE, logact = NULL, autobalance = TRUE, use.polymorphs = TRUE, IS = 0) {
# Revise the call if the states are the second argument
if(!is.null(coeff[1])) {
@@ -118,9 +118,8 @@
# Get species information
thermo <- get("thermo", CHNOSZ)
- # Pre-20110808, we sent numeric species argument through info() to get species name and state(s)
+ # Before 20110808, we sent numeric species argument through info() to get species name and state(s)
# But why slow things down if we already have a species index?
- # ... so now polymorph stuff will only work for character species names
if(is.numeric(species[1])) {
ispecies <- species
species <- as.character(thermo$OBIGT$name[ispecies])
@@ -127,28 +126,28 @@
state <- as.character(thermo$OBIGT$state[ispecies])
newstate <- as.character(thermo$OBIGT$state[ispecies])
} else {
- # From names, get species indices and states and possibly
- # keep track of polymorphs (cr,cr2 ...)
+ # Species are named ... look up the index
ispecies <- numeric()
newstate <- character()
for(i in 1:length(species)) {
# Get the species index for a named species
- if(!can.be.numeric(species[i])) si <- info.character(species[i], state[i])
+ if(!can.be.numeric(species[i])) sindex <- info.character(species[i], state[i])
else {
- # Check that a numeric argument is a rownumber of thermo()$OBIGT
- si <- as.numeric(species[i])
- if(!si %in% 1:nrow(thermo$OBIGT)) stop(paste(species[i], "is not a row number of thermo()$OBIGT"))
+ # Check that a coerced-to-numeric argument is a rownumber of thermo()$OBIGT
+ sindex <- as.numeric(species[i])
+ if(!sindex %in% 1:nrow(thermo$OBIGT)) stop(paste(species[i], "is not a row number of thermo()$OBIGT"))
}
- # That could have the side-effect of adding a protein; re-read thermo
+ # info.character() has the possible side-effect of adding a protein; re-read thermo to use the (possible) additions
thermo <- get("thermo", CHNOSZ)
- if(is.na(si[1])) stop("no info found for ", species[i], " ",state[i])
+ if(is.na(sindex[1])) stop("no info found for ", species[i], " ",state[i])
if(!is.null(state[i])) is.cr <- state[i]=="cr" else is.cr <- FALSE
- if(thermo$OBIGT$state[si[1]] == "cr" & (is.null(state[i]) | is.cr)) {
+ # If we found multiple species, take the first one (useful for minerals with polymorphs)
+ if(thermo$OBIGT$state[sindex[1]] == "cr" & (is.null(state[i]) | is.cr)) {
newstate <- c(newstate, "cr")
- ispecies <- c(ispecies, si[1])
+ ispecies <- c(ispecies, sindex[1])
} else {
- newstate <- c(newstate, as.character(thermo$OBIGT$state[si[1]]))
- ispecies <- c(ispecies, si[1])
+ newstate <- c(newstate, as.character(thermo$OBIGT$state[sindex[1]]))
+ ispecies <- c(ispecies, sindex[1])
}
}
}
@@ -162,17 +161,19 @@
state <- as.character(thermo$OBIGT$state[ispecies])
name <- as.character(thermo$OBIGT$name[ispecies])
# A counter of all species considered
- # iphases is longer than ispecies if cr,cr2 ... phases are present
+ # iphases is longer than ispecies if multiple polymorphs (cr, cr2, ...) are present
# polymorph.species shows which of ispecies correspond to iphases
- # pre-20091114: the success of this depends on there not being duplicated aqueous or other
- # non-mineral-polymorphs (i.e., two entries in OBIGT for Cu+ screw this up when running the skarn example).
- # after 20091114: we can deal with duplicated species (aqueous at least)
+ # Before 20091114: the success of this depends on there not being duplicated aqueous or other
+ # non-mineral species (i.e., two entries in OBIGT for Cu+ mess this up when running the skarn example).
+ # After 20091114: we can deal with duplicated species (aqueous at least)
iphases <- polymorph.species <- coeff.new <- numeric()
for(i in 1:length(ispecies)) {
- if(newstate[i] == "cr") {
- searchstates <- c("cr", "cr2", "cr3", "cr4", "cr5", "cr6", "cr7", "cr8", "cr9")
- tghs <- thermo$OBIGT[(thermo$OBIGT$name %in% name[i]) & thermo$OBIGT$state %in% searchstates, ]
- # we only take one if they are in fact duplicated species and not polymorphs
+ # Add check for use.polymorphs argument 20230620
+ if(newstate[i] == "cr" & use.polymorphs) {
+ # Check for available polymorphs in OBIGT
+ polymorph.states <- c("cr", "cr2", "cr3", "cr4", "cr5", "cr6", "cr7", "cr8", "cr9")
+ tghs <- thermo$OBIGT[(thermo$OBIGT$name %in% name[i]) & thermo$OBIGT$state %in% polymorph.states, ]
+ # Only take the first one if they are duplicated non-mineral species (i.e., not polymorphs)
if(all(tghs$state == tghs$state[1])) tghs <- thermo$OBIGT[ispecies[i], ]
} else tghs <- thermo$OBIGT[ispecies[i], ]
iphases <- c(iphases, as.numeric(rownames(tghs)))
@@ -354,12 +355,12 @@
# Name and state
myname <- reaction$name[i]
mystate <- reaction$state[i]
- # If we are considering multiple phases, and if this phase is cr2 or higher, check if we're below the transition temperature
+ # If we are considering multiple polymorphs, and if this polymorph is cr2 or higher, check if we're below the transition temperature
if(length(iphases) > length(ispecies) & i > 1) {
if(!(reaction$state[i] %in% c("liq", "cr", "gas")) & reaction$name[i-1] == reaction$name[i]) {
# After add.OBIGT("SUPCRT92"), quartz cr and cr2 are not next to each other in thermo()$OBIGT,
# so use iphases[i-1] here, not iphases[i]-1 20181107
- Ttr <- Ttr(iphases[i-1], iphases[i], P=P, dPdT = dPdTtr(iphases[i-1], iphases[i]))
+ Ttr <- Ttr(iphases[i-1], iphases[i], P = P, dPdT = dPdTtr(iphases[i-1], iphases[i]))
if(all(is.na(Ttr))) next
if(any(T <= Ttr)) {
status.Ttr <- "(extrapolating G)"
@@ -443,8 +444,8 @@
npolymorphs <- length(are.polymorphs) / ndups
are.polymorphs <- are.polymorphs[1:npolymorphs]
}
- if(length(are.polymorphs)>1) {
- message(paste("subcrt:", length(are.polymorphs), "phases for", thermo$OBIGT$name[ispecies[i]], "... "), appendLF = FALSE)
+ if(length(are.polymorphs) > 1) {
+ message(paste("subcrt:", length(are.polymorphs), "polymorphs for", thermo$OBIGT$name[ispecies[i]], "... "), appendLF = FALSE)
# Assemble the Gibbs energies for each species
for(j in 1:length(are.polymorphs)) {
G.this <- outprops[[are.polymorphs[j]]]$G
@@ -451,7 +452,7 @@
if(sum(is.na(G.this)) > 0 & exceed.Ttr) warning(paste("subcrt: NAs found for G of ",
reaction$name[are.polymorphs[j]], " ", reaction$state[are.polymorphs[j]], " at T-P point(s) ",
c2s(which(is.na(G.this)), sep = " "), sep = ""), call. = FALSE)
- if(j==1) G <- as.data.frame(G.this)
+ if(j == 1) G <- as.data.frame(G.this)
else G <- cbind(G, as.data.frame(G.this))
}
# Find the minimum-energy polymorph at each T-P point
@@ -461,11 +462,11 @@
ps <- which.min(as.numeric(G[j, ]))
if(length(ps)==0) {
# minimum not found (we have NAs)
- # - no non-NA value of G to begin with, e.g. aegerine) --> probably should use lowest-T phase
+ # - no non-NA value of G to begin with (e.g. aegerine) --> probably should use lowest-T phase
#ps <- 1
# - above temperature limit for the highest-T phase (subcrt.Rd skarn example) --> use highest-T phase 20171110
ps <- ncol(G)
- if(exceed.Ttr) warning("subcrt: stable phase for ", reaction$name[are.polymorphs[ps]], " at T-P point ", j,
+ if(exceed.Ttr) warning("subcrt: stable polymorph for ", reaction$name[are.polymorphs[ps]], " at T-P point ", j,
" undetermined (using ", reaction$state[are.polymorphs[ps]], ")", call. = FALSE)
}
stable.polymorph <- c(stable.polymorph, ps)
Modified: pkg/CHNOSZ/inst/NEWS.Rd
===================================================================
--- pkg/CHNOSZ/inst/NEWS.Rd 2023-06-20 06:44:10 UTC (rev 789)
+++ pkg/CHNOSZ/inst/NEWS.Rd 2023-06-20 11:47:58 UTC (rev 790)
@@ -12,12 +12,15 @@
% links to vignettes 20220723
\newcommand{\viglink}{\ifelse{html}{\out{<a href="../CHNOSZ/doc/#1.html"><strong>#1.Rmd</strong></a>}}{\bold{#1.Rmd}}}
-\section{Changes in CHNOSZ version 2.0.0-5 (2023-05-15)}{
+\section{Changes in CHNOSZ version 2.0.0-10 (2023-06-20)}{
\itemize{
- \item Restore EOSregress.R, eos-regress.Rmd, and demo/adenine.R.
+ \item Restore \file{EOSregress.R}, \file{eos-regress.Rmd}, and \file{demo/adenine.R}.
+ \item Add \code{use.polymorphs} argument to \code{subcrt()} to allow
+ turning off automatic identification of stable polymorphs.
+
}
}
Modified: pkg/CHNOSZ/inst/tinytest/test-subcrt.R
===================================================================
--- pkg/CHNOSZ/inst/tinytest/test-subcrt.R 2023-06-20 06:44:10 UTC (rev 789)
+++ pkg/CHNOSZ/inst/tinytest/test-subcrt.R 2023-06-20 11:47:58 UTC (rev 790)
@@ -202,6 +202,23 @@
# Before version 1.1.3-63, having more than one invalid property gave a mangled error message
expect_error(subcrt("H2O", -1, "liq", c(1, 2)), "invalid property names: 1 2", info = info)
+# Added on 20230620
+info <- "Polymorphs are used by default"
+sres_poly <- subcrt("pyrrhotite")
+expect_equal(unique(sres_poly$out[[1]]$polymorph), c(1, 2, 3), info = info)
+info <- "Polymorphs work for named species or numeric indices"
+iPo <- info("pyrrhotite")
+sres_poly1 <- subcrt(iPo)
+expect_identical(sres_poly, sres_poly1, info = info)
+info <- "Automatic identificatio of polymorphs can be turned off"
+sres_nopoly <- subcrt("pyrrhotite", use.polymorphs = FALSE)
+expect_null(sres_nopoly$out[[1]]$polymorph, info = info)
+info <- "Gibbs energy is NA beyond the transition temperature"
+expect_true(anyNA(sres_nopoly$out[[1]]$G))
+info <- "Gibbs energy can be extrapolated beyond the transition temperature"
+sres_nopoly_extrap <- subcrt("pyrrhotite", use.polymorphs = FALSE, exceed.Ttr = TRUE)
+expect_false(anyNA(sres_nopoly_extrap$out[[1]]$G))
+
# References
# Amend, J. P. and Shock, E. L. (2001)
Modified: pkg/CHNOSZ/man/subcrt.Rd
===================================================================
--- pkg/CHNOSZ/man/subcrt.Rd 2023-06-20 06:44:10 UTC (rev 789)
+++ pkg/CHNOSZ/man/subcrt.Rd 2023-06-20 11:47:58 UTC (rev 790)
@@ -11,7 +11,7 @@
property = c("logK","G","H","S","V","Cp"),
T = seq(273.15,623.15,25), P = "Psat", grid = NULL,
convert = TRUE, exceed.Ttr = FALSE, exceed.rhomin = FALSE,
- logact = NULL, autobalance = TRUE, IS = 0)
+ logact = NULL, autobalance = TRUE, use.polymorphs = TRUE, IS = 0)
}
\arguments{
@@ -27,6 +27,7 @@
\item{logact}{numeric, logarithms of activities of species in reaction}
\item{convert}{logical, are units of T, P, and energy settable by the user (default) (see \code{\link{T.units}})?}
\item{autobalance}{logical, attempt to automatically balance reaction with basis species?}
+ \item{use.polymorphs}{logical, automatically identify available polymorphs in OBIGT and use the stable one at each value of \code{T}?}
\item{IS}{numeric, ionic strength(s) at which to calculate adjusted molal properties, mol kg\eqn{^{-1}}{^-1}}
}
@@ -101,12 +102,12 @@
}
\section{Note on phase transitions}{
-Minerals with polymorphic transitions (denoted by having states \samp{cr} (lowest T phase), \samp{cr2}, \samp{cr3} etc.) can be defined generically by \samp{cr} in the \code{state} argument with a character value for \code{species}.
-\code{subcrt} in this case simultaneously calculates the requested properties of all the phases of each such mineral (phase species) and, using the values of the transition temperatures calculated from those at P = 1 bar given in the thermodynamic database together with functions of the entropies and volumes of transitions (see \code{\link{dPdTtr}}), determines the stable phase of the mineral at any grid point and substitutes the properties of this phase at that point for further calculations.
-If individual phase species of minerals are specified (by \samp{cr}, \samp{cr2} etc. in \code{state}), and \code{exceed.Ttr} is \code{FALSE} (the default), the Gibbs energies of these minerals are assigned values of NA at conditions beyond their transition temperature, where another of the phases is stable.
-If you set \code{exceed.Ttr} to \code{TRUE} to calculate the properties of mineral polymorphs (i.e., using \samp{cr}) the function will identify the stable polymorph using the calculated Gibbs energies of the phase species instead of the tabulated transition temperatures.
-This is not generally advised, as the computed standard molal properties of a phase species lose their physical meaning beyond the transition temperatures of the phase.
-
+Minerals with polymorphic transitions (denoted in OBIGT by having states \samp{cr} (lowest-\T phase), \samp{cr2}, etc.) can be specified by name with \samp{cr} for the \code{state} or by using a numeric species index for the lowest-\T polymorph.
+If \code{use.polymorphs} is TRUE, \code{subcrt} uses the transition temperatures calculated from those at P = 1 bar given in OBIGT together with functions of the entropies and volumes of transitions (see \code{\link{dPdTtr}}) to determine the stable polymorph at each grid point and uses the properties of that polymorph in the output.
+A \code{polymorph} column is added to the output to indicate the stable polymorph at each \T-\P condition.
+If \code{exceed.Ttr} is \code{FALSE} (the default), output values of Gibbs energy are assigned NA beyond the transition temperature of the highest-temperature polymorph.
+Set \code{exceed.Ttr} to \code{TRUE} to identify the stable polymorphs by comparing their extrapolated Gibbs energies instead of the tabulated transition temperatures.
+This is generally not advised, as extrapolated Gibbs energies may not reliably determine the stable polymorph at extreme temperatures.
}
\value{
@@ -150,12 +151,12 @@
subcrt(c("HS-", "SO4-2"), c(-1, 1), T = 100)
## Mineral polymorphs
-# Properties of the stable polymorph
+# Properties of the stable polymorph at each temperature
subcrt("pyrrhotite")
-# Properties of one of the polymorphs (metastable at other temperatures)
-subcrt(c("pyrrhotite"), state = "cr2")
# Reactions automatically use stable polymorph
subcrt(c("pyrite", "pyrrhotite", "H2O", "H2S", "O2"), c(-1, 1, -1, 1, 0.5))
+# Extrapolated properties of the lowest-T polymorph (metastable at higher temperatures)
+subcrt(c("pyrrhotite"), use.polymorphs = FALSE, exceed.Ttr = TRUE)
## Messages about problems with the calculation
# Above the T, P limits for the H2O equations of state
More information about the CHNOSZ-commits
mailing list