[CHNOSZ-commits] r161 - in pkg/CHNOSZ: . inst inst/extdata/thermo man vignettes

Sat Feb 18 11:15:38 CET 2017

Author: jedick
Date: 2017-02-18 11:15:38 +0100 (Sat, 18 Feb 2017)
New Revision: 161

Modified:
   pkg/CHNOSZ/DESCRIPTION
   pkg/CHNOSZ/inst/NEWS
   pkg/CHNOSZ/inst/extdata/thermo/OBIGT-2.csv
   pkg/CHNOSZ/man/data.Rd
   pkg/CHNOSZ/man/util.data.Rd
   pkg/CHNOSZ/vignettes/eos-regress.Rmd
Log:
eos-regress.Rmd: try to discuss H4SiO4_Ste01 carefully


Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================

--- pkg/CHNOSZ/DESCRIPTION	2017-02-18 07:17:18 UTC (rev 160)
+++ pkg/CHNOSZ/DESCRIPTION	2017-02-18 10:15:38 UTC (rev 161)
@@ -1,6 +1,6 @@
 Date: 2017-02-18
 Package: CHNOSZ
-Version: 1.0.8-50
+Version: 1.0.8-51
 Title: Chemical Thermodynamics and Activity Diagrams
 Author: Jeffrey Dick
 Maintainer: Jeffrey Dick <j3ffdick at gmail.com>

Modified: pkg/CHNOSZ/inst/NEWS
===================================================================
--- pkg/CHNOSZ/inst/NEWS	2017-02-18 07:17:18 UTC (rev 160)
+++ pkg/CHNOSZ/inst/NEWS	2017-02-18 10:15:38 UTC (rev 161)
@@ -1,4 +1,4 @@
-CHANGES IN CHNOSZ 1.0.8-50 (2017-02-18)
+CHANGES IN CHNOSZ 1.0.8-51 (2017-02-18)
 ---------------------------------------
 
 DOCUMENTATION:
@@ -80,7 +80,7 @@
 - OBIGT.csv: Add H4SiO4 pseudospecies (pseudo-H4SiO4); see calculations
   in eos-regress.Rmd.
 
