[CHNOSZ-commits] r243 - in pkg/CHNOSZ: . R inst man tests/testthat vignettes
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Oct 5 15:55:46 CEST 2017
Author: jedick
Date: 2017-10-05 15:55:46 +0200 (Thu, 05 Oct 2017)
New Revision: 243
Added:
pkg/CHNOSZ/R/util.test.R
pkg/CHNOSZ/man/util.test.Rd
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/NAMESPACE
pkg/CHNOSZ/R/basis.R
pkg/CHNOSZ/R/cgl.R
pkg/CHNOSZ/R/info.R
pkg/CHNOSZ/R/makeup.R
pkg/CHNOSZ/inst/NEWS
pkg/CHNOSZ/tests/testthat/test-berman.R
pkg/CHNOSZ/tests/testthat/test-makeup.R
pkg/CHNOSZ/tests/testthat/test-subcrt.R
pkg/CHNOSZ/vignettes/obigt.Rmd
Log:
add expect_maxdiff() for testing maximum differences between values
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2017-10-04 15:51:38 UTC (rev 242)
+++ pkg/CHNOSZ/DESCRIPTION 2017-10-05 13:55:46 UTC (rev 243)
@@ -1,6 +1,6 @@
-Date: 2017-10-04
+Date: 2017-10-05
Package: CHNOSZ
-Version: 1.1.0-41
+Version: 1.1.0-42
Title: Thermodynamic Calculations for Geobiochemistry
Author: Jeffrey Dick
Maintainer: Jeffrey Dick <j3ffdick at gmail.com>
Modified: pkg/CHNOSZ/NAMESPACE
===================================================================
--- pkg/CHNOSZ/NAMESPACE 2017-10-04 15:51:38 UTC (rev 242)
+++ pkg/CHNOSZ/NAMESPACE 2017-10-05 13:55:46 UTC (rev 243)
@@ -58,7 +58,8 @@
"nonideal", "anim.TCA", "uniprot.aa", "run.guess",
# added 20170301 or later
"GHS_Tr", "calculateDensity", "calculateGibbsOfWater",
- "calculateEpsilon", "calculateQ", "water.DEW", "berman"
+ "calculateEpsilon", "calculateQ", "water.DEW", "berman",
+ "maxdiff", "expect_maxdiff"
)
# Load shared objects
@@ -77,5 +78,5 @@
"example", "head", "installed.packages", "read.csv", "tail",
"write.csv", "write.table")
-# Imports from colorspace
+# Imports from other packages
importFrom("colorspace", "diverge_hcl")
Modified: pkg/CHNOSZ/R/basis.R
===================================================================
--- pkg/CHNOSZ/R/basis.R 2017-10-04 15:51:38 UTC (rev 242)
+++ pkg/CHNOSZ/R/basis.R 2017-10-05 13:55:46 UTC (rev 243)
@@ -150,7 +150,7 @@
ispecies <- suppressMessages(info(as.character(thermo$buffers$species)[ibuff[k]],
as.character(thermo$buffers$state)[ibuff[k]], check.it=FALSE))
bufmakeup <- makeup(ispecies)
- inbasis <- rownames(bufmakeup) %in% colnames(basis())
+ inbasis <- names(bufmakeup) %in% colnames(basis())
if(FALSE %in% inbasis) {
stop(paste("the elements '",c2s(rownames(bufmakeup)[!inbasis]),
"' of species '",thermo$buffers$species[ibuff[k]],"' in buffer '",state[i],
Modified: pkg/CHNOSZ/R/cgl.R
===================================================================
--- pkg/CHNOSZ/R/cgl.R 2017-10-04 15:51:38 UTC (rev 242)
+++ pkg/CHNOSZ/R/cgl.R 2017-10-05 13:55:46 UTC (rev 243)
@@ -2,6 +2,7 @@
# calculate standard thermodynamic properties of non-aqueous species
# 20060729 jmd
+
cgl <- function(property = NULL, parameters = NULL, T = 298.15, P = 1) {
# calculate properties of crystalline, liquid (except H2O) and gas species
Tr <- 298.15
@@ -20,6 +21,10 @@
iprop <- match(property, colnames(properties))
values <- properties[, iprop, drop=FALSE]
} else {
+ # in CHNOSZ, we have
+ # 1 cm^3 bar --> convert(1, "calories") == 0.02390057 cal
+ # but SUPCRT92 in REAC92D.F uses
+ cm3bar_to_cal <- 0.023901488 # cal
# start with NA values
values <- data.frame(matrix(NA, ncol = length(property), nrow=ncond))
colnames(values) <- property
@@ -66,7 +71,7 @@
# entropy and volume terms
if(!is.na(PAR$S)) p <- p - PAR$S * (T - Tr)
if(isqtz) p <- p + qtz$G
- else if(!is.na(PAR$V)) p <- p + convert(PAR$V * (P - Pr), "calories")
+ else if(!is.na(PAR$V)) p <- p + PAR$V * (P - Pr) * cm3bar_to_cal
p <- PAR$G + p
}
if(PROP == "H") {
@@ -81,8 +86,7 @@
else p <- p - PAR$f * ( T^(PAR$lambda + 1) - Tr^(PAR$lambda + 1) ) / (PAR$lambda + 1)
}
if(isqtz) p <- p + qtz$H
- ## SUPCRT seems to ignore this term? ... 20070802
- #else p <- p + convert(PAR$V*(P-Pr),'calories')
+ else if(!is.na(PAR$V)) p <- p + PAR$V * (P-Pr) * cm3bar_to_cal
p <- PAR$H + p
}
if(PROP=="S") {
@@ -147,11 +151,12 @@
VPrTtb <- VPrTtb - Vdiff
V <- V - Vdiff
}
- GVterm <- convert(1, "calories") * (VPrTra * (P - Pstar) + VPrTtb * (Pstar - Pr) -
+ cm3bar_to_cal <- 0.023901488
+ GVterm <- cm3bar_to_cal * (VPrTra * (P - Pstar) + VPrTtb * (Pstar - Pr) -
0.5 * ca * (2 * Pr * (P - Pstar) - (P^2 - Pstar^2)) -
ca * k * (T - Tr) * (P - Pstar) +
k * (ba + aa * ca * k) * (T - Tr) * log((aa + P/k) / (aa + Pstar/k)))
- SVterm <- convert(1, "calories") * (-k * (ba + aa * ca * k) *
+ SVterm <- cm3bar_to_cal * (-k * (ba + aa * ca * k) *
log((aa + P/k) / (aa + Pstar/k)) + ca * k * (P - Pstar)) - Sstar
list(G=GVterm, S=SVterm, H=GVterm + T*SVterm, V=V)
}
Modified: pkg/CHNOSZ/R/info.R
===================================================================
--- pkg/CHNOSZ/R/info.R 2017-10-04 15:51:38 UTC (rev 242)
+++ pkg/CHNOSZ/R/info.R 2017-10-05 13:55:46 UTC (rev 243)
@@ -132,18 +132,22 @@
ispecies.out <- ispecies[1]
# let user know if there is more than one state for this species
if(length(ispecies) > length(ispecies.out)) {
+ mystate <- thermo$obigt$state[ispecies.out]
ispecies.other <- ispecies[!ispecies %in% ispecies.out]
otherstates <- thermo$obigt$state[ispecies.other]
transtext <- othertext <- ""
# we count, but don't show the states for phase transitions (cr2, cr3, etc)
istrans <- otherstates %in% c("cr2", "cr3", "cr4", "cr5", "cr6", "cr7", "cr8", "cr9")
- ntrans <- sum(istrans)
- if(ntrans == 1) transtext <- paste(" with", ntrans, "phase transition")
- else if(ntrans > 1) transtext <- paste(" with", ntrans, "phase transitions")
+ if(mystate=="cr") {
+ # if we are "cr" we show the number of phase transitions
+ ntrans <- sum(istrans)
+ if(ntrans == 1) transtext <- paste(" with", ntrans, "phase transition")
+ else if(ntrans > 1) transtext <- paste(" with", ntrans, "phase transitions")
+ }
otherstates <- otherstates[!istrans]
if(length(otherstates) > 0) othertext <- paste0(", also available in ", paste(otherstates, collapse=", "))
if(transtext != "" | othertext != "") {
- starttext <- paste0("info.character: found ", species, "(", thermo$obigt$state[ispecies.out], ")")
+ starttext <- paste0("info.character: found ", species, "(", mystate, ")")
message(starttext, transtext, othertext)
}
}
Modified: pkg/CHNOSZ/R/makeup.R
===================================================================
--- pkg/CHNOSZ/R/makeup.R 2017-10-04 15:51:38 UTC (rev 242)
+++ pkg/CHNOSZ/R/makeup.R 2017-10-05 13:55:46 UTC (rev 243)
@@ -68,7 +68,7 @@
mcc <- lapply(seq_along(cf), function(i) ce[[i]]*cf[i])
# unlist the subformula counts and sum them together by element
um <- unlist(mcc)
- out <- tapply(um, names(um), sum)
+ out <- unlist(tapply(um, names(um), sum, simplify=FALSE))
}
# all done with the counting, now apply the multiplier
out <- out * multiplier
@@ -123,7 +123,8 @@
count <- c(count, mycount)
# in case there are repeated elements, sum all of their counts
# (tapply hint from https://stat.ethz.ch/pipermail/r-help/2011-January/265341.html)
- out <- tapply(count, element, sum)
+ # use simplify=FALSE followed by unlist to get a vector, not array 20171005
+ out <- unlist(tapply(count, element, sum, simplify=FALSE))
# tapply returns alphabetical sorted list. keep the order appearing in the formula
out <- out[match(unique(element), names(out))]
return(out)
Added: pkg/CHNOSZ/R/util.test.R
===================================================================
--- pkg/CHNOSZ/R/util.test.R (rev 0)
+++ pkg/CHNOSZ/R/util.test.R 2017-10-05 13:55:46 UTC (rev 243)
@@ -0,0 +1,25 @@
+# CHNOSZ/util.test.R
+# functions for writing tests 20171005
+
+# the maximum (absolute) pairwise difference between x and y
+maxdiff <- function(x, y) max(abs(y - x))
+
+# modelled after the "expect_" functions in testthat ...
+# is the maximum of the pairwise differences between two objects less than some value?
+expect_maxdiff <- function(object, expected, maxdiff = 0) {
+ if(!"testthat" %in% row.names(installed.packages())) {
+ stop("please install the 'testthat' package to use this function")
+ } else {
+ # get the names for the object and expected values
+ # we use a double substitute here because
+ # as.character(substitute(x$a$b)) == c("$", "x$a", "b")
+ # as.character(substitute(substitute(x$a$b)))[2] == "x$a$b" --> better for printing
+ lab_act <- as.character(substitute(substitute(object)))[2]
+ lab_exp <- as.character(substitute(substitute(expected)))[2]
+ truemd <- maxdiff(object, expected)
+ failmsg <- sprintf("maxdiff(%s, %s) not less than %s", lab_act, lab_exp, maxdiff)
+ truemsg <- sprintf("actual value: %s", truemd)
+ testthat::expect(maxdiff(object, expected) <= maxdiff, sprintf("%s\n%s", failmsg, truemsg))
+ invisible(object)
+ }
+}
Modified: pkg/CHNOSZ/inst/NEWS
===================================================================
--- pkg/CHNOSZ/inst/NEWS 2017-10-04 15:51:38 UTC (rev 242)
+++ pkg/CHNOSZ/inst/NEWS 2017-10-05 13:55:46 UTC (rev 243)
@@ -1,4 +1,4 @@
-CHANGES IN CHNOSZ 1.1.0-41 (2017-10-04)
+CHANGES IN CHNOSZ 1.1.0-42 (2017-10-05)
---------------------------------------
MAJOR CHANGES:
@@ -98,6 +98,9 @@
- Components of subcrt() output indicating the stable polymorph of a
mineral are now named 'polymorph', not 'state'.
+- Add maxdiff() and expect_maxdiff() for calculating and testing the
+ maximum absolute pairwise difference between two objects.
+
CLEANUP:
- To save space, taxid_names.csv has been trimmed to hold only those
Added: pkg/CHNOSZ/man/util.test.Rd
===================================================================
--- pkg/CHNOSZ/man/util.test.Rd (rev 0)
+++ pkg/CHNOSZ/man/util.test.Rd 2017-10-05 13:55:46 UTC (rev 243)
@@ -0,0 +1,31 @@
+\encoding{UTF-8}
+\name{util.test}
+\alias{util.test}
+\alias{maxdiff}
+\alias{expect_maxdiff}
+\title{Functions for Writing Tests}
+\description{
+Functions modelled after the \code{expect_} functions in \CRANpkg{testthat}.
+}
+
+\usage{
+ maxdiff(x, y)
+ expect_maxdiff(object, expected, maxdiff = 0)
+}
+
+\arguments{
+ \item{x}{numeric object}
+ \item{y}{numeric object}
+ \item{object}{numeric, object to test}
+ \item{expected}{numeric, expected value}
+ \item{maxdiff}{numeric, maximum pairwise difference between object and expected value}
+}
+
+\details{
+\code{maxdiff} computes the maximum (absolute) pairwise difference between x and y, i.e. \code{max(abs(y - x))}.
+
+\code{expect_maxdiff} tests that the maximum of the pairwise differences between two objects is less than the value of the argument \code{maxdiff}.
+The function uses \code{\link[testthat]{expect}} to generate an expectation in the \CRANpkg{testthat} framework.
+}
+
+\keyword{utilities}
Modified: pkg/CHNOSZ/tests/testthat/test-berman.R
===================================================================
--- pkg/CHNOSZ/tests/testthat/test-berman.R 2017-10-04 15:51:38 UTC (rev 242)
+++ pkg/CHNOSZ/tests/testthat/test-berman.R 2017-10-05 13:55:46 UTC (rev 243)
@@ -62,8 +62,7 @@
expect_match(mineral[idiffV], "anthophyllite|antigorite|chrysotile|merwinite")
})
-
-test_that("Berman and Helgeson calculated properties are similar", {
+test_that("high-T,P calculated properties are similar to precalculated ones", {
# Reference values for G were taken from the spreadsheet Berman_Gibbs_Free_Energies.xlsx
# (http://www.dewcommunity.org/uploads/4/1/7/6/41765907/sunday_afternoon_sessions__1_.zip accessed on 2017-10-03)
# The spreadsheet also has values for some minerals using the Helgeson data
@@ -74,30 +73,32 @@
# Helgeson akermanite
Ak_Hel_G <- c(-872485, -772662, -970152, -870328)
Ak_Hel <- subcrt("akermanite", T=T, P=P)$out[[1]]
- expect_equal(Ak_Hel_G, Ak_Hel$G, tol=1e-5)
+ expect_maxdiff(Ak_Hel_G, Ak_Hel$G, 5)
# Berman andalusite
And_Ber_G <- c(-579368, -524987, -632421, -576834)
And_Ber <- subcrt("andalusite", "cr_Berman", T=T, P=P)$out[[1]]
- expect_equal(And_Ber_G, And_Ber$G, tol=1e-4)
+ expect_maxdiff(And_Ber_G, And_Ber$G, 10)
## Now a more complicated case with polymorphic transitions
# Helgeson quartz
Qz_Hel_G <- c(-202778, -179719, -223906, -199129)
Qz_Hel <- subcrt("quartz", T=T, P=P)$out[[1]]
# they're very close, but less so at the extremest condition
- expect_equal(Qz_Hel_G[-4], Qz_Hel$G[-4], tol=1e-5)
- expect_equal(Qz_Hel_G[4], Qz_Hel$G[4], tol=1e-2)
+ expect_maxdiff(Qz_Hel_G[-4], Qz_Hel$G[-4], 1)
+ expect_maxdiff(Qz_Hel_G[4], Qz_Hel$G[4], 250)
# Berman alpha-quartz
aQz_Ber_G <- c(-202800, -179757, -223864, -200109)
aQz_Ber <- subcrt("quartz", "cr_Berman", T=T, P=P)$out[[1]]
# here, the high-P, low-T point suffers
- expect_equal(aQz_Ber_G[-2], aQz_Ber$G[-2], tol=1e-5)
- expect_equal(aQz_Ber_G[2], aQz_Ber$G[2], tol=1e-2)
+ expect_maxdiff(aQz_Ber_G[-2], aQz_Ber$G[-2], 1.2)
+ expect_maxdiff(aQz_Ber_G[2], aQz_Ber$G[2], 1250)
## This one has disordering effects
+ # (see test-subcrt.R for test with Helgeson/SUPCRT92 values)
+ # Berman K-feldspar
Kfs_Ber_G <- c(-888115, -776324, -988950, -874777)
Kfs_Ber <- subcrt("K-feldspar", "cr_Berman", T=T, P=P)$out[[1]]
# we lose some accuracy at high T
- expect_equal(Kfs_Ber_G[1:2], Kfs_Ber$G[1:2], tol=1e-5)
- expect_equal(Kfs_Ber_G[3:4], Kfs_Ber$G[3:4], tol=1e-2)
+ expect_maxdiff(Kfs_Ber_G[1:2], Kfs_Ber$G[1:2], 7)
+ expect_maxdiff(Kfs_Ber_G[3:4], Kfs_Ber$G[3:4], 1600)
})
Modified: pkg/CHNOSZ/tests/testthat/test-makeup.R
===================================================================
--- pkg/CHNOSZ/tests/testthat/test-makeup.R 2017-10-04 15:51:38 UTC (rev 242)
+++ pkg/CHNOSZ/tests/testthat/test-makeup.R 2017-10-05 13:55:46 UTC (rev 243)
@@ -14,11 +14,11 @@
test_that("numeric species indices, and coefficients indicating charge can be parsed", {
# these are all equivalent formulas for the electron
expect_equal(makeup("-1"), makeup("Z0-1"))
- expect_equal(makeup("-1"), makeup("(Z-1)")[1])
+ expect_equal(makeup("-1"), makeup("(Z-1)"))
expect_equal(makeup("-1"), makeup("Z-1+0"))
# the species index of the electron in thermo$obigt
ie <- info("e-")
- expect_equal(makeup("-1"), makeup(ie)[1])
+ expect_equal(makeup("-1"), makeup(ie))
})
test_that("signed and fractional coefficients can be parsed", {
@@ -56,7 +56,7 @@
pf <- protein.formula(protein) # a matrix with elements on the columns
basis(protein) # yup, a basis set made of proteins, just for fun
bmat <- basis.elements() # a matrix with elements on the columns
- expect_equal(as.array(makeup(pf)[[1]]), makeup(as.chemical.formula(pf)[1]))
+ expect_equal(makeup(pf)[[1]], makeup(as.chemical.formula(pf)[1]))
expect_equal(makeup(pf), makeup(bmat))
})
Modified: pkg/CHNOSZ/tests/testthat/test-subcrt.R
===================================================================
--- pkg/CHNOSZ/tests/testthat/test-subcrt.R 2017-10-04 15:51:38 UTC (rev 242)
+++ pkg/CHNOSZ/tests/testthat/test-subcrt.R 2017-10-05 13:55:46 UTC (rev 243)
@@ -3,32 +3,21 @@
# delete the basis definition in case there is one
basis(delete=TRUE)
-test_that("unbalanced reactions give a warning", {
+test_that("unbalanced reactions give a warning or are balanced given sufficient basis species", {
expect_warning(subcrt(c("glucose", "ethanol"), c(-1, 3)), "reaction was unbalanced, missing H-6O3")
-})
-
-test_that("unbalanced reactions are balanced given sufficient basis species", {
basis("CHNOS")
s <- subcrt(c("malic acid", "citric acid"), c(-1, 1))
expect_equal(s$reaction$coeff, c(-1, 1, -2, -1, 1.5))
expect_equal(s$reaction$name, c("malic acid", "citric acid", "CO2", "water", "oxygen"))
})
-test_that("phase transitions of minerals give expected messages and results", {
- iacanthite <- info("acanthite", "cr2")
- #expect_message(subcrt(iacanthite), "subcrt: some points below transition temperature for acanthite cr2 \\(using NA for G\\)")
- expect_message(subcrt(iacanthite), "subcrt: some points above temperature limit for acanthite cr2 \\(using NA for G\\)")
- expect_equal(subcrt("acanthite")$out$acanthite$polymorph, c(1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3))
-})
+#test_that("heat capacity of wollastonite is consistent with the literature", {
+# # from Helgeson et al., 1978 Table 6
+# # here, set P to NA so that water() is not called to get Psat
+# expect_maxdiff(subcrt("wollastonite", T=c(200, 500, 800), P=NA)$out$wollastonite$Cp, c(25.4, 28.3, 29.9), 0.05)
+#})
-test_that("heat capacity of minerals are consistent with literature", {
- # from Helgeson et al., 1978 Table 6
- # here, set P to NA so that water() is not called to get Psat
- expect_equal(subcrt("wollastonite", T=c(200, 500, 800), P=NA)$out$wollastonite$Cp, c(25.4, 28.3, 29.9), tolerance=1e-2)
-})
-
test_that("standard Gibbs energies of reactions involving aqueous species are consistent with the literature", {
-
# from Amend and Shock, 2001 [AS01] Table 3
T <- c(2, 18, 25, 37, 45, 55, 70, 85, 100, 115, 150, 200)
# standard Gibbs energies in kJ/mol
@@ -38,40 +27,22 @@
AS01.H2O <- c(78.25, 79.34, 79.89, 80.90, 81.63, 82.59, 84.13, 85.78, 87.55, 89.42, 94.22, 102.21)
sout.H2O <- subcrt(c("H2O", "H+", "OH-"), c(-1, 1, 1), T=T)$out
# tolerances set to lowest order of magnitute to pass
- expect_equal(sout.H2O$G/1000, AS01.H2O, tolerance=1e-4)
+ expect_maxdiff(sout.H2O$G/1000, AS01.H2O, 0.01)
# AS01 Table 4.3 Reaction A1: H2(aq) + 0.5O2(aq) = H2O(l)
AS01.A1 <- c(-263.94, -263.45, -263.17, -262.62, -262.20, -261.63, -260.67, -259.60, -258.44, -257.18, -253.90, -248.44)
sout.A1 <- subcrt(c("H2", "O2", "H2O"), "aq", c(-1, -0.5, 1), T=T)$out
- expect_equal(sout.A1$G/1000, AS01.A1, tolerance=1e-4)
+ expect_maxdiff(sout.A1$G/1000, AS01.A1, 0.01)
- # AS01 Table 5.1 NO(g) -> NO(aq)
- DG0.NO.g <- c(91.39, 88.04, 86.57, 84.03, 82.33, 80.20, 76.99, 73.75, 70.50, 67.23, 59.53, 48.38)
- DG0.NO.aq <- c(104.56, 102.87, 102.06, 100.58, 99.53, 98.15, 95.98, 93.67, 91.25, 88.72, 82.46, 72.77)
- sout.NO <- subcrt("NO", c("gas", "aq"), c(-1, 1), T=T)$out
- # higher tolerance, our values for NO(aq) differ slightly from AS01
- expect_equal(sout.NO$G/1000, DG0.NO.aq - DG0.NO.g, tolerance=1e-2)
-
- # AS01 Table 5.1 Reaction B10: NH3(aq) + 1.5O2(aq) = H+ + NO2- + H2O(l)
- AS01.B10 <- c(-268.85, -268.01, -267.50, -266.46, -265.66, -264.55, -262.66, -260.54, -258.19, -255.63, -248.81, -237.06)
- sout.B10 <- subcrt(c("NH3", "O2", "H+", "NO2-", "H2O"), "aq", c(-1, -1.5, 1, 1, 1), T=T)$out
- expect_equal(sout.B10$G/1000, AS01.B10, tolerance=1e-3)
-
# AS01 Table 6.3 Reaction C7: 5S2O3-2 + H2O(l) + 4O2(aq) = 6SO4-2 + 2H+ + 4S(s)
AS01.C7 <- c(-1695.30, -1686.90, -1682.80, -1675.30, -1670.00, -1663.10, -1652.00, -1640.30, -1628.00, -1615.20, -1583.50, -1533.00)
s.C7 <- subcrt(c("S2O3-2", "H2O", "O2", "SO4-2", "H+", "S"), c("aq", "liq", "aq", "aq", "aq", "cr"), c(-5, -1, -4, 6, 2, 4), T=T)
sout.C7 <- s.C7$out
- expect_equal(sout.C7$G/1000, AS01.C7, tolerance=1e-4)
+ expect_maxdiff(sout.C7$G/1000, AS01.C7, 0.05)
# we can also check that sulfur has expected phase transitions
expect_equal(s.C7$polymorphs$sulfur, c(1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 3, 3))
- ## AS01 Table 8.3 Reaction E12: 4(2-)propanol(aq) + 3CO2(aq) + 2H2O(l) = 3CH4(aq) + 4lactic acid(aq)
- #DG0.E12 <- c(132.52, 132.26, 132.29, 132.49, 132.74, 133.15, 133.98, 135.04, 136.31, 137.79, 141.97, 149.53)
- #sout.E12 <- subcrt(c("2-propanol", "CO2", "H2O", "CH4", "lactic acid"), c(-4, -3, -2, 3, 4), T=T)$out
- ## this one is on hold until the HKF parameters of 2-propanol can be located
- #expect_equal(sout.E12$G/1000, DG0.E12, 1e-4)
-
- # return to our favourite units
+ # return to our favorite units
E.units("cal")
})
@@ -82,7 +53,7 @@
expect_equal(s.H2O$G[1], NA_real_)
# we should get something at -20 deg C
expect_equal(floor(s.H2O$G[2]), -56001)
- # for historical reasons, an input temperature of 0 was converted to 0.01
+ # following SUPCRT92, an input temperature of 0 is converted to 0.01
expect_equal(s.H2O$T[22], 0.01)
})
@@ -95,40 +66,72 @@
water(oldwat)
})
+test_that("phase transitions of minerals give expected messages and results", {
+ iacanthite <- info("acanthite", "cr2")
+ #expect_message(subcrt(iacanthite), "subcrt: some points below transition temperature for acanthite cr2 \\(using NA for G\\)")
+ expect_message(subcrt(iacanthite), "subcrt: some points above temperature limit for acanthite cr2 \\(using NA for G\\)")
+ expect_equal(subcrt("acanthite")$out$acanthite$polymorph, c(1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3))
+})
+
+test_that("calculations for K-feldspar are consistent with SUPCRT92", {
+ T <- c(100, 100, 1000, 1000)
+ P <- c(5000, 50000, 5000, 50000)
+ SUPCRT_G <- c(-886628, -769531, -988590, -871493)
+ SUPCRT_H <- c(-932344, -815247, -865868, -748771)
+ SUPCRT_S <- c(62.6, 62.6, 150.6, 150.6)
+ SUPCRT_V <- c(108.9, 108.9, 108.9, 108.9)
+ SUPCRT_Cp <- c(56.7, 56.7, 80.3, 80.3)
+ CHNOSZ <- subcrt("K-feldspar", T=T, P=P)$out[[1]]
+ expect_equal(round(CHNOSZ$G), SUPCRT_G)
+ expect_equal(round(CHNOSZ$H), SUPCRT_H)
+ expect_equal(round(CHNOSZ$S, 1), SUPCRT_S)
+ expect_equal(round(CHNOSZ$V, 1), SUPCRT_V)
+ expect_equal(round(CHNOSZ$Cp, 1), SUPCRT_Cp)
+})
+
test_that("calculations for quartz are consistent with SUPCRT92", {
# output from SUPCRT92 for reaction specified as "1 QUARTZ" run at 1 bar
# (SUPCRT shows phase transition at 574.850 deg C, and does not give Cp values around the transition)
S92_1bar <- read.table(header = TRUE, text = "
T G H S V
- 572 -214482 -209535 24.7 23.3
573 -214507 -209517 24.7 23.3
574 -214532 -209499 24.8 23.3
575 -214557 -209192 25.1 23.7
576 -214582 -209176 25.1 23.7
- 577 -214607 -209159 25.2 23.7
")
- CHNOSZ_1bar <- subcrt("quartz", T=seq(572, 577), P=1)$out[[1]]
- expect_equal(CHNOSZ_1bar$G, S92_1bar$G, tolerance = 1e-5)
- expect_equal(CHNOSZ_1bar$H, S92_1bar$H, tolerance = 1e-5)
- expect_equal(CHNOSZ_1bar$S, S92_1bar$S, tolerance = 1e-2)
- expect_equal(CHNOSZ_1bar$V, S92_1bar$V, tolerance = 1e-2)
+ CHNOSZ_1bar <- subcrt("quartz", T=seq(573, 576), P=1)$out[[1]]
+ expect_maxdiff(CHNOSZ_1bar$G, S92_1bar$G, 0.5)
+ expect_maxdiff(CHNOSZ_1bar$H, S92_1bar$H, 0.5)
+ expect_maxdiff(CHNOSZ_1bar$S, S92_1bar$S, 0.05)
+ expect_maxdiff(CHNOSZ_1bar$V, S92_1bar$V, 0.05)
- # output from SUPCRT92 for reaction specified as "1 QUARTZ" run at 500 bar
- # (SUPCRT shows phase transition at 587.811 deg C)
+ # 500 bar: SUPCRT shows phase transition at 587.811 deg C
S92_500bar <- read.table(header = TRUE, text = "
T G H S V
- 585 -214523 -209335 24.6 23.3
586 -214548 -209318 24.7 23.3
587 -214573 -209301 24.7 23.3
588 -214602 -208700 25.4 23.7
589 -214627 -208684 25.4 23.7
- 590 -214653 -208668 25.4 23.7
")
- CHNOSZ_500bar <- subcrt("quartz", T=seq(585, 590), P=500)$out[[1]]
- expect_equal(CHNOSZ_500bar$G, S92_500bar$G, tolerance = 1e-5)
- expect_equal(CHNOSZ_500bar$H, S92_500bar$H, tolerance = 1e-4)
- expect_equal(CHNOSZ_500bar$S, S92_500bar$S, tolerance = 1e-3)
- expect_equal(CHNOSZ_500bar$V, S92_500bar$V, tolerance = 1e-2)
+ CHNOSZ_500bar <- subcrt("quartz", T=seq(586, 589), P=500)$out[[1]]
+ expect_maxdiff(CHNOSZ_500bar$G, S92_500bar$G, 0.5)
+ expect_maxdiff(CHNOSZ_500bar$H, S92_500bar$H, 25)
+ expect_maxdiff(CHNOSZ_500bar$S, S92_500bar$S, 0.05)
+ expect_maxdiff(CHNOSZ_500bar$V, S92_500bar$V, 0.05)
+
+ # 5000 bar: SUPCRT shows phase transition at 704.694 deg C
+ S92_5000bar <- read.table(header = TRUE, text = "
+ T G H S V
+ 703 -215044 -204913 26.7 23.3
+ 704 -215071 -204895 26.7 23.3
+ 705 -215142 -204254 27.4 23.7
+ 706 -215170 -204238 27.5 23.7
+ ")
+ CHNOSZ_5000bar <- subcrt("quartz", T=seq(703, 706), P=5000)$out[[1]]
+ expect_maxdiff(CHNOSZ_5000bar$G, S92_5000bar$G, 20)
+ expect_maxdiff(CHNOSZ_5000bar$H, S92_5000bar$H, 300)
+ expect_maxdiff(CHNOSZ_5000bar$S, S92_5000bar$S, 0.5)
+ expect_maxdiff(CHNOSZ_5000bar$V, S92_5000bar$V, 0.05)
})
test_that("duplicated species yield correct phase transitions", {
Modified: pkg/CHNOSZ/vignettes/obigt.Rmd
===================================================================
--- pkg/CHNOSZ/vignettes/obigt.Rmd 2017-10-04 15:51:38 UTC (rev 242)
+++ pkg/CHNOSZ/vignettes/obigt.Rmd 2017-10-05 13:55:46 UTC (rev 243)
@@ -148,8 +148,8 @@
cat('Data in this this column are used in CHNOSZ only to specify the charge that is input to the "*g*-function" ([Tanger and Helgeson, 1988](https://doi.org/10.2475/ajs.288.1.19); [Shock and Helgeson, 1988](https://doi.org/10.1016/0016-7037(88)90181-0)).\n')
cat('Setting it to zero prevents activation of the *g*-function, which would result in non-zero contributions to thermodynamic properties, conflicting with the conventions mentioned above.\n')
cat('All other calculations in CHNOSZ obtain the elemental makeup, including the correct charge for the species, by parsing the chemical formulas stored in the database.<span style="color:red">^^</span>\n\n')
-cat('<span style="color:red">**</span>Likewise, [GEM-Selektor](http://gems.web.psi.ch/) defines "independent components" to be stoichiometric units usually consisting of elements and charge; the latter, [which is named Zz](http://gems.web.psi.ch/tests/TestNaCl-dep.html) and has a standard molal entropy of -65.34 J/mol/K and heat capacity of -14.418 J/mol/K (negative one-half those of gaseous hydrogen), is negated in the formula of the hypothetical "aqueous electron" ([Kulik, 2006](https://doi.org/10.1016/j.chemgeo.2005.08.014)).\n\n')
-cat('<span style="color:red">^^</span> Relatedly, charged amino acid sidechain groups have a charge that is tabulated as zero, because other values would be incompatible with group additivity of cations and anions (for which the *g*-function produces derivatives of the omega parameter (ω) in the revised HKF equations of state that are not opposite each other) to give a neutral species (for which the derivatives of ω are taken to be zero) (cf. [Dick et al., 2006](https://doi.org/10.5194/bg-3-311-2006)).\n')
+cat('<span style="color:red">**</span>Likewise, [GEM-Selektor](http://gems.web.psi.ch/) defines "independent components" to be stoichiometric units usually consisting of elements and charge; the latter, [which is named Zz](http://gems.web.psi.ch/tests/TestNaCl-dep.html) and has a standard molal entropy of -65.34 J/mol/K and heat capacity of -14.418 J/mol/K (negative one-half those of gaseous hydrogen), is negated in the formula of the fictive "aqueous electron" ([Kulik, 2006](https://doi.org/10.1016/j.chemgeo.2005.08.014)).\n\n')
+cat('<span style="color:red">^^</span> Relatedly, charged amino acid sidechain groups have a charge that is tabulated as zero, because other values would be incompatible with group additivity of cations and anions (which have derivatives of the omega parameter (ω) in the revised HKF equations of state that are not opposites of each other) to give a neutral species (for which the derivatives of ω are taken to be zero) (cf. [Dick et al., 2006](https://doi.org/10.5194/bg-3-311-2006)).\n')
```
```{r used, include=FALSE}
# initialize the list of used species
More information about the CHNOSZ-commits
mailing list