[CHNOSZ-commits] r358 - in pkg/CHNOSZ: . R inst

noreply at r-forge.r-project.org noreply at r-forge.r-project.org
Sun Jan 20 09:24:28 CET 2019


Author: jedick
Date: 2019-01-20 09:24:26 +0100 (Sun, 20 Jan 2019)
New Revision: 358

Modified:
   pkg/CHNOSZ/DESCRIPTION
   pkg/CHNOSZ/R/solubility.R
   pkg/CHNOSZ/inst/NEWS
Log:
solubility(): calculate ionic strength for ions in basis and species


Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION	2019-01-17 15:40:17 UTC (rev 357)
+++ pkg/CHNOSZ/DESCRIPTION	2019-01-20 08:24:26 UTC (rev 358)
@@ -1,6 +1,6 @@
-Date: 2019-01-17
+Date: 2019-01-20
 Package: CHNOSZ
-Version: 1.1.3-65
+Version: 1.1.3-66
 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/solubility.R
===================================================================
--- pkg/CHNOSZ/R/solubility.R	2019-01-17 15:40:17 UTC (rev 357)
+++ pkg/CHNOSZ/R/solubility.R	2019-01-20 08:24:26 UTC (rev 358)
@@ -38,7 +38,7 @@
   bout <- balance(aout)
   n.balance <- bout$n.balance
   balance <- bout$balance
