[CHNOSZ-commits] r360 - in pkg/CHNOSZ: . R inst man tests/testthat

Mon Jan 21 04:48:18 CET 2019

Author: jedick
Date: 2019-01-21 04:48:17 +0100 (Mon, 21 Jan 2019)
New Revision: 360

Modified:
   pkg/CHNOSZ/DESCRIPTION
   pkg/CHNOSZ/R/mosaic.R
   pkg/CHNOSZ/inst/NEWS
   pkg/CHNOSZ/man/solubility.Rd
   pkg/CHNOSZ/tests/testthat/test-mosaic.R
Log:
mosaic(): with blend=TRUE, include Gibbs energy of mixing


Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================

--- pkg/CHNOSZ/DESCRIPTION	2019-01-20 18:53:02 UTC (rev 359)
+++ pkg/CHNOSZ/DESCRIPTION	2019-01-21 03:48:17 UTC (rev 360)
@@ -1,6 +1,6 @@
-Date: 2019-01-20
+Date: 2019-01-21
 Package: CHNOSZ
-Version: 1.1.3-67
+Version: 1.1.3-68
 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/mosaic.R
===================================================================
--- pkg/CHNOSZ/R/mosaic.R	2019-01-20 18:53:02 UTC (rev 359)
+++ pkg/CHNOSZ/R/mosaic.R	2019-01-21 03:48:17 UTC (rev 360)
@@ -93,8 +93,10 @@
         for(i in seq_along(A.species$values)) {
           # start with zero affinity
           if(j==1) A.species$values[[i]][] <- 0
-          # add affinity scaled by __relative__ abundance of this basis species
-          A.species$values[[i]] <- A.species$values[[i]] + affs[[j]]$values[[i]] * 10^e$loga.equil[[j]]/a.tot
+          # add affinity scaled by relative abundance of this basis species
+          # and include mixing term (-x*log10(x)) 20190121
+          x <- 10^e$loga.equil[[j]]/a.tot
+          A.species$values[[i]] <- A.species$values[[i]] + affs[[j]]$values[[i]] * x - x * log10(x)
         }
       }
     }

Modified: pkg/CHNOSZ/inst/NEWS
===================================================================
--- pkg/CHNOSZ/inst/NEWS	2019-01-20 18:53:02 UTC (rev 359)
+++ pkg/CHNOSZ/inst/NEWS	2019-01-21 03:48:17 UTC (rev 360)
@@ -1,6 +1,25 @@
-CHANGES IN CHNOSZ 1.1.3-67 (2019-01-20)
+CHANGES IN CHNOSZ 1.1.3-68 (2019-01-21)
 ---------------------------------------
 
+BUG FIXES
+
+- Fix a bug where subcrt()$reaction$coeffs was incorrect for reactions
+  involving minerals with phase transitions. Also ensure that the output
+  reaction stoichiometry is correct for duplicated species in reactions.
+  Thanks to Grayson Boyer for the bug report.
+
+- Fix bug in nonideal() where "Zn" in formula was identified as charge.
+  Thanks to Feng Lai for reporting the incorrect behavior caused by
+  this bug.
+
+- For species in the revised HKF model, subcrt() now sets properties to
+  NA where the density of H2O is less than 0.35 g/cm3, avoiding the
+  output of bogus values in this region. Thanks to Evgeniy Bastrakov.
+
+- In mosaic() with blend = TRUE, add a previously missing term for
+  Gibbs energy of mixing. An example using mosaic() to calculate the
+  pH-dependent solubility of calcite has been added to solubility.Rd.
+
 NEW FEATURE: SOLUBILITY CALCULATIONS
 
 - Add solubility(). Run this after affinity() to calculate the
@@ -147,21 +166,6 @@
   coordinates of lines (field boundaries) on 2-D diagrams (these are
   taken from the output of contourLines()).
 
-BUG FIXES
-
-- Fix a bug where subcrt()$reaction$coeffs was incorrect for reactions
-  involving minerals with phase transitions. Also ensure that the output
-  reaction stoichiometry is correct for duplicated species in reactions.
-  Thanks to Grayson Boyer for the bug report.
-
-- Fix bug in nonideal() where "Zn" in formula was identified as charge.
-  Thanks to Feng Lai for reporting the incorrect behavior caused by
-  this bug.
-
-- For species in the revised HKF model, subcrt() now sets properties to
-  NA where the density of H2O is less than 0.35 g/cm3, avoiding the
-  output of bogus values in this region. Thanks to Evgeniy Bastrakov.
-
 COMPUTATIONAL OPTIONS
 
 - Add 'exceed.rhomin' argument to subcrt() and affinity() to enable

