[CHNOSZ-commits] r244 - in pkg/CHNOSZ: . R inst man tests/testthat
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sat Oct 7 05:48:12 CEST 2017
Author: jedick
Date: 2017-10-07 05:48:11 +0200 (Sat, 07 Oct 2017)
New Revision: 244
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/R/cgl.R
pkg/CHNOSZ/R/util.misc.R
pkg/CHNOSZ/inst/NEWS
pkg/CHNOSZ/man/util.expression.Rd
pkg/CHNOSZ/tests/testthat/test-berman.R
pkg/CHNOSZ/tests/testthat/test-eos.R
pkg/CHNOSZ/tests/testthat/test-subcrt.R
Log:
clean up and simplify code in cgl()
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2017-10-05 13:55:46 UTC (rev 243)
+++ pkg/CHNOSZ/DESCRIPTION 2017-10-07 03:48:11 UTC (rev 244)
@@ -1,6 +1,6 @@
-Date: 2017-10-05
+Date: 2017-10-07
Package: CHNOSZ
-Version: 1.1.0-42
+Version: 1.1.0-43
Title: Thermodynamic Calculations for Geobiochemistry
Author: Jeffrey Dick
Maintainer: Jeffrey Dick <j3ffdick at gmail.com>
Modified: pkg/CHNOSZ/R/cgl.R
===================================================================
--- pkg/CHNOSZ/R/cgl.R 2017-10-05 13:55:46 UTC (rev 243)
+++ pkg/CHNOSZ/R/cgl.R 2017-10-07 03:48:11 UTC (rev 244)
@@ -2,7 +2,6 @@
# 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
@@ -23,86 +22,61 @@
} else {
# in CHNOSZ, we have
# 1 cm^3 bar --> convert(1, "calories") == 0.02390057 cal
- # but SUPCRT92 in REAC92D.F uses
+ # but REAC92D.F in SUPCRT92 uses
cm3bar_to_cal <- 0.023901488 # cal
# start with NA values
values <- data.frame(matrix(NA, ncol = length(property), nrow=ncond))
colnames(values) <- property
- # additional calculations for quartz and coesite
- qtz <- quartz_coesite(PAR, T, P)
- isqtz <- !identical(qtz$V, 0)
- for(i in 1:length(property)) {
- PROP <- property[i]
- # a test for availability of the EoS parameters
- # here we assume that the parameters are in the same columns as in thermo$obigt
- # leave T transition (in 20th column) alone
- hasEOS <- any(!is.na(PAR[, 13:19]))
- # if at least one of the EoS parameters is available, zero out any NA's in the rest
- if(hasEOS) PAR[, 13:19][, is.na(PAR[, 13:19])] <- 0
- # equations for lambda adapted from HOK+98
- if(PROP == "Cp") {
- # use constant Cp if the EoS parameters are not available
- if(!hasEOS) p <- PAR$Cp
- else p <- PAR$a + PAR$b * T + PAR$c * T^-2 + PAR$d * T^-0.5 + PAR$e * T^2 + PAR$f * T^PAR$lambda
+ # a test for availability of heat capacity coefficients (a, b, c, d, e, f)
+ # based on the column assignments in thermo$obigt
+ if(any(!is.na(PAR[, 13:18]))) {
+ # we have at least one of the heat capacity coefficients;
+ # zero out any NA's in the rest (leave lambda and T of transition (columns 19-20) alone)
+ PAR[, 13:18][, is.na(PAR[, 13:18])] <- 0
+ # calculate the heat capacity and its integrals
+ Cp <- PAR$a + PAR$b*T + PAR$c*T^-2 + PAR$d*T^-0.5 + PAR$e*T^2 + PAR$f*T^PAR$lambda
+ intCpdT <- PAR$a*(T - Tr) + PAR$b*(T^2 - Tr^2)/2 + PAR$c*(1/T - 1/Tr)/-1 + PAR$d*(T^0.5 - Tr^0.5)/0.5 + PAR$e*(T^3-Tr^3)/3
+ intCpdlnT <- PAR$a*log(T / Tr) + PAR$b*(T - Tr) + PAR$c*(T^-2 - Tr^-2)/-2 + PAR$d*(T^-0.5 - Tr^-0.5)/-0.5 + PAR$e*(T^2 - Tr^2)/2
+ # do we also have the lambda parameter (Cp term with adjustable exponent on T)?
+ if(!is.na(PAR$lambda) & !identical(PAR$lambda, 0)) {
+ # equations for lambda adapted from Helgeson et al., 1998 (doi:10.1016/S0016-7037(97)00219-6)
+ if(PAR$lambda == -1) intCpdT <- intCpdT + PAR$f*log(T/Tr)
+ else intCpdT <- intCpdT - PAR$f*( T^(PAR$lambda + 1) - Tr^(PAR$lambda + 1) ) / (PAR$lambda + 1)
+ intCpdlnT <- intCpdlnT + PAR$f*(T^PAR$lambda - Tr^PAR$lambda) / PAR$lambda
}
- if(PROP == "V") {
- if(isqtz) p <- qtz$V
- else p <- rep(PAR$V, ncond)
+ } else {
+ # use constant heat capacity if the coefficients are not available
+ Cp <- PAR$Cp
+ CpdT <- PAR$Cp*(T - Tr)
+ intCpdlnT <- PAR$Cp*log(T / Tr)
+ }
+ # volume and its integrals
+ if(PAR$name %in% c("quartz", "coesite")) {
+ # volume calculations for quartz and coesite
+ qtz <- quartz_coesite(PAR, T, P)
+ V <- qtz$V
+ intVdP <- qtz$intVdP
+ intdVdTdP <- qtz$intdVdTdP
+ } else {
+ # for other minerals, volume is constant (Helgeson et al., 1978)
+ V <- rep(PAR$V, ncond)
+ # if the volume is NA, set its integrals to zero
+ if(is.na(PAR$V)) intVdP <- intdVdTdP <- numeric(ncond)
+ else {
+ intVdP <- PAR$V*(P - Pr) * cm3bar_to_cal
+ intdVdTdP <- 0
}
- if(PROP %in% c("E", "kT")) {
- p <- rep(NA, ncond)
- warning("cgl: E and/or kT of cr, gas and/or liq species are NA.")
- }
- if(PROP == "G") {
- # use constant Cp if the EoS parameters are not available
- if(!hasEOS) p <- PAR$Cp * (T - Tr - T * log(T/Tr)) else {
- # Gibbs energy integral: the value at Tref plus heat capacity terms
- p <- PAR$a * (T - Tr - T * log(T/Tr)) -
- PAR$b * (T - Tr)^2 / 2 - PAR$c * (1/T + T/Tr^2 - 2/Tr) / 2 -
- PAR$d * (T^0.5 - 0.5 * T * Tr^-0.5 - 0.5 * Tr^0.5) / -0.25 -
- PAR$e * (T^3 - 3 * T * Tr^2 + 2 * Tr^3) / 6
- }
- # use additional heat capacity term if it's defined
- if(!is.na(PAR$f) & !is.na(PAR$lambda)) if(PAR$f != 0) {
- if(PAR$lambda == -1) p <- p + PAR$f * (log(T/Tr) - T * (1/Tr - 1/T))
- else p <- p + PAR$f * ( T^(PAR$lambda + 1) - (PAR$lambda + 1) * T * Tr^PAR$lambda +
- PAR$lambda * Tr^(PAR$lambda + 1) ) / ( PAR$lambda * (PAR$lambda + 1) )
- }
- # 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 + PAR$V * (P - Pr) * cm3bar_to_cal
- p <- PAR$G + p
- }
- if(PROP == "H") {
- # use constant Cp if the EoS parameters are not available
- if(!hasEOS) p <- PAR$Cp * (T - Tr) else {
- p <- PAR$a * (T - Tr) + PAR$b * (T^2 - Tr^2) / 2 +
- PAR$c * (1/T - 1/Tr) / -1 + PAR$d * (T^0.5 - Tr^0.5) / 0.5 +
- PAR$e * (T^3 - Tr^3) / 3
- }
- if(!is.na(PAR$f) & !is.na(PAR$lambda)) if(PAR$f != 0) {
- if(PAR$lambda == -1) p <- p + PAR$f * log(T/Tr)
- else p <- p - PAR$f * ( T^(PAR$lambda + 1) - Tr^(PAR$lambda + 1) ) / (PAR$lambda + 1)
- }
- if(isqtz) p <- p + qtz$H
- else if(!is.na(PAR$V)) p <- p + PAR$V * (P-Pr) * cm3bar_to_cal
- p <- PAR$H + p
- }
- if(PROP=="S") {
- # use constant Cp if the EoS parameters are not available
- if(!hasEOS) p <- PAR$Cp * log(T/Tr) else {
- p <- PAR$a * log(T / Tr) + PAR$b * (T - Tr) +
- PAR$c * (T^-2 - Tr^-2) / -2 + PAR$e * (T^2 - Tr^2) / 2 +
- PAR$d * (T^-0.5 - Tr^-0.5) / -0.5
- }
- if(!is.na(PAR$f) & !is.na(PAR$lambda)) if(PAR$f != 0) {
- p <- p + PAR$f * (T^PAR$lambda - Tr^PAR$lambda) / PAR$lambda
- }
- p <- PAR$S + p + qtz$S
- }
- values[, i] <- p
}
+ # get the values of each of the requested thermodynamic properties
+ for(i in 1:length(property)) {
+ if(property[i] == "Cp") values[, i] <- Cp
+ if(property[i] == "V") values[, i] <- V
+ if(property[i] == "E") values[, i] <- rep(NA, ncond)
+ if(property[i] == "kT") values[, i] <- rep(NA, ncond)
+ if(property[i] == "G") values[, i] <- PAR$G - PAR$S*(T - Tr) + intCpdT - T*intCpdlnT + intVdP
+ if(property[i] == "H") values[, i] <- PAR$H + intCpdT + intVdP - T*intdVdTdP
+ if(property[i] == "S") values[, i] <- PAR$S + intCpdlnT - intdVdTdP
+ }
} # end calculations using parameters from thermo$obigt
out[[k]] <- values
} # end loop over species
@@ -112,7 +86,7 @@
### unexported function ###
# calculate GHS and V corrections for quartz and coesite 20170929
-# (these are the only mineral phases which SUPCRT92 applies a variable volume)
+# (these are the only mineral phases for which SUPCRT92 uses an inconstant volume)
quartz_coesite <- function(PAR, T, P) {
# the corrections are 0 for anything other than quartz and coesite
if(!PAR$name %in% c("quartz", "coesite")) return(list(G=0, H=0, S=0, V=0))
@@ -131,8 +105,8 @@
# constants from REAC92D.f
VPrTra <- 22.688 # VPrTr(a-quartz)
Vdiff <- 2.047 # VPrTr(a-quartz) - VPrTr(coesite)
- #k <- 38.5 # dPdTtr(a/b-quartz)
- k <- 38.45834 # calculated in CHNOSZ: dPdTtr(info("quartz"))
+ k <- 38.5 # dPdTtr(a/b-quartz)
+ #k <- 38.45834 # calculated in CHNOSZ: dPdTtr(info("quartz"))
# code adapted from REAC92D.f
qphase <- gsub("cr", "", PAR$state)
if(qphase == 2) {
@@ -158,6 +132,7 @@
k * (ba + aa * ca * k) * (T - Tr) * log((aa + P/k) / (aa + Pstar/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)
+ # note the minus sign on "SVterm" in order that intdVdTdP has the correct sign
+ list(intVdP=GVterm, intdVdTdP=-SVterm, V=V)
}
Modified: pkg/CHNOSZ/R/util.misc.R
===================================================================
--- pkg/CHNOSZ/R/util.misc.R 2017-10-05 13:55:46 UTC (rev 243)
+++ pkg/CHNOSZ/R/util.misc.R 2017-10-07 03:48:11 UTC (rev 244)
@@ -21,11 +21,12 @@
}
Ttr <- function(ispecies,P=1,dPdT=NULL) {
- # calculate a phase transition temperature
- # taking account of P (from dP/dT)
- T <- get("thermo")$obigt$z.T[ispecies]
+ # calculate a phase transition temperature for given P
+ TtrPr <- get("thermo")$obigt$z.T[ispecies]
+ # the constant slope, dP/dT
if(is.null(dPdT)) dPdT <- dPdTtr(ispecies)
- return(T + P / dPdT)
+ Pr <- 1
+ TtrPr + (P - Pr) / dPdT
}
GHS_Tr <- function(ispecies, Htr) {
Modified: pkg/CHNOSZ/inst/NEWS
===================================================================
--- pkg/CHNOSZ/inst/NEWS 2017-10-05 13:55:46 UTC (rev 243)
+++ pkg/CHNOSZ/inst/NEWS 2017-10-07 03:48:11 UTC (rev 244)
@@ -1,4 +1,4 @@
-CHANGES IN CHNOSZ 1.1.0-42 (2017-10-05)
+CHANGES IN CHNOSZ 1.1.0-42 (2017-10-06)
---------------------------------------
MAJOR CHANGES:
Modified: pkg/CHNOSZ/man/util.expression.Rd
===================================================================
--- pkg/CHNOSZ/man/util.expression.Rd 2017-10-05 13:55:46 UTC (rev 243)
+++ pkg/CHNOSZ/man/util.expression.Rd 2017-10-07 03:48:11 UTC (rev 244)
@@ -61,8 +61,10 @@
\tabular{ll}{
\samp{D} \tab Delta \cr
\samp{A} \tab bold A (chemical affinity) \cr
- \samp{p} \tab subscript italic P (isobaric heat capacity) \cr
- \samp{0} \tab degree sign (standard-state property) \cr }
+ \samp{p} \tab subscript italic P (for isobaric heat capacity) \cr
+ \samp{0} \tab degree sign (for a standard-state property) \cr
+ \samp{l} \tab subscript lambda \cr
+ }
A \samp{0} gets interpreted as a degree sign only if it does not immediately follow a number (so that e.g. \samp{2.303} can be included in an expression).
Modified: pkg/CHNOSZ/tests/testthat/test-berman.R
===================================================================
--- pkg/CHNOSZ/tests/testthat/test-berman.R 2017-10-05 13:55:46 UTC (rev 243)
+++ pkg/CHNOSZ/tests/testthat/test-berman.R 2017-10-07 03:48:11 UTC (rev 244)
@@ -65,40 +65,27 @@
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
T <- c(100, 100, 1000, 1000)
P <- c(5000, 50000, 5000, 50000)
- ## Start with uncomplicated minerals (no transitions)
- # Helgeson akermanite
- Ak_Hel_G <- c(-872485, -772662, -970152, -870328)
- Ak_Hel <- subcrt("akermanite", T=T, P=P)$out[[1]]
- 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_maxdiff(And_Ber_G, And_Ber$G, 10)
+ # anadalusite: an uncomplicated mineral (no transitions)
+ And_G <- c(-579368, -524987, -632421, -576834)
+ And <- subcrt("andalusite", "cr_Berman", T=T, P=P)$out[[1]]
+ expect_maxdiff(And_G, And$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_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]]
+ # quartz: a mineral with polymorphic transitions
+ aQz_G <- c(-202800, -179757, -223864, -200109)
+ aQz <- subcrt("quartz", "cr_Berman", T=T, P=P)$out[[1]]
# here, the high-P, low-T point suffers
- expect_maxdiff(aQz_Ber_G[-2], aQz_Ber$G[-2], 1.2)
- expect_maxdiff(aQz_Ber_G[2], aQz_Ber$G[2], 1250)
+ expect_maxdiff(aQz_G[-2], aQz$G[-2], 1.2)
+ expect_maxdiff(aQz_G[2], aQz$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]]
+ # K-feldspar: this one has disordering effects
+ Kfs_G <- c(-888115, -776324, -988950, -874777)
+ Kfs <- subcrt("K-feldspar", "cr_Berman", T=T, P=P)$out[[1]]
# we lose some accuracy at high T
- 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)
+ expect_maxdiff(Kfs_G[1:2], Kfs$G[1:2], 7)
+ expect_maxdiff(Kfs_G[3:4], Kfs$G[3:4], 1600)
+
+ # (see test-subcrt.R for feldspar and quartz tests with Helgeson/SUPCRT92 values)
})
Modified: pkg/CHNOSZ/tests/testthat/test-eos.R
===================================================================
--- pkg/CHNOSZ/tests/testthat/test-eos.R 2017-10-05 13:55:46 UTC (rev 243)
+++ pkg/CHNOSZ/tests/testthat/test-eos.R 2017-10-07 03:48:11 UTC (rev 244)
@@ -1,5 +1,19 @@
context("eos")
+test_that("cgl() with NA volume produces reasonable results", {
+ # 20171006: after rewriting much of the code in cgl(), melanterite and hydronium jarosite
+ # disappeared from the example diagram after Majzlan et al., 2006.
+ # Because the volumes are NA, the integral properties became NA,
+ # but it makes more sense to set them to zero.
+ ispecies <- info("melanterite")
+ expect_equal(info(ispecies)$V, NA_real_)
+ sout <- subcrt(ispecies, T=c(25, 25, 100, 100), P=c(1, 100, 1, 100))$out[[1]]
+ expect_false(any(is.na(sout$H)))
+ # for melanterite, which is listed in the database with zero heat capacity,
+ # all Cp and V integrals evaluate to zero
+ expect_length(unique(sout$H), 1)
+})
+
test_that("gfun() gives expected results", {
# calculate values of g and its derivatives
# up to 350 degrees C at Psat
Modified: pkg/CHNOSZ/tests/testthat/test-subcrt.R
===================================================================
--- pkg/CHNOSZ/tests/testthat/test-subcrt.R 2017-10-05 13:55:46 UTC (rev 243)
+++ pkg/CHNOSZ/tests/testthat/test-subcrt.R 2017-10-07 03:48:11 UTC (rev 244)
@@ -11,12 +11,6 @@
expect_equal(s$reaction$name, c("malic acid", "citric acid", "CO2", "water", "oxygen"))
})
-#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("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)
@@ -89,7 +83,33 @@
expect_equal(round(CHNOSZ$Cp, 1), SUPCRT_Cp)
})
-test_that("calculations for quartz are consistent with SUPCRT92", {
+test_that("calculations for quartz are nearly consistent with SUPCRT92", {
+ # using SUPCRT's equations, the alpha-beta transition occurs at
+ # 705 degC at 5000 bar and 1874 degC at 50000 bar,
+ # so here beta-quartz is stable only at T=1000, P=5000
+ T <- c(100, 100, 1000, 1000)
+ P <- c(5000, 50000, 5000, 50000)
+ SUPCRT_G <- c(-202778, -179719, -223906, -199129)
+ SUPCRT_H <- c(-214133, -191708, -199359, -177118)
+ SUPCRT_S <- c(12.3, 10.6, 31.8, 29.8)
+ SUPCRT_V <- c(22.5, 20.3, 23.7, 21.9)
+ SUPCRT_Cp <- c(12.3, 12.3, 16.9, 16.9)
+ CHNOSZ <- subcrt("quartz", T=T, P=P)$out[[1]]
+ # NOTE: Testing has shown that, where alpha-quartz is stable above Ttr(Pr) but below Ttr(P),
+ # SUPCRT92 computes the heat capacity and its integrals using parameters of beta-quartz.
+ # (see e.g. the equation for CprdT under the (Cpreg .EQ. 2) case in the Cptrms subroutine of SUPCRT).
+ # ... is that incorrect?
+ # CHNOSZ's behavior is different - it only uses the lower-T polymorph below Ttr(P);
+ # so we get different values from SUPCRT at T=1000, P=50000 (4th condition) where alpha-quartz is stable
+ # (above Ttr at 1 bar (575 degC), but below Ttr at 50000 bar)
+ expect_equal(round(CHNOSZ$G)[-4], SUPCRT_G[-4])
+ expect_equal(round(CHNOSZ$H)[-4], SUPCRT_H[-4])
+ expect_equal(round(CHNOSZ$S, 1)[-4], SUPCRT_S[-4])
+ expect_equal(round(CHNOSZ$Cp, 1)[-4], SUPCRT_Cp[-4])
+ expect_equal(round(CHNOSZ$V, 1), SUPCRT_V)
+})
+
+test_that("more calculations for quartz are nearly 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 = "
@@ -100,25 +120,11 @@
576 -214582 -209176 25.1 23.7
")
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)
+ expect_equal(round(CHNOSZ_1bar$G), S92_1bar$G)
+ expect_equal(round(CHNOSZ_1bar$H), S92_1bar$H)
+ expect_equal(round(CHNOSZ_1bar$S, 1), S92_1bar$S)
+ expect_equal(round(CHNOSZ_1bar$V, 1), S92_1bar$V)
- # 500 bar: SUPCRT shows phase transition at 587.811 deg C
- S92_500bar <- read.table(header = TRUE, text = "
- T G H S V
- 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
- ")
- 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
@@ -128,6 +134,7 @@
706 -215170 -204238 27.5 23.7
")
CHNOSZ_5000bar <- subcrt("quartz", T=seq(703, 706), P=5000)$out[[1]]
+ # NOTE: calculated values *below* the transition are different
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)
More information about the CHNOSZ-commits
mailing list