-- OBIGT-2.csv: Rename H4SiO4 to pseudo-H4SiO4_Ste01 (Steffánson, 2001).
+- OBIGT-2.csv: Rename H4SiO4 to H4SiO4_Ste01 (Steffánson, 2001).
 
 - extdata: Add fasta/rubisco.fasta and cpetc/*.csv files SS98_Fig5a,
   SS98_Fig5b, PM90 (all used by anintro.Rmd).

Modified: pkg/CHNOSZ/inst/extdata/thermo/OBIGT-2.csv
===================================================================
--- pkg/CHNOSZ/inst/extdata/thermo/OBIGT-2.csv	2017-02-18 07:17:18 UTC (rev 160)
+++ pkg/CHNOSZ/inst/extdata/thermo/OBIGT-2.csv	2017-02-18 10:15:38 UTC (rev 161)
@@ -229,7 +229,7 @@
 YbCl2+,NA,YbCl2+,aq,MWW09,NA,03.Jul.10,-216360,NA,-19.22,-4.55,-0.6,1.9985,-2.9017,6.8906,-2.659,11.1191,-3.9614,0.8351,1
 LuCl2+,NA,LuCl2+,aq,MWW09,NA,03.Jul.10,-222750,NA,-30.31,11.26,-1.23,1.9757,-2.9572,6.9124,-2.6567,21.9175,-0.7409,1.0031,1
 Ti(OH)4,NA,Ti(OH)4,aq,Ste01,NA,31.Aug.06,-313690.248565966,-361331.022944551,5.23900573613767,48.4942638623327,40.6,7.36687380497132,10.210086042065,1.15296367112811,-3.20100382409178,51.2762906309751,-1.27652963671128,0.0149856596558317,0
-pseudo-H4SiO4_Ste01,NA,H4SiO4,aq,Ste01,NA,31.Aug.06,-312919.933078394,-348676.147227533,45.1003824091778,15.0095602294455,52.3,18.7299235181644,-21.2600382409178,18.6199808795411,-1.20004780114723,58.0305927342256,-20.789913957935,0.0869024856596558,0
+H4SiO4_Ste01,NA,H4SiO4,aq,Ste01,NA,31.Aug.06,-312919.933078394,-348676.147227533,45.1003824091778,15.0095602294455,52.3,18.7299235181644,-21.2600382409178,18.6199808795411,-1.20004780114723,58.0305927342256,-20.789913957935,0.0869024856596558,0
 arsenolite,NA,As2O3,cr,PGS+96,NA,25.Aug.06,-137640,-156930,25.67,23.15,51.12,20.9,17.73,-2.7,0,0,0,0,548
 claudetite,NA,As2O3,cr,PGS+96,NA,25.Aug.06,-137730,-156600,27.09,23.19,47.26,23.86,10.74,-3.44,0,0,0,0,588
 orpiment,NA,As2S3,cr,PGS+96,NA,25.Aug.06,-20550,-20700,39.1,27.53,70.5,22.73,13.82,0.605,0,0,0,0,580

Modified: pkg/CHNOSZ/man/data.Rd
===================================================================
--- pkg/CHNOSZ/man/data.Rd	2017-02-18 07:17:18 UTC (rev 160)
+++ pkg/CHNOSZ/man/data.Rd	2017-02-18 10:15:38 UTC (rev 161)
@@ -90,15 +90,9 @@
      \item The properties and parameters of some reference unfolded proteins used by Dick et al., 2006 are included here. Their names have dashes, instead of underscores, so that they are not confused with proteins whose properties are generated at runtime.
      \item The standard molal Gibbs energies and enthalpies of formation of the elements and entropies at 25 \eqn{^{\circ}}{°}C and 1 bar of aqueous metal-amino acid (alanate or glycinate) complexes reported by Shock and Koretsky, 1995 were recalculated by adding to their values the differences in the corresponding properties between the values for aqueous alanate and glycinate used by Shock and Koretsky, 1995, and those used by Amend and Helgeson, 1997b and Dick et al., 2006.
      \item The standard molal properties and equations-of-state parameters of four phase species (see below) of Fe(cr) were generated from heat capacity data given by Robie and Hemingway, 1995. 
+     \item \samp{pseudo-H4SiO4} is an approximate pseudospecies for H4SiO4. \samp{pseudo-H4SiO4} has thermodynamic properties that give zero properties for the reaction H4SiO4 = SiO2 + 2H2O at 25 \eqn{^{\circ}}{°}C, and HKF coefficients that are applicable to moderate temperatures (see the vignette \emph{Regressing thermodynamic data}).
   } 
 
-  Notes on the secondary database added by \code{add.obigt} (OBIGT-2.csv):
-
-  \itemize{
-     \item \samp{pseudo-H4SiO4} and \samp{pseudo-H4SiO4_Ste01} are approximate pseudospecies for H4SiO4. \samp{pseudo-H4SiO4} has thermodynamic properties that give zero properties for the reaction H4SiO4 = SiO2 + 2H2O at 25 \eqn{^{\circ}}{°}C, and HKF coefficients that are applicable to moderate temperatures (see the vignette `eos-regress.Rmd`). \samp{pseudo-H4SiO4_Ste01}, taken from Steffánson, 2001, is applicable near 300 \eqn{^{\circ}}{°}C (see Figure 3 of that paper). \samp{pseudo-H4SiO4} can be used to make diagrams with H4SiO4 as a basis species, instead of SiO2(aq).
-  } 
-
-
   These modifications are indicated in \code{OBIGT.csv} by having \samp{CHNOSZ} as one of the sources of data. Note also that some data appearing in the \acronym{slop98.dat} file (Shock et al., 1998) were corrected or modified as noted in that file, and are indicated in \code{OBIGT.csv} by having \samp{SLOP98} as one of the sources of data.
 
   In order to represent thermodynamic data for minerals with phase transitions, the different phases of these minerals are represented as phase species that have states denoted by \samp{cr1}, \samp{cr2}, etc. The standard molar thermodynamic properties at 25 \eqn{^{\circ}}{°}C and 1 bar (\eqn{T_r}{Pr} and \eqn{P_r}{Pr}) of the \samp{cr2} phase species of minerals were generated by first calculating those of the \samp{cr1} phase species at the transition temperature (\eqn{T_{tr}}{Ttr}) and 1 bar then taking account of the volume and entropy of transition (the latter can be retrieved by combining the former with the Clausius-Clapeyron equation and values of \eqn{(dP/dT)} of transitions taken from the \acronym{SUPCRT92} data file) to calculate the standard molar entropy of the \samp{cr2} phase species at \eqn{T_{tr}}{Ttr}, and taking account of the enthalpy of transition (\eqn{{\Delta}H^{\circ}}{DeltaH0}, taken from the \acronym{SUPCRT92} data file) to calculate the standard molar enthalpy of the \samp{cr2} phase species at \eqn{T_{tr}}{Ttr}. The standard molar properties of the \samp{cr2} phase species at \eqn{T_{tr}}{Ttr} and 1 bar calculated in this manner were combined with the equations-of-state parameters of the species to generate values of the standard molar properties at 25 \eqn{^{\circ}}{°}C and 1 bar. This process was repeated as necessary to generate the standard molar properties of phase species represented by \samp{cr3} and \samp{cr4}, referencing at each iteration the previously calculated values of the standard molar properties of the lower-temperature phase species (i.e., \samp{cr2} and \samp{cr3}). A consequence of tabulating the standard molar thermodynamic properties of the phase species is that the values of \eqn{(dP/dT)} and \eqn{{\Delta}H^{\circ}}{DeltaH0} of phase transitions can be calculated using the equations of state and therefore do not need to be stored in the thermodynamic database. However, the transition temperatures (\eqn{T_{tr}}{Ttr}) generally can not be assessed by comparing the Gibbs energies of phase species and are tabulated in the database.
@@ -253,8 +247,6 @@
 
   Shock, E. L. et al. {1998} \emph{SLOP98.dat} (computer data file). http://geopig.asu.edu/supcrt92_data/slop98.dat, accessed on 2005-11-05. Current location: \url{http://geopig.asu.edu/sites/default/files/slop98.dat}.
 
-  Stefánsson, Andri (2001) Dissolution of primary minerals of basalt in natural waters. I. Calculation of mineral solubilities from 0°C to 350°C \emph{Chem. Geol.} \bold{172}, 225--250. \url{http://dx.doi.org/10.1016/S0009-2541(00)00263-1}
-
   Wagman, D. D., Evans, W. H., Parker, V. B., Schumm, R. H., Halow, I., Bailey, S. M., Churney, K. L. and Nuttall, R. L. (1982) The NBS tables of chemical thermodynamic properties. Selected values for inorganic and C\eqn{_1}{1} and C\eqn{_2}{2} organic substances in SI units. \emph{J. Phys. Chem. Ref. Data} \bold{11} (supp. 2), 1--392. \url{https://srd.nist.gov/JPCRD/jpcrdS2Vol11.pdf}
 
 }

Modified: pkg/CHNOSZ/man/util.data.Rd
===================================================================
--- pkg/CHNOSZ/man/util.data.Rd	2017-02-18 07:17:18 UTC (rev 160)
+++ pkg/CHNOSZ/man/util.data.Rd	2017-02-18 10:15:38 UTC (rev 161)
@@ -49,6 +49,8 @@
 
   \code{add.obigt} reads data from the specified \code{file} and adds it to \code{\link{thermo}$obigt}. Set \code{force} to TRUE to replace species that exist in the thermodynamic database (each unique combination of a name and a state in the database is considered to be a unique species). \code{force}, if not specified, reverts to TRUE if the \code{file} is left at its default. Given the default setting of \code{E.units}, the function does not perform any unit conversions. If \code{E.units} is set to \samp{J}, then the thermodynamic parameters are converted from units of Joules to calories, as used in the CHNOSZ database. 
 
+  \code{add.obigt} affects the order of entries in \code{thermo}$obigt; therefore, it should be called before any \code{basis} or \code{species} definition. Also, there is no attempt made to keep the order of physical states in the database (aq-cr-liq-gas); the function simply adds new rows to the end of \code{thermo}$obigt. As a result, retrieving the properties of an added aqueous species using \code{\link{info}} requires an explicit \code{state="aq"} argument to that function if a species with the same name is in one of the cr, liq or gas states.
+
 \code{mod.obigt} changes one or more of the properties of species or adds species to the thermodynamic database.
 These changes are lost if you reload the database by calling \code{\link{data}(thermo)} or if you quit the \R session without saving it.
 The name of the species to add or change must be supplied as the first argument of \code{...} or as a named argument (named \samp{name}).
@@ -82,7 +84,11 @@
 }
 
 \section{Warning}{
-  \code{add.obigt} affects the order of entries in \code{thermo}$obigt; therefore, it should be called before any \code{basis} or \code{species} definition. Also, there is no attempt made to keep the order of physical states in the database (aq-cr-liq-gas); the function simply adds new rows to the end of \code{thermo}$obigt. As a result, retrieving the properties of an added aqueous species using \code{\link{info}} requires an explicit \code{state="aq"} argument to that function if a species with the same name is in one of the cr, liq or gas states.
+  Cautionary advice about the secondary database added by \code{add.obigt} (OBIGT-2.csv):
+
+  \itemize{
+     \item \samp{H4SiO4_Ste01}, taken from Steffánson, 2001, is used in the vignette \emph{Regressing thermodynamic data} for comparison with the H4SiO4 pseudospecies made in CHNOSZ (\samp{pseudo-H4SiO4}). As with other species in OBIGT-2.csv, it is likely to be inconsistent with specific minerals and other species in the default database, so use it with care.
+  } 
 }
 
 \value{
@@ -130,6 +136,9 @@
 
   Richard, L. and Helgeson, H. C. (1998) Calculation of the thermodynamic properties at elevated temperatures and pressures of saturated and aromatic high molecular weight solid and liquid hydrocarbons in kerogen, bitumen, petroleum, and other organic matter of biogeochemical interest. \emph{Geochim. Cosmochim. Acta} \bold{62}, 3591--3636. \url{http://dx.doi.org/10.1016/S0016-7037(97)00345-1}
 
+  Stefánsson, Andri (2001) Dissolution of primary minerals of basalt in natural waters. I. Calculation of mineral solubilities from 0°C to 350°C \emph{Chem. Geol.} \bold{172}, 225--250. \url{http://dx.doi.org/10.1016/S0009-2541(00)00263-1}
+
+
 }
 
 

Modified: pkg/CHNOSZ/vignettes/eos-regress.Rmd
===================================================================
--- pkg/CHNOSZ/vignettes/eos-regress.Rmd	2017-02-18 07:17:18 UTC (rev 160)
+++ pkg/CHNOSZ/vignettes/eos-regress.Rmd	2017-02-18 10:15:38 UTC (rev 161)
@@ -7,6 +7,7 @@
   tufte::tufte_html:
     tufte_features: ["background"]
     css: "vig.css"
+    toc: true
     mathjax: null
   tufte::tufte_handout:
     citation_package: natbib
@@ -95,7 +96,7 @@
 
 Applying a linear model (`lm`), the *coefficients* in the model can be obtained.
 In this case, they correspond to the HKF parameters, e.g. `c1` (the intercept), `c2` and ω.
-The coefficients are constants, meaning that **the model is applicable only to neutral (uncharged) species**, for which the high-temperature derivatives of ω are 0.
+The coefficients **in this model** are constants, restricing application of the model to neutral (uncharged) species, for which the high-temperature derivatives of ω are 0.
 Further below, a procedure is described to iteratively solve the equations for charged species.
 
 # An example for neutral species
@@ -147,7 +148,7 @@
 Note the second call to <span style="color:green">`EOSplot()`</span>, where the coefficients are not provided.
 This causes the data to be re-regressed, this time without the *T* < 600 K limitation.
 Although the overall fit, represented by the dotted line, is better, the fit to the low-temperature data is considerably worse.
-For nonelectroyltes, consider excluding extreme near-critical values from the regression (or modify the procedure to give them lower weight) in order to obtain more generally useful results.
+For nonelectrolytes, consider excluding extreme near-critical values from the regression (or modify the procedure to give them lower weight) in order to obtain more generally useful results.
 
 <span style="color:green">`EOScoeffs()`</span> is a small function to retrieve the HKF parameters from the database in CHNOSZ.
 The `stopifnot` expression tests that the regressed values are within 10% of the database values of `r c1`, `r c2` and ω for `r ch4`.
@@ -237,7 +238,7 @@
 # An example for charged species
 
 For this example, let's generate synthetic data for Na<sup>+</sup> using its parameters in the database.
-In the call to <span style="color:green">`subcrt()`</span> below, `convert=FALSE` means to take *T* in units of K.
+In the call to <span style="color:green">`subcrt()`</span> below, `convert = FALSE` means to take *T* in units of K.
 
 ```{r Nadat}
 T <- convert(seq(0, 600, 50), "K")
@@ -260,9 +261,9 @@
 The fits using constant ω are clearly not acceptable.
 
 There is, however, a way out.
-A different variable, `Cp_s_var`, can be used to specify the calculation of the "solvation" heat capacity in the HKF equations **using the temperature- and pressure-dependent corrections for charged species**. 
+A different variable, `Cp_s_var`, can be used to specify the calculation of the "solvation" heat capacity in the HKF equations using the temperature- and pressure-dependent corrections for charged species.
 To use this variable, the values of `r wPrTr` (omega at the reference temperature and pressure) and *Z* (charge) must be given, in addition to *T* and *P*.
-Of course, right now we **don't know** the value of `r wPrTr` -- it is the purpose of the regression to find it!
+Of course, right now we *don't know* the value of `r wPrTr`---it is the purpose of the regression to find it!
 But we can make a first guess using the value of ω found above.
 
 ```{r Navars1}
@@ -325,14 +326,14 @@
 
 Some mineral stability diagrams use the activity of `r h4sio4` as a variable.
 However, the primary species for dissolved silica in CHNOSZ is `r sio2`(aq) [@SHS89].
-The reaction with zero properties, `r h4sio4` = `r sio2`  + 2 `r h2o`, defines the properties of the pseudospecies `r h4sio4`.
+The pseudo-reaction with zero properties, `r h4sio4` = `r sio2`  + 2 `r h2o`, defines the properties of the pseudospecies `r h4sio4`.
 
 First we go about calculating the properties of `r sio2` + 2 `r h2o`.
 We do this over a range of *T* and *P*, but include many points near 25 °C to improve the fit of the regression in that region:
 `r op <- options(warn = -1)`
 ```{r SiO2_2H2O, message=FALSE}
 s_25C <- subcrt(c("SiO2", "H2O"), c(1, 2), T = 25)$out
-s_near25 <- subcrt(c("SiO2", "H2O"), c(1, 2), T = seq(20, 30, length.out = 50))$out
+s_near25 <- subcrt(c("SiO2", "H2O"), c(1, 2), T = seq(20, 30, length.out=50))$out
 s_lowT <- subcrt(c("SiO2", "H2O"), c(1, 2), T = seq(0, 100, 10))$out
 s_Psat <- subcrt(c("SiO2", "H2O"), c(1, 2))$out
 s_P500 <- subcrt(c("SiO2", "H2O"), c(1, 2), T = seq(0, 1000, 100), P = 500)$out
@@ -342,8 +343,9 @@
 
 Now we can start making the new species, with thermodynamic properties calculated at 25 °C:
 ```{r new_H4SiO4}
-mod.obigt("new-H4SiO4", formula = "H4SiO4", ref1 = "this_vignette", date = today(),
-          G = s_25C$G, H = s_25C$H, S = s_25C$S, Cp = s_25C$Cp, V = s_25C$V, z = 0)
+mod.obigt("new-H4SiO4", formula = "H4SiO4", ref1 = "this_vignette",
+          date = today(), G = s_25C$G, H = s_25C$H, S = s_25C$S,
+          Cp = s_25C$Cp, V = s_25C$V, z = 0)
 ```
 
 To prepare for the regression, combine the calculated data and convert °C to K:
@@ -385,36 +387,38 @@
 ```
 
 We just made a new `r h4sio4` pseudospecies.
-The database in CHNOSZ was previously updated with a `pseudo-H4SiO4` that was made using this procedure.
-The secondary database has a `pseudo-H4SiO4_Ste01` with properties of `r h4sio4` from @Ste01:
+In February 2017, the database in CHNOSZ was updated with a `pseudo-H4SiO4` that was made using this procedure.
+In addition, the secondary database has a `H4SiO4_Ste01` with properties of `r h4sio4` from @Ste01:
 ```{r width180, include=FALSE}
 options(width = 180)
 ```
 ```{r info_H4SiO4, message=FALSE}
 add.obigt()
-info(info(c("new-H4SiO4", "pseudo-H4SiO4", "pseudo-H4SiO4_Ste01")))
+info(info(c("new-H4SiO4", "pseudo-H4SiO4", "H4SiO4_Ste01")))
 ```
 ```{r width80, include=FALSE}
 ```
 
 ```{r subcrt_H4SiO4, fig.margin=TRUE, fig.width=4, fig.height=4, small.mar=TRUE, echo=FALSE, results="hide", message=FALSE, dpi=dpi, out.width="100%", cache=TRUE, fig.cap="Comparison of H<sub>4</sub>SiO<sub>4</sub> pseudospecies.", pngquant=pngquant}
 s1 <- subcrt(c("pseudo-H4SiO4", "SiO2", "H2O"), c(-1, 1, 2))
-plot(s1$out$T, s1$out$G, type = "l", ylim = c(-200, 500),
+plot(s1$out$T, s1$out$G, type = "l", ylim = c(-100, 600),
   xlab = axis.label("T"), ylab = axis.label("DG0"))
 add.obigt()
-s2 <- subcrt(c("pseudo-H4SiO4_Ste01", "SiO2", "H2O"), c(-1, 1, 2))
+s2 <- subcrt(c("H4SiO4_Ste01", "SiO2", "H2O"), c(-1, 1, 2))
 lines(s2$out$T, s2$out$G, lty = 2)
 abline(h = 0, lty = 3)
 legend("topright", legend = c("pseudo-H4SiO4 (CHNOSZ)",
-  "pseudo-H4SiO4_Ste01", "\t(Stefánsson, 2001)"), lty = c(1, 2, NA), bty = "n")
+  "H4SiO4_Ste01", "\t(Stefánsson, 2001)"), lty = c(1, 2, NA), bty = "n")
 text(225, 250, describe.reaction(s1$reaction))
 ```
 
-Let's compare `pseudo-H4SiO4` (with the same properties as `new-H4SiO4`) and `pseudo-H4SiO4_Ste01` with the actual properties of `r sio2` + 2 `r h2o` as a function of temperature:
+Let's compare CHNOSZ's `pseudo-H4SiO4` (with the same properties as the `new-H4SiO4` we just made) and `H4SiO4_Ste01` from Stefánsson (2001) with the calculated properties of `r sio2` + 2 `r h2o` as a function of temperature:
 ```{r subcrt_H4SiO4, eval=FALSE}
 ```
 
-Neither pseudospecies is perfect, but `pseudo-H4SiO4` is preferable for making low-temperature mineral activity diagrams---see <span style="color:blue">`?transfer`</span> and <span style="color:blue">`demo(activity_ratios)`</span>.
+The large differences for `H4SiO4_Ste01` at low temperature agree with the comparison to Shock et al. (1989) shown in Figure 3 of Stefánsson (2001).
+Therefore, `H4SiO4_Ste01` should not be used _in CHNOSZ_ for making low-temperature mineral activity diagrams.
+Instead, the pseudospecies with properties calculated here, `pseudo-H4SiO4`, is preferable for use _in CHNOSZ_---see <span style="color:blue">`?transfer`</span> and <span style="color:blue">`demo(activity_ratios)`</span>.
 
 # Other possibilities
 

