[CHNOSZ-commits] r759 - in pkg/CHNOSZ: . R inst inst/tinytest man vignettes
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Dec 8 10:00:48 CET 2022
Author: jedick
Date: 2022-12-08 10:00:48 +0100 (Thu, 08 Dec 2022)
New Revision: 759
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/man/logB.to.OBIGT.Rd
pkg/CHNOSZ/vignettes/customizing.Rmd
Log:
Change logB.to.OBIGT() to use HKF parameters c1, c2, and omega
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2022-12-01 06:51:13 UTC (rev 758)
+++ pkg/CHNOSZ/DESCRIPTION 2022-12-08 09:00:48 UTC (rev 759)
@@ -1,6 +1,6 @@
-Date: 2022-12-01
+Date: 2022-12-08
Package: CHNOSZ
-Version: 1.9.9-50
+Version: 1.9.9-51
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-12-01 06:51:13 UTC (rev 758)
+++ pkg/CHNOSZ/R/logB.to.OBIGT.R 2022-12-08 09:00:48 UTC (rev 759)
@@ -1,11 +1,15 @@
# CHNOSZ/logB.to.OBIGT.R
-# Fit G, S, and Cp to formation constants (logB) as a function of temperature
-# 20220324 jmd
+# Get thermodynamic parameters from formation constants (logB) as a function of temperature
+# 20220324 v1 Regress three parameters (G, S, and Cp)
+# 20221202 v2 Regress HKF parameters (assume constant pressure and optimize omega parameter for charged species)
-logB.to.OBIGT <- function(logB, species, coeff, T, P, tolerance = 0.05) {
+logB.to.OBIGT <- function(logB, species, coeffs, T, P, optimize.omega = TRUE, tolerance = 0.05) {
- # We need at least three temperature values
- if(length(unique(T)) < 3) stop("at least 3 unique temperature values are needed")
+ # We need at least five temperature values to optimize omega
+ if(optimize.omega & length(unique(T)) < 5) {
+ message("logB.to.OBIGT: forcing optimize.omega = FALSE because there are < 5 T values")
+ optimize.omega <- FALSE
+ }
## 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
@@ -14,7 +18,7 @@
# Be sure to zap properties of a formed species that is already in the database
suppressMessages(mod.OBIGT(name, formula = name, E_units = E.units(), G = 0, S = 0, Cp = 0, zap = TRUE))
# Calculate logK of the formation reaction with ΔG°f = 0 for the formed species
- logK0 <- suppressMessages(subcrt(species, coeff, T = T, P = P)$out$logK)
+ logK0 <- suppressMessages(subcrt(species, coeffs, T = T, P = P)$out$logK)
## Get Gibbs energy of species from logB of reaction
# Calculate T in Kelvin
@@ -26,30 +30,73 @@
# Calculate ΔG°f of the formed species
Gf <- Gr - Gr0
- ## Solve for G, S, and Cp
- # Make an 'lm' model object for given Cp
- Gfun <- function(Cp = 0) {
- Tr <- 298.15
- TTr <- TK - Tr
- # Subtract Cp term from Gf
- GCp <- Cp * (TK - Tr - TK * log(TK / Tr))
- GCp[is.na(GCp)] <- 0
- GfCp <- Gf - GCp
- # Write linear model in Ttr -- slope is -S
- lm(GfCp ~ TTr)
+ ## Fit to HKF 20221202
+
+ # Get blank set of HKF parameters
+ PAR <- info(info("H+"))
+ # DON'T keep the name H+, becuase it tells hkf() to not calculate g function
+ PAR$name <- name
+ # Insert charge of formed species (for calculating g function and variable omega)
+ Z <- makeup(c(name, "Z0"), sum = TRUE)["Z"]
+ # Force charge to be 0 when optimize.omega is FALSE 20221208
+ if(!optimize.omega) Z <- 0
+ PAR$Z <- Z
+ # Get values of T and P (so "Psat" becomes numeric) 20221208
+ tpargs <- TP.args(T = TK, P = P)
+ # Calculate variable terms for S, c1, c2, and omega
+ # NOTE: Units of Kelvin for HKF
+ S_var <- hkf("G", T = tpargs$T, P = tpargs$P, parameters = transform(PAR, S = 1))$aq[[1]]$G
+ c1_var <- hkf("G", T = tpargs$T, P = tpargs$P, parameters = transform(PAR, c1 = 1))$aq[[1]]$G
+ c2_var <- hkf("G", T = tpargs$T, P = tpargs$P, parameters = transform(PAR, c2 = 1))$aq[[1]]$G
+ omega_var <- hkf("G", T = tpargs$T, P = tpargs$P, parameters = transform(PAR, omega = 1))$aq[[1]]$G
+
+ if(!optimize.omega) {
+
+ # Perform linear regression (constant initial omega)
+ Glm <- lm(Gf ~ S_var + c1_var + c2_var + omega_var)
+ omega <- Glm$coefficients["omega_var"]
+
+ } else {
+
+ # Function to perform linear regression for a given omega(Pr,Tr)
+ Gfun <- function(omega) {
+ PAR$omega <- omega
+ omega_var <- hkf("G", T = tpargs$T, P = tpargs$P, parameters = PAR)$aq[[1]]$G
+ # Subtract omega term from Gf
+ Gf_omega <- Gf - omega_var
+ # Perform linear regression with all terms except omega
+ lm(Gf_omega ~ S_var + c1_var + c2_var)
+ }
+ # Calculate the sum of squares of residuals for given omega(Pr,Tr)
+ Sqfun <- function(omega) sum(Gfun(omega)$residuals ^ 2)
+ # Find the omega(Pr,Tr) that minimizes the sum of squares of residuals
+ omega <- optimize(Sqfun, c(-1e10, 1e10))$minimum
+ # Construct the linear model with this omega
+ Glm <- Gfun(omega)
+
}
- # Calculate the sum of squares of residuals for given Cp
- Sqfun <- function(Cp) sum(Gfun(Cp)$residuals ^ 2)
- # Find the Cp with minimum sum of squares of residuals
- Cp <- optimize(Sqfun, c(-100, 100))$minimum
- # Calculate the fitted G and S for this Cp
- G <- Gfun(Cp)$coefficients[[1]]
- S <- - Gfun(Cp)$coefficients[[2]]
+ # Get coefficients
+ G <- Glm$coefficients["(Intercept)"]
+ S <- Glm$coefficients["S_var"]
+ c1 <- Glm$coefficients["c1_var"]
+ c2 <- Glm$coefficients["c2_var"]
+ # NAs propagate as zero in the EOS
+ if(is.na(S)) S <- 0
+ 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())
+ PAR$c1 <- c1
+ PAR$c2 <- c2
+ PAR$omega <- omega
+ Cp <- hkf("Cp", parameters = PAR)$aq[[1]]$Cp
+
## Update the thermodynamic parameters of the formed species
- ispecies <- do.call(mod.OBIGT, list(name = name, G = G, S = S, Cp = Cp))
+ # NOTE: we have to apply OOM scaling
+ ispecies <- do.call(mod.OBIGT, list(name = name, G = G, S = S, Cp = Cp, c1 = c1, c2 = c2 * 1e-4, omega = omega * 1e-5, z = Z))
# Calculate logK of the formation reaction with "real" ΔG°f for the formed species
- logK <- suppressMessages(subcrt(species, coeff, T = T, P = P)$out$logK)
+ logK <- suppressMessages(subcrt(species, coeffs, T = T, P = P)$out$logK)
# Calculate the mean absolute difference
mad <- mean(abs(logK - logB))
message(paste("logB.to.OBIGT: mean absolute difference between logB (experimental) and logK (calculated) is", round(mad, 4)))
Modified: pkg/CHNOSZ/inst/NEWS.Rd
===================================================================
--- pkg/CHNOSZ/inst/NEWS.Rd 2022-12-01 06:51:13 UTC (rev 758)
+++ pkg/CHNOSZ/inst/NEWS.Rd 2022-12-08 09:00:48 UTC (rev 759)
@@ -12,7 +12,7 @@
% 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 1.9.9-50 (2022-12-01)}{
+\section{Changes in CHNOSZ version 1.9.9-51 (2022-12-08)}{
\subsection{MAJOR USER-VISIBLE CHANGES}{
\itemize{
@@ -73,9 +73,10 @@
\subsection{NEW FEATURES: THERMODYNAMIC DATA}{
\itemize{
- \item Add \code{logB.to.OBIGT()} to fit three thermodynamic parameters
- (\samp{G}, \samp{S}, and \samp{Cp} at 25 \degC) to formation constants of
- aqueous species as a function of temperature.
+ \item Add \code{logB.to.OBIGT()} to fit selected thermodynamic parameters
+ (\samp{G} and \samp{S} at 25 \degC and \samp{c1}, \samp{c2}, and
+ \samp{omega} HKF coefficients) to formation constants of aqueous species
+ as a function of temperature.
\item Add vignette \viglink{customizing} with description of database
format, data-entry conventions, and examples of customizing the
Modified: pkg/CHNOSZ/inst/tinytest/test-mod.OBIGT.R
===================================================================
--- pkg/CHNOSZ/inst/tinytest/test-mod.OBIGT.R 2022-12-01 06:51:13 UTC (rev 758)
+++ pkg/CHNOSZ/inst/tinytest/test-mod.OBIGT.R 2022-12-08 09:00:48 UTC (rev 759)
@@ -27,15 +27,6 @@
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))
-
# Test added 20220920
info <- "Can add > 1 species; different states get different default models; info() works after mod.OBIGT"
iC12C13 <- mod.OBIGT(c("X", "Y"), formula = c("C12", "C13"), state = c("aq", "cr"))
Modified: pkg/CHNOSZ/man/logB.to.OBIGT.Rd
===================================================================
--- pkg/CHNOSZ/man/logB.to.OBIGT.Rd 2022-12-01 06:51:13 UTC (rev 758)
+++ pkg/CHNOSZ/man/logB.to.OBIGT.Rd 2022-12-08 09:00:48 UTC (rev 759)
@@ -7,76 +7,94 @@
}
\usage{
- logB.to.OBIGT(logB, species, coeff, T, P, tolerance = 0.05)
+ logB.to.OBIGT(logB, species, coeffs, T, P, optimize.omega = TRUE, tolerance = 0.05)
}
\arguments{
\item{logB}{Values of \logB}
\item{species}{Species in the formation reaction (the formed species must be last)}
- \item{coeff}{Coefficients of the formation reaction}
+ \item{coeffs}{Coefficients of the formation reaction}
\item{T}{Temperature (units of \code{\link{T.units}})}
\item{P}{Temperature (\samp{Psat} or numerical values with units of \code{\link{P.units}})}
+ \item{optimize.omega}{logical, optimize the omega parameter (relevant for charged species)?}
\item{tolerance}{Tolerance for checking equivalence of input and calculated \logB values}
}
\details{
-\code{logB.to.OBIGT} starts with the decimal logarithms of equilibrium constants (\logB) for a formation reaction.
-The formation reaction is defined by \code{species} and \code{coeff}.
-All species in the formation reaction must be available in OBIGT except for the \emph{last} species, which must be the newly formed species.
-The properties of the known species are combined with \logB to calculate the standard Gibbs energy (G[T]) of the formed species as a function of \T and \P.
+This function can be used to update the \code{\link{OBIGT}} thermodynamic database with parameters fit to formation constants of aqueous species as a function of temperature.
+The \code{logB} argument should have decimal logarithm of equilibrium constants (\logB) for a formation reaction.
+The formation reaction is defined by \code{species} and \code{coeffs}.
+All species in the formation reaction must be available in OBIGT except for the \emph{last} species, which is the newly formed species.
-The values of G[T] are fit using the thermodynamic parameters \code{G}, \code{S}, and \code{Cp} (i.e., standard-state values at 25 °C and 1 bar).
-The equation is this: G[T] = G[Tr] + -S[Tr] * (T - Tr) + Cp[Tr] (T - Tr - T * ln(T / Tr)), where G, S, and Cp are standard-state Gibbs energy, entropy, and isobaric heat capacity, T is the temperature in Kelvin, Tr is the reference temperature (298.15 K), and ln stands for the natural logarithm.
-Note that the equation assumes a constant heat capacity and no volume change, and it may be unsuitable for extrapolation beyond the original \T and \P range.
+The function works as follows.
+First, the properties of the known species are combined with \logB to calculate the standard Gibbs energy (G[T]) of the formed species at each value of \code{T} and \code{P}.
+Then, the values of G[T] are used to retrieve the thermodynamic parameters \code{G}, \code{S}, \code{c1}, \code{c2}, and \code{omega}.
+The first two of these are standard-state values at 25 \degC and 1 bar, and the last three are parameters in the revised HKF equations (e.g. Eq. B25 of Shock et al., 1992).
+Volumetric parameters (\a1 to \a4) in the HKF equations currently aren't included, so the resulting fits are valid only at the pressure corresponding to the given values of \logB.
+The fitted parameters for the formed species are then added to OBIGT.
+\code{\link{all.equal}} is used to test for approximate equivalence of the input and calculated equilibrium constants; if the mean absolute difference exceeds \code{tolerance}, an error occurs.
-The fitted \code{G}, \code{S}, and \code{Cp} for the formed species are added to OBIGT.
-\code{\link{all.equal}} is used to test the near equivalence of the input and calculated equilibrium constants; if the mean absolute difference exceeds \code{tolerance}, an error occurs.
+The number of parameters used in the fit is not more than the number of data values (\emph{n}), and the parameters are added in the order stated above.
+If \emph{n} is less than 5, then the values of the unused parameters are set to 0 for addition to OBIGT.
+For instance, a single value of \logB would be fit only with \code{G}, implying that computed values of G[T] have no temperature dependence.
+For uncharged species, the value of \omega is a constant, but for charged species, it is a function of \T and \P as described by the \dQuote{g function} (Shock et al., 1992).
+By default an optimization step is used to refine the value of \code{omega} at 25 \degC and 1 bar for charged species.
+In case the data are not sufficient to constrain \code{omega} (see first example), this step can be skipped by setting \code{optimize.omega} to FALSE; then, it is assumed that \omega is constant during the fitting.
+In OBIGT, a constant omega for a charged species can be enforced by setting the \code{z} parameter to 0, but that is not done by this function.
+
}
\seealso{
\code{logB.to.OBIGT} calls \code{\link{mod.OBIGT}} with \code{zap = TRUE} to clear the parameters of a formed species if it already exists in the OBIGT database.
-If preexisting parameters (e.g. HKF coefficients) weren't cleared, they would interfere with the fitting done here, which uses only three parameters (\code{G}, \code{S}, and \code{Cp}).
+If preexisting parameters (e.g. volumetric HKF coefficients) weren't cleared, they would interfere with the fitting done here, which uses only selected parameters.
}
+\value{
+ The species index of the new species in OBIGT.
+}
+
\examples{
## CoHS+ from Migdisov et al. (2011)
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"
-inew <- logB.to.OBIGT(logB, species, coeff, T, P)
-# Print the new database entry
-# (Note: H is calculated from G and S)
-info(inew)
-# Make a plot of experimental logB
-xlim <- c()
-plot(T, logB, "n", c(100, 320), c(5.8, 6.8), xlab = axis.label("T"), ylab = quote(log~beta))
-points(T, logB, pch = 19, cex = 2)
-# Plot fitted values
-Tfit <- seq(100, 320, 10)
-sres <- subcrt(species, coeff, T = Tfit)
-lines(sres$out$T, sres$out$logK, col = 4)
-title(describe.reaction(sres$reaction))
-legend <- c("Migdisov et al. (2011)", "logB.to.OBIGT()")
-legend("topleft", legend, pch = c(19, NA), lty = c(0, 1), col = c(1, 4), pt.cex = 2)
+species <- c("Co+2", "HS-", "CoHS+")
+coeffs <- c(-1, -1, 1)
+opar <- par(mfrow = c(2, 1))
+for(optimize.omega in c(TRUE, FALSE)) {
+ # Fit the parameters with or without variable omega
+ inew <- logB.to.OBIGT(logB, species, coeffs, T, P, optimize.omega = optimize.omega)
+ # Print the new database entry
+ info(inew)
+ # Plot experimental logB
+ plot(T, logB, "n", c(100, 320), c(5.8, 6.8), xlab = axis.label("T"), ylab = quote(log~beta))
+ points(T, logB, pch = 19, cex = 2)
+ # Plot calculated values
+ Tfit <- seq(100, 320, 10)
+ sres <- subcrt(species, coeffs, T = Tfit)
+ lines(sres$out$T, sres$out$logK, col = 4)
+ title(describe.reaction(sres$reaction))
+ legend <- c("Migdisov et al. (2011)", paste0("logB.to.OBIGT(optimize.omega =",optimize.omega,")"))
+ legend("top", legend, pch = c(19, NA), lty = c(0, 1), col = c(1, 4), pt.cex = 2, bg = "#FFFFFFB0")
+}
+par(opar)
## ZnCl+ from Mei et al. (2015)
-# These are values for 5000 bar calculated from the modified Ryzhenko-Bryzgalin model
+# These are values for 5000 bar calculated from the modified Ryzhenko-Bryzgalin (RB) model
logB <- c(-1.93, -1.16, -0.38, 0.45, 1.15, 1.76, 2.30, 2.80, 3.26, 3.70, 4.12, 4.53, 4.92)
species <- c("Zn+2", "Cl-", "ZnCl+")
-coeff <- c(-1, -1, 1)
+coeffs <- c(-1, -1, 1)
T <- c(25, 60, 100, 150, 200, 250, 300, 350, 400, 450, 500, 550, 600)
P <- 5000
# Note: ZnCl+ is in the default OBIGT database;
# logB.to.OBIGT() "zaps" the previous parameters before adding the fitted ones
-inew <- logB.to.OBIGT(logB, species, coeff, T, P)
-# Plot MRB and logB.to.OBIGT fitted values
+inew <- logB.to.OBIGT(logB, species, coeffs, T, P)
+# Plot RB and logB.to.OBIGT values
plot(T, logB, xlab = axis.label("T"), ylab = axis.label("logB"), pch = 19, cex = 2)
Tfit <- seq(25, 600, 25)
-sres <- subcrt(species, coeff, T = Tfit, P = P)
+sres <- subcrt(species, coeffs, T = Tfit, P = P)
lines(sres$out$T, sres$out$logK, col = 4)
title(describe.reaction(sres$reaction), line = 3)
title("5000 bar", font.main = 1, line = 1)
@@ -88,4 +106,6 @@
Migdisov, Art. A., Zezin, D. and Williams-Jones, A. E. (2011) An experimental study of Cobalt (II) complexation in Cl\S{-} and H\s{2}S-bearing hydrothermal solutions. \emph{Geochim. Cosmochim. Acta} \bold{75}, 4065--4079. \doi{10.1016/j.gca.2011.05.003}
Mei, Y., Sherman, D. M., Liu, W., Etschmann, B., Testemale, D. and Brugger, J. (2015) Zinc complexation in chloride-rich hydrothermal fluids (25--600 \degC): A thermodynamic model derived from \emph{ab initio} molecular dynamics. \emph{Geochim. Cosmochim. Acta} \bold{150}, 264--284. \doi{10.1016/j.gca.2014.09.023}
+
+Shock, E. L., Oelkers, E. H., Johnson, J. W., Sverjensky, D. A. and Helgeson, H. C. (1992) Calculation of the thermodynamic properties of aqueous species at high pressures and temperatures: Effective electrostatic radii, dissociation constants and standard partial molal properties to 1000 \degC and 5 kbar. \emph{J. Chem. Soc. Faraday Trans.} \bold{88}, 803--826. \doi{10.1039/FT9928800803}
}
Modified: pkg/CHNOSZ/vignettes/customizing.Rmd
===================================================================
--- pkg/CHNOSZ/vignettes/customizing.Rmd 2022-12-01 06:51:13 UTC (rev 758)
+++ pkg/CHNOSZ/vignettes/customizing.Rmd 2022-12-08 09:00:48 UTC (rev 759)
@@ -40,6 +40,9 @@
e <- '<tt style="color: blue;">e</tt>'
f <- '<tt style="color: blue;">f</tt>'
lambda <- '<tt style="color: blue;">lambda</tt>'
+c1 <- '<tt style="color: blue;">c1</tt>'
+c2 <- '<tt style="color: blue;">c2</tt>'
+omega <- '<tt style="color: blue;">omega</tt>'
G_ <- '<tt style="color: blue;">G</tt>'
H_ <- '<tt style="color: blue;">H</tt>'
S_ <- '<tt style="color: blue;">S</tt>'
@@ -409,7 +412,7 @@
<!-- ######## SECTION MARKER ######## -->
## Case study: Formation constants for aqueous tungsten species
-Here we use `r logB.to.OBIGT_` to fit formation constants to thermodynamic parameters using a constant heat capacity.
+Here we use `r logB.to.OBIGT_` to fit to thermodynamic parameters to experimental formation constants.
Some additional steps are shown to refine a thermodynamic model to generate a speciation diagram as a function of pH.
### Fitting formation constants
@@ -416,7 +419,7 @@
`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);
+- Experimental decimal logarithms of formation constants (`r logB`) as a function of temperature;
- The stoichiometry of the formation reaction in terms of known species (the new species must be last);
- The experimental temperature and pressure.
@@ -423,8 +426,8 @@
`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.
+- Fits `r DG_0` of the new species using 25 °C thermodynamic properties and selected HKF model parameters (i.e., `r G_`, `r S_`, `r c1`, `r c2`, and `r omega` parameters in OBIGT);
+- Adds the parameters to OBIGT for use by other functions in CHNOSZ.
First we set the pressure for all `r logB` data.
```{r Psat}
@@ -432,6 +435,7 @@
```
Add first species: `r HWO4_` [@WTW_19].
+The fit to 3 temperature values involves only 3 parameters (`r G_`, `r S_`, and `r c1`), not including `r omega`.
```{r HWO4_}
T <- c(250, 300, 350)
logB <- c(5.58, 6.51, 7.99)
@@ -441,6 +445,7 @@
```
Add second species: `r H3WO4F2_` [@WWH_21].
+We get a better fit by setting `optimize.omega` to FALSE (although this is a charged species, it is treated as having a constant $\omega$).
```{r H3WO4F2-}
T <- seq(100, 250, 25)
logB <- c(17.00, 17.11, 17.46, 17.75, 18.17, 18.71, 19.23)
@@ -447,19 +452,19 @@
# 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)
+logB.to.OBIGT(logB, species, coeff, T, P, optimize.omega = FALSE)
```
Add third species: `r H2WO4` [@WWH_21].
-We have to increase the tolerance because there is considerable scatter in the experimental values.
+We 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)
+logB.to.OBIGT(logB, species, coeff, T, P, tolerance = 0.2)
```
-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.
+After running, `r logB.to.OBIGT_` returns just the species indices; the low values for `r HWO4_` (`r info("HWO4-")`) and `r H2WO4` (`r info("H2WO4")`) indicate that the function replaced parameters for these species that were already present in OBIGT.
### Diagram 1: Constant molality of `r F_`
@@ -515,7 +520,7 @@
@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_`.
+`r NOTE`: The coefficients in the model used by `r logB.to.OBIGT_` include 25 °C values of `r G_` and `r S_`.
In a conservative interpretation, these should be regarded only as *fitting parameters*, and should not be used to compute thermodynamic properties close to 25 °C unless they were fit to experimental data in that temperature range.
## References
More information about the CHNOSZ-commits
mailing list