[CHNOSZ-commits] r710 - in pkg/CHNOSZ: . R man man/macros vignettes
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Mar 24 05:21:32 CET 2022
Author: jedick
Date: 2022-03-24 05:21:31 +0100 (Thu, 24 Mar 2022)
New Revision: 710
Added:
pkg/CHNOSZ/R/logB_to_OBIGT.R
pkg/CHNOSZ/man/logB_to_OBIGT.Rd
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/NAMESPACE
pkg/CHNOSZ/R/add.OBIGT.R
pkg/CHNOSZ/man/add.OBIGT.Rd
pkg/CHNOSZ/man/macros/macros.Rd
pkg/CHNOSZ/vignettes/customizing.Rmd
Log:
Add logB_to_OBIGT() (fit formation constants of aqueous species)
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2022-03-24 00:37:09 UTC (rev 709)
+++ pkg/CHNOSZ/DESCRIPTION 2022-03-24 04:21:31 UTC (rev 710)
@@ -1,6 +1,6 @@
Date: 2022-03-24
Package: CHNOSZ
-Version: 1.4.3-2
+Version: 1.4.3-3
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 2022-03-24 00:37:09 UTC (rev 709)
+++ pkg/CHNOSZ/NAMESPACE 2022-03-24 04:21:31 UTC (rev 710)
@@ -57,7 +57,9 @@
"CHNOSZ", "thermo", "reset", "OBIGT", "retrieve", "AD", "moles",
"lNaCl", "lS", "lT", "lP", "lTP", "lex",
# added 20200716 or later
- "mash", "mix", "rebalance"
+ "mash", "mix", "rebalance",
+# added 20220324
+ "logB_to_OBIGT"
)
# Load shared objects
@@ -73,7 +75,7 @@
"plot.new", "plot.window", "points", "rect", "text", "title")
importFrom("stats", "D", "as.formula", "cor", "lm", "loess.smooth",
"na.omit", "predict.lm", "qqline", "qqnorm", "sd", "splinefun",
- "uniroot", "median", "predict", "smooth.spline")
+ "uniroot", "median", "predict", "smooth.spline", "optimize")
importFrom("utils", "browseURL", "capture.output", "combn", "demo",
"example", "head", "installed.packages", "read.csv", "tail",
"write.csv", "write.table", "read.table", "packageDescription")
Modified: pkg/CHNOSZ/R/add.OBIGT.R
===================================================================
--- pkg/CHNOSZ/R/add.OBIGT.R 2022-03-24 00:37:09 UTC (rev 709)
+++ pkg/CHNOSZ/R/add.OBIGT.R 2022-03-24 04:21:31 UTC (rev 710)
@@ -6,23 +6,23 @@
#source("util.data.R")
mod.OBIGT <- function(...) {
- # add or modify species in thermo()$OBIGT
+ # Add or modify species in thermo()$OBIGT
thermo <- get("thermo", CHNOSZ)
- # the names and values are in the arguments
- # this works for providing arguments via do.call
+ # The names and values are in the arguments
+ # This works for providing arguments via do.call
args <- list(...)
- # this is needed if we are called with a list as the actual argument
+ # This is needed if we are called with a list as the actual argument
if(is.list(args[[1]])) args <- args[[1]]
if(length(args) < 2) stop("please supply at least a species name and a property to update")
if(is.null(names(args))) stop("all arguments after the first should be named")
if(any(tail(nchar(names(args)), -1)==0)) stop("all arguments after the first should be named")
- # if the first argument is numeric, it's the species index
+ # If the first argument is numeric, it's the species index
if(is.numeric(args[[1]][1])) {
ispecies <- args[[1]]
args <- args[-1]
speciesname <- info(ispecies, check.it = FALSE)$name
} else {
- # if the name of the first argument is missing, assume it's the species name
+ # If the name of the first argument is missing, assume it's the species name
if(names(args)[1]=="") names(args)[1] <- "name"
speciesname <- args$name
# search for this species, use check.protein=FALSE to avoid infinite loop when adding proteins
@@ -32,56 +32,56 @@
else ispecies <- suppressMessages(mapply(info.character,
species=args$name, check.protein=FALSE, SIMPLIFY=TRUE, USE.NAMES=FALSE))
}
- # the column names of thermo()$OBIGT, split at the "."
+ # The column names of thermo()$OBIGT, split at the "."
cnames <- c(do.call(rbind, strsplit(colnames(thermo$OBIGT), ".", fixed=TRUE)), colnames(thermo$OBIGT))
- # the columns we are updating
+ # The columns we are updating
icol <- match(names(args), cnames)
if(any(is.na(icol))) stop(paste("properties not in thermo$OBIGT:", paste(names(args)[is.na(icol)], collapse=" ")) )
- # the column numbers for properties that matched after the split
+ # The column numbers for properties that matched after the split
icol[icol > 42] <- icol[icol > 42] - 42
icol[icol > 21] <- icol[icol > 21] - 21
- # which species are new and which are old
+ # Which species are new and which are old
inew <- which(is.na(ispecies))
iold <- which(!is.na(ispecies))
- # the arguments as data frame
+ # The arguments as data frame
args <- data.frame(args, stringsAsFactors=FALSE)
if(length(inew) > 0) {
- # the right number of blank rows of thermo()$OBIGT
+ # The right number of blank rows of thermo()$OBIGT
newrows <- thermo$OBIGT[1:length(inew), ]
- # if we don't know something it's NA
+ # If we don't know something it's NA
newrows[] <- NA
- # put in a default state
+ # Put in a default state
newrows$state <- thermo$opt$state
- # the formula defaults to the name
+ # The formula defaults to the name
newrows$formula <- args$name[inew]
- # the units should also be set 20190530
+ # The units should also be set 20190530
newrows$E_units <- thermo$opt$E.units
- # fill in the columns
+ # Fill in the columns
newrows[, icol] <- args[inew, ]
- # now check the formulas
+ # Now check the formulas
e <- tryCatch(makeup(newrows$formula), error=function(e) e)
if(inherits(e, "error")) {
warning("please supply a valid chemical formula as the species name or in the 'formula' argument")
- # transmit the error from makeup
+ # Transmit the error from makeup
stop(e)
}
- # assign to thermo()$OBIGT
+ # Assign to thermo()$OBIGT
thermo$OBIGT <- rbind(thermo$OBIGT, newrows)
rownames(thermo$OBIGT) <- NULL
assign("thermo", thermo, CHNOSZ)
- # update ispecies
+ # Update ispecies
ntotal <- nrow(thermo$OBIGT)
ispecies[inew] <- (ntotal-length(inew)+1):ntotal
- # inform user
+ # Inform user
message(paste("mod.OBIGT: added ", newrows$name, "(", newrows$state, ")", " with energy units of ", newrows$E_units, sep="", collapse="\n"))
}
if(length(iold) > 0) {
- # loop over species
+ # Loop over species
for(i in 1:length(iold)) {
- # the old values and the state
+ # The old values and the state
oldprop <- thermo$OBIGT[ispecies[iold[i]], icol]
state <- thermo$OBIGT$state[ispecies[iold[i]]]
- # tell user if they're the same, otherwise update the data entry
+ # Tell user if they're the same, otherwise update the data entry
if(isTRUE(all.equal(oldprop, args[iold[i], ], check.attributes=FALSE)))
message("mod.OBIGT: no change for ", speciesname[iold[i]], "(", state, ")")
else {
@@ -95,29 +95,29 @@
}
add.OBIGT <- function(file, species=NULL, force=TRUE) {
- # add/replace entries in thermo$OBIGT from values saved in a file
- # only replace if force==TRUE
+ # Add/replace entries in thermo$OBIGT from values saved in a file
+ # Only replace if force==TRUE
thermo <- get("thermo", CHNOSZ)
to1 <- thermo$OBIGT
id1 <- paste(to1$name,to1$state)
- # we match system files with the file suffixes (.csv) removed
+ # We match system files with the file suffixes (.csv) removed
sysfiles <- dir(system.file("extdata/OBIGT/", package="CHNOSZ"))
sysnosuffix <- sapply(strsplit(sysfiles, "\\."), "[", 1)
isys <- match(file, sysnosuffix)
if(!is.na(isys)) file <- system.file(paste0("extdata/OBIGT/", sysfiles[isys]), package="CHNOSZ")
# else {
-# # we also match single system files with the state suffix removed
+# # We also match single system files with the state suffix removed
# # (e.g. "DEW" for "DEW_aq", but not "organic" because we have "organic_aq", "organic_cr", etc.)
# sysnostate <- sapply(strsplit(sysnosuffix, "_"), "[", 1)
# isys <- which(file==sysnostate)
# if(length(isys)==1) file <- system.file(paste0("extdata/OBIGT/", sysfiles[isys]), package="CHNOSZ")
# }
- # read data from the file
+ # Tead data from the file
to2 <- read.csv(file, as.is=TRUE)
- # add E_units column if it's missing 20190529
+ # Add E_units column if it's missing 20190529
if(!"E_units" %in% colnames(to2)) to2 <- data.frame(to2[, 1:7], E_units = "cal", to2[, 8:20], stringsAsFactors = FALSE)
Etxt <- paste(unique(to2$E_units), collapse = " and ")
- # load only selected species if requested
+ # Load only selected species if requested
if(!is.null(species)) {
idat <- match(species, to2$name)
ina <- is.na(idat)
@@ -125,17 +125,17 @@
else stop(paste("file", file, "doesn't have", paste(species[ina], collapse=", ")))
}
id2 <- paste(to2$name,to2$state)
- # check if the data is compatible with thermo$OBIGT
+ # Check if the data is compatible with thermo$OBIGT
tr <- tryCatch(rbind(to1, to2), error = identity)
if(inherits(tr, "error")) stop(paste(file, "is not compatible with thermo$OBIGT data table."))
- # match the new species to existing ones
+ # Match the new species to existing ones
does.exist <- id2 %in% id1
ispecies.exist <- na.omit(match(id2, id1))
nexist <- sum(does.exist)
- # keep track of the species we've added
+ # Keep track of the species we've added
inew <- numeric()
if(force) {
- # replace existing entries
+ # Replace existing entries
if(nexist > 0) {
to1[ispecies.exist, ] <- to2[does.exist, ]
to2 <- to2[!does.exist, ]
@@ -142,20 +142,20 @@
inew <- c(inew, ispecies.exist)
}
} else {
- # ignore any new entries that already exist
+ # Ignore any new entries that already exist
to2 <- to2[!does.exist, ]
nexist <- 0
}
- # add new entries
+ # Add new entries
if(nrow(to2) > 0) {
to1 <- rbind(to1, to2)
inew <- c(inew, (length(id1)+1):nrow(to1))
}
- # commit the change
+ # Commit the change
thermo$OBIGT <- to1
rownames(thermo$OBIGT) <- 1:nrow(thermo$OBIGT)
assign("thermo", thermo, CHNOSZ)
- # give the user a message
+ # Give the user a message
message("add.OBIGT: read ", length(does.exist), " rows; made ",
nexist, " replacements, ", nrow(to2), " additions [energy units: ", Etxt, "]")
#message("add.OBIGT: use OBIGT() or reset() to restore default database")
Added: pkg/CHNOSZ/R/logB_to_OBIGT.R
===================================================================
--- pkg/CHNOSZ/R/logB_to_OBIGT.R (rev 0)
+++ pkg/CHNOSZ/R/logB_to_OBIGT.R 2022-03-24 04:21:31 UTC (rev 710)
@@ -0,0 +1,66 @@
+# CHNOSZ/logB_to_OBIGT.R
+# Fit G, S, and Cp to formation constants (logB) as a function of temperature
+# 20220324 jmd
+
+logB_to_OBIGT <- function(logB, species, coeff, T, P, tolerance = 0.05) {
+
+ ## 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)
+ # Check if the species is already present in OBIGT and delete its properties
+ iname <- suppressMessages(info(name))
+ if(!is.na(iname)) {
+ message(paste("logB_to_OBIGT: deleting existing", name, "from OBIGT"))
+ thermo <- get("thermo", CHNOSZ)
+ thermo$OBIGT <- thermo$OBIGT[-iname, ]
+ assign("thermo", thermo, CHNOSZ)
+ }
+ # Add the formed species to the thermodynamic database with ΔG°f = 0 at all temperatures
+ suppressMessages(mod.OBIGT(name, G = 0, S = 0, Cp = 0))
+ # 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)
+
+ ## 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
+
+ ## 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)
+ }
+ # 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]]
+
+ ## Update the thermodynamic parameters of the formed species
+ ispecies <- do.call(mod.OBIGT, list(name = name, G = G, S = S, Cp = Cp))
+ # 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)
+ # 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)))
+ # 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
+
+}
+
Modified: pkg/CHNOSZ/man/add.OBIGT.Rd
===================================================================
--- pkg/CHNOSZ/man/add.OBIGT.Rd 2022-03-24 00:37:09 UTC (rev 709)
+++ pkg/CHNOSZ/man/add.OBIGT.Rd 2022-03-24 04:21:31 UTC (rev 710)
@@ -55,7 +55,7 @@
The values returned (\code{\link{invisible}}-y) are the indices of the added and/or modified species.
}
-\seealso{ \code{\link{thermo}}, \code{\link{util.data}}, \code{\link{mod.buffer}} }
+\seealso{ \code{\link{thermo}}, \code{\link{util.data}}, \code{\link{mod.buffer}}, \code{\link{logB_to_OBIGT}} }
\examples{\dontshow{reset()}
## Modify an existing species (not real properties)
Added: pkg/CHNOSZ/man/logB_to_OBIGT.Rd
===================================================================
--- pkg/CHNOSZ/man/logB_to_OBIGT.Rd (rev 0)
+++ pkg/CHNOSZ/man/logB_to_OBIGT.Rd 2022-03-24 04:21:31 UTC (rev 710)
@@ -0,0 +1,68 @@
+\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, coeff, T, P, 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{T}{Temperature (units of \code{\link{T.units}})}
+ \item{P}{Temperature (\samp{Psat} or numerical values with units of \code{\link{P.units}})}
+ \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 last species, which is 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.
+
+The values of G[T] are fit using the thermodynamic parameters \emph{G}, \emph{S}, and \emph{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 * log(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 log 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 fitted \emph{G}, \emph{S}, and \emph{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 differences exceed \code{tolerance}, the code generates an error.
+
+}
+
+\seealso{
+See \code{\link{add.OBIGT}} for more information about modifying the database.
+}
+
+
+\examples{
+# Formation constants of 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"))
+points(T, logB, pch = 19, cex = 2)
+# Add 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)
+}
+
+\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}
+}
Modified: pkg/CHNOSZ/man/macros/macros.Rd
===================================================================
--- pkg/CHNOSZ/man/macros/macros.Rd 2022-03-24 00:37:09 UTC (rev 709)
+++ pkg/CHNOSZ/man/macros/macros.Rd 2022-03-24 04:21:31 UTC (rev 710)
@@ -63,3 +63,6 @@
\newcommand{\c1}{\ifelse{latex}{\eqn{c_1}}{\ifelse{html}{\out{<I>c</I><sub>1</sub>}}{c1}}}
\newcommand{\c2}{\ifelse{latex}{\eqn{c_2}}{\ifelse{html}{\out{<I>c</I><sub>2</sub>}}{c2}}}
\newcommand{\omega}{\ifelse{latex}{\eqn{\omega}}{\ifelse{html}{\out{ω}}{ω}}}
+
+% For logB_to_OBIGT.Rd
+\newcommand{\logB}{\ifelse{latex}{\eqn{\log \beta}}{\ifelse{html}{\out{log β}}{log β}}}
Modified: pkg/CHNOSZ/vignettes/customizing.Rmd
===================================================================
--- pkg/CHNOSZ/vignettes/customizing.Rmd 2022-03-24 00:37:09 UTC (rev 709)
+++ pkg/CHNOSZ/vignettes/customizing.Rmd 2022-03-24 04:21:31 UTC (rev 710)
@@ -170,7 +170,7 @@
Note that a missing (blank) value in the file is treated as NA.
- Unknown values for character values (Columns 1-8) should be NA.
-- Also, if you have only two of G, H, and S, then the missing one should be NA, not 0. This is because CHNOSZ "knows" the equation DG = DGH - TDS, and will handle a single missing one of those parameters.
+- 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 HKF or CGL parameters aren't present: thermodynamic properties *can't* be calculated at T > 25 °C (NAs will propagate to higher T).
@@ -187,7 +187,7 @@
OOM scaling is also applied to the CGL parameters.
See `?thermo` for details of the OOM scaling.
-`r info_` provides a way for the user to query the OBIGT database.
+`r info_` provides a simple user interface to the OBIGT database and is called by other functions in CHNOSZ to retrieve unscaled values from the database.
This is a summary of its main features:
- 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.
@@ -210,7 +210,7 @@
<summary><span style="background-color: #DDDDDD">* Click to expand *</span></summary>
Let's look at some minerals first.
-Use `r info_` to get the species indices (i.e. rownumbers) in OBIGT, and then the "raw" data (including OOM scaling and any NA values).
+First use `r info_` to get the species indices (i.e. rownumbers) in OBIGT, then pull out the "raw" data (including OOM scaling and any NA values).
```{r icr, message = FALSE}
icr <- info(c("orpiment,amorphous", "arsenic,alpha", "tin"))
thermo()$OBIGT[icr, ]
@@ -245,7 +245,7 @@
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.
-`r info_` actually does this to fills 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:
+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:
```{r info_.tin}
info(info("tin"))
More information about the CHNOSZ-commits
mailing list