[CHNOSZ-commits] r345 - in pkg/CHNOSZ: . R man tests/testthat
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Nov 7 14:14:21 CET 2018
Author: jedick
Date: 2018-11-07 14:14:21 +0100 (Wed, 07 Nov 2018)
New Revision: 345
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/R/NaCl.R
pkg/CHNOSZ/R/examples.R
pkg/CHNOSZ/R/solubility.R
pkg/CHNOSZ/R/subcrt.R
pkg/CHNOSZ/R/util.misc.R
pkg/CHNOSZ/man/NaCl.Rd
pkg/CHNOSZ/man/util.misc.Rd
pkg/CHNOSZ/tests/testthat/test-subcrt.R
Log:
subcrt(): fix some phase transition bugs
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2018-11-07 02:18:53 UTC (rev 344)
+++ pkg/CHNOSZ/DESCRIPTION 2018-11-07 13:14:21 UTC (rev 345)
@@ -1,6 +1,6 @@
Date: 2018-11-06
Package: CHNOSZ
-Version: 1.1.3-52
+Version: 1.1.3-53
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/NaCl.R
===================================================================
--- pkg/CHNOSZ/R/NaCl.R 2018-11-07 02:18:53 UTC (rev 344)
+++ pkg/CHNOSZ/R/NaCl.R 2018-11-07 13:14:21 UTC (rev 345)
@@ -6,7 +6,7 @@
# 20181105 use activity coefficient of Na+
# 20181106 use activity coefficient of NaCl
-NaCl <- function(T=seq(100, 500, 100), P=1000, m_tot=2) {
+NaCl <- function(T=seq(100, 500, 100), P=1000, m_tot=2, ...) {
# define a function for the reaction quotient
logQ <- function(m_Cl, gam_NaCl, gam_Na, gam_Cl) {
# starting with Q = a_NaCl / (a_Na+ * a_Cl-),
@@ -17,7 +17,7 @@
# define a function for affinity = log(K / Q)
A <- function(m_Cl, gam_NaCl, gam_Na, gam_Cl, logK) logK - logQ(m_Cl, gam_NaCl, gam_Na, gam_Cl)
# calculate equilibrium constant at all temperatures (standard conditions: IS = 0)
- logK <- subcrt(c("Na+", "Cl-", "NaCl"), c(-1, -1, 1), T = T, P = P)$out$logK
+ logK <- subcrt(c("Na+", "Cl-", "NaCl"), c(-1, -1, 1), T = T, P = P, ...)$out$logK
# calculate Debye-Huckel parameters at all temperatures
wout <- water(c("A_DH", "B_DH"), T = convert(T, "K"), P = P)
# initialize output variables
Modified: pkg/CHNOSZ/R/examples.R
===================================================================
--- pkg/CHNOSZ/R/examples.R 2018-11-07 02:18:53 UTC (rev 344)
+++ pkg/CHNOSZ/R/examples.R 2018-11-07 13:14:21 UTC (rev 345)
@@ -8,7 +8,7 @@
topics <- c("thermo", "examples",
"util.array", "util.blast", "util.data", "util.expression",
"util.fasta", "util.formula", "util.matrix", "util.misc", "util.seq", "util.units",
- "util.water", "taxonomy", "info", "protein.info", "hkf", "water", "IAPWS95", "subcrt",
+ "util.water", "taxonomy", "info", "protein.info", "hkf", "water", "IAPWS95", "subcrt", "berman",
"makeup", "basis", "swap.basis", "species", "affinity", "solubility", "equilibrate",
"diagram", "buffer", "nonideal", "NaCl", "add.protein", "protein", "ionize.aa", "yeast.aa",
"objective", "revisit", "EOSregress", "wjd")
Modified: pkg/CHNOSZ/R/solubility.R
===================================================================
--- pkg/CHNOSZ/R/solubility.R 2018-11-07 02:18:53 UTC (rev 344)
+++ pkg/CHNOSZ/R/solubility.R 2018-11-07 13:14:21 UTC (rev 345)
@@ -14,16 +14,16 @@
## unit activities of species, so we have to take away the activites
Astar <- function(i) aout$values[[i]] + aout$species$logact[i]
loga.equil <- lapply(1:length(aout$values), Astar)
+
## for a dissociation (split) on a *per reaction* (not system) basis,
## apply the divisor here and skip the if(split){} part below
## (can be used to reproduce Fig. 4 of Manning et al., 2013)
if(is.numeric(split)) loga.equil <- lapply(loga.equil, "/", split)
- # get the balancing coefficients
+ ## to output loga.balance we need the balancing coefficients
bout <- balance(aout, balance)
n.balance <- bout$n.balance
balance <- bout$balance
-
# get logarithm of total activity of the balancing basis species
logabfun <- function(loga.equil, n.balance) {
# exponentiate, multiply by n.balance, sum, logarithm
Modified: pkg/CHNOSZ/R/subcrt.R
===================================================================
--- pkg/CHNOSZ/R/subcrt.R 2018-11-07 02:18:53 UTC (rev 344)
+++ pkg/CHNOSZ/R/subcrt.R 2018-11-07 13:14:21 UTC (rev 345)
@@ -328,20 +328,22 @@
# name and state
myname <- reaction$name[i]
mystate <- reaction$state[i]
-# # don't proceed if the state is cr_Berman
-# if(mystate=="cr_Berman") next
- # if this phase is cr2 or higher, check if we're below the transition temperature
- if(!(reaction$state[i] %in% c('liq','cr','gas'))) {
- Ttr <- Ttr(iphases[i]-1,P=P,dPdT=dPdTtr(iphases[i]-1))
- if(all(is.na(Ttr))) next
- if(any(T < Ttr)) {
- status.Ttr <- "(extrapolating G)"
- if(!exceed.Ttr) {
- # put NA into the value of G
- p.cgl[[ncgl[i]]]$G[T<Ttr] <- NA
- status.Ttr <- "(using NA for G)"
- }
- #message(paste('subcrt: some points below transition temperature for',myname, mystate, status.Ttr))
+ # if we are considering multiple phases, and if this phase is cr2 or higher, check if we're below the transition temperature
+ if(length(iphases) > length(ispecies) & i > 1) {
+ if(!(reaction$state[i] %in% c('liq','cr','gas')) & reaction$name[i-1] == reaction$name[i]) {
+ # after add.obigt("SUPCRT92"), quartz cr and cr2 are not next to each other in thermo$obigt,
+ # so use iphases[i-1] here, not iphases[i]-1 20181107
+ Ttr <- Ttr(iphases[i-1], iphases[i], P=P, dPdT = dPdTtr(iphases[i-1], iphases[i]))
+ if(all(is.na(Ttr))) next
+ if(any(T < Ttr)) {
+ status.Ttr <- "(extrapolating G)"
+ if(!exceed.Ttr) {
+ # put NA into the value of G
+ p.cgl[[ncgl[i]]]$G[T<Ttr] <- NA
+ status.Ttr <- "(using NA for G)"
+ }
+ #message(paste('subcrt: some points below transition temperature for',myname, mystate, status.Ttr))
+ }
}
}
# check if we're above the temperature limit or transition temperature
@@ -351,19 +353,21 @@
# calculate Ttr at higher P if a phase transition is present
if(i < nrow(reaction)) {
# if the next one is cr2, cr3, etc we have a transition
- if(reaction$state[i+1] %in% c("cr1", "cr2", "cr3", "cr4", "cr5", "cr6", "cr7", "cr8", "cr9"))
- Ttr <- Ttr(iphases[i],P=P,dPdT=dPdTtr(iphases[i]))
- # we don't warn here about the transition
- warn.above <- FALSE
+ if(reaction$state[i+1] %in% c("cr1", "cr2", "cr3", "cr4", "cr5", "cr6", "cr7", "cr8", "cr9") & reaction$name[i+1] == reaction$name[i]) {
+ Ttr <- Ttr(iphases[i], iphases[i+1], P = P, dPdT = dPdTtr(iphases[i], iphases[i+1]))
+ # we don't warn here about the transition
+ warn.above <- FALSE
+ }
}
if(any(is.na(Ttr))) next
- if(!all(Ttr==0) & any(T >= Ttr)) {
+ if(!all(Ttr == 0) & any(T >= Ttr)) {
status.Ttr <- "(extrapolating G)"
if(!exceed.Ttr) {
- p.cgl[[ncgl[i]]]$G[T>=Ttr] <- NA
+ p.cgl[[ncgl[i]]]$G[T >= Ttr] <- NA
status.Ttr <- "(using NA for G)"
}
- if(warn.above) message(paste('subcrt: some points above temperature limit for',myname, mystate, status.Ttr))
+ Tmax <- min(T[T >= Ttr])
+ if(warn.above) message(paste("subcrt: temperature(s) of", Tmax, "K and above exceed limit for", myname, mystate, status.Ttr))
}
}
}
@@ -430,7 +434,7 @@
}
# find the minimum-energy phase at each T-P point
phasestate <- numeric()
- out.new.entry <- outprops[[1]]
+ out.new.entry <- outprops[[arephases[1]]]
for(j in 1:nrow(G)) {
ps <- which.min(as.numeric(G[j,]))
if(length(ps)==0) {
@@ -458,7 +462,7 @@
up <- unique(phasestate)
if(length(up)>1) { word <- 'are'; p.word <- 'phases' }
else { word <- 'is'; p.word <- 'phase' }
- message(paste(p.word,c2s(unique(phasestate)),word,'stable'))
+ message(paste(p.word,paste(unique(phasestate), collapse=","),word,'stable'))
} else {
# multiple phases aren't involved ... things stay the same
out.new[[i]] <- outprops[[arephases]]
@@ -469,6 +473,7 @@
isH2O.new <- c(isH2O.new,isH2O[arephases])
}
}
+
outprops <- out.new
# remove the rows that were added to keep track of phase transitions
reaction <- reaction.new[1:length(ispecies),]
Modified: pkg/CHNOSZ/R/util.misc.R
===================================================================
--- pkg/CHNOSZ/R/util.misc.R 2018-11-07 02:18:53 UTC (rev 344)
+++ pkg/CHNOSZ/R/util.misc.R 2018-11-07 13:14:21 UTC (rev 345)
@@ -2,14 +2,14 @@
# some utility functions for the CHNOSZ package
# speciate/thermo.R 20051021 jmd
-dPdTtr <- function(ispecies) {
+dPdTtr <- function(ispecies, ispecies2 = NULL) {
# calculate dP/dT for a phase transition
# (argument is index of the lower-T phase)
thermo <- get("thermo")
- pars <- info(c(ispecies, ispecies+1), check.it=FALSE)
- # if these aren't the same mineral all we can say is zero
- # actually, should be infinity ... the volume change is zero
- if(as.character(pars$name[1]) != as.character(pars$name[2])) return(Inf)
+ if(is.null(ispecies2)) ispecies2 <- ispecies + 1
+ pars <- info(c(ispecies, ispecies2), check.it=FALSE)
+ # if these aren't the same mineral, we shouldn't be here
+ if(as.character(pars$name[1]) != as.character(pars$name[2])) stop("different names for species ", ispecies, " and ", ispecies2)
# the special handling for quartz and coesite interfere with this function,
# so we convert to uppercase names to prevent cgl() from calling quartz_coesite()
pars$name <- toupper(pars$name)
@@ -20,11 +20,11 @@
return(dP.dT)
}
-Ttr <- function(ispecies,P=1,dPdT=NULL) {
+Ttr <- function(ispecies, ispecies2 = NULL, P = 1, dPdT = NULL) {
# 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)
+ if(is.null(dPdT)) dPdT <- dPdTtr(ispecies, ispecies2)
Pr <- 1
TtrPr + (P - Pr) / dPdT
}
Modified: pkg/CHNOSZ/man/NaCl.Rd
===================================================================
--- pkg/CHNOSZ/man/NaCl.Rd 2018-11-07 02:18:53 UTC (rev 344)
+++ pkg/CHNOSZ/man/NaCl.Rd 2018-11-07 13:14:21 UTC (rev 345)
@@ -7,13 +7,14 @@
}
\usage{
- NaCl(T = seq(100, 500, 100), P = 1000, m_tot = 2)
+ NaCl(T = seq(100, 500, 100), P = 1000, m_tot = 2, ...)
}
\arguments{
\item{T}{numeric, temperature in \degC}
\item{P}{numeric, pressure in bar (single value)}
\item{m_tot}{numeric, total molality of NaCl (single value)}
+ \item{...}{additional arguments for \code{\link{subcrt}}}
}
\details{
Modified: pkg/CHNOSZ/man/util.misc.Rd
===================================================================
--- pkg/CHNOSZ/man/util.misc.Rd 2018-11-07 02:18:53 UTC (rev 344)
+++ pkg/CHNOSZ/man/util.misc.Rd 2018-11-07 13:14:21 UTC (rev 345)
@@ -11,14 +11,15 @@
}
\usage{
- dPdTtr(ispecies)
- Ttr(ispecies, P = 1, dPdT = NULL)
+ dPdTtr(ispecies, ispecies2 = NULL)
+ Ttr(ispecies, ispecies2 = NULL, P = 1, dPdT = NULL)
GHS_Tr(ispecies, Htr)
unitize(logact = NULL, length = NULL, logact.tot = 0)
}
\arguments{
\item{ispecies}{numeric, species index of a mineral phase}
+ \item{ispecies2}{numeric, species index of next mineral phase (the default is ispecies + 1)}
\item{P}{numeric, pressure (bar)}
\item{dPdT}{numeric, values of (\eqn{dP/dT}{dP/dT}) of phase transitions (\code{Ttr})}
\item{Htr}{numeric, enthalpy(ies) of transition (cal/mol)}
@@ -28,7 +29,7 @@
}
\details{
-\code{dPdTtr} returns values of \eqn{(dP/dT)_{Ttr}}{(dP/dT)Ttr}, where \eqn{Ttr}{Ttr} represents the transition temperature, of the phase transition at the high-\eqn{T}{T} stability limit of the \code{ispecies} in \code{thermo$obigt} (no checking is done to verify that the species represents in fact one phase of a mineral with phase transitions).
+\code{dPdTtr} returns values of \eqn{(dP/dT)_{Ttr}}{(dP/dT)Ttr}, where \eqn{Ttr}{Ttr} represents the transition temperature, of the phase transition at the high-\eqn{T}{T} stability limit of the \code{ispecies} in \code{thermo$obigt} (other than checking that the names match, the function does not check that the species in fact represent different phases of the same mineral).
\code{dPdTtr} takes account of the Clapeyron equation, \eqn{(dP/dT)_{Ttr}}{(dP/dT)Ttr}=\eqn{{\Delta}S/{\Delta}V}{DS/DV}, where \eqn{{\Delta}S}{DS} and \eqn{{\Delta}V}{DV} represent the changes in entropy and volume of phase transition, and are calculated using \code{subcrt} at Ttr from the standard molal entropies and volumes of the two phases involved.
Using values of \code{dPdT} calculated using \code{dPdTtr} or supplied in the arguments, \code{Ttr} returns as a function of \code{P} values of the upper transition temperature of the mineral phase represented by \code{ispecies}.
@@ -42,17 +43,21 @@
}
\examples{\dontshow{data(thermo)}
-### the first example is commented because after removing most of the
-### Helgeson et al. minerals, we don't have a suitable mineral in the
-### database to demonstrate this calculation (but ice might work)
-## properties of phase transitions
-#si <- info("enstatite")
-## (dP/dT) of transitions
-#dPdTtr(si) # first transition
-#dPdTtr(si+1) # second transition
-## temperature of transitions (Ttr) as a function of P
-#Ttr(si,P=c(1,10,100,1000))
-#Ttr(si,P=c(1,10,100,1000))
+# we need the Helgeson et al., 1978 minerals for this example
+add.obigt("SUPCRT92")
+# that replaces the existing enstatite with the first phase;
+# the other phases are appended to the end of thermo$obigt
+i1 <- info("enstatite")
+i2 <- info("enstatite", "cr2")
+i3 <- info("enstatite", "cr3")
+# (dP/dT) of transitions
+dPdTtr(i1, i2) # first transition
+dPdTtr(i2, i3) # second transition
+# temperature of transitions (Ttr) as a function of P
+Ttr(i1, i2, P=c(1,10,100,1000))
+Ttr(i2, i3, P=c(1,10,100,1000))
+# restore default database
+data(OBIGT)
# calculate the GHS at Tr for the high-temperature phases of iron
# using transition enthalpies from the SUPCRT92 database (sprons92.dat)
Modified: pkg/CHNOSZ/tests/testthat/test-subcrt.R
===================================================================
--- pkg/CHNOSZ/tests/testthat/test-subcrt.R 2018-11-07 02:18:53 UTC (rev 344)
+++ pkg/CHNOSZ/tests/testthat/test-subcrt.R 2018-11-07 13:14:21 UTC (rev 345)
@@ -62,8 +62,7 @@
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_message(subcrt(iacanthite), "subcrt: temperature\\(s\\) of 623.15 K and above exceed 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))
# the reaction coefficients in the output should be unchanged 20171214
expect_equal(subcrt(c("bunsenite", "nickel", "oxygen"), c(-1, 1, 0.5))$reaction$coeff, c(-1, 1, 0.5))
@@ -89,8 +88,6 @@
})
test_that("calculations for quartz are nearly consistent with SUPCRT92", {
- # remove existing quartz so that SUPCRT92 quartz gets added with cr and cr2 next to each other
- thermo$obigt <<- thermo$obigt[thermo$obigt$name!="quartz", ]
add.obigt("SUPCRT92")
# using SUPCRT's equations, the alpha-beta transition occurs at
# 705 degC at 5000 bar and 1874 degC at 50000 bar,
@@ -119,8 +116,6 @@
})
test_that("more calculations for quartz are nearly consistent with SUPCRT92", {
- # remove existing quartz so that SUPCRT92 quartz gets added with cr and cr2 next to each other
- thermo$obigt <<- thermo$obigt[thermo$obigt$name!="quartz", ]
add.obigt("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)
@@ -194,6 +189,20 @@
expect_equal(sum(is.na(s2$out$`Na+`$logK)), 0)
})
+test_that("combining minerals with phase transitions and aqueous species with IS > 0 does not mangle output", {
+ # s2 was giving quartz an extraneous loggam column and incorrect G and logK 20181107
+ add.obigt("SUPCRT92")
+ s1 <- subcrt(c("quartz", "K+"), T=25, IS=1)
+ s2 <- subcrt(c("K+", "quartz"), T=25, IS=1)
+ expect_true(identical(colnames(s1$out[[1]]), c("T", "P", "rho", "logK", "G", "H", "S", "V", "Cp", "polymorph")))
+ expect_true(identical(colnames(s2$out[[2]]), c("T", "P", "rho", "logK", "G", "H", "S", "V", "Cp", "polymorph")))
+ expect_true(identical(colnames(s1$out[[2]]), c("T", "P", "rho", "logK", "G", "H", "S", "V", "Cp", "loggam", "IS")))
+ expect_true(identical(colnames(s2$out[[1]]), c("T", "P", "rho", "logK", "G", "H", "S", "V", "Cp", "loggam", "IS")))
+ # another one ... pyrrhotite was getting a loggam
+ expect_true(identical(colnames(subcrt(c("iron", "Na+", "Cl-", "OH-", "pyrrhotite"), T=25, IS=1)$out$pyrrhotite),
+ c("T", "P", "rho", "logK", "G", "H", "S", "V", "Cp", "polymorph")))
+})
+
# references
# Amend, J. P. and Shock, E. L. (2001)
More information about the CHNOSZ-commits
mailing list