[CHNOSZ-commits] r869 - in pkg/CHNOSZ: . R demo inst man man/macros vignettes
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Jan 22 07:16:22 CET 2025
Author: jedick
Date: 2025-01-22 07:16:21 +0100 (Wed, 22 Jan 2025)
New Revision: 869
Added:
pkg/CHNOSZ/R/logK.to.OBIGT.R
pkg/CHNOSZ/man/logK.to.OBIGT.Rd
Removed:
pkg/CHNOSZ/R/logB.to.OBIGT.R
pkg/CHNOSZ/man/logB.to.OBIGT.Rd
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/NAMESPACE
pkg/CHNOSZ/R/examples.R
pkg/CHNOSZ/demo/00Index
pkg/CHNOSZ/demo/yttrium.R
pkg/CHNOSZ/inst/NEWS.Rd
pkg/CHNOSZ/man/add.OBIGT.Rd
pkg/CHNOSZ/man/examples.Rd
pkg/CHNOSZ/man/macros/macros.Rd
pkg/CHNOSZ/vignettes/custom_data.Rmd
pkg/CHNOSZ/vignettes/mklinks.sh
Log:
Rename logB.to.OBIGT() to logK.to.OBIGT()
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2025-01-17 10:26:48 UTC (rev 868)
+++ pkg/CHNOSZ/DESCRIPTION 2025-01-22 06:16:21 UTC (rev 869)
@@ -1,6 +1,6 @@
-Date: 2025-01-17
+Date: 2025-01-22
Package: CHNOSZ
-Version: 2.1.0-40
+Version: 2.1.0-41
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/NAMESPACE
===================================================================
--- pkg/CHNOSZ/NAMESPACE 2025-01-17 10:26:48 UTC (rev 868)
+++ pkg/CHNOSZ/NAMESPACE 2025-01-22 06:16:21 UTC (rev 869)
@@ -56,7 +56,7 @@
# added 20200716 or later
"mash", "mix", "rebalance",
# added 20220324
- "logB.to.OBIGT",
+ "logK.to.OBIGT",
# added 20220416
"rank.affinity",
# added 20220620
Modified: pkg/CHNOSZ/R/examples.R
===================================================================
--- pkg/CHNOSZ/R/examples.R 2025-01-17 10:26:48 UTC (rev 868)
+++ pkg/CHNOSZ/R/examples.R 2025-01-22 06:16:21 UTC (rev 869)
@@ -14,7 +14,7 @@
"diagram", "mosaic", "mix",
"buffer", "nonideal", "NaCl",
"add.protein", "ionize.aa", "EOSregress", "rank.affinity",
- "DEW", "logB.to.OBIGT", "stack_mosaic"
+ "DEW", "logK.to.OBIGT", "stack_mosaic"
)
plot.it <- FALSE
if(is.character(save.png))
Deleted: pkg/CHNOSZ/R/logB.to.OBIGT.R
===================================================================
--- pkg/CHNOSZ/R/logB.to.OBIGT.R 2025-01-17 10:26:48 UTC (rev 868)
+++ pkg/CHNOSZ/R/logB.to.OBIGT.R 2025-01-22 06:16:21 UTC (rev 869)
@@ -1,142 +0,0 @@
-# CHNOSZ/logB.to.OBIGT.R
-# 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, coeffs, T, P, 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("logB.to.OBIGT: forcing optimize.omega = FALSE because there are < 5 T values")
- optimize.omega <- FALSE
- }
- # We need five parameters to optimize omega 20221209
- if(optimize.omega & npar != 5) {
- message("logB.to.OBIGT: forcing optimize.omega = FALSE because npar != 5")
- 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
- name <- tail(species, 1)
- # 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))
- # 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)
-
- ## Get Gibbs energy of species from logB of reaction
- # Calculate T in Kelvin
- TK <- convert(T, "K")
- # logB gives values for ΔG°r of the reaction
- Gr <- convert(logB, "G", TK)
- # logK0 gives values for ΔG°r of the reaction with ΔG°f = 0 for the formed species
- Gr0 <- convert(logK0, "G", TK)
- # Calculate ΔG°f of the formed species
- Gf <- Gr - Gr0
-
- ## 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) {
-
- # 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
- 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, coeffs, T = T, P = P)$out$logK)
- # Calculate the mean absolute difference
- mad <- mean(abs(logK - logB))
- message(paste("logB.to.OBIGT: mean difference between logB (experimental) and logK (calculated) is", round(mad, 4)))
- # Check that calculated values are close to input values
- stopifnot(all.equal(logK, logB, tolerance = tolerance, scale = 1))
- # Return the species index in OBIGT
- ispecies
-
-}
Copied: pkg/CHNOSZ/R/logK.to.OBIGT.R (from rev 868, pkg/CHNOSZ/R/logB.to.OBIGT.R)
===================================================================
--- pkg/CHNOSZ/R/logK.to.OBIGT.R (rev 0)
+++ pkg/CHNOSZ/R/logK.to.OBIGT.R 2025-01-22 06:16:21 UTC (rev 869)
@@ -0,0 +1,142 @@
+# CHNOSZ/logK.to.OBIGT.R
+# Get thermodynamic parameters from formation constants (logK) 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)
+
+logK.to.OBIGT <- function(logK, species, coeffs, T, P, 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")
+ optimize.omega <- FALSE
+ }
+ # We need five parameters to optimize omega 20221209
+ if(optimize.omega & npar != 5) {
+ message("logK.to.OBIGT: forcing optimize.omega = FALSE because npar != 5")
+ 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
+ name <- tail(species, 1)
+ # 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))
+ # 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)
+
+ ## Get Gibbs energy of species from logK of reaction
+ # Calculate T in Kelvin
+ TK <- convert(T, "K")
+ # logK gives values for ΔG°r of the reaction
+ Gr <- convert(logK, "G", TK)
+ # logK0 gives values for ΔG°r of the reaction with ΔG°f = 0 for the formed species
+ Gr0 <- convert(logK0, "G", TK)
+ # Calculate ΔG°f of the formed species
+ Gf <- Gr - Gr0
+
+ ## 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) {
+
+ # 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
+ 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)
+ # 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
+
+}
Modified: pkg/CHNOSZ/demo/00Index
===================================================================
--- pkg/CHNOSZ/demo/00Index 2025-01-17 10:26:48 UTC (rev 868)
+++ pkg/CHNOSZ/demo/00Index 2025-01-22 06:16:21 UTC (rev 869)
@@ -32,5 +32,5 @@
Pourbaix Eh-pH diagram for Fe-O-H with equisolubility lines
E_coli Gibbs energy of biomass synthesis in E. coli
rank.affinity Affinity ranking for proteins in yeast nutrient limitation
-yttrium logB.to.OBIGT fits at 800 and 1000 bar and Y speciation in NaCl solution at varying pH
+yttrium logK.to.OBIGT fits at 800 and 1000 bar and Y speciation in NaCl solution at varying pH
uranyl Total (carbonate|sulfate)-pH diagrams for uranyl species
Modified: pkg/CHNOSZ/demo/yttrium.R
===================================================================
--- pkg/CHNOSZ/demo/yttrium.R 2025-01-17 10:26:48 UTC (rev 868)
+++ pkg/CHNOSZ/demo/yttrium.R 2025-01-22 06:16:21 UTC (rev 869)
@@ -5,12 +5,12 @@
library(CHNOSZ)
-# Function to add Y(III)-Cl species using logB values from Table 9 of Guan et al. (2020)
+# Function to add Y(III)-Cl species using logK values from Table 9 of Guan et al. (2020)
add.Y.species <- function(P, plot.it = FALSE) {
# Temperature
T <- c(25, seq(50, 500, 50))
if(P == 800) {
- logB <- list(
+ logK <- list(
c(1.04, 0.48, 0.13, 0.43, 1.12, 2.06, 3.21, 4.58, 6.26, 8.39, 11.14),
c(-9.14, -7.34, -4.14, -1.37, 1.12, 3.46, 5.78, 8.24, 11.05, 14.48, 18.77),
c(-14, -11.48, -7.06, -3.25, 0.14, 3.32, 6.45, 9.74, 13.47, 18.01, 23.65),
@@ -17,13 +17,13 @@
c(-15.94, -13.2, -8.39, -4.27, -0.61, 2.79, 6.15, 9.67, 13.63, 18.45, 24.41)
)
} else if(P == 1000) {
- logB <- list(
+ logK <- list(
c(1.13, 0.54, 0.16, 0.43, 1.09, 2, 3.1, 4.39, 5.9, 7.69, 9.85),
c(-9.33, -7.51, -4.3, -1.55, 0.9, 3.18, 5.4, 7.68, 10.15, 12.94, 16.16),
c(-14.24, -11.71, -7.27, -3.49, -0.14, 2.95, 5.95, 9.01, 12.3, 16, 20.25),
c(-16.19, -13.43, -8.62, -4.52, -0.91, 2.41, 5.62, 8.89, 12.4, 16.33, 20.84)
)
- } else stop("logB values for P =", P, "are not available here")
+ } else stop("logK values for P =", P, "are not available here")
# Define species and coefficients in formation reactions
species <- list(
c("Y+3", "Cl-", "YCl+2"),
@@ -38,7 +38,7 @@
c(-1, -4, 1)
)
# Fit the formation constants to thermodynamic parameters and add them to OBIGT
- for(i in 1:4) logB.to.OBIGT(logB[[i]], species[[i]], coeffs[[i]], T = T, P = P, tolerance = 0.6, npar = 5)
+ for(i in 1:4) logK.to.OBIGT(logK[[i]], species[[i]], coeffs[[i]], T = T, P = P, tolerance = 0.6, npar = 5)
# Plot the given and fitted values
if(plot.it) {
par(mfrow = c(2, 2))
@@ -45,9 +45,9 @@
for(i in 1:4) {
sres <- subcrt(species[[i]], coeffs[[i]], T = T, P = P)
plot(T, sres$out$logK, type = "l", xlab = axis.label("T"), ylab = axis.label("logK"))
- points(T, logB[[i]], pch = 19)
+ points(T, logK[[i]], pch = 19)
title(describe.reaction(sres$reaction))
- if(i == 1) legend("topleft", c("Guan et al. (2020)", "logB.to.OBIGT()"), pch = c(19, NA), lty = c(NA, 1))
+ if(i == 1) legend("topleft", c("Guan et al. (2020)", "logK.to.OBIGT()"), pch = c(19, NA), lty = c(NA, 1))
legend("bottomright", paste(P, "bar"), bty = "n")
}
}
Modified: pkg/CHNOSZ/inst/NEWS.Rd
===================================================================
--- pkg/CHNOSZ/inst/NEWS.Rd 2025-01-17 10:26:48 UTC (rev 868)
+++ pkg/CHNOSZ/inst/NEWS.Rd 2025-01-22 06:16:21 UTC (rev 869)
@@ -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-40 (2025-01-17)}{
+\section{Changes in CHNOSZ version 2.1.0-41 (2025-01-22)}{
\subsection{OBIGT DEFAULT DATA}{
\itemize{
@@ -122,6 +122,8 @@
\samp{m_Clminus}, \samp{m_NaCl0}, and \samp{m_HCl0} (molalities of
aqueous species).
+ \item Rename \code{logB.to.OBIGT()} to \code{logK.to.OBIGT()}.
+
}
}
Modified: pkg/CHNOSZ/man/add.OBIGT.Rd
===================================================================
--- pkg/CHNOSZ/man/add.OBIGT.Rd 2025-01-17 10:26:48 UTC (rev 868)
+++ pkg/CHNOSZ/man/add.OBIGT.Rd 2025-01-22 06:16:21 UTC (rev 869)
@@ -63,7 +63,7 @@
The values returned (\code{\link{invisible}}-y) are the indices of the added and/or modified species.
}
-\seealso{ \code{\link{thermo}} (description of OBIGT), \code{\link{mod.buffer}} (modify buffer definitions), \code{\link{logB.to.OBIGT}} (fit thermodynamic parameters to formation constants) }
+\seealso{ \code{\link{thermo}} (description of OBIGT), \code{\link{mod.buffer}} (modify buffer definitions), \code{\link{logK.to.OBIGT}} (fit thermodynamic parameters to formation constants) }
\examples{\dontshow{reset()}
## Modify an existing species (not real properties)
Modified: pkg/CHNOSZ/man/examples.Rd
===================================================================
--- pkg/CHNOSZ/man/examples.Rd 2025-01-17 10:26:48 UTC (rev 868)
+++ pkg/CHNOSZ/man/examples.Rd 2025-01-22 06:16:21 UTC (rev 869)
@@ -65,7 +65,7 @@
\item{Pourbaix}{Eh-pH diagram for Fe-O-H with equisolubility lines (Pourbaix, 1974)}
\item{E_coli}{Gibbs energy of biomass synthesis in \emph{E. coli} (LaRowe and Amend, 2016)}
\item{rank.affinity}{Affinity ranking for proteins in yeast nutrient limitation (data from Tai et al., 2005)}
- \item{yttrium}{\code{\link{logB.to.OBIGT}} fits at 800 and 1000 bar and Y speciation in \code{\link{NaCl}} solution at varying pH (Guan et al., 2020)}
+ \item{yttrium}{\code{\link{logK.to.OBIGT}} fits at 800 and 1000 bar and Y speciation in \code{\link{NaCl}} solution at varying pH (Guan et al., 2020)}
\item{uranyl}{Total (carbonate|sulfate)-pH diagrams for uranyl species (Migdisov et al., 2024)}
\item{sum_S}{Summed molality of S species and solubility contours for iron and gold (Skirrow and Walshe, 2002)}
}
Deleted: pkg/CHNOSZ/man/logB.to.OBIGT.Rd
===================================================================
--- pkg/CHNOSZ/man/logB.to.OBIGT.Rd 2025-01-17 10:26:48 UTC (rev 868)
+++ pkg/CHNOSZ/man/logB.to.OBIGT.Rd 2025-01-22 06:16:21 UTC (rev 869)
@@ -1,120 +0,0 @@
-\encoding{UTF-8}
-\name{logB.to.OBIGT}
-\alias{logB.to.OBIGT}
-\title{Fit thermodynamic parameters to formation constants (\logB)}
-\description{
-Fit thermodynamic parameters to experimental formation constants for an aqueous species and add the parameters to OBIGT.
-}
-
-\usage{
- logB.to.OBIGT(logB, species, coeffs, T, P, npar = 3,
- optimize.omega = FALSE, tolerance = 0.05)
-}
-
-\arguments{
- \item{logB}{Values of \logB}
- \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{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 \logB values}
-}
-
-\details{
-
-This function updates 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 formation constants for an aqueous complex (\logB).
-The formation reaction is defined by \code{species} and \code{coeffs}.
-All species in the formation reaction must be present in OBIGT except for the \emph{last} species, which is the newly formed species.
-
-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, \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).
-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 \logB 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}.
-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}).
-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.
-
-The value of \omega is a constant in the revised HKF equations for uncharged species, but for charged species, it is a function of \T and \P as described by the \dQuote{g function} (Shock et al., 1992).
-An optimization step is available to refine the value of \code{omega} at 25 \degC and 1 bar for charged species taking its temperature dependence into account for the fitting.
-However, representative datasets are not well represented by variable \code{omega} (see first example), so this step is skipped by default.
-Furthermore, \code{logB.to.OBIGT} by default also sets the \code{z} parameter in OBIGT to 0; this enforces a constant \omega for charged species in calculations with \code{\link{subcrt}} (note that this is a parameter for the HKF equations and does not affect reaction balancing).
-Set \code{optimize.omega} to TRUE to override the defaults and activate variable \omega for charged species; this takes effect only if \code{npar == 5} and \emph{n} > 5.
-
-}
-
-\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. 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)
-T <- c(120, 150, 200, 250, 300)
-P <- "Psat"
-species <- c("Co+2", "HS-", "CoHS+")
-coeffs <- c(-1, -1, 1)
-opar <- par(mfrow = c(2, 1))
-for(o.o in c(TRUE, FALSE)) {
- # Fit the parameters with or without variable omega
- inew <- logB.to.OBIGT(logB, species, coeffs, T, P, npar = 5, optimize.omega = o.o)
- # 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 = ",o.o,")"))
- legend("top", legend, pch = c(19, NA), lty = c(0, 1), col = c(1, 4),
- pt.cex = 2, bg = "#FFFFFFB0")
-}
-par(opar)
-# NB. Optimizing omega leads to unphysical oscillations in the logK (first plot)
-
-## ZnCl+ from Mei et al. (2015)
-# Values for 5000 bar calculated with 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+")
-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, coeffs, T, P, npar = 5)
-# 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, 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)
-legend <- c("Mei et al. (2015)", "logB.to.OBIGT()")
-legend("topleft", legend, pch = c(19, NA), lty = c(0, 1), col = c(1, 4), pt.cex = 2)
-}
-
-\references{
-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}
-}
Copied: pkg/CHNOSZ/man/logK.to.OBIGT.Rd (from rev 868, pkg/CHNOSZ/man/logB.to.OBIGT.Rd)
===================================================================
--- pkg/CHNOSZ/man/logK.to.OBIGT.Rd (rev 0)
+++ pkg/CHNOSZ/man/logK.to.OBIGT.Rd 2025-01-22 06:16:21 UTC (rev 869)
@@ -0,0 +1,120 @@
+\encoding{UTF-8}
+\name{logK.to.OBIGT}
+\alias{logK.to.OBIGT}
+\title{Fit thermodynamic parameters to formation constants (\logK)}
+\description{
+Fit thermodynamic parameters to experimental formation constants for an aqueous species and add the parameters to OBIGT.
+}
+
+\usage{
+ logK.to.OBIGT(logK, species, coeffs, T, P, 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{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}
+}
+
+\details{
+
+This function updates the \code{\link{OBIGT}} thermodynamic database with parameters fit to formation constants of aqueous species as a function of temperature.
+The \code{logK} argument should have decimal logarithm of formation constants for an aqueous complex (\logK).
+The formation reaction is defined by \code{species} and \code{coeffs}.
+All species in the formation reaction must be present in OBIGT except for the \emph{last} species, which is the newly formed species.
+
+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).
+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}.
+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}).
+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 \logK would be fit only with \code{G}, implying that computed values of G[T] have no temperature dependence.
+
+The value of \omega is a constant in the revised HKF equations for uncharged species, but for charged species, it is a function of \T and \P as described by the \dQuote{g function} (Shock et al., 1992).
+An optimization step is available to refine the value of \code{omega} at 25 \degC and 1 bar for charged species taking its temperature dependence into account for the fitting.
+However, representative datasets are not well represented by variable \code{omega} (see first example), so this step is skipped by default.
+Furthermore, \code{logK.to.OBIGT} by default also sets the \code{z} parameter in OBIGT to 0; this enforces a constant \omega for charged species in calculations with \code{\link{subcrt}} (note that this is a parameter for the HKF equations and does not affect reaction balancing).
+Set \code{optimize.omega} to TRUE to override the defaults and activate variable \omega for charged species; this takes effect only if \code{npar == 5} and \emph{n} > 5.
+
+}
+
+\seealso{
+\code{logK.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. 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)
+logK <- c(6.24, 6.02, 5.84, 5.97, 6.52)
+T <- c(120, 150, 200, 250, 300)
+P <- "Psat"
+species <- c("Co+2", "HS-", "CoHS+")
+coeffs <- c(-1, -1, 1)
+opar <- par(mfrow = c(2, 1))
+for(o.o in c(TRUE, FALSE)) {
+ # Fit the parameters with or without variable omega
+ inew <- logK.to.OBIGT(logK, species, coeffs, T, P, npar = 5, optimize.omega = o.o)
+ # Print the new database entry
+ info(inew)
+ # Plot experimental logK
+ plot(T, logK, "n", c(100, 320), c(5.8, 6.8),
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/chnosz -r 869
More information about the CHNOSZ-commits
mailing list