Modified: pkg/CHNOSZ/man/solubility.Rd
===================================================================
--- pkg/CHNOSZ/man/solubility.Rd	2019-01-20 18:53:02 UTC (rev 359)
+++ pkg/CHNOSZ/man/solubility.Rd	2019-01-21 03:48:17 UTC (rev 360)
@@ -64,16 +64,14 @@
 }
 
 \examples{\dontshow{data(thermo)}
-# solubility of CO2 and calcite as a function of pH
+## solubility of CO2 and calcite as a function of pH
 opar <- par(mfrow = c(1, 2))
-
-# set pH range and resolution, constant temperature and ionic strength
+## set pH range and resolution, constant temperature and ionic strength
 pH <- c(0, 14)
 res <- 100
 T <- 25
 IS <- 0
-
-# start with CO2 (not a dissociation reaction)
+## start with CO2 (not a dissociation reaction)
 basis(c("carbon dioxide", "H2O", "O2", "H+"))
 # ca. atmospheric PCO2
 basis("CO2", -3.5)
@@ -92,8 +90,7 @@
 title(main = substitute("Solubility of"~what~"at"~T~degree*"C",
   list(what = expr.species("CO2"), T = T)), line = 1.5)
 mtext("cf. Fig. 4.5 of Stumm and Morgan, 1996")
-
-# now do calcite (a dissociation reaction)
+## now do calcite (a dissociation reaction)
 basis(c("calcite", "Ca+2", "H2O", "O2", "H+"))
 species(c("CO2", "HCO3-", "CO3-2"))
 a <- affinity(pH = c(pH, res), T = T, IS = IS)
@@ -105,8 +102,38 @@
 title(main = substitute("Solubility of"~what~"at"~T~degree*"C",
   list(what = "calcite", T = T)))
 mtext("cf. Fig. 4A of Manning et al., 2013")
+par(opar)
 
-par(opar)
+## two ways to calculate pH-dependent solubility of calcite
+## with ionic strength determination
+## method 1: CO2 and carbonate species as formed species
+basis(c("calcite", "Ca+2", "H2O", "O2", "H+"))
+species(c("CO2", "HCO3-", "CO3-2"))
+# ionic strength calculations don't converge below around pH=3
+a <- affinity(pH = c(3, 14))
+sa0 <- solubility(a)
+saI <- solubility(a, find.IS = TRUE)
+## method 2: CO2 and carbonate species as basis species
+basis(c("calcite", "CO2", "H2O", "O2", "H+"))
+species(c("Ca+2"))
+m <- mosaic(c("CO2", "HCO3-", "CO3-2"), pH = c(3, 14), blend = TRUE)
+sm0 <- solubility(m)
+smI <- solubility(m, find.IS = TRUE)
+## plot the results
+plot(0, 0, xlab="pH", ylab="solubility, log mol", xlim = c(3, 14), ylim = c(-5, 2))
+# method 1 with/without ionic strength
+lines(a$vals[[1]], saI$loga.balance, lwd=5, col="lightblue")
+lines(a$vals[[1]], sa0$loga.balance, lwd=5, col="pink")
+# method 2 with/without ionic strength
+lines(a$vals[[1]], smI$loga.balance, lty=2)
+lines(a$vals[[1]], sm0$loga.balance, lty=2)
+legend("topright", c("I = 0", "I = calculated", "mosaic method"),
+       col = c("pink", "lightblue", "black"), lwd = c(5, 5, 1), lty = c(1, 1, 2))
+# the two methods give nearly equivalent results
+stopifnot(all.equal(sa0$loga.balance, sm0$loga.balance))
+stopifnot(all.equal(saI$loga.balance, smI$loga.balance, tolerance = 0.003))
+## NOTE: the second method (using mosaic) takes longer, but is
+## more flexible; complexes with Ca+2 could be included
 }
 
 \references{

Modified: pkg/CHNOSZ/tests/testthat/test-mosaic.R
===================================================================
--- pkg/CHNOSZ/tests/testthat/test-mosaic.R	2019-01-20 18:53:02 UTC (rev 359)
+++ pkg/CHNOSZ/tests/testthat/test-mosaic.R	2019-01-21 03:48:17 UTC (rev 360)
@@ -36,10 +36,13 @@
   pH <- c(0, 14, 29)
   m1 <- mosaic(bases, pH=pH)
   m2 <- mosaic(bases, pH=pH, blend=TRUE)
-  # these species have no S so the results should be the same
-  expect_equal(m1$A.species$values, m2$A.species$values)
+  # these species have no S so the results should be similar,
+  # 20190121 except for a negative free energy of mixing (positive affinity)
+  expect_true(all(m2$A.species$values[[1]] - m1$A.species$values[[1]] > 0))
+  # the differences increase, then decrease
+  expect_equal(unique(sign(diff(as.numeric(m2$A.species$values[[1]] - m1$A.species$values[[1]])))), c(1, -1))
+  # now with S-bearing species ...
   species(c("pyrrhotite", "pyrite"))
-  # now with S-bearing species ...
   m3 <- mosaic(bases, pH=pH)
   m4 <- mosaic(bases, pH=pH, blend=TRUE)
   # the results are different ...

