[CHNOSZ-commits] r696 - in pkg/CHNOSZ: . R inst inst/extdata/OBIGT inst/tinytest man src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sun Feb 6 12:33:15 CET 2022
Author: jedick
Date: 2022-02-06 12:33:15 +0100 (Sun, 06 Feb 2022)
New Revision: 696
Added:
pkg/CHNOSZ/inst/tinytest/test-AkDi.R
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/NAMESPACE
pkg/CHNOSZ/R/AkDi.R
pkg/CHNOSZ/R/basis.R
pkg/CHNOSZ/R/findit.R
pkg/CHNOSZ/inst/NEWS.Rd
pkg/CHNOSZ/inst/TODO
pkg/CHNOSZ/inst/extdata/OBIGT/AkDi.csv
pkg/CHNOSZ/inst/tinytest/test-eos.R
pkg/CHNOSZ/man/eos.Rd
pkg/CHNOSZ/src/H2O92D.f
pkg/CHNOSZ/src/init.c
Log:
AkDi() now calculates entropy; more tests added to test-AkDi.R
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2022-02-04 04:10:05 UTC (rev 695)
+++ pkg/CHNOSZ/DESCRIPTION 2022-02-06 11:33:15 UTC (rev 696)
@@ -1,6 +1,6 @@
-Date: 2022-02-03
+Date: 2022-02-06
Package: CHNOSZ
-Version: 1.4.1-21
+Version: 1.4.1-22
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-02-04 04:10:05 UTC (rev 695)
+++ pkg/CHNOSZ/NAMESPACE 2022-02-06 11:33:15 UTC (rev 696)
@@ -74,7 +74,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")
+ "uniroot", "median", "predict", "smooth.spline")
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/AkDi.R
===================================================================
--- pkg/CHNOSZ/R/AkDi.R 2022-02-04 04:10:05 UTC (rev 695)
+++ pkg/CHNOSZ/R/AkDi.R 2022-02-06 11:33:15 UTC (rev 696)
@@ -1,46 +1,52 @@
# CHNOSZ/AkDi.R
# Akinfiev-Diamond model for aqueous species
-# 20190219 first version
+# 20190219 First version
+# 20220206 Calculate S, Cp, and V
AkDi <- function(property = NULL, parameters = NULL, T = 298.15, P = 1, isPsat = TRUE) {
# Some constants (from Akinfiev and Diamond, 2004 doi:10.1016/j.fluid.2004.06.010)
MW <- 18.0153 # g mol-1
- NW <- 1000/MW # mol kg-1
+ NW <- 1000 / MW # mol kg-1
R <- 8.31441 # J K-1 mol-1
+ # R expressed in volume units
+ RV <- 10 * R # cm3 bar K-1 mol-1
- # A list for the output
+ # Calculate properties of H2O
+ # Moved these calculations to 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)
+ S1_g <- sapply(T, .S1_g)
+ S1 <- mapply(.S1, T, P, MoreArgs = list(isPsat = isPsat))
+
+#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)) {
PAR <- parameters[i, ]
+ # Send a message
+ message("AkDi: Akinfiev-Diamond model for ", PAR$name, " gas to aq")
# Start with an NA-filled data frame
myprops <- as.data.frame(matrix(NA, ncol = length(property), nrow = length(T)))
colnames(myprops) <- property
+
# Loop over properties
for(j in seq_along(property)) {
- # Just calculate G for now
+
if(property[[j]] == "G") {
- # Send a message
- message("AkDi: Akinfiev-Diamond model for ", PAR$name, " gas to aq")
# 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)
- # Get H2O fugacity (bar)
- GH2O_P <- water("G", T = T, P = P)$G
- GH2O_1 <- water("G", T = T, P = 1)$G
- f1 <- exp ( (GH2O_P - GH2O_1) / (1.9872 * T) )
- # For Psat, calculate the real liquid-vapor curve (not 1 bar below 100 degC)
- if(isPsat) {
- P <- water("Psat", T = T, P = "Psat", P1 = FALSE)$Psat
- f1[P < 1] <- P[P < 1]
- }
- # Density (g cm-3)
- rho1 <- water("rho", T = T, P = P)$rho / 1000
# Calculate G_hyd (J mol-1)
- G_hyd <- R*T * ( -log(NW) + (1 - PAR$xi) * log(f1) + PAR$xi * log(10 * R * T * rho1 / MW) + rho1 * (PAR$a + PAR$b * (1000/T)^0.5) )
+ 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)
G <- G_gas + G_hyd
# Convert J to cal
@@ -48,11 +54,92 @@
# Insert into data frame of properties
myprops$G <- G
}
+
if(property[[j]] == "S") {
- # https://stackoverflow.com/questions/11081069/calculate-the-derivative-of-a-data-function-in-r
+ # 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)
+ # Calculate S_hyd
+ S_hyd <- R * (
+ (1 - PAR$xi) * (S1 - S1_g)
+ + log(NW)
+ - (PAR$xi + PAR$xi * log(RV * T / MW) + PAR$xi * log(rho1) + PAR$xi * T / rho1 * drho1_dT)
+ - (
+ PAR$a * (rho1 + T * drho1_dT)
+ + PAR$b * (0.5 * 10^1.5 * T^-0.5 * rho1 + 10^1.5 * T^0.5 * drho1_dT)
+ )
+ )
+ S <- S_gas + S_hyd
+ S <- convert(S, "cal", T = T)
+ myprops$S <- S
}
+
}
out[[i]] <- myprops
}
out
}
+
+### UNEXPORTED FUNCTIONS ###
+
+.f1 <- function(T, P, isPsat) {
+ # Get H2O fugacity (bar)
+ GH2O_P <- water("G", T = T, P = P)$G
+ GH2O_1 <- water("G", T = T, P = 1)$G
+ f1 <- exp ( (GH2O_P - GH2O_1) / (1.9872 * T) )
+ # For Psat, calculate the real liquid-vapor curve (not 1 bar below 100 degC)
+ if(isPsat) {
+ P <- water("Psat", T = T, P = "Psat", P1 = FALSE)$Psat
+ f1[P < 1] <- P[P < 1]
+ }
+ f1
+}
+
+.rho1 <- function(T, P) {
+ # Density of H2O (g cm-3)
+ water("rho", T = T, P = P)$rho / 1000
+}
+
+.drho1_dT <- function(T, P) {
+ # Partial derivative of density with respect to temperature at constant pressure 20220206
+ dT <- 0.1
+ T1 <- T - dT
+ T2 <- T + dT
+ # Add 1 bar to P so the derivative doesn't blow up when P = Psat at T > 100 degC
+ # TODO: Is there a better way?
+ rho1 <- .rho1(c(T1, T2), P + 1)
+ diff(rho1) / (T2 - T1)
+}
+
+.drho1_dP <- function(T, P) {
+ # Partial derivative of density with respect to pressure at constant temperature 20220206
+ dP <- 0.1
+ P1 <- P - dP
+ P2 <- P + dP
+ rho1 <- .rho1(T, c(P1, P2))
+ diff(rho1) / (P2 - P1)
+}
+
+.d2rho1_dT2 <- function(T, P) {
+ # Second partial derivative of density with respect to temperature at constant pressure 20220206
+ dT <- 0.1
+ # Calculate density at five temperature values
+ Tval <- seq(T - 2 * dT, T + 2 * dT, dT)
+ rho1 <- .rho1(Tval, P)
+ # 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)
+ predict(spl, deriv = 2)$y[3]
+}
+
+.S1_g <- function(T) {
+ # Entropy of the ideal gas (dimensionless)
+ .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
+}
Modified: pkg/CHNOSZ/R/basis.R
===================================================================
--- pkg/CHNOSZ/R/basis.R 2022-02-04 04:10:05 UTC (rev 695)
+++ pkg/CHNOSZ/R/basis.R 2022-02-06 11:33:15 UTC (rev 696)
@@ -54,7 +54,7 @@
if(!is.null(logact)) species[species=="pe"] <- "e-"
}
if("Eh" %in% species) {
- # 20090209 should be careful with this conversion as it's only for 25 deg C
+ # 20090209 should be careful with this conversion as it's only for 25 degC
# to be sure, just don't call species("Eh")
if(!is.null(logact)) logact[species=="Eh"] <- -convert(logact[species=="Eh"],"pe")
species[species=="Eh"] <- "e-"
Modified: pkg/CHNOSZ/R/findit.R
===================================================================
--- pkg/CHNOSZ/R/findit.R 2022-02-04 04:10:05 UTC (rev 695)
+++ pkg/CHNOSZ/R/findit.R 2022-02-06 11:33:15 UTC (rev 696)
@@ -11,7 +11,7 @@
# we iteratively move toward a higher/lower value of the objective
# within these limits
- # fun stuff: when running on either side of 100 deg C,
+ # fun stuff: when running on either side of 100 degC,
# set P=1 to force an extreme of some functions near that temperature
# the number of dimensions
Modified: pkg/CHNOSZ/inst/NEWS.Rd
===================================================================
--- pkg/CHNOSZ/inst/NEWS.Rd 2022-02-04 04:10:05 UTC (rev 695)
+++ pkg/CHNOSZ/inst/NEWS.Rd 2022-02-06 11:33:15 UTC (rev 696)
@@ -10,7 +10,7 @@
\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-16 (2022-02-03)}{
+\section{Changes in CHNOSZ version 1.4.1-22 (2022-02-06)}{
\subsection{OBIGT DATABASE}{
\itemize{
@@ -28,17 +28,28 @@
from \href{https://doi.org/10.1016/j.chemgeo.2021.120488}{Liu et al.
(2021)}.
- \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}.
+ \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}.
}
}
- \subsection{TESTING}{
+ \subsection{THERMODYNAMIC CALCULATIONS}{
\itemize{
- \item Tests are now run using the \CRANpkg{tinytest} package.
+ \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}).
- \item \code{maxdiff()} and \code{expect_maxdiff()} have been removed.
+ \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.
}
}
@@ -46,16 +57,16 @@
\subsection{OTHER CHANGES}{
\itemize{
+ \item Tests are now run using the \CRANpkg{tinytest} package.
+
+ \item \code{maxdiff()} and \code{expect_maxdiff()} have been removed.
+
\item \strong{multi-metal.Rmd} now has a link to the associated paper
(\href{https://doi.org/10.1016/j.acags.2021.100059}{Dick, 2021}).
- \item \code{info()} with a numeric argument now includes values of
- \code{G}, \code{H}, \code{S}, and \code{V} at 25 \degC and 1 bar for
- minerals with thermodynamic parameters in the Berman equations.
-
\item All function, variables, and data file names now use capitalized
\samp{Berman}, not \samp{berman}.
-
+
}
}
Modified: pkg/CHNOSZ/inst/TODO
===================================================================
--- pkg/CHNOSZ/inst/TODO 2022-02-04 04:10:05 UTC (rev 695)
+++ pkg/CHNOSZ/inst/TODO 2022-02-06 11:33:15 UTC (rev 696)
@@ -50,4 +50,4 @@
- Allow numeric values for bases in mosaic() (test-mix.R)
- Add a 'model' column to thermo()$OBIGT to indicate the thermodynamic data
- model (e.g. HKF, CGL, Berman, AkDi, or DEW)
+ model (e.g. H2O, HKF, CGL, Berman, AkDi, or DEW)
Modified: pkg/CHNOSZ/inst/extdata/OBIGT/AkDi.csv
===================================================================
--- pkg/CHNOSZ/inst/extdata/OBIGT/AkDi.csv 2022-02-04 04:10:05 UTC (rev 695)
+++ pkg/CHNOSZ/inst/extdata/OBIGT/AkDi.csv 2022-02-06 11:33:15 UTC (rev 696)
@@ -20,6 +20,6 @@
B(OH)3,AkDi,B(OH)3,aq,AP14,NA,2019-02-22,cal,NA,NA,NA,NA,NA,-4.2561,4.0194,-1.057,NA,NA,NA,NA,NA
Si(OH)4,AkDi,Si(OH)4,aq,AP14,NA,2019-02-22,cal,NA,NA,NA,NA,NA,0.9285,-0.9409,-1.8933,NA,NA,NA,NA,NA
As(OH)3,AkDi,As(OH)3,aq,AP14,NA,2019-02-22,cal,NA,NA,NA,NA,NA,-9.903,7.6818,-1.23,NA,NA,NA,NA,NA
-B(OH)3,AkDi,B(OH)3,gas,AP14,NA,2019-02-22,cal,-222523.9,NA,66.589,NA,0,46.7734,-0.097,0.5932,-522.4665,0,0,0,1000
-Si(OH)4,AkDi,Si(OH)4,gas,AP14,NA,2019-02-22,cal,-296285.9,NA,82.784,NA,0,22.9685,37.6434,-4.4407,0,-2.7782,7.923E-09,3,1000
-As(OH)3,AkDi,As(OH)3,gas,AP14,NA,2019-02-22,cal,-143112.8,NA,73.599,NA,0,22.6338,-16.3623,-3.7213,0,-1.0194,2.7414E-09,3,1000
+B(OH)3,NA,B(OH)3,gas,AP14,NA,2019-02-22,cal,-222523.9,NA,66.589,NA,0,46.7734,-0.097,0.5932,-522.4665,0,0,0,1000
+Si(OH)4,NA,Si(OH)4,gas,AP14,NA,2019-02-22,cal,-296285.9,NA,82.784,NA,0,22.9685,37.6434,-4.4407,0,-2.7782,7.923E-09,3,1000
+As(OH)3,NA,As(OH)3,gas,AP14,NA,2019-02-22,cal,-143112.8,NA,73.599,NA,0,22.6338,-16.3623,-3.7213,0,-1.0194,2.7414E-09,3,1000
Added: pkg/CHNOSZ/inst/tinytest/test-AkDi.R
===================================================================
--- pkg/CHNOSZ/inst/tinytest/test-AkDi.R (rev 0)
+++ pkg/CHNOSZ/inst/tinytest/test-AkDi.R 2022-02-06 11:33:15 UTC (rev 696)
@@ -0,0 +1,73 @@
+# Load default settings for CHNOSZ
+reset()
+# Set subcrt() to output values in Joules
+E.units("J")
+
+# 20220206
+info <- "Database is set up correctly"
+# Activate the AkDi model for aqueous nonelectrolytes
+iAkDi <- add.OBIGT("AkDi")
+# Get the indices of the modified aqueous species
+iaq <- iAkDi[info(iAkDi)$state == "aq"]
+# Each aqueous species should be associated with the AkDi model
+expect_equal(unique(info(iaq)$abbrv), "AkDi", info = info)
+# Each aqueous species should have a gaseous counterpart
+igas <- info(info(iAkDi)$name, "gas")
+expect_true(!any(is.na(igas)), info = info)
+
+# First version of tests used a circular reference (values previously calculated in CHNOSZ) 20190220
+# 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)
+# 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)
+# 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)
+
+info <- "AkDi produces correct values for CO2 at 1000 bar"
+P <- 1000
+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)
+# 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)
+
+# 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
+expect_equal(log(CHNOSZ:::.f1(298.15, 1, TRUE)), -3.4773, tolerance = 0.026, scale = 1, info = info)
+# This is the value of fugacity output by AD_full
+expect_equal(log(CHNOSZ:::.f1(298.15, 1, TRUE)), -3.4524, tolerance = 0.0007, scale = 1, info = info)
+# The following values are from the paper
+# g / cm3
+expect_equal(CHNOSZ:::.rho1(298.15, 1), 0.9970, tolerance = 0.0001, scale = 1, info = info)
+# 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)
+# g / cm3 / K^2
+expect_equal(CHNOSZ:::.d2rho1_dT2(298.15, 1), -0.000009503, tolerance = 0.0000014, scale = 1, info = info)
+
+# 20190220 Compare Gibbs energies at 25 degrees calculated with AkDi model to default OBIGT database
+info <- "Gibbs energies at 25 degree C are comparable between AkDi and default OBIGT database"
+# Remove some hydroxides that aren't in the default database
+iaq <- iaq[!info(iaq)$name %in% c("Si(OH)4", "As(OH)3")]
+# This would produce an error if any calculations failed
+# (e.g. because gases corresponding to any aqueous species were unavailable)
+sAkDi <- subcrt(iaq, T = 25)
+GAkDi <- do.call(rbind, sAkDi$out)$G
+# Now calculate the parameters using default OBIGT database
+OBIGT()
+sOBIGT <- subcrt(iaq, T = 25)
+GOBIGT <- do.call(rbind, sOBIGT$out)$G
+# Mean absolute differences are less than 280 J / mol
+expect_equal(GAkDi, GOBIGT, tolerance = 280, scale = 1, info = info)
+# The largest differences are for HCl, ethane, and B(OH)3
+expect_equal(sort(info(iaq[abs(GAkDi - GOBIGT) > 900])$name), sort(c("HCl", "ethane", "B(OH)3")))
Modified: pkg/CHNOSZ/inst/tinytest/test-eos.R
===================================================================
--- pkg/CHNOSZ/inst/tinytest/test-eos.R 2022-02-04 04:10:05 UTC (rev 695)
+++ pkg/CHNOSZ/inst/tinytest/test-eos.R 2022-02-06 11:33:15 UTC (rev 696)
@@ -77,42 +77,6 @@
expect_equal(gfun.1000$dgdP * 1e6, dgdP.1000.ref, tolerance = 1e-1, info = info)
expect_equal(gfun.4000$dgdP * 1e6, dgdP.4000.ref, tolerance = 1e-3, info = info)
-info <- "AkDi produces expected results"
-# 20190220
-# modify aqueous CO2 to use the AkDi model
-iCO2 <- mod.OBIGT("CO2", abbrv = "AkDi", a = -8.8321, b = 11.2684, c = -0.0850)
-# do the properties we calculate match previously calculated values?
-P <- "Psat"
-T <- seq(50, 350, 100)
-# J mol-1
-G_ref <- c(-389122.3, -405138.4, -425410.7, -450573.2)
-G_calc <- subcrt(iCO2, T = T, P = P)$out[[1]]$G
-# convert to J mol-1
-G_calc <- convert(G_calc, "J", T = convert(T, "K"))
-expect_equal(round(G_calc, 1), G_ref, info = info)
-
-P <- 500
-T <- seq(200, 1000, 200)
-G_ref <- c(-412767.4, -459654.1, -515231.6, -565736.3, -617927.9)
-G_calc <- subcrt(iCO2, T = T, P = P)$out[[1]]$G
-G_calc <- convert(G_calc, "J", T = convert(T, "K"))
-expect_equal(round(G_calc, 1), G_ref, info = info)
-
-# compare Gibbs energies at 25 degrees calculated with AkDi model to database values
-iAkDi <- add.OBIGT("AkDi")
-# remove hydroxides because they aren't in the default database (except B(OH)3(aq))
-iAkDi <- iAkDi[-grep("OH", info(iAkDi, check.it = FALSE)$name)]
-# this would produce an error if any calculations failed
-# (e.g. because gases corresponding to any aqueous species were unavailable)
-sAkDi <- subcrt(iAkDi, T = 25)
-GAkDi <- do.call(rbind, sAkDi$out)$G
-# now get the parameters from default OBIGT
-reset()
-GOBIGT <- info(iAkDi, check.it = FALSE)$G
-# the differences are not that big, except for HCl(aq)
-maxdiff <- function(x, y) max(abs(y - x))
-expect_true(maxdiff(GAkDi, GOBIGT) < 300, info = info)
-
# reference
# Shock, E. L., Oelkers, E. H., Johnson, J. W., Sverjensky, D. A. and Helgeson, H. C. (1992)
Modified: pkg/CHNOSZ/man/eos.Rd
===================================================================
--- pkg/CHNOSZ/man/eos.Rd 2022-02-04 04:10:05 UTC (rev 695)
+++ pkg/CHNOSZ/man/eos.Rd 2022-02-06 11:33:15 UTC (rev 696)
@@ -50,7 +50,7 @@
\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) is calculated.
+Currently, only the standard chemical potential (Gibbs energy) and entropy are calculated.
}
Modified: pkg/CHNOSZ/src/H2O92D.f
===================================================================
--- pkg/CHNOSZ/src/H2O92D.f 2022-02-04 04:10:05 UTC (rev 695)
+++ pkg/CHNOSZ/src/H2O92D.f 2022-02-06 11:33:15 UTC (rev 696)
@@ -1106,6 +1106,50 @@
RETURN
END
+************************************************************************
+
+*** Modified version of ideal with no COMMON block 20220206 jmd
+*** ideal2 - Computes thermodynamic properties for H2O in the
+* ideal gas state using equations given by Woolley (1979).
+
+ SUBROUTINE ideal2(t, ai, gi, si, ui, hi, cvi, cpi)
+
+ IMPLICIT DOUBLE PRECISION (a-h,o-z)
+
+ DOUBLE PRECISION c(18), ai, gi, si, ui, hi, cvi, cpi
+
+ SAVE
+
+ DATA c / .19730271018d2, .209662681977d2, -.483429455355d0,
+ 1 .605743189245d1, .2256023885d2, -.987532442d1,
+ 2 -.43135538513d1, .458155781d0, -.47754901883d-1,
+ 3 .41238460633d-2, -.27929052852d-3, .14481695261d-4,
+ 4 -.56473658748d-6, .16200446d-7, -.3303822796d-9,
+ 5 .451916067368d-11, -.370734122708d-13,
+ 6 .137546068238d-15/
+
+
+ tt = t / 1.0d2
+ tl = DLOG(tt)
+ gi = -(c(1)/tt + c(2)) * tl
+ hi = (c(2) + c(1)*(1.0d0 - tl)/tt)
+ cpi = c(2) - c(1)/tt
+
+ DO i=3,18
+ emult = power(tt,DBLE(i-6))
+ gi = gi - c(i) * emult
+ hi = hi + c(i) * (i-6) * emult
+ cpi = cpi + c(i) * (i-6) * (i-5) * emult
+ END DO
+
+ ai = gi - 1.0d0
+ ui = hi - 1.0d0
+ cvi = cpi - 1.0d0
+ si = ui - ai
+
+ RETURN
+ END
+
******************************************************************************
*** dalHGK - Computes/returns (d(alpha)/dt)p(d,t,alpha)
Modified: pkg/CHNOSZ/src/init.c
===================================================================
--- pkg/CHNOSZ/src/init.c 2022-02-04 04:10:05 UTC (rev 695)
+++ pkg/CHNOSZ/src/init.c 2022-02-06 11:33:15 UTC (rev 696)
@@ -5,8 +5,11 @@
/* .Fortran calls */
extern void F77_NAME(h2o92)(void *, void *, void *, void *);
+extern void F77_NAME(ideal2)(void *, void *, void *, void *, void *, void *, void *, void *);
+
static const R_FortranMethodDef FortranEntries[] = {
{"h2o92", (DL_FUNC) &F77_NAME(h2o92), 4},
+ {"ideal2", (DL_FUNC) &F77_NAME(ideal2), 8},
{NULL, NULL, 0}
};
More information about the CHNOSZ-commits
mailing list