-  # get logarithm of total activity of the balancing basis species
+  # get logarithm of total activity of the conserved basis species
   logabfun <- function(loga.equil, n.balance) {
     # exponentiate, multiply by n.balance, sum, logarithm
     a.equil <- mapply("^", 10, loga.equil, SIMPLIFY = FALSE)
@@ -79,16 +79,36 @@
 
     # get the old IS
     IS.old <- rep(IS, length.out = length(aout$values[[1]]))
-    # calculate the ionic strength (assuming no ion pairing)
-    # i.e. for Sr+2 and SO4-2
-    # TODO: determine the correct values for the actual species in the system
-    mol.cation <- mol.anion <- 10^loga.balance
-    Z <- 2
-    IS <- (mol.cation * Z^2 + mol.anion * Z^2) / 2
+    # calculate the ionic strength for unpaired ions: IS = sum(molality * Z^2) / 2
+    sum.mZ2 <- 0
+    ion.names <- character()
+    # is there an ion in the basis species?
+    if(dissociation) {
+      basis.ion <- rownames(aout$basis)[2]
+      Z.basis.ion <- makeup(basis.ion)["Z"]
+      if(!is.na(Z.basis.ion)) {
+        sum.mZ2 <- sum.mZ2 + 10^loga.balance * Z.basis.ion^2
+        ion.names <- c(ion.names, basis.ion)
+      }
+    }
+    # add ions present in the species of interest
+    for(i in 1:length(loga.equil)) {
+      species.ion <- aout$species$name[i]
+      Z.species.ion <- makeup(species.ion)["Z"]
+      if(!is.na(Z.species.ion)) {
+        sum.mZ2 <- sum.mZ2 + 10^loga.equil[[i]] * Z.species.ion^2
+        ion.names <- c(ion.names, species.ion)
+      }
+    }
+    IS <- sum.mZ2 / 2
+    # report ions used in the ionic strength calculation
+    if(find.IS & niter==1) message("solubility: ionic strength calculated for ", paste(ion.names, collapse=" "))
     # report the current ionic strength
     if(find.IS) message("solubility: (iteration ", niter, ") ionic strength range is ", paste(round(range(IS), 4), collapse=" "))
     # stop iterating if we reached the tolerance (or find.IS=FALSE)
     if(!find.IS | all(IS - IS.old < 1e-4)) break
+    # expand argument values for affinity()
+    for(i in 1:length(aout$vals)) aout$args[[i]] <- aout$vals[[i]]
     # recalculate the affinity using the new IS
     aout <- suppressMessages(do.call(aout$fun, list(aout, IS = IS)))
     niter <- niter + 1

Modified: pkg/CHNOSZ/inst/NEWS
===================================================================
--- pkg/CHNOSZ/inst/NEWS	2019-01-17 15:40:17 UTC (rev 357)
+++ pkg/CHNOSZ/inst/NEWS	2019-01-20 08:24:26 UTC (rev 358)
@@ -1,23 +1,24 @@
-CHANGES IN CHNOSZ 1.1.3-65 (2019-01-17)
+CHANGES IN CHNOSZ 1.1.3-66 (2019-01-20)
 ---------------------------------------
 
-NEW FEATURES
+NEW FEATURE: SOLUBILITY CALCULATIONS
 
 - Add solubility(). Run this after affinity() to calculate the
   solubility of a solid or gas defined as the conserved basis species,
   which is involved in the formation of one or more dissolved species.
-  Features in development include automatic detection of dissociation
-  reactions and finding the final ionic strength for dissolution of a
-  mineral into pure water.
 
+- Features include automatic detection of dissociation reactions and
+  finding the final ionic strength for dissolution of a mineral into
+  pure water (find.IS argument).
+
+- find.IS depends on the new argument recall feature of affinity().
+  This allows a calculation to be re-run with the same settings except
+  for particular additions or modifications.
+
 - Revise demo/solubility.R to show solubility calculations for CO2(gas)
   and calcite as a function of T and pH.
 
-- Add demo/gold.R for calculations of Au solubility in hydrothermal
-  chloride and sulfide solutions (based on diagrams from Akinfiev and
-  Zotov, 2001, Stefánsson and Seward, 2004, and Williams-Jones et al.,
-  2009). This depends on the revised nonideal() and new NaCl() functions
-  described next.
+NEW FEATURE: EXPANDED ACTIVITY COEFFICIENT CALCULATIONS
 
 - Reorganize and expand options for activity coefficient calculations
   (set in thermo$opt$nonideal: Bdot, Bdot0, bgamma, bgamma0, or Alberty).
@@ -41,17 +42,26 @@
   of NaCl in water, taking account of activity coefficients and the
   reaction Na+ + Cl- = NaCl(aq).
 
-- Add dumpdata() for returning/writing all packaged thermodynamic data
-  (including default database and optional data files).
+DOCUMENTATION
 
+- Add demo/gold.R for calculations of Au solubility in hydrothermal
+  chloride and sulfide solutions (based on diagrams from Akinfiev and
+  Zotov, 2001, Stefánsson and Seward, 2004, and Williams-Jones et al.,
+  2009). This depends on the revised nonideal() and new NaCl() functions
+  described above.
+
+- anintro.Rmd: add cuprite to mosaic diagram example, and note about
+  implications of changing balance to 1.
+
 - Add demo/bison.R (average oxidation state of carbon of metagenome-
   derived proteins in different microbial phyla at Bison Pool)
 
-- affinity(): implement argument recall, to re-run a calculation with
-  the same settings except for particular additions or modifications.
-
 THERMODYNAMIC DATA
 
+- Add dumpdata() for returning/writing all packaged thermodynamic data
+  (including default database and optional data files). The file is
+  also available on the website (chnosz.net/download/alldata.csv).
+
 - The Berman data (Berman, 1988 and later additions) have replaced the
   SUPCRT92 data (based on Helgeson et al., 1978) for most minerals in
   the default database (i.e. the one loaded by data(thermo)). Only
@@ -136,11 +146,6 @@
   coordinates of lines (field boundaries) on 2-D diagrams (these are
   taken from the output of contourLines()).
 
-DOCUMENTATION
-
-- anintro.Rmd: add cuprite to mosaic diagram example, and note about
-  implications of changing balance to 1.
-
 BUG FIXES
 
 - Fix a bug where subcrt()$reaction$coeffs was incorrect for reactions
@@ -178,7 +183,7 @@
 
 USABILITY ENHANCEMENTS
 
-- To provide better diagnostics for potential web apps, warning
+- To provide better diagnostics for other apps using CHNOSZ, warning
   messages produced by subcrt() are now available in the output of
   affinity(), under 'sout$warnings'.
 



More information about the CHNOSZ-commits mailing list