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

[noreply at r-forge.r-project.org]

Author: jedick
Date: 2018-01-28 08:16:39 +0100 (Sun, 28 Jan 2018)
New Revision: 302

Modified:
   pkg/CHNOSZ/DESCRIPTION
   pkg/CHNOSZ/R/equilibrate.R
   pkg/CHNOSZ/inst/NEWS
   pkg/CHNOSZ/tests/testthat/test-equilibrate.R
Log:
equilibrate(): accept 'normalize' argument with length > 1


Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION	2017-12-22 01:54:12 UTC (rev 301)
+++ pkg/CHNOSZ/DESCRIPTION	2018-01-28 07:16:39 UTC (rev 302)
@@ -1,6 +1,6 @@
-Date: 2017-12-21
+Date: 2018-01-28
 Package: CHNOSZ
-Version: 1.1.3-9
+Version: 1.1.3-10
 Title: Thermodynamic Calculations for Geobiochemistry
 Authors at R: c(
     person("Jeffrey", "Dick", , "j3ffdick at gmail.com", role = c("aut", "cre"),

Modified: pkg/CHNOSZ/R/equilibrate.R
===================================================================
--- pkg/CHNOSZ/R/equilibrate.R	2017-12-22 01:54:12 UTC (rev 301)
+++ pkg/CHNOSZ/R/equilibrate.R	2018-01-28 07:16:39 UTC (rev 302)
@@ -54,12 +54,14 @@
   ## normalize -- normalize the molar formula by the balance coefficients
   m.balance <- n.balance
   isprotein <- grepl("_", as.character(aout$species$name))
-  if(normalize | as.residue) {
+  if(any(normalize) | as.residue) {
     if(any(n.balance < 0)) stop("one or more negative balancing coefficients prohibit using normalized molar formulas")
-    n.balance <- rep(1, nspecies)
+    n.balance[normalize|as.residue] <- 1
     if(as.residue) message(paste("equilibrate: using 'as.residue' for molar formulas"))
     else message(paste("equilibrate: using 'normalize' for molar formulas"))
-  } else m.balance <- rep(1, nspecies)
+    # set the formula divisor (m.balance) to 1 for species whose formulas are *not* normalized
+    m.balance[!(normalize|as.residue)] <- 1
+  } else m.balance[] <- 1
   ## Astar: the affinities/2.303RT of formation reactions with
   ## formed species in their standard-state activities
   Astar <- lapply(1:nspecies, function(i) { 
@@ -76,7 +78,7 @@
   if(method[1]=="boltzmann") loga.equil <- equil.boltzmann(Astar, n.balance, loga.balance)
   else if(method[1]=="reaction") loga.equil <- equil.reaction(Astar, n.balance, loga.balance, tol)
   ## if we normalized the formulas, get back to activities to species
-  if(normalize & !as.residue) {
+  if(any(normalize) & !as.residue) {
     loga.equil <- lapply(1:nspecies, function(i) {
       loga.equil[[i]] - log10(m.balance[i])
     })

Modified: pkg/CHNOSZ/inst/NEWS
===================================================================
--- pkg/CHNOSZ/inst/NEWS	2017-12-22 01:54:12 UTC (rev 301)
+++ pkg/CHNOSZ/inst/NEWS	2018-01-28 07:16:39 UTC (rev 302)
@@ -1,5 +1,5 @@
-CHANGES IN CHNOSZ 1.1.3-8 (2017-12-21)
---------------------------------------
+CHANGES IN CHNOSZ 1.1.3-10 (2018-01-28)
+---------------------------------------
 
 - Lines in 1-D diagram()s can optionally be drawn as splines using the
   method for splinefun() given in the 'spline.method' argument (the
@@ -18,8 +18,13 @@
 - Fix a bug where subcrt()$reaction$coeffs was incorrect for reactions
   involving minerals with phase transitions.
 
-- Fix bug in nonideal() where 'Zn' in formula was identified as charge.
+- Fix bug in nonideal() where "Zn" in formula was identified as charge.
 
+- TODO: fix overly long message for info("SiO2").
+
+- In equilibrate(), accept a length > 1 'normalize' argument in order
+  normalize the chemical formulas of only the selected species.
+
 CHANGES IN CHNOSZ 1.1.3 (2017-11-13)
 ------------------------------------
 

Modified: pkg/CHNOSZ/tests/testthat/test-equilibrate.R
===================================================================
--- pkg/CHNOSZ/tests/testthat/test-equilibrate.R	2017-12-22 01:54:12 UTC (rev 301)
+++ pkg/CHNOSZ/tests/testthat/test-equilibrate.R	2018-01-28 07:16:39 UTC (rev 302)
@@ -155,6 +155,40 @@
   expect_equal(list2array(eO2$loga.equil)[21, ], e60.8)
 })
 
+test_that("normalizing formulas of only selected species works as expected", {
+  iC6 <- info("n-hexane", "liq")
+  iC12 <- info("n-dodecane", "liq")
+  `n-alkane` <- iC6:iC12
+  i2C6 <- info("2-methylpentane", "liq")
+  i2C9 <- info("2-methyloctane", "liq")
+  `2-isoalkane` <- i2C6:i2C9
+  basis("CHNOS")
+  basis("O2", -49.5)
+  species(`n-alkane`)
+  species(`2-isoalkane`)
+  # approximate conditions of Computer Experiment 27 (Helgeson et al., 2009, GCA)
+  a <- affinity(T=150, P=830, exceed.Ttr=TRUE)
+  # using full chemical formulas
+  efull <- equilibrate(a)
+  dloga_isoalkane_full <- diff(unlist(efull$loga.equil[c(8, 11)]))
+  # normalize all the formulas
+  enorm <- equilibrate(a, normalize=TRUE)
+  dloga_nalkane_norm <- diff(unlist(enorm$loga.equil[c(1, 7)]))
+  dloga_isoalkane_norm <- diff(unlist(enorm$loga.equil[c(8, 11)]))
+  # normalize only the n-alkane formulas
+  isalk <- grepl("n-", species()$name)
+  emix <- equilibrate(a, normalize=isalk)
+  # the activity ratios for the normalized formulas should be the same in both calculations
+  dloga_nalkane_mix <- diff(unlist(emix$loga.equil[c(1, 7)]))
+  expect_equal(dloga_nalkane_mix, dloga_nalkane_norm)
+  # the actvitity ratios for the not-normalized formulas should be similar in both calculations
+  # (not identical becuase they are affected by total activity, unlike normalized formulas)
+  dloga_isoalkane_mix <- diff(unlist(emix$loga.equil[c(8, 11)]))
+  expect_maxdiff(dloga_isoalkane_mix, dloga_isoalkane_full, 0.07)
+  # however, the difference between normalized and not-normalized formulas is much greater
+  expect_true(maxdiff(dloga_isoalkane_mix, dloga_isoalkane_norm) > maxdiff(dloga_isoalkane_mix, dloga_isoalkane_full))
+})
+
 # references
 
 # Seewald, J. S. (2001)