[CHNOSZ-commits] r349 - in pkg/CHNOSZ: . R demo man tests/testthat
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sun Nov 11 03:22:40 CET 2018
Author: jedick
Date: 2018-11-11 03:22:09 +0100 (Sun, 11 Nov 2018)
New Revision: 349
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/R/basis.R
pkg/CHNOSZ/R/mosaic.R
pkg/CHNOSZ/R/swap.basis.R
pkg/CHNOSZ/demo/gold.R
pkg/CHNOSZ/man/berman.Rd
pkg/CHNOSZ/tests/testthat/test-swap.basis.R
Log:
demo/gold.R: estimate molality of Cl-
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2018-11-08 13:59:57 UTC (rev 348)
+++ pkg/CHNOSZ/DESCRIPTION 2018-11-11 02:22:09 UTC (rev 349)
@@ -1,6 +1,6 @@
-Date: 2018-11-08
+Date: 2018-11-11
Package: CHNOSZ
-Version: 1.1.3-56
+Version: 1.1.3-57
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/basis.R
===================================================================
--- pkg/CHNOSZ/R/basis.R 2018-11-08 13:59:57 UTC (rev 348)
+++ pkg/CHNOSZ/R/basis.R 2018-11-11 02:22:09 UTC (rev 349)
@@ -185,7 +185,11 @@
}
}
# then modify the logact
- if(!is.null(logact)) thermo$basis$logact[ib] <- as.numeric(logact[i])
+ if(!is.null(logact)) {
+ # allow this to be non-numeric in case we're called by swap.basis() while a buffer is active 20181109
+ if(can.be.numeric(logact[i])) thermo$basis$logact[ib] <- as.numeric(logact[i])
+ else thermo$basis$logact[ib] <- logact[i]
+ }
# assign the result to the CHNOSZ environment
assign("thermo", thermo, "CHNOSZ")
}
Modified: pkg/CHNOSZ/R/mosaic.R
===================================================================
--- pkg/CHNOSZ/R/mosaic.R 2018-11-08 13:59:57 UTC (rev 348)
+++ pkg/CHNOSZ/R/mosaic.R 2018-11-11 02:22:09 UTC (rev 349)
@@ -33,6 +33,7 @@
A.bases2 <- mcall$A.bases
}
# change the basis species; restore the original at the end of the loop
+ if(can.be.numeric(logact.swap)) logact.swap <- as.numeric(logact.swap)
if(i < length(bases)) {
swap.basis(bases[i], bases[i+1])
# TODO: basis() requires the formula to identify the basis species,
Modified: pkg/CHNOSZ/R/swap.basis.R
===================================================================
--- pkg/CHNOSZ/R/swap.basis.R 2018-11-08 13:59:57 UTC (rev 348)
+++ pkg/CHNOSZ/R/swap.basis.R 2018-11-11 02:22:09 UTC (rev 349)
@@ -94,13 +94,19 @@
ispecies <- oldbasis$ispecies
ispecies[ib] <- ispecies2
newbasis <- put.basis(ispecies)
- # if put.basis didn't stop with an error, we're good to go!
- # what were the original chemical potentials of the elements?
- emu <- element.mu(oldbasis, T=T)
- # the corresponding logarithms of activities of the new basis species
- bl <- basis.logact(emu, newbasis, T=T)
+ # now deal with the activities
+ if(!all(can.be.numeric(oldbasis$logact))) {
+ # if there are any buffers, just set the old activities
+ bl <- oldbasis$logact
+ } else {
+ # no buffers, so we can recalculate activities to maintain the chemical potentials of the elements
+ # what were the original chemical potentials of the elements?
+ emu <- element.mu(oldbasis, T=T)
+ # the corresponding logarithms of activities of the new basis species
+ bl <- basis.logact(emu, newbasis, T=T)
+ }
# update the basis with these logacts
- mb <- mod.basis(ispecies, logact=bl)
+ mb <- mod.basis(ispecies, state = newbasis$state, logact = bl)
# delete, then restore species if they were defined
species(delete=TRUE)
if(!is.null(ts)) {
Modified: pkg/CHNOSZ/demo/gold.R
===================================================================
--- pkg/CHNOSZ/demo/gold.R 2018-11-08 13:59:57 UTC (rev 348)
+++ pkg/CHNOSZ/demo/gold.R 2018-11-11 02:22:09 UTC (rev 349)
@@ -1,5 +1,6 @@
# CHNOSZ/demo/gold.R: Au solubility calculations
# 20181101 jmd first version
+# 20181109 add calculation of K+ molality
## additions to OBIGT:
# Au(HS) from Akinfiev and Zotov, 2010
@@ -19,12 +20,13 @@
## modifications to OBIGT:
# AuCl2- from Akinfiev and Zotov, 2001 (reported in AZ10)
# (http://pleiades.online/cgi-perl/search.pl/?type=abstract&name=geochem&number=10&year=1&page=990)
-mod.obigt("AuCl2-", formula = "AuCl2-", state = "aq", ref1 = "AZ01", date = today(),
+mod.obigt("AuCl2-", ref1 = "AZ01", date = today(),
G = -36795, H = -46664, S = 47.16, Cp = -26.4, V = 68.6,
a1 = 11.4774, a2 = 20.2425, a3 = -2.2063, a4 = -3.6158,
c1 = 27.0677, c2 = -22.240, omega = 0.8623, z = -1)
# Au(HS)2- from Pokrovski et al., 2014
-mod.obigt("Au(HS)2-", G = 3487, H = 4703, S = 77.46, Cp = 3.3, V = 75.1,
+mod.obigt("Au(HS)2-", ref1 = "PAB+14", date = today(),
+ G = 3487, H = 4703, S = 77.46, Cp = 3.3, V = 75.1,
a1 = 12.3373, a2 = 22.3421, a3 = 3.0317, a4 = -3.7026,
c1 = -53.6010, c2 = 31.4030, omega = 0.7673, z = -1)
@@ -123,6 +125,18 @@
title(main=("After Stef\u00e1nsson and Seward, 2004, Fig. 12b"), font.main = 1, cex.main = 1.1)
}
+# estimate the Cl- molality and ionic strength for a hypothetical
+# NaCl solution with total chloride equal to specified NaCl + KCl solution,
+# then estimate the molality of K+ in that solution 20181109
+chloride <- function(T, P, m_NaCl, m_KCl) {
+ NaCl <- NaCl(T = T, P = P, m_tot = m_NaCl + m_KCl)
+ # calculate logK of K+ + Cl- = KCl, adjusted for ionic strength
+ logKadj <- subcrt(c("K+", "Cl-", "KCl"), c(-1, -1, 1), T = T, P = P, IS = NaCl$IS)$out$logK
+ # what is the molality of K+ from 0.5 mol/kg KCl, assuming total chloride from above
+ m_K <- m_KCl / (10^logKadj * NaCl$m_Cl + 1)
+ list(IS = NaCl$IS, m_Cl = NaCl$m_Cl, m_K = m_K)
+}
+
# log(m_Au)-T diagram like Fig. 2B of Williams-Jones et al., 2009
# (doi:10.2113/gselements.5.5.281)
Au_T1 <- function() {
@@ -132,14 +146,13 @@
basis("H2S", "PPM")
# apply QMK buffer for pH
basis("H+", "QMK")
- basis("K+", log10(0.5))
- # calculate solution composition for 2 mol/kg NaCl
- NaCl <- NaCl(T = seq(150, 550, 10), P = 1000, m_tot=2)
+ # estimate solution composition for 1.5 m NaCl and 0.5 m KCl
+ chl <- chloride(T = seq(150, 550, 10), P = 1000, m_NaCl = 1.5, m_KCl = 0.5)
# calculate affinity and solubility
- a <- affinity(T = seq(150, 550, 10), `Cl-` = log10(NaCl$m_Cl), P = 1000, IS = NaCl$IS)
+ a <- affinity(T = seq(150, 550, 10), `Cl-` = log10(chl$m_Cl), `K+` = log10(chl$m_K), P = 1000, IS = chl$IS)
s <- solubility(a)
# make diagram and show total log molality
- diagram(s, ylim = c(-10, -4), col = col, lwd = 2, lty = 1)
+ diagram(s, ylim = c(-10, -3), col = col, lwd = 2, lty = 1)
diagram(s, add = TRUE, type = "loga.balance", lwd = 3, lty = 2)
# make legend and title
dP <- describe.property("P", 1000)
@@ -147,7 +160,7 @@
dK <- describe.basis(ibasis=5, use.molality=TRUE)
legend("topleft", c(dP, dNaCl, dK), bty = "n")
dbasis <- describe.basis(ibasis = c(9, 7, 10))
- legend(320, -4, dbasis, bty = "n")
+ legend(320, -3, dbasis, bty = "n")
title(main=("After Williams-Jones et al., 2009, Fig. 2B"), font.main = 1)
}
@@ -162,11 +175,13 @@
basis("O2", "HM")
# apply QMK buffer for pH
basis("H+", "QMK")
- basis("K+", log10(0.5))
- # calculate solution composition for 2 mol/kg NaCl
- NaCl <- NaCl(T = seq(150, 550, 10), P = 1000, m_tot=2)
- # calculate affinity and solubility
- a <- affinity(T = seq(150, 550, 10), `Cl-` = log10(NaCl$m_Cl), P = 1000, IS = NaCl$IS)
+ # estimate solution composition for 1.5 m NaCl and 0.5 m KCl
+ chl <- chloride(T = seq(150, 550, 10), P = 1000, m_NaCl = 1.5, m_KCl = 0.5)
+# # calculate affinity and solubility, considering speciation of sulfur
+# bases <- c("H2S", "HS-", "SO4-2", "HSO4-")
+# m <- mosaic(bases, T = seq(150, 550, 10), `Cl-` = log10(chl$m_Cl), `K+` = log10(chl$m_K), P = 1000, IS = chl$IS)
+# s <- solubility(m$A.species)
+ a <- affinity(T = seq(150, 550, 10), `Cl-` = log10(chl$m_Cl), `K+` = log10(chl$m_K), P = 1000, IS = chl$IS)
s <- solubility(a)
# make diagram and show total log molality
diagram(s, ylim = c(-10, -3), col = col, lwd = 2, lty = 1)
Modified: pkg/CHNOSZ/man/berman.Rd
===================================================================
--- pkg/CHNOSZ/man/berman.Rd 2018-11-08 13:59:57 UTC (rev 348)
+++ pkg/CHNOSZ/man/berman.Rd 2018-11-11 02:22:09 UTC (rev 349)
@@ -140,6 +140,8 @@
title(main = syslab(c("K2O", "Al2O3", "SiO2", "H2O", "HCl")), line = 1.8)
title(main = "Helgeson and Berman minerals, after Sverjensky et al., 1991",
line = 0.3, font.main = 1)
+# cleanup for next example
+data(thermo)
# make a P-T diagram for SiO2 minerals (Ber88 Fig. 4)
basis(c("SiO2", "O2"), c("cr", "gas"))
Modified: pkg/CHNOSZ/tests/testthat/test-swap.basis.R
===================================================================
--- pkg/CHNOSZ/tests/testthat/test-swap.basis.R 2018-11-08 13:59:57 UTC (rev 348)
+++ pkg/CHNOSZ/tests/testthat/test-swap.basis.R 2018-11-11 02:22:09 UTC (rev 349)
@@ -44,3 +44,13 @@
basis(names(bl99), bl99)
expect_equal(element.mu(T=99), ep99)
})
+
+# 20181111
+test_that("swapping works with a buffer (no recalculation of activities)", {
+ basis("FeCHNOS+")
+ oldb <- basis("O2", "PPM")
+ # before version 1.1.3-57, this gave Error in value * (-log(10) * R * T) : non-numeric argument to binary operator
+ newb <- swap.basis("O2", "hydrogen")
+ # note: logact includes "PPM" for O2 (old) and H2 (new)
+ expect_identical(oldb$logact, newb$logact)
+})
More information about the CHNOSZ-commits
mailing list