[CHNOSZ-commits] r400 - in pkg/CHNOSZ: . R inst man tests/testthat
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Feb 19 23:57:50 CET 2019
Author: jedick
Date: 2019-02-19 23:57:50 +0100 (Tue, 19 Feb 2019)
New Revision: 400
Added:
pkg/CHNOSZ/R/AkDi.R
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/NAMESPACE
pkg/CHNOSZ/R/hkf.R
pkg/CHNOSZ/R/info.R
pkg/CHNOSZ/R/subcrt.R
pkg/CHNOSZ/R/util.data.R
pkg/CHNOSZ/inst/NEWS
pkg/CHNOSZ/man/eos.Rd
pkg/CHNOSZ/man/thermo.Rd
pkg/CHNOSZ/tests/testthat/test-eos.R
pkg/CHNOSZ/tests/testthat/test-info.R
pkg/CHNOSZ/tests/testthat/test-thermo.R
Log:
add AkDi(): new function implementing Akinfiev-Diamond model
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2019-02-19 11:12:12 UTC (rev 399)
+++ pkg/CHNOSZ/DESCRIPTION 2019-02-19 22:57:50 UTC (rev 400)
@@ -1,6 +1,6 @@
-Date: 2019-02-19
+Date: 2019-02-20
Package: CHNOSZ
-Version: 1.2.0-6
+Version: 1.2.0-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 2019-02-19 11:12:12 UTC (rev 399)
+++ pkg/CHNOSZ/NAMESPACE 2019-02-19 22:57:50 UTC (rev 400)
@@ -58,7 +58,7 @@
# added 20171121 or later
"dumpdata", "thermo.axis", "solubility", "NaCl",
# added 20190213 or later
- "CHNOSZ", "thermo", "reset", "obigt", "retrieve"
+ "CHNOSZ", "thermo", "reset", "obigt", "retrieve", "AkDi"
)
# Load shared objects
Added: pkg/CHNOSZ/R/AkDi.R
===================================================================
--- pkg/CHNOSZ/R/AkDi.R (rev 0)
+++ pkg/CHNOSZ/R/AkDi.R 2019-02-19 22:57:50 UTC (rev 400)
@@ -0,0 +1,54 @@
+# CHNOSZ/AkDi.R
+# Akinfiev-Diamond model for aqueous species
+# 20190219 first version
+
+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
+
+ # a list for the output
+ out <- list()
+ # loop over species
+ nspecies <- nrow(parameters)
+ for(i in seq_len(nspecies)) {
+ PAR <- parameters[i, ]
+ # start with an NA-filled data frame
+ myprops <- as.data.frame(matrix(NA, ncol=length(property), nrow=length(T)))
+ colnames(myprops) <- property
+ # just calculate G for now
+ for(j in seq_along(property)) {
+ 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) )
+ # 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
+ }
+ }
+ out[[i]] <- myprops
+ }
+ out
+}
Modified: pkg/CHNOSZ/R/hkf.R
===================================================================
--- pkg/CHNOSZ/R/hkf.R 2019-02-19 11:12:12 UTC (rev 399)
+++ pkg/CHNOSZ/R/hkf.R 2019-02-19 22:57:50 UTC (rev 400)
@@ -48,7 +48,7 @@
# a list to store the result
aq.out <- list()
nspecies <- nrow(parameters)
- for(k in 1:nspecies) {
+ for(k in seq_len(nspecies)) {
# loop over each species
PAR <- parameters[k, ]
# substitute Cp and V for missing EoS parameters
@@ -71,7 +71,7 @@
dwdP <- dwdT <- d2wdT2 <- numeric(length(T))
Z <- PAR$Z
omega.PT <- rep(PAR$omega, length(T))
- if(!identical(Z, 0) & !identical(PAR$name, "H+")) {
+ if(!identical(Z, 0) & !is.na(Z) & !identical(PAR$name, "H+")) {
# compute derivatives of omega: g and f functions (Shock et al., 1992; Johnson et al., 1992)
rhohat <- H2O.PT$rho/1000 # just converting kg/m3 to g/cm3
g <- gfun(rhohat, convert(T, "C"), P, H2O.PT$alpha, H2O.PT$daldT, H2O.PT$beta)
Modified: pkg/CHNOSZ/R/info.R
===================================================================
--- pkg/CHNOSZ/R/info.R 2019-02-19 11:12:12 UTC (rev 399)
+++ pkg/CHNOSZ/R/info.R 2019-02-19 22:57:50 UTC (rev 400)
@@ -194,7 +194,8 @@
this[, which(naGHS)+7] <- GHS[naGHS]
}
# now perform consistency checks for GHS and EOS parameters if check.it=TRUE
- if(check.it) {
+ # don't do it for the AkDi species 20190219
+ if(check.it & !"xi" %in% colnames(this)) {
# check GHS if they were all present
if(sum(naGHS)==0) calcG <- checkGHS(this)
# check tabulated heat capacities against EOS parameters
Modified: pkg/CHNOSZ/R/subcrt.R
===================================================================
--- pkg/CHNOSZ/R/subcrt.R 2019-02-19 11:12:12 UTC (rev 399)
+++ pkg/CHNOSZ/R/subcrt.R 2019-02-19 22:57:50 UTC (rev 400)
@@ -9,6 +9,7 @@
#source("util.units.R")
#source("util.data.R")
#source("species.R")
+#source("AkDi.R")
subcrt <- function(species, coeff = 1, state = NULL, property = c("logK", "G", "H", "S", "V", "Cp"),
T = seq(273.15, 623.15, 25), P = "Psat", grid = NULL, convert = TRUE, exceed.Ttr = FALSE,
@@ -75,6 +76,7 @@
}
# gridding?
+ isPsat <- FALSE
do.grid <- FALSE
if(!is.null(grid)) if(!is.logical(grid)) do.grid <- TRUE
newIS <- IS
@@ -101,6 +103,8 @@
P <- rep(tpargs$P,length.out=length(newIS))
}
} else {
+ # for AkDi, 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)
T <- tpargs$T; P <- tpargs$P
@@ -277,17 +281,22 @@
if(TRUE %in% isaq) {
# 20110808 get species parameters using obigt2eos() (faster than using info())
param <- obigt2eos(thermo$obigt[iphases[isaq],], "aq", fixGHS = TRUE)
+ # aqueous species with NA for Z use the AkDi model
+ isAkDi <- is.na(param$Z)
# always get density
H2O.props <- "rho"
# calculate A_DH and B_DH if we're using the B-dot (Helgeson) equation
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)
- hkfstuff <- hkf(eosprop, parameters = param, T = T, P = P, H2O.props=H2O.props)
+ # in case everything is AkDi, run hkf (for water properties) but exclude all species
+ hkfpar <- param
+ if(all(isAkDi)) 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) {
+ if(!exceed.rhomin & !all(isAkDi)) {
ilowrho <- H2O.PT$rho < 350
ilowrho[is.na(ilowrho)] <- FALSE
if(any(ilowrho)) {
@@ -296,6 +305,12 @@
warnings <- c(warnings, paste0("below minimum density for applicability of revised HKF equations (", sum(ilowrho), " T,P ", ptext, ")"))
}
}
+ # calculate properties using Akinfiev-Diamond model 20190219
+ if(any(isAkDi)) {
+ # get the parameters with the right names
+ param <- obigt2eos(param[isAkDi, ], "aq")
+ p.aq[isAkDi] <- AkDi(eosprop, parameters = param, T = T, P = P, isPsat = isPsat)
+ }
# calculate activity coefficients if ionic strength is not zero
if(any(IS != 0)) {
if(thermo$opt$nonideal %in% c("Bdot", "Bdot0", "bgamma", "bgamma0")) p.aq <- nonideal(iphases[isaq], p.aq, newIS, T, P, H2O.PT$A_DH, H2O.PT$B_DH)
Modified: pkg/CHNOSZ/R/util.data.R
===================================================================
--- pkg/CHNOSZ/R/util.data.R 2019-02-19 11:12:12 UTC (rev 399)
+++ pkg/CHNOSZ/R/util.data.R 2019-02-19 22:57:50 UTC (rev 400)
@@ -397,8 +397,15 @@
# remove scaling factors from EOS parameters
# and apply column names depending on the EOS
if(identical(state, "aq")) {
- obigt[,13:20] <- t(t(obigt[,13:20]) * 10^c(-1,2,0,4,0,4,5,0))
- colnames(obigt)[13:20] <- c('a1','a2','a3','a4','c1','c2','omega','Z')
+ # species in the Akinfiev-Diamond model (AkDi) have NA for Z 20190219
+ isAkDi <- is.na(obigt$z.T)
+ # remove scaling factors for the HKF species, but not for the AkDi species
+ obigt[!isAkDi, 13:20] <- t(t(obigt[!isAkDi, 13:20]) * 10^c(-1,2,0,4,0,4,5,0))
+ # for AkDi specie, set NA values in remaining columns (for display only)
+ obigt[isAkDi, 16:19] <- NA
+ # if all of the species are AkDi, change the variable names
+ if(all(isAkDi)) colnames(obigt)[13:20] <- c('a','b','xi','XX1','XX2','XX3','XX4','Z')
+ else colnames(obigt)[13:20] <- c('a1','a2','a3','a4','c1','c2','omega','Z')
} else {
obigt[,13:20] <- t(t(obigt[,13:20]) * 10^c(0,-3,5,0,-5,0,0,0))
colnames(obigt)[13:20] <- c('a','b','c','d','e','f','lambda','T')
Modified: pkg/CHNOSZ/inst/NEWS
===================================================================
--- pkg/CHNOSZ/inst/NEWS 2019-02-19 11:12:12 UTC (rev 399)
+++ pkg/CHNOSZ/inst/NEWS 2019-02-19 22:57:50 UTC (rev 400)
@@ -1,4 +1,4 @@
-CHANGES IN CHNOSZ 1.2.0-6 (2019-02-19)
+CHANGES IN CHNOSZ 1.2.0-8 (2019-02-20)
--------------------------------------
CRAN COMPLIANCE
@@ -18,9 +18,9 @@
- Add retrieve() to retrieve all the species with given elements. Thanks
to Evgeniy Bastrakov for the suggestion.
-- water() and water.SUPCRT92(): add 'P1' argument to choose whether to
- output 1 bar for Psat at temperatures less than 100 degrees C
- (default is TRUE).
+- Add AkDi() to calculate thermodynamic properties of aqueous
+ nonelectrolytes using the Akinfiev-Diamond model. Thanks to Evgeniy
+ Bastrakov for guidance and advice.
THERMODYNAMIC DATA
@@ -30,18 +30,23 @@
- Move SiO2(aq) from Apps and Spycher, 2004 and recalculated HSiO3- to
new optional data file, OBIGT/AS04.csv.
-- Move H2AsO3- from SLOP98.csv to SUPCRT92.csv.
+- Move H2AsO3- from OBIGT/SLOP98.csv to OBIGT/SUPCRT92.csv.
DOCUMENTATION
- In demo/NaCl.R, indicate region not considered by Shock et al., 1992
- in developing the revised HKF model and presence of the associated
- discontinuities in SUPCRT92.
+ in developing the "g function" applicable to electrolytes in the
+ revised HKF model, and note presence of resulting discontinuities (see
+ man/examples.Rd).
OTHER CHANGES
- In describe.reaction(), change equals sign to reaction double arrow.
+- water() and water.SUPCRT92(): add 'P1' argument to choose whether to
+ output 1 bar for Psat at temperatures less than 100 degrees C
+ (default is TRUE).
+
CHANGES IN CHNOSZ 1.2.0 (2019-02-09)
------------------------------------
Modified: pkg/CHNOSZ/man/eos.Rd
===================================================================
--- pkg/CHNOSZ/man/eos.Rd 2019-02-19 11:12:12 UTC (rev 399)
+++ pkg/CHNOSZ/man/eos.Rd 2019-02-19 22:57:50 UTC (rev 400)
@@ -3,15 +3,17 @@
\alias{eos}
\alias{hkf}
\alias{cgl}
+\alias{AkDi}
\title{Equations of State}
\description{
-Calculate thermodynamic properties using the revised Helgeson-Kirkham-Flowers (HKF) equations of state for aqueous species, or using a generic heat capacity equation for crystalline, gas, and liquid species.
+Calculate thermodynamic properties using the revised Helgeson-Kirkham-Flowers (HKF) or Akinfiev-Diamond (AkDi) 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")
+ AkDi(property = NULL, parameters = NULL, T = 298.15, P = 1, isPsat = TRUE)
}
\arguments{
@@ -21,6 +23,7 @@
\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{
@@ -43,11 +46,16 @@
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{AkDi} provides the Akinfiev-Diamond model for aqueous species (Akinfiev and Diamond, 2003).
+To run this code, the database must also include the corresponding gasesous species (with the same name or chemical formula).
+
}
\section{Warning}{
-The temperature and pressure range of validity of the revised HKF equations of state for aqueous species corresponds to the stability region of liquid water or the supercritical fluid at conditions between 0 to 1000 \degC and 1 to 5000 bar (Tanger and Helgeson, 1988; Shock and Helgeson, 1988).
-The \code{hkf} function does not check these limits and will compute properties as long as the requisite electrostatic properties of water are available. There are conceptually no temperature limits (other than 0 Kelvin) for the validity of the \code{cgl} equations of state.
+The range of applicability of the revised HKF equations of state for aqueous species corresponds to the stability region of liquid water or the supercritical fluid with density greater than 0.35 g/cm3, and between 0 to 1000 \degC and 1 to 5000 bar (Tanger and Helgeson, 1988; Shock and Helgeson, 1988).
+The \code{hkf} function does not check these limits and will compute properties as long as the requisite electrostatic properties of water are available.
+There are conceptually no temperature limits (other than 0 Kelvin) for the validity of the \code{cgl} equations of state.
However, the actual working upper temperature limits correspond to the temperatures of phase transitions of minerals or to those temperatures beyond which extrapolations from experimental data become highly uncertain.
These temperature limits are stored in the thermodynamic database for some minerals, but \code{cgl} ignores them; however, \code{\link{subcrt}} warns if they are exceeded.
}
@@ -86,6 +94,8 @@
\references{
+ Akinfiev, N. N. and Diamond, L. W. (2003) Thermodynamic description of aqueous nonelectrolytes at infinite dilution over a wide range of state parameters. \emph{Geochim. Cosmochim. Acta} \bold{67}, 613--629. \url{https://doi.org/10.1016/S0016-7037(02)01141-9}
+
Helgeson, H. C., Kirkham, D. H. and Flowers, G. C. (1981) Theoretical prediction of the thermodynamic behavior of aqueous electrolytes at high pressures and temperatures. IV. Calculation of activity coefficients, osmotic coefficients, and apparent molal and standard and relative partial molal properties to 600\degC and 5 Kb. \emph{Am. J. Sci.} \bold{281}, 1249--1516. \url{https://doi.org/10.2475/ajs.281.10.1249}
Helgeson, H. C., Owens, C. E., Knox, A. M. and Richard, L. (1998) Calculation of the standard molal thermodynamic properties of crystalline, liquid, and gas organic molecules at high temperatures and pressures. \emph{Geochim. Cosmochim. Acta} \bold{62}, 985--1081. \url{https://doi.org/10.1016/S0016-7037(97)00219-6}
Modified: pkg/CHNOSZ/man/thermo.Rd
===================================================================
--- pkg/CHNOSZ/man/thermo.Rd 2019-02-19 11:12:12 UTC (rev 399)
+++ pkg/CHNOSZ/man/thermo.Rd 2019-02-19 22:57:50 UTC (rev 400)
@@ -83,19 +83,21 @@
Note the following database conventions:
\itemize{
\item The combination of \code{name} and \code{state} defines a species in \code{thermo$obigt}. A species can not be duplicated (this is checked when running \code{reset()}).
- \item English names of gases are used only for the gas state. The dissolved species is named with the chemical formula. Therefore, \code{info("oxygen")} refers to the gas, and \code{info("O2")} refers to the aqueous species.
+ \item English names of inorganic gases are used only for the gas state. The dissolved species is named with the chemical formula. Therefore, \code{info("oxygen")} refers to the gas, and \code{info("O2")} refers to the aqueous species.
+ \item Properties of most aqueous species (\code{state} = \samp{aq}) are calculated using the revised Helgeson-Kirkham-Flowers (HKF) model (see \code{\link{hkf}}).
+ \item Properties of aqueous species with an NA value of \code{Z} (the final column of thermo$obigt) are calculated using the Akinfiev-Diamond model (see \code{\link{AkDi}}).
+ \item Properties of most non-aqueous species (liquids, gases, and minerals) are calculated using a heat capacity polynomial expression with up to six terms (see \code{\link{cgl}}).
+ \item Properties of minerals with NA values of all heat capacity parameters are calculated using the Berman model (see \code{\link{berman}}).
}
\samp{OrganoBioGeoTherm} is the name of a GUI program to use SUPCRT in Windows, produced in Harold C. Helgeson's Laboratory of Theoretical Geochemistry and Biogeochemistry at the University of California, Berkeley.
- The \acronym{OBIGT} database was originally developed for that program, and has been ported to CHNOSZ, with additional modifications.
+ The \acronym{OBIGT} database was originally developed for that program, and was the original basis for the database in CHNOSZ.
There may be an additional meaning for the acronym: \dQuote{One BIG Table} of thermodynamic data.
Each entry is referenced to one or two literature sources listed in \code{thermo$refs}.
Use \code{\link{thermo.refs}} to look up the citation information for the references.
- \acronym{OBIGT} was initially built from the \acronym{SUPCRT92} (Johnson et al., 1992) and \acronym{slop98.dat} data files (Shock et al., 1998), and the references in those files are included here.
- Some data in \acronym{slop98.dat} were corrected or modified as noted in that file; these modifications are indicated in OBIGT by having \samp{SLOP98} as one of the sources of data.
- Other additions or modifications used in CHNOSZ are indicated by having \samp{CHNOSZ} as one of the sources of data.
See the vignette \emph{Thermodynamic data in CHNOSZ} for a complete description of the sources of data.
+ The original \acronym{OBIGT} database was influenced by the \acronym{SUPCRT92} (Johnson et al., 1992) and \acronym{slop98.dat} data files (Shock et al., 1998), and the references in those files are included here.
In order to represent thermodynamic data for minerals with phase transitions, the higher-temperature phases of these minerals are represented as phase species that have states denoted by \samp{cr2}, \samp{cr3}, etc.
The standard molar thermodynamic properties at 25 \degC and 1 bar (\eqn{T_r}{Pr} and \eqn{P_r}{Pr}) of the \samp{cr2} phase species of minerals were generated by first calculating those of the \samp{cr} (lowest-T) phase species at the transition temperature (\eqn{T_{tr}}{Ttr}) and 1 bar then taking account of the volume and entropy of transition (the latter can be retrieved by combining the former with the Clausius-Clapeyron equation and values of \eqn{(dP/dT)} of transitions taken from the \acronym{SUPCRT92} data file) to calculate the standard molar entropy of the \samp{cr2} phase species at \eqn{T_{tr}}{Ttr}, and taking account of the enthalpy of transition (\eqn{{\Delta}H^{\circ}}{DeltaH0}, taken from the \acronym{SUPCRT92} data file) to calculate the standard molar enthalpy of the \samp{cr2} phase species at \eqn{T_{tr}}{Ttr}.
@@ -124,9 +126,8 @@
}
-The meanings of the remaining columns depend on the physical state of a particular species.
-If it is aqueous, the values in these columns represent parameters in the revised HKF equations of state (see \code{\link{hkf}}), otherwise they denote parameters in a general equations for crystalline, gas and liquid species (see \code{\link{cgl}}).
-The names of these columns are compounded from those of the parameters in each of the equations of state (for example, column 13 is named \code{a1.a}).
+The meanings of the remaining columns depend on the model used for a particular species (see database conventions above).
+The names of these columns are compounded from those of the parameters in the HKF equations of state and general heat capacity polynomial; for example, column 13 is named \code{a1.a}.
Scaling of the values by orders of magnitude is adopted for some of the parameters, following common usage in the literature.
Columns 13-20 for aqueous species (parameters in the revised HKF equations of state):
@@ -156,6 +157,19 @@
\tab \tab temperature limit of validity of extrapolation (K)
}
+Columns 13-20 for aqueous species using the Akinfiev-Diamond model. Note that the \code{c} column is used to store the \eqn{\xi}{xi} parameter, and that \code{Z} must be NA to activate the code for this model.
+
+ \tabular{lll}{
+ \code{a} \tab numeric \tab \eqn{a} \cr
+ \code{b} \tab numeric \tab \eqn{b} \cr
+ \code{c} \tab numeric \tab \eqn{\xi}{xi} \cr
+ \code{d} \tab numeric \tab \eqn{XX1} NA \cr
+ \code{e} \tab numeric \tab \eqn{XX2} NA \cr
+ \code{f} \tab numeric \tab \eqn{XX3} NA \cr
+ \code{lambda} \tab numeric \tab \eqn{XX4} NA \cr
+ \code{Z} \tab numeric \tab \eqn{Z} NA \cr
+ }
+
\item \code{thermo$refs}
Dataframe of references to sources of thermodynamic data.
\tabular{lll}{
Modified: pkg/CHNOSZ/tests/testthat/test-eos.R
===================================================================
--- pkg/CHNOSZ/tests/testthat/test-eos.R 2019-02-19 11:12:12 UTC (rev 399)
+++ pkg/CHNOSZ/tests/testthat/test-eos.R 2019-02-19 22:57:50 UTC (rev 400)
@@ -78,6 +78,30 @@
expect_equal(gfun.4000$dgdP * 1e6, dgdP.4000.ref, tolerance=1e-3)
})
+test_that("AkDi produces expected results", {
+ # 20190220
+ # add an aqueous species conforming to the AkDi model: it has NA for Z
+ iCO2 <- mod.obigt("CO2", a=-8.8321, b=11.2684, c=-0.0850, z=NA)
+ # 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)
+
+ 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)
+
+ reset()
+})
+
# reference
# Shock, E. L., Oelkers, E. H., Johnson, J. W., Sverjensky, D. A. and Helgeson, H. C. (1992)
Modified: pkg/CHNOSZ/tests/testthat/test-info.R
===================================================================
--- pkg/CHNOSZ/tests/testthat/test-info.R 2019-02-19 11:12:12 UTC (rev 399)
+++ pkg/CHNOSZ/tests/testthat/test-info.R 2019-02-19 22:57:50 UTC (rev 400)
@@ -33,3 +33,12 @@
expect_equal(thermo()$obigt$state[i2], c("cr", "aq"))
expect_equal(info(i2)[1, ], info(i2[1]), check.attributes=FALSE)
})
+
+test_that("info() gives correct column names for species using the AkDi model", {
+ # add an aqueous species conforming to the AkDi model: it has NA for Z
+ iCO2 <- mod.obigt("CO2", a = -8.8321, b = 11.2684, c = -0.0850, z = NA)
+ params <- info(iCO2)
+ expect_equal(params$a, -8.8321)
+ expect_equal(params$b, 11.2684)
+ expect_equal(params$xi, -0.0850)
+})
Modified: pkg/CHNOSZ/tests/testthat/test-thermo.R
===================================================================
--- pkg/CHNOSZ/tests/testthat/test-thermo.R 2019-02-19 11:12:12 UTC (rev 399)
+++ pkg/CHNOSZ/tests/testthat/test-thermo.R 2019-02-19 22:57:50 UTC (rev 400)
@@ -30,6 +30,8 @@
"please supply a valid chemical formula")
# the default state is aq
expect_message(itest <- mod.obigt("test", formula="Z0", date=today()), "added test\\(aq\\)")
+ # set the charge so following test use hkf() rather than AkDi()
+ mod.obigt("test", z = 0)
# we should get NA values of G for a species with NA properties
expect_true(all(is.na(subcrt(itest)$out[[1]]$G)))
# a single value of G comes through to subcrt
More information about the CHNOSZ-commits
mailing list