[CHNOSZ-commits] r241 - in pkg/CHNOSZ: . R demo inst man tests/testthat
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Oct 3 10:41:50 CEST 2017
Author: jedick
Date: 2017-10-03 10:41:50 +0200 (Tue, 03 Oct 2017)
New Revision: 241
Added:
pkg/CHNOSZ/demo/berman.R
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/R/berman.R
pkg/CHNOSZ/inst/NEWS
pkg/CHNOSZ/man/CHNOSZ-package.Rd
pkg/CHNOSZ/man/berman.Rd
pkg/CHNOSZ/man/subcrt.Rd
pkg/CHNOSZ/tests/testthat/test-subcrt.R
Log:
get tests and examples to work with 'polymorph' output from subcrt()
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2017-10-03 07:49:57 UTC (rev 240)
+++ pkg/CHNOSZ/DESCRIPTION 2017-10-03 08:41:50 UTC (rev 241)
@@ -1,6 +1,6 @@
Date: 2017-10-03
Package: CHNOSZ
-Version: 1.1.0-39
+Version: 1.1.0-40
Title: Thermodynamic Calculations for Geobiochemistry
Author: Jeffrey Dick
Maintainer: Jeffrey Dick <j3ffdick at gmail.com>
Modified: pkg/CHNOSZ/R/berman.R
===================================================================
--- pkg/CHNOSZ/R/berman.R 2017-10-03 07:49:57 UTC (rev 240)
+++ pkg/CHNOSZ/R/berman.R 2017-10-03 08:41:50 UTC (rev 241)
@@ -4,7 +4,7 @@
# in the system Na2O-K2O-CaO-MgO-FeO-Fe2O3-Al2O3-SiO2-TiO2-H2O-CO2.
# J. Petrol. 29, 445-522. https://doi.org/10.1093/petrology/29.2.445
-berman <- function(name, T = 298.15, P = 1, thisinfo=NULL, check.G=FALSE, calc.transition=TRUE, calc.disorder=TRUE, units="cal") {
+berman <- function(name, T = 298.15, P = 1, thisinfo=NULL, check.G=FALSE, calc.transition=TRUE, calc.disorder=FALSE, units="cal") {
# reference temperature and pressure
Pr <- 1
Tr <- 298.15
@@ -143,17 +143,14 @@
S <- S + Sds
V <- V + Vds
Cp <- Cp + Cpds
- } else {
-
- # FIXME: for now, we skip this check if disorder properties are calculated
-
- ### (for testing) use G = H - TS to check that integrals for H and S are written correctly
- Ga_fromHminusTS <- Ha - T * S
- # (fails with with berman("K-feldspar", T=convert(600, "K"), P=10000))
- if(!isTRUE(all.equal(Ga_fromHminusTS, Ga))) stop(paste0(name, ": incorrect integrals detected using DG = DH - T*S"))
-
}
+ ### (for testing) use G = H - TS to check that integrals for H and S are written correctly
+ Ga_fromHminusTS <- Ha - T * S
+ # FIXME: this check fails if disorder properties are calculated:
+ # berman("K-feldspar", T=convert(600, "K"), P=10000, calc.disorder=TRUE)
+ if(!isTRUE(all.equal(Ga_fromHminusTS, Ga))) stop(paste0(name, ": incorrect integrals detected using DG = DH - T*S"))
+
### thermodynamic and unit conventions used in SUPCRT ###
# use entropy of the elements in calculation of G --> Benson-Helgeson convention (DG = DH - T*DS)
Gf <- Ga + Tr * SPrTr_elements
Added: pkg/CHNOSZ/demo/berman.R
===================================================================
--- pkg/CHNOSZ/demo/berman.R (rev 0)
+++ pkg/CHNOSZ/demo/berman.R 2017-10-03 08:41:50 UTC (rev 241)
@@ -0,0 +1,28 @@
+# CHNOSZ/demo/berman.R 20171003
+# make some mineral activity diagrams using Berman (1988) and related data
+
+# using the Helgeson data
+# set up basis species
+basis(c("K+", "Al+3", "quartz", "H2O", "O2", "H+"))
+# use pH = 0 so that aK+ = aK+/aH+
+basis("pH", 0)
+# load the species
+species(c("K-feldspar", "muscovite", "kaolinite", "pyrophyllite", "andalusite"), "cr")
+# calculate affinities in aK+ - temperature space
+a <- affinity(`K+`=c(0, 6), T=c(25, 650), P=1000)
+# note that we go just past the quartz transition, but it has no effect on the diagram
+diagram(a)
+
+# now using the Berman data
+basis("SiO2", "cr_Berman")
+# it might be good to check that we have Berman's quartz and not coesite or some other SiO2 phase
+info(basis()$ispecies[3])
+# remove the Helgeson minerals
+species(delete=TRUE)
+# load the Berman minerals
+species(c("K-feldspar", "muscovite", "kaolinite", "pyrophyllite", "andalusite"), "cr_Berman")
+# calculate affinities in aK+ - temperature space
+a <- affinity(`K+`=c(0, 6), T=c(25, 650), P=1000)
+diagram(a, add=TRUE, names="", col="blue", lwd=2)
+
+legend("topleft", lty=c(1, 1), lwd=c(1, 2), col=c("black", "blue"), legend=c("Helgeson et al., 1978", "Berman, 1988"))
Modified: pkg/CHNOSZ/inst/NEWS
===================================================================
--- pkg/CHNOSZ/inst/NEWS 2017-10-03 07:49:57 UTC (rev 240)
+++ pkg/CHNOSZ/inst/NEWS 2017-10-03 08:41:50 UTC (rev 241)
@@ -1,4 +1,4 @@
-CHANGES IN CHNOSZ 1.1.0-38 (2017-10-03)
+CHANGES IN CHNOSZ 1.1.0-40 (2017-10-03)
---------------------------------------
MAJOR CHANGES:
@@ -90,6 +90,9 @@
- In water.* functions, rename the "diel" variable to "epsilon".
+- Components of subcrt() output indicating the stable polymorph of a
+ mineral are now named 'polymorph', not 'state'.
+
CLEANUP:
- To save space, taxid_names.csv has been trimmed to hold only those
Modified: pkg/CHNOSZ/man/CHNOSZ-package.Rd
===================================================================
--- pkg/CHNOSZ/man/CHNOSZ-package.Rd 2017-10-03 07:49:57 UTC (rev 240)
+++ pkg/CHNOSZ/man/CHNOSZ-package.Rd 2017-10-03 08:41:50 UTC (rev 241)
@@ -83,7 +83,7 @@
}
\section{Known Bugs}{
- \code{\link{subcrt}} does not correctly identify the stable polymorph of some minerals at high temperature.
+ \code{\link{subcrt}} does not correctly identify the stable polymorph of some minerals at high temperature (see skarn example).
The values generated by \code{\link{buffer}} may not be applied correctly by \code{\link{affinity}} in calculating the affinities of the formation reactions. (The values returned by \code{affinity(..., return.buffer=TRUE)} do appear to be correct in the examples).
}
Modified: pkg/CHNOSZ/man/berman.Rd
===================================================================
--- pkg/CHNOSZ/man/berman.Rd 2017-10-03 07:49:57 UTC (rev 240)
+++ pkg/CHNOSZ/man/berman.Rd 2017-10-03 08:41:50 UTC (rev 241)
@@ -8,7 +8,7 @@
\usage{
berman(name, T = 298.15, P = 1, thisinfo = NULL, check.G = FALSE,
- calc.transition = TRUE, calc.disorder = TRUE, units = "cal")
+ calc.transition = TRUE, calc.disorder = FALSE, units = "cal")
}
\arguments{
@@ -34,6 +34,8 @@
A warning is produced if the absolute value of the difference between tabulated and calculated DGfTrPr is greater than 1000 J/mol.
Providing \code{thisinfo} means that the mineral name is not searched in \code{thermo$obigt}, potentially saving some running time.
+
+FIXME: the default for \code{calc.disorder} is temporarily FALSE, as the code produces incorrect Gibbs energies.
}
\examples{
Modified: pkg/CHNOSZ/man/subcrt.Rd
===================================================================
--- pkg/CHNOSZ/man/subcrt.Rd 2017-10-03 07:49:57 UTC (rev 240)
+++ pkg/CHNOSZ/man/subcrt.Rd 2017-10-03 08:41:50 UTC (rev 241)
@@ -204,13 +204,10 @@
coeff=c(-1, -3, -4, -2, 3, 3, 2, 2, 3),
state=c("cr", "g", "aq", "aq", "cr", "cr", "cr", "aq", "liq"),
P=P, T=T, grid="P")
-# the results are not identical to SUPCRT92, at least because CHNOSZ
-# doesn't have volume changes for quartz, and has updated
-# parameters for species e.g. Cu+ from Shock et al., 1997
-# the following is to help detect unintended changes to the code
-# across revisions; the code *should* be fixed sometime so that
-# the last 1 becomes a 3
-stopifnot(all.equal(s$state$chalcopyrite[1:11],
+# The results are not identical to SUPCRT92, as CHNOSZ has updated
+# parameters for species e.g. Cu+ from Shock et al., 1997.
+# FIXME: the last 1 should be a 3
+stopifnot(all.equal(s$polymorphs$chalcopyrite[1:11],
c(1, 1, 1, 1, 1, 1, 2, 3, 3, 3, 1)))
## Standard Gibbs energy of reactions with HCN and
Modified: pkg/CHNOSZ/tests/testthat/test-subcrt.R
===================================================================
--- pkg/CHNOSZ/tests/testthat/test-subcrt.R 2017-10-03 07:49:57 UTC (rev 240)
+++ pkg/CHNOSZ/tests/testthat/test-subcrt.R 2017-10-03 08:41:50 UTC (rev 241)
@@ -18,7 +18,7 @@
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$state, c(1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3))
+ 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 minerals are consistent with literature", {
@@ -63,7 +63,7 @@
sout.C7 <- s.C7$out
expect_equal(sout.C7$G/1000, AS01.C7, tolerance=1e-4)
# we can also check that sulfur has expected phase transitions
- expect_equal(s.C7$state$sulfur, c(1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 3, 3))
+ 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)
More information about the CHNOSZ-commits
mailing list