[CHNOSZ-commits] r715 - in pkg/CHNOSZ: . demo inst inst/tinytest man vignettes
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sat Mar 26 09:11:31 CET 2022
Author: jedick
Date: 2022-03-26 09:11:30 +0100 (Sat, 26 Mar 2022)
New Revision: 715
Removed:
pkg/CHNOSZ/man/eos.Rd
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/NAMESPACE
pkg/CHNOSZ/demo/DEW.R
pkg/CHNOSZ/demo/copper.R
pkg/CHNOSZ/inst/NEWS.Rd
pkg/CHNOSZ/inst/tinytest/test-eos.R
pkg/CHNOSZ/man/CHNOSZ-package.Rd
pkg/CHNOSZ/man/subcrt.Rd
pkg/CHNOSZ/man/thermo.Rd
pkg/CHNOSZ/man/water.Rd
pkg/CHNOSZ/vignettes/anintro.Rmd
pkg/CHNOSZ/vignettes/customizing.Rmd
pkg/CHNOSZ/vignettes/vig.bib
Log:
Remove exports of cgl(), hkf(), and AD()
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2022-03-25 11:13:38 UTC (rev 714)
+++ pkg/CHNOSZ/DESCRIPTION 2022-03-26 08:11:30 UTC (rev 715)
@@ -1,6 +1,6 @@
-Date: 2022-03-25
+Date: 2022-03-26
Package: CHNOSZ
-Version: 1.4.3-7
+Version: 1.4.3-8
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-25 11:13:38 UTC (rev 714)
+++ pkg/CHNOSZ/NAMESPACE 2022-03-26 08:11:30 UTC (rev 715)
@@ -23,7 +23,7 @@
"count.aa", "nucleic.complement", "nucleic.formula",
"rho.IAPWS95", "IAPWS95", "water.AW90", "WP02.auxiliary", "water.IAPWS95",
"getrank", "parent", "sciname", "allparents", "getnodes", "getnames",
- "protein.OBIGT", "hkf", "cgl", "which.pmax",
+ "protein.OBIGT", "which.pmax",
"equil.boltzmann", "equil.reaction", "find.tp",
"ionize.aa", "MP90.cp", "aasum",
"qqr", "RMSD", "CVRMSD", "spearman", "DGmix", "DDGmix", "DGtr",
@@ -54,7 +54,7 @@
# added 20171121 or later
"dumpdata", "thermo.axis", "solubility", "NaCl",
# added 20190213 or later
- "CHNOSZ", "thermo", "reset", "OBIGT", "retrieve", "AD", "moles",
+ "CHNOSZ", "thermo", "reset", "OBIGT", "retrieve", "moles",
"lNaCl", "lS", "lT", "lP", "lTP", "lex",
# added 20200716 or later
"mash", "mix", "rebalance",
Modified: pkg/CHNOSZ/demo/DEW.R
===================================================================
--- pkg/CHNOSZ/demo/DEW.R 2022-03-25 11:13:38 UTC (rev 714)
+++ pkg/CHNOSZ/demo/DEW.R 2022-03-26 08:11:30 UTC (rev 715)
@@ -1,24 +1,24 @@
-# demo for the Deep Earth Water (DEW) model in CHNOSZ 20170927
+# Demo for the Deep Earth Water (DEW) model in CHNOSZ 20170927
library(CHNOSZ)
-# set up subplots
+# Set up subplots
opar <- par(no.readonly = TRUE)
par(mfrow = c(2, 2), mar=c(3.0, 3.5, 2.5, 1.0), mgp=c(1.7, 0.3, 0), las=1, tcl=0.3, xaxs="i", yaxs="i")
-# activate DEW model
+# Activate DEW model
oldwat <- water("DEW")
###########
-#### plot 1: quartz solubility at high pressure
+#### Plot 1: quartz solubility at high pressure
## after Figure 7D of Sverjensky et al., 2014a [SHA14]
## (Geochim. Cosmochim. Acta, https://doi.org/10.1016/j.gca.2013.12.019)
###########
-# load SiO2 and Si2O4 data taken from DEW spreadsheet
+# Load SiO2 and Si2O4 data taken from DEW spreadsheet
iSi <- add.OBIGT("DEW", c("SiO2", "Si2O4"))
-# print the data references to confirm we got the right ones
+# Print the data references to confirm we got the right ones
thermo.refs(iSi)
-# set temperature ranges for different pressures
+# Set temperature ranges for different pressures
# data.frame is used to make P and T the same length
PT0.5 <- data.frame(P=500, T=seq(200, 550, 10))
PT1.0 <- data.frame(P=1000, T=seq(200, 700, 10))
@@ -27,11 +27,11 @@
PT10.0 <- data.frame(P=10000, T=seq(200, 825, 10))
PT20.0 <- data.frame(P=20000, T=seq(200, 800, 10))
PT <- rbind(PT0.5, PT1.0, PT2.0, PT5.0, PT10.0, PT20.0)
-# reaction 1: quartz = SiO2(aq) [equivalent to quartz + 3 H2O = Si(OH)4]
+# Reaction 1: quartz = SiO2(aq) [equivalent to quartz + 3 H2O = Si(OH)4]
SiO2_logK <- suppressWarnings(subcrt(c("quartz", "SiO2"), c("cr", "aq"), c(-1, 1), P=PT$P, T=PT$T))$out$logK
-# reaction 2: 2 quartz = Si2O4(aq) [equivalent to 2 quartz + 3 H2O = Si2O(OH)6]
+# Reaction 2: 2 quartz = Si2O4(aq) [equivalent to 2 quartz + 3 H2O = Si2O(OH)6]
Si2O4_logK <- suppressWarnings(subcrt(c("quartz", "Si2O4"), c("cr", "aq"), c(-2, 1), P=PT$P, T=PT$T))$out$logK
-# plot the sum of molalities (== activities) for each pressure
+# Plot the sum of molalities (== activities) for each pressure
plot(c(200, 1000), c(-2.5, 0.5), type="n", xlab=axis.label("T"), ylab="log molality")
for(P in unique(PT$P)) {
icond <- PT$P == P
@@ -50,54 +50,59 @@
# TODO: lines are a little low at highest P and T ...
###########
-#### plot 2: correlations between non-solvation volume and HKF a1 parameter
+#### Plot 2: correlations between non-solvation volume and HKF a1 parameter
## after Figures 12B and 12C of Sverjensky et al., 2014a [SHA14]
###########
-# load the fitted parameters for species as used by SHA14
+# For this plot we don't need the DEW water model
+reset()
+# Load the fitted parameters for species as used by SHA14
# TODO: also use their Ca+2??
# NOTE: don't load NaCl, NH4+, or HS- here because the DEW spreadsheet lists a1 from the correlation
add.OBIGT("DEW", c("CO3-2", "BO2-", "MgCl+", "SiO2", "HCO3-", "Si2O4"))
-# set up the plot
+# Set up the plot
V0nlab <- expression(Delta * italic(V) * degree[n]~~(cm^3~mol^-1))
a1lab <- expression(italic(a)[1]%*%10~~(cal~mol~bar^-1))
plot(c(-25, 50), c(-4, 12), type="n", xlab=V0nlab, ylab=a1lab)
-# a function to get the HKF parameters, calculate nonsolvation volume, plot points, labels, error bars, and correlation lines
+# A function to get the HKF parameters, calculate nonsolvation volume, plot points, labels, error bars, and correlation lines
plotfun <- function(species, col, pch, cex, dy, error, xlim, corrfun) {
- # get HKF parameters
- par <- info(info(species), check.it = FALSE)
- a1 <- par$a1 * 10
- # get the nonsolvation volume
- Vn <- unlist(hkf("V", par, contrib="n")$aq)
+ # Get HKF a1 parameter
+ a1 <- info(info(species), check.it = FALSE)$a1 * 10
+ # Set omega to zero to calculate non-solvation volume 20220326
+ mod.OBIGT(species, omega = 0)
+ Vn <- do.call(rbind, subcrt(species, T = 25)$out)$V
points(Vn, a1, col=col, pch=pch, cex=cex)
for(i in 1:length(species)) text(Vn[i], a1[i]+dy, expr.species(species[i]))
arrows(Vn, a1 - error, Vn, a1 + error, length = 0.03, angle = 90, code = 3, col=col)
lines(xlim, corrfun(xlim), col=col)
}
-# monovalent ions: Na+, K+, Cl-, Br-
+# Monovalent ions: Na+, K+, Cl-, Br-
monofun <- function(Vn) 2.0754 + 0.10871 * Vn
-# for easier reading, set y-offset to NA so the labels aren't plotted
+# For easier reading, set y-offset to NA so the labels aren't plotted
plotfun(c("Na+", "K+", "Cl-", "Br-"), "red", 19, 1, NA, 0.5, c(-7, 35), monofun)
-# divalent ions: Mg+2, Ca+2, CO3-2, SO4-2
+# Divalent ions: Mg+2, Ca+2, CO3-2, SO4-2
difun <- function(Vn) 3.5321 + 0.23911 * Vn
plotfun(c("Mg+2", "Ca+2", "CO3-2", "SO4-2"), "black", 15, 1, 1.2, 0.7, c(-20, 25), difun)
-# complexes and neutral molecules: BO2-, MgCl+, SiO2, NaCl, HCO3-, Si2O4, NH4+, HS-
+# Complexes and neutral molecules: BO2-, MgCl+, SiO2, NaCl, HCO3-, Si2O4, NH4+, HS-
compfun <- function(Vn) 1.5204 + 0.19421 * Vn
plotfun(c("MgCl+", "SiO2", "NaCl", "HCO3-", "Si2O4"), "blue1", 18, 1.5, 1, 0.5, c(-20, 50), compfun)
-# for easier reading, put some labels below the points
+# For easier reading, put some labels below the points
plotfun(c("BO2-", "NH4+", "HS-"), "blue1", 18, 1.5, -1.2, 0.5, c(-20, 50), compfun)
-# include an empty subscript for better spacing between the lines
+# Include an empty subscript for better spacing between the lines
t1 <- quote("Correlations between non-solvation"[])
t2 <- quote("volume and HKF "*italic(a)[1]*" parameter")
mtitle(as.expression(c(t1, t2)))
+# Clean up for next plot
+OBIGT()
+
###########
-#### plot 3: logfO2-pH diagram for aqueous inorganic and organic carbon species at high pressure
+#### Plot 3: logfO2-pH diagram for aqueous inorganic and organic carbon species at high pressure
## after Figure 1b of Sverjensky et al., 2014b [SSH14]
## (Nature Geoscience, https://doi.org/10.1038/NGEO2291)
###########
-# define system
+# Define system
basis("CHNOS+")
inorganic <- c("CO2", "HCO3-", "CO3-2", "CH4")
organic <- c("formic acid", "formate", "acetic acid", "acetate", "propanoic acid", "propanoate")
@@ -104,10 +109,10 @@
species(c(inorganic, organic), 0)
fill <- c(rep("aliceblue", length(inorganic)), rep("aquamarine", length(organic)))
-# a function to make the diagrams
+# A function to make the diagrams
dfun <- function(T = 600, P = 50000, res = 300) {
a <- affinity(pH = c(0, 10, res), O2 = c(-24, -12, res), T = T, P = P)
- # set total C concentration to 0.03 molal
+ # Set total C concentration to 0.03 molal
# (from EQ3NR model for eclogite [Supporting Information of SSH14])
e <- equilibrate(a, loga.balance = log10(0.03))
diagram(e, fill = fill)
@@ -115,11 +120,11 @@
legend("bottomleft", legend=dp, bty="n")
}
-OBIGT()
-## (not run: make diagram using CHNOSZ default database;
-## not recommended for high P)
+water("DEW")
+## Not run: make diagram using CHNOSZ default database
+## (not recommended for high P)
#dfun()
-# make diagram using CO2, HCO3-, CO3-2, CH4, and acetic acid data from DEW spreadsheet
+# Make diagram using CO2, HCO3-, CO3-2, CH4, and acetic acid data from DEW spreadsheet
# (the acetate field disappears if we also use the data for acetate from the spreadsheet 20200629)
#add.OBIGT("DEW", c("CO2", "HCO3-", "CO3-2", "CH4", "acetic acid"))
add.OBIGT("DEW")
@@ -127,11 +132,11 @@
mtitle(c("Inorganic and organic species", "C[total] = 0.03 molal"))
###########
-#### plot 4: speciation of carbon as a function T, logfO2 and pH (added 20171008)
+#### Plot 4: speciation of carbon as a function T, logfO2 and pH (added 20171008)
## after SSH14 Fig. 3
###########
-# conditions:
+# Conditions:
# T = 600, 700, 800, 900, 1000 degC
# P = 5.0GPa (50000 bar)
# fO2 = QFM - 2
@@ -140,33 +145,33 @@
# dissolved carbon: 0.03, 0.2, 1, 4, 20 molal
# true ionic strength: 0.39, 0.57, 0.88, 1.45, 2.49
# pH: 3.80, 3.99, 4.14, 4.25, 4.33
-## activate DEW model
+## Activate DEW model
reset()
water("DEW")
-# add species data for DEW
+# Add species data for DEW
inorganics <- c("CH4", "CO2", "HCO3-", "CO3-2")
organics <- c("formic acid", "formate", "acetic acid", "acetate", "propanoic acid", "propanoate")
add.OBIGT("DEW")
-## set basis species
+## Set basis species
basis(c("Fe", "SiO2", "CO3-2", "H2O", "oxygen", "H+"))
-## calculate logfO2 in QFM buffer
+## Calculate logfO2 in QFM buffer
basis("O2", "QFM")
T <- seq(600, 1000, 100)
buf <- affinity(T=T, P=50000, return.buffer=TRUE)
-## add species
+## Add species
species(c(inorganics, organics))
-## generate spline functions from IS, pH, and molC values at every 100 degC
+## Generate spline functions from IS, pH, and molC values at every 100 degC
IS <- c(0.39, 0.57, 0.88, 1.45, 2.49)
pH <- c(3.80, 3.99, 4.14, 4.25, 4.33)
molC <- c(0.03, 0.2, 1, 4, 20)
-## use extended Debye-Huckel equation with b_gamma set to zero
+## Use extended Debye-Huckel equation with b_gamma set to zero
nonideal("bgamma0")
-## calculate affinities on the T-logfO2-pH-IS transect
+## Calculate affinities on the T-logfO2-pH-IS transect
a <- affinity(T = T, O2 = buf$O2 - 2, IS = IS, pH = pH, P = 50000)
-## calculate metastable equilibrium activities using the total
-## carbon molality as an approximation of total activity
+## Calculate metastable equilibrium activities using the total
+## Carbon molality as an approximation of total activity
e <- equilibrate(a, loga.balance = log10(molC))
-## make the diagram; don't plot names of low-abundance species
+## Make the diagram; don't plot names of low-abundance species
names <- c(inorganics, organics)
names[c(4, 5, 7, 9)] <- ""
col <- rep("black", length(names))
@@ -177,7 +182,7 @@
diagram(e, alpha = "balance", names = names, col = col, ylim = c(0, 0.8), ylab="carbon fraction")
}
-## add legend and title
+## Add legend and title
ltxt1 <- "P = 50000 bar"
ltxt2 <- substitute(logfO2=="QFM-2", list(logfO2 = axis.label("O2")))
pH <- seq(3.8, 4.3, length.out = length(T))
@@ -186,11 +191,11 @@
t2 <- "after Sverjensky et al., 2014b"
mtitle(c(t1, t2))
-### additional checks
+### Additional checks
# The maximum absolute pairwise difference between x and y
maxdiff <- function(x, y) max(abs(y - x))
-## check that we're within 0.1 of the QFM-2 values used by SSH14
+## Check that we're within 0.1 of the QFM-2 values used by SSH14
stopifnot(maxdiff((buf$O2-2), c(-17.0, -14.5, -12.5, -10.8, -9.4)) < 0.1)
# Here are the logKs of aqueous species dissociation reactions at 600 degC and 50000 bar,
@@ -197,22 +202,22 @@
# values from EQ3NR output in Supporting Information of the paper (p. 103-109):
inorganic.logK <- c(24.4765, -9.0784, -5.3468, 0)
organic.logK <- c(1.7878, 2.5648, 15.3182, 16.9743, 30.4088, 28.9185)
-# calculate equilibrium constants of the reactions in CHNOSZ; use a negative sign to change from formation to dissociation
+# Calculate equilibrium constants of the reactions in CHNOSZ; use a negative sign to change from formation to dissociation
logK.calc <- -unlist(affinity(T = 600, P = 50000, property = "logK")$values)
logK.calc - c(inorganic.logK, organic.logK)
-## except for acetate, we're within 0.021 of the logK values used by SSH14
+## Except for acetate, we're within 0.021 of the logK values used by SSH14
stopifnot(maxdiff(logK.calc[-8], c(inorganic.logK, organic.logK)[-8]) < 0.021)
-## check that we get similar activity coefficients
-# activity coefficients for monovalent species from EQ3NR output
+## Check that we get similar activity coefficients
+# Activity coefficients for monovalent species from EQ3NR output
loggamma <- c(-0.15, -0.18, -0.22, -0.26, -0.31)
-# activity coefficients calculated in CHNOSZ
+# Activity coefficients calculated in CHNOSZ
sres <- subcrt("propanoate", T = seq(600, 1000, 100), P = 50000, IS = c(0.39, 0.57, 0.88, 1.45, 2.49))
stopifnot(maxdiff(sres$out[[1]]$loggam, loggamma) < 0.023)
-# if m_star in nonideal() was zero, we could decrease the tolerance here
+# If m_star in nonideal() was zero, we could decrease the tolerance here
#stopifnot(maxdiff(sres$out[[1]]$loggam, loggamma) < 0.004)
-# reset OBIGT database
+# Reset OBIGT database
reset()
par(opar)
Modified: pkg/CHNOSZ/demo/copper.R
===================================================================
--- pkg/CHNOSZ/demo/copper.R 2022-03-25 11:13:38 UTC (rev 714)
+++ pkg/CHNOSZ/demo/copper.R 2022-03-26 08:11:30 UTC (rev 715)
@@ -4,14 +4,14 @@
## solutions. J. Electrochem. Soc., 148, B51-B57. doi:10.1149/1.1344532)
library(CHNOSZ)
-# we need data for Cu-Gly complexes 20190206
+# We need data for Cu-Gly complexes 20190206
add.OBIGT(system.file("extdata/adds/SK95.csv", package = "CHNOSZ"))
-# add some new species to thermo()$OBIGT
-m1 <- makeup(info(c("Cu+", "glycinate", "glycinate")), sum=TRUE)
-mod.OBIGT(name="Cu(Gly)2-", formula=as.chemical.formula(m1))
-m2 <- makeup(info(c("Cu+2", "glycinate", "H+")), sum=TRUE)
-mod.OBIGT(name="HCu(Gly)+2", formula=as.chemical.formula(m2))
-# Gibbs energies from A&D's Table 1 and Table II
+# Add some new species to thermo()$OBIGT
+m1 <- makeup(info(c("Cu+", "glycinate", "glycinate")), sum = TRUE)
+mod.OBIGT(name = "Cu(Gly)2-", formula = as.chemical.formula(m1))
+m2 <- makeup(info(c("Cu+2", "glycinate", "H+")), sum = TRUE)
+mod.OBIGT(name = "HCu(Gly)+2", formula = as.chemical.formula(m2))
+# Names of species in AD01 Table 1 and Table II
Cu_s <- c("copper", "cuprite", "tenorite")
Gly <- c("glycinium", "glycine", "glycinate")
Cu_aq <- c("Cu+", "Cu+2", "CuO2-2", "HCuO2-")
@@ -18,52 +18,57 @@
CuGly <- c("Cu(Gly)+", "Cu(Gly)2", "Cu(Gly)2-", "HCu(Gly)+2")
names <- c(Cu_s, Gly, Cu_aq, CuGly)
G <- c(
- convert(c(0, -146, -129.7,
+ # Table I: Gibbs energies in kJ/mol
+ c(0, -146, -129.7,
-384.061, -370.647, -314.833,
- 49.98, 65.49, -183.6, -258.5, -298.2)*1000, "cal"),
- convert(c(15.64, 10.1, 2.92), "G"))
-# run updates in order so later species take account of prev. species' values
+ 49.98, 65.49, -183.6, -258.5, -298.2)*1000,
+ # Table II: Association constants, converted to Gibbs energy
+ convert(c(15.64, 10.1, 2.92), "G")
+)
+
+# Run updates in order so later species take account of prev. species' values
getG <- function(x) info(info(x))$G
for(i in 1:length(G)) {
myG <- G[i]
- if(i==12) myG <- myG + getG("Cu+2") + 2*getG("glycinate")
- if(i==13) myG <- myG + getG("Cu+") + 2*getG("glycinate")
- if(i==14) myG <- myG + getG("Cu(Gly)+")
- mod.OBIGT(names[i], G=myG)
+ if(i == 12) myG <- myG + getG("Cu+2") + 2*getG("glycinate")
+ if(i == 13) myG <- myG + getG("Cu+") + 2*getG("glycinate")
+ if(i == 14) myG <- myG + getG("Cu(Gly)+")
+ # Energies are in Joules, so we have to change units of species in default OBIGT 20220326
+ mod.OBIGT(names[i], G = myG, E_units = "J")
}
-# in Fig. 2b, total log activities of Cu (Cu_T) and glycine (L_T) are -4 and -1
+# In Fig. 2b, total log activities of Cu (Cu_T) and glycine (L_T) are -4 and -1
basis(c("Cu+2", "H2O", "H+", "e-", "glycinium", "CO2"), c(999, 0, 999, 999, -1, 999))
-# add solids and aqueous species
+# Add solids and aqueous species
species(Cu_s)
species(c(Cu_aq, CuGly), -4, add = TRUE)
names <- c(Cu_s, Cu_aq, CuGly)
-# mosaic diagram with speciate glycine as a function of pH
-m <- mosaic(bases=Gly, pH=c(0, 16, 500), Eh=c(-0.6, 1.0, 500))
+# Mosaic diagram with glycine speciation as a function of pH
+m <- mosaic(bases = Gly, pH = c(0, 16, 500), Eh = c(-0.6, 1.0, 500))
fill <- c(rep("lightgrey", 3), rep("white", 4), rep("lightblue", 4))
-d <- diagram(m$A.species, fill=fill, names=FALSE, xaxs="i", yaxs="i", fill.NA="pink2", limit.water = TRUE)
-# to make the labels look nicer
+d <- diagram(m$A.species, fill = fill, names = FALSE, xaxs = "i", yaxs = "i", fill.NA = "pink2", limit.water = TRUE)
+# Adjustments for labels
names <- names[sort(unique(as.numeric(d$predominant)))]
for(i in 1:length(names)) {
if(i %in% 1:3) lab <- names[i] else lab <- expr.species(names[i])
- # some manual adjustment so labels don't collide
+ # Some manual adjustment so labels don't collide
srt <- dy <- dx <- 0
- if(names[i]=="tenorite") dy <- -0.1
- if(names[i]=="CuO2-2") dy <- -0.1
- if(names[i]=="HCu(Gly)+2") srt <- 90
- if(names[i]=="HCu(Gly)+2") dx <- -0.2
- if(names[i]=="Cu(Gly)+") srt <- 90
- text(na.omit(d$namesx)[i]+dx, na.omit(d$namesy)[i]+dy, lab, srt=srt)
+ if(names[i] == "tenorite") dy <- -0.1
+ if(names[i] == "CuO2-2") dy <- -0.1
+ if(names[i] == "HCu(Gly)+2") srt <- 90
+ if(names[i] == "HCu(Gly)+2") dx <- -0.2
+ if(names[i] == "Cu(Gly)+") srt <- 90
+ text(na.omit(d$namesx)[i]+dx, na.omit(d$namesy)[i]+dy, lab, srt = srt)
}
-# add glycine ionization lines
-d <- diagram(m$A.bases, add=TRUE, col="darkblue", lty=3, names=FALSE)
-text(d$namesx, -0.5, Gly, col="darkblue")
+# Add glycine ionization lines
+d <- diagram(m$A.bases, add = TRUE, col = "darkblue", lty = 3, names = FALSE)
+text(d$namesx, -0.5, Gly, col = "darkblue")
-# add water lines and title
+# Add water lines and title
water.lines(d)
mtitle(expression("Copper-water-glycine at 25"~degree*"C and 1 bar",
"After Aksu and Doyle, 2001 (Fig. 2b)"))
-# done!
+# Done!
reset()
Modified: pkg/CHNOSZ/inst/NEWS.Rd
===================================================================
--- pkg/CHNOSZ/inst/NEWS.Rd 2022-03-25 11:13:38 UTC (rev 714)
+++ pkg/CHNOSZ/inst/NEWS.Rd 2022-03-26 08:11:30 UTC (rev 715)
@@ -18,13 +18,12 @@
\item Units of Joules instead of calories are now used by default for the
thermodynamic properties output by \code{subcrt()}. That is,
\code{E.units("J")} is the default setting. User scripts that implicitly
- depend on the previous default setting of \code{E.units("cal")} need to
- be modified to produce expected output.
+ depend on the previous default setting of \code{E.units("cal")} may need
+ to be modified to produce expected output.
- \item The switch to Joules is work in progress. As much as possible,
- functions in CHNOSZ have been modified to use Joules in internal
- variables, and the gas constant has been changed to Joules in
- \code{convert()} and other functions.
+ \item As much as possible, functions in CHNOSZ have been modified to use
+ Joules in internal variables, and the gas constant has been changed to
+ Joules in \code{convert()} and other functions.
\item The documentation is still catching up ...
@@ -53,6 +52,10 @@
\item The eos-regress.Rmd vignette has been removed to debug problems
associated with the switch to Joules.
+ \item Remove exports of \code{cgl}, \code{hkf}, and \code{AD}.
+ \code{subcrt} should be used for all calculations of thermodynamic
+ properties.
+
}
}
Modified: pkg/CHNOSZ/inst/tinytest/test-eos.R
===================================================================
--- pkg/CHNOSZ/inst/tinytest/test-eos.R 2022-03-25 11:13:38 UTC (rev 714)
+++ pkg/CHNOSZ/inst/tinytest/test-eos.R 2022-03-26 08:11:30 UTC (rev 715)
@@ -1,8 +1,8 @@
# Load default settings for CHNOSZ
reset()
-info <- "cgl() with NA volume produces reasonable results"
-# 20171006: after rewriting much of the code in cgl(), melanterite and hydronium jarosite
+info <- "CGL species with NA volume produces reasonable results"
+# 20171006: After rewriting much of the code in cgl(), melanterite and hydronium jarosite
# disappeared from the example diagram after Majzlan et al., 2006.
# Because the volumes are NA, the integral properties became NA,
# but it makes more sense to set them to zero.
@@ -10,44 +10,43 @@
expect_equal(info(ispecies)$V, NA_real_, info = info)
sout <- subcrt(ispecies, T = c(25, 25, 100, 100), P = c(1, 100, 1, 100))$out[[1]]
expect_false(any(is.na(sout$H)), info = info)
-# for melanterite, which is listed in the database with zero heat capacity,
+# For melanterite, which is listed in the database with zero heat capacity,
# all Cp and V integrals evaluate to zero
expect_true(length(unique(sout$H)) == 1, info = info)
info <- "gfun() gives expected results"
-# calculate values of g and its derivatives
-# up to 350 degrees C at Psat
+# Calculate values of g and its derivatives up to350 degrees C at Psat
Tc <- c(0, 25, 50, 100, 150, 200, 250, 300, 350)
-# get the required properties of water
+# Get the required properties of water
w <- water(c("rho", "alpha", "daldT", "beta", "Psat"), T = convert(Tc, "K"), P = "Psat")
-# calculate g and its derivatives
+# Calculate g and its derivatives
gfun.Psat <- CHNOSZ:::gfun(w$rho/1000, Tc, w$Psat, w$alpha, w$daldT, w$beta)
-# up to 450 degrees C at 500 bar
+# Up to 450 degrees C at 500 bar
Tc <- c(Tc, 400, 450)
w <- water(c("rho", "alpha", "daldT", "beta"), T = convert(Tc, "K"), P = 500)
gfun.500 <- CHNOSZ:::gfun(w$rho/1000, Tc, rep(500, length(Tc)), w$alpha, w$daldT, w$beta)
-# up to 600 degrees C at 1000 bar
+# Up to 600 degrees C at 1000 bar
Tc <- c(Tc, 500, 550, 600)
w <- water(c("rho", "alpha", "daldT", "beta"), T = convert(Tc, "K"), P = 1000)
gfun.1000 <- CHNOSZ:::gfun(w$rho/1000, Tc, rep(1000, length(Tc)), w$alpha, w$daldT, w$beta)
-# up to 1000 degrees C at 4000 bar
+# Up to 1000 degrees C at 4000 bar
Tc <- c(Tc, 700, 800, 900, 1000)
w <- water(c("rho", "alpha", "daldT", "beta"), T = convert(Tc, "K"), P = 4000)
gfun.4000 <- CHNOSZ:::gfun(w$rho/1000, Tc, rep(4000, length(Tc)), w$alpha, w$daldT, w$beta)
-# values from table 5 of Shock et al., 1992
+# Values from table 5 of Shock et al., 1992
g.Psat.ref <- c(0, 0, 0, 0, -0.09, -1.40, -8.05, -35.23, -192.05)
g.500.ref <- c(0, 0, 0, 0, -0.02, -0.43, -3.55, -16.81, -56.78, -287.16, -1079.75)
g.1000.ref <- c(0, 0, 0, 0, 0, -0.12, -1.49, -8.44, -30.59, -85.09, -201.70, -427.17, -803.53, -1312.67)
g.4000.ref <- c(0, 0, 0, 0, 0, 0, 0, -0.01, -0.20, -1.14, -3.52, -7.84, -14.19, -22.21, -40.15, -54.28, -59.35, -55.17)
-# compare the calculations with the reference values
-# note: tolerance is set as low as possible (order of magnitude) for successful tests
+# Compare the calculations with the reference values
+# Note: tolerance is set as low as possible (order of magnitude) for successful tests
expect_equal(gfun.Psat$g * 1e4, g.Psat.ref, tolerance = 1e-4, info = info)
expect_equal(gfun.500$g * 1e4, g.500.ref, tolerance = 1e-4, info = info)
expect_equal(gfun.1000$g * 1e4, g.1000.ref, tolerance = 1e-5, info = info)
expect_equal(gfun.4000$g * 1e4, g.4000.ref, tolerance = 1e-4, info = info)
-# values from table 6 of Shock et al., 1992
+# Values from table 6 of Shock et al., 1992
dgdT.Psat.ref <- c(0, 0, 0, -0.01, -0.62, -6.13, -25.33, -123.20, -1230.59)
dgdT.500.ref <- c(0, 0, 0, 0, -0.13, -2.21, -12.54, -46.50, -110.45, -734.92, -2691.41)
dgdT.1000.ref <- c(0, 0, 0, 0, -0.02, -0.75, -6.12, -24.83, -69.29, -158.38, -324.13, -594.43, -904.41, -1107.84)
@@ -57,7 +56,7 @@
expect_equal(gfun.1000$dgdT * 1e6, dgdT.1000.ref, tolerance = 1e-5, info = info)
expect_equal(gfun.4000$dgdT * 1e6, dgdT.4000.ref, tolerance = 1e-3, info = info)
-# values from table 7 of Shock et al., 1992
+# Values from table 7 of Shock et al., 1992
d2gdT2.Psat.ref <- c(0, 0, 0, -0.1, -3.7, -20.7, -74.3, -527.8, -15210.4)
d2gdT2.500.ref <- c(0, 0, 0, 0, -1.0, -9.3, -36.8, -109.4, -26.3, -2043.5, -4063.1)
d2gdT2.1000.ref <- c(0, 0, 0, 0, -0.2, -4.0, -20.5, -58.4, -125.4, -241.8, -433.9, -628.0, -550.0, -265.2)
@@ -67,7 +66,7 @@
expect_equal(gfun.1000$d2gdT2 * 1e8, d2gdT2.1000.ref, tolerance = 1e-4, info = info)
expect_equal(gfun.4000$d2gdT2 * 1e8, d2gdT2.4000.ref, tolerance = 1e-2, info = info)
-# values from table 8 of Shock et al., 1992
+# Values from table 8 of Shock et al., 1992
dgdP.Psat.ref <- c(0, 0, 0, 0, 0.03, 0.36, 1.92, 12.27, 227.21)
dgdP.500.ref <- c(0, 0, 0, 0, 0.01, 0.53, 1.75, 5.69, 107.61, 704.88, 1110.20)
dgdP.1000.ref <- c(0, 0, 0, 0, 0, 0.03, 0.31, 1.50, 5.25, 15.53, 4.78, 101.39, 202.25, 313.90)
@@ -77,8 +76,31 @@
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)
-# reference
+# load SiO2 and Si2O4 data taken from DEW spreadsheet
+iSi <- add.OBIGT("DEW", c("SiO2", "Si2O4"))
+Vn1 <- Vn2 <- numeric()
+species <- c("CO3-2", "BO2-", "MgCl+", "SiO2", "HCO3-", "Si2O4")
+# First method: use hkf() function
+for(i in 1:length(species)) {
+ # Get HKF parameters
+ par <- info(info(species[i]), check.it = FALSE)
+ # Get the nonsolvation volume from hkf()
+ Vn1 <- c(Vn1, CHNOSZ:::hkf("V", par, contrib="n")$aq[[1]]$V)
+}
+# In version 2.0.0, hkf() assumes the parameters are Joules,
+# but we gave it calorie-based parameters, so we need to convert to Joules
+Vn1 <- convert(Vn1, "J")
+# Second method: subcrt() with omega = 0
+for(i in 1:length(species)) {
+ mod.OBIGT(species[i], omega = 0)
+ Vn2 <- c(Vn2, subcrt(species[i], T = 25)$out[[1]]$V)
+}
+expect_equal(Vn1, Vn2)
+
+
+# Reference
+
# 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 degrees C and 5 kbar.
Modified: pkg/CHNOSZ/man/CHNOSZ-package.Rd
===================================================================
--- pkg/CHNOSZ/man/CHNOSZ-package.Rd 2022-03-25 11:13:38 UTC (rev 714)
+++ pkg/CHNOSZ/man/CHNOSZ-package.Rd 2022-03-26 08:11:30 UTC (rev 715)
@@ -22,7 +22,7 @@
\item Main workflow: \code{\link{info}}, \code{\link{subcrt}}, \code{\link{basis}}, \code{\link{species}}, \code{\link{affinity}}, \code{\link{equilibrate}}, \code{\link{diagram}}
\item Extended workflow: \code{\link{swap.basis}}, \code{\link{buffer}}, \code{\link{mosaic}}, \code{\link{objective}}, \code{\link{revisit}}, \code{\link{findit}}, \code{\link{EOSregress}}
\item Thermodynamic data: \code{\link{data}}, \code{\link{extdata}}, \code{\link{add.OBIGT}}, \code{\link{util.data}}
- \item Thermodynamic calculations: \code{\link{util.formula}}, \code{\link{makeup}}, \code{\link{util.units}}, \code{\link{eos}}, \code{\link{Berman}}, \code{\link{nonideal}}, \code{\link{util.misc}}
+ \item Thermodynamic calculations: \code{\link{util.formula}}, \code{\link{makeup}}, \code{\link{util.units}}, \code{\link{Berman}}, \code{\link{nonideal}}, \code{\link{util.misc}}
\item Water properties: \code{\link{water}}, \code{\link{util.water}}, \code{\link{DEW}}, \code{\link{IAPWS95}}
\item Protein properties: \code{\link{protein.info}}, \code{\link{add.protein}}, \code{\link{util.fasta}}, \code{\link{util.protein}}, \code{\link{util.seq}}, \code{\link{ionize.aa}}
\item Other tools: \code{\link{examples}}, \code{\link{taxonomy}}
Deleted: pkg/CHNOSZ/man/eos.Rd
===================================================================
--- pkg/CHNOSZ/man/eos.Rd 2022-03-25 11:13:38 UTC (rev 714)
+++ pkg/CHNOSZ/man/eos.Rd 2022-03-26 08:11:30 UTC (rev 715)
@@ -1,116 +0,0 @@
-\encoding{UTF-8}
-\name{eos}
-\alias{eos}
-\alias{hkf}
-\alias{cgl}
-\alias{AD}
-\title{Equations of State}
-\description{
-Calculate thermodynamic properties using the revised Helgeson-Kirkham-Flowers (HKF) or Akinfiev-Diamond (AD) equations of state for aqueous species, or using a generic heat capacity equation for crystalline, gas, and liquid species.
-}
-
-\usage{
- cgl(property = NULL, parameters = NULL, T = 298.15, P = 1)
- hkf(property = NULL, parameters = NULL, T = 298.15, P = 1,
- contrib = c("n", "s", "o"), H2O.props = "rho")
- AD(property = NULL, parameters = NULL, T = 298.15, P = 1, isPsat = TRUE)
-}
-
-\arguments{
- \item{property}{character, name(s) of properties to calculate}
- \item{parameters}{dataframe, species parameters as one or more rows from \code{thermo()$OBIGT}}
- \item{T}{numeric, temperature(s) at which to calculate properties (K)}
- \item{P}{numeric, pressure(s) at which to calculate properties (bar)}
- \item{contrib}{character, which contributions to consider in the revised HKF equations equations of state: (\code{n})onsolvation, (\code{s})olvation (the \eqn{\omega}{omega} terms), or (o)rigination contributions (i.e., the property itself at 25 \degC and 1 bar). Default is \code{c("n","s","o")}, for all contributions}
- \item{H2O.props}{character, properties to calculate for water}
- \item{isPsat}{logical, is this a calculation along the liquid-vapor saturation curve (Psat)?}
-}
-
-\details{
-The equations of state permit the calculation of the standard molal properties of species as a function of temperature and pressure.
-The \code{property} argument is required and refers to one or more of \samp{G}, \samp{H}, \samp{S}, \samp{Cp} and \samp{V}, and for aqueous species only, \samp{kT} and \samp{E}.
-The units of these properties are the first ones shown in the description for \code{\link{subcrt}}.
-The names of the properties are matched without regard to case.
-
-\code{hkf} implements the revised HKF equations of state (Helgeson et al., 1981; Tanger and Helgeson, 1988; Shock and Helgeson, 1988).
-The equations-of-state parameters are \code{a1}, \code{a2}, \code{a3}, \code{a4}, \code{c1}, \code{c2}, \code{omega} and \code{Z}; the units of these parameters are as indicated for \code{\link{thermo}$OBIGT}, without the order of magnitude multipliers.
-Note that the equation-of-state parameter \code{Z} (appearing in the \eqn{g}{g}-function for the temperature derivatives of the omega parameter; Shock et al., 1992) is taken from \code{thermo()$OBIGT} and \emph{not} from the \code{\link{makeup}} of the species.
-\code{H2O.props} is an optional argument that lists the properties of water that should be returned; it is used by \code{\link{subcrt}} so that the time-consuming water calculations are only performed once.
-
-The temperature and pressure derivatives of the \code{omega} parameter for charged species (where \code{Z != 0}, but not for the aqueous proton, H+) are calculated using the \eqn{g}{g}- and \eqn{f}{f}-functions described by Shock et al., 1992 and Johnson et al., 1992.
-If the IAPWS-95 or DEW equations are activated (see \code{\link{water}}), only the \eqn{g}{g}-function (applicable to \samp{G}), but not its derivatives (needed for \samp{H}, \samp{S}, \samp{Cp}, and \samp{V}), is calculated.
-
-The parameters in the \code{cgl} equations of state for crystalline, gas and liquid species (except liquid water) include \code{V}, \code{a}, \code{b}, \code{c}, \code{d}, \code{e}, \code{f} and \code{lambda}.
-The terms denoted by \code{a}, \code{b} and \code{c} correspond to the Maier-Kelley equation for heat capacity (Maier and Kelley, 1932); the additional terms are useful for representing heat capacities of minerals (Robie and Hemingway, 1995) and gaseous or liquid organic species (Helgeson et al., 1998).
-The standard molal volumes (\samp{V}) of species in these calculations are taken to be independent of temperature and pressure.
-
-For both \code{hkf} and \code{cgl}, if at least one equations-of-state parameter for a species is provided, any NA values of the other parameters are reset to zero.
-If all equations-of-state parameters are NA, but values of \samp{Cp} and/or \samp{V} are available, those values are used in the integration of \samp{G}, \samp{H} and \samp{S} as a function of temperature.
-
-\code{AD} implements the Akinfiev-Diamond model for aqueous species (Akinfiev and Diamond, 2003).
-This model is activated by setting \code{abbrv = "AD"} 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).
-
-}
-
-\section{Warning}{
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/chnosz -r 715
More information about the CHNOSZ-commits
mailing list