[CHNOSZ-commits] r329 - in pkg/CHNOSZ: . R demo inst man tests/testthat
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Sep 26 02:57:25 CEST 2018
Author: jedick
Date: 2018-09-26 02:57:25 +0200 (Wed, 26 Sep 2018)
New Revision: 329
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/R/affinity.R
pkg/CHNOSZ/R/berman.R
pkg/CHNOSZ/R/diagram.R
pkg/CHNOSZ/R/mosaic.R
pkg/CHNOSZ/R/subcrt.R
pkg/CHNOSZ/R/util.affinity.R
pkg/CHNOSZ/demo/mosaic.R
pkg/CHNOSZ/inst/NEWS
pkg/CHNOSZ/man/CHNOSZ-package.Rd
pkg/CHNOSZ/man/affinity.Rd
pkg/CHNOSZ/man/subcrt.Rd
pkg/CHNOSZ/tests/testthat/test-berman.R
pkg/CHNOSZ/tests/testthat/test-mosaic.R
pkg/CHNOSZ/tests/testthat/test-subcrt.R
Log:
subcrty(), affinity(): add 'exceed.rhomin' argument
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2018-09-23 08:17:42 UTC (rev 328)
+++ pkg/CHNOSZ/DESCRIPTION 2018-09-26 00:57:25 UTC (rev 329)
@@ -1,6 +1,6 @@
-Date: 2018-09-23
+Date: 2018-09-26
Package: CHNOSZ
-Version: 1.1.3-36
+Version: 1.1.3-37
Title: Thermodynamic Calculations and Diagrams for Geo(bio)chemistry
Authors at R: c(
person("Jeffrey", "Dick", , "j3ffdick at gmail.com", role = c("aut", "cre"),
Modified: pkg/CHNOSZ/R/affinity.R
===================================================================
--- pkg/CHNOSZ/R/affinity.R 2018-09-23 08:17:42 UTC (rev 328)
+++ pkg/CHNOSZ/R/affinity.R 2018-09-26 00:57:25 UTC (rev 329)
@@ -12,8 +12,8 @@
#source("util.data.R")
#source("species.R")
-affinity <- function(...,property=NULL,sout=NULL,exceed.Ttr=FALSE,
- return.buffer=FALSE,balance="PBB",iprotein=NULL,loga.protein=-3) {
+affinity <- function(..., property=NULL, sout=NULL, exceed.Ttr=FALSE, exceed.rhomin=FALSE,
+ return.buffer=FALSE, balance="PBB", iprotein=NULL, loga.protein=-3) {
# ...: variables over which to calculate
# property: what type of energy
# (G.basis,G.species,logact.basis,logK,logQ,A)
@@ -30,7 +30,7 @@
# the argument list
args <- energy.args(list(...))
- args <- c(args,list(sout=sout,exceed.Ttr=exceed.Ttr))
+ args <- c(args, list(sout=sout, exceed.Ttr=exceed.Ttr, exceed.rhomin=exceed.rhomin))
# the species we're given
thermo <- get("thermo")
Modified: pkg/CHNOSZ/R/berman.R
===================================================================
--- pkg/CHNOSZ/R/berman.R 2018-09-23 08:17:42 UTC (rev 328)
+++ pkg/CHNOSZ/R/berman.R 2018-09-26 00:57:25 UTC (rev 329)
@@ -110,6 +110,8 @@
Tprime <- T + Td
# with the condition that Tref < Tprime < Tlambda(1bar)
iTprime <- Tref < Tprime & Tprime < Tlambda
+ # handle NA values (arising from NA in input P values e.g. Psat above Tcritical) 20180925
+ iTprime[is.na(iTprime)] <- FALSE
Tprime <- Tprime[iTprime]
Cptr[iTprime] <- Tprime * (l1 + l2 * Tprime)^2
# we got Cp, now calculate the integrations for H and S
@@ -118,6 +120,8 @@
Ttr <- T[iTtr]
Tlambda_P <- Tlambda_P[iTtr]
Td <- Td[iTtr]
+ # handle NA values 20180925
+ Tlambda_P[is.na(Tlambda_P)] <- Inf
# the upper integration limit is Tlambda_P
Ttr[Ttr >= Tlambda_P] <- Tlambda_P[Ttr >= Tlambda_P]
# derived variables
Modified: pkg/CHNOSZ/R/diagram.R
===================================================================
--- pkg/CHNOSZ/R/diagram.R 2018-09-23 08:17:42 UTC (rev 328)
+++ pkg/CHNOSZ/R/diagram.R 2018-09-26 00:57:25 UTC (rev 329)
@@ -560,8 +560,10 @@
zs <- t(predominant[, ncol(predominant):1])
if(!is.null(fill)) fill.color(xs, ys, zs, fill, ngroups)
pn <- plot.names(zs, xs, ys, names)
- # only draw the lines if there is more than one field (avoid warnings from contour)
- if(length(unique(as.vector(zs))) > 1) {
+ # only draw the lines if there is more than one field 20180923
+ # (to avoid warnings from contour, which seem to be associated with weird
+ # font metric state and subsequent errors adding e.g. subscripted text to plot)
+ if(length(na.omit(unique(as.vector(zs)))) > 1) {
if(!is.null(dotted)) plot.line(zs, xlim, ylim, dotted, col, lwd, xrange=xrange)
else contour.lines(predominant, xlim, ylim, lty=lty, col=col, lwd=lwd)
}
Modified: pkg/CHNOSZ/R/mosaic.R
===================================================================
--- pkg/CHNOSZ/R/mosaic.R 2018-09-23 08:17:42 UTC (rev 328)
+++ pkg/CHNOSZ/R/mosaic.R 2018-09-26 00:57:25 UTC (rev 329)
@@ -60,15 +60,18 @@
A.species <- affs[[1]]
if(blend) {
# calculate affinities using relative abundances of basis species
- e <- equilibrate(A.bases)
- # what is the total activity of the basis species?
- a.tot <- Reduce("+", lapply(e$loga.equil, function(x) 10^x))
- for(j in seq_along(affs)) {
- for(i in seq_along(A.species$values)) {
- # start with zero affinity
- if(j==1) A.species$values[[i]][] <- 0
- # add affinity scaled by __relative__ abundance of this basis species
- A.species$values[[i]] <- A.species$values[[i]] + affs[[j]]$values[[i]] * 10^e$loga.equil[[j]]/a.tot
+ # this isn't needed (and doesn't work) if all the affinities are NA 20180925
+ if(any(!sapply(A.species$values, is.na))) {
+ e <- equilibrate(A.bases)
+ # what is the total activity of the basis species?
+ a.tot <- Reduce("+", lapply(e$loga.equil, function(x) 10^x))
+ for(j in seq_along(affs)) {
+ for(i in seq_along(A.species$values)) {
+ # start with zero affinity
+ if(j==1) A.species$values[[i]][] <- 0
+ # add affinity scaled by __relative__ abundance of this basis species
+ A.species$values[[i]] <- A.species$values[[i]] + affs[[j]]$values[[i]] * 10^e$loga.equil[[j]]/a.tot
+ }
}
}
} else {
Modified: pkg/CHNOSZ/R/subcrt.R
===================================================================
--- pkg/CHNOSZ/R/subcrt.R 2018-09-23 08:17:42 UTC (rev 328)
+++ pkg/CHNOSZ/R/subcrt.R 2018-09-26 00:57:25 UTC (rev 329)
@@ -12,7 +12,7 @@
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,
- logact = NULL, action.unbalanced = "warn", IS = 0) {
+ exceed.rhomin = FALSE, logact = NULL, action.unbalanced = "warn", IS = 0) {
# revise the call if the states have
# come as the second argument
@@ -287,12 +287,14 @@
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
- ilowrho <- H2O.PT$rho < 350
- ilowrho[is.na(ilowrho)] <- FALSE
- if(any(ilowrho)) {
- for(i in 1:length(p.aq)) p.aq[[i]][ilowrho, ] <- NA
- if(sum(ilowrho)==1) ptext <- "pair" else ptext <- "pairs"
- warnings <- c(warnings, paste0("below minimum density for applicability of revised HKF equations (", sum(ilowrho), " T,P ", ptext, ")"))
+ if(!exceed.rhomin) {
+ ilowrho <- H2O.PT$rho < 350
+ ilowrho[is.na(ilowrho)] <- FALSE
+ if(any(ilowrho)) {
+ for(i in 1:length(p.aq)) p.aq[[i]][ilowrho, ] <- NA
+ if(sum(ilowrho)==1) ptext <- "pair" else ptext <- "pairs"
+ warnings <- c(warnings, paste0("below minimum density for applicability of revised HKF equations (", sum(ilowrho), " T,P ", ptext, ")"))
+ }
}
# calculate activity coefficients if ionic strength is not zero
if(any(IS != 0)) {
Modified: pkg/CHNOSZ/R/util.affinity.R
===================================================================
--- pkg/CHNOSZ/R/util.affinity.R 2018-09-23 08:17:42 UTC (rev 328)
+++ pkg/CHNOSZ/R/util.affinity.R 2018-09-26 00:57:25 UTC (rev 329)
@@ -17,7 +17,7 @@
### unexported functions ###
-energy <- function(what,vars,vals,lims,T=298.15,P="Psat",IS=0,sout=NULL,exceed.Ttr=FALSE,transect=FALSE) {
+energy <- function(what,vars,vals,lims,T=298.15,P="Psat",IS=0,sout=NULL,exceed.Ttr=FALSE,exceed.rhomin=FALSE,transect=FALSE) {
# 20090329 extracted from affinity() and made to
# deal with >2 dimensions (variables)
@@ -136,7 +136,7 @@
if("T" %in% vars) T <- vals[[which(vars=="T")]]
if("P" %in% vars) P <- vals[[which(vars=="P")]]
if("IS" %in% vars) IS <- vals[[which(vars=="IS")]]
- s.args <- list(species=species,property=property,T=T,P=P,IS=IS,grid=grid,convert=FALSE,exceed.Ttr=exceed.Ttr)
+ s.args <- list(species=species,property=property,T=T,P=P,IS=IS,grid=grid,convert=FALSE,exceed.Ttr=exceed.Ttr,exceed.rhomin=exceed.rhomin)
return(do.call("subcrt",s.args))
}
}
Modified: pkg/CHNOSZ/demo/mosaic.R
===================================================================
--- pkg/CHNOSZ/demo/mosaic.R 2018-09-23 08:17:42 UTC (rev 328)
+++ pkg/CHNOSZ/demo/mosaic.R 2018-09-26 00:57:25 UTC (rev 329)
@@ -22,13 +22,13 @@
# calculate affinities using the predominant basis species
# using blend=TRUE we get curvy lines, particularly at the boundaries with siderite
# compare with the plot in Garrels and Christ, 1965
-m1 <- mosaic(bases, bases2, TRUE, pH=pH, Eh=Eh, T=T)
+m1 <- mosaic(bases, bases2, blend=TRUE, pH=pH, Eh=Eh, T=T)
# make a diagram and add water stability lines
diagram(m1$A.species, lwd=2)
water.lines(m1$A.species, col="seagreen", lwd=1.5)
# show lines for Fe(aq) = 10^-4 M
species(c("Fe+2", "Fe+3"), -4)
-m2 <- mosaic(bases, bases2, TRUE, pH=pH, Eh=Eh, T=T)
+m2 <- mosaic(bases, bases2, blend=TRUE, pH=pH, Eh=Eh, T=T)
diagram(m2$A.species, add=TRUE, names=NULL)
title(main=paste("Iron oxides, sulfides and carbonate in water, log(total S) = -6,",
"log(total C)=0, after Garrels and Christ, 1965", sep="\n"))
Modified: pkg/CHNOSZ/inst/NEWS
===================================================================
--- pkg/CHNOSZ/inst/NEWS 2018-09-23 08:17:42 UTC (rev 328)
+++ pkg/CHNOSZ/inst/NEWS 2018-09-26 00:57:25 UTC (rev 329)
@@ -1,4 +1,4 @@
-CHANGES IN CHNOSZ 1.1.3-36 (2018-09-23)
+CHANGES IN CHNOSZ 1.1.3-37 (2018-09-26)
---------------------------------------
THERMODYNAMIC DATA
@@ -62,8 +62,34 @@
plot saturation lines for minerals (where their affinity equals
zero).
+BUG FIXES
+
+- Fix a bug where subcrt()$reaction$coeffs was incorrect for reactions
+ involving minerals with phase transitions. Also ensure that the output
+ reaction stoichiometry is correct for duplicated species in reactions.
+ Thanks to Grayson Boyer for the bug report.
+
+- Fix bug in nonideal() where "Zn" in formula was identified as charge.
+ Thanks to Feng Lai for reporting the incorrect behavior caused by
+ this bug.
+
+- For species in the revised HKF model, subcrt() now sets properties to
+ NA where the density of H2O is less than 0.35 g/cm3, avoiding the
+ output of bogus values in this region. Thanks to Evgeniy Bastrakov.
+
OTHER CHANGES
+- Add 'exceed.rhomin' argument to subcrt() and affinity() to enable
+ output of properties for species in the revised HKF model below 0.35
+ g/cm3.
+
+- To provide betters diagnostics for potential web apps, warning
+ messages produced by subcrt() are now available in the output of
+ affinity(), under 'sout$warnings'.
+
+- Change internal variable names in subcrt() for better readability
+ (sinfo -> ispecies, inpho -> iphases, sinph -> phasespecies).
+
- Add dumpdata() for returning/writing all packaged thermodynamic data
(including default database and optional data files).
@@ -87,28 +113,6 @@
- Keywords in basis(): Change 'CHNOPS+' to use O2 instead of e-, and add
'CHNOPSe' and 'MgCHNOPSe' for sets of basis species that have e-.
-- Change internal variable names in subcrt() for better readability
- (sinfo -> ispecies, inpho -> iphases, sinph -> phasespecies).
-
-- To provide betters diagnostics for potential web apps, warning
- messages produced by subcrt() are now available in the output of
- affinity(), under 'sout$warnings'.
-
-BUG FIXES
-
-- Fix a bug where subcrt()$reaction$coeffs was incorrect for reactions
- involving minerals with phase transitions. Also ensure that the output
- reaction stoichiometry is correct for duplicated species in reactions.
- Thanks to Grayson Boyer for the bug report.
-
-- Fix bug in nonideal() where "Zn" in formula was identified as charge.
- Thanks to Feng Lai for reporting the incorrect behavior caused by
- this bug.
-
-- For species in the revised HKF model, subcrt() now sets properties to
- NA where the density of H2O is less than 0.35 g/cm3, avoiding the
- output of bogus values in this region. Thanks to Evgeniy Bastrakov.
-
CHANGES IN CHNOSZ 1.1.3 (2017-11-13)
------------------------------------
Modified: pkg/CHNOSZ/man/CHNOSZ-package.Rd
===================================================================
--- pkg/CHNOSZ/man/CHNOSZ-package.Rd 2018-09-23 08:17:42 UTC (rev 328)
+++ pkg/CHNOSZ/man/CHNOSZ-package.Rd 2018-09-26 00:57:25 UTC (rev 329)
@@ -2,7 +2,7 @@
\name{CHNOSZ-package}
\alias{CHNOSZ-package}
\docType{package}
-\title{Thermodynamic Calculations for Geobiochemistry}
+\title{Thermodynamic Calculations and Diagrams for Geo(bio)chemistry}
\description{
CHNOSZ is a package for thermodynamic calculations, primarily with applications in geochemistry and compositional biology.
It can be used to calculate the standard molal thermodynamic properties and chemical affinities of reactions relevant to geobiochemical processes, and to visualize the equilibrium activities of species on chemical speciation and predominance diagrams.
Modified: pkg/CHNOSZ/man/affinity.Rd
===================================================================
--- pkg/CHNOSZ/man/affinity.Rd 2018-09-23 08:17:42 UTC (rev 328)
+++ pkg/CHNOSZ/man/affinity.Rd 2018-09-26 00:57:25 UTC (rev 329)
@@ -7,7 +7,7 @@
}
\usage{
- affinity(..., property=NULL, sout=NULL, exceed.Ttr=FALSE,
+ affinity(..., property=NULL, sout=NULL, exceed.Ttr=FALSE, exceed.rhomin = FALSE,
return.buffer=FALSE, balance="PBB", iprotein=NULL, loga.protein=-3)
}
@@ -16,6 +16,7 @@
\item{property}{character, the property to be calculated. Default is \samp{A}, for chemical affinity of formation reactions of species of interest}
\item{sout}{list, output from \code{\link{subcrt}}}
\item{exceed.Ttr}{logical, allow \code{\link{subcrt}} to compute properties for phases beyond their transition temperature?}
+ \item{exceed.rhomin}{logical, allow \code{\link{subcrt}} to compute properties of species in the HKF model below 0.35 g cm\S{-3}?}
\item{return.buffer}{logical. If \code{TRUE}, and a \code{\link{buffer}} has been associated with one or more basis species in the system, return the values of the activities of the basis species calculated using the buffer. Default is \code{FALSE}.}
\item{balance}{character. This argument is used to identify a conserved basis species (or \samp{PBB}) in a chemical activity buffer. Default is \samp{PBB}.}
\item{iprotein}{numeric, indices of proteins in \code{\link{thermo}$protein} for which to calculate properties}
Modified: pkg/CHNOSZ/man/subcrt.Rd
===================================================================
--- pkg/CHNOSZ/man/subcrt.Rd 2018-09-23 08:17:42 UTC (rev 328)
+++ pkg/CHNOSZ/man/subcrt.Rd 2018-09-26 00:57:25 UTC (rev 329)
@@ -10,8 +10,8 @@
subcrt(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, logact = NULL,
- action.unbalanced = "warn", IS = 0)
+ convert = TRUE, exceed.Ttr = FALSE, exceed.rhomin = FALSE,
+ logact = NULL, action.unbalanced = "warn", IS = 0)
}
\arguments{
@@ -23,6 +23,7 @@
\item{P}{numeric, pressure(s) of the calculation, or character, \samp{Psat}}
\item{grid}{character, type of \code{P}\eqn{\times}{x}\code{T} grid to produce (NULL, the default, means no gridding)}
\item{exceed.Ttr}{logical, calculate Gibbs energies of mineral phases and other species beyond their transition temperatures?}
+ \item{exceed.rhomin}{logical, return properties of species in the HKF model below 0.35 g cm\S{-3}?}
\item{logact}{numeric, logarithms of activities of species in reaction}
\item{convert}{logical, are input and output units of T and P those of the user (\code{TRUE}) (see \code{\link{T.units}}), or are they Kelvin and bar (\code{FALSE})?}
\item{action.unbalanced}{character \samp{warn} or NULL, what action to take if unbalanced reaction is provided}
Modified: pkg/CHNOSZ/tests/testthat/test-berman.R
===================================================================
--- pkg/CHNOSZ/tests/testthat/test-berman.R 2018-09-23 08:17:42 UTC (rev 328)
+++ pkg/CHNOSZ/tests/testthat/test-berman.R 2018-09-26 00:57:25 UTC (rev 329)
@@ -88,3 +88,14 @@
thermo$opt$Berman <<- NA
expect_error(berman("xxx"), "Data for xxx not available. Please add it to your_data_file.csv")
})
+
+test_that("NA values of P are handled", {
+ sresult <- subcrt("H2O", T = seq(0, 500, 100))
+ T <- sresult$out$water$T
+ P <- sresult$out$water$P
+ # this stopped with a error prior to version 1.1.3-37
+ bresult <- berman("quartz", T = convert(T, "K"), P = P)
+ expect_equal(sum(is.na(bresult$G)), 2)
+ # this also now works (producing the same NA values)
+ #subcrt("quartz", T = seq(0, 500, 100))
+})
Modified: pkg/CHNOSZ/tests/testthat/test-mosaic.R
===================================================================
--- pkg/CHNOSZ/tests/testthat/test-mosaic.R 2018-09-23 08:17:42 UTC (rev 328)
+++ pkg/CHNOSZ/tests/testthat/test-mosaic.R 2018-09-26 00:57:25 UTC (rev 329)
@@ -3,15 +3,22 @@
test_that("results are consistent with affinity()", {
basis(c("CO2", "H2O", "NH3", "O2"), c(0, 0, 0, 0))
species(c("alanine", "glycine"))
- a <- affinity()
+ a25 <- affinity()
# this is a degenerate case because we only allow NH3 to swap for NH3, and CO2 for CO2;
# however it still exercises the affinity scaling and summing code
- m1 <- mosaic("NH3", "CO2", blend=TRUE)
+ m1_25 <- mosaic("NH3", "CO2", blend=TRUE)
# this failed before we divided by loga.tot to get _relative_ abundances of basis species in mosaic.R
- expect_equal(a$values, m1$A.species$values)
+ expect_equal(a25$values, m1_25$A.species$values)
# the next call failed when which.pmax(), called by diagram(), choked on a list of length one
- m2 <- mosaic("NH3", "CO2")
- expect_equal(a$values, m2$A.species$values)
+ m2_25 <- mosaic("NH3", "CO2")
+ expect_equal(a25$values, m2_25$A.species$values)
+ # make sure the function works when all affinities are NA
+ a500 <- affinity(T=500)
+ # using blend=TRUE was failing prior to version 1.1.3-37
+ m1_500 <- mosaic("NH3", "CO2", blend=TRUE, T=500)
+ expect_equal(a500$values, m1_500$A.species$values)
+ m2_500 <- mosaic("NH3", "CO2", T=500)
+ expect_equal(a500$values, m2_500$A.species$values)
})
test_that("blend=TRUE produces reasonable values", {
Modified: pkg/CHNOSZ/tests/testthat/test-subcrt.R
===================================================================
--- pkg/CHNOSZ/tests/testthat/test-subcrt.R 2018-09-23 08:17:42 UTC (rev 328)
+++ pkg/CHNOSZ/tests/testthat/test-subcrt.R 2018-09-26 00:57:25 UTC (rev 329)
@@ -186,9 +186,12 @@
test_that("properties of HKF species below 0.35 g/cm3 are NA and give a warning", {
wtext <- "below minimum density for applicability of revised HKF equations \\(2 T,P pairs\\)"
- expect_warning(s1 <- subcrt(c("Na+", "quartz"), T=450, P=c(400, 450, 500)), wtext)
+ expect_warning(s1 <- subcrt(c("Na+", "quartz"), T=450, P=c(400, 450, 500)), wtext)
expect_equal(sum(is.na(s1$out$`Na+`$logK)), 2)
expect_equal(sum(is.na(s1$out$quartz$logK)), 0)
+ # use exceed.rhomin to go below the minimum density
+ s2 <- subcrt(c("Na+", "quartz"), T=450, P=c(400, 450, 500), exceed.rhomin=TRUE)
+ expect_equal(sum(is.na(s2$out$`Na+`$logK)), 0)
})
# references
More information about the CHNOSZ-commits
mailing list