[CHNOSZ-commits] r712 - in pkg/CHNOSZ: . R inst inst/tinytest vignettes
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Mar 25 00:57:51 CET 2022
Author: jedick
Date: 2022-03-25 00:57:51 +0100 (Fri, 25 Mar 2022)
New Revision: 712
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/R/logB_to_OBIGT.R
pkg/CHNOSZ/inst/NEWS.Rd
pkg/CHNOSZ/inst/tinytest/test-mod.OBIGT.R
pkg/CHNOSZ/vignettes/customizing.Rmd
pkg/CHNOSZ/vignettes/vig.bib
Log:
Add tungsten case study to customizing.Rmd
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2022-03-24 07:28:21 UTC (rev 711)
+++ pkg/CHNOSZ/DESCRIPTION 2022-03-24 23:57:51 UTC (rev 712)
@@ -1,6 +1,6 @@
Date: 2022-03-24
Package: CHNOSZ
-Version: 1.4.3-4
+Version: 1.4.3-5
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/R/logB_to_OBIGT.R
===================================================================
--- pkg/CHNOSZ/R/logB_to_OBIGT.R 2022-03-24 07:28:21 UTC (rev 711)
+++ pkg/CHNOSZ/R/logB_to_OBIGT.R 2022-03-24 23:57:51 UTC (rev 712)
@@ -4,6 +4,9 @@
logB_to_OBIGT <- function(logB, species, coeff, T, P, tolerance = 0.05) {
+ # We need at least three temperature values
+ if(length(unique(T)) < 3) stop("at least 3 unique temperature values are needed")
+
## Get the formula of the formed species (used as the species name in OBIGT)
### NOTE: the formed species has to be *last* in this list
name <- tail(species, 1)
Modified: pkg/CHNOSZ/inst/NEWS.Rd
===================================================================
--- pkg/CHNOSZ/inst/NEWS.Rd 2022-03-24 07:28:21 UTC (rev 711)
+++ pkg/CHNOSZ/inst/NEWS.Rd 2022-03-24 23:57:51 UTC (rev 712)
@@ -10,6 +10,26 @@
\newcommand{\s}{\ifelse{latex}{\eqn{_{#1}}}{\ifelse{html}{\out{<sub>#1</sub>}}{#1}}}
\newcommand{\S}{\ifelse{latex}{\eqn{^{#1}}}{\ifelse{html}{\out{<sup>#1</sup>}}{^#1}}}
+\section{Changes in CHNOSZ version 1.4.3-4 (2022-03-24)}{
+
+ \subsection{NEW FEATURES}{
+ \itemize{
+
+ \item Add logB_to_OBIGT() to fit three thermodynamic parameters (G, S,
+ and Cp) to formation constants of aqueous species as a function of
+ temperature.
+
+ \item Add \code{zap} argument to mod.OBIGT() to clear parameters of
+ preexisting species.
+
+ \item Add vignette customizing.Rmd about customizing the thermodynamic
+ database.
+
+ }
+ }
+
+}
+
\section{Changes in CHNOSZ version 1.4.3 (2022-02-20)}{
\subsection{UNIT TESTING}{
Modified: pkg/CHNOSZ/inst/tinytest/test-mod.OBIGT.R
===================================================================
--- pkg/CHNOSZ/inst/tinytest/test-mod.OBIGT.R 2022-03-24 07:28:21 UTC (rev 711)
+++ pkg/CHNOSZ/inst/tinytest/test-mod.OBIGT.R 2022-03-24 23:57:51 UTC (rev 712)
@@ -2,27 +2,36 @@
reset()
info <- "minimal usage of mod.OBIGT() creates usable data entries"
-# we need at least a name and some property
+# We need at least a name and some property
expect_error(mod.OBIGT("test"), "species name and a property", info = info)
-# a valid formula is needed
+# A valid formula is needed
expect_warning(expect_error(mod.OBIGT("test", date = as.character(Sys.Date())), "is not a simple chemical formula", info = info),
"please supply a valid chemical formula", info = info)
-# the default state is aq
+# The default state is aq
expect_message(itest <- mod.OBIGT("test", formula = "Z0", date = as.character(Sys.Date())), "added test\\(aq\\)", info = info)
-# we should get NA values of G for a species with NA properties
+# We should get NA values of G for a species with NA properties
expect_true(all(is.na(subcrt(itest)$out[[1]]$G)), info = info)
-# a single value of G comes through to subcrt
+# A single value of G comes through to subcrt
mod.OBIGT("test", G = 100)
expect_equal(subcrt("test", T = 25, P = 1)$out[[1]]$G, 100, info = info)
-# values for Cp and c1 integrate to the same values of G
+# Values for Cp and c1 integrate to the same values of G
G.Cp <- subcrt(mod.OBIGT(list(name = "test", S = 0, Cp = 100)))$out[[1]]$G
G.c1 <- subcrt(mod.OBIGT(list(name = "test", S = 0, c1 = 100)))$out[[1]]$G
expect_equal(G.Cp, G.c1, info = info)
info <- "mod.OBIGT() works with numeric argument (species index)"
-# test added 20200707
+# Test added 20200707
i1 <- info("CH4")
-# previously this failed with an error
+# Previously this failed with an error
i2 <- mod.OBIGT(i1, name = "methane")
expect_identical(i1, i2, info = info)
expect_identical(info(i2)$name, "methane", info = info)
+
+# Test added 20220325
+info <- "logK_to_OBIGT() errors for < 3 temperature values"
+logB <- c(6.24, 6.02, 5.84, 5.97, 6.52)
+species <- c("Co+2", "HS-", "CoHS+")
+coeff <- c(-1, -1, 1)
+T <- c(120, 150, 200, 250, 300)
+P <- "Psat"
+expect_error(logB_to_OBIGT(logB[1:2], species, coeff, T[1:2], P))
Modified: pkg/CHNOSZ/vignettes/customizing.Rmd
===================================================================
--- pkg/CHNOSZ/vignettes/customizing.Rmd 2022-03-24 07:28:21 UTC (rev 711)
+++ pkg/CHNOSZ/vignettes/customizing.Rmd 2022-03-24 23:57:51 UTC (rev 712)
@@ -14,7 +14,12 @@
```{r setup, include = FALSE}
library(CHNOSZ)
-oldopt <- options(width = 72)
+options(width = 80)
+## 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
```
```{r HTML, include = FALSE}
@@ -34,7 +39,13 @@
e <- '<tt style="color: blue;">e</tt>'
f <- '<tt style="color: blue;">f</tt>'
lambda <- '<tt style="color: blue;">lambda</tt>'
-Cp <- '<tt style="color: blue;">Cp</tt>'
+G_ <- '<tt style="color: blue;">G</tt>'
+H_ <- '<tt style="color: blue;">H</tt>'
+S_ <- '<tt style="color: blue;">S</tt>'
+Cp_ <- '<tt style="color: blue;">Cp</tt>'
+V_ <- '<tt style="color: blue;">V</tt>'
+Cp_0 <- "<i>C<sub>p</sub></i>°"
+DG_0 <- "Δ<i>G</i>°"
# CHNOSZ functions
reset_ <- '<code style="color: red;">reset()</code>'
OBIGT_ <- '<code style="color: red;">OBIGT()</code>'
@@ -49,6 +60,14 @@
affinity_ <- '<code style="color: green;">affinity()</code>'
thermo.refs_ <- '<code style="color: green;">thermo.refs()</code>'
thermo_ <- '<code style="color: green;">thermo()</code>'
+# Math stuff
+logK <- "log <i>K</i>"
+logB <- "log β"
+F_ <- "F<sup>-</sup>"
+Hplus <- "H<sup>+</sup>"
+HWO4_ <- "HWO<sub>4</sub><sup>-</sup>"
+H2WO4 <- "H<sub>2</sub>WO<sub>4</sub>"
+H3WO4F2_ <- "H<sub>3</sub>WO<sub>4</sub>F<sub>2</sub><sup>-</sup>"
```
This vignette was compiled on `r Sys.Date()` with CHNOSZ version `r sessionInfo()$otherPkgs$CHNOSZ$Version`.
@@ -55,7 +74,7 @@
This vignette will cover some topics about using custom thermodynamic data in CHNOSZ.
The two main functions to remember are `r add.OBIGT_` to add data from a CSV file and `r mod.OBIGT_` to add or modify data through a function interface.
-A third function, `logB_to_OBIGT()`, has been provided to fit thermodynamic parameters to experimental formation constants (β).
+A third function, `r logB_to_OBIGT_`, is provided to fit thermodynamic parameters to experimental formation constants (`r logB`).
Before describing the methods to add or modify data, some notes on the basic structure of the database and data entry conventions are given.
Column names (or parts thereof) are formatted in blue (e.g. `r formula`), and important notes are highlighted in yellow (`r NOTE`).
@@ -110,7 +129,7 @@
`r thermo_` is a convenience function to access or modify parts of the `thermo` list object.
To see the first few entries in OBIGT, do this:
-```{r thermo_.OBIGT}
+```{r thermo.OBIGT}
head(thermo()$OBIGT)
```
@@ -154,7 +173,7 @@
OPTIONAL: Everything else.
Really, it depends on what you need.
-For instance, if you just want to use `r subcrt_` to calculate log*K* of a reaction from ΔG°f values at 25 °C, then G is the only parameter that is needed.
+For instance, if you just want to use `r subcrt_` to calculate `r logK` of a reaction from `r DG_0` of species at 25 °C, then `r G_` is the only parameter that is needed.
OPTIONAL but useful:
@@ -166,19 +185,19 @@
### NA or 0?
-If a character value or thermo_dynamic parameter (Columns 1-13) is unknown, use `NA`.
+If a character value or thermodynamic parameter (Columns 1-13) is unknown, use `NA`.
Note that a missing (blank) value in the file is treated as NA.
- Unknown values for character values (Columns 1-8) should be NA.
-- If you have only two of G, H, and S, then the missing one should be NA (definitely not 0). `r NOTE`: `r info_` and dependent functions (particularly `r subcrt_`) "know" the equation DG = DGH - TDS; together with the entropies of the elements, this equation is used to compute a single missing one of G, H, or S from the other two.
-- If you don't have Cp or V, then set it to NA.
- - If HKF or CGL parameters are present: they will be used to calculate Cp, so thermodynamic properties *can* be calculated at T > 25 °C.
+- If you have only two of G, H, and S, then the missing one should be NA (definitely not 0). `r NOTE`: `r info_` and dependent functions (particularly `r subcrt_`) "know" the equation DG = DGH - TDS; together with the entropies of the elements, this equation is used to compute a single missing one of `r G_`, `r H_`, or `r S_` from the other two.
+- If you don't have `r Cp_` or `r V_`, then set it to NA.
+ - If HKF or CGL parameters are present: they will be used to calculate `r Cp_`, so thermodynamic properties *can* be calculated at T > 25 °C.
- If HKF or CGL parameters aren't present: thermodynamic properties *can't* be calculated at T > 25 °C (NAs will propagate to higher T).
If an "equation-of-state" or heat capacity coefficient (Columns 14-20) is unknown, use 0.
-- Furthermore, if you would like to assume that Cp or V is 0, then set it to 0.
- - Thermodynamic properties will be extrapolated to T > 25 °C and P > 1 bar assuming that Cp and V are 0.
+- Furthermore, if you would like to assume that `r Cp_` or V is 0, then set it to 0.
+ - Thermodynamic properties will be extrapolated to T > 25 °C and P > 1 bar assuming that `r Cp_` and `r V_` are 0.
### OOM scaling and `r info_`
@@ -192,10 +211,10 @@
- Remove OOM scaling. This is used primarily by other functions in CHNOSZ to get a set of unscaled `r model` parameters for calculating thermodynamic properties as a function of T and P.
- Extract the HKF or CGL parts of column names (only if all matching species have the same `r model`).
-- Calculate a missing one of G, H, or S if two of them are present.
-- Cross-check G, H, and S if all of them are present, and print a message if the difference is above a threshold (see `?checkGHS`).
-- Calculate a missing Cp or V from the `r model` parameters, if possible.
-- Cross-check Cp or V (if present) against the `r model` parameters, if possible, and print a message if the difference is above a threshold (see `?checkEOS`).
+- Calculate a missing one of `r G_`, `r H_`, or `r S_` if two of them are present.
+- Cross-check `r G_`, `r H_`, and `r S_` if all of them are present, and print a message if the difference is above a threshold (see `?checkGHS`).
+- Calculate a missing `r Cp_` or `r V_` from the `r model` parameters, if possible.
+- Cross-check `r Cp_` or `r V_` (if present) against the `r model` parameters, if possible, and print a message if the difference is above a threshold (see `?checkEOS`).
`r NOTE`: `r info_` does **NOT** change the units; the values it displays (including possibly calculated ones) correspond to the `r E.units_` for that species. On the other hand, `r subcrt_` outputs values in the units selected with `r E.units_`.
@@ -216,7 +235,7 @@
thermo()$OBIGT[icr, ]
```
-Based on the values in the `r Cp` column, would you predict that thermodynamic properties at T > 25 °C could be calculated for all of these minerals?
+Based on the values in the `r Cp_` column, would you predict that thermodynamic properties at T > 25 °C could be calculated for all of these minerals?
Let's see ...
For conciseness we'll consider a relatively small temperature range and display only the `out` part of the `r subcrt_` output.
@@ -224,7 +243,7 @@
subcrt("orpiment,amorphous", T = c(25, 50, 75))$out[[1]]
```
-That makes sense; integrating Cp to calculate Gibbs energy and other thermo_dynamic properties would propagate NA, and that is what appears in the output.
+That makes sense; integrating NA `r Cp_` to calculate Gibbs energy and other thermodynamic properties would propagate NA, and that is what appears in the output.
Now let's run the calculation for the alpha phase of arsenic.
```{r arsenic}
@@ -232,7 +251,7 @@
```
What happened here?
-Even though there are no heat capacity coefficients (see above), there is a non-NA value of Cp, and that value is used together with the entropy for calculating Gibbs energy at T > 25 °C.
+Even though there are no heat capacity coefficients (see above), there is a non-NA value of `r Cp_`, and that value is used together with the entropy for calculating Gibbs energy at T > 25 °C.
Note that zero for the 25 °C values of G and H is not a mistake or a placeholder for unknown values (unknown values should be represented by NA).
This is the reference state for the element, for which G and H are by convention equal to zero.
@@ -243,9 +262,9 @@
```
Are you surprised?
-You might be if you only noticed the NA value for Cp in OBIGT.
-However, there are non-NA values for the heat capacity coefficients, which are used to calculate Cp as a function of temperature.
-When supplied with a numeric argument (species index), `r info_` actually does this to fill in missing 25 °C values of Cp, V, and G, H, or S if possible, in addition to removing OOM scaling and simplifying column names:
+You might be if you only noticed the NA value for `r Cp_` in OBIGT.
+However, there are non-NA values for the heat capacity coefficients, which are used to calculate `r Cp_0` as a function of temperature.
+When supplied with a numeric argument (species index), `r info_` actually does this to fill in missing 25 °C values of `r Cp_`, `r V_`, and `r G_`, `r H_`, or `r S_` if possible, in addition to removing OOM scaling and simplifying column names:
```{r info_.tin}
info(info("tin"))
@@ -404,12 +423,111 @@
<details>
<summary><span style="background-color: #DDDDDD">* Click to expand *</span></summary>
-TODO
+### Fitting formation constants
-- Refinement 1: Account for variable activity of F- as a function of pH
-- Refinement 2: Set ionic strength
-- Refinement 3: Use extrapolation based on Bryzgalin-Ryzhenko model
+`r logB_to_OBIGT_` requires three things:
+- Experimental decimal logarithms of formation constants (`r logB`) as a function of temperature (at least three different temperatures);
+- The stoichiometry of the formation reaction in terms of known species (the new species must be last);
+- The experimental temperature and pressure.
+
+`r logB_to_OBIGT_` does three things:
+
+- Combines the formation constants with standard Gibbs energies (`r DG_0`) of the known species to calculate `r DG_0` of the new species;
+- Fits `r DG_0` of the new species using a three-term equation with constant heat capacity (i.e., 25 °C `r G_`, `r S_`, and `r Cp_` parameters in OBIGT);
+- Adds the fitted parameters to OBIGT for use by other functions in CHNOSZ.
+
+First we set the pressure for all `r logB` data.
+```{r Psat}
+P <- "Psat"
+```
+
+Add first species: `r HWO4_` [@WTW_19].
+```{r HWO4_}
+T <- c(250, 300, 350)
+logB <- c(5.58, 6.51, 7.99)
+species <- c("WO4-2", "H+", "HWO4-")
+coeff <- c(-1, -1, 1)
+logB_to_OBIGT(logB, species, coeff, T, P)
+```
+
+Add second species: `r H3WO4F2_` [@WWH_21].
+```{r H3WO4F2-}
+T <- seq(100, 250, 25)
+logB <- c(17.00, 17.11, 17.46, 17.75, 18.17, 18.71, 19.23)
+# Species and coefficients in the formation reaction
+species <- c("H+", "WO4-2", "F-", "H3WO4F2-")
+coeff <- c(-3, -1, -2, 1)
+logB_to_OBIGT(logB, species, coeff, T, P)
+```
+
+Add third species: `r H2WO4` [@WWH_21].
+We have to increase the tolerance because there is considerable scatter in the experimental values.
+```{r H2WO4}
+logB <- c(7.12, 7.82, 7.07, 7.76, 7.59, 7.98, 8.28)
+species <- c("H+", "WO4-2", "H2WO4")
+coeff <- c(-2, -1, 1)
+logB_to_OBIGT(logB, species, coeff, T, P, tolerance = 0.3)
+```
+
+We can tell from the low species indices of `r HWO4_` (`r info("HWO4-")`) and `r H2WO4` (`r info("H2WO4")`) that `r logB_to_OBIGT_` replaced parameters for these species that were already present in OBIGT.
+
+### Diagram 1: Constant activity of `r F_`
+
+Now we're ready to make a speciation diagram.
+Our target is to reproduce Fig. 7b of @WWH_21, which is made for 300 °C.
+A constant activity of `r F_` is based on the assumption of complete dissociation of 0.1 m NaF (we'll change this later).
+An ionic strength of 1.2 mol/kg is estimated for a solution with 1.8 m NaCl (use `NaCl(300, m_tot = 1.8)`).
+
+```{r diagram1, message = FALSE, results = "hide", fig.width = 6, fig.height = 5, out.width = "75%", fig.align = "center", pngquant = pngquant}
+basis(c("H+", "WO4-2", "F-", "H2O", "O2"))
+basis("F-", log10(0.1))
+iaq <- retrieve("W", c("O", "H", "F"), "aq")
+species(iaq)
+a <- affinity(pH = c(2, 7), T = 300, IS = 1.2)
+e <- equilibrate(a)
+col <- c(1, 4, 5, 2)
+diagram(e, alpha = TRUE, col = col, lty = 1, lwd = 2, ylab = "Fraction total W")
+```
+
+This isn't quite the diagram we were looking for.
+The published diagram shows a broad region of coexistence of `r H3WO4F2_` and `r HWO4_` at pH < 5 and increasing abundance of `r H2WO4` at lower pH.
+
+### Diagram 2: Variable activity of `r F_`
+
+In reality, the activity of `r F_` depends strongly on pH according to the reaction `r Hplus` + `r F_` = HF.
+With a little algebra, we can calculate the activity of `r F_` (`a_F` in the code below) from the equilbrium constant of this reaction for a given total F concentration (`F_tot`).
+`r NOTE`: It is important to call `r subcrt_` with a non-zero `IS` so that it returns effective equilibrium constants corrected for ionic strength (see what hapens to the diagram if you set `IS = 0`).
+
+```{r a_F}
+T <- 300
+pH <- seq(2, 7, 0.1)
+logK_HF <- subcrt(c("H+", "F-", "HF"), c(-1, -1, 1), T = T, IS = 1.2)$out$logK
+F_tot <- 0.1
+a_F <- F_tot / (1 + 10^(logK_HF - pH))
+```
+
+Now that we have the activity of `r F_` as a function of pH, we can provide it in the call to `r affinity_`.
+```{r diagram2, message = FALSE, results = "hide", results = "hide", fig.width = 6, fig.height = 5, out.width = "75%", fig.align = "center", pngquant = pngquant}
+basis(c("H+", "WO4-2", "F-", "H2O", "O2"))
+iaq <- retrieve("W", c("O", "H", "F"), "aq")
+species(iaq)
+a <- affinity(pH = pH, "F-" = log10(a_F), T = T, IS = 1.2)
+e <- equilibrate(a)
+diagram(e, alpha = TRUE, col = col, lty = 1, lwd = 2, ylab = "Fraction total W")
+```
+
+That's more like it.
+We have captured the basic geometry of Fig. 7b in @WWH_21.
+For instance, in accord with the published diagram, `r HWO4_` plateaus at around 40% of total W, and `r H2WO4` and `r H3WO4F2_` are nearly equally abundant at pH = 2.
+
+The highest experimental temperature for the formation constants of `r H2WO4` and `r H3WO4F2_` is 250 °C, but this diagram is drawn for 300 °C.
+ at WWH_21 used the modified Ryzhenko-Bryzgalin (MRB) model to extrapolate to 300 °C.
+In contrast, we used a different model but obtained quite similar results.
+
+`r NOTE`: The coefficients in the model used by `r logB_to_OBIGT_` are 25 °C values of `r G_`, `r S_`, and `r Cp_`.
+In a conservative interpretation, these should be regarded only as *fitting parameters*, and should not be used to report thermodynamic properties close to 25 °C unless they were fit to experimental data in that temperature range.
+
</details>
## References
Modified: pkg/CHNOSZ/vignettes/vig.bib
===================================================================
--- pkg/CHNOSZ/vignettes/vig.bib 2022-03-24 07:28:21 UTC (rev 711)
+++ pkg/CHNOSZ/vignettes/vig.bib 2022-03-24 23:57:51 UTC (rev 712)
@@ -178,14 +178,14 @@
}
@Article{LH06a,
- author = {LaRowe, Douglas E. and Helgeson, Harold C.},
- journal = {Geochimica et Cosmochimica Acta},
- title = {{B}iomolecules in hydrothermal systems: {C}alculation of the standard molal thermodynamic properties of nucleic-acid bases, nucleosides, and nucleotides at elevated temperatures and pressures},
- year = {2006},
- volume = {70},
- number = {18},
- pages = {4680--4724},
- doi = {10.1016/j.gca.2006.04.010},
+ author = {LaRowe, Douglas E. and Helgeson, Harold C.},
+ journal = {Geochimica et Cosmochimica Acta},
+ title = {{B}iomolecules in hydrothermal systems: {C}alculation of the standard molal thermodynamic properties of nucleic-acid bases, nucleosides, and nucleotides at elevated temperatures and pressures},
+ year = {2006},
+ number = {18},
+ pages = {4680--4724},
+ volume = {70},
+ doi = {10.1016/j.gca.2006.04.010},
}
@Article{LH06b,
@@ -723,3 +723,23 @@
volume = {317},
doi = {10.2475/07.2017.01},
}
+
+ at Article{WWH_21,
+ author = {Xin-Song Wang and A. E. Williams-Jones and Rui-Zhong Hu and Lin-Bo Shang and Xian-Wu Bi},
+ journal = {Geochimica et Cosmochimica Acta},
+ title = {The role of fluorine in granite-related hydrothermal tungsten ore genesis: Results of experiments and modeling},
+ year = {2021},
+ pages = {170--187},
+ volume = {292},
+ doi = {10.1016/j.gca.2020.09.032},
+}
+
+ at Article{WTW_19,
+ author = {Xin-Song Wang and Alexander Timofeev and A. E. Williams-Jones and Lin-Bo Shang and Xian-Wu Bi},
+ journal = {Geochimica et Cosmochimica Acta},
+ title = {An experimental study of the solubility and speciation of tungsten in {NaCl}-bearing aqueous solutions at 250, 300, and 350 °{C}},
+ year = {2019},
+ pages = {313--329},
+ volume = {265},
+ doi = {10.1016/j.gca.2019.09.013},
+}
More information about the CHNOSZ-commits
mailing list