[CHNOSZ-commits] r256 - in pkg/CHNOSZ: . R demo inst man

noreply at r-forge.r-project.org noreply at r-forge.r-project.org
Fri Oct 13 09:31:39 CEST 2017 

Author: jedick
Date: 2017-10-13 09:31:38 +0200 (Fri, 13 Oct 2017)
New Revision: 256

Modified:
   pkg/CHNOSZ/DESCRIPTION
   pkg/CHNOSZ/R/affinity.R
   pkg/CHNOSZ/R/nonideal.R
   pkg/CHNOSZ/demo/DEW.R
   pkg/CHNOSZ/inst/NEWS
   pkg/CHNOSZ/man/nonideal.Rd
Log:
affinity(): variables in arguments override buffers in thermo$basis


Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION	2017-10-13 05:25:56 UTC (rev 255)
+++ pkg/CHNOSZ/DESCRIPTION	2017-10-13 07:31:38 UTC (rev 256)
@@ -1,6 +1,6 @@
 Date: 2017-10-13
 Package: CHNOSZ
-Version: 1.1.0-54
+Version: 1.1.0-55
 Title: Thermodynamic Calculations for Geobiochemistry
 Author: Jeffrey Dick
 Maintainer: Jeffrey Dick <j3ffdick at gmail.com>

Modified: pkg/CHNOSZ/R/affinity.R
===================================================================
--- pkg/CHNOSZ/R/affinity.R	2017-10-13 05:25:56 UTC (rev 255)
+++ pkg/CHNOSZ/R/affinity.R	2017-10-13 07:31:38 UTC (rev 256)
@@ -7,6 +7,10 @@
 #source("util.character.R")
 #source("util.list.R")
 #source("subcrt.R")
+#source("buffer.R")
+#source("util.args.R")
+#source("util.data.R")
+#source("species.R")
 
 affinity <- function(...,property=NULL,sout=NULL,exceed.Ttr=FALSE,
   return.buffer=FALSE,balance="PBB",iprotein=NULL,loga.protein=-3) {
@@ -70,7 +74,10 @@
 
     # buffer stuff
     buffer <- FALSE
-    ibufbasis <- which(!can.be.numeric(mybasis$logact))
+    hasbuffer <- !can.be.numeric(mybasis$logact)
+    # however, variables specified in the arguments take precedence over the buffer 20171013
+    isnotvar <- !rownames(mybasis) %in% args$vars
+    ibufbasis <- which(hasbuffer & isnotvar)
     if(!is.null(mybasis) & length(ibufbasis) > 0) {
       buffer <- TRUE
       message('affinity: loading buffer species')

Modified: pkg/CHNOSZ/R/nonideal.R
===================================================================
--- pkg/CHNOSZ/R/nonideal.R	2017-10-13 05:25:56 UTC (rev 255)
+++ pkg/CHNOSZ/R/nonideal.R	2017-10-13 07:31:38 UTC (rev 256)
@@ -62,7 +62,7 @@
     R <- 1.9872  # gas constant, cal K^-1 mol^-1
     if(prop=="loggamma") return(loggamma)
     else if(prop=="G") return(R * T * log(10) * loggamma)
-    # note the log(10) (=2.303) ... use natural logarithm to calculate G!!!
+    # note the log(10) (=2.303) ... use natural logarithm to calculate G
   }
 
   # get B-dot if we're using the Helgeson method

Modified: pkg/CHNOSZ/demo/DEW.R
===================================================================
--- pkg/CHNOSZ/demo/DEW.R	2017-10-13 05:25:56 UTC (rev 255)
+++ pkg/CHNOSZ/demo/DEW.R	2017-10-13 07:31:38 UTC (rev 256)
@@ -131,71 +131,52 @@
 # T = 600, 700, 800, 900, 1000 degC
 # P = 5.0GPa (50000 bar)
 # fO2 = QFM - 2
-# pH set by jadeite + kyanite + coesite (approximated here as constant)
-# output from EQ3 calculations:
+# pH set by jadeite + kyanite + coesite
+# output from EQ3NR calculations (SSH14 Supporting Information)
 # dissolved carbon: 0.03, 0.2, 1, 4, 20 molal
 # ionic strength: 0.39, 0.57, 0.88, 1.45, 2.49
-# activity coefficient (log gamma for singly charged species): -0.15, -0.18, -0.22, -0.26, -0.31
-
+# pH: 3.80, 3.99, 4.14, 4.25, 4.33
 T <- seq(600, 1000, 5)
-
 ## activate DEW model
 data(thermo)
 water("DEW")
 # add species data for DEW
 inorganics <- c("methane", "CO2", "HCO3-", "CO3-2")
 organics <- c("formic acid", "formate", "acetic acid", "acetate", "propanoic acid", "propanoate")
-
 # skip updating acetate because the new data from the DEW spreadsheet give different logK
 add.obigt("DEW", c(inorganics, organics[-4]))
 ## set basis species
 basis(c("Fe", "SiO2", "CO3-2", "H2O", "oxygen", "H+"))
-# for the time being we use a constant pH
-basis("H+", -4)
-
 ## define a QFM buffer using Berman's equations for minerals
 mod.buffer("QFM_Berman", c("quartz", "fayalite", "magnetite"), "cr_Berman", 0)
-
-## calculate logfO2 as QFM minus 2
+## calculate logfO2 in QFM buffer
 basis("O2", "QFM_Berman")
-a <- affinity(T=T, P=50000, return.buffer=TRUE)
-QFM_2 <- a$O2 - 2
-
-## now that we have QFM-2, remove QFM buffer for calculations below
-basis("O2", 0)
-
+buf <- affinity(T=T, P=50000, return.buffer=TRUE)
 ## add species
 species(c(inorganics, organics))
-
-## generate spline function for ionic strength
-IS <- splinefun(seq(600, 1000, 100), c(0.39, 0.57, 0.88, 1.45, 2.49))
-
+## generate spline functions from IS, pH, and molC values at every 100 degC
+IS <- splinefun(T[!T%%100], c(0.39, 0.57, 0.88, 1.45, 2.49))
+pH <- splinefun(T[!T%%100], c(3.80, 3.99, 4.14, 4.25, 4.33))
+molC <- splinefun(T[!T%%100], c(0.03, 0.2, 1, 4, 20))
 ## use Debye-Huckel equation with B-dot set to zero
 nonideal("Helgeson0")
-## and calculate activity coefficient for the proton
-thermo$opt$ideal.H <<- FALSE
-
-## calculate affinities on the T-logfO2 transect
-a <- affinity(T=T, O2=QFM_2, P=50000, IS=IS(T))
-
-## use the total carbon molality as an approximation of total activity
-molC <- splinefun(seq(600, 1000, 100), c(0.03, 0.2, 1, 4, 20))
-loga.C <- log10(molC(T))
-
-## calculate metastable equilibrium activities
-e <- equilibrate(a, loga.balance=loga.C)
-
+## calculate affinities on the T-logfO2-pH-IS transect
+a <- affinity(T=T, O2=buf$O2 - 2, IS=IS(T), pH=pH(T), P=50000)
+## calculate metastable equilibrium activities using the total
+## carbon molality as an approximation of total activity
+e <- equilibrate(a, loga.balance=log10(molC(T)))
 ## make the diagram; don't plot names of low-abundance species
 names <- c(inorganics, organics)
 names[c(4, 5, 7, 9)] <- ""
 col <- rep("black", length(names))
 col[c(1, 3, 6, 8, 10)] <- c("red", "darkgreen", "purple", "orange", "navyblue")
-diagram(e, alpha="balance", names=names, col=col, ylim=c(0, 0.8))
+diagram(e, alpha="balance", ylab="carbon fraction", names=names, col=col, ylim=c(0, 0.8))
 
 ## add legend and title
 ltxt1 <- "P = 50000 bar"
 ltxt2 <- substitute(logfO2=="QFM-2", list(logfO2=axis.label("O2")))
 ltxt3 <- "pH = 4"
+pH <- seq(3.8, 4.3, length.out=length(T))
 legend("left", legend=as.expression(c(ltxt1, ltxt2, ltxt3)), bty="n")
 t1 <- "Aqueous carbon speciation"
 t2 <- "after Sverjensky et al., 2014b"
@@ -217,8 +198,9 @@
 stopifnot(maxdiff(logK.calc, c(inorganic.logK, organic.logK)) < 0.021)
 
 ## check that we get similar activity coefficients
-# values for monovalent species from EQ3NR output
+# activity coefficients for monovalent species from EQ3NR output
 loggamma <- c(-0.15, -0.18, -0.22, -0.26, -0.31)
+# activity coefficients calculated in CHNOSZ
 sres <- subcrt("propanoate", T=seq(600, 1000, 100), P=50000, IS=c(0.39, 0.57, 0.88, 1.45, 2.49))
 stopifnot(maxdiff(sres$out[[1]]$loggam, loggamma) < 0.004)
 

Modified: pkg/CHNOSZ/inst/NEWS
===================================================================
--- pkg/CHNOSZ/inst/NEWS	2017-10-13 05:25:56 UTC (rev 255)
+++ pkg/CHNOSZ/inst/NEWS	2017-10-13 07:31:38 UTC (rev 256)
@@ -1,4 +1,4 @@
-CHANGES IN CHNOSZ 1.1.0-54 (2017-10-13)
+CHANGES IN CHNOSZ 1.1.0-55 (2017-10-13)
 ---------------------------------------
 
 MAJOR CHANGES:
@@ -68,6 +68,9 @@
   logfH2, or logaH2 vs pH, T, or P. It is possible to have T or P on
   either the x- or y-axis.
 
+- Variables that are in the aguments to affinity() now override any
+  buffers previously specified using basis().
+
 DATABASE UPDATES:
 
 - Add data(OBIGT) command to reset only the thermodynamic database

Modified: pkg/CHNOSZ/man/nonideal.Rd
===================================================================
--- pkg/CHNOSZ/man/nonideal.Rd	2017-10-13 05:25:56 UTC (rev 255)
+++ pkg/CHNOSZ/man/nonideal.Rd	2017-10-13 07:31:38 UTC (rev 256)
@@ -29,7 +29,8 @@
 \code{nonideal} takes a list of dataframes (in \code{speciesprops}) containing the standard molal properties of the identified \code{species}.
 The function calculates the *apparent* properties for given ionic strength (\code{IS}); they are equal to the *standard* values at IS=0.
 The function bypasses (leaves unchanged) properties of all species whose charge (determined by the number of Z in their \code{\link{makeup}}) is equal to zero.
-The proton (H+) and electron (e-) are also bypassed by default; to apply the calculations to H+ and/or e-, change \code{\link{thermo}$opt$ideal.H} or \code{ideal.e} to FALSE.
+The proton (\Hplus) and electron (\eminus) are also bypassed by default; this makes sense if you are setting the pH, i.e. activity of \Hplus, to some value.
+To apply the calculations to H+ and/or e-, change \code{\link{thermo}$opt$ideal.H} or \code{ideal.e} to FALSE.
 The lengths of \code{IS} and \code{T} supplied in the arguments should be equal to the number of rows of each dataframe in \code{speciesprops}, or length one to use single values throughout.
 
 If \code{method} is \samp{Alberty}, the values of \code{IS} are combined with Alberty's (2003) equation 3.6-1 (extended Debye-Hückel equation) and its derivatives, to calculate apparent molal properties at the specified ionic strength(s) and temperature(s).

More information about the CHNOSZ-commits mailing list