[CHNOSZ-commits] r205 - in pkg/CHNOSZ: . R data demo inst inst/extdata/OBIGT man tests/testthat vignettes
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Jul 27 02:01:57 CEST 2017
Author: jedick
Date: 2017-07-27 02:01:56 +0200 (Thu, 27 Jul 2017)
New Revision: 205
Added:
pkg/CHNOSZ/R/add.obigt.R
pkg/CHNOSZ/demo/adenine.R
pkg/CHNOSZ/inst/extdata/OBIGT/CHNOSZ_aq.csv
pkg/CHNOSZ/man/add.obigt.Rd
Removed:
pkg/CHNOSZ/inst/extdata/OBIGT/CHNOSZ_aq.csv.xz
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/R/examples.R
pkg/CHNOSZ/R/transfer.R
pkg/CHNOSZ/R/util.data.R
pkg/CHNOSZ/data/refs.csv
pkg/CHNOSZ/data/thermo.R
pkg/CHNOSZ/demo/00Index
pkg/CHNOSZ/demo/activity_ratios.R
pkg/CHNOSZ/inst/NEWS
pkg/CHNOSZ/man/EOSregress.Rd
pkg/CHNOSZ/man/examples.Rd
pkg/CHNOSZ/man/transfer.Rd
pkg/CHNOSZ/man/util.data.Rd
pkg/CHNOSZ/tests/testthat/test-util.data.R
pkg/CHNOSZ/vignettes/anintro.Rmd
pkg/CHNOSZ/vignettes/eos-regress.Rmd
pkg/CHNOSZ/vignettes/obigt.Rmd
pkg/CHNOSZ/vignettes/obigt.bib
pkg/CHNOSZ/vignettes/vig.bib
Log:
add adenine (Lowe et al., 2017) to CHNOSZ_aq.csv
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2017-06-06 14:35:05 UTC (rev 204)
+++ pkg/CHNOSZ/DESCRIPTION 2017-07-27 00:01:56 UTC (rev 205)
@@ -1,14 +1,14 @@
-Date: 2017-06-06
+Date: 2017-07-27
Package: CHNOSZ
-Version: 1.1.0-3
+Version: 1.1.0-4
Title: Chemical Thermodynamics and Activity Diagrams
Author: Jeffrey Dick
Maintainer: Jeffrey Dick <j3ffdick at gmail.com>
Depends: R (>= 3.1.0)
Suggests: limSolve, testthat, knitr, rmarkdown, tufte, RSVGTipsDevice
Imports: grDevices, graphics, stats, utils, colorspace
-Description: An integrated set of tools for thermodynamic calculations in compositional
- biology and geochemistry. Thermodynamic properties are taken from a database for minerals
+Description: An integrated set of tools for thermodynamic calculations in geochemistry and
+ compositional biology. Thermodynamic properties are taken from a database for minerals
and inorganic and organic aqueous species including biomolecules, or from amino acid
group additivity for proteins. High-temperature properties are calculated using the
revised Helgeson-Kirkham-Flowers equations of state for aqueous species. Functions are
Added: pkg/CHNOSZ/R/add.obigt.R
===================================================================
--- pkg/CHNOSZ/R/add.obigt.R (rev 0)
+++ pkg/CHNOSZ/R/add.obigt.R 2017-07-27 00:01:56 UTC (rev 205)
@@ -0,0 +1,162 @@
+# CHNOSZ/add.obigt.R
+# add or change entries in the thermodynamic database
+
+today <- function() {
+ # write today's date in the format used in SUPCRT data files
+ # e.g. 13.May.12 for 2012-05-13
+ t <- date()
+ tt <- unlist(strsplit(t, " "))
+ # for single-digit days there is an extra space
+ tt <- tt[!tt==""]
+ tday <- tt[3]
+ tmonth <- tt[2]
+ tyear <- substr(tt[5], start=3, stop=4)
+ return(paste(tday, tmonth, tyear, sep="."))
+}
+
+mod.obigt <- function(...) {
+ # add or modify species in thermo$obigt
+ thermo <- get("thermo")
+ # the names and values are in the arguments
+ # this works for providing arguments via do.call
+ args <- list(...)
+ # this is needed if we are called with a list as the actual argument
+ if(is.list(args[[1]])) args <- args[[1]]
+ if(length(args) < 2) stop("please supply at least a species name and a property to update")
+ if(is.null(names(args))) stop("please supply named arguments")
+ # if the first argument is numeric, it's the species index
+ if(is.numeric(args[[1]][1])) {
+ ispecies <- args[[1]]
+ } else {
+ # if the name of the first argument is missing, assume it's the species name
+ if(names(args)[1]=="") names(args)[1] <- "name"
+ # search for this species, use check.protein=FALSE to avoid infinite loop when adding proteins
+ # and suppressMessages to not show messages about matches of this name to other states
+ if("state" %in% names(args)) ispecies <- suppressMessages(mapply(info.character,
+ species=args$name, state=args$state, check.protein=FALSE, SIMPLIFY=TRUE, USE.NAMES=FALSE))
+ else ispecies <- suppressMessages(mapply(info.character,
+ species=args$name, check.protein=FALSE, SIMPLIFY=TRUE, USE.NAMES=FALSE))
+ }
+ # the column names of thermo$obigt, split at the "."
+ cnames <- c(do.call(rbind, strsplit(colnames(thermo$obigt), ".", fixed=TRUE)), colnames(thermo$obigt))
+ # the columns we are updating
+ icol <- match(names(args), cnames)
+ if(any(is.na(icol))) stop(paste("properties not in thermo$obigt:", paste(names(args)[is.na(icol)], collapse=" ")) )
+ # the column numbers for properties that matched after the split
+ icol[icol > 40] <- icol[icol > 40] - 40
+ icol[icol > 20] <- icol[icol > 20] - 20
+ # which species are new and which are old
+ inew <- which(is.na(ispecies))
+ iold <- which(!is.na(ispecies))
+ # the arguments as data frame
+ args <- data.frame(args, stringsAsFactors=FALSE)
+ if(length(inew) > 0) {
+ # the right number of blank rows of thermo$obigt
+ newrows <- thermo$obigt[1:length(inew), ]
+ # if we don't know something it's NA
+ newrows[] <- NA
+ # put in a default state
+ newrows$state <- thermo$opt$state
+ # the formula defaults to the name
+ newrows$formula <- args$name[inew]
+ # fill in the columns
+ newrows[, icol] <- args[inew, ]
+ # now check the formulas
+ e <- tryCatch(makeup(newrows$formula), error=function(e) e)
+ if(inherits(e, "error")) {
+ warning("please supply a valid chemical formula as the species name or in the 'formula' argument")
+ # transmit the error from makeup
+ stop(e)
+ }
+ # assign to thermo$obigt
+ thermo$obigt <- rbind(thermo$obigt, newrows)
+ rownames(thermo$obigt) <- NULL
+ assign("thermo", thermo, "CHNOSZ")
+ # update ispecies
+ ntotal <- nrow(thermo$obigt)
+ ispecies[inew] <- (ntotal-length(inew)+1):ntotal
+ # inform user
+ message(paste("mod.obigt: added ", newrows$name, "(", newrows$state, ")", sep="", collapse="\n"))
+ }
+ if(length(iold) > 0) {
+ # loop over species
+ for(i in 1:length(iold)) {
+ # the old values and the state
+ oldprop <- thermo$obigt[ispecies[iold[i]], icol]
+ state <- thermo$obigt$state[ispecies[iold[i]]]
+ # tell user if they're the same, otherwise update the data entry
+ if(isTRUE(all.equal(oldprop, args[iold[i], ], check.attributes=FALSE)))
+ message("mod.obigt: no change for ", args$name[iold[i]], "(", state, ")")
+ else {
+ thermo$obigt[ispecies[iold[i]], icol] <- args[iold[i], ]
+ assign("thermo", thermo, "CHNOSZ")
+ message("mod.obigt: updated ", args$name[iold[i]], "(", state, ")")
+ }
+ }
+ }
+ return(ispecies)
+}
+
+add.obigt <- function(file, force=TRUE, E.units="cal") {
+ # add/replace entries in thermo$obigt from values saved in a file
+ # only replace if force==TRUE
+ thermo <- get("thermo")
+ to1 <- thermo$obigt
+ id1 <- paste(to1$name,to1$state)
+ # is the source the name of a species or a file?
+ CHNOSZfile <- system.file("extdata/OBIGT/CHNOSZ_aq.csv", package="CHNOSZ")
+ dat <- read.csv(CHNOSZfile, as.is=TRUE)
+ idat <- match(file, dat$name)
+ newdat <- dat[idat, ]
+ if(!all(is.na(newdat))) to2 <- newdat # a species
+ else to2 <- read.csv(file,as.is=TRUE) # a file
+ id2 <- paste(to2$name,to2$state)
+ # check if the data is compatible with thermo$obigt
+ tr <- try(rbind(to1,to2),silent=TRUE)
+ if(identical(class(tr),'try-error')) stop(paste(file,"is not compatible with thermo$obigt data table."))
+ # match the new species to existing ones
+ does.exist <- id2 %in% id1
+ ispecies.exist <- na.omit(match(id2, id1))
+ nexist <- sum(does.exist)
+ # convert from J if necessary
+ if(tolower(E.units)=="j") {
+ # loop over each row
+ for(i in 1:nrow(to2)) {
+ # GHS and EOS parameters
+ icol <- (8:18)
+ # if it's aqueous, also include omega
+ if(to2$state[i]=="aq") icol <-(8:19)
+ # don't touch volume (in column 12)
+ icol <- icol[icol!=12]
+ # convert to calories
+ to2[i,icol] <- convert(to2[i,icol],"cal")
+ }
+ }
+ # keep track of the species we've added
+ inew <- numeric()
+ if(force) {
+ # replace existing entries
+ if(nexist > 0) {
+ to1[ispecies.exist, ] <- to2[does.exist, ]
+ to2 <- to2[!does.exist, ]
+ inew <- c(inew, ispecies.exist)
+ }
+ } else {
+ # ignore any new entries that already exist
+ to2 <- to2[!does.exist, ]
+ nexist <- 0
+ }
+ # add new entries
+ if(nrow(to2) > 0) {
+ to1 <- rbind(to1, to2)
+ inew <- c(inew, (length(id1)+1):nrow(to1))
+ }
+ # commit the change
+ thermo$obigt <- to1
+ rownames(thermo$obigt) <- 1:nrow(thermo$obigt)
+ assign("thermo", thermo, "CHNOSZ")
+ message("add.obigt: read ", length(does.exist), " rows; made ",
+ nexist, " replacements, ", nrow(to2), " additions, units = ", E.units)
+ message("add.obigt: use data(thermo) to restore default database")
+ return(invisible(inew))
+}
Modified: pkg/CHNOSZ/R/examples.R
===================================================================
--- pkg/CHNOSZ/R/examples.R 2017-06-06 14:35:05 UTC (rev 204)
+++ pkg/CHNOSZ/R/examples.R 2017-07-27 00:01:56 UTC (rev 205)
@@ -28,7 +28,8 @@
demos <- function(which=c("sources", "protein.equil", "affinity", "NaCl", "density",
"ORP", "revisit", "findit", "ionize", "buffer", "protbuff", "yeastgfp", "mosaic",
- "copper", "solubility", "wjd", "dehydration", "bugstab", "Shh", "activity_ratios"), to.file=FALSE) {
+ "copper", "solubility", "wjd", "dehydration", "bugstab", "Shh", "activity_ratios",
+ "adenine"), to.file=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/transfer.R
===================================================================
--- pkg/CHNOSZ/R/transfer.R 2017-06-06 14:35:05 UTC (rev 204)
+++ pkg/CHNOSZ/R/transfer.R 2017-07-27 00:01:56 UTC (rev 205)
@@ -806,6 +806,8 @@
# or feldspar("open")
# setup conditions for feldspar reaction
basis(delete=TRUE)
+ # make pseudo-H4SiO4 available
+ add.obigt("pseudo-H4SiO4")
basis(c('Al+3','pseudo-H4SiO4','K+','H2O','H+','O2'))
# some of SLS94's initial conditions
basis(c('K+','H4SiO4'),c(-6,-6))
Modified: pkg/CHNOSZ/R/util.data.R
===================================================================
--- pkg/CHNOSZ/R/util.data.R 2017-06-06 14:35:05 UTC (rev 204)
+++ pkg/CHNOSZ/R/util.data.R 2017-07-27 00:01:56 UTC (rev 205)
@@ -1,165 +1,6 @@
# CHNOSZ/util.data.R
-# add or change entries in the thermodynamic database
+# check entries in the thermodynamic database
-today <- function() {
- # write today's date in the format used in SUPCRT data files
- # e.g. 13.May.12 for 2012-05-13
- t <- date()
- tt <- unlist(strsplit(t, " "))
- # for single-digit days there is an extra space
- tt <- tt[!tt==""]
- tday <- tt[3]
- tmonth <- tt[2]
- tyear <- substr(tt[5], start=3, stop=4)
- return(paste(tday, tmonth, tyear, sep="."))
-}
-
-mod.obigt <- function(...) {
- # add or modify species in thermo$obigt
- thermo <- get("thermo")
- # the names and values are in the arguments
- # this works for providing arguments via do.call
- args <- list(...)
- # this is needed if we are called with a list as the actual argument
- if(is.list(args[[1]])) args <- args[[1]]
- if(length(args) < 2) stop("please supply at least a species name and a property to update")
- if(is.null(names(args))) stop("please supply named arguments")
- # if the first argument is numeric, it's the species index
- if(is.numeric(args[[1]][1])) {
- ispecies <- args[[1]]
- } else {
- # if the name of the first argument is missing, assume it's the species name
- if(names(args)[1]=="") names(args)[1] <- "name"
- # search for this species, use check.protein=FALSE to avoid infinite loop when adding proteins
- # and suppressMessages to not show messages about matches of this name to other states
- if("state" %in% names(args)) ispecies <- suppressMessages(mapply(info.character,
- species=args$name, state=args$state, check.protein=FALSE, SIMPLIFY=TRUE, USE.NAMES=FALSE))
- else ispecies <- suppressMessages(mapply(info.character,
- species=args$name, check.protein=FALSE, SIMPLIFY=TRUE, USE.NAMES=FALSE))
- }
- # the column names of thermo$obigt, split at the "."
- cnames <- c(do.call(rbind, strsplit(colnames(thermo$obigt), ".", fixed=TRUE)), colnames(thermo$obigt))
- # the columns we are updating
- icol <- match(names(args), cnames)
- if(any(is.na(icol))) stop(paste("properties not in thermo$obigt:", paste(names(args)[is.na(icol)], collapse=" ")) )
- # the column numbers for properties that matched after the split
- icol[icol > 40] <- icol[icol > 40] - 40
- icol[icol > 20] <- icol[icol > 20] - 20
- # which species are new and which are old
- inew <- which(is.na(ispecies))
- iold <- which(!is.na(ispecies))
- # the arguments as data frame
- args <- data.frame(args, stringsAsFactors=FALSE)
- if(length(inew) > 0) {
- # the right number of blank rows of thermo$obigt
- newrows <- thermo$obigt[1:length(inew), ]
- # if we don't know something it's NA
- newrows[] <- NA
- # put in a default state
- newrows$state <- thermo$opt$state
- # the formula defaults to the name
- newrows$formula <- args$name[inew]
- # fill in the columns
- newrows[, icol] <- args[inew, ]
- # now check the formulas
- e <- tryCatch(makeup(newrows$formula), error=function(e) e)
- if(inherits(e, "error")) {
- warning("please supply a valid chemical formula as the species name or in the 'formula' argument")
- # transmit the error from makeup
- stop(e)
- }
- # assign to thermo$obigt
- thermo$obigt <- rbind(thermo$obigt, newrows)
- rownames(thermo$obigt) <- NULL
- assign("thermo", thermo, "CHNOSZ")
- # update ispecies
- ntotal <- nrow(thermo$obigt)
- ispecies[inew] <- (ntotal-length(inew)+1):ntotal
- # inform user
- message(paste("mod.obigt: added ", newrows$name, "(", newrows$state, ")", sep="", collapse="\n"))
- }
- if(length(iold) > 0) {
- # loop over species
- for(i in 1:length(iold)) {
- # the old values and the state
- oldprop <- thermo$obigt[ispecies[iold[i]], icol]
- state <- thermo$obigt$state[ispecies[iold[i]]]
- # tell user if they're the same, otherwise update the data entry
- if(isTRUE(all.equal(oldprop, args[iold[i], ], check.attributes=FALSE)))
- message("mod.obigt: no change for ", args$name[iold[i]], "(", state, ")")
- else {
- thermo$obigt[ispecies[iold[i]], icol] <- args[iold[i], ]
- assign("thermo", thermo, "CHNOSZ")
- message("mod.obigt: updated ", args$name[iold[i]], "(", state, ")")
- }
- }
- }
- return(ispecies)
-}
-
-add.obigt <- function(file, force=FALSE, E.units="cal") {
- # add/replace entries in thermo$obigt from values saved in a file
- # only replace if force==TRUE
- thermo <- get("thermo")
- to1 <- thermo$obigt
- id1 <- paste(to1$name,to1$state)
- to2 <- read.csv(file,as.is=TRUE)
- id2 <- paste(to2$name,to2$state)
- # check if the file is compatible with thermo$obigt
- tr <- try(rbind(to1,to2),silent=TRUE)
- if(identical(class(tr),'try-error')) stop(paste(file,"is not compatible with thermo$obigt data table."))
- # match the new species to existing ones
- does.exist <- id2 %in% id1
- ispecies.exist <- na.omit(match(id2, id1))
- nexist <- sum(does.exist)
- # convert from J if necessary
- if(tolower(E.units)=="j") {
- # loop over each row
- for(i in 1:nrow(to2)) {
- # GHS and EOS parameters
- icol <- (8:18)
- # if it's aqueous, also include omega
- if(to2$state[i]=="aq") icol <-(8:19)
- # don't touch volume (in column 12)
- icol <- icol[icol!=12]
- # convert to calories
- to2[i,icol] <- convert(to2[i,icol],"cal")
- }
- }
- # keep track of the species we've added
- inew <- numeric()
- if(force) {
- # replace existing entries
- if(nexist > 0) {
- to1[ispecies.exist, ] <- to2[does.exist, ]
- to2 <- to2[!does.exist, ]
- inew <- c(inew, ispecies.exist)
- }
- } else {
- # ignore any new entries that already exist
- to2 <- to2[!does.exist, ]
- nexist <- 0
- }
- # add new entries
- if(nrow(to2) > 0) {
- to1 <- rbind(to1, to2)
- inew <- c(inew, (length(id1)+1):nrow(to1))
- }
- # commit the change
- thermo$obigt <- to1
- rownames(thermo$obigt) <- 1:nrow(thermo$obigt)
- assign("thermo", thermo, "CHNOSZ")
- # message about file, if file argument is missing (default)
- if(missing(file)) {
- message("add.obigt: using default file:")
- message(file)
- }
- message("add.obigt: read ", length(does.exist), " rows; made ",
- nexist, " replacements, ", nrow(to2), " additions, units = ", E.units)
- message("add.obigt: use data(thermo) to restore default database")
- return(invisible(inew))
-}
-
thermo.refs <- function(key=NULL) {
## return references for thermodynamic data.
## 20110615 browse.refs() first version
Modified: pkg/CHNOSZ/data/refs.csv
===================================================================
--- pkg/CHNOSZ/data/refs.csv 2017-06-06 14:35:05 UTC (rev 204)
+++ pkg/CHNOSZ/data/refs.csv 2017-07-27 00:01:56 UTC (rev 205)
@@ -116,3 +116,4 @@
CHNOSZ.6,"J. M. Dick",2017,"CHNOSZ package documentation","dipeptides not included in slop files after slop98.dat",http://chnosz.net
CHNOSZ.7,"J. M. Dick",2017,"CHNOSZ package documentation","charge of NpO2(Oxal), La(Succ)+, NH4(Succ)-, and NpO2(Succ) as listed by @PSK99",http://chnosz.net
CHNOSZ.8,"J. M. Dick",2017,"CHNOSZ package documentation","Incorrect values of HKF a<sub>1</sub>--a<sub>4</sub> parameters for [-CH2NH2] were printed in Table 6 of @DLH06; corrected values are used here.",http://chnosz.net
+LCT17,"A. R. Lowe, J. S. Cox and P. R. Tremaine",2017,"J. Chem. Thermodynamics 112, 129-145","adenine Cp and V",https://doi.org/10.1016/j.jct.2017.04.005
Modified: pkg/CHNOSZ/data/thermo.R
===================================================================
--- pkg/CHNOSZ/data/thermo.R 2017-06-06 14:35:05 UTC (rev 204)
+++ pkg/CHNOSZ/data/thermo.R 2017-07-27 00:01:56 UTC (rev 205)
@@ -9,7 +9,7 @@
local({
# create obigt data frame
- sources_aq <- paste0(c("H2O", "inorganic", "organic", "biotic", "CHNOSZ"), "_aq")
+ sources_aq <- paste0(c("H2O", "inorganic", "organic", "biotic"), "_aq")
sources_cr <- paste0(c("inorganic", "organic"), "_cr")
sources_liq <- paste0(c("organic"), "_liq")
sources_gas <- paste0(c("inorganic", "organic"), "_gas")
Modified: pkg/CHNOSZ/demo/00Index
===================================================================
--- pkg/CHNOSZ/demo/00Index 2017-06-06 14:35:05 UTC (rev 204)
+++ pkg/CHNOSZ/demo/00Index 2017-07-27 00:01:56 UTC (rev 205)
@@ -18,3 +18,4 @@
bugstab formation potential of microbial proteins in colorectal cancer
Shh affinities of transcription factors relative to Sonic hedgehog
activity_ratios mineral stability plots with activity ratios on the axes
+adenine HKF parameters regressed from heat capacity and volume of aqueous adenine
Modified: pkg/CHNOSZ/demo/activity_ratios.R
===================================================================
--- pkg/CHNOSZ/demo/activity_ratios.R 2017-06-06 14:35:05 UTC (rev 204)
+++ pkg/CHNOSZ/demo/activity_ratios.R 2017-07-27 00:01:56 UTC (rev 205)
@@ -10,8 +10,8 @@
## Steinmann et al., 1994 (http://ccm.geoscienceworld.org/content/42/2/197)
## Garrels and Christ, p. 361 (http://www.worldcat.org/oclc/517586)
## https://wwwbrr.cr.usgs.gov/projects/GWC_coupled/phreeqc/html/final-75.html
-# the pseudospecies H4SiO4 is used as a basis species
-# (see the vignette "Regressing thermodynamic data")
+# make the pseudospecies H4SiO4 available for use as a basis species
+add.obigt("pseudo-H4SiO4")
basis(c("Al+3", "pseudo-H4SiO4", "K+", "H2O", "H+", "O2"))
species(c("gibbsite", "muscovite", "kaolinite", "pyrophyllite", "k-feldspar"))
a <- affinity(H4SiO4 = c(-6, -2, res), `K+` = c(-3, 6, res))
Added: pkg/CHNOSZ/demo/adenine.R
===================================================================
--- pkg/CHNOSZ/demo/adenine.R (rev 0)
+++ pkg/CHNOSZ/demo/adenine.R 2017-07-27 00:01:56 UTC (rev 205)
@@ -0,0 +1,117 @@
+# CHNOSZ/demos/adenine.R
+# plot thermodynamic data and model fits for aqueous adenine 20170725
+
+# notable functions in this demo:
+# EOSregress() - to regress HKF coefficients from Cp data
+# add.obigt() - to add thermodynamic data that are not in the default database
+
+# LCT17 = Lowe et al., 2017 (J. Chem. Thermodyn., doi:10.1016/j.jct.2017.04.005)
+# LH06 = LaRowe and Helgeson, 2006 (Geochim. Cosmochim. Acta, doi:10.1016/j.gca.2006.04.010)
+# HKF = Helgeson-Kirkham-Flowers equations (e.g. Am. J. Sci., doi:10.2475/ajs.281.10.1249)
+
+# adenine Cp0 and V0 from LCT17 Table 4
+AdH <- data.frame(
+ T = seq(288.15, 363.15, 5),
+ V = c(89.1, 89.9, 90.6, 91.3, 92, 92.7, 93.1, 93.6, 94.1, 94.9, 95.4, 95.9, 96.3, 96.9, 97.1, 97.8),
+ V_SD = c(1.1, 1.3, 1.1, 1, 1.1, 1, 0.9, 0.9, 0.8, 0.6, 0.7, 0.7, 0.7, 0.4, 1.1, 0.7),
+ Cp = c(207, 212, 216, 220, 224, 227, 230, 234, 237, 241, 243, 245, 248, 250, 252, 255),
+ Cp_SD = c(5, 7, 8, 7, 8, 7, 6, 7, 6, 5, 6, 6, 5, 4, 7, 5)
+)
+# functions to calculate V and Cp using density model (LCT17 Eq. 28)
+Vfun <- function(v1, v2, k, T) {
+ # gas constant (cm3 bar K-1 mol-1)
+ R <- 83.144598
+ # isothermal compressibility (bar-1)
+ beta <- water("beta", TK)$beta
+ v1 + v2 / (T - 228) - k * R * beta
+}
+Cpfun <- function(c1, c2, k, T) {
+ # gas constant (J K-1 mol-1)
+ R <- 8.3144598
+ # isobaric temperature derivative of expansivity (K-2)
+ daldT <- water("daldT", TK)$daldT
+ c1 + c2 / (T - 228) ^ 2 - k * R * T * daldT
+}
+# set up units (used for axis labels and HKF calculations)
+E.units("J")
+T.units("K")
+# temperature and pressure points for calculations
+TK <- seq(275, 425)
+P <- water("Psat", TK)$Psat
+# set up plots
+layout(matrix(1:3), heights=c(1, 8, 8))
+# title at top
+par(mar=c(0, 0, 0, 0), cex=1)
+plot.new()
+text(0.5, 0.5, "Heat capacity and volume of aqueous adenine\n(Lowe et al., 2017)", font=2)
+# settings for plot areas
+par(mar = c(4, 4, 0.5, 1), mgp = c(2.4, 0.5, 0))
+# location of x-axis tick marks
+xaxp <- c(275, 425, 6)
+
+### Cp plot (LCT17 Figures 4 and 12) ###
+plot(AdH$T, AdH$Cp, type = "p", xlim = range(TK), ylim = c(150, 350),
+ xlab = axis.label("T"), ylab = axis.label("Cp0"),
+ pch = 5, tcl = 0.3, xaxs = "i", yaxs = "i", las = 1, xaxp = xaxp
+)
+# error bars (arrows trick from https://stackoverflow.com/questions/13032777/scatter-plot-with-error-bars)
+arrows(AdH$T, AdH$Cp - AdH$Cp_SD, AdH$T, AdH$Cp + AdH$Cp_SD, length = 0.05, angle = 90, code = 3)
+# get LH06 predictions using HKF model
+LH06 <- subcrt("adenine", T = TK)$out$adenine
+lines(TK, LH06$Cp, lty = 3)
+# density model (parameters from LCT17 Table 11)
+lines(TK, Cpfun(160.4, -653239, -7930.3, TK), lty = 2)
+
+# regress HKF parameters
+# specify the terms in the HKF equations
+var <- c("invTTheta2", "TXBorn")
+# build a data frame with T, P, and Cp columns
+Cpdat <- data.frame(T = AdH[, "T"], P = 1, Cp = AdH[, "Cp"])
+# convert Cp data from J to cal
+Cpdat$Cp <- convert(Cpdat$Cp, "cal")
+# regress HKF parameters from Cp data
+HKFcoeffs <- EOSregress(Cpdat, var)$coefficients
+# get predictions from the fitted model
+Cpfit <- EOScalc(HKFcoeffs, TK, P)
+# plot the fitted model
+lines(TK, convert(Cpfit, "J"), lwd = 2, col = "green3")
+# format coefficients for legend; use scientific notation for c2 and omega
+coeffs <- format(signif(HKFcoeffs, 4), scientific = TRUE)
+# keep c1 out of scientific notation
+coeffs[1] <- signif(HKFcoeffs[[1]], 4)
+ipos <- which(coeffs >= 0)
+coeffs[ipos] <- paste("+", coeffs[ipos], sep = "")
+fun.lab <- as.expression(lapply(1:length(coeffs), function(x) {
+ EOSlab(names(HKFcoeffs)[x], coeffs[x])
+}))
+# add legend: regressed HKF coefficients
+legend("topleft", legend = fun.lab, pt.cex = 0.1, box.col = "green3")
+# add legend: lines
+legend("bottomright", lty = c(3, 2, 1), lwd = c(1, 1, 2), col = c("black", "black", "green3"), bty = "n",
+ legend = c("HKF model (LaRowe and Helgeson, 2006)",
+ "density model (Lowe et al., 2017)", "HKF model (fit using CHNOSZ)")
+)
+
+### V plot (LCT17 Figures 3 and 11) ###
+plot(AdH$T, AdH$V, type = "p", xlim = range(TK), ylim = c(85, 105),
+ xlab = axis.label("T"), ylab = axis.label("V0"),
+ pch = 5, tcl = 0.3, xaxs = "i", yaxs = "i", las = 1, xaxp = xaxp
+)
+axis(3, labels = FALSE, tcl = 0.3, xaxp = xaxp)
+axis(4, labels = FALSE, tcl = 0.3)
+arrows(AdH$T, AdH$V - AdH$V_SD, AdH$T, AdH$V + AdH$V_SD, length = 0.05, angle = 90, code = 3)
+# HKF model with coefficients from LH06
+lines(TK, LH06$V, lty = 3)
+# density model with coefficients from LCT17
+lines(TK, Vfun(73.9, -917.0, -7930.3, TK), lty = 2)
+# load HKF coefficients for adenine reported by LCT17
+# note that the Cp coefficients here are very close to those regressed above!
+add.obigt("adenine")
+LCT17 <- subcrt("adenine", T = TK)$out$adenine
+lines(TK, LCT17$V, lwd = 2, col = "royalblue")
+legend("bottomright", lty = c(3, 2, 1), lwd = c(1, 1, 2), col = c("black", "black", "royalblue"), bty = "n",
+ legend=c("HKF model (LaRowe and Helgeson, 2006)",
+ "density model (Lowe et al., 2017)", "HKF model (fit by Lowe et al., 2017 using CHNOSZ)")
+)
+# reset database
+data(thermo)
Modified: pkg/CHNOSZ/inst/NEWS
===================================================================
--- pkg/CHNOSZ/inst/NEWS 2017-06-06 14:35:05 UTC (rev 204)
+++ pkg/CHNOSZ/inst/NEWS 2017-07-27 00:01:56 UTC (rev 205)
@@ -1,10 +1,19 @@
-CHANGES IN CHNOSZ 1.1.0-3 (2017-06-06)
+CHANGES IN CHNOSZ 1.1.0-4 (2017-07-27)
--------------------------------------
- water.lines() now works for diagrams of Eh, pe, logfO2, logaO2,
logfH2, or logaH2 vs pH, T, or P. It is possible to have T or P on
either the x- or y-axis.
+- Data file CHNOSZ_aq.csv is no longer added to the database by default.
+ It is now intended to hold provisional updates. Updated heat capacity
+ and volume HKF parameters for adenine (Lowe et al., 2017) have been
+ added to this file.
+
+- Modify add.obigt() to search for species from CHNOSZ_aq.csv to add to
+ the current database. If the species is not found, add.obigt() treats
+ the argument as a file name, as before.
+
CHANGES IN CHNOSZ 1.1.0 (2017-05-04)
------------------------------------
Copied: pkg/CHNOSZ/inst/extdata/OBIGT/CHNOSZ_aq.csv (from rev 204, pkg/CHNOSZ/inst/extdata/OBIGT/CHNOSZ_aq.csv.xz)
===================================================================
--- pkg/CHNOSZ/inst/extdata/OBIGT/CHNOSZ_aq.csv (rev 0)
+++ pkg/CHNOSZ/inst/extdata/OBIGT/CHNOSZ_aq.csv 2017-07-27 00:01:56 UTC (rev 205)
@@ -0,0 +1,3 @@
+name,abbrv,formula,state,ref1,ref2,date,G,H,S,Cp,V,a1.a,a2.b,a3.c,a4.d,c1.e,c2.f,omega.lambda,z.T
+pseudo-H4SiO4,NA,H4SiO4,aq,CHNOSZ.4,NA,18.Feb.17,-312565,-346409,51.4246,-40.0964,52.1998,89.2031,-176.5071,-452.1431,101.36051,67.0854,-52.0776,0.1215745,0
+adenine,NA,C5H5N5,aq,LCT17,LH06a,26.Jul.17,74770,31235,53.41,51.63,90.6,23.53,0,-17.75,0,48.54,-3.318,-1.093,0
Deleted: pkg/CHNOSZ/inst/extdata/OBIGT/CHNOSZ_aq.csv.xz
===================================================================
--- pkg/CHNOSZ/inst/extdata/OBIGT/CHNOSZ_aq.csv.xz 2017-06-06 14:35:05 UTC (rev 204)
+++ pkg/CHNOSZ/inst/extdata/OBIGT/CHNOSZ_aq.csv.xz 2017-07-27 00:01:56 UTC (rev 205)
@@ -1 +0,0 @@
-ý7zXZ æÖ´F ! t/å£à þ Í] 7Iýúb‹¿ñ9²¶Ë¦¨•TB;qý"½qù ¦cL; ù"ß“<²ëbï¨gò;ø¥K䇥Šfœ‘1Æw?oz†ÆCþaì]IØŽ¨Ý5i‡Rw¼lù¼ßg ˆnÞÛ¿¡õö|É\¿ë%}Bèà>—œ)~)ÉÛMa=~Ô+ @ Ü,†è\²ëªZÜ"Ôþ”•edŠèí¦…È»C €C¡Ã4—È• ØÍ»å!ž- ¹ÉT‡C%Ò:Ô®€ ’›UPÚUÖô éÿ áÛ¹±Ägû YZ
\ No newline at end of file
Modified: pkg/CHNOSZ/man/EOSregress.Rd
===================================================================
--- pkg/CHNOSZ/man/EOSregress.Rd 2017-06-06 14:35:05 UTC (rev 204)
+++ pkg/CHNOSZ/man/EOSregress.Rd 2017-07-27 00:01:56 UTC (rev 205)
@@ -49,6 +49,7 @@
\details{
\code{EOSregress} uses a linear model (\code{\link{lm}}) to regress the experimental heat capacity or volume data in \code{exptdata}, which is a data frame with columns \samp{T} (temperature in degrees Kelvin), \samp{P} (pressure in bars), and \samp{Cp} or \samp{V} (heat capacity in cal/mol.K or volume in cm3/mol).
+The \samp{Cp} or \samp{V} data must be in the third column.
Only data below the temperature of \code{T.max} are included in the regression.
The regression formula is specified by a vector of names in \code{var}.
The names of the variables can be any combination of the following (listed in the order of search): variables listed in the following table, any available property of \code{\link{water}} (e.g. \samp{V}, \samp{alpha}, \samp{QBorn}), or the name of a function that can be found using \code{\link{get}} in the default environment.
@@ -130,6 +131,8 @@
var <- c("invTTheta2", "TXBorn")
# perform regression, with a temperature limit
EOSlm <- EOSregress(d, var, T.max=600)
+# calculate the Cp at some temperature and pressure
+EOScalc(EOSlm$coefficients, 298.15, 1)
# get the database values of c1, c2 and omega for CH4(aq)
CH4coeffs <- EOScoeffs("CH4", "Cp")
## make plots comparing the regressions
Added: pkg/CHNOSZ/man/add.obigt.Rd
===================================================================
--- pkg/CHNOSZ/man/add.obigt.Rd (rev 0)
+++ pkg/CHNOSZ/man/add.obigt.Rd 2017-07-27 00:01:56 UTC (rev 205)
@@ -0,0 +1,107 @@
+\encoding{UTF-8}
+\name{add.obigt}
+\alias{add.obigt}
+\alias{mod.obigt}
+\alias{today}
+\title{Functions to Work with the Thermodynamic Database}
+\description{
+Add or modify species in the thermodynamic database.
+}
+
+\usage{
+ add.obigt(file, force = TRUE, E.units = "cal")
+ mod.obigt(...)
+ today()
+}
+
+\arguments{
+ \item{file}{character, path to a file}
+ \item{force}{logical, force replacement of already existing species?}
+ \item{E.units}{character, units of energy, \samp{cal} or \samp{J}}
+ \item{...}{character or numeric, properties of species to modify in the thermodynamic database}
+}
+
+\details{
+
+\code{\link{add.obigt}} is available to update the thermodynamic database in use in the running session.
+If the \code{file} argument matches the name of a species in \code{extdata/OBIGT/CHNOSZ_aq.csv}, this species is added to the database (possibly replacing an existing species with the same name).
+Otherwise, \code{add.obigt} reads data from the specified \code{file} and adds it to \code{\link{thermo}$obigt}.
+
+The format of the file must be the same as the OBIGT files provided with CHNOSZ.
+Set \code{force} to FALSE to avoid replacing species that exist in the thermodynamic database (each unique combination of a name and a state in the database is considered to be a different species).
+Given the default setting of \code{E.units}, the function does not perform any unit conversions.
+If \code{E.units} is set to \samp{J}, then the thermodynamic parameters are converted from units of Joules to calories, as used in the CHNOSZ database.
+
+When adding species, there is no attempt made to keep the order of physical states in the database (aq-cr-liq-gas); the function simply adds new rows to the end of \code{thermo}$obigt.
+As a result, retrieving the properties of an added aqueous species using \code{\link{info}} requires an explicit \code{state="aq"} argument to that function if a species with the same name is in one of the cr, liq or gas states.
+
+\code{mod.obigt} changes one or more of the properties of species or adds species to the thermodynamic database.
+These changes are lost if you reload the database by calling \code{\link{data}(thermo)} or if you quit the \R session without saving it.
+The name of the species to add or change must be supplied as the first argument of \code{...} or as a named argument (named \samp{name}).
+When adding new species, a chemical formula should be included along with the values of any of the thermodynamic properties.
+The formula is taken from the \samp{formula} argument, or if that is missing, is taken to be the same as the \samp{name} of the species.
+An error results if the formula is not valid (i.e. can not be parsed by\code{\link{makeup}}).
+Additional arguments refer to the name of the property(s) to be updated and are matched to any part of compound column names in \code{\link{thermo}$obigt}, such as \samp{z} or \samp{T} in \samp{z.T}.
+Unless \samp{state} is specified as one of the properties, its value is taken from \code{thermo$opt$state}.
+When adding species, properties that are not specified become NA (except for \samp{state}).
+The values provided should be in the units specifed in the documentation for the \code{thermo} data object, including any order-of-magnitude scaling factors.
+
+\code{today} returns the current date in the format adopted for \code{thermo$obigt} (inherited from SUPCRT-format data files) e.g. \samp{13.May.12} for May 13, 2012.
+}
+
+\value{
+The values returned (\code{\link{invisible}}-y) by \code{mod.obigt} are the rownumbers of the affected species.
+}
+
+\seealso{ \code{\link{thermo}}, \code{\link{util.data}}, \code{\link{mod.buffer}} }
+
+\examples{\dontshow{data(thermo)}
+## modify an existing species (example only)
+ialanine <- mod.obigt("alanine", state="cr", G=0, H=0, S=0)
+# we have made the values of G, H, and S inconsistent
+# with the elemental composition of alanine, so the following
+# now produces a message about that
+info(ialanine)
+## add a species
+iCl2O <- mod.obigt("Cl2O", G=20970)
+info(iCl2O)
+# add a species with a name that is different from the formula
+mod.obigt("buckminsterfullerene", formula="C60", state="cr", date=today())
+# retrieve the species data (thermodynamic properties in this toy example are NA)
+info(info("C60"))
+# reset database
+data(thermo)
+
+# using add.obigt():
+# compare stepwise stability constants of cadmium chloride complexes
+# using data from Sverjensky et al., 1997 and Bazarkina et al., 2010
+Cdspecies <- c("Cd+2", "CdCl+", "CdCl2", "CdCl3-", "CdCl4-2")
+P <- c(1, seq(25, 1000, 25))
+SSH97 <- lapply(1:4, function(i) {
+ subcrt(c(Cdspecies[i], "Cl-", Cdspecies[i+1]),
+ c(-1, -1, 1), T=25, P=P)$out$logK
+})
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/chnosz -r 205
More information about the CHNOSZ-commits
mailing list