From 5bdffe1108e6b62c3d9a3158305606e0af642e3c Mon Sep 17 00:00:00 2001 From: Jeffrey Dick Date: Sun, 11 Feb 2024 14:48:58 +0000 Subject: [PATCH] Submit version 2.1.0 to CRAN git-svn-id: svn://scm.r-forge.r-project.org/svnroot/chnosz/pkg/CHNOSZ@827 edb9625f-4e0d-4859-8d74-9fd3b1da38cb --- DESCRIPTION | 4 ++-- R/cgl.R | 9 ++++---- R/util.misc.R | 5 +++-- demo/aluminum.R | 2 +- inst/CHECKLIST | 9 +------- inst/NEWS.Rd | 16 +++++++++------ inst/tinytest/test-AD.R | 4 ++-- inst/tinytest/test-subcrt.R | 41 +++++++++++++++++++------------------ inst/tinytest/test-util.R | 9 ++++++++ vignettes/FAQ.Rmd | 3 ++- vignettes/OBIGT.bib | 2 +- vignettes/vig.bib | 11 ++++++++++ 12 files changed, 67 insertions(+), 48 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 66bd3ed..8588de0 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ -Date: 2024-02-10 +Date: 2024-02-11 Package: CHNOSZ -Version: 2.0.0-46 +Version: 2.1.0 Title: Thermodynamic Calculations and Diagrams for Geochemistry Authors@R: c( person("Jeffrey", "Dick", , "j3ffdick@gmail.com", role = c("aut", "cre"), diff --git a/R/cgl.R b/R/cgl.R index 8223601..0273a61 100644 --- a/R/cgl.R +++ b/R/cgl.R @@ -111,7 +111,7 @@ quartz_coesite <- function(PAR, T, P) { ca <- -0.4973e-4 VPtTta <- 23.348 VPrTtb <- 23.72 - Stran <- convert(0.342, "J") + Stran <- 0.342 # Constants from REAC92D.f VPrTra <- 22.688 # VPrTr(a-quartz) Vdiff <- 2.047 # VPrTr(a-quartz) - VPrTr(coesite) @@ -136,13 +136,12 @@ quartz_coesite <- function(PAR, T, P) { V <- V - Vdiff } cm3bar_to_cal <- 0.023901488 - cm3bar_to_J <- convert(cm3bar_to_cal, "J") - GVterm <- cm3bar_to_J * (VPrTra * (P - Pstar) + VPrTtb * (Pstar - Pr) - + GVterm <- cm3bar_to_cal * (VPrTra * (P - Pstar) + VPrTtb * (Pstar - Pr) - 0.5 * ca * (2 * Pr * (P - Pstar) - (P^2 - Pstar^2)) - ca * k * (T - Tr) * (P - Pstar) + k * (ba + aa * ca * k) * (T - Tr) * log((aa + P/k) / (aa + Pstar/k))) - SVterm <- cm3bar_to_J * (-k * (ba + aa * ca * k) * + SVterm <- cm3bar_to_cal * (-k * (ba + aa * ca * k) * log((aa + P/k) / (aa + Pstar/k)) + ca * k * (P - Pstar)) - Sstar # Note the minus sign on "SVterm" in order that intdVdTdP has the correct sign - list(intVdP = GVterm, intdVdTdP = -SVterm, V = V) + list(intVdP = convert(GVterm, "J"), intdVdTdP = convert(-SVterm, "J"), V = V) } diff --git a/R/util.misc.R b/R/util.misc.R index 284b60b..342d16c 100644 --- a/R/util.misc.R +++ b/R/util.misc.R @@ -7,7 +7,8 @@ dPdTtr <- function(ispecies, ispecies2 = NULL) { # (argument is index of the lower-T phase) thermo <- get("thermo", CHNOSZ) if(is.null(ispecies2)) ispecies2 <- ispecies + 1 - pars <- info(c(ispecies, ispecies2), check.it = FALSE) + # Use OBIGT2eos to ensure that all parameters are in Joules 20240211 + pars <- OBIGT2eos(thermo$OBIGT[c(ispecies, ispecies2), ], "cgl", fixGHS = TRUE, toJoules = TRUE) # If these aren't the same mineral, we shouldn't be here if(as.character(pars$name[1]) != as.character(pars$name[2])) stop("different names for species ", ispecies, " and ", ispecies2) # The special handling for quartz and coesite interfere with this function, @@ -15,7 +16,7 @@ dPdTtr <- function(ispecies, ispecies2 = NULL) { pars$name <- toupper(pars$name) props <- cgl(c("G", "S", "V"), pars, P = 0, T = thermo$OBIGT$z.T[ispecies]) # The G's should be the same ... - #if(abs(props[[2]]$G - props[[1]]$G) > 0.1) warning('dP.dT: inconsistent values of G for different phases of ',ispecies,call. = FALSE) + #if(abs(props[[2]]$G - props[[1]]$G) > 0.1) warning('dP.dT: inconsistent values of G for different polymorphs of ',ispecies,call. = FALSE) dP.dT <- convert( ( props[[2]]$S - props[[1]]$S ) / ( props[[2]]$V - props[[1]]$V ), 'cm3bar' ) return(dP.dT) } diff --git a/demo/aluminum.R b/demo/aluminum.R index 12d7860..e4b3a96 100644 --- a/demo/aluminum.R +++ b/demo/aluminum.R @@ -55,7 +55,7 @@ legend("topleft", c("Boehmite - Kaolinite", "After Zhu and Lu, 2009 Fig. A1"), b # Shock et al., 1989 (SHS89): doi:10.1016/0016-7037(89)90341-4 # Berman, 1988 (Ber88): doi:10.1093/petrology/29.2.445 # Holland and Powell, 2011 (HP11): 10.1111/j.1525-1314.2010.00923.x -# Hemingway et al., 1991 (HRA91): http://pubs.er.usgs.gov/publication/70016664 +# Hemingway et al., 1991 (HRA91): https://pubs.usgs.gov/publication/70016664 # Apps and Spycher, 2004 (AS04): Bechtel SAIC Company, LLC ANL-NBS-HS-000043 REV 00 (DOC.20041118.0004) ########### diff --git a/inst/CHECKLIST b/inst/CHECKLIST index 8b848f1..f1e8b63 100644 --- a/inst/CHECKLIST +++ b/inst/CHECKLIST @@ -1,6 +1,6 @@ **************************** Release checklist for CHNOSZ - (updated 2023-11-29) + (updated 2024-02-11) **************************** - Run examples() and demos() and inspect their output (especially plots) @@ -77,10 +77,3 @@ Making vignettes for website (https://chnosz.net) - Install the package and run doc/mklinks.sh within the installation directory. (this adds links to the HTML renditions of Rd files) - -*************** -Temporary Items -*************** - -- Replace reference to "development version of CHNOSZ (to be 2.0.1)" in FAQ.Rmd. - diff --git a/inst/NEWS.Rd b/inst/NEWS.Rd index bcb2e3a..9f1550a 100644 --- a/inst/NEWS.Rd +++ b/inst/NEWS.Rd @@ -15,7 +15,7 @@ \newcommand{\Cp}{\ifelse{latex}{\eqn{C_P}}{\ifelse{html}{\out{CP}}{Cp}}} \newcommand{\DG0}{\ifelse{latex}{\eqn{{\Delta}G^{\circ}}}{\ifelse{html}{\out{ΔG°}}{ΔG°}}} -\section{Changes in CHNOSZ version 2.0.0-44 (2023-02-06)}{ +\section{Changes in CHNOSZ version 2.1.0 (2024-02-11)}{ \subsection{NEW FEATURES}{ \itemize{ @@ -40,15 +40,14 @@ \item \code{mosaic()} now handles the situation where the alternate basis species include one of the axis variables on a diagram (by changing the - argument names in its internal calls to \code{affinity()}) and - adjusts the labels for the diagram accordingly (e.g. \dQuote{total C}, - \dQuote{total S}, etc.). Thanks to Evgeniy Bastrakov for the feature - request. + argument names in its internal calls to \code{affinity()}) and adjusts + the labels for the diagram accordingly (e.g. \dQuote{total C}, + \dQuote{total S}, etc.). Thanks to Evgeniy Bastrakov for the suggestion. \item The environment variable CHNOSZ_BUILD_LARGE_VIGNETTES is used to control \strong{dpi} in \pkg{knitr} chunk options. Setting this variable results in larger vignettes that are used for the CHNOSZ website; if this - variable is unset (as on CRAN), a smaller package is built. + variable is unset (as in CRAN checks), a smaller package is built. \item \code{rank.affinity()} now returns average group rankings as percentage values. @@ -101,6 +100,11 @@ \strong{IS} arguments were not applied to automatically balanced reactions. + \item Fix bug in \code{dPdTtr()} where OBIGT database parameters were not + converted to Joules. This allows re-activating the tests for quartz + (which depend on the computed transition temperature) in + \file{test-subcrt.R}. + } } diff --git a/inst/tinytest/test-AD.R b/inst/tinytest/test-AD.R index c0560fc..096a7eb 100644 --- a/inst/tinytest/test-AD.R +++ b/inst/tinytest/test-AD.R @@ -81,8 +81,8 @@ expect_equal(GAD, GOBIGT, tolerance = 280, scale = 1, info = info) # The largest differences are for HCl, ethane, and B(OH)3 expect_equal(sort(info(iaq[abs(GAD - GOBIGT) > 900])$name), sort(c("HCl", "ethane", "B(OH)3"))) -# This line should be commented for a released package -exit_file("Skipping tests so development builds on R-Forge work") +## This line should be commented for a released package +#exit_file("Skipping tests so development builds on R-Forge work") ## The following tests work on JMD's Linux machine "at home" but not on some CRAN machines 20220210 if(!at_home()) exit_file("Skipping tests on CRAN") diff --git a/inst/tinytest/test-subcrt.R b/inst/tinytest/test-subcrt.R index d9e2f02..03c4df8 100644 --- a/inst/tinytest/test-subcrt.R +++ b/inst/tinytest/test-subcrt.R @@ -34,16 +34,6 @@ expect_true(maxdiff(sout.C7$G/1000, AS01.C7) < 0.06, info = info) # We can also check that sulfur has expected polymorphic transitions expect_equal(s.C7$polymorphs$sulfur, c(1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 3, 3), info = info) -info <- "Subzero degree C calculations are possible" -## Start with H2O -s.H2O <- subcrt("H2O", T = c(-20.1, seq(-20, 0)), P = 1)$out$water -# We shouldn't get anything at -20.1 deg C -expect_true(is.na(s.H2O$G[1]), info = info) -# We should get something at -20 deg C -expect_equal(s.H2O$G[2], convert(-56001, "J"), tolerance = 1, scale = 1, info = info) -# Following SUPCRT92, an input temperature of 0 is converted to 0.01 -expect_equal(s.H2O$T[22], 0.01, info = info) - info <- "Calculations using IAPWS-95 are possible" oldwat <- water("IAPWS95") # The tests pass if we get numeric values and no error @@ -81,9 +71,7 @@ expect_equal(round(CHNOSZ$V, 1), SUPCRT_V, info = info) expect_equal(round(CHNOSZ$Cp, 1), SUPCRT_Cp, info = info) OBIGT() -# TODO: fix quartz_coesite() for switch to Joules 20220325 -if(FALSE) { - +# Quartz tests re-activated after fixing dPdTtr to use parameters converted to Joules 20240211 info <- "Calculations for quartz are nearly consistent with SUPCRT92" add.OBIGT("SUPCRT92") # Using SUPCRT's equations, the alpha-beta transition occurs at @@ -97,7 +85,7 @@ SUPCRT_S <- c(12.3, 10.6, 31.8, 29.8) SUPCRT_V <- c(22.5, 20.3, 23.7, 21.9) SUPCRT_Cp <- c(12.3, 12.3, 16.9, 16.9) CHNOSZ <- subcrt("quartz", T = T, P = P)$out[[1]] -# NOTE: Testing has shown that, where alpha-quartz is stable above Ttr(Pr) but below Ttr(P), +# NOTE: It appears that, where alpha-quartz is stable above Ttr(Pr) but below Ttr(P), # SUPCRT92 computes the heat capacity and its integrals using parameters of beta-quartz. # (see e.g. the equation for CprdT under the (Cpreg .EQ. 2) case in the Cptrms subroutine of SUPCRT). # ... is that incorrect? @@ -138,14 +126,12 @@ S92_5000bar <- read.table(header = TRUE, text = " ") CHNOSZ_5000bar <- subcrt("quartz", T = seq(703, 706), P = 5000)$out[[1]] # NOTE: calculated values *below* the transition are different -expect_true(maxdiff(CHNOSZ_5000bar$G, S92_5000bar$G) < 20, info = info) -expect_true(maxdiff(CHNOSZ_5000bar$H, S92_5000bar$H) < 300, info = info) -expect_true(maxdiff(CHNOSZ_5000bar$S, S92_5000bar$S) < 0.5, info = info) -expect_true(maxdiff(CHNOSZ_5000bar$V, S92_5000bar$V) < 0.05, info = info) +expect_equal(CHNOSZ_5000bar$G, S92_5000bar$G, tolerance = 20, scale = 1, info = info) +expect_equal(CHNOSZ_5000bar$H, S92_5000bar$H, tolerance = 300, scale = 1, info = info) +expect_equal(CHNOSZ_5000bar$S, S92_5000bar$S, tolerance = 0.5, scale = 1, info = info) +expect_equal(CHNOSZ_5000bar$V, S92_5000bar$V, tolerance = 0.05, scale = 1, info = info) OBIGT() -} # end if(FALSE) - info <- "Duplicated species yield correct polymorphic transitions" # If a mineral with polymorphic transitions is in both the basis and species lists, # energy()'s call to subcrt() will have duplicated species. @@ -260,6 +246,21 @@ T <- c(25, 50, 103, 104) sout <- subcrt("carrollite", T = T, P = 1)$out[[1]] expect_equal(sout$polymorph, c(1, 1, 1, 2), info = info) +info <- "Subzero degree C calculations are possible" +# Set default units +E.units("J") +# Start with H2O +s.H2O <- subcrt("H2O", T = c(-20.1, seq(-20, 0)), P = 1)$out$water +# We should get something at -20 deg C +expect_equal(s.H2O$G[2], convert(-56001, "J"), tolerance = 1, scale = 1, info = info) +# Following SUPCRT92, an input temperature of 0 is converted to 0.01 +expect_equal(s.H2O$T[22], 0.01, info = info) +# The following test fails CRAN checks with Intel oneAPI 2023.x compilers +# (Expected TRUE, got FALSE) 20240211 +if(!at_home()) exit_file("Skipping tests on CRAN") +# We shouldn't get anything at -20.1 deg C +expect_true(is.na(s.H2O$G[1]), info = info) + # References # Amend, J. P. and Shock, E. L. (2001) diff --git a/inst/tinytest/test-util.R b/inst/tinytest/test-util.R index b98a532..334adde 100644 --- a/inst/tinytest/test-util.R +++ b/inst/tinytest/test-util.R @@ -41,3 +41,12 @@ expect_error(E.units("X"), "units of energy must be either cal or J", info = inf info <- "describe.property() does not accept NULL values" expect_error(describe.property(), "property or value is NULL", info = info) + +# Test added on 20240211 +# Incorrect result was causing quartz tests in test-subcrt.R to fail +info <- "dPdTtr gives expected result" +add.OBIGT("SUPCRT92") +dPdTtr.calc <- round(dPdTtr(info("quartz", "cr"), info("quartz", "cr2")), 5) +# The reference value was calculated with CHNOSZ_1.4.3 +# (prior to bug in dPdTtr introduced by switch to Joules) +expect_equal(dPdTtr.calc, 38.45834, info = info) diff --git a/vignettes/FAQ.Rmd b/vignettes/FAQ.Rmd index 5eaadb2..8f90293 100644 --- a/vignettes/FAQ.Rmd +++ b/vignettes/FAQ.Rmd @@ -747,7 +747,7 @@ reset() The warning is similar to that produced by SUPCRT92 ("CAUTION: BEYOND T LIMIT OF CP COEFFS FOR A MINERAL OR GAS") at temperatures above maximum temperature of validity of the Maier-Kelley equation (Tmax). Notably, SUPCRT92 outputs Δ*G*° and other standard thermodynamic properties at temperatures higher than Tmax despite the warning. -This is a new feature in the development version of CHNOSZ (to be 2.0.1). +This is a new feature in CHNOSZ version 2.1.0. In previous versions of CHNOSZ, values of Δ*G*° above the *C~p~* equation limit were set to NA without a warning, as with the phase stability limit described above. **4. Finally, if `T` is NA or 0, then no upper temerature limit is imposed by `subcrt()`.** @@ -820,6 +820,7 @@ The parameters for these minerals in the default OBIGT database come from @Ber88 thermo.refs(species()$ispecies) ``` +CHNOSZ doesn't implement the thermodynamic model for minerals from @HP98, which is one of the sources cited by @HC14. If we use the thermodynamic parameters for minerals from @HDNB78 [these do not include the revisions for aluminosilicates described by @SHD91], we get the lines shown in the second plot above, representing a larger stability field for muscovite. This moves the KMQ buffer closer to the value shown by @HC14, but the MC buffer further away, so this still doesn't explain why we get a different result. diff --git a/vignettes/OBIGT.bib b/vignettes/OBIGT.bib index e19c58c..710e6cb 100644 --- a/vignettes/OBIGT.bib +++ b/vignettes/OBIGT.bib @@ -1274,7 +1274,7 @@ @Article{HRA91 number = {3-4}, pages = {445--457}, volume = {76}, - url = {https://pubs.er.usgs.gov/publication/70016664}, + url = {https://pubs.usgs.gov/publication/70016664}, } @Article{BPAH07, diff --git a/vignettes/vig.bib b/vignettes/vig.bib index 3bbd6e3..2e62d95 100644 --- a/vignettes/vig.bib +++ b/vignettes/vig.bib @@ -886,3 +886,14 @@ @Article{Yar05 volume = {100}, doi = {10.2113/gsecongeo.100.4.613}, } + +@Article{HP98, + author = {Holland, T. J. B. and Powell, R.}, + journal = {Journal of Metamorphic Geology}, + title = {An internally consistent thermodynamic data set for phases of petrological interest}, + year = {1998}, + number = {3}, + pages = {309--343}, + volume = {16}, + doi = {10.1111/j.1525-1314.1998.00140.x}, +}