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

noreply at r-forge.r-project.org noreply at r-forge.r-project.org
Mon Dec 19 05:17:19 CET 2016 

Author: jedick
Date: 2016-12-19 05:17:18 +0100 (Mon, 19 Dec 2016)
New Revision: 117

Added:
   pkg/CHNOSZ/tests/testthat/test-nonideal.R
Modified:
   pkg/CHNOSZ/DESCRIPTION
   pkg/CHNOSZ/R/nonideal.R
   pkg/CHNOSZ/R/subcrt.R
   pkg/CHNOSZ/demo/solubility.R
   pkg/CHNOSZ/inst/NEWS
   pkg/CHNOSZ/man/data.Rd
   pkg/CHNOSZ/man/nonideal.Rd
   pkg/CHNOSZ/man/util.units.Rd
   pkg/CHNOSZ/man/water.Rd
Log:
subcrt() uses common logarithm for loggam


Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION	2016-11-29 02:50:22 UTC (rev 116)
+++ pkg/CHNOSZ/DESCRIPTION	2016-12-19 04:17:18 UTC (rev 117)
@@ -1,6 +1,6 @@
-Date: 2016-11-29
+Date: 2016-12-19
 Package: CHNOSZ
-Version: 1.0.8-4
+Version: 1.0.8-5
 Title: Chemical Thermodynamics and Activity Diagrams
 Author: Jeffrey Dick
 Maintainer: Jeffrey Dick <j3ffdick at gmail.com>

Modified: pkg/CHNOSZ/R/nonideal.R
===================================================================
--- pkg/CHNOSZ/R/nonideal.R	2016-11-29 02:50:22 UTC (rev 116)
+++ pkg/CHNOSZ/R/nonideal.R	2016-12-19 04:17:18 UTC (rev 117)
@@ -25,7 +25,7 @@
     }
     # Alberty, 2003 Eq. 3.6-1
     loggamma <- function(a,Z,I,B) { - a * Z^2 * I^(1/2) / (1 + B * I^(1/2)) }
-    # XXX check the following equations 20080208 jmd
+    # TODO: check the following equations 20080208 jmd
     if(prop=='log') return(loggamma(eval(A),Z,I,B))
     else if(prop=='G') return(thermo$opt$R * T * loggamma(eval(A),Z,I,B))
     else if(prop=='H') return(thermo$opt$R * T^2 * loggamma(eval(DD(A,'T',1)),Z,I,B))

Modified: pkg/CHNOSZ/R/subcrt.R
===================================================================
--- pkg/CHNOSZ/R/subcrt.R	2016-11-29 02:50:22 UTC (rev 116)
+++ pkg/CHNOSZ/R/subcrt.R	2016-12-19 04:17:18 UTC (rev 117)
@@ -392,14 +392,10 @@
   isaq.new <- logical()
   iscgl.new <- logical()
   isH2O.new <- logical()
-  #print(sinfo)
-  #print(sinph)
-  #print(reaction)
   for(i in 1:length(sinfo)) {
     iphases <- which(sinfo[i]==sinph)
     # deal with repeated species here ... divide iphases 
     # by the number of duplicates
-    #print(iphases)
     if(TRUE %in% duplicated(inpho[iphases])) {
       iphases <- iphases[length(which(sinfo==sinfo[i]))]
     }
@@ -467,11 +463,8 @@
   iscgl <- iscgl.new
   isH2O <- isH2O.new
 
-  #print(out)
-
   newprop <- eprop[eprop!='rho']
   # the order of the properties
-  #if(ncol(out[[1]])>1) for(i in 1:length(out)) {
   if(length(newprop)>1) for(i in 1:length(out)) {
     # keep state/loggam columns if they exists
     ipp <- match(newprop,tolower(colnames(out[[i]])))
@@ -548,8 +541,14 @@
       }
     }
   }
+  # convert loggam to common logarithm and
   # put ionic strength next to any loggam columns
-  for(i in 2:length(out)) if('loggam' %in% colnames(out[[i]])) out[[i]] <- cbind(out[[i]],IS=newIS)
+  for(i in 2:length(out)) {
+    if('loggam' %in% colnames(out[[i]])) {
+      out[[i]] <- cbind(out[[i]],IS=newIS)
+      out[[i]][, "loggam"] <- out[[i]][, "loggam"]/log(10)
+    }
+  }
   # more fanagling for species
   if(!do.reaction) {
     out <- list(species=out$species,out=out[2:length(out)])

Modified: pkg/CHNOSZ/demo/solubility.R
===================================================================
--- pkg/CHNOSZ/demo/solubility.R	2016-11-29 02:50:22 UTC (rev 116)
+++ pkg/CHNOSZ/demo/solubility.R	2016-12-19 04:17:18 UTC (rev 117)
@@ -31,7 +31,7 @@
   species(1:3, unlist(e$loga.equil))
   a <- affinity(T=T)
   # check they're actually equal
-  stopifnot(all(abs(unlist(a$values) - a$values[[1]]) < 1e-10))
+  stopifnot(all(abs(unlist(a$values) - as.vector(a$values[[1]])) < 1e-10))
   return(a$values[[1]])
 }
 

