Add 'use.polymorphs' argument to subcrt()

git-svn-id: svn://scm.r-forge.r-project.org/svnroot/chnosz/pkg/CHNOSZ@790 edb9625f-4e0d-4859-8d74-9fd3b1da38cb
jedick · Jun 20, 2023 · 4895ad2 · 4895ad2
1 parent 3eafb75
commit 4895ad2
Show file tree

Hide file tree

Showing 5 changed files with 65 additions and 43 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -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@R: c(
     person("Jeffrey", "Dick", , "[email protected]", role = c("aut", "cre"),

diff --git a/R/subcrt.R b/R/subcrt.R
@@ -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,37 +118,36 @@ subcrt <- function(species, coeff = 1, state = NULL, property = c("logK", "G", "
 
   # 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])
     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 @@ subcrt <- function(species, coeff = 1, state = NULL, property = c("logK", "G", "
   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 @@ subcrt <- function(species, coeff = 1, state = NULL, property = c("logK", "G", "
         # 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,15 +444,15 @@ subcrt <- function(species, coeff = 1, state = NULL, property = c("logK", "G", "
       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
         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 @@ subcrt <- function(species, coeff = 1, state = NULL, property = c("logK", "G", "
         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)

diff --git a/inst/NEWS.Rd b/inst/NEWS.Rd
@@ -12,11 +12,14 @@
 % 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.
 
     }
 

diff --git a/inst/tinytest/test-subcrt.R b/inst/tinytest/test-subcrt.R
@@ -202,6 +202,23 @@ expect_error(subcrt("H2O", -1, "liq", "xxx"), "invalid property name: xxx", info
 # 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) 

diff --git a/man/subcrt.Rd b/man/subcrt.Rd
@@ -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 @@ On the other hand, if the \code{water} option is \samp{SUPCRT} (the default), \c
 }
 
 \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 @@ basis(c("H2O", "H2S", "O2", "H+"))
 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