[CHNOSZ-commits] r644 - in pkg/CHNOSZ: . R demo inst man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Mar 1 10:22:59 CET 2021
Author: jedick
Date: 2021-03-01 10:22:59 +0100 (Mon, 01 Mar 2021)
New Revision: 644
Added:
pkg/CHNOSZ/demo/pourbaix.R
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/R/diagram.R
pkg/CHNOSZ/R/examples.R
pkg/CHNOSZ/demo/00Index
pkg/CHNOSZ/demo/zinc.R
pkg/CHNOSZ/inst/NEWS.Rd
pkg/CHNOSZ/man/diagram.Rd
pkg/CHNOSZ/man/examples.Rd
Log:
Add demo/pourbaix.R
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2021-03-01 00:30:30 UTC (rev 643)
+++ pkg/CHNOSZ/DESCRIPTION 2021-03-01 09:22:59 UTC (rev 644)
@@ -1,6 +1,6 @@
Date: 2021-03-01
Package: CHNOSZ
-Version: 1.4.0-13
+Version: 1.4.0-14
Title: Thermodynamic Calculations and Diagrams for Geochemistry
Authors at R: c(
person("Jeffrey", "Dick", , "j3ffdick at gmail.com", role = c("aut", "cre"),
Modified: pkg/CHNOSZ/R/diagram.R
===================================================================
--- pkg/CHNOSZ/R/diagram.R 2021-03-01 00:30:30 UTC (rev 643)
+++ pkg/CHNOSZ/R/diagram.R 2021-03-01 09:22:59 UTC (rev 644)
@@ -24,7 +24,7 @@
# sizes
cex=par("cex"), cex.names=1, cex.axis=par("cex"),
# line styles
- lty=NULL, lwd=par("lwd"), dotted=NULL,
+ lty=NULL, lty.cr=NULL, lty.aq=NULL, lwd=par("lwd"), dotted=NULL,
spline.method=NULL, contour.method="edge", levels=NULL,
# colors
col=par("col"), col.names=par("col"), fill=NULL,
@@ -532,27 +532,66 @@
}
# the categories (species/groups/etc) on the plot
zvals <- na.omit(unique(as.vector(predominant)))
- # loop over species
- for(i in 1:(length(zvals)-1)) {
- # get the "z" values
- z <- predominant
- # assign values to get one contour line between this species and all others
- i0 <- z==zvals[i]
- i1 <- z!=zvals[i]
- z[i0] <- 0
- z[i1] <- 1
- # use contourLines() instead of contour() in order to get line coordinates 20181029
- cLines <- contourLines(xs, ys, z, levels=0.5)
- if(length(cLines) > 0) {
- # loop in case contourLines returns multiple lines
- for(k in 1:length(cLines)) {
- # draw the lines
- lines(cLines[[k]][2:3], lty=lty[i], col=col[i], lwd=lwd[i])
+
+ if(is.null(lty.aq) & is.null(lty.cr)) {
+
+ # DEFAULT method: loop over species
+ for(i in 1:(length(zvals)-1)) {
+ # get the "z" values
+ z <- predominant
+ # assign values to get one contour line between this species and all others
+ i0 <- z==zvals[i]
+ i1 <- z!=zvals[i]
+ z[i0] <- 0
+ z[i1] <- 1
+ # use contourLines() instead of contour() in order to get line coordinates 20181029
+ cLines <- contourLines(xs, ys, z, levels = 0.5)
+ if(length(cLines) > 0) {
+ # loop in case contourLines returns multiple lines
+ for(k in 1:length(cLines)) {
+ # draw the lines
+ lines(cLines[[k]][2:3], lty = lty[i], col = col[i], lwd = lwd[i])
+ }
}
+ # mask species to prevent double-plotting contour lines
+ predominant[i0] <- NA
}
- # mask species to prevent double-plotting contour lines
- predominant[i0] <- NA
- }
+
+ } else {
+
+ # ALTERNATE (slower) method to handle lty.aq and lty.cr: take each possible pair of species
+ # Reinstated on 20210301
+ for(i in 1:(length(zvals)-1)) {
+ for(j in (i+1):length(zvals)) {
+ z <- predominant
+ # draw contours only for this pair
+ z[!z %in% c(zvals[i], zvals[j])] <- NA
+ # give them neighboring values (so we get one contour line)
+ z[z==zvals[i]] <- 0
+ z[z==zvals[j]] <- 1
+ # use contourLines() instead of contour() in order to get line coordinates 20181029
+ cLines <- contourLines(xs, ys, z, levels=0.5)
+ if(length(cLines) > 0) {
+ # loop in case contourLines returns multiple lines
+ for(k in 1:length(cLines)) {
+ # draw the lines
+ mylty <- lty[i]
+ if(!is.null(lty.cr)) {
+ # use lty.cr for cr-cr boundaries 20190530
+ if(all(grepl("cr", eout$species$state[c(zvals[i], zvals[j])]))) mylty <- lty.cr
+ }
+ if(!is.null(lty.aq)) {
+ # use lty.aq for aq-aq boundaries 20190531
+ if(all(grepl("aq", eout$species$state[c(zvals[i], zvals[j])]))) mylty <- lty.aq
+ }
+ lines(cLines[[k]][2:3], lty = mylty, col = col[i], lwd = lwd[i])
+ }
+ }
+ }
+ }
+
+ }
+
}
## label plot function
plot.names <- function(out, xs, ys, xlim, ylim, names, srt, min.area) {
@@ -679,19 +718,20 @@
}
pn <- list(namesx=NULL, namesy=NULL)
} else {
- # with a predominance matrix, color fields and make field boundaries
+
+ # Put predominance matrix in the right order for image()
+ zs <- t(predominant[, ncol(predominant):1])
if(!is.null(H2O.predominant)) {
- # isTRUE(limit.water): clip diagram to H2O stability region
- if(isTRUE(limit.water)) predominant <- H2O.predominant
- else {
- # is.null(limit.water): overlay diagram on H2O stability region
- zs <- t(H2O.predominant[, ncol(H2O.predominant):1])
- if(!is.null(fill)) fill.color(xs, ys, zs, fill, ngroups)
- }
+ zsH2O <- t(H2O.predominant[, ncol(H2O.predominant):1])
+ # This colors in fields and greyed-out H2O non-stability regions
+ if(!is.null(fill)) fill.color(xs, ys, zsH2O, fill, ngroups)
+ # Clip entire diagram to H2O stability region?
+ if(isTRUE(limit.water)) zs <- zsH2O
}
- # put predominance matrix in the right order for image() etc
- zs <- t(predominant[, ncol(predominant):1])
+ # This colors in fields (possibly a second time, to overlay on H2O regions)
if(!is.null(fill)) fill.color(xs, ys, zs, fill, ngroups)
+
+ # Plot field labels
pn <- plot.names(zs, xs, ys, xlim, ylim, names, srt, min.area)
# only draw the lines if there is more than one field 20180923
# (to avoid warnings from contour, which seem to be associated with weird
Modified: pkg/CHNOSZ/R/examples.R
===================================================================
--- pkg/CHNOSZ/R/examples.R 2021-03-01 00:30:30 UTC (rev 643)
+++ pkg/CHNOSZ/R/examples.R 2021-03-01 09:22:59 UTC (rev 644)
@@ -34,7 +34,7 @@
"ORP", "findit", "ionize", "buffer", "protbuff", "glycinate",
"mosaic", "copper", "arsenic", "solubility", "gold", "contour", "sphalerite", "zinc",
"Shh", "saturation", "adenine", "DEW", "lambda", "berman", "TCA", "aluminum",
- "AkDi", "comproportionation"), save.png=FALSE) {
+ "AkDi", "comproportionation", "pourbaix"), save.png=FALSE) {
# run one or more demos from CHNOSZ with ask=FALSE, and return the value of the last one
for(i in 1:length(which)) {
# say something so the user sees where we are
Modified: pkg/CHNOSZ/demo/00Index
===================================================================
--- pkg/CHNOSZ/demo/00Index 2021-03-01 00:30:30 UTC (rev 643)
+++ pkg/CHNOSZ/demo/00Index 2021-03-01 09:22:59 UTC (rev 644)
@@ -29,3 +29,4 @@
carboxylase Rank abundance distribution for RuBisCO and acetyl-CoA carboxylase
AkDi Henry's constant of dissolved gases
comproportionation Gibbs energy of sulfur comproportionation
+pourbaix Eh-pH diagram for Fe-O-H with equisolubility lines
Added: pkg/CHNOSZ/demo/pourbaix.R
===================================================================
--- pkg/CHNOSZ/demo/pourbaix.R (rev 0)
+++ pkg/CHNOSZ/demo/pourbaix.R 2021-03-01 09:22:59 UTC (rev 644)
@@ -0,0 +1,148 @@
+# CHNOSZ/demos/pourbaix.R
+# Eh-pH diagram for Fe-O-H with equisolubility lines
+# After Pourbaix (1974, p. 312)
+# 20210301 jmd first version
+
+### PARAMETERS (to be changed by the user) ###
+
+# Choose an element
+# Some mostly working ones: Fe, Cu, Au, Rh, Mn
+# Incomplete: Al (no native metal), Ni, ...
+# Not working: C, Cr, ...
+# (C gives lots of organic species, probably getting an NA affinity somewhere)
+# (Cr has no solids in OBIGT)
+element <- "Fe"
+
+# Set temperature and pressure
+T <- 25
+P <- "Psat"
+
+# Set plot limits and resolution
+pH <- c(-2, 16)
+Eh <- c(-2, 2)
+res <- 400
+
+# Assign levels for equisolubility lines
+levels <- seq(-6, 0, 2)
+
+# Switches for using colors
+color.fields <- TRUE
+color.lines <- TRUE
+color.water <- TRUE
+
+# Adjustment for most aqueous species labels
+# (e.g. to get them out of the way of the upper water line in the Fe diagram)
+dy <- -0.2
+
+### SCRIPT (can also be changed by the user if wanted!) ###
+
+# Find a species with this element
+# (to be used as a basis species)
+elem_basis <- element
+if(is.na(suppressMessages(info(elem_basis)))) {
+ elem_basis <- paste0(element, "+")
+ if(is.na(suppressMessages(info(elem_basis)))) {
+ elem_basis <- paste0(element, "+2")
+ if(is.na(suppressMessages(info(elem_basis)))) {
+ elem_basis <- paste0(element, "+3")
+ }
+ }
+}
+
+# Define system
+basis(c(elem_basis, "H2O", "H+", "e-"))
+
+# Find species
+i_cr <- retrieve(element, c("O", "H"), "cr")
+i_aq <- retrieve(element, c("O", "H"), "aq")
+
+# Add solids with unit activity
+species(i_cr, 0)
+# Add aqueous species with activity for first equisolubility line
+species(i_aq, levels[1], add = TRUE)
+
+# Calculate affinities of formation of species
+# from basis species as a function of Eh and pH
+a_all <- affinity(pH = c(pH, res), Eh = c(Eh, res), T = T, P = P)
+
+# Create labels using chemical formulas instead of name of minerals
+names_all_a <- names_all_b <- info(a_all$species$ispecies)$formula
+
+# Plot diagram (LAYER 1: colors for all fields and labels for large-enough aqueous fields)
+# (it's better to label smaller aqueous fields relative to each other than to the minerals)
+names_all_a[a_all$species$state == "cr"] <- ""
+limit.water <- fill <- NULL
+if(!color.water) limit.water <- FALSE
+if(!color.fields) fill <- NA
+d_all <- diagram(a_all, names = names_all_a, lty = 0, min.area = 0.1, limit.water = limit.water, fill = fill)
+
+# Calculate affinities for minerals
+species(i_cr)
+a_cr <- affinity(pH = c(pH, res), Eh = c(Eh, res), T = T, P = P)
+
+# Find all stable minerals across diagram
+d_cr <- diagram(a_cr, plot.it = FALSE)
+d_cr.stable <- d_cr$species$name[unique(as.vector(d_cr$predominant))]
+
+# Make a list to store the calculated solubilities for each mineral
+slist <- list()
+# Loop over stable minerals
+for(i in seq_along(d_cr.stable)) {
+ # Define basis species with mineral to dissolve
+ basis(c(d_cr.stable[i], "H2O", "H+", "e-"))
+ # Add aqueous species (no need to define activities here - they will be calculated)
+ species(i_aq)
+ # Calculate affinities of formation reactions
+ a <- affinity(pH = c(pH, res), Eh = c(Eh, res), T = T, P = P)
+ # Calculate solubility of this mineral
+ # FIXME: what to do about 'dissociation' argument?
+ s <- solubility(a, in.terms.of = element, dissociation = FALSE)
+ # Store the solubilities in the list
+ slist[[i]] <- s$loga.balance
+}
+
+# The overall solubility is the *minimum* among all the minerals
+smin <- do.call(pmin, slist)
+# Put this into the last-computed 'solubility' object
+s$loga.balance <- smin
+
+# Plot diagram (LAYER 3: equisolubility lines)
+diagram(s, type = "loga.balance", levels = levels, contour.method = "flattest", add = TRUE, lwd = 1.7)
+
+# Calculate affinities for aqueous species
+# FIXME: should be able to remove cr species from previous affinity object
+species(i_aq, 0)
+a_aq <- affinity(pH = c(pH, res), Eh = c(Eh, res), T = T, P = P)
+
+# Don't add existing labels
+d_all.ispecies_plotted <- d_all$species$ispecies[!is.na(d_all$namesx)]
+names_aq <- info(a_aq$species$ispecies)$formula
+names_aq[a_aq$species$ispecies %in% d_all.ispecies_plotted] <- ""
+
+# Plot diagram (LAYER 4: equal-activity lines and remaining labels for aqueous species)
+col <- ifelse(color.lines, 4, 1)
+# Use a white base to improve contrast
+# (so the line remains visible if it coincides with an equisolubility contour)
+diagram(a_aq, add = TRUE, col = "white", names = FALSE)
+diagram(a_aq, add = TRUE, lty = 2, col = col, names = names_aq, dy = dy)
+
+# Add solids with unit activity
+species(i_cr, 0)
+# Add aqueous species with activity for last equisolubility line
+# (i.e. unit activity)
+species(i_aq, levels[length(levels)], add = TRUE)
+
+# Calculate affinities of formation of species
+# from basis species as a function of Eh and pH
+a_all_0 <- affinity(pH = c(pH, res), Eh = c(Eh, res), T = T, P = P)
+
+# Plot diagram (LAYER 5: mineral stability boundaries and water lines)
+d_all_0 <- diagram(a_all_0, fill = NA, names = FALSE, lty = 0, lty.cr = 1, lwd = 3, add = TRUE)
+water.lines(d_all_0, lty = 5, lwd = 1.3)
+
+# Plot diagram (LAYER 6: large bold labels for all mineral fields)
+# (this is the last layer, so names are above equal-activity lines,
+# but we use positions calculated with the first equisolubility line
+# so that names are within the shrunken parts of the mineral fields)
+names_all_b[a_all$species$state == "aq"] <- ""
+diagram(a_all, fill = NA, names = names_all_b, bold = TRUE, cex.names = 1.2, lty = 0, add = TRUE)
Modified: pkg/CHNOSZ/demo/zinc.R
===================================================================
--- pkg/CHNOSZ/demo/zinc.R 2021-03-01 00:30:30 UTC (rev 643)
+++ pkg/CHNOSZ/demo/zinc.R 2021-03-01 09:22:59 UTC (rev 644)
@@ -59,7 +59,7 @@
pch = c("A", "", "B", "", "C", "", ""), inset = c(-0.1, 0), cex = 0.95)
par(xpd = FALSE)
-# Make diagram for Co minerals only 20201007
+# Make diagram for Zn minerals only 20201007
if(packageVersion("CHNOSZ") <= "1.3.6") species(delete = TRUE)
species(icr)
mcr <- mosaic(bases, pH = pH, O2 = O2, T = T, P = P, IS = sat$IS)
Modified: pkg/CHNOSZ/inst/NEWS.Rd
===================================================================
--- pkg/CHNOSZ/inst/NEWS.Rd 2021-03-01 00:30:30 UTC (rev 643)
+++ pkg/CHNOSZ/inst/NEWS.Rd 2021-03-01 09:22:59 UTC (rev 644)
@@ -10,7 +10,7 @@
\newcommand{\s}{\ifelse{latex}{\eqn{_{#1}}}{\ifelse{html}{\out{<sub>#1</sub>}}{#1}}}
\newcommand{\S}{\ifelse{latex}{\eqn{^{#1}}}{\ifelse{html}{\out{<sup>#1</sup>}}{^#1}}}
-\section{Changes in CHNOSZ version 1.4.0-11 (2021-02-25)}{
+\section{Changes in CHNOSZ version 1.4.0-14 (2021-03-01)}{
\subsection{CHANGES}{
\itemize{
@@ -35,6 +35,13 @@
Stacking 2 (minerals and aqueous species); add \emph{K}\S{eff} calculation
(\href{https://doi.org/10.1016/j.gca.2021.01.038}{Robinson et al., 2021}).
+ \item Restore \samp{lty.aq} and \samp{lty.cr} arguments to
+ \code{diagram()} to control plotting of aq-aq and cr-cr field boundaries.
+
+ \item Add \samp{demo/pourbaix.R} (Eh-pH diagrams with equisolubility
+ lines, after \href{https://www.worldcat.org/oclc/563921897}{Pourbaix,
+ 1974}).
+
}
}
Modified: pkg/CHNOSZ/man/diagram.Rd
===================================================================
--- pkg/CHNOSZ/man/diagram.Rd 2021-03-01 00:30:30 UTC (rev 643)
+++ pkg/CHNOSZ/man/diagram.Rd 2021-03-01 09:22:59 UTC (rev 644)
@@ -21,7 +21,7 @@
# character sizes
cex = par("cex"), cex.names = 1, cex.axis = par("cex"),
# line styles
- lty = NULL, lwd = par("lwd"), dotted = NULL,
+ lty = NULL, lty.cr = NULL, lty.aq = NULL, lwd = par("lwd"), dotted = NULL,
spline.method = NULL, contour.method = "edge", levels = NULL,
# colors
col = par("col"), col.names = par("col"), fill = NULL,
@@ -58,6 +58,8 @@
\item{cex.names}{numeric, character expansion factor to be used for names of species on plots}
\item{cex.axis}{numeric, character expansion factor for names of axes}
\item{lty}{numeric, line types to be used in plots}
+ \item{lty.cr}{numeric, line types for cr-cr boundaries (between two minerals)}
+ \item{lty.aq}{numeric, line types for aq-aq boundaries (between two aqueous species)}
\item{lwd}{numeric, line width}
\item{dotted}{numeric, how often to skip plotting points on predominance field boundaries (to gain the effect of dotted or dashed boundary lines)}
\item{spline.method}{character, method used in \code{\link{splinefun}}}
@@ -160,6 +162,8 @@
The default line-drawing algorithm uses \code{\link{contourLines}} to obtain smooth-looking diagonal and curved lines, at the expense of not coinciding exactly with the rectangular grid that is used for drawing colors.
\code{lty}, \code{col}, and \code{lwd} can be specified, but limiting the lines via \code{xrange} is not currently supported.
+Set \code{lty.cr} or \code{lty.aq} to 0 to suppress boundary lines between minerals or aqueous species.
+
To go back to the old behavior for drawing lines, set \code{dotted} to \samp{0}.
The old behavior does not respect \code{lty}; instead, the style of the boundary lines on 2-D diagrams can be altered by supplying one or more non-zero integers in \code{dotted}, which indicates the fraction of line segments to omit; a value of \samp{1} or NULL for \code{dotted} has the effect of not drawing the boundary lines.
Modified: pkg/CHNOSZ/man/examples.Rd
===================================================================
--- pkg/CHNOSZ/man/examples.Rd 2021-03-01 00:30:30 UTC (rev 643)
+++ pkg/CHNOSZ/man/examples.Rd 2021-03-01 09:22:59 UTC (rev 644)
@@ -18,7 +18,7 @@
"glycinate", "mosaic", "copper", "arsenic", "solubility", "gold",
"contour", "sphalerite", "zinc", "Shh", "saturation",
"adenine", "DEW", "lambda", "berman", "TCA", "aluminum", "AkDi",
- "comproportionation"),
+ "comproportionation", "pourbaix"),
save.png=FALSE)
}
@@ -62,6 +62,7 @@
\code{carboxylase} \tab Rank abundance distribution for RuBisCO and acetyl-CoA carboxylase \cr
\code{AkDi} \tab Henry's constant of dissolved gases (Akinfiev and Diamond, 2003) \cr
\code{comproportionation} \tab Gibbs energy of sulfur comproportionation (Amend et al., 2020) \cr
+ \code{pourbaix} \tab Eh-pH diagram for Fe-O-H with equisolubility lines (Pourbaix, 1974) \cr
}
For either function, if \code{save.png} is TRUE, the plots are saved in \code{\link{png}} files whose names begin with the names of the help topics or demos.
@@ -121,6 +122,8 @@
Manning, C. E., Shock, E. L. and Sverjensky, D. A. (2013) The chemistry of carbon in aqueous fluids at crustal and upper-mantle conditions: Experimental and theoretical constraints. \emph{Rev. Mineral. Geochem.} \bold{75}, 109--148. \url{https://doi.org/10.2138/rmg.2013.75.5}
+Pourbaix, M. (1974) \emph{Atlas of Electrochemical Equilibria in Aqueous Solutions}, NACE, Houston, TX and CEBELCOR, Brussels. \url{https://www.worldcat.org/oclc/563921897}
+
Schulte, M. D. and Shock, E. L. (1995) Thermodynamics of Strecker synthesis in hydrothermal systems. \emph{Orig. Life Evol. Biosph.} \bold{25}, 161--173. \url{https://doi.org/10.1007/BF01581580}
Shock, E. L. and Koretsky, C. M. (1995) Metal-organic complexes in geochemical processes: Estimation of standard partial molal thermodynamic properties of aqueous complexes between metal cations and monovalent organic acid ligands at high pressures and temperatures. \emph{Geochim. Cosmochim. Acta} \bold{59}, 1497--1532. \url{https://doi.org/10.1016/0016-7037(95)00058-8}
More information about the CHNOSZ-commits
mailing list