OBIGT: Use model = "CGL_Ttr" for species with a phase stability limit

git-svn-id: svn://scm.r-forge.r-project.org/svnroot/chnosz/pkg/CHNOSZ@864 edb9625f-4e0d-4859-8d74-9fd3b1da38cb

DESCRIPTION

@@ -1,6 +1,6 @@
-Date: 2024-12-28
+Date: 2024-12-29
 Package: CHNOSZ
-Version: 2.1.0-35
+Version: 2.1.0-36
 Title: Thermodynamic Calculations and Diagrams for Geochemistry
 Authors@R: c(
     person("Jeffrey", "Dick", , "[email protected]", role = c("aut", "cre"),

R/subcrt.R

@@ -191,8 +191,8 @@ subcrt <- function(species, coeff = 1, state = NULL, property = c("logK", "G", "
     ispecies = iphases, model = model, stringsAsFactors = FALSE)
   # Make the rownames readable ... but they have to be unique
   if(length(unique(iphases))==length(iphases)) rownames(reaction) <- as.character(iphases)
-  # Which species use models for aqueous species
-  # This breaks things if we have state = "aq" and model = "CGL" 20230220
+  # Identify aqueous species from their 'model' parameter
+  # NOTE: checking for state == "aq" would break using CGL for aqueous species 20230220
   #isaq <- reaction$state == "aq"
   isaq <- toupper(reaction$model) %in% c("HKF", "AD", "DEW")
 
@@ -334,7 +334,7 @@ subcrt <- function(species, coeff = 1, state = NULL, property = c("logK", "G", "
   }
 
   # Crystalline, gas, or liquid (except water) species
-  iscgl <- reaction$model %in% c("CGL", "Berman")
+  iscgl <- reaction$model %in% c("CGL", "CGL_Ttr", "Berman")
 
   if(TRUE %in% iscgl) {
     param <- OBIGT2eos(thermo$OBIGT[iphases[iscgl],], "cgl", fixGHS = TRUE, toJoules = TRUE)
@@ -387,26 +387,23 @@ subcrt <- function(species, coeff = 1, state = NULL, property = c("logK", "G", "
         }
 
         if(!is.polymorphic.transition) {
-          # Check if we're above the T limit for a Cp equation or T for phase change (e.g. melting or vaporization) listed in OBIGT
-          is.Cp.equation <- FALSE
+          # Check if we're above the T limit for a Cp equation or a phase change (e.g. melting or vaporization)
           if(all(is.na(Ttr))) next
           if(all(Ttr == 0)) next
-          ineg <- Ttr < 0
-          if(any(ineg)) {
-            # This is the T limit for a Cp equation (not a phase change)
-            Ttr[ineg] <- -Ttr[ineg]
-            is.Cp.equation <- TRUE
-          }
+          is.phase_change <- any(thermo$OBIGT$model[iphases[i]] == "CGL_Ttr")
           if(any(T > Ttr)) {
-            if(!exceed.Ttr) {
-                if(is.Cp.equation) {
-                  warning(paste0("above T limit of ", Ttr, " K for the Cp equation for ", myname, "(", mystate, ")"))
-                } else {
-                  message(paste0("subcrt: G is set to NA for ", myname, "(", mystate, ") above its stability limit of ", Ttr, " K (use exceed.Ttr = TRUE to output G)"))
-                  p.cgl[[ncgl[i]]]$G[T > Ttr] <- NA
-                }
+            if(is.phase_change) {
+              if(exceed.Ttr) {
+                message(paste0("subcrt: showing G for ", myname, "(", mystate, ") above its stability limit of ", Ttr, " K (use exceed.Ttr = FALSE to prevent this)"))
+              } else {
+                message(paste0("subcrt: setting G to NA for ", myname, "(", mystate, ") above its stability limit of ", Ttr, " K (use exceed.Ttr = TRUE to output G)"))
+                p.cgl[[ncgl[i]]]$G[T > Ttr] <- NA
+              }
             } else {
-              message(paste0("subcrt: showing G for ", myname, "(", mystate, ") above its stability limit of ", Ttr, " K (use exceed.Ttr = FALSE to prevent this)"))
+              if(! exceed.Ttr) {
+                # Warn if we're above a Cp limit (but don't change the output)
+                warning(paste0("above T limit of ", Ttr, " K for the Cp equation for ", myname, "(", mystate, ")"))
+              }
             }
           }
         }

R/util.data.R

@@ -413,7 +413,7 @@ OBIGT2eos <- function(OBIGT, state, fixGHS = FALSE, toJoules = FALSE) {
   # Figure out the model for each species 20220929
   model <- OBIGT$model
   model[is.na(model)] <- ""
-  isCGL <- model == "CGL"
+  isCGL <- model %in% c("CGL", "CGL_Ttr")
   isHKF <- model == "HKF"
   isDEW <- model == "DEW"
   isAD <- model == "AD"

inst/NEWS.Rd

@@ -15,7 +15,7 @@
 \newcommand{\Cp}{\ifelse{latex}{\eqn{C_P}}{\ifelse{html}{\out{<I>C<sub>P</sub></I>}}{Cp}}}
 \newcommand{\DG0}{\ifelse{latex}{\eqn{{\Delta}G^{\circ}}}{\ifelse{html}{\out{&Delta;<I>G</I>&deg;}}{ΔG°}}}
 
-\section{Changes in CHNOSZ version 2.1.0-35 (2024-12-28)}{
+\section{Changes in CHNOSZ version 2.1.0-36 (2024-12-29)}{
 
   \subsection{OBIGT DEFAULT DATA}{
     \itemize{
@@ -80,6 +80,11 @@
   \subsection{OTHER CHANGES}{
     \itemize{
 
+      \item \code{OBIGT} now lists \samp{CGL_Ttr} in \code{model} for species
+      whose \code{T} value denotes a phase stability limit (as opposed to a Cp
+      equation limit).  This replaces the (overloaded) use of signed values of
+      \code{T}.
+
       \item \code{check.EOS()} now uses values of Born coefficients \emph{X}
       and \emph{Q} that are consistent with either SUPCRT92 or DEW, depending
       on the \code{model} defined in OBIGT (HKF or DEW, respectively).

inst/extdata/OBIGT/AD.csv

@@ -20,6 +20,6 @@ SO2,NA,SO2,aq,AD03.2,NA,2019-02-21,AD,J,-71980,NA,38.7,NA,38.5,-14.5223,14.3512,
 B(OH)3,NA,B(OH)3,aq,AP14,NA,2019-02-22,AD,J,NA,NA,NA,NA,NA,-4.2561,4.0194,-1.057,NA,NA,NA,NA,NA
 Si(OH)4,NA,Si(OH)4,aq,AP14,NA,2019-02-22,AD,J,NA,NA,NA,NA,NA,0.9285,-0.9409,-1.8933,NA,NA,NA,NA,NA
 As(OH)3,NA,As(OH)3,aq,AP14,NA,2019-02-22,AD,J,NA,NA,NA,NA,NA,-9.903,7.6818,-1.23,NA,NA,NA,NA,NA
-B(OH)3,NA,B(OH)3,gas,AP14,NA,2019-02-22,CGL,cal,-222523.9,NA,66.589,NA,0,46.7734,-9.7E-05,59320,-522.4665,0,0,0,-1000
-Si(OH)4,NA,Si(OH)4,gas,AP14,NA,2019-02-22,CGL,cal,-296285.9,NA,82.784,NA,0,22.9685,0.0376434,-444070,0,-2.7782E-05,7.923E-09,3,-1000
-As(OH)3,NA,As(OH)3,gas,AP14,NA,2019-02-22,CGL,cal,-143112.8,NA,73.599,NA,0,22.6338,-0.0163623,-372130,0,-1.0194E-05,2.7414E-09,3,-1000
+B(OH)3,NA,B(OH)3,gas,AP14,NA,2019-02-22,CGL,cal,-222523.9,NA,66.589,NA,0,46.7734,-9.7E-05,59320,-522.4665,0,0,0,1000
+Si(OH)4,NA,Si(OH)4,gas,AP14,NA,2019-02-22,CGL,cal,-296285.9,NA,82.784,NA,0,22.9685,0.0376434,-444070,0,-2.7782E-05,7.923E-09,3,1000
+As(OH)3,NA,As(OH)3,gas,AP14,NA,2019-02-22,CGL,cal,-143112.8,NA,73.599,NA,0,22.6338,-0.0163623,-372130,0,-1.0194E-05,2.7414E-09,3,1000