Modified: pkg/CHNOSZ/inst/NEWS
===================================================================
--- pkg/CHNOSZ/inst/NEWS	2016-11-29 02:50:22 UTC (rev 116)
+++ pkg/CHNOSZ/inst/NEWS	2016-12-19 04:17:18 UTC (rev 117)
@@ -1,4 +1,4 @@
-CHANGES IN CHNOSZ 1.0.8-4 (2016-11-29)
+CHANGES IN CHNOSZ 1.0.8-5 (2016-12-19)
 --------------------------------------
 
 - Add "AA" as a keyword for preset species in basis() (cysteine,
@@ -9,6 +9,12 @@
 - More flexible parsing of chemical formulas for ZC() and other
   functions; e.g. ZC(colMeans(protein.formula(1:4))) now works.
 
+- supcrt() now returns 'loggam' using the common logarithm; add
+  test-nonideal.R to check for consistency between loggam and logK
+  values returned by supcrt(). Thanks to David T. Wang for the bug
+  report and test. (This also fixes the bug previously noted for the
+  first example in ?nonideal.)
+
 CHANGES IN CHNOSZ 1.0.8 (2016-05-28)
 ------------------------------------
 

Modified: pkg/CHNOSZ/man/data.Rd
===================================================================
--- pkg/CHNOSZ/man/data.Rd	2016-11-29 02:50:22 UTC (rev 116)
+++ pkg/CHNOSZ/man/data.Rd	2016-12-19 04:17:18 UTC (rev 117)
@@ -246,7 +246,7 @@
 
   Shock, E. L. et al. {1998} \emph{SLOP98.dat} (computer data file). http://geopig.asu.edu/supcrt92_data/slop98.dat, accessed on 2005-11-05. Current location: \url{http://geopig.asu.edu/sites/default/files/slop98.dat}.
 
-  Wagman, D. D., Evans, W. H., Parker, V. B., Schumm, R. H., Halow, I., Bailey, S. M., Churney, K. L. and Nuttall, R. L. (1982) The NBS tables of chemical thermodynamic properties. Selected values for inorganic and C\eqn{_1}{1} and C\eqn{_2}{2} organic substances in SI units. \emph{J. Phys. Chem. Ref. Data} \bold{11} (supp. 2), 1--392. \url{http://www.nist.gov/data/PDFfiles/jpcrdS2Vol11.pdf}
+  Wagman, D. D., Evans, W. H., Parker, V. B., Schumm, R. H., Halow, I., Bailey, S. M., Churney, K. L. and Nuttall, R. L. (1982) The NBS tables of chemical thermodynamic properties. Selected values for inorganic and C\eqn{_1}{1} and C\eqn{_2}{2} organic substances in SI units. \emph{J. Phys. Chem. Ref. Data} \bold{11} (supp. 2), 1--392. \url{https://srd.nist.gov/JPCRD/jpcrdS2Vol11.pdf}
 
 }
 

Modified: pkg/CHNOSZ/man/nonideal.Rd
===================================================================
--- pkg/CHNOSZ/man/nonideal.Rd	2016-11-29 02:50:22 UTC (rev 116)
+++ pkg/CHNOSZ/man/nonideal.Rd	2016-12-19 04:17:18 UTC (rev 117)
@@ -26,8 +26,9 @@
 }
 
 \section{Warning}{
-The activity coefficients (gamma) in the first example below are lower than those shown on p. 274-276 of Alberty, 2003 (compare gamma near 0.75 in for charge = -1, IS = 0.25, T = 25 on p. 274 with gamma near 0.5 here).
-This may be due to an unidentified bug in the function.
+The logarithms of activity coefficients (\code{loggam}) returned by \code{nonideal} use the natural logarithm (cf. Alberty, 2003 Eq. 3.6-1).
+To maintain consistency with the conventions used elsewhere in the package (i.e. for logarithms of equilibrium constants and of chemical activities), the values of \code{loggam} returned by \code{\link{subcrt}} are expressed using the common (base 10) logarithm.
+Note that the first example below uses \code{loggam} returned by \code{subcrt}, therefore requiring a base of 10 for calculating gamma.
 }
 
 \examples{\dontshow{data(thermo)}
@@ -43,9 +44,10 @@
   s <- subcrt(c("H2PO4-", "HADP-2", "HATP-3", "ATP-4"), IS=IS, grid="IS", T=T[j])
   for(i in 1:4) lines(IS, 10^s$out[[i]]$loggam, col=col[j])
 }
-text(0.1, 0.62, "Z = -1")
-text(0.075, 0.18, "Z = -2")
-text(0.05, 0.06, "Z = -3")
+text(0.125, 0.8, "Z = -1")
+text(0.1, 0.42, "Z = -2")
+text(0.075, 0.18, "Z = -3")
+text(0.05, 0.08, "Z = -4")
 title(main=paste("activity coefficient (gamma) of -1,-2,-3,-4",
   "charged species at 0, 25, 40 deg C, after Alberty, 2003",
   sep="\n"), cex.main=0.95)

Modified: pkg/CHNOSZ/man/util.units.Rd
===================================================================
--- pkg/CHNOSZ/man/util.units.Rd	2016-11-29 02:50:22 UTC (rev 116)
+++ pkg/CHNOSZ/man/util.units.Rd	2016-12-19 04:17:18 UTC (rev 117)
@@ -33,7 +33,9 @@
 
 \details{
 
-   The units settings are used by \code{\link{subcrt}}, \code{\link{affinity}}, and \code{\link{diagram}} to accept input in or convert output to the units desired by the user. The settings, which can be queried or changed with \code{T.units}, \code{E.units} and \code{P.units}, refer to the units of temperature (\code{K} or \code{C}), energy (\code{cal} or \code{J}), and pressure (\code{bar}, \code{MPa}). (The first value in each of those pairs refers to the default units).
+The units settings are used by \code{\link{subcrt}}, \code{\link{affinity}}, and \code{\link{diagram}} to accept input in or convert output to the units desired by the user.
+The settings, which can be queried or changed with \code{T.units}, \code{E.units} and \code{P.units}, refer to the units of temperature (\code{C} or \code{K}), energy (\code{cal} or \code{J}), and pressure (\code{bar}, \code{MPa}).
+(The first value in each of those pairs refers to the default units).
 
 The actual units conversions are handled by \code{convert}, through which \code{values} are transformed into destination \code{units} (names not case sensitive).
 The possible conversions and settings for the \code{units} argument are shown in the following table.

Modified: pkg/CHNOSZ/man/water.Rd
===================================================================
--- pkg/CHNOSZ/man/water.Rd	2016-11-29 02:50:22 UTC (rev 116)
+++ pkg/CHNOSZ/man/water.Rd	2016-12-19 04:17:18 UTC (rev 117)
@@ -142,7 +142,7 @@
 
 \references{
 
-   Archer, D. G. and Wang, P. M. (1990) The dielectric constant of water and Debye-Huckel limiting law slopes. \emph{J. Phys. Chem. Ref. Data} \bold{19}, 371--411. \url{http://www.nist.gov/data/PDFfiles/jpcrd383.pdf}
+   Archer, D. G. and Wang, P. M. (1990) The dielectric constant of water and Debye-Huckel limiting law slopes. \emph{J. Phys. Chem. Ref. Data} \bold{19}, 371--411. \url{http://dx.doi.org/10.1063/1.555853}
 
   Haar, L., Gallagher, J. S. and Kell, G. S. (1984) \emph{NBS/NRC Steam Tables}. Hemisphere, Washington, D. C., 320 p. \url{http://www.worldcat.org/oclc/301304139}
 
@@ -152,7 +152,7 @@
 
   Johnson, J. W., Oelkers, E. H. and Helgeson, H. C. (1992) SUPCRT92: A software package for calculating the standard molal thermodynamic properties of minerals, gases, aqueous species, and reactions from 1 to 5000 bar and 0 to 1000\eqn{^{\circ}}{°}C. \emph{Comp. Geosci.} \bold{18}, 899--947. \url{http://dx.doi.org/10.1016/0098-3004(92)90029-Q}
 
-  Levelt-Sengers, J. M. H., Kamgarparsi, B., Balfour, F. W. and Sengers, J. V. (1983) Thermodynamic properties of steam in the critical region. \emph{J. Phys. Chem. Ref. Data} \bold{12}, 1--28. \url{http://www.nist.gov/data/PDFfiles/jpcrd214.pdf}
+  Levelt-Sengers, J. M. H., Kamgarparsi, B., Balfour, F. W. and Sengers, J. V. (1983) Thermodynamic properties of steam in the critical region. \emph{J. Phys. Chem. Ref. Data} \bold{12}, 1--28. \url{http://dx.doi.org/10.1063/1.555676}
 
   Wagner, W. and Pruss, A. (2002) The IAPWS formulation 1995 for the thermodynamic properties of ordinary water substance for general and scientific use. \emph{J. Phys. Chem. Ref. Data} \bold{31}, 387--535. \url{http://dx.doi.org/10.1063/1.1461829}
 

Added: pkg/CHNOSZ/tests/testthat/test-nonideal.R
===================================================================
--- pkg/CHNOSZ/tests/testthat/test-nonideal.R	                        (rev 0)
+++ pkg/CHNOSZ/tests/testthat/test-nonideal.R	2016-12-19 04:17:18 UTC (rev 117)
@@ -0,0 +1,7 @@
+context("nonideal")
+
+test_that("loggam and logK values are consistent", {
+  rxn1 <- subcrt(c("anhydrite", "Ca+2", "SO4-2"), c(-1, 1, 1), P=1, T=25, I=0)
+  rxn2 <- subcrt(c("anhydrite", "Ca+2", "SO4-2"), c(-1, 1, 1), P=1, T=25, I=0.7)
+  expect_equal(rxn1$out$logK, rxn2$out$loggam + rxn2$out$logK)
+})

More information about the CHNOSZ-commits mailing list