[CHNOSZ-commits] r257 - in pkg/CHNOSZ: . demo inst
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sat Oct 14 09:02:09 CEST 2017
Author: jedick
Date: 2017-10-14 09:02:09 +0200 (Sat, 14 Oct 2017)
New Revision: 257
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/demo/DEW.R
pkg/CHNOSZ/demo/lambda.R
pkg/CHNOSZ/inst/NEWS
Log:
demo/lambda.R: calculate V using finite difference derivative (dGlambda/dP)
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2017-10-13 07:31:38 UTC (rev 256)
+++ pkg/CHNOSZ/DESCRIPTION 2017-10-14 07:02:09 UTC (rev 257)
@@ -1,6 +1,6 @@
-Date: 2017-10-13
+Date: 2017-10-14
Package: CHNOSZ
-Version: 1.1.0-55
+Version: 1.1.0-56
Title: Thermodynamic Calculations for Geobiochemistry
Author: Jeffrey Dick
Maintainer: Jeffrey Dick <j3ffdick at gmail.com>
Modified: pkg/CHNOSZ/demo/DEW.R
===================================================================
--- pkg/CHNOSZ/demo/DEW.R 2017-10-13 07:31:38 UTC (rev 256)
+++ pkg/CHNOSZ/demo/DEW.R 2017-10-14 07:02:09 UTC (rev 257)
@@ -185,7 +185,7 @@
### additional checks
## check that we're within 0.1 of the QFM-2 values used by SSH14
-stopifnot(maxdiff(QFM_2[T %% 100 == 0], c(-17.0, -14.5, -12.5, -10.8, -9.4)) < 0.1)
+stopifnot(maxdiff((buf$O2-2)[!T%%100], c(-17.0, -14.5, -12.5, -10.8, -9.4)) < 0.1)
# Here are the logKs of aqueous species dissociation reactions at 600 degC and 50000 bar,
# values from EQ3NR output in Supporting Information of the paper (p. 103-109):
Modified: pkg/CHNOSZ/demo/lambda.R
===================================================================
--- pkg/CHNOSZ/demo/lambda.R 2017-10-13 07:31:38 UTC (rev 256)
+++ pkg/CHNOSZ/demo/lambda.R 2017-10-14 07:02:09 UTC (rev 257)
@@ -7,7 +7,8 @@
text(0.5, 0.5, "Effects of lambda transition in quartz, after Berman (1988) Figs. 1 and 2", cex=1.8)
opar <- par(mar=c(4, 4.5, 1, 0.5), cex=0.8)
-T <- convert(seq(0, 1400, 1), "K")
+TC <- 0:1200
+T <- convert(TC, "K")
labplot <- function(x) label.plot(x, xfrac=0.9, yfrac=0.1, paren=TRUE)
# this sets the units used for making the axis labels
@@ -19,25 +20,49 @@
# calculate properties at 1 kbar with and without transition
Qz_1bar <- berman("quartz", T=T, units="J")
Qz_1bar_notrans <- berman("quartz", T=T, calc.transition=FALSE, units="J")
-# Fig. 1a: volume
-plot(T, Qz_1bar$V, type="l", xlab=Tlab, ylab=Vlab, ylim=c(22.5, 24))
+
+## Fig. 1a: volume
+plot(TC, Qz_1bar$V, type="l", xlab=Tlab, ylab=Vlab, ylim=c(22.5, 24))
+# FIXME: why don't we get the curvature his plot for V shows?
+# Should it be in the v4 parameter (but it's zero)??
+## add data points digitized from Fig. 3B of Helgeson et al., 1978
+## and Fig. 1a of Berman, 1988
+skinner <- list(T=c(550, 560, 570, 580, 590), V=c(23.44, 23.48, 23.54, 23.72, 23.72))
+robie <- list(T=575, V=23.72)
+ghioroso <- list(T=c(0, 100, 200, 300, 400, 500, 575, 575, 675, 775, 875, 975),
+ V=c(22.7, 22.77, 22.85, 22.94, 23.1, 23.3, 23.5, 23.71, 23.69, 23.69, 23.71, 23.73))
+points(skinner$T, skinner$V, pch=19)
+points(robie$T, robie$V)
+points(ghioroso$T, ghioroso$V, pch=0)
+
+## calculate finite difference derivative of dGlambda/dP = Vlambda between 1 and 1.001 bar
+Glambda_1bar <- Qz_1bar_notrans$G - Qz_1bar$G
+Qz_1.001bar <- berman("quartz", T=T, P=1.001, units="J")
+Qz_1.001bar_notrans <- berman("quartz", T=T, P=1.001, calc.transition=FALSE, units="J")
+Glambda_1.001bar <- Qz_1.001bar_notrans$G - Qz_1.001bar$G
+dGlambdadP <- (Glambda_1.001bar - Glambda_1bar) / 0.001
+# we're using Joules, so multiply by ten to get cm^3
+Vlambda <- -10 * dGlambdadP
+VQz <- Qz_1bar$V + Vlambda
+# above 848K, we have beta quartz
+Qz_beta <- berman("quartz,beta", T=T, P=1, units="J")
+VQz[T >= 848.15] <- Qz_beta$V[T >= 848.14]
+lines(TC, VQz, lty=2)
legend("topleft", legend="1 bar", bty="n")
labplot("a")
-# TODO: why don't we get the curvature his plot for V shows?
-# Should it be in the v4 parameter (but it's zero)??
-# Fig. 1b: heat capacity
-plot(T, Qz_1bar$Cp, type="l", xlab=Tlab, ylab=Cplab)
-lines(T, Qz_1bar_notrans$Cp, lty=3)
+## Fig. 1b: heat capacity
+plot(TC, Qz_1bar$Cp, type="l", xlab=Tlab, ylab=Cplab)
+lines(TC, Qz_1bar_notrans$Cp, lty=3)
legend("topleft", legend="1 bar", bty="n")
labplot("b")
-# calculate properties at 10 kbar with and without transition
+## calculate properties at 10 kbar with and without transition
Qz_10bar <- berman("quartz", T=T, P=10000, units="J")
Qz_10bar_notrans <- berman("quartz", T=T, P=10000, calc.transition=FALSE, units="J")
-# Fig. 1c: heat capacity
-plot(T, Qz_10bar$Cp, type="l", xlab=Tlab, ylab=Cplab)
-lines(T, Qz_10bar_notrans$Cp, lty=3)
+## Fig. 1c: heat capacity
+plot(TC, Qz_10bar$Cp, type="l", xlab=Tlab, ylab=Cplab)
+lines(TC, Qz_10bar_notrans$Cp, lty=3)
legend("topleft", legend="10 kb", bty="n")
labplot("c")
Modified: pkg/CHNOSZ/inst/NEWS
===================================================================
--- pkg/CHNOSZ/inst/NEWS 2017-10-13 07:31:38 UTC (rev 256)
+++ pkg/CHNOSZ/inst/NEWS 2017-10-14 07:02:09 UTC (rev 257)
@@ -1,4 +1,4 @@
-CHANGES IN CHNOSZ 1.1.0-55 (2017-10-13)
+CHANGES IN CHNOSZ 1.1.0-56 (2017-10-14)
---------------------------------------
MAJOR CHANGES:
@@ -7,7 +7,7 @@
thermodynamic properties of minerals using equations of Berman, 1988.
- Calculations related to Berman's (1988) Figs. 1 and 2 for the lambda
- transition of quartz are presented in the new demo lambda.R.
+ transition of quartz are available in the new demo lambda.R.
- Add functions implementing the Deep Earth Water (DEW) model
(Sverjensky et al., 2014): water.DEW() and its supporting functions
@@ -22,11 +22,11 @@
the demo, the following *four* NEWS items:
- In equilibrate(), it is now possible to combine affinity calculations
- with changing activity of the balancing basis species (loga.balance).
+ with variable activity of the balancing basis species (loga.balance).
For example, in the last plot of the DEW demo, the affinity transect
involves simultaneously varying temperature and logfO2 (given as
arguments to affinity()) as well as total concentration of carbon
- (approximated by the loga.balance argument in equilibrate()).
+ (given by the loga.balance argument in equilibrate()).
COMPUTATIONAL IMPROVEMENTS:
@@ -35,15 +35,16 @@
"percent carbon" plots for systems where the species have different
carbon numbers.
-- nonideal() now has three methods for calculating activity coefficients:
- 'Helgeson' and 'Helgeson0', which utilize the "B-dot" equation of
- Helgeson, 1969, and 'Alberty' (the only method previously available).
- The 'Helgeson' method depends on the following *two* NEWS items (but
- the 'Helgeson0' method omits the B-dot term):
+- nonideal() now has three methods for calculating activity
+ coefficients: 'Helgeson' and 'Helgeson0', which utilize the
+ Debye-Huckel equation with parameters as described by Helgeson, 1969,
+ and 'Alberty' (the only method previously available). The 'Helgeson'
+ method depends on the following *two* NEWS items (but the 'Helgeson0'
+ method omits the B-dot term):
- All three water options (SUPCRT92, IAPWS95, DEW) can be used to
calculate 'A_DH' and 'B_DH', i.e. A and B coefficients in the extended
- Debye-Huckel equation of Helgeson, 1969.
+ Debye-Huckel equation as described by Helgeson, 1969.
- Add Bdot() to calculate the "B-dot" (or b_gamma) extended term
parameter for activity coefficients in NaCl solutions at high
@@ -84,13 +85,14 @@
- add.obigt() gets new argument 'species' for selecting species to add.
- Add DEW_aq.csv with aqueous species data from the DEW spreadsheet
- (May 2017 version). Species with data unchanged from the default
- database in CHNOSZ are not included. For the time being, only
- selected species are fully referenced.
+ (May 2017 version). Species with data that do not differ from the
+ default database in CHNOSZ are not included. Currently, detailed
+ reference keys and descriptions (in thermo$refs) are given for
+ selected species only.
- Add Berman_cr.csv with names and formulas for minerals whose
- thermodynamic properties are calculated using the Berman 1988
- formulation. The state "cr_Berman" is used to distinguish these
+ thermodynamic properties are calculated using the Berman, 1988
+ equations. The state "cr_Berman" is used to distinguish these
minerals from the Helgeson data. The parameters themselves are
stored in extdata/Berman/*.csv.
@@ -114,24 +116,25 @@
OTHER CHANGES:
- Add 'tol' argument to equil.reaction() (convergence tolerance for
- uniroot).
+ uniroot()).
- In hkf() and cgl(), combine 'ghs' and 'eos' arguments into single
argument named 'parameters'.
- Remove 'H2O.PrTr' and 'domega' arguments from hkf(). Properties of
- water at Pr and Tr and the logic for when to calculate derivatives of
- omega are now calculated within the function, thereby simplifying the
+ water at Pr and Tr and the logic for when to find derivatives of omega
+ are now calculated within the function, thereby simplifying the
function call from subcrt().
- Add 'H2O.props' argument to hkf(). This lists the properties of water
- to calculate (at P, T) so that the calculations aren't duplicated by
- subcrt().
+ to calculate (at P, T); these are now included in the results, in
+ order to avoid duplicating the calculations in subcrt().
- 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'.
+- Components of subcrt() output indicating the stable polymorph of
+ minerals with phase transitions are now named 'polymorph', not
+ 'state'.
- Add maxdiff() and expect_maxdiff() for calculating and testing the
maximum absolute pairwise difference between two objects.
More information about the CHNOSZ-commits
mailing list