[CHNOSZ-commits] r792 - in pkg/CHNOSZ: . R inst inst/tinytest man vignettes
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Jun 23 16:00:30 CEST 2023
Author: jedick
Date: 2023-06-23 16:00:30 +0200 (Fri, 23 Jun 2023)
New Revision: 792
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/NAMESPACE
pkg/CHNOSZ/R/info.R
pkg/CHNOSZ/R/logB.to.OBIGT.R
pkg/CHNOSZ/R/subcrt.R
pkg/CHNOSZ/R/util.data.R
pkg/CHNOSZ/inst/NEWS.Rd
pkg/CHNOSZ/inst/tinytest/test-subcrt.R
pkg/CHNOSZ/inst/tinytest/test-util.data.R
pkg/CHNOSZ/man/extdata.Rd
pkg/CHNOSZ/man/info.Rd
pkg/CHNOSZ/man/thermo.Rd
pkg/CHNOSZ/man/util.data.Rd
pkg/CHNOSZ/vignettes/FAQ.Rmd
pkg/CHNOSZ/vignettes/anintro.Rmd
pkg/CHNOSZ/vignettes/custom_data.Rmd
pkg/CHNOSZ/vignettes/mklinks.sh
pkg/CHNOSZ/vignettes/vig.bib
Log:
Add FAQ for polymorphs and rename functions to check.GHS() and check.EOS()
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2023-06-21 06:24:29 UTC (rev 791)
+++ pkg/CHNOSZ/DESCRIPTION 2023-06-23 14:00:30 UTC (rev 792)
@@ -1,6 +1,6 @@
-Date: 2023-06-21
+Date: 2023-06-23
Package: CHNOSZ
-Version: 2.0.0-11
+Version: 2.0.0-12
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/NAMESPACE
===================================================================
--- pkg/CHNOSZ/NAMESPACE 2023-06-21 06:24:29 UTC (rev 791)
+++ pkg/CHNOSZ/NAMESPACE 2023-06-23 14:00:30 UTC (rev 792)
@@ -39,7 +39,7 @@
"count.elements",
# (no other functions are used in the tests)
# other exported functions that are not used above
- "checkEOS", "checkGHS", "check.OBIGT",
+ "check.EOS", "check.GHS", "check.OBIGT",
"V_s_var", "Cp_s_var",
"EOSlab", "EOScalc",
"basis.elements", "element.mu", "ibasis",
Modified: pkg/CHNOSZ/R/info.R
===================================================================
--- pkg/CHNOSZ/R/info.R 2023-06-21 06:24:29 UTC (rev 791)
+++ pkg/CHNOSZ/R/info.R 2023-06-23 14:00:30 UTC (rev 792)
@@ -231,9 +231,9 @@
# Don't do it for the AD species 20190219
if(check.it & this$model != "AD") {
# Check GHS if they are all present
- if(sum(naGHS)==0) calcG <- checkGHS(this)
- # Check tabulated heat capacities against EOS parameters
- calcCp <- checkEOS(this, this$model, "Cp")
+ if(sum(naGHS)==0) calcG <- check.GHS(this, ret.diff = FALSE)
+ # Check heat capacities in database against EOS parameters
+ calcCp <- check.EOS(this, this$model, "Cp", ret.diff = FALSE)
# Fill in NA heat capacity
if(!is.na(calcCp) & is.na(this$Cp)) {
message("info.numeric: Cp of ", this$name, "(", this$state, ") is NA; set by EOS parameters to ", round(calcCp, 2), " ", this$E_units, " K-1 mol-1")
@@ -243,9 +243,9 @@
calcCp <- Berman(this$name)$Cp
this$Cp <- calcCp
}
- # Check tabulated volumes - only for HKF model (aq species)
+ # Check volumes in database - only for HKF model (aq species)
if(this$model %in% c("HKF", "DEW")) {
- calcV <- checkEOS(this, this$model, "V")
+ calcV <- check.EOS(this, this$model, "V", ret.diff = FALSE)
# Fill in NA volume
if(!is.na(calcV) & is.na(this$V)) {
message("info.numeric: V of ", this$name, "(", this$state, ") is NA; set by EOS parameters to ", round(calcV, 2), " cm3 mol-1")
Modified: pkg/CHNOSZ/R/logB.to.OBIGT.R
===================================================================
--- pkg/CHNOSZ/R/logB.to.OBIGT.R 2023-06-21 06:24:29 UTC (rev 791)
+++ pkg/CHNOSZ/R/logB.to.OBIGT.R 2023-06-23 14:00:30 UTC (rev 792)
@@ -120,7 +120,7 @@
if(is.na(c1)) c1 <- 0
if(is.na(c2)) c2 <- 0
if(is.na(omega)) omega <- 0
- # Calculate Cp at 25 °C (not used in HKF - just for info() and checkEOS())
+ # Calculate Cp at 25 °C (not used in HKF - just for info() and check.EOS())
PAR$c1 <- c1
PAR$c2 <- c2
PAR$omega <- omega
Modified: pkg/CHNOSZ/R/subcrt.R
===================================================================
--- pkg/CHNOSZ/R/subcrt.R 2023-06-21 06:24:29 UTC (rev 791)
+++ pkg/CHNOSZ/R/subcrt.R 2023-06-23 14:00:30 UTC (rev 792)
@@ -20,22 +20,19 @@
# Revise the call if the states are the second argument
if(!is.null(coeff[1])) {
- if(is.numeric(state[1])) newcoeff <- state else newcoeff <- 1
- if(is.character(coeff[1])) newstate <- coeff else newstate <- NULL
if(is.character(coeff[1])) {
- if(missing(T)) {
- if(identical(newcoeff, 1) & !(identical(newcoeff, state)))
- return(subcrt(species, state = coeff, property = property, P = P, grid = grid,
- convert = convert, exceed.Ttr = exceed.Ttr, logact = logact))
- else return(subcrt(species, coeff = newcoeff, state = coeff, property = property,
- P = P, grid = grid, convert = convert, exceed.Ttr = exceed.Ttr, logact = logact))
- } else {
- if(identical(newcoeff,1) & !(identical(newcoeff, state)))
- return(subcrt(species, state = coeff, property = property, T = T, P = P, grid = grid,
- convert = convert, exceed.Ttr = exceed.Ttr, logact = logact))
- else return(subcrt(species, coeff = newcoeff, state = coeff, property = property,
- T = T, P = P, grid = grid, convert = convert, exceed.Ttr = exceed.Ttr, logact = logact))
+ newstate <- coeff
+ # This is missing coeff and T in order that missing values are correctly detected further below 20230621
+ newargs <- list(species = species, state = newstate,
+ property = property, P = P, grid = grid, convert = convert,
+ exceed.Ttr = exceed.Ttr, exceed.rhomin = exceed.rhomin, logact = logact,
+ autobalance = autobalance, use.polymorphs = use.polymorphs, IS = IS)
+ if(!missing(state)) {
+ if(is.numeric(state[1])) newcoeff <- state else stop("If they are both given, one of arguments 2 and 3 should be numeric reaction coefficients")
+ newargs <- c(list(coeff = newcoeff), newargs)
}
+ if(!missing(T)) newargs <- c(list(T = T), newargs)
+ return(do.call(subcrt, newargs))
}
}
Modified: pkg/CHNOSZ/R/util.data.R
===================================================================
--- pkg/CHNOSZ/R/util.data.R 2023-06-21 06:24:29 UTC (rev 791)
+++ pkg/CHNOSZ/R/util.data.R 2023-06-23 14:00:30 UTC (rev 792)
@@ -159,8 +159,8 @@
}
}
-checkEOS <- function(eos, model, prop, ret.diff = FALSE) {
- # Compare calculated properties from thermodynamic parameters with reference (tabulated) values.
+check.EOS <- function(eos, model, prop, ret.diff = TRUE) {
+ # Compare calculated properties from thermodynamic parameters with given (database) values.
# Print message and return the calculated value if tolerance is exceeded
# or NA if the difference is within the tolerance.
# 20110808 jmd
@@ -209,8 +209,8 @@
if(!is.na(calcval)) {
if(!is.na(refval)) {
if(abs(diff) > tol) {
- message(paste("checkEOS: ", prop, " of ", eos$name, "(", eos$state,
- ") differs by ", round(diff,2), " ", units, " from tabulated value", sep=""))
+ message(paste("check.EOS: calculated ", prop, " of ", eos$name, "(", eos$state,
+ ") differs by ", round(diff,2), " ", units, " from database value", sep=""))
return(calcval)
}
} else return(calcval)
@@ -220,8 +220,8 @@
return(NA)
}
-checkGHS <- function(ghs, ret.diff=FALSE) {
- # Compare calculated G from H and S with reference (tabulated) values
+check.GHS <- function(ghs, ret.diff = TRUE) {
+ # Compare G calculated from H and S with given (database) values
# Print message and return the calculated value if tolerance is exceeded
# or NA if the difference is within the tolerance
# 20110808 jmd
@@ -243,8 +243,8 @@
if(!is.na(refval)) {
diff <- calcval - refval
if(abs(diff) > thermo$opt$G.tol) {
- message(paste("checkGHS: G of ", ghs$name, "(", ghs$state,
- ") differs by ", round(diff), " ", ghs$E_units, " mol-1 from tabulated value", sep=""))
+ message(paste("check.GHS: calculated G of ", ghs$name, "(", ghs$state,
+ ") differs by ", round(diff), " ", ghs$E_units, " mol-1 from database value", sep=""))
return(calcval)
}
} else return(calcval)
@@ -278,10 +278,10 @@
isHKF <- tdata$model %in% c("HKF", "DEW")
if(any(isHKF)) {
eos.HKF <- OBIGT2eos(tdata[isHKF,], "HKF")
- DCp.HKF <- checkEOS(eos.HKF, "HKF", "Cp", ret.diff = TRUE)
- DV.HKF <- checkEOS(eos.HKF, "HKF", "V", ret.diff = TRUE)
+ DCp.HKF <- check.EOS(eos.HKF, "HKF", "Cp")
+ DV.HKF <- check.EOS(eos.HKF, "HKF", "V")
cat(paste("check.OBIGT: GHS for", sum(isHKF), "species with HKF model in", what, "\n"))
- DG.HKF <- checkGHS(eos.HKF, ret.diff = TRUE)
+ DG.HKF <- check.GHS(eos.HKF)
# Store the results
DCp[isHKF] <- DCp.HKF
DV[isHKF] <- DV.HKF
@@ -290,9 +290,9 @@
# Then other species, if they are present
if(sum(!isHKF) > 0) {
eos.cgl <- OBIGT2eos(tdata[!isHKF,], "cgl")
- DCp.cgl <- checkEOS(eos.cgl, "CGL", "Cp", ret.diff = TRUE)
+ DCp.cgl <- check.EOS(eos.cgl, "CGL", "Cp")
cat(paste("check.OBIGT: GHS for", sum(!isHKF), "cr,gas,liq species in", what, "\n"))
- DG.cgl <- checkGHS(eos.cgl, ret.diff = TRUE)
+ DG.cgl <- check.GHS(eos.cgl)
DCp[!isHKF] <- DCp.cgl
DG[!isHKF] <- DG.cgl
}
Modified: pkg/CHNOSZ/inst/NEWS.Rd
===================================================================
--- pkg/CHNOSZ/inst/NEWS.Rd 2023-06-21 06:24:29 UTC (rev 791)
+++ pkg/CHNOSZ/inst/NEWS.Rd 2023-06-23 14:00:30 UTC (rev 792)
@@ -12,15 +12,21 @@
% 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-10 (2023-06-20)}{
+\section{Changes in CHNOSZ version 2.0.0-12 (2023-06-23)}{
\itemize{
- \item Restore \file{EOSregress.R}, \file{eos-regress.Rmd}, and \file{demo/adenine.R}.
+ \item Add FAQ vignette.
+ \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.
+ \item Rename \code{checkGHS()} and \code{checkEOS()} to
+ \code{check.GHS()} and \code{check.EOS()}.
+
}
}
Modified: pkg/CHNOSZ/inst/tinytest/test-subcrt.R
===================================================================
--- pkg/CHNOSZ/inst/tinytest/test-subcrt.R 2023-06-21 06:24:29 UTC (rev 791)
+++ pkg/CHNOSZ/inst/tinytest/test-subcrt.R 2023-06-23 14:00:30 UTC (rev 792)
@@ -219,6 +219,10 @@
sres_nopoly_extrap <- subcrt("pyrrhotite", use.polymorphs = FALSE, exceed.Ttr = TRUE)
expect_false(anyNA(sres_nopoly_extrap$out[[1]]$G))
+# Added on 20230621
+info <- "Arguments 2 and 3 can't both be character"
+expect_error(subcrt(c("hydrogen", "H2"), c("gas", "aq"), "G"), info = info)
+
# References
# Amend, J. P. and Shock, E. L. (2001)
Modified: pkg/CHNOSZ/inst/tinytest/test-util.data.R
===================================================================
--- pkg/CHNOSZ/inst/tinytest/test-util.data.R 2023-06-21 06:24:29 UTC (rev 791)
+++ pkg/CHNOSZ/inst/tinytest/test-util.data.R 2023-06-23 14:00:30 UTC (rev 792)
@@ -7,19 +7,19 @@
d2 <- thermo()$OBIGT
expect_equal(d1, d2, info = info)
-info <- "checkGHS() and checkEOS() (via info()) produce messages"
+info <- "check.GHS() and check.EOS() (via info()) produce messages"
i1 <- info("S2O3-2")
-expect_message(info(i1), "G of S2O3-2\\(aq\\) differs by 939 cal mol-1 from tabulated value", info = info)
+expect_message(info(i1), "calculated G of S2O3-2\\(aq\\) differs by 939 cal mol-1 from database value", info = info)
i2 <- info("Cu+2")
-expect_message(info(i2), "Cp of Cu\\+2\\(aq\\) differs by 3.62 cal K-1 mol-1 from tabulated value", info = info)
+expect_message(info(i2), "calculated Cp of Cu\\+2\\(aq\\) differs by 3.62 cal K-1 mol-1 from database value", info = info)
-info <- "checkGHS() and checkEOS() respond to thermo()$opt$*.tol"
+info <- "check.GHS() and check.EOS() respond to thermo()$opt$*.tol"
i1 <- info("SO4-2")
thermo("opt$Cp.tol" = 0.5)
-expect_message(info(i1), "checkEOS", info = info)
+expect_message(info(i1), "check.EOS", info = info)
i2 <- info("a,w-dicarboxytetracosane")
thermo("opt$G.tol" = 50)
-expect_message(info(i2), "checkGHS", info = info)
+expect_message(info(i2), "check.GHS", info = info)
info <- "RH2OBIGT() gives group additivity results consistent with database values (from Richard and Helgeson, 1998)"
file <- system.file("extdata/adds/RH98_Table15.csv", package = "CHNOSZ")
@@ -64,7 +64,7 @@
expect_message(add.OBIGT(system.file("extdata/adds/LA19_test.csv", package = "CHNOSZ")), "energy units: J and cal", info = info)
expect_message(info(info("DMA_cal")), "-1.92 cal", info = info)
expect_message(info(info("DMA_J")), "-8.02 J", info = info)
-# For TMA, only a checkGHS message for the entry in J is produced,
+# For TMA, only a check.GHS message for the entry in J is produced,
# because it's above the threshold of 100 set in thermo()$opt$G.tol
expect_silent(info(info("TMA_cal")), info = info)
expect_message(info(info("TMA_J")), "-102 J", info = info)
@@ -120,7 +120,7 @@
# Test added 20230220
icr <- mod.OBIGT("fake_cr", formula = "Na2Cl2", state = "cr", model = "CGL", G = -1000, H = -1000, S = 10, Cp = 10, V = 10)
iaq <- mod.OBIGT("fake_aq", formula = "Na2Cl2", state = "aq", model = "CGL", G = -1000, H = -1000, S = 10, Cp = 10, V = 10)
-# Make sure info() runs (message is from checkEOS())
+# Make sure info() runs (message is from check.EOS())
expect_message(info(iaq), "differs", info = info)
expect_silent(info(iaq), info = info)
# Make sure subcrt() runs
Modified: pkg/CHNOSZ/man/extdata.Rd
===================================================================
--- pkg/CHNOSZ/man/extdata.Rd 2023-06-21 06:24:29 UTC (rev 791)
+++ pkg/CHNOSZ/man/extdata.Rd 2023-06-23 14:00:30 UTC (rev 792)
@@ -60,7 +60,7 @@
\item \code{OBIGT_check.csv} contains the results of running \code{\link{check.OBIGT}} to check the internal consistency of entries in the default and optional datafiles.
\item \code{RH98_Table15.csv} Group stoichiometries for high molecular weight crystalline and liquid organic compounds taken from Table 15 of Richard and Helgeson, 1998. The first three columns have the \code{compound} name, \code{formula} and physical \code{state} (\samp{cr} or \samp{liq}). The remaining columns have the numbers of each group in the compound; the names of the groups (columns) correspond to species in \code{\link{thermo}$OBIGT}. The compound named \samp{5a(H),14a(H)-cholestane} in the paper has been changed to \samp{5a(H),14b(H)-cholestane} here to match the group stoichiometry given in the table. See \code{\link{RH2OBIGT}} for a function that uses this file.
\item \code{SK95.csv} contains thermodynamic data for alanate, glycinate, and their complexes with metals, taken from Amend and Helgeson (1997) and Shock and Koretsky (1995) as corrected in slop98.dat. These data are used in \code{demo("copper")} and \code{demo("glycinate")}.
- \item \code{LA19_test.csv} contains thermodynamic data for dimethylamine and trimethylamine from LaRowe and Amend (2019) in energy units of both J and cal. This file is used in \code{test-util.data.R}) to check the messages produced by \code{\link{checkGHS}} and \code{\link{checkEOS}}.
+ \item \code{LA19_test.csv} contains thermodynamic data for dimethylamine and trimethylamine from LaRowe and Amend (2019) in energy units of both J and cal. This file is used in \code{test-util.data.R}) to check the messages produced by \code{\link{check.GHS}} and \code{\link{check.EOS}}.
}
Modified: pkg/CHNOSZ/man/info.Rd
===================================================================
--- pkg/CHNOSZ/man/info.Rd 2023-06-21 06:24:29 UTC (rev 791)
+++ pkg/CHNOSZ/man/info.Rd 2023-06-23 14:00:30 UTC (rev 792)
@@ -38,7 +38,7 @@
With a numeric argument, the rows of \code{thermo()$OBIGT} indicated by \code{ispecies} are returned, after removing any order-of-magnitude scaling factors (see \code{\link{thermo}}).
If these species are all aqueous or are all not aqueous, the compounded column names used in \code{thermo()$OBIGT} are replaced with names appropriate for the corresponding equations of state.
A missing value of one of the standard molal Gibbs energy (\code{G}) or enthalpy (\code{H}) of formation from the elements or entropy (\code{S}) is calculated from the other two, if available.
-If \code{check.it} is TRUE, several checks of self consistency among the thermodynamic properties and parameters are performed using \code{\link{checkGHS}} and \code{\link{checkEOS}}.
+If \code{check.it} is TRUE, several checks of self consistency among the thermodynamic properties and parameters are performed using \code{\link{check.GHS}} and \code{\link{check.EOS}}.
}
Modified: pkg/CHNOSZ/man/thermo.Rd
===================================================================
--- pkg/CHNOSZ/man/thermo.Rd 2023-06-21 06:24:29 UTC (rev 791)
+++ pkg/CHNOSZ/man/thermo.Rd 2023-06-23 14:00:30 UTC (rev 792)
@@ -50,9 +50,9 @@
\code{P.units} \tab character \tab The user's units of pressure ([\samp{bar}] or \samp{MPa})\cr
\code{state} \tab character \tab The default physical state for searching species [\samp{aq}] (see \code{\link{info}})\cr
\code{water} \tab character \tab Computational option for properties of water ([\samp{SUPCRT}] or \samp{IAPWS}; see \code{\link{water}})\cr
- \code{G.tol} \tab numeric \tab Difference in G above which \code{\link{checkGHS}} produces a message (cal mol\ifelse{latex}{\eqn{^{-1}}}{\ifelse{html}{\out{<sup>-1</sup>}}{^-1}}) [100]\cr
- \code{Cp.tol} \tab numeric \tab Difference in Cp above which \code{\link{checkEOS}} produces a message (cal K\ifelse{latex}{\eqn{^{-1}}}{\ifelse{html}{\out{<sup>-1</sup>}}{^-1}} mol\ifelse{latex}{\eqn{^{-1}}}{\ifelse{html}{\out{<sup>-1</sup>}}{^-1}}) [1]\cr
- \code{V.tol} \tab numeric \tab Difference in V above which \code{\link{checkEOS}} produces a message (cm\ifelse{latex}{\eqn{^{3}}}{\ifelse{html}{\out{<sup>3</sup>}}{^3}} mol\ifelse{latex}{\eqn{^{-1}}}{\ifelse{html}{\out{<sup>-1</sup>}}{^-1}}) [1]\cr
+ \code{G.tol} \tab numeric \tab Difference in G above which \code{\link{check.GHS}} produces a message (cal mol\ifelse{latex}{\eqn{^{-1}}}{\ifelse{html}{\out{<sup>-1</sup>}}{^-1}}) [100]\cr
+ \code{Cp.tol} \tab numeric \tab Difference in Cp above which \code{\link{check.EOS}} produces a message (cal K\ifelse{latex}{\eqn{^{-1}}}{\ifelse{html}{\out{<sup>-1</sup>}}{^-1}} mol\ifelse{latex}{\eqn{^{-1}}}{\ifelse{html}{\out{<sup>-1</sup>}}{^-1}}) [1]\cr
+ \code{V.tol} \tab numeric \tab Difference in V above which \code{\link{check.EOS}} produces a message (cm\ifelse{latex}{\eqn{^{3}}}{\ifelse{html}{\out{<sup>3</sup>}}{^3}} mol\ifelse{latex}{\eqn{^{-1}}}{\ifelse{html}{\out{<sup>-1</sup>}}{^-1}}) [1]\cr
\code{varP} \tab logical \tab Use variable-pressure standard state for gases? [\code{FALSE}] (see \code{\link{subcrt}})\cr
\code{IAPWS.sat} \tab character \tab State of water for saturation properties [\samp{liquid}] (see \code{\link{util.water}})\cr
\code{paramin} \tab integer \tab Minimum number of calculations to launch parallel processes [1000] (see \code{\link{palply}}) \cr
Modified: pkg/CHNOSZ/man/util.data.Rd
===================================================================
--- pkg/CHNOSZ/man/util.data.Rd 2023-06-21 06:24:29 UTC (rev 791)
+++ pkg/CHNOSZ/man/util.data.Rd 2023-06-23 14:00:30 UTC (rev 792)
@@ -2,8 +2,8 @@
\name{util.data}
\alias{util.data}
\alias{thermo.refs}
-\alias{checkEOS}
-\alias{checkGHS}
+\alias{check.EOS}
+\alias{check.GHS}
\alias{check.OBIGT}
\alias{dumpdata}
\alias{RH2OBIGT}
@@ -15,8 +15,8 @@
\usage{
thermo.refs(key = NULL, keep.duplicates = FALSE)
- checkEOS(eos, model, prop, ret.diff = FALSE)
- checkGHS(ghs, ret.diff = FALSE)
+ check.EOS(eos, model, prop, ret.diff = TRUE)
+ check.GHS(ghs, ret.diff = TRUE)
check.OBIGT()
dumpdata(file)
RH2OBIGT(compound = NULL, state = "cr",
@@ -29,7 +29,7 @@
\item{eos}{dataframe, equations-of-state parameters in the format of \code{thermo()$OBIGT}}
\item{model}{character, thermodynamic model (see \code{\link{thermo}})}
\item{prop}{character, property of interest (\samp{Cp} or \samp{V})}
- \item{ret.diff}{logical, return the difference between calculated and tabulated values?}
+ \item{ret.diff}{logical, return the difference between database and calculated values?}
\item{ghs}{dataframe, containing G, H and S, in the format of \code{thermo()$OBIGT}}
\item{file}{character, path to a file}
\item{compound}{character, name of compound(s) in group additivity calculation}
@@ -45,17 +45,20 @@
Only unique references are returned, unless \code{keep.duplicates} is TRUE.
In that case, a single reference for each species is returned, ignoring anything in \code{thermo()$OBIGT$ref2}.
-\code{checkEOS} compares heat capacity and volume calculated from equation-of-state parameters with reference (tabulated) values at 25 \degC and 1 bar and prints a message and returns the calculated value if tolerance is exceeded.
-The thermodynamic parameters should be provided in \code{eos}, which is a data frame with columns (and column names) in the same format as \code{\link{thermo}$OBIGT}.
-The \code{prop}erty can be one of \samp{Cp} or \samp{V}.
-The default tolerances, given in \code{thermo()$opt$Cp.tol} and \code{thermo()$opt$V.tol}, are 1 J/K.mol or 1 cal/K.mol for Cp and 1 cm3/mol for V.
-If \code{ret.diff} is TRUE, the differences are returned irrespective of their values, and no messages are printed.
+\code{check.EOS} calculates heat capacity (\code{prop = "Cp"}) or volume (\code{prop = "V"}) from equation-of-state parameters at 25 \degC and 1 bar.
+\code{check.GHS} calculates G (standard molal Gibbs energy of formation from the elements) from H (standard molal enthalpy of formation) and S (standard molal entropy) at 25 \degC and 1 bar.
+The calculated values of Cp, V, or G are then compared with the given values (i.e., database values).
+If \code{ret.diff} is TRUE (the default), the difference between the database and calculated values is returned.
+If \code{ret.diff} is FALSE, the difference is compared with a tolerance setting (see below).
+If the absolute value of the difference exceeds the tolerance, the function prints a message and returns the calculated value (not the difference) of the property.
+If the absolute value of the difference is less than the tolerance, the function returns NA with no message.
-\code{checkGHS} compares G (standard molal Gibbs energy of formation from the elements) calculated from H (standard molal enthalpy of formation) and S (standard molal entropy) with reference (tabulated) values of G at 25 \degC and 1 bar.
-A message is printed and the calculated difference is returned if it exceeds the value given in \code{thermo()$opt$G.tol}, which has a default value of 100 cal/mol.
-The calculation requires that G, H and S, and the chemical formula of the species all be present.
+For \code{check.EOS}, the thermodynamic parameters should be provided in \code{eos}, which is a data frame with column names in the same format as \code{\link{thermo}$OBIGT}.
+For \code{check.GHS}, the data frame should include G, H, S, and the chemical formula of the species.
+The default tolerances are 1 J/K.mol or 1 cal/K.mol for Cp (depending on the \code{E_units} for the species), 1 cm3/mol for V, and 100 cal/mol for G.
+These can be changed by setting \code{thermo()$opt$Cp.tol}, \code{thermo()$opt$V.tol}, and \code{thermo()$opt$G.tol}.
-\code{check.OBIGT} is a function to check self-consistency of each entry in the thermodynamic database, using \code{checkEOS} and \code{checkGHS}.
+\code{check.OBIGT} is a function to check self-consistency of each entry in the thermodynamic database, using \code{check.EOS} and \code{check.GHS}.
The output is a table listing only species that exceed at least one of the tolerance limits, giving the species index (rownumber in `thermo()$OBIGT`), species name and state, and DCp, DV and DG, for the calculated differences (only those above the tolerances are given).
Values of DCp and DG are given in the units present in the data files.
This function is used to generate the file found at \code{extdata/thermo/OBIGT_check.csv}.
Modified: pkg/CHNOSZ/vignettes/FAQ.Rmd
===================================================================
--- pkg/CHNOSZ/vignettes/FAQ.Rmd 2023-06-21 06:24:29 UTC (rev 791)
+++ pkg/CHNOSZ/vignettes/FAQ.Rmd 2023-06-23 14:00:30 UTC (rev 792)
@@ -13,6 +13,47 @@
csl: elementa.csl
---
+<style>
+/* CSS snippet from boostrap.css modified by jmd on 2023-06-23 */
+/*!
+ * Bootstrap v5.3.0 (https://getbootstrap.com/)
+ * Copyright 2011-2023 The Bootstrap Authors
+ * Licensed under MIT (https://github.com/twbs/bootstrap/blob/main/LICENSE)
+ */
+:root,
+[data-bs-theme=light] {
+ --bs-info-text-emphasis: #31708f;
+ --bs-info-bg-subtle: #d9edf7;
+ --bs-info-border-subtle: #9eeaf9;
+ --bs-border-width: 1px;
+ --bs-border-radius: 0.375rem;
+}
+.alert {
+ --bs-alert-bg: transparent;
+ --bs-alert-padding-x: 0.5rem;
+ --bs-alert-padding-y: 0.5rem;
+ --bs-alert-margin-bottom: 1rem;
+ --bs-alert-color: inherit;
+ --bs-alert-border-color: transparent;
+ --bs-alert-border: var(--bs-border-width) solid var(--bs-alert-border-color);
+ --bs-alert-border-radius: var(--bs-border-radius);
+ --bs-alert-link-color: inherit;
+ position: relative;
+ padding: var(--bs-alert-padding-y) var(--bs-alert-padding-x);
+ margin-bottom: var(--bs-alert-margin-bottom);
+ color: var(--bs-alert-color);
+ background-color: var(--bs-alert-bg);
+ border: var(--bs-alert-border);
+ border-radius: var(--bs-alert-border-radius);
+}
+.alert-info {
+ --bs-alert-color: var(--bs-info-text-emphasis);
+ --bs-alert-bg: var(--bs-info-bg-subtle);
+ --bs-alert-border-color: var(--bs-info-border-subtle);
+ --bs-alert-link-color: var(--bs-info-text-emphasis);
+}
+</style>
+
<script>
function ToggleDiv(ID) {
var D = document.getElementById("D-" + ID);
@@ -32,11 +73,30 @@
```{r setup, include = FALSE}
library(CHNOSZ)
options(width = 80)
-## use pngquant to optimize PNG images
+# Use pngquant to optimize PNG images
library(knitr)
knit_hooks$set(pngquant = hook_pngquant)
pngquant <- "--speed=1 --quality=0-25"
if (!nzchar(Sys.which("pngquant"))) pngquant <- NULL
+# To make warnings appear in text box 20230619
+# https://selbydavid.com/2017/06/18/rmarkdown-alerts/
+knitr::knit_hooks$set(
+ error = function(x, options) {
+ paste('\n\n<div class="alert alert-danger">',
+ gsub('##', '\n', gsub('^##\ Error', '**Error:**', x)),
+ '</div>', sep = '\n')
+ },
+ warning = function(x, options) {
+ paste('\n\n<div class="alert alert-warning">',
+ gsub('##', '\n', gsub('^##\ Warning:', '**Warning:**', x)),
+ '</div>', sep = '\n')
+ },
+ message = function(x, options) {
+ paste('\n\n<div class="alert alert-info">',
+ gsub('##', '\n', x),
+ '</div>', sep = '\n')
+ }
+)
```
```{r HTML, include = FALSE}
@@ -69,7 +129,7 @@
## How is 'CHNOSZ' pronounced?
-As one syllable, starting with an *sh* sound and [rhyming with *Oz*](https://en.wiktionary.org/wiki/Rhymes:English/%C9%92z).
+As one syllable that starts with an *sh* sound and [rhymes with *Oz*](https://en.wiktionary.org/wiki/Rhymes:English/%C9%92z).
CHNOSZ and [schnoz](https://en.wiktionary.org/wiki/schnozz) are homophones.
*Answered on 2023-05-22.*
@@ -185,4 +245,164 @@
*Answered on 2023-05-17.*
+## How can minerals with phase transitions be added to the database?
+
+Many minerals undergo phase transitions upon heating.
+The stable phases at different temperatures are called polymorphs.
+Each polymorph for a given mineral should have its own entry in OBIGT, including values of the standard thermodynamic properties (Δ*G*°, Δ*H*°, and *S*°) at 25 °C.
+The 25 °C (or 298.15 K) values for high-temperature polymorphs are often not listed in thermodynamic tables, but they can be calculated.
+This thermodynamic cycle shows how: we calculate the changes of a thermodyamic property (pictured here as `DS1` and `DS2`) between 298.15 K and the transition temperature (`Ttr`) for two polymorphs, then combine those with the property of the polymorphic transition (`DStr`) to obtain the difference of the property between the polymorphs at 298.15 K (`DS298`).
+
+```text
+ DStr DStr: entropy of transition between polymorphs 1 and 2
+ Ttr o---------->o Ttr: temperature of transition
+ ^ |
+ | |
+ DS1 | | DS2 DS1: entropy change of polymorph 1 from 298.15 K to Ttr
+ | | DS2: entropy change of polymorph 2 from 298.15 K to Ttr
+ | v
+ 298.15 K o==========>o DS298: entropy difference between polymorphs 1 and 2 at 298.15 K
+ DS298 DS298 = DS1 + DStr - DS2
+Polymorph 1 2
+```
+
+As an example, let's add pyrrhotite (Fe0.877S) from @PMW87.
+The formula and thermodynamic properties of this pyrrhotite differ from those of FeS from @HDNB78, which is already in OBIGT.
+We begin by defining all the input values in the next code block.
+In addition to `G`, `H`, `S`, and the heat capacity coefficients, non-NA values of volume (`V`) must be provided for the polymorph transitions to be calculated correctly by `subcrt()`.
+
+```{r pyrrhotite_values, message = FALSE}
+# The formula of the new mineral
+formula <- "Fe0.877S"
+# Because the properties from Pankratz et al. (1987) are listed in calories,
+# we need to change the output of subcrt() to calories (the default is Joules)
+E.units("cal")
+# Use temperature in Kelvin for the calculations below
+T.units("K")
+# Thermodynamic properties of polymorph 1 at 25 °C (298.15 K)
+G1 <- -25543
+H1 <- -25200
+S1 <- 14.531
+Cp1 <- 11.922
+# Heat capacity coefficients for polymorph 1
+a1 <- 7.510
+b1 <- 0.014798
+# For volume, use the value from Helgeson et al. (1978)
+V1 <- V2 <- 18.2
+# Transition temperature
+Ttr <- 598
+# Transition enthalpy (cal/mol)
+DHtr <- 95
+# Heat capacity coefficients for polymorph 2
+a2 <- -1.709
+b2 <- 0.011746
+c2 <- 3073400
+# Maximum temperature of polymorph 2
+T2 <- 1800
+```
+
+Use the temperature (`Ttr`) and enthalpy of transition (`DHtr`) to calculate the entropy of transition (`DStr`).
+Note that the Gibbs energy of transition (`DGtr`) is zero at `Ttr`.
+```{r pyrrhotite_Ttr, message = FALSE}
+DGtr <- 0 # DON'T CHANGE THIS
+TDStr <- DHtr - DGtr # TΔS° = ΔH° - ΔG°
+DStr <- TDStr / Ttr
+```
+
+Put the heat capacity coefficients and volume in the database.
+There is no tabulated value of `Cp` of the high-temperature polymorph extrapolated to 298.15 K, so leave it out.
+If the properties were in Joules, we would omit `E_units = "cal"` or change it to `E_units = "J"`.
+```{r pyrrhotite_Cp, results = "hide", message = FALSE}
+mod.OBIGT("pyrrhotite_new", formula = formula, state = "cr", E_units = "cal",
+ G = 0, H = 0, S = 0, V = V1, Cp = Cp1,
+ a = a1, b = b1, c = 0, d = 0, e = 0, f = 0, lambda = 0, T = Ttr)
+mod.OBIGT("pyrrhotite_new", formula = formula, state = "cr2", E_units = "cal",
+ G = 0, H = 0, S = 0, V = V2,
+ a = a2, b = b2, c = c2, d = 0, e = 0, f = 0, lambda = 0, T = T2)
+```
+
+For the time being, we set `G`, `H`, and `S` (i.e., the properties at 25 °C) to zero in order to easily calculate the temperature integrals of the properties from 298.15 K to `Ttr`.
+Values of zero are placeholders that don't satisfy Δ*G*° = Δ*H*° − *T*Δ*S*° for the formation properties from the elements for either polymorph, as indicated by the following messages.
+We will check again for self-consistency at the end of the example.
+```{r pyrrhotite_info, results = "hide", collapse = TRUE}
+info(info("pyrrhotite_new", "cr"))
+info(info("pyrrhotite_new", "cr2"))
+```
+
+In order to calculate the temperature integrals of Δ*G*° and Δ*H*°, we need not only the heat capacity coefficients but also actual values of *S*°.
+Therefore, we start by calculating the entropy changes of each polymorph from 298.15 to `r Ttr` K (`DS1` and `DS2`) and combining those with the entropy of transition to obtain the difference of entropy between the polymorphs at 298.15 K.
+For polymorph 1 (with `state = "cr"`) it's advisable to include `use.polymorphs = FALSE` to prevent `subcrt()` from trying to identify the most stable polymorph at the indicated temperature.
+```{r pyrrhotite_S, message = FALSE}
+DS1 <- subcrt("pyrrhotite_new", "cr", P = 1, T = Ttr, use.polymorphs = FALSE)$out[[1]]$S
+DS2 <- subcrt("pyrrhotite_new", "cr2", P = 1, T = Ttr)$out[[1]]$S
+DS25 <- DS1 + DStr - DS2
+```
+
+Put the values of *S*° at 298.15 in OBIGT, then calculate the changes of all thermodynamic properties of each polymorph between 298.15 K and `Ttr`.
+```{r pyrrhotite_D1_D2, message = FALSE, results = "hide"}
+mod.OBIGT("pyrrhotite_new", state = "cr", S = S1)
+mod.OBIGT("pyrrhotite_new", state = "cr2", S = S1 + DS25)
+D1 <- subcrt("pyrrhotite_new", "cr", P = 1, T = Ttr, use.polymorphs = FALSE)$out[[1]]
+D2 <- subcrt("pyrrhotite_new", "cr2", P = 1, T = Ttr)$out[[1]]
+```
+
+It's a good idea to check that the entropy of transition is calculated correctly.
+```{r pyrrhotite_check_S, results = "hide"}
+stopifnot(all.equal(D2$S - D1$S, DStr))
+```
+
+Now we're ready to add up the contributions to Δ*G*° and Δ*H*° from the left, top, and right sides of the cycle.
+This gives us the differences between the polymorphs at 298.15 K, which we use to make the final changes to the database.
+```{r pyrrhotite_DG25_DH25, results = "hide", message = FALSE}
+DG25 <- D1$G + DGtr - D2$G
+DH25 <- D1$H + DHtr - D2$H
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/chnosz -r 792
More information about the CHNOSZ-commits
mailing list