[CHNOSZ-commits] r870 - in pkg/CHNOSZ: . R inst inst/tinytest man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Jan 22 13:13:21 CET 2025
Author: jedick
Date: 2025-01-22 13:13:20 +0100 (Wed, 22 Jan 2025)
New Revision: 870
Added:
pkg/CHNOSZ/inst/tinytest/test-logK.to.OBIGT.R
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/R/logK.to.OBIGT.R
pkg/CHNOSZ/inst/NEWS.Rd
pkg/CHNOSZ/man/logK.to.OBIGT.Rd
Log:
Add 'state' argument to logK.to.OBIGT()
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2025-01-22 06:16:21 UTC (rev 869)
+++ pkg/CHNOSZ/DESCRIPTION 2025-01-22 12:13:20 UTC (rev 870)
@@ -1,6 +1,6 @@
Date: 2025-01-22
Package: CHNOSZ
-Version: 2.1.0-41
+Version: 2.1.0-42
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/logK.to.OBIGT.R
===================================================================
--- pkg/CHNOSZ/R/logK.to.OBIGT.R 2025-01-22 06:16:21 UTC (rev 869)
+++ pkg/CHNOSZ/R/logK.to.OBIGT.R 2025-01-22 12:13:20 UTC (rev 870)
@@ -3,8 +3,13 @@
# 20220324 v1 Regress three parameters (G, S, and Cp)
# 20221202 v2 Regress HKF parameters (assume constant pressure and optimize omega parameter for charged species)
-logK.to.OBIGT <- function(logK, species, coeffs, T, P, npar = 3, optimize.omega = FALSE, tolerance = 0.05) {
+## If this file is interactively sourced, the following are also needed to provide unexported functions:
+#source("util.args.R")
+#source("hkf.R")
+#source("cgl.R")
+logK.to.OBIGT <- function(logK, species, coeffs, T, P, name = NULL, state = "aq", V = NA, npar = 3, optimize.omega = FALSE, tolerance = 0.05) {
+
# We need at least five temperature values to optimize omega
if(optimize.omega & length(unique(T)) < 5) {
message("logK.to.OBIGT: forcing optimize.omega = FALSE because there are < 5 T values")
@@ -18,12 +23,17 @@
## 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)
+ formula <- tail(species, 1)
+ if(is.null(name)) name <- formula
# Add the formed species to the thermodynamic database with ΔG°f = 0 at all temperatures
- # 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))
+ # - but with a non-zero G if V != 0 (for minerals)
+ # - zap properties of a formed species that is already in the database
+ suppressMessages(inew <- mod.OBIGT(name, formula = formula, state = state, E_units = E.units(), G = 0, S = 0, Cp = 0, V = V, zap = TRUE))
+ # Use species indices in case formed species has same formula and state (but different name) as a species in the database
+ ispecies <- info(species)
+ ispecies[length(ispecies)] <- inew
# Calculate logK of the formation reaction with ΔG°f = 0 for the formed species
- logK0 <- suppressMessages(subcrt(species, coeffs, T = T, P = P)$out$logK)
+ logK0 <- suppressMessages(subcrt(ispecies, coeffs, T = T, P = P)$out$logK)
## Get Gibbs energy of species from logK of reaction
# Calculate T in Kelvin
@@ -35,51 +45,141 @@
# Calculate ΔG°f of the formed species
Gf <- Gr - Gr0
- ## Fit to HKF 20221202
+ if(state == "aq") {
- # 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
+ ## Fit to HKF 20221202
- if(!optimize.omega) {
+ # 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(formula, "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) {
+
+ # Default values
+ omega <- NA
+ c2 <- NA
+ c1 <- NA
+ S <- NA
+ if(npar == 5) {
+ # Perform linear regression with constant omega
+ Glm <- lm(Gf ~ S_var + c1_var + c2_var + omega_var)
+ omega <- Glm$coefficients["omega_var"]
+ c2 <- Glm$coefficients["c2_var"]
+ c1 <- Glm$coefficients["c1_var"]
+ S <- Glm$coefficients["S_var"]
+ G <- Glm$coefficients["(Intercept)"]
+ } else if(npar == 4) {
+ # Now with fewer parameters
+ Glm <- lm(Gf ~ S_var + c1_var + c2_var)
+ c2 <- Glm$coefficients["c2_var"]
+ c1 <- Glm$coefficients["c1_var"]
+ S <- Glm$coefficients["S_var"]
+ G <- Glm$coefficients["(Intercept)"]
+ } else if(npar == 3) {
+ Glm <- lm(Gf ~ S_var + c1_var)
+ c1 <- Glm$coefficients["c1_var"]
+ S <- Glm$coefficients["S_var"]
+ G <- Glm$coefficients["(Intercept)"]
+ } else if(npar == 2) {
+ Glm <- lm(Gf ~ S_var)
+ S <- Glm$coefficients["S_var"]
+ G <- Glm$coefficients["(Intercept)"]
+ } else if(npar == 1) {
+ G <- mean(Gf)
+ } else stop("invalid value for npar (should be from 1 to 5)")
+
+ } 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)
+ # Get coefficients other than omega
+ 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 HKF equations
+ 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 check.EOS())
+ 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
+ # NOTE: we have to apply OOM scaling for HKF
+ do.call(mod.OBIGT, list(name = name, formula = formula, state = state, G = G, S = S, Cp = Cp, c1 = c1, c2 = c2 * 1e-4, omega = omega * 1e-5, z = Z))
+
+ } else if(state == "cr") {
+
+ ## Fit to CGL (Maier-Kelley Cp equation) 20250122
+
+ # Get blank set of CGL parameters
+ PAR <- data.frame(name = NA, abbrv = NA, formula = NA, state = "cr", ref1 = NA, ref2 = NA, date = NA, model = "CGL", E_units = "J",
+ G = 0, H = 0, S = 0, Cp = 0, V = 0, a = 0, b = 0, c = 0, d = 0, e = 0, f = 0, lambda = 0, T = 0)
+ # Get values of T and P (so "Psat" becomes numeric) 20221208
+ tpargs <- TP.args(T = TK, P = P)
+ # Calculate variable terms for S, a, b, and c
+ S_var <- cgl("G", T = tpargs$T, P = tpargs$P, parameters = transform(PAR, S = 1))[[1]]$G
+ a_var <- cgl("G", T = tpargs$T, P = tpargs$P, parameters = transform(PAR, a = 1))[[1]]$G
+ b_var <- cgl("G", T = tpargs$T, P = tpargs$P, parameters = transform(PAR, b = 1))[[1]]$G
+ c_var <- cgl("G", T = tpargs$T, P = tpargs$P, parameters = transform(PAR, c = 1))[[1]]$G
+
# Default values
- omega <- NA
- c2 <- NA
- c1 <- NA
- S <- NA
+ c <- 0
+ b <- 0
+ a <- 0
+ S <- 0
if(npar == 5) {
- # Perform linear regression with constant omega
- Glm <- lm(Gf ~ S_var + c1_var + c2_var + omega_var)
- omega <- Glm$coefficients["omega_var"]
- c2 <- Glm$coefficients["c2_var"]
- c1 <- Glm$coefficients["c1_var"]
+ # Perform linear regression with all parameters
+ Glm <- lm(Gf ~ S_var + a_var + b_var + c_var)
+ c <- Glm$coefficients["c_var"]
+ b <- Glm$coefficients["b_var"]
+ a <- Glm$coefficients["a_var"]
S <- Glm$coefficients["S_var"]
G <- Glm$coefficients["(Intercept)"]
} else if(npar == 4) {
# Now with fewer parameters
- Glm <- lm(Gf ~ S_var + c1_var + c2_var)
- c2 <- Glm$coefficients["c2_var"]
- c1 <- Glm$coefficients["c1_var"]
+ Glm <- lm(Gf ~ S_var + a_var + b_var)
+ b <- Glm$coefficients["b_var"]
+ a <- Glm$coefficients["a_var"]
S <- Glm$coefficients["S_var"]
G <- Glm$coefficients["(Intercept)"]
} else if(npar == 3) {
- Glm <- lm(Gf ~ S_var + c1_var)
- c1 <- Glm$coefficients["c1_var"]
+ Glm <- lm(Gf ~ S_var + a_var)
+ a <- Glm$coefficients["a_var"]
S <- Glm$coefficients["S_var"]
G <- Glm$coefficients["(Intercept)"]
} else if(npar == 2) {
@@ -90,53 +190,32 @@
G <- mean(Gf)
} else stop("invalid value for npar (should be from 1 to 5)")
- } else {
+ # Calculate Cp at 25 °C (for info() and check.EOS())
+ PAR$a <- a
+ PAR$b <- b
+ PAR$c <- c
+ print(PAR)
+ Cp <- cgl("Cp", parameters = PAR)[[1]]$Cp
- # 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)
- # Get coefficients other than omega
- G <- Glm$coefficients["(Intercept)"]
- S <- Glm$coefficients["S_var"]
- c1 <- Glm$coefficients["c1_var"]
- c2 <- Glm$coefficients["c2_var"]
+ ## Update the thermodynamic parameters of the formed species
+ # NOTE: do NOT apply OOM scaling for CGL
+ # Use maximum provided T as T limit for Cp equation
+ do.call(mod.OBIGT, list(name = name, formula = formula, state = state, G = G, S = S, Cp = Cp, V = V,
+ a = a, b = b, c = c, d = 0, e = 0, f = 0, lambda = 0, T = max(TK)))
+ } else {
+ stop(paste0("Unrecognized state (", state, "); should be aq or cr"))
}
- # NAs propagate as zero in the HKF equations
- 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 check.EOS())
- 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
- # 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_calc <- suppressMessages(subcrt(species, coeffs, T = T, P = P)$out$logK)
+ logK_calc <- suppressMessages(subcrt(ispecies, coeffs, T = T, P = P)$out$logK)
+ print(inew)
# Calculate the mean absolute difference
mad <- mean(abs(logK_calc - logK))
message(paste("logK.to.OBIGT: mean absolute difference between experimental and calculated logK is", round(mad, 4)))
# Check that calculated values are close to input values
stopifnot(all.equal(logK_calc, logK, tolerance = tolerance, scale = 1))
- # Return the species index in OBIGT
- ispecies
+ # Return the index of the new species in OBIGT
+ inew
}
Modified: pkg/CHNOSZ/inst/NEWS.Rd
===================================================================
--- pkg/CHNOSZ/inst/NEWS.Rd 2025-01-22 06:16:21 UTC (rev 869)
+++ pkg/CHNOSZ/inst/NEWS.Rd 2025-01-22 12:13:20 UTC (rev 870)
@@ -15,7 +15,7 @@
\newcommand{\Cp}{\ifelse{latex}{\eqn{C_P}}{\ifelse{html}{\out{<I>C<sub>P</sub></I>}}{Cp}}}
\newcommand{\DG0}{\ifelse{latex}{\eqn{{\Delta}G^{\circ}}}{\ifelse{html}{\out{Δ<I>G</I>°}}{ΔG°}}}
-\section{Changes in CHNOSZ version 2.1.0-41 (2025-01-22)}{
+\section{Changes in CHNOSZ version 2.1.0-42 (2025-01-22)}{
\subsection{OBIGT DEFAULT DATA}{
\itemize{
@@ -81,71 +81,80 @@
}
}
- \subsection{OTHER CHANGES}{
+ \subsection{DOCUMENTATION}{
\itemize{
+ \item Add FAQ question: Why are mineral stability boundaries curved on
+ mosaic diagrams?
+
+ \item Add \file{demo/total_S.R}: total activity of S--pH diagram for
+ Fe-S-O-C minerals.
+
+ \item Add \file{demo/uranyl.R}: total (carbonate|sulfate)-pH diagrams
+ for uranyl species, after
+ \href{https://doi.org/10.1016/j.gca.2024.04.023}{Migdisov et al. (2024)}.
+
+ \item Revise \file{demo/sum_S.R} for log \emph{f}O\s{2} - log \emph{m} ΣS
+ diagram and add solubility contours for Fe and Au.
+
+ }
+ }
+
+ \subsection{NEW FEATURES}{
+ \itemize{
+
+ \item Rename \code{logB.to.OBIGT()} to \code{logK.to.OBIGT()} and add a
+ \strong{state} argument which accepts \samp{aq} (aqueous species using
+ heat capacity parameters in the \samp{HKF} model; the previous method and
+ current default) or \samp{cr} (minerals using the Maier-Kelley equation
+ in the \samp{CGL} model; a new method)).
+
\item \code{OBIGT} now uses \samp{CGL_Ttr} as the \code{model} for
species whose \code{T} value denotes a phase stability limit above which
- \code{subcrt()} returns NA values (in contrast to \samp{CGL}, for which
+ \code{subcrt()} returns NA values, in contrast to \samp{CGL}, for which
the value of \code{T} is a Cp equation limit that only produces a
- warning). This replaces the use of signed values of \code{T}.
+ warning. This replaces the use of signed values of \code{T}.
+ \item Adjust \code{axis.label()} to show sum symbol (Σ) for total
+ activity or molality of basis species used as a plotting variable in
+ \code{mosaic()}.
+
\item \code{check.EOS()} now uses values of Born coefficients \emph{X}
and \emph{Q} that are consistent with either SUPCRT92 or DEW, depending
on the \code{model} defined in OBIGT (HKF or DEW, respectively).
- \item Remove \code{ispecies} from output of \code{check.OBIGT()} to avoid
- superfluous diffs.
+ }
+ }
+
+ \subsection{OTHER CHANGES}{
+ \itemize{
+
\item Restore \code{lines} to the output of \code{diagram()} for the x
and y values of predominance field boundaries.
- \item Adjust \code{axis.label()} to show sum symbol (Σ) for total
- activity or molality of basis species used as a plotting variable in
- \code{mosaic()}.
-
- \item For better readability, formulas of species are no longer
- subscripted by \code{axis.label()} (via \code{expr.species()}). An
- example of the new formatting is log \emph{f}O\s{2}. Similar formatting
- for other expressions has been used in the vignettes.
-
- \item Merge \file{extdata/adds} and \file{extdata/cpetc} into
- \file{extdata/misc}.
-
\item In \code{diagram()}, the default for the \strong{type} argument
when using the output from \code{solubility()} is now \code{loga.balance}
(sum of activities of aqueous species) rather than \code{loga.equil}
(activities of individual aqueous species).
- \item In \code{NaCl()}, rename \strong{m_tot} to \strong{m_NaCl} (moles
- of of NaCl added to 1 kg H\s{2}O) and rename output with \samp{m_Naplus},
- \samp{m_Clminus}, \samp{m_NaCl0}, and \samp{m_HCl0} (molalities of
- aqueous species).
+ \item For better readability, formulas of species are no longer
+ subscripted by \code{axis.label()} (via \code{expr.species()}). An
+ example of the new formatting is log \emph{f}O\s{2}. Similar formatting
+ for other expressions has been used in the vignettes.
- \item Rename \code{logB.to.OBIGT()} to \code{logK.to.OBIGT()}.
+ \item In \code{NaCl()}, rename argument \strong{m_tot} to \strong{m_NaCl}
+ (moles of of NaCl added to 1 kg H\s{2}O) and rename output with
+ \samp{m_Naplus}, \samp{m_Clminus}, \samp{m_NaCl0}, and \samp{m_HCl0}
+ (molalities of aqueous species).
- }
- }
+ \item Remove \code{ispecies} from output of \code{check.OBIGT()} to avoid
+ superfluous diffs.
- \subsection{DOCUMENTATION}{
+ \item Merge \file{extdata/adds} and \file{extdata/cpetc} into
+ \file{extdata/misc}.
- \itemize{
-
- \item Add FAQ question: Why are mineral stability boundaries curved on
- mosaic diagrams?
-
- \item Add \file{demo/total_S.R}: total activity of S--pH diagram for
- Fe-S-O-C minerals.
-
- \item Add \file{demo/uranyl.R}: total (carbonate|sulfate)-pH diagrams
- for uranyl species, after
- \href{https://doi.org/10.1016/j.gca.2024.04.023}{Migdisov et al. (2024)}.
-
- \item Revise \file{demo/sum_S.R} for log \emph{f}O\s{2} - log \emph{m} ΣS
- diagram and add solubility contours for Fe and Au.
-
}
-
}
\subsection{REMOVED FEATURES}{
Added: pkg/CHNOSZ/inst/tinytest/test-logK.to.OBIGT.R
===================================================================
--- pkg/CHNOSZ/inst/tinytest/test-logK.to.OBIGT.R (rev 0)
+++ pkg/CHNOSZ/inst/tinytest/test-logK.to.OBIGT.R 2025-01-22 12:13:20 UTC (rev 870)
@@ -0,0 +1,19 @@
+# Test added 20250122
+info <- "Parameters are fit correctly for a mineral"
+species <- c("iron", "sulfur", "FeS2")
+coeffs <- c(-1, -2, 1)
+P <- 1000
+T <- seq(120, 400, 10)
+# Use calculated logK for formation of pyrite
+logK <- subcrt(c("iron", "sulfur", "pyrite"), c(-1, -2, 1), T = T, P = P)$out$logK
+# Use V for pyrite to reproduce value of G at 25 degC
+inew <- logK.to.OBIGT(logK, species, coeffs, T, P, name = "newpyrite", state = "cr", V = info(info("pyrite"))$V, npar = 5)
+# Get parameters of pyrite (from OBIGT database) and newpyrite (created by logK.to.OBIGT)
+oldpar <- info(info("pyrite"))
+newpar <- info(info("newpyrite"))
+expect_equal(convert(oldpar$G, "J"), newpar$G, info = info)
+expect_equal(convert(oldpar$S, "J"), newpar$S, info = info)
+expect_equal(convert(oldpar$a, "J"), newpar$a, info = info)
+expect_equal(convert(oldpar$b, "J"), newpar$b, info = info)
+expect_equal(convert(oldpar$c, "J"), newpar$c, info = info)
+
Modified: pkg/CHNOSZ/man/logK.to.OBIGT.Rd
===================================================================
--- pkg/CHNOSZ/man/logK.to.OBIGT.Rd 2025-01-22 06:16:21 UTC (rev 869)
+++ pkg/CHNOSZ/man/logK.to.OBIGT.Rd 2025-01-22 12:13:20 UTC (rev 870)
@@ -7,16 +7,19 @@
}
\usage{
- logK.to.OBIGT(logK, species, coeffs, T, P, npar = 3,
- optimize.omega = FALSE, tolerance = 0.05)
+ logK.to.OBIGT(logK, species, coeffs, T, P, name = NULL, state = "aq",
+ V = NA, npar = 3, optimize.omega = FALSE, tolerance = 0.05)
}
\arguments{
- \item{logK}{Values of \logK}
- \item{species}{Species in the formation reaction (the formed species must be last)}
- \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{logK}{numeric, values of \logK}
+ \item{species}{character, formulas of species in the formation reaction (the newly formed species must be last)}
+ \item{coeffs}{numeric, coefficients of the formation reaction}
+ \item{T}{numeric, temperature (units of \code{\link{T.units}})}
+ \item{P}{numeric, pressure (\samp{Psat} or numerical values with units of \code{\link{P.units}})}
+ \item{name}{character, name of newly formed species; if NULL it is set to the formula from the last member of \code{species}}
+ \item{state}{character, state: \samp{aq} or \samp{cr}}
+ \item{V}{numeric, molar volume of the formed species, if it is a mineral (for \code{state == "cr"})}
\item{npar}{numeric, number of parameters to use in fitting}
\item{optimize.omega}{logical, optimize the omega parameter (relevant for charged species)?}
\item{tolerance}{Tolerance for checking equivalence of input and calculated \logK values}
@@ -31,13 +34,15 @@
The function works as follows.
First, the properties of the known species are combined with \logK to calculate the standard Gibbs energy (G[T]) of the formed species at each value of \code{T} and \code{P}.
-Then, \code{\link{lm}} is used to fit one or more of the thermodynamic parameters \code{G}, \code{S}, \code{c1}, \code{c2}, and \code{omega} to the values of G[T].
-The first two of these parameters 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).
+Then, for \code{state == "aq"}, \code{\link{lm}} is used to fit one or more of the thermodynamic parameters \code{G}, \code{S}, \code{c1}, \code{c2}, and \code{omega} to the values of G[T].
+The first two of these parameters 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) (\samp{HKF} model in OBIGT).
+For \code{state == "cr"}, the supported fitting parameters are \code{G}, \code{S}, \code{a}, \code{b}, and \code{c}; the last three of these are the coefficient in the Maier-Kelley heat capacity equation (\samp{CGL} model in OBIGT).
The fitted parameters for the formed species are then added to OBIGT.
Finally, \code{\link{all.equal}} is used to test for approximate equivalence of the input values of \logK and calculated equilibrium constants; if the mean absolute difference exceeds \code{tolerance}, an error occurs.
To avoid overfitting, only the first three parameters (\code{G}, \code{S}, and \code{c1}) are used by default.
More parameters (up to 5) or fewer (down to 1) can be selected by changing \code{npar}.
+The default of 3 correspond to a constant heat capacity for both \code{aq} and \code{cr} states.
Volumetric parameters (\a1 to \a4) in the HKF equations currently aren't included, so the resulting fits are valid only at the input pressure values.
Independent of \code{npar}, the number of parameters used in the fit is not more than the number of data values (\emph{n}).
More information about the CHNOSZ-commits
mailing list