[CHNOSZ-commits] r697 - in pkg/CHNOSZ: . R inst inst/tinytest man vignettes
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sun Feb 6 15:27:27 CET 2022
Author: jedick
Date: 2022-02-06 15:27:27 +0100 (Sun, 06 Feb 2022)
New Revision: 697
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/R/AkDi.R
pkg/CHNOSZ/inst/NEWS.Rd
pkg/CHNOSZ/inst/tinytest/test-AkDi.R
pkg/CHNOSZ/man/eos.Rd
pkg/CHNOSZ/vignettes/multi-metal.Rmd
Log:
AkDi() now calculates heat capacity and volume
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2022-02-06 11:33:15 UTC (rev 696)
+++ pkg/CHNOSZ/DESCRIPTION 2022-02-06 14:27:27 UTC (rev 697)
@@ -1,6 +1,6 @@
Date: 2022-02-06
Package: CHNOSZ
-Version: 1.4.1-22
+Version: 1.4.1-23
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/AkDi.R
===================================================================
--- pkg/CHNOSZ/R/AkDi.R 2022-02-06 11:33:15 UTC (rev 696)
+++ pkg/CHNOSZ/R/AkDi.R 2022-02-06 14:27:27 UTC (rev 697)
@@ -12,25 +12,37 @@
# R expressed in volume units
RV <- 10 * R # cm3 bar K-1 mol-1
- # Calculate properties of H2O
- # Moved these calculations to unexported functions
+ # Calculate H2O fugacity and derivatives of density
+ # These calculations are done through (unexported) functions
# to be able to test their output in test-AkDi.R 20220206
f1 <- .f1(T, P, isPsat)
- rho1 <- .rho1(T, P)
drho1_dT <- mapply(.drho1_dT, T, P)
drho1_dP <- mapply(.drho1_dP, T, P)
d2rho1_dT2 <- mapply(.d2rho1_dT2, T, P)
+
+ # Calculate other properties of H2O solvent
+ waterTP <- water(c("rho", "S", "Cp", "V"), T = T, P = P)
+ # Density (g cm-3)
+ rho1 <- waterTP$rho / 1000
+ # Entropy (dimensionless)
+ S1 <- waterTP$S / 1.9872
+ # Heat capacity (dimensionless)
+ Cp1 <- waterTP$Cp / 1.9872
+ # Volume (cm3 mol-1)
+ V1 <- waterTP$V
+ # Calculate properties of ideal H2O gas
S1_g <- sapply(T, .S1_g)
- S1 <- mapply(.S1, T, P, MoreArgs = list(isPsat = isPsat))
+ Cp1_g <- sapply(T, .Cp1_g)
-#ADstuff <<- list(f1 = f1, drho1_dT = drho1_dT, drho1_dP = drho1_dP, d2rho1_dT2 = d2rho1_dT2, S1_g = S1_g, S1 = S1)
-
# Initialize a list for the output
out <- list()
# Loop over species
nspecies <- nrow(parameters)
for(i in seq_len(nspecies)) {
+
+ # Get thermodynamic parameters for the gas and calculate properties at T, P
PAR <- parameters[i, ]
+ gasprops <- subcrt(PAR$name, "gas", T = T, P = P, convert = FALSE)$out[[1]]
# Send a message
message("AkDi: Akinfiev-Diamond model for ", PAR$name, " gas to aq")
# Start with an NA-filled data frame
@@ -42,9 +54,7 @@
if(property[[j]] == "G") {
# Get gas properties (J mol-1)
- G_gas <- subcrt(PAR$name, "gas", T = T, P = P, convert = FALSE)$out[[1]]$G
- # TODO: Does this work if E.units is cal or J?
- G_gas <- convert(G_gas, "J", T = T)
+ G_gas <- convert(gasprops$G, "J", T = T)
# Calculate G_hyd (J mol-1)
G_hyd <- R*T * ( -log(NW) + (1 - PAR$xi) * log(f1) + PAR$xi * log(RV * T * rho1 / MW) + rho1 * (PAR$a + PAR$b * (1000/T)^0.5) )
# Calculate the chemical potential (J mol-1)
@@ -57,8 +67,7 @@
if(property[[j]] == "S") {
# Get S_gas
- S_gas <- subcrt(PAR$name, "gas", T = T, P = P, convert = FALSE)$out[[1]]$S
- S_gas <- convert(S_gas, "J", T = T)
+ S_gas <- convert(gasprops$S, "J", T = T)
# Calculate S_hyd
S_hyd <- R * (
(1 - PAR$xi) * (S1 - S1_g)
@@ -74,7 +83,34 @@
myprops$S <- S
}
+ if(property[[j]] == "Cp") {
+ # Get Cp_gas
+ Cp_gas <- convert(gasprops$Cp, "J", T = T)
+ # Calculate Cp_hyd
+ Cp_hyd <- R * (
+ (1 - PAR$xi) * (Cp1 - Cp1_g)
+ - (PAR$xi + 2 * PAR$xi * T / rho1 * drho1_dT - PAR$xi * T^2 / rho1^2 * drho1_dT^2 + PAR$xi * T^2 / rho1 * d2rho1_dT2)
+ ) - R*T * (
+ PAR$a * (2 * drho1_dT + T * d2rho1_dT2)
+ + PAR$b * (-0.25 * 10^1.5 * T^-1.5 * rho1 + 10^1.5 * T^-0.5 * drho1_dT + 10^1.5 * T^0.5 * d2rho1_dT2)
+ )
+ Cp <- Cp_gas + Cp_hyd
+ Cp <- convert(Cp, "cal", T = T)
+ myprops$Cp <- Cp
+ }
+
+ if(property[[j]] == "V") {
+ # Get V_gas
+ V_gas <- 0
+ # Calculate V_hyd
+ V_hyd <- V1 * (1 - PAR$xi) + PAR$xi * RV * T / rho1 * drho1_dP + RV * T * drho1_dP * (PAR$a + PAR$b * (1000/T)^0.5)
+ V <- V_gas + V_hyd
+ myprops$V <- V
+ }
+
}
+ # Calculate enthalpy (NOTE: this is in calories) 20220206
+ myprops$H <- myprops$G - 298.15 * entropy(PAR$formula) + T * myprops$S
out[[i]] <- myprops
}
out
@@ -116,7 +152,9 @@
dP <- 0.1
P1 <- P - dP
P2 <- P + dP
- rho1 <- .rho1(T, c(P1, P2))
+ # Subtract 1 degC from T so the derivative doesn't blow up when P = Psat at T > 100 degC
+ # TODO: Is there a better way?
+ rho1 <- .rho1(T - 1, c(P1, P2))
diff(rho1) / (P2 - P1)
}
@@ -125,7 +163,8 @@
dT <- 0.1
# Calculate density at five temperature values
Tval <- seq(T - 2 * dT, T + 2 * dT, dT)
- rho1 <- .rho1(Tval, P)
+ # TODO: Is there a better way to calculate the partial derivative for P = Psat?
+ rho1 <- .rho1(Tval, P + 1)
# https://stackoverflow.com/questions/11081069/calculate-the-derivative-of-a-data-function-in-r
spl <- smooth.spline(Tval, rho1)
# The second derivative of the fitted spline function at the third point (i.e., T)
@@ -137,9 +176,7 @@
.Fortran(C_ideal2, T, 0, 0, 0, 0, 0, 0, 0)[[4]]
}
-.S1 <- function(T, P, isPsat) {
- # Entropy of solvent H2O (dimensionless)
- if(isPsat) S1 <- water("S", T = T, P = "Psat")$S / 1.9872
- else S1 <- water("S", T = T, P = P)$S / 1.9872
- S1
+.Cp1_g <- function(T) {
+ # Heat capacity of the ideal gas (dimensionless)
+ .Fortran(C_ideal2, T, 0, 0, 0, 0, 0, 0, 0)[[8]]
}
Modified: pkg/CHNOSZ/inst/NEWS.Rd
===================================================================
--- pkg/CHNOSZ/inst/NEWS.Rd 2022-02-06 11:33:15 UTC (rev 696)
+++ pkg/CHNOSZ/inst/NEWS.Rd 2022-02-06 14:27:27 UTC (rev 697)
@@ -10,19 +10,19 @@
\newcommand{\s}{\ifelse{latex}{\eqn{_{#1}}}{\ifelse{html}{\out{<sub>#1</sub>}}{#1}}}
\newcommand{\S}{\ifelse{latex}{\eqn{^{#1}}}{\ifelse{html}{\out{<sup>#1</sup>}}{^#1}}}
-\section{Changes in CHNOSZ version 1.4.1-22 (2022-02-06)}{
+\section{Changes in CHNOSZ version 1.4.1-23 (2022-02-06)}{
\subsection{OBIGT DATABASE}{
\itemize{
- \item Move Ar, Xe, CH\s{4}, and CO\s{2} from \samp{organic_aq.csv} to
- \samp{inorganic_aq.csv}.
+ \item Move Ar, Xe, CH\s{4}, and CO\s{2} from \file{organic_aq.csv} to
+ \file{inorganic_aq.csv}.
- \item Remove \samp{OldAA.csv} (superseded thermodynamic parameters for
+ \item Remove \file{OldAA.csv} (superseded thermodynamic parameters for
amino acids). This file is now available in the
\href{https://github.com/jedick/JMDplots}{JMDplots} package.
- \item Fix formula of CaCl\s{2} in \samp{DEW.csv}. Thanks to Grayson Boyer.
+ \item Fix formula of CaCl\s{2} in \file{DEW.csv}. Thanks to Grayson Boyer.
\item Add H\s{2}WO\s{4}(aq) and Cp coefficients of scheelite (CaWO\s{4})
from \href{https://doi.org/10.1016/j.chemgeo.2021.120488}{Liu et al.
@@ -31,7 +31,7 @@
\item Add aqueous species in the system Ca-Mg-Na-K-Al-Si-O-H-C-Cl from
\href{https://doi.org/10.1016/j.gca.2016.04.026}{Miron et al. (2016)} and
\href{https://doi.org/10.2475/07.2017.01}{Miron et al. (2017)} to
- \samp{GEMSFIT.csv}.
+ \file{GEMSFIT.csv}.
}
}
@@ -42,14 +42,15 @@
\item For minerals with thermodynamic parameters in the equations of
\href{https://doi.org/10.1093/petrology/29.2.445}{Berman (1988)},
\code{info()} with a numeric argument now lists values of \code{G},
- \code{H}, \code{S}, and \code{V} at 25 \degC and 1 bar (previously shown
- as \code{NA}).
+ \code{H}, \code{S}, and \code{V} at 25 \degC and 1 bar (these were
+ previously shown as \code{NA}).
\item The \code{AkDi()} function, which provides an implementation of the
\href{https://doi.org/10.1016/S0016-7037(02)01141-9}{Akinfiev and Diamond
(2003)} equation of state for aqueous nonelectrolytes, now calculates
- entropy in addition to the previously available Gibbs energy. Thanks to
- Evgeniy Bastrakov for suggesting both of these changes.
+ entropy, heat capacity, and volume in addition to the previously
+ available Gibbs energy. Thanks to Evgeniy Bastrakov for suggesting both
+ of these changes.
}
}
@@ -59,12 +60,13 @@
\item Tests are now run using the \CRANpkg{tinytest} package.
- \item \code{maxdiff()} and \code{expect_maxdiff()} have been removed.
+ \item Remove \code{maxdiff()} and \code{expect_maxdiff()}.
- \item \strong{multi-metal.Rmd} now has a link to the associated paper
+ \item In \file{vignettes/multi-metal.Rmd} (\dQuote{Diagrams with multiple
+ metals}), add a link to the associated paper
(\href{https://doi.org/10.1016/j.acags.2021.100059}{Dick, 2021}).
- \item All function, variables, and data file names now use capitalized
+ \item Names of functions, variables, and data files now use capitalized
\samp{Berman}, not \samp{berman}.
}
Modified: pkg/CHNOSZ/inst/tinytest/test-AkDi.R
===================================================================
--- pkg/CHNOSZ/inst/tinytest/test-AkDi.R 2022-02-06 11:33:15 UTC (rev 696)
+++ pkg/CHNOSZ/inst/tinytest/test-AkDi.R 2022-02-06 14:27:27 UTC (rev 697)
@@ -19,15 +19,22 @@
# Tests now use values calculated with the AD_full program provided by N. Akinfiev and E. Bastrakov 20210206
info <- "AkDi produces correct values for CO2 along saturation curve"
P <- "Psat"
-T <- seq(50, 350, 100)
+T <- c(50, 150, 250, 350)
# J mol-1
G_ref <- c(-389128.8, -405167.1, -425423.6, -450572.9)
# J K-1 mol-1
S_ref <- c(135.082, 183.295, 226.32, 366.626)
+Cp_ref <- c(189.791, 178.452, 273.662, 7231.723)
+V_ref <- c(32.57, 38.125, 58.269, 298.659)
# Calculate values using AkDi model in CHNOSZ
sout <- subcrt("CO2", T = T, P = P)$out[[1]]
expect_equal(sout$G, G_ref, tolerance = 13, scale = 1, info = info)
expect_equal(sout$S, S_ref, tolerance = 0.8, scale = 1, info = info)
+expect_equal(sout$Cp[1:3], Cp_ref[1:3], tolerance = 14, scale = 1, info = info)
+expect_equal(sout$V[1:3], V_ref[1:3], tolerance = 0.11, scale = 1, info = info)
+# Cp and V get much larger, and so do the differences, near the critical point
+expect_equal(sout$Cp[4], Cp_ref[4], tolerance = 1950, scale = 1, info = info)
+expect_equal(sout$V[4], V_ref[4], tolerance = 24, scale = 1, info = info)
info <- "AkDi produces correct values for CO2 at 1000 bar"
P <- 1000
@@ -34,11 +41,20 @@
T <- c(300, 400, 500, 600)
G_ref <- c(-431672.1, -455065.1, -480692, -507005.3)
S_ref <- c(221.255, 246.139, 263.759, 259.812)
+Cp_ref <- c(150.613, 157.014, 54.489, -76.384)
+V_ref <- c(45.831, 67.408, 107.656, 129.264)
# Calculate values using AkDi model in CHNOSZ
sout <- subcrt("CO2", T = T, P = P)$out[[1]]
expect_equal(sout$G, G_ref, tolerance = 11, scale = 1, info = info)
expect_equal(sout$S, S_ref, tolerance = 0.09, scale = 1, info = info)
+expect_equal(sout$Cp, Cp_ref, tolerance = 10, scale = 1, info = info)
+expect_equal(sout$V, V_ref, tolerance = 0.4, scale = 1, info = info)
+info <- "AkDi gives consistent values of G, H, and S"
+T_in_Kelvin <- convert(T, "K")
+S_of_elements_in_Joules <- convert(entropy("CO2"), "J")
+expect_equal(sout$H - T_in_Kelvin * sout$S + 298.15 * S_of_elements_in_Joules, sout$G, tolerance = 0.24, scale = 1)
+
# 20220206
info <- "Fugacity, density, and density derivatives of H2O are close to values in Akinfiev and Diamond (2003)"
# This is the value of fugacity given in the paper
@@ -51,7 +67,7 @@
# g / cm3 / K
expect_equal(CHNOSZ:::.drho1_dT(298.15, 1), -0.0002571, tolerance = 0.000002, scale = 1, info = info)
# g / cm3 / bar
-expect_equal(CHNOSZ:::.drho1_dP(298.15, 1), 0.00004511, tolerance = 0.00000003, scale = 1, info = info)
+expect_equal(CHNOSZ:::.drho1_dP(298.15, 1), 0.00004511, tolerance = 0.0000001, scale = 1, info = info)
# g / cm3 / K^2
expect_equal(CHNOSZ:::.d2rho1_dT2(298.15, 1), -0.000009503, tolerance = 0.0000014, scale = 1, info = info)
Modified: pkg/CHNOSZ/man/eos.Rd
===================================================================
--- pkg/CHNOSZ/man/eos.Rd 2022-02-06 11:33:15 UTC (rev 696)
+++ pkg/CHNOSZ/man/eos.Rd 2022-02-06 14:27:27 UTC (rev 697)
@@ -50,7 +50,6 @@
\code{AkDi} implements the Akinfiev-Diamond model for aqueous species (Akinfiev and Diamond, 2003).
This model is activated by setting \code{abbrv = "AkDi"} in \code{thermo()$OBIGT} for a given aqueous species.
The database must also include the corresponding gasesous species (with the same name or chemical formula).
-Currently, only the standard chemical potential (Gibbs energy) and entropy are calculated.
}
Modified: pkg/CHNOSZ/vignettes/multi-metal.Rmd
===================================================================
--- pkg/CHNOSZ/vignettes/multi-metal.Rmd 2022-02-06 11:33:15 UTC (rev 696)
+++ pkg/CHNOSZ/vignettes/multi-metal.Rmd 2022-02-06 14:27:27 UTC (rev 697)
@@ -85,7 +85,7 @@
This vignette was compiled on `r Sys.Date()` with CHNOSZ version `r sessionInfo()$otherPkgs$CHNOSZ$Version`.
-An accompanying paper to this vignette has been published in *Applied Computing and Geosciences* ([Dick, 2021](https://doi.org/10.1016/j.acags.2021.100059 "Diagrams with multiple metals in CHNOSZ")).
+This vignette is associated with a paper that has been published in *Applied Computing and Geosciences* ([Dick, 2021](https://doi.org/10.1016/j.acags.2021.100059 "Diagrams with multiple metals in CHNOSZ")).
Please consider citing that paper if you use the functions or diagrams described here.
The plots in this vignette were made using the following resolution settings, which can be changed if desired (low resolutions are used to save time in CRAN checks):
More information about the CHNOSZ-commits
mailing list