[CHNOSZ-commits] r698 - in pkg/CHNOSZ: . R demo inst inst/extdata/OBIGT inst/tinytest man vignettes
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Feb 7 04:46:09 CET 2022
Author: jedick
Date: 2022-02-07 04:46:08 +0100 (Mon, 07 Feb 2022)
New Revision: 698
Added:
pkg/CHNOSZ/R/AD.R
pkg/CHNOSZ/demo/AD.R
pkg/CHNOSZ/inst/extdata/OBIGT/AD.csv
pkg/CHNOSZ/inst/tinytest/test-AD.R
Removed:
pkg/CHNOSZ/R/AkDi.R
pkg/CHNOSZ/demo/AkDi.R
pkg/CHNOSZ/inst/extdata/OBIGT/AkDi.csv
pkg/CHNOSZ/inst/tinytest/test-AkDi.R
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/NAMESPACE
pkg/CHNOSZ/R/examples.R
pkg/CHNOSZ/R/info.R
pkg/CHNOSZ/R/subcrt.R
pkg/CHNOSZ/R/util.data.R
pkg/CHNOSZ/demo/00Index
pkg/CHNOSZ/demo/aluminum.R
pkg/CHNOSZ/demo/sources.R
pkg/CHNOSZ/inst/NEWS.Rd
pkg/CHNOSZ/inst/TODO
pkg/CHNOSZ/inst/tinytest/test-info.R
pkg/CHNOSZ/man/eos.Rd
pkg/CHNOSZ/man/examples.Rd
pkg/CHNOSZ/man/extdata.Rd
pkg/CHNOSZ/man/thermo.Rd
pkg/CHNOSZ/vignettes/OBIGT.Rmd
pkg/CHNOSZ/vignettes/anintro.Rmd
Log:
Rename AkDi to AD
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2022-02-06 14:27:27 UTC (rev 697)
+++ pkg/CHNOSZ/DESCRIPTION 2022-02-07 03:46:08 UTC (rev 698)
@@ -1,6 +1,6 @@
-Date: 2022-02-06
+Date: 2022-02-07
Package: CHNOSZ
-Version: 1.4.1-23
+Version: 1.4.1-24
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-06 14:27:27 UTC (rev 697)
+++ pkg/CHNOSZ/NAMESPACE 2022-02-07 03:46:08 UTC (rev 698)
@@ -55,7 +55,7 @@
# added 20171121 or later
"dumpdata", "thermo.axis", "solubility", "NaCl",
# added 20190213 or later
- "CHNOSZ", "thermo", "reset", "OBIGT", "retrieve", "AkDi", "moles",
+ "CHNOSZ", "thermo", "reset", "OBIGT", "retrieve", "AD", "moles",
"lNaCl", "lS", "lT", "lP", "lTP", "lex",
# added 20200716 or later
"mash", "mix", "rebalance"
Copied: pkg/CHNOSZ/R/AD.R (from rev 697, pkg/CHNOSZ/R/AkDi.R)
===================================================================
--- pkg/CHNOSZ/R/AD.R (rev 0)
+++ pkg/CHNOSZ/R/AD.R 2022-02-07 03:46:08 UTC (rev 698)
@@ -0,0 +1,182 @@
+# CHNOSZ/AD.R
+# Akinfiev-Diamond model for aqueous species
+# 20190219 First version
+# 20220206 Calculate S, Cp, and V
+
+AD <- 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
+ R <- 8.31441 # J K-1 mol-1
+ # R expressed in volume units
+ RV <- 10 * R # cm3 bar K-1 mol-1
+
+ # Calculate H2O fugacity and derivatives of density
+ # These calculations are done through (unexported) functions
+ # to be able to test their output in test-AD.R 20220206
+ f1 <- .f1(T, P, isPsat)
+ 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)
+ Cp1_g <- sapply(T, .Cp1_g)
+
+ # 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("AD: 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)) {
+
+ if(property[[j]] == "G") {
+ # Get gas properties (J mol-1)
+ 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)
+ G <- G_gas + G_hyd
+ # Convert J to cal
+ G <- convert(G, "cal", T = T)
+ # Insert into data frame of properties
+ myprops$G <- G
+ }
+
+ if(property[[j]] == "S") {
+ # Get S_gas
+ S_gas <- convert(gasprops$S, "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
+ }
+
+ 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
+}
+
+### 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
+ # 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)
+}
+
+.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)
+ # 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)
+ 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]]
+}
+
+.Cp1_g <- function(T) {
+ # Heat capacity of the ideal gas (dimensionless)
+ .Fortran(C_ideal2, T, 0, 0, 0, 0, 0, 0, 0)[[8]]
+}
Deleted: pkg/CHNOSZ/R/AkDi.R
===================================================================
--- pkg/CHNOSZ/R/AkDi.R 2022-02-06 14:27:27 UTC (rev 697)
+++ pkg/CHNOSZ/R/AkDi.R 2022-02-07 03:46:08 UTC (rev 698)
@@ -1,182 +0,0 @@
-# CHNOSZ/AkDi.R
-# Akinfiev-Diamond model for aqueous species
-# 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
- R <- 8.31441 # J K-1 mol-1
- # R expressed in volume units
- RV <- 10 * R # cm3 bar K-1 mol-1
-
- # 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)
- 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)
- Cp1_g <- sapply(T, .Cp1_g)
-
- # 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
- myprops <- as.data.frame(matrix(NA, ncol = length(property), nrow = length(T)))
- colnames(myprops) <- property
-
- # Loop over properties
- for(j in seq_along(property)) {
-
- if(property[[j]] == "G") {
- # Get gas properties (J mol-1)
- 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)
- G <- G_gas + G_hyd
- # Convert J to cal
- G <- convert(G, "cal", T = T)
- # Insert into data frame of properties
- myprops$G <- G
- }
-
- if(property[[j]] == "S") {
- # Get S_gas
- S_gas <- convert(gasprops$S, "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
- }
-
- 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
-}
-
-### 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
- # 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)
-}
-
-.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)
- # 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)
- 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]]
-}
-
-.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/R/examples.R
===================================================================
--- pkg/CHNOSZ/R/examples.R 2022-02-06 14:27:27 UTC (rev 697)
+++ pkg/CHNOSZ/R/examples.R 2022-02-07 03:46:08 UTC (rev 698)
@@ -34,7 +34,7 @@
"ORP", "findit", "ionize", "buffer", "protbuff", "glycinate",
"mosaic", "copper", "arsenic", "solubility", "gold", "contour", "sphalerite", "zinc",
"Shh", "saturation", "adenine", "DEW", "lambda", "potassium", "TCA", "aluminum",
- "AkDi", "comproportionation", "Pourbaix", "E_coli"), save.png=FALSE) {
+ "AD", "comproportionation", "Pourbaix", "E_coli"), save.png=FALSE) {
# run one or more demos from CHNOSZ with ask=FALSE, and return the value of the last one
for(i in 1:length(which)) {
# say something so the user sees where we are
Modified: pkg/CHNOSZ/R/info.R
===================================================================
--- pkg/CHNOSZ/R/info.R 2022-02-06 14:27:27 UTC (rev 697)
+++ pkg/CHNOSZ/R/info.R 2022-02-07 03:46:08 UTC (rev 698)
@@ -214,7 +214,7 @@
}
# Perform consistency checks for GHS and EOS parameters if check.it = TRUE
- # Don't do it for the AkDi species 20190219
+ # Don't do it for the AD species 20190219
if(check.it & !"xi" %in% colnames(this)) {
# check GHS if they are all present
if(sum(naGHS)==0) calcG <- checkGHS(this)
@@ -269,7 +269,7 @@
}
# if we got here there were no approximate matches
# 20190127 look for the species in optional data files
- for(opt in c("SLOP98", "SUPCRT92", "AkDi")) {
+ for(opt in c("SLOP98", "SUPCRT92", "AD")) {
optdat <- read.csv(system.file(paste0("extdata/OBIGT/", opt, ".csv"), package="CHNOSZ"), as.is=TRUE)
if(species %in% optdat$name) {
message('info.approx: ', species, ' is in an optional database; use add.OBIGT("', opt, '", "', species, '") to load it')
Modified: pkg/CHNOSZ/R/subcrt.R
===================================================================
--- pkg/CHNOSZ/R/subcrt.R 2022-02-06 14:27:27 UTC (rev 697)
+++ pkg/CHNOSZ/R/subcrt.R 2022-02-07 03:46:08 UTC (rev 698)
@@ -9,7 +9,7 @@
#source("util.units.R")
#source("util.data.R")
#source("species.R")
-#source("AkDi.R")
+#source("AD.R")
#source("nonideal.R")
subcrt <- function(species, coeff = 1, state = NULL, property = c("logK", "G", "H", "S", "V", "Cp"),
@@ -106,7 +106,7 @@
P <- rep(tpargs$P,length.out=length(newIS))
}
} else {
- # for AkDi, remember if P = "Psat" 20190219
+ # for AD, remember if P = "Psat" 20190219
if(identical(P, "Psat")) isPsat <- TRUE
# expansion of Psat and equivalence of argument lengths
tpargs <- TP.args(T=T,P=P)
@@ -284,10 +284,10 @@
if(TRUE %in% isaq) {
# 20110808 get species parameters using OBIGT2eos() (faster than using info())
param <- OBIGT2eos(thermo$OBIGT[iphases[isaq],], "aq", fixGHS = TRUE, tocal = TRUE)
- # Aqueous species with abbrv = "AkDi" use the AkDi model 20210407
+ # Aqueous species with abbrv = "AD" use the AD model 20210407
abbrv <- thermo$OBIGT$abbrv[iphases[isaq]]
abbrv[is.na(abbrv)] <- ""
- isAkDi <- abbrv == "AkDi"
+ isAD <- abbrv == "AD"
# always get density
H2O.props <- "rho"
# calculate A_DH and B_DH if we're using the B-dot (Helgeson) equation
@@ -294,14 +294,14 @@
if(any(IS != 0) & thermo$opt$nonideal %in% c("Bdot", "Bdot0", "bgamma", "bgamma0")) H2O.props <- c(H2O.props, "A_DH", "B_DH")
# get other properties for H2O only if it's in the reaction
if(any(isH2O)) H2O.props <- c(H2O.props, eosprop)
- # in case everything is AkDi, run hkf (for water properties) but exclude all species
+ # in case everything is AD, run hkf (for water properties) but exclude all species
hkfpar <- param
- if(all(isAkDi)) hkfpar <- param[0, ]
+ if(all(isAD)) hkfpar <- param[0, ]
hkfstuff <- hkf(eosprop, parameters = hkfpar, T = T, P = P, H2O.props=H2O.props)
p.aq <- hkfstuff$aq
H2O.PT <- hkfstuff$H2O
# set properties to NA for density below 0.35 g/cm3 (a little above the critical isochore, threshold used in SUPCRT92) 20180922
- if(!exceed.rhomin & !all(isAkDi)) {
+ if(!exceed.rhomin & !all(isAD)) {
ilowrho <- H2O.PT$rho < 350
ilowrho[is.na(ilowrho)] <- FALSE
if(any(ilowrho)) {
@@ -311,10 +311,10 @@
}
}
# calculate properties using Akinfiev-Diamond model 20190219
- if(any(isAkDi)) {
+ if(any(isAD)) {
# get the parameters with the right names
- param <- OBIGT2eos(param[isAkDi, ], "aq", tocal = TRUE)
- p.aq[isAkDi] <- AkDi(eosprop, parameters = param, T = T, P = P, isPsat = isPsat)
+ param <- OBIGT2eos(param[isAD, ], "aq", tocal = TRUE)
+ p.aq[isAD] <- AD(eosprop, parameters = param, T = T, P = P, isPsat = isPsat)
}
# calculate activity coefficients if ionic strength is not zero
if(any(IS != 0)) {
Modified: pkg/CHNOSZ/R/util.data.R
===================================================================
--- pkg/CHNOSZ/R/util.data.R 2022-02-06 14:27:27 UTC (rev 697)
+++ pkg/CHNOSZ/R/util.data.R 2022-02-07 03:46:08 UTC (rev 698)
@@ -270,7 +270,7 @@
else if(what=="SLOP98") tdata <- read.csv(system.file("extdata/OBIGT/SLOP98.csv", package="CHNOSZ"), as.is=TRUE)
else if(what=="SUPCRT92") tdata <- read.csv(system.file("extdata/OBIGT/SUPCRT92.csv", package="CHNOSZ"), as.is=TRUE)
else if(what=="AS04") tdata <- read.csv(system.file("extdata/OBIGT/AS04.csv", package="CHNOSZ"), as.is=TRUE)
- else if(what=="AkDi") tdata <- read.csv(system.file("extdata/OBIGT/AkDi.csv", package="CHNOSZ"), as.is=TRUE)
+ else if(what=="AD") tdata <- read.csv(system.file("extdata/OBIGT/AD.csv", package="CHNOSZ"), as.is=TRUE)
else if(what=="GEMSFIT") tdata <- read.csv(system.file("extdata/OBIGT/GEMSFIT.csv", package="CHNOSZ"), as.is=TRUE)
ntot <- nrow(tdata)
# where to keep the results
@@ -406,18 +406,18 @@
# remove scaling factors from EOS parameters
# and apply column names depending on the EOS
if(identical(state, "aq")) {
- # Aqueous species with abbrv = "AkDi" use the AkDi model 20210407
+ # Aqueous species with abbrv = "AD" use the AD model 20210407
abbrv <- OBIGT$abbrv
abbrv[is.na(abbrv)] <- ""
- isAkDi <- abbrv == "AkDi"
- # remove scaling factors for the HKF species, but not for the AkDi species
+ isAD <- abbrv == "AD"
+ # remove scaling factors for the HKF species, but not for the AD species
# protect this by an if statement to workaround error in subassignment to empty subset of data frame in R < 3.6.0
# (https://bugs.r-project.org/bugzilla/show_bug.cgi?id=17483) 20190302
- if(any(!isAkDi)) OBIGT[!isAkDi, 14:21] <- t(t(OBIGT[!isAkDi, 14:21]) * 10^c(-1,2,0,4,0,4,5,0))
- # for AkDi species, set NA values in remaining columns (for display only)
- if(any(isAkDi)) OBIGT[isAkDi, 17:20] <- NA
- # if all of the species are AkDi, change the variable names
- if(all(isAkDi)) colnames(OBIGT)[14:21] <- c('a','b','xi','XX1','XX2','XX3','XX4','Z')
+ if(any(!isAD)) OBIGT[!isAD, 14:21] <- t(t(OBIGT[!isAD, 14:21]) * 10^c(-1,2,0,4,0,4,5,0))
+ # for AD species, set NA values in remaining columns (for display only)
+ if(any(isAD)) OBIGT[isAD, 17:20] <- NA
+ # if all of the species are AD, change the variable names
+ if(all(isAD)) colnames(OBIGT)[14:21] <- c('a','b','xi','XX1','XX2','XX3','XX4','Z')
else colnames(OBIGT)[14:21] <- c('a1','a2','a3','a4','c1','c2','omega','Z')
} else {
OBIGT[,14:21] <- t(t(OBIGT[,14:21]) * 10^c(0,-3,5,0,-5,0,0,0))
Modified: pkg/CHNOSZ/demo/00Index
===================================================================
--- pkg/CHNOSZ/demo/00Index 2022-02-06 14:27:27 UTC (rev 697)
+++ pkg/CHNOSZ/demo/00Index 2022-02-07 03:46:08 UTC (rev 698)
@@ -27,7 +27,7 @@
TCA Standard Gibbs energies of steps of the citric acid cycle
aluminum Reactions involving Al-bearing minerals
carboxylase Rank abundance distribution for RuBisCO and acetyl-CoA carboxylase
-AkDi Henry's constant of dissolved gases
+AD Henry's constant of dissolved gases
comproportionation Gibbs energy of sulfur comproportionation
Pourbaix Eh-pH diagram for Fe-O-H with equisolubility lines
E_coli Gibbs energy of biomass synthesis in E. coli
Copied: pkg/CHNOSZ/demo/AD.R (from rev 697, pkg/CHNOSZ/demo/AkDi.R)
===================================================================
--- pkg/CHNOSZ/demo/AD.R (rev 0)
+++ pkg/CHNOSZ/demo/AD.R 2022-02-07 03:46:08 UTC (rev 698)
@@ -0,0 +1,90 @@
+# CHNOSZ/demo/AD.R
+# calculations using the Akinfiev-Diamond model 20190220
+# after Fig. 1 of Akinfiev and Diamond, 2003
+library(CHNOSZ)
+
+# Function to plot natural logarithm of Henry's constant
+lines.KH <- function(name = "CO2", T = 1:373, P = "Psat", HKF = FALSE, altH2S = FALSE) {
+ # use AD or HKF model?
+ if(!HKF) add.OBIGT("AD")
+ # use alternative parameters for H2S? (AD03 Table 1)
+ if(altH2S) mod.OBIGT("H2S", state = "aq", a = -11.2303, b = 12.6104, c = -0.2102)
+ # get properties of aq - gas reaction
+ sres <- suppressWarnings(subcrt(c(name, name), c("aq", "gas"), c(-1, 1), T = T, P = P))
+ # calculate natural logarithm of Henry's constant in mole-fraction units
+ ln_KH <- log(1000/18.0153) + log(10) * sres$out$logK
+ # plot with units of reciprocal temperature (1000/K)
+ TK <- convert(T, "K")
+ lty <- 1
+ if(altH2S) lty <- 2
+ if(HKF) lty <- 3
+ if(HKF) col <- "red" else col <- "black"
+ lines(1000/TK, ln_KH, lty = lty, col = col)
+ reset()
+}
+
+# Set up plot
+opar <- par(no.readonly = TRUE)
+par(mfrow = c(2, 2))
+par(mar = c(3.5, 3.5, 2.5, 1))
+par(mgp = c(2.4, 1, 0))
+
+ylab <- quote(ln~italic(K[H]))
+xlab <- quote(1000 / list(italic(T), K))
+
+# CO2 (Fig. 1a of AD03)
+plot(0, 0, xlim = c(1, 4), ylim = c(4, 10), xlab = xlab, ylab = ylab)
+lines.KH("CO2", 1:373, "Psat")
+lines.KH("CO2", seq(100, 650, 10), 500)
+lines.KH("CO2", 1:373, "Psat", HKF = TRUE)
+lines.KH("CO2", seq(100, 650, 10), 500, HKF = TRUE)
+dat <- read.csv(system.file("extdata/cpetc/AD03_Fig1a.csv", package = "CHNOSZ"))
+points(dat$x, dat$y, pch = dat$pch)
+text(3.5, 7.8, quote(italic(P)[sat]))
+text(3.05, 9.2, "500 bar")
+legend("bottom", c("Data (AD03, Fig. 1a)", "AD model", "HKF model"),
+ lty = c(0, 1, 3), pch = c(1, NA, NA), col = c(1, 1, 2), bty = "n")
+title(main = syslab(c("CO2", "H2O"), dash = " - "))
+
+# H2 (Fig. 1b of AD03)
+plot(0, 0, xlim = c(1, 4), ylim = c(8, 12), xlab = xlab, ylab = ylab)
+lines.KH("H2", 1:373, "Psat")
+lines.KH("H2", seq(100, 650, 10), 1000)
+lines.KH("H2", 1:373, "Psat", HKF = TRUE)
+lines.KH("H2", seq(100, 650, 10), 1000, HKF = TRUE)
+text(3.4, 11.4, quote(italic(P)[sat]))
+text(1.5, 11, "1000 bar")
+dat <- read.csv(system.file("extdata/cpetc/AD03_Fig1b.csv", package = "CHNOSZ"))
+points(dat$x, dat$y, pch = dat$pch)
+legend("bottomright", c("Data (AD03, Fig. 1b)", "AD model", "HKF model"),
+ lty = c(0, 1, 3), pch = c(1, NA, NA), col = c(1, 1, 2), bty = "n")
+title(main = syslab(c("H2", "H2O"), dash = " - "))
+
+# H2S (Fig. 1c of AD03)
+plot(0, 0, xlim = c(1, 4), ylim = c(4, 9), xlab = xlab, ylab = ylab)
+lines.KH("H2S", 1:373, "Psat")
+lines.KH("H2S", seq(100, 650, 10), 1000)
+lines.KH("H2S", 1:373, "Psat", altH2S = TRUE)
+lines.KH("H2S", seq(100, 650, 10), 1000, altH2S = TRUE)
+lines.KH("H2S", 1:373, "Psat", HKF = TRUE)
+lines.KH("H2S", seq(100, 650, 10), 1000, HKF = TRUE)
+dat <- read.csv(system.file("extdata/cpetc/AD03_Fig1c.csv", package = "CHNOSZ"))
+points(dat$x, dat$y, pch = dat$pch)
+text(3.4, 6.9, quote(italic(P)[sat]))
+text(3.1, 8.6, "1000 bar")
+legend("bottom", c("Data (AD03, Fig. 1c)", "AD model", "AD model (alt. H2S)", "HKF model"),
+ lty = c(0, 1, 2, 3), pch = c(1, NA, NA, NA), col = c(1, 1, 1, 2), bty = "n")
+title(main = syslab(c("H2S", "H2O"), dash = " - "))
+
+# CH4 (Fig. 1d of AD03)
+plot(0, 0, xlim = c(1.5, 4), ylim = c(8, 12), xlab = xlab, ylab = ylab)
+lines.KH("CH4", 1:350, "Psat")
+lines.KH("CH4", 1:350, "Psat", HKF = TRUE)
+dat <- read.csv(system.file("extdata/cpetc/AD03_Fig1d.csv", package = "CHNOSZ"))
+points(dat$x, dat$y, pch = dat$pch)
+text(3.4, 11, quote(italic(P)[sat]))
+legend("bottomright", c("Data (AD03, Fig. 1d)", "AD model", "HKF model"),
+ lty = c(0, 1, 3), pch = c(1, NA, NA), col = c(1, 1, 2), bty = "n")
+title(main = syslab(c("CH4", "H2O"), dash = " - "))
+
+par(opar)
Deleted: pkg/CHNOSZ/demo/AkDi.R
===================================================================
--- pkg/CHNOSZ/demo/AkDi.R 2022-02-06 14:27:27 UTC (rev 697)
+++ pkg/CHNOSZ/demo/AkDi.R 2022-02-07 03:46:08 UTC (rev 698)
@@ -1,86 +0,0 @@
-# CHNOSZ/demo/AkDi.R
-# calculations using the Akinfiev-Diamond model 20190220
-# after Fig. 1 of Akinfiev and Diamond, 2003
-library(CHNOSZ)
-
-# function to plot natural logarithm of Henry's constant
-lines.KH <- function(name = "CO2", T = 1:373, P = "Psat", HKF = FALSE, altH2S = FALSE) {
- # use AkDi or HKF model?
- if(!HKF) add.OBIGT("AkDi")
- # use alternative parameters for H2S? (AD03 Table 1)
- if(altH2S) mod.OBIGT("H2S", state="aq", a=-11.2303, b=12.6104, c=-0.2102)
- # get properties of aq - gas reaction
- sres <- suppressWarnings(subcrt(c(name, name), c("aq", "gas"), c(-1, 1), T = T, P = P))
- # calculate natural logarithm of Henry's constant in mole-fraction units
- ln_KH <- log(1000/18.0153) + log(10) * sres$out$logK
- # plot with units of reciprocal temperature (1000/K)
- TK <- convert(T, "K")
- lty <- 1
- if(altH2S) lty <- 2
- if(HKF) lty <- 3
- if(HKF) col <- "red" else col <- "black"
- lines(1000/TK, ln_KH, lty = lty, col = col)
- reset()
-}
-
-# set up plot
-opar <- par(no.readonly = TRUE)
-par(mfrow = c(2, 2))
-par(mar = c(3.5, 3.5, 2.5, 1))
-par(mgp = c(2.4, 1, 0))
-
-ylab <- quote(ln~italic(K[H]))
-xlab <- quote(1000 / list(italic(T), K))
-
-# CO2 (Fig. 1a of AD03)
-plot(0, 0, xlim=c(1, 4), ylim=c(4, 10), xlab=xlab, ylab=ylab)
-lines.KH("CO2", 1:373, "Psat")
-lines.KH("CO2", seq(100, 650, 10), 500)
-lines.KH("CO2", 1:373, "Psat", HKF = TRUE)
-lines.KH("CO2", seq(100, 650, 10), 500, HKF = TRUE)
-dat <- read.csv(system.file("extdata/cpetc/AD03_Fig1a.csv", package="CHNOSZ"))
-points(dat$x, dat$y, pch=dat$pch)
-text(3.5, 7.8, quote(italic(P)[sat]))
-text(3.05, 9.2, "500 bar")
-legend("bottom", c("Data (AD03, Fig. 1a)", "AkDi model", "HKF model"), lty=c(0, 1, 3), pch=c(1, NA, NA), col=c(1, 1, 2), bty="n")
-title(main=syslab(c("CO2", "H2O"), dash = " - "))
-
-# H2 (Fig. 1b of AD03)
-plot(0, 0, xlim=c(1, 4), ylim=c(8, 12), xlab=xlab, ylab=ylab)
-lines.KH("H2", 1:373, "Psat")
-lines.KH("H2", seq(100, 650, 10), 1000)
-lines.KH("H2", 1:373, "Psat", HKF = TRUE)
-lines.KH("H2", seq(100, 650, 10), 1000, HKF = TRUE)
-text(3.4, 11.4, quote(italic(P)[sat]))
-text(1.5, 11, "1000 bar")
-dat <- read.csv(system.file("extdata/cpetc/AD03_Fig1b.csv", package="CHNOSZ"))
-points(dat$x, dat$y, pch=dat$pch)
-legend("bottomright", c("Data (AD03, Fig. 1b)", "AkDi model", "HKF model"), lty=c(0, 1, 3), pch=c(1, NA, NA), col=c(1, 1, 2), bty="n")
-title(main=syslab(c("H2", "H2O"), dash = " - "))
-
-# H2S (Fig. 1c of AD03)
-plot(0, 0, xlim=c(1, 4), ylim=c(4, 9), xlab=xlab, ylab=ylab)
-lines.KH("H2S", 1:373, "Psat")
-lines.KH("H2S", seq(100, 650, 10), 1000)
-lines.KH("H2S", 1:373, "Psat", altH2S = TRUE)
-lines.KH("H2S", seq(100, 650, 10), 1000, altH2S = TRUE)
-lines.KH("H2S", 1:373, "Psat", HKF = TRUE)
-lines.KH("H2S", seq(100, 650, 10), 1000, HKF = TRUE)
-dat <- read.csv(system.file("extdata/cpetc/AD03_Fig1c.csv", package="CHNOSZ"))
-points(dat$x, dat$y, pch=dat$pch)
-text(3.4, 6.9, quote(italic(P)[sat]))
-text(3.1, 8.6, "1000 bar")
-legend("bottom", c("Data (AD03, Fig. 1c)", "AkDi model", "AkDi model (alt. H2S)", "HKF model"), lty=c(0, 1, 2, 3), pch=c(1, NA, NA, NA), col=c(1, 1, 1, 2), bty="n")
-title(main=syslab(c("H2S", "H2O"), dash = " - "))
-
-# CH4 (Fig. 1d of AD03)
-plot(0, 0, xlim=c(1.5, 4), ylim=c(8, 12), xlab=xlab, ylab=ylab)
-lines.KH("CH4", 1:350, "Psat")
-lines.KH("CH4", 1:350, "Psat", HKF = TRUE)
-dat <- read.csv(system.file("extdata/cpetc/AD03_Fig1d.csv", package="CHNOSZ"))
-points(dat$x, dat$y, pch=dat$pch)
-text(3.4, 11, quote(italic(P)[sat]))
-legend("bottomright", c("Data (AD03, Fig. 1d)", "AkDi model", "HKF model"), lty=c(0, 1, 3), pch=c(1, NA, NA), col=c(1, 1, 2), bty="n")
-title(main=syslab(c("CH4", "H2O"), dash = " - "))
-
-par(opar)
Modified: pkg/CHNOSZ/demo/aluminum.R
===================================================================
--- pkg/CHNOSZ/demo/aluminum.R 2022-02-06 14:27:27 UTC (rev 697)
+++ pkg/CHNOSZ/demo/aluminum.R 2022-02-07 03:46:08 UTC (rev 698)
@@ -29,7 +29,7 @@
r2 <- subcrt(c("boehmite", "H2O", "SiO2", "kaolinite"), c(-1, -0.5, -1, 0.5), T = T, P = 1000, exceed.Ttr = TRUE)
reset()
## third calculation: get Si(OH)4 from Akinfiev and Plyasunov, 2014
-add.OBIGT("AkDi")
+add.OBIGT("AD")
r3 <- subcrt(c("boehmite", "Si(OH)4", "H2O", "kaolinite"), c(-1, -1, 1.5, 0.5), T = T, P = 1000, exceed.Ttr = TRUE)
reset()
## fourth calculation: minerals as in SUPCRT92
@@ -48,7 +48,7 @@
title(main = describe.reaction(r4$reaction), cex.main = 1.1)
legend("bottomright", lty = c(0, 0, 1, 2, 3, 0), pch = c(1, 4, NA, NA, NA, NA), lwd = c(1, 1, 1.5, 1, 1, 0),
col = c("black", "red", "black", "red", "purple", NA), bty = "n", cex = 0.9,
- legend = c("Hemley et al., 1980", "SUPCRTBL", "CHNOSZ", 'add.OBIGT("AS04")', 'add.OBIGT("AkDi")', ""))
+ legend = c("Hemley et al., 1980", "SUPCRTBL", "CHNOSZ", 'add.OBIGT("AS04")', 'add.OBIGT("AD")', ""))
legend("bottomright", lty = 4, pch = NA, lwd = 1, col = "blue", legend = 'add.OBIGT("SUPCRT92")', bty = "n", cex = 0.9)
legend("topleft", c("Boehmite - Kaolinite", "After Zhu and Lu, 2009 Fig. A1"), bty = "n")
# Helgeson et al., 1978 (HDNB78): http://www.worldcat.org/oclc/13594862
@@ -111,7 +111,7 @@
lines(invTK, pK, col = "red", lty = 2)
reset()
# plot line: Si(OH)4 from Akinfiev and Plyasunov, 2014
-add.OBIGT("AkDi")
+add.OBIGT("AD")
sres <- subcrt(c("kaolinite", "OH-", "H2O", "Al(OH)4-", "Si(OH)4"), c(-1, -2, -5, 2, 2), T = T)
pK <- -sres$out$logK
lines(invTK, pK, col = "purple", lty = 3)
@@ -133,7 +133,7 @@
legend("topright", c("Kaolinite solubility", "After Tutolo et al., 2014 Fig. 2"), bty = "n")
legend("bottomleft", lty = c(0, 0, 0, 1, 2, 3, 4), pch = c(1, NA, 4, NA, NA, NA, NA),
lwd = c(1, 1, 1, 1.5, 1, 1, 1), col = c("black", "black", "red", "black", "red", "purple", "blue"),
- legend = c("Various sources \u2013", " see Tutolo et al., 2014", "SUPCRTBL", "CHNOSZ", 'add.OBIGT("AS04")', 'add.OBIGT("AkDi")', 'add.OBIGT("SUPCRT92")'),
+ legend = c("Various sources \u2013", " see Tutolo et al., 2014", "SUPCRTBL", "CHNOSZ", 'add.OBIGT("AS04")', 'add.OBIGT("AD")', 'add.OBIGT("SUPCRT92")'),
bty = "n", cex = 0.9)
###########
Modified: pkg/CHNOSZ/demo/sources.R
===================================================================
--- pkg/CHNOSZ/demo/sources.R 2022-02-06 14:27:27 UTC (rev 697)
+++ pkg/CHNOSZ/demo/sources.R 2022-02-07 03:46:08 UTC (rev 698)
@@ -21,7 +21,7 @@
tdata <- read.csv(system.file("extdata/OBIGT/AS04.csv", package="CHNOSZ"), as.is=TRUE)
os7 <- gsub("\ .*", "", tdata$ref1)
os8 <- gsub("\ .*", "", tdata$ref2)
-tdata <- read.csv(system.file("extdata/OBIGT/AkDi.csv", package="CHNOSZ"), as.is=TRUE)
+tdata <- read.csv(system.file("extdata/OBIGT/AD.csv", package="CHNOSZ"), as.is=TRUE)
os9 <- gsub("\ .*", "", tdata$ref1)
os10 <- gsub("\ .*", "", tdata$ref2)
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/chnosz -r 698
More information about the CHNOSZ-commits
mailing list