[CHNOSZ-commits] r245 - in pkg/CHNOSZ: . R data inst man tests/testthat vignettes

noreply at r-forge.r-project.org noreply at r-forge.r-project.org
Sat Oct 7 08:08:38 CEST 2017 

Author: jedick
Date: 2017-10-07 08:08:37 +0200 (Sat, 07 Oct 2017)
New Revision: 245

Removed:
   pkg/CHNOSZ/data/CHNOSZ.R
Modified:
   pkg/CHNOSZ/DESCRIPTION
   pkg/CHNOSZ/R/basis.R
   pkg/CHNOSZ/R/berman.R
   pkg/CHNOSZ/inst/NEWS
   pkg/CHNOSZ/man/basis.Rd
   pkg/CHNOSZ/man/berman.Rd
   pkg/CHNOSZ/man/data.Rd
   pkg/CHNOSZ/tests/testthat/test-berman.R
   pkg/CHNOSZ/vignettes/obigt.Rmd
Log:
berman(): improve disorder calculations


Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION	2017-10-07 03:48:11 UTC (rev 244)
+++ pkg/CHNOSZ/DESCRIPTION	2017-10-07 06:08:37 UTC (rev 245)
@@ -1,6 +1,6 @@
 Date: 2017-10-07
 Package: CHNOSZ
-Version: 1.1.0-43
+Version: 1.1.0-44
 Title: Thermodynamic Calculations for Geobiochemistry
 Author: Jeffrey Dick
 Maintainer: Jeffrey Dick <j3ffdick at gmail.com>

Modified: pkg/CHNOSZ/R/basis.R
===================================================================
--- pkg/CHNOSZ/R/basis.R	2017-10-07 03:48:11 UTC (rev 244)
+++ pkg/CHNOSZ/R/basis.R	2017-10-07 06:08:37 UTC (rev 245)
@@ -3,11 +3,13 @@
 
 basis <- function(species=NULL, state=NULL, logact=NULL, delete=FALSE) {
   thermo <- get("thermo")
+  oldbasis <- thermo$basis
   ## delete the basis species if requested
-  oldbasis <- thermo$basis
-  if(delete) {
+  if(delete | identical(species, "")) {
     thermo$basis <- NULL
+    thermo$species <- NULL
     assign("thermo", thermo, "CHNOSZ")
+    return(invisible(oldbasis))
   }
   ## return the basis definition if requested
   if(is.null(species)) return(oldbasis)
@@ -182,7 +184,7 @@
   # just list the keywords if none is specified
   if(is.null(key)) return(basis.key)
   # delete any previous basis definition
-  basis(delete=TRUE)
+  basis("")
   # match the keyword to the available ones
   ibase <- match(key, basis.key)
   if(is.na(ibase)) stop(paste(key, "is not a keyword for preset basis species"))

Modified: pkg/CHNOSZ/R/berman.R
===================================================================
--- pkg/CHNOSZ/R/berman.R	2017-10-07 03:48:11 UTC (rev 244)
+++ pkg/CHNOSZ/R/berman.R	2017-10-07 06:08:37 UTC (rev 245)
@@ -4,7 +4,7 @@
 #      in the system Na2O-K2O-CaO-MgO-FeO-Fe2O3-Al2O3-SiO2-TiO2-H2O-CO2.
 #      J. Petrol. 29, 445-522. https://doi.org/10.1093/petrology/29.2.445
 
-berman <- function(name, T = 298.15, P = 1, thisinfo=NULL, check.G=FALSE, calc.transition=TRUE, calc.disorder=FALSE, units="cal") {
+berman <- function(name, T = 298.15, P = 1, thisinfo=NULL, check.G=FALSE, calc.transition=TRUE, calc.disorder=TRUE, units="cal") {
   # reference temperature and pressure
   Pr <- 1
   Tr <- 298.15
@@ -124,19 +124,28 @@
     Tds <- T[iTds]
     # the upper integration limit is Tmax
     Tds[Tds > Tmax] <- Tmax
-    # Ber88 Eqs. 15, 16, 17, 18, 19
-    Cpds[iTds] <- d0 + d1 * Tds^-0.5 + d2 * Tds^-2 + d3 * Tds + d4 * Tds^2
-    Hds[iTds] <- d0 * (Tds - Tmin) + 2 * d1 * (Tds^-0.5 - Tmin^-0.5) -
-      d2 * (Tds^-1 - Tmin^-1) + d3 * (Tds^2 - Tmin^2) / 2 + d4 * (Tds^3 - Tmin^3) / 3
-    Sds[iTds] <- d0 * (log(Tds) - log(Tmin)) - 2 * d1 * (Tds^-0.5 - Tmin^-0.5) -
-      d2 * (Tds^-2 - Tmin^-2) / 2 + d3 * (Tds - Tmin) + d4 * (Tds^2 - Tmin^2) / 2
-    # we can't do this if d5 == 0 (dolomite and gehlenite)
+    # Ber88 Eqs. 15, 16, 17
+    Cpds[iTds] <- d0 + d1*Tds^-0.5 + d2*Tds^-2 + d3*Tds + d4*Tds^2
+    Hds[iTds] <- d0*(Tds - Tmin) + d1*(Tds^0.5 - Tmin^0.5)/0.5 +
+      d2*(Tds^-1 - Tmin^-1)/-1 + d3*(Tds^2 - Tmin^2)/2 + d4*(Tds^3 - Tmin^3)/3
+    Sds[iTds] <- d0*(log(Tds) - log(Tmin)) + d1*(Tds^-0.5 - Tmin^-0.5)/-0.5 +
+      d2*(Tds^-2 - Tmin^-2)/-2 + d3*(Tds - Tmin) + d4*(Tds^2 - Tmin^2)/2
+    # Eq. 18; we can't do this if d5 == 0 (dolomite and gehlenite)
+    # "d5 is a constant computed in such as way as to scale the disordring enthalpy to the volume of disorder" (Berman, 1988)
     if(d5 != 0) Vds <- Hds / d5
-    Gds <- Hds - T * Sds + Vds * (P - Pr)
-    # Gds above Tmax (Eq. 20)
+    # Berman puts the Vds term directly into Eq. 19 (commented below), but that necessarily makes Gds != Hds - T * Sds
+    #Gds <- Hds - T * Sds + Vds * (P - Pr)
+    # instead, we include the Vds term with Hds
+    Hds <- Hds + Vds * (P - Pr)
+    # disordering properties above Tmax (Eq. 20)
     ihigh <- T > Tmax
-    # note that Gds[ihigh] and Sds[ihigh] on the rhs were both calculated at Tmax (above)
-    Gds[ihigh] <- Gds[ihigh] - (T[ihigh] - Tmax) * Sds[ihigh]
+    # again, Berman put the Sds term (for T > Tmax) into Eq. 20 for Gds (commented below), which would also make Gds != Hds - T * Sds
+    #Gds[ihigh] <- Gds[ihigh] - (T[ihigh] - Tmax) * Sds[ihigh]
+    # instead, we add the Sds[ihigh] term to Hds
+    Hds[ihigh] <- Hds[ihigh] - (T[ihigh] - Tmax) * Sds[ihigh]
+    # by writing Gds = Hds - T * Sds, the above two changes w.r.t. Berman's
+    # equations affect the computed values only for Hds, not Gds
+    Gds <- Hds - T * Sds
     # apply the disorder contributions
     Ga <- Ga + Gds
     Ha <- Ha + Hds
@@ -147,8 +156,6 @@
 
   ### (for testing) use G = H - TS to check that integrals for H and S are written correctly
   Ga_fromHminusTS <- Ha - T * S
-  # FIXME: this check fails if disorder properties are calculated:
-  # berman("K-feldspar", T=convert(600, "K"), P=10000, calc.disorder=TRUE)
   if(!isTRUE(all.equal(Ga_fromHminusTS, Ga))) stop(paste0(name, ": incorrect integrals detected using DG = DH - T*S"))
 
   ### thermodynamic and unit conventions used in SUPCRT ###

Deleted: pkg/CHNOSZ/data/CHNOSZ.R
===================================================================
--- pkg/CHNOSZ/data/CHNOSZ.R	2017-10-07 03:48:11 UTC (rev 244)
+++ pkg/CHNOSZ/data/CHNOSZ.R	2017-10-07 06:08:37 UTC (rev 245)
@@ -1,21 +0,0 @@
-# CHNOSZ/data/CHNOSZ.R
-# clear system settings (basis/species)
-
-# we only work if the CHNOSZ environment exists
-if(!"CHNOSZ" %in% search()) {
-  message("data(CHNOSZ): please run data(thermo) first")
-} else {
-
-  local({
-    # get thermo from CHNOSZ environment
-    thermo <- get("thermo", "CHNOSZ")
-    # set basis,species components
-    thermo$basis <- NULL
-    thermo$species <- NULL
-    # place it in CHNOSZ environment
-    assign("thermo", thermo, "CHNOSZ")
-  })
-
-  message("data(CHNOSZ): cleared basis and species settings")
-
-}

Modified: pkg/CHNOSZ/inst/NEWS
===================================================================
--- pkg/CHNOSZ/inst/NEWS	2017-10-07 03:48:11 UTC (rev 244)
+++ pkg/CHNOSZ/inst/NEWS	2017-10-07 06:08:37 UTC (rev 245)
@@ -1,4 +1,4 @@
-CHANGES IN CHNOSZ 1.1.0-42 (2017-10-06)
+CHANGES IN CHNOSZ 1.1.0-44 (2017-10-07)
 ---------------------------------------
 
 MAJOR CHANGES:
@@ -40,9 +40,8 @@
 
 DATABASE UPDATES:
 
-- Add data(OBIGT) and data(CHNOSZ) commands to reset only the
-  thermodynamic database (OBIGT) or system settings (CHNOSZ: basis and
-  species) to their defaults. data(thermo) is still used to reset
+- Add data(OBIGT) command to reset only the thermodynamic database
+  (OBIGT) to its default. data(thermo) is still used to reset
   *everything* (the database and all computational and system settings).
 
 - In add.obigt(), using 'file' without a file suffix now can be used to

Modified: pkg/CHNOSZ/man/basis.Rd
===================================================================
--- pkg/CHNOSZ/man/basis.Rd	2017-10-07 03:48:11 UTC (rev 244)
+++ pkg/CHNOSZ/man/basis.Rd	2017-10-07 06:08:37 UTC (rev 245)
@@ -39,8 +39,8 @@
 If this occurs, any existing species definition (created by the \code{species} function) is deleted.
 However, \code{\link{swap.basis}} can be used to change the species (the compositions and/or physical states thereof) in the basis set while maintaining the list of species of interest, with the added benefit of equivalence of the chemical potentials of the elements before and after the swap.
 
-Call \code{basis} with \code{delete} set to TRUE to clear the basis definition.
-Any current basis definition (before being deleted) is returned by this call or by calling \code{basis} with all default arguments.
+Call \code{basis} with \code{delete} set to TRUE or code{species} set to \code{""} to clear the basis definition.
+This also deletes the \code{\link{species}} list, if any.
 
 If the value of \code{basis} is one of the keywords in the following table, the corresponding set of basis species is loaded, and their activities are given preset values.
 This approach is used by many of the examples in the package.
@@ -63,7 +63,7 @@
 }
 
 \value{
-Returns the value of \code{thermo$basis} after any modifications; or, if \code{delete} is TRUE, its value before deletion.
+Returns the value of \code{thermo$basis} after any modifications; or, if \code{delete} is TRUE, its value before deletion (invisibly).
 }
 
 \section{Side Effects}{
@@ -91,7 +91,7 @@
 basis(c("H2O", "O2"))
 basis(c("H2O", "O2", "H+"))
 ## clear the basis species
-basis(delete=TRUE)
+basis("")
 
 \dontrun{
 ## marked dontrun because they produce errors

Modified: pkg/CHNOSZ/man/berman.Rd
===================================================================
--- pkg/CHNOSZ/man/berman.Rd	2017-10-07 03:48:11 UTC (rev 244)
+++ pkg/CHNOSZ/man/berman.Rd	2017-10-07 06:08:37 UTC (rev 245)
@@ -8,7 +8,7 @@
 
 \usage{
   berman(name, T = 298.15, P = 1, thisinfo = NULL, check.G = FALSE,
-         calc.transition = TRUE, calc.disorder = FALSE, units = "cal")
+         calc.transition = TRUE, calc.disorder = TRUE, units = "cal")
 }
 
 \arguments{
@@ -34,8 +34,6 @@
 A warning is produced if the absolute value of the difference between tabulated and calculated DGfTrPr is greater than 1000 J/mol.
 
 Providing \code{thisinfo} means that the mineral name is not searched in \code{thermo$obigt}, potentially saving some running time.
-
-FIXME: the default for \code{calc.disorder} is temporarily FALSE, as the code produces incorrect Gibbs energies.
 }
 
 \examples{

Modified: pkg/CHNOSZ/man/data.Rd
===================================================================
--- pkg/CHNOSZ/man/data.Rd	2017-10-07 03:48:11 UTC (rev 244)
+++ pkg/CHNOSZ/man/data.Rd	2017-10-07 06:08:37 UTC (rev 245)
@@ -3,7 +3,6 @@
 \docType{data}
 \alias{thermo}
 \alias{OBIGT}
-\alias{CHNOSZ}
 \alias{opt}
 \alias{element}
 \alias{refs}
@@ -15,7 +14,7 @@
 This includes the computational settings, thermodynamic database, and system settings (chemical species).
 
 The system settings are changed using \code{\link{basis}} and \code{\link{species}}.
-To restore the default system settings (no species loaded), use \code{data(CHNOSZ)}.
+To restore the default system settings (no species loaded), use \code{basis("")}.
 
 The thermodynamic database is changed using \code{\link{add.obigt}} and \code{\link{mod.obigt}}.
 To restore the default database, use \code{data(OBIGT)}.
@@ -24,7 +23,7 @@
 To restore the default computational settings, thermodynamic database, and system settings, use \code{data(thermo)}.
 
 In an interactive session, use \code{data(thermo)} to restore the thermodynamic database and computational and system settings to their initial state.
-Or, use \code{data(CHNOSZ)} to clear the system settings (basis and species), or \code{data(OBIGT)} to restore the default database.
+Or, use \code{basis("")} to clear the system settings (basis and species), or \code{data(OBIGT)} to restore the default database.
 
 The data files provided with CHNOSZ are in the \code{data} and \code{extdata/OBIGT} directories of the package.
 The \code{*.csv} files in these directories are used to build the \code{thermo} data object in an environment named \code{CHNOSZ}.
@@ -34,7 +33,6 @@
 \usage{
   data(thermo)
   data(OBIGT)
-  data(CHNOSZ)
 }
 
 \format{

Modified: pkg/CHNOSZ/tests/testthat/test-berman.R
===================================================================
--- pkg/CHNOSZ/tests/testthat/test-berman.R	2017-10-07 03:48:11 UTC (rev 244)
+++ pkg/CHNOSZ/tests/testthat/test-berman.R	2017-10-07 06:08:37 UTC (rev 245)
@@ -71,21 +71,19 @@
   # anadalusite: an uncomplicated mineral (no transitions)
   And_G <- c(-579368, -524987, -632421, -576834)
   And <- subcrt("andalusite", "cr_Berman", T=T, P=P)$out[[1]]
-  expect_maxdiff(And_G, And$G, 10)
+  expect_maxdiff(And$G, And_G, 7.5)
 
   # quartz: a mineral with polymorphic transitions
   aQz_G <- c(-202800, -179757, -223864, -200109)
   aQz <- subcrt("quartz", "cr_Berman", T=T, P=P)$out[[1]]
+  expect_maxdiff(aQz$G[-2], aQz_G[-2], 1.2)
   # here, the high-P, low-T point suffers
-  expect_maxdiff(aQz_G[-2], aQz$G[-2], 1.2)
-  expect_maxdiff(aQz_G[2], aQz$G[2], 1250)
+  expect_maxdiff(aQz$G[2], aQz_G[2], 1250)
 
   # K-feldspar: this one has disordering effects
   Kfs_G <- c(-888115, -776324, -988950, -874777)
   Kfs <- subcrt("K-feldspar", "cr_Berman", T=T, P=P)$out[[1]]
-  # we lose some accuracy at high T
-  expect_maxdiff(Kfs_G[1:2], Kfs$G[1:2], 7)
-  expect_maxdiff(Kfs_G[3:4], Kfs$G[3:4], 1600)
-
-  # (see test-subcrt.R for feldspar and quartz tests with Helgeson/SUPCRT92 values)
+  expect_maxdiff(Kfs$G[1:2], Kfs_G[1:2], 5)
+  # we are less consistent with the reference values at high T
+  expect_maxdiff(Kfs$G[3:4], Kfs_G[3:4], 350)
 })

Modified: pkg/CHNOSZ/vignettes/obigt.Rmd
===================================================================
--- pkg/CHNOSZ/vignettes/obigt.Rmd	2017-10-07 03:48:11 UTC (rev 244)
+++ pkg/CHNOSZ/vignettes/obigt.Rmd	2017-10-07 06:08:37 UTC (rev 245)
@@ -14,7 +14,7 @@
 bibliography: obigt.bib
 # so that these appear in the bibliography
 nocite: | 
-  @SPRONS92, @SLOP98, @SLOP07, @SLOP15, @CHNOSZ, @JOH92, @WP02, @CWM89, @PRPG97, @TH88, @Kul06
+  @SPRONS92, @SLOP98, @SLOP07, @SLOP15, @CHNOSZ, @JOH92, @WP02, @CWM89, @PRPG97, @TH88, @Kul06, @Sho09
 # 20170212: comment the csl line to build package on R-Forge, to avoid
 # getting an error there (pandoc-citeproc: error while parsing the XML string)
 csl: elementa.csl
@@ -123,9 +123,9 @@
 
 # Recent additions (late 2017)
 
-* Mineral data using the [Berman (1988)]() equations are listed under **Solids**.
+* Mineral data using the [Berman (1988)]() equations are listed under **Solids** / **Berman**.
 
-* Aqueous species data compatible with the [Deep Earth Water](http://www.dewcommunity.org/) model are listed under **Optional Data**.
+* Aqueous species data intended for use with the [Deep Earth Water](http://www.dewcommunity.org/) model are listed under **Optional Data** / **DEW**.
 
 # Sources of data {.tabset .tabset-fade}
 
@@ -184,7 +184,7 @@
 ### `r setfile("inorganic_cr.csv")`
 ```{r inorganic_cr, results="asis", echo=FALSE}
 cat("Chamosite,7A and witherite were present in sprons92.dat but not in slop98.dat or later files, and are not included in CHNOSZ.\n\n")
-cat("Parameters used here for goethite differ slightly from those listed in the slop files [@Sho09].<hr>")
+cat("The source of parameters used here for goethite is different from that in the slop files ([Shock, 2009](https://doi.org/10.2113/gsecongeo.104.8.1235)).<hr>")
 ```
 
 ```{r reflist, results="asis", echo=FALSE}
@@ -196,8 +196,8 @@
 
 ### `r setfile("Berman_cr.csv")`
 ```{r Berman_cr, results="asis", echo=FALSE}
-cat("This file gives the identifiying information for minerals whose properties are calculated using the formulation of @Ber88.\n")
-cat("To distinguish these minerals from the original set of mineral data in CHNOSZ [based on the compliation of @HDNB78], the physical states are listed as `cr_Berman`.\n")
+cat("This file gives the identifiying information for minerals whose properties are calculated using the formulation of [Berman (1988)](https://doi.org/10.1093/petrology/29.2.445).\n")
+cat("To distinguish these minerals from the original set of mineral data in CHNOSZ (based on the compliation of [Helgeson et al., 1978](http://www.worldcat.org/oclc/13594862)), the physical states are listed as `cr_Berman`.\n")
 cat("The actual data are stored separately (`extdata/Berman/*.csv`).<hr>")
 ```
 
@@ -223,7 +223,7 @@
 ## Optional Data {.tabset .tabset-pills}
 
 ```{r Optional_Data, results="asis", echo=FALSE}
-cat("These files contain optional data updates that are not loaded by default (using `data(thermo)`). To load these data, run `add.obigt('CHNOSZ_aq')` or `add.obigt('DEW_aq')`.\n\n")
+cat("These files contain optional data updates that are not loaded by default (using `data(thermo)`). To load these data, run `add.obigt('CHNOSZ')` or `add.obigt('DEW')`.\n\n")
 ```
 
 ### `r setfile("DEW_aq.csv")`
@@ -236,8 +236,8 @@
 
 ### `r setfile("CHNOSZ_aq.csv")`
 ```{r CHNOSZ_aq, results="asis", echo=FALSE}
-cat("The primary aqueous silica species in CHNSOZ is SiO<sub>2</sub> [@SHS89]. The pseudospecies H<sub>4</sub>SiO<sub>4</sub> is used to make activity diagrams with *a*<sub>H<sub>4</sub>SiO<sub>4</sub></sub> as a variable. The GHS and HKF parameters for this pseudospecies were calculated using CHNOSZ; see the vignette [*Regressing thermodynamic data*](eos-regress.html) for more information.\n\n")
-cat("HKF coefficients for adenine regressed from experimental Cp and V data were reported by Lowe et al. (2017). See [`demo(adenine)`](../demo/adenine.R).<hr>")
+cat("The primary aqueous silica species in CHNSOZ is SiO<sub>2</sub> ([Shock et al., 1989](https://doi.org/10.1016/0016-7037(89)90341-4)). The pseudospecies H<sub>4</sub>SiO<sub>4</sub> is used to make activity diagrams with *a*<sub>H<sub>4</sub>SiO<sub>4</sub></sub> as a variable. The GHS and HKF parameters for this pseudospecies were calculated using CHNOSZ; see the vignette [*Regressing thermodynamic data*](eos-regress.html) for more information.\n\n")
+cat("HKF coefficients for adenine regressed from experimental Cp and V data were reported by [Lowe et al. (2017)](https://doi.org/10.1016/j.jct.2017.04.005). See [`demo(adenine)`](../demo/adenine.R).<hr>")
 ```
 
 ```{r noreflist, results="asis", echo=FALSE}

More information about the CHNOSZ-commits mailing list