[CHNOSZ-commits] r572 - in pkg/CHNOSZ: . R man vignettes
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Jul 22 19:20:44 CEST 2020
Author: jedick
Date: 2020-07-22 19:20:44 +0200 (Wed, 22 Jul 2020)
New Revision: 572
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/R/diagram.R
pkg/CHNOSZ/R/mix.R
pkg/CHNOSZ/man/mix.Rd
pkg/CHNOSZ/vignettes/mklinks.sh
pkg/CHNOSZ/vignettes/multi-metal.Rmd
pkg/CHNOSZ/vignettes/vig.bib
Log:
Add mixing example (Fe-V-O-H) to multi-metal.Rmd
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2020-07-22 03:11:31 UTC (rev 571)
+++ pkg/CHNOSZ/DESCRIPTION 2020-07-22 17:20:44 UTC (rev 572)
@@ -1,6 +1,6 @@
Date: 2020-07-22
Package: CHNOSZ
-Version: 1.3.6-45
+Version: 1.3.6-46
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 2020-07-22 03:11:31 UTC (rev 571)
+++ pkg/CHNOSZ/R/diagram.R 2020-07-22 17:20:44 UTC (rev 572)
@@ -259,7 +259,7 @@
## apply formatting to chemical formulas 20170204
if(all(grepl("_", names))) is.pname <- TRUE
if(format.names & !is.pname) {
- # check if names are a deparsed expression (used in combine()) 20200718
+ # check if names are a deparsed expression (used in mix()) 20200718
parsed <- FALSE
if(any(grepl("paste\\(", names))) {
exprnames <- parse(text = names)
@@ -268,7 +268,12 @@
} else {
exprnames <- as.expression(names)
# get formatted chemical formulas
- for(i in seq_along(exprnames)) exprnames[[i]] <- expr.species(exprnames[[i]])
+ for(i in seq_along(exprnames)) {
+ # don't try to format the names if they have "+" followed by a character
+ # (created by mix()ing diagrams with format.names = FALSE);
+ # expr.species() can't handle it 20200722
+ if(!grepl("\\+[a-zA-Z]", names[i])) exprnames[[i]] <- expr.species(exprnames[[i]])
+ }
}
# apply bold or italic
bold <- rep(bold, length.out = length(exprnames))
Modified: pkg/CHNOSZ/R/mix.R
===================================================================
--- pkg/CHNOSZ/R/mix.R 2020-07-22 03:11:31 UTC (rev 571)
+++ pkg/CHNOSZ/R/mix.R 2020-07-22 17:20:44 UTC (rev 572)
@@ -2,6 +2,9 @@
# Combine diagrams for two metals
# 20200713 first version jmd
+## if this file is interactively sourced, the following are also needed to provide unexported functions:
+#source("equilibrate.R")
+
# Function to combine two diagrams (simple overlay, no interaction) 20200717
# -- makes new "species" from all combinations of those in d1 and d2
flatten <- function(d1, d2) {
@@ -12,21 +15,23 @@
# mix the systems to include bimetal species 20200721
# 'parts' gives number of moles of each metal in 1 mole of the mixture
mix <- function(d1, d2, d3 = NULL, parts = c(1, 1), .balance = NULL) {
- check_d1_d2(d1, d2, .balance)
+ if(is.null(.balance)) check_d1_d2(d1, d2)
+ else check_d1_d3(d1, d3)
# Handle mixing here
- if(!is.null(d3)) {
+ if(!is.null(d3) & is.null(.balance)) {
# Mix d1 (e.g. Fe only) and d2 (e.g. V only)
m12 <- mix(d1, d2, parts = parts)
# Mix d1 (Fe) and d3 (FeV bimetallic species)
# - '.balance' is defined to add d3 to get required amount of V
- m13 <- mix(d1, d3, parts = parts, .balance = d2$balance)
+ # - second 'd3' in argument list is used only for check_d1_d3()
+ m13 <- mix(d1, d3, d3, parts = parts, .balance = d2$balance)
# Mix d2 (V) and d3 (FeV)
# - '.balance' is defined to add d3 to get required amount of Fe
# - d2 (V) is first, so we need to reverse the 'parts' values
- m23 <- mix(d2, d3, parts = rev(parts), .balance = d1$balance)
+ m23 <- mix(d2, d3, d3, parts = rev(parts), .balance = d1$balance)
# Mix d3 with itself (combinations of bimetallic species)
- m33 <- mix(d3, d3, parts = parts, .balance = c(d1$balance, d2$balance))
+ m33 <- mix(d3, d3, d3, parts = parts, .balance = c(d1$balance, d2$balance))
# Merge all the species and affinity values
species <- rbind(m12$species, m13$species, m23$species, m33$species)
values <- c(m12$values, m13$values, m23$values, m33$values)
@@ -68,11 +73,11 @@
if(!is.null(.balance)) {
if(length(.balance)==2) {
# For combinations of bimetallic species (m33)
- b1 <- balance(d1, .balance[1])$n.balance[combs[, 1]]
- b2 <- balance(d2, .balance[2])$n.balance[combs[, 2]]
- ibal <- match(balance, colnames(s1))
+ b1 <- suppressMessages(balance(d1, .balance[1])$n.balance[combs[, 1]])
+ b2 <- suppressMessages(balance(d2, .balance[2])$n.balance[combs[, 2]])
+ ibal <- match(.balance, colnames(s1))
} else {
- b2 <- balance(d2, .balance)$n.balance[combs[, 2]]
+ b2 <- suppressMessages(balance(d2, .balance)$n.balance[combs[, 2]])
ibal <- match(c(d1$balance, .balance), colnames(s1))
}
}
@@ -85,14 +90,22 @@
frac1 <- c(frac1, x[1])
frac2 <- c(frac2, x[2])
}
- # Note that some of frac1 might be < 0 ... we remove these combinations below
+ # Note that some of frac1 or frac2 might be < 0 ... we remove these combinations below
# Make a new species data frame starting with sum of scaled formation reactions
nbasis <- nrow(d1$basis)
species <- s1[, 1:nbasis] * frac1 + s2[, 1:nbasis] * frac2
+ # Concatenate ispecies, logact, state
ispecies <- paste(s1$ispecies, s2$ispecies, sep = ",")
logact <- paste(s1$logact, s2$logact, sep = ",")
state <- paste(s1$state, s2$state, sep = ",")
+ # Omit species with zero mole fraction from concatenated values 20200722
+ ispecies[frac1 == 0] <- s2$ispecies[frac1 == 0]
+ logact[frac1 == 0] <- s2$logact[frac1 == 0]
+ state[frac1 == 0] <- s2$state[frac1 == 0]
+ ispecies[frac2 == 0] <- s1$ispecies[frac2 == 0]
+ logact[frac2 == 0] <- s1$logact[frac2 == 0]
+ state[frac2 == 0] <- s1$state[frac2 == 0]
# Use names from diagram()
if(is.expression(d1$names) | is.expression(d2$names)) {
# Convert non-expressions to lists so we can use [[ indexing below
@@ -99,10 +112,15 @@
if(!is.expression(d1$names)) d1$names <- as.list(d1$names)
if(!is.expression(d2$names)) d2$names <- as.list(d2$names)
name <- lapply(1:nrow(combs), function(i) {
- # Don't include names of species that are not present (zero or negative mole fractions)
- if(frac1[i] <= 0) bquote(.(d2$names[[combs[i, 2]]]))
- else if(frac2[i] <= 0) bquote(.(d1$names[[combs[i, 1]]]))
- else bquote(.(d1$names[[combs[i, 1]]])+.(d2$names[[combs[i, 2]]]))
+ # Don't include names of species that are not present (zero mole fractions)
+ if(frac1[i] == 0 & frac2[i] == 0) bquote("")
+ else if(frac1[i] == 0) bquote(.(d2$names[[combs[i, 2]]]))
+ else if(frac2[i] == 0) bquote(.(d1$names[[combs[i, 1]]]))
+ else {
+ # Put the names together, with the species from d2 first if it is has a higher mole fraction 20200722
+ if(frac2[i] > frac1[i]) bquote(.(d2$names[[combs[i, 2]]])+.(d1$names[[combs[i, 1]]]))
+ else bquote(.(d1$names[[combs[i, 1]]])+.(d2$names[[combs[i, 2]]]))
+ }
})
name <- unlist(lapply(name, deparse, width.cutoff = 500, control = NULL))
if(length(name) != nrow(combs)) stop("deparse()-ing expressions gives unequal length; try diagram(., format.names = FALSE)")
@@ -127,7 +145,7 @@
ipredominant <- logical(length(values))
if(length(.balance)==2) {
# For combinations of bimetallic species (m33), don't do predominance masking
- # (there are not mono-metallic species to look at)
+ # (there are no mono-metallic species to look at)
ipredominant[] <- TRUE
} else {
# Loop over combinations to find predominant species in the single-metal diagram(s)
@@ -136,7 +154,8 @@
i1 <- combs[i, 1]
ip1 <- d1$predominant == i1
# If the mole fraction is zero, it is predominant by definition
- # (this allows a single bimetallic species to be formed)
+ # (this allows a single bimetallic species (that is paired with
+ # the zero-mole-fraction mono-metallic species) to be formed)
if(frac1[i] == 0) ip1[] <- TRUE
if(is.null(.balance)) {
# Get predominant species in second diagram
@@ -143,7 +162,7 @@
i2 <- combs[i, 2]
ip2 <- d2$predominant == i2
# If the mole fraction is zero, it is predominant by definition
- # (this allows us to recover the single-element diagram with frac = c(0, 1) or c(1, 0))
+ # (this allows us to recover the first single-element diagram with frac = c(1, 0))
if(frac2[i] == 0) ip2[] <- TRUE
# Assign -Inf affinity where any species isn't predominant in the corresponding single-metal diagram
ip12 <- ip1 & ip2
@@ -157,9 +176,9 @@
}
}
# Remove combinations that:
- # 1) involve a species with no predominance field in the single-metal diagram or
- # 2) have a negative mole fraction of a single-metal species
- inotnegative <- frac1 >= 0
+ # 1) involve a mono-metallic species with no predominance field in the corresponding single-metal diagram or
+ # 2) have a negative mole fraction of any species
+ inotnegative <- frac1 >= 0 & frac2 >= 0
values <- values[ipredominant & inotnegative]
species <- species[ipredominant & inotnegative, ]
@@ -234,21 +253,35 @@
}
-### unexported function ###
+### unexported functions ###
# Check that d1 and d2 can be combined
# Extracted from duplex() (now rebalance()) 20200717
-check_d1_d2 <- function(d1, d2, .balance = NULL) {
- if(is.null(.balance)) d2txt <- "d2" else d2txt <- "d3"
+check_d1_d2 <- function(d1, d2) {
# Check that the basis species are the same
- if(!identical(d1$basis, d2$basis)) stop(paste0("basis species in objects 'd1' and '", d2txt, "' are not identical"))
+ if(!identical(d1$basis, d2$basis)) stop(paste0("basis species in objects 'd1' and 'd2' are not identical"))
# Check that the variables and their values are the same
- if(!identical(d1$vars, d2$vars)) stop(paste0("variable names in objects 'd1' and '", d2txt, "' are not identical"))
- if(!identical(d1$vals, d2$vals)) stop(paste0("variable values in objects 'd1' and '", d2txt, "' are not identical"))
+ if(!identical(d1$vars, d2$vars)) stop(paste0("plot variables in objects 'd1' and 'd2' are not identical"))
+ if(!identical(d1$vals, d2$vals)) stop(paste0("plot ranges in objects 'd1' and 'd2' are not identical"))
# Check that T and P are the same
- if(!identical(d1$T, d2$T)) stop(paste0("temperatures in objects 'd1' and '", d2txt, "' are not identical"))
- if(!identical(d1$P, d2$P)) stop(paste0("pressures in objects 'd1' and '", d2txt, "' are not identical"))
+ if(!identical(d1$T, d2$T)) stop(paste0("temperatures in objects 'd1' and 'd2' are not identical"))
+ if(!identical(d1$P, d2$P)) stop(paste0("pressures in objects 'd1' and 'd2' are not identical"))
# Check that we have plotvals and predominant (from diagram())
if(is.null(d1$plotvals) | is.null(d1$predominant)) stop("object 'd1' is missing 'plotvals' or 'predominant' components (not made by diagram()?)")
- if(is.null(d2$plotvals) | is.null(d2$predominant)) stop(paste0("object '", d2txt, "' is missing 'plotvals' or 'predominant' components (not made by diagram()?)"))
+ if(is.null(d2$plotvals) | is.null(d2$predominant)) stop(paste0("object 'd2' is missing 'plotvals' or 'predominant' components (not made by diagram()?)"))
}
+
+# Check that d1 and d3 can be combined
+check_d1_d3 <- function(d1, d3) {
+ # Check that the basis species are the same
+ if(!identical(d1$basis, d3$basis)) stop(paste0("basis species in objects 'd1' and 'd3' are not identical"))
+ # Check that the variables and their values are the same
+ if(!identical(d1$vars, d3$vars)) stop(paste0("plot variables in objects 'd1' and 'd3' are not identical"))
+ if(!identical(d1$vals, d3$vals)) stop(paste0("plot ranges in objects 'd1' and 'd3' are not identical"))
+ # Check that T and P are the same
+ if(!identical(d1$T, d3$T)) stop(paste0("temperatures in objects 'd1' and 'd3' are not identical"))
+ if(!identical(d1$P, d3$P)) stop(paste0("pressures in objects 'd1' and 'd3' are not identical"))
+ # Check that we have plotvals and predominant (from diagram())
+ if(is.null(d1$plotvals) | is.null(d1$predominant)) stop("object 'd1' is missing 'plotvals' or 'predominant' components (not made by diagram()?)")
+ if(is.null(d3$plotvals) | is.null(d3$predominant)) stop(paste0("object 'd3' is missing 'plotvals' or 'predominant' components (not made by diagram()?)"))
+}
Modified: pkg/CHNOSZ/man/mix.Rd
===================================================================
--- pkg/CHNOSZ/man/mix.Rd 2020-07-22 03:11:31 UTC (rev 571)
+++ pkg/CHNOSZ/man/mix.Rd 2020-07-22 17:20:44 UTC (rev 572)
@@ -36,14 +36,14 @@
\code{mix} is an expanded form of \code{flatten} that allows combinations not only between two single-metal diagrams (\code{d1} and \code{d2}) but also between each of those diagrams and third diagram for bimetallic species (\code{d3}).
All combinations of species in all crosses between the diagrams (\code{d1-d2}, \code{d1-d3}, \code{d2-d3}, \code{d3-d3}) are identified.
The mole fractions of species in each combination are computed to satisfy the ratio of metals defined in \code{parts}.
-For example, if \code{d1} and \code{d2} are balanced on Fe\s{+2} and VO\S{4}\s{-3}, the species are combined by default to give equal parts of Fe and V.
-Note that the \code{d3-d3} combinations include pairs of distinct bimetallic species as well as single bimetallic species that satisfy the ratios (e.g. FeV for \code{c(1, 1)} or Fe3V for \code{c(3, 1)}).
+For example, if \code{d1} and \code{d2} are balanced on Fe\S{+2} and VO\s{4}\S{-3}, the species are combined by default to give equal parts of Fe and V.
+Note that pairs of distinct bimetallic species in \code{d3} are included as well as single bimetallic species that satisfy the composition in \code{parts} (e.g. FeV for \code{c(1, 1)} or Fe\s{3}V for \code{c(3, 1)}).
(\code{parts} can be expressed as integers or fractions, so \code{c(0.75, 0.25)} is equivalent to \code{c(3, 1)}.)
-From the possible combinations of species, combinations are removed that involve any species that has no predominance field in the respective single-metal diagram or that appears in the combination with a negative mole fraction.
+From the possible combinations of species, combinations are removed that have a negative mole fraction of any species or that involve any mono-metallic species that has no predominance field in the corresponding single-metal diagram.
The output consists of each unique combination of species, including the combined formation reactions and affinities (in the \code{species} and \code{values} elements of the output list),
The affinities are assigned values of -Inf wherever one of the species is not predominant in the respective single-metal diagram.
-Therefore, either the single-metal diagrams can be recovered by setting \code{parts} to \code{c(1, 0)} (for \code{d1}) or \code{c(0, 1)} (for \code{d2}).
+Therefore, either the single-metal diagrams (\code{d1} or \code{d2}) can be recovered by setting \code{parts} to \code{c(1, 0)} or \code{c(0, 1)}, respectively.
NOTE: Unlike the \code{diagram} calls used to make \code{d1} and \code{d2}, which by themselves should produce reasonable diagrams for a single-metal system, the \code{d3} diagram by itself probably has no useful interpretation.
It is only used in \code{mix} as a way to transmit the results of \code{\link{affinity}} for the bimetmal system and the formatted names that are made by \code{diagram}.
Modified: pkg/CHNOSZ/vignettes/mklinks.sh
===================================================================
--- pkg/CHNOSZ/vignettes/mklinks.sh 2020-07-22 03:11:31 UTC (rev 571)
+++ pkg/CHNOSZ/vignettes/mklinks.sh 2020-07-22 17:20:44 UTC (rev 572)
@@ -116,3 +116,5 @@
sed -i 's/equilibrate()/<a href="..\/html\/equilibrate.html">equilibrate()<\/a>/g' multi-metal.html
sed -i 's/rebalance()/<a href="..\/html\/mix.html">rebalance()<\/a>/g' multi-metal.html
sed -i 's/ratlab()/<a href="..\/html\/util.expression.html">ratlab()<\/a>/g' multi-metal.html
+sed -i 's/mix()/<a href="..\/html\/mix.html">mix()<\/a>/g' multi-metal.html
+sed -i 's/mod.OBIGT()/<a href="..\/html\/add.OBIGT.html">mod.OBIGT()<\/a>/g' multi-metal.html
Modified: pkg/CHNOSZ/vignettes/multi-metal.Rmd
===================================================================
--- pkg/CHNOSZ/vignettes/multi-metal.Rmd 2020-07-22 03:11:31 UTC (rev 571)
+++ pkg/CHNOSZ/vignettes/multi-metal.Rmd 2020-07-22 17:20:44 UTC (rev 572)
@@ -74,12 +74,11 @@
Basic diagrams in CHNOSZ are made for reactions that are *balanced on an element* (see [Equilibrium in CHNOSZ](equilibrium.html)) and therefore represent minerals or aqueous species that all have one element, often a metal, in common.
The package documentation has many examples of diagrams for a single metal appearing in different minerals or complexed with different ligands, but a common request is to make diagrams for multiple metals.
This vignette describes some methods for constructing diagrams for multi-metal minerals and other multi-element systems.
-The methods are **flattening**, **mosaic series**, and **secondary balancing**.
+The methods are **flattening**, **mixing**, **mosaic stacking**, and **secondary balancing**.
## Flattening
Flattening or simple overlay refers to independent calculations for two different systems that are displayed on the same diagram.
-Although it is easy to make such a diagram, there is no interaction between the systems.
This example starts with a log*f*~O<sub>2</sub>~--pH base diagram for the C-O-H system then overlays a diagram for S-O-H.
The second call to `affinity()` uses the argument recall feature, where the arguments are taken from the previous command.
@@ -110,13 +109,183 @@
Tip: the names of the fields in the second diagram come from `aCS$species$name`, which are expressions made by combining `aC$names` and `aS$names`.
If you prefer plain text names without formatting, add `format.names = FALSE` to all of the `diagram()` calls.
-## Mosaic Series 1
+## Mixing
+In simple terms, flattening diagrams portrays the result of mixing two single-metal systems.
+Although it is easy to make such a diagram, there is no interaction between the systems.
+If there is a possibility of forming bimetallic species, then additional considerations are needed to account for the stoichiometry of the mixture.
+The stoichiometry can be given as a fixed composition of both metals; then, all combinations of (mono- and/or bimetallic) species that satisfy this compositional constraint are used as the candidate "species" in the system.
+This is the same type of calculation as that described for [binary Pourbaix diagrams in the Materials Project](https://matgenb.materialsvirtuallab.org/2017/12/15/Plotting-a-Pourbaix-Diagram.html#Plotting-k-nary-systems).
+
+This example makes a Pourbaix diagram for the Fe-V-O-H system that is similar to Figure 1 of @SZS_17.
+First we need to assemble the energies of the solids and aqueous species.
+For solids, values of formation energy (in eV / atom) computed using density functional theory (DFT) are taken from the Materials Project website and are converted to units of J / mol.
+For aqueous species, values of standard Gibbs energy of formation at 25 °C (in J / mol) are taken mostly from @WEP_82 augmented with data for FeO~2~^-^ from @SSWS97 and FeO~4~^-2^ from @Mis73.
+Adapting the method described by @PWLC12, a reference-state correction for each metal is calculated from the difference between the DFT-based formation energy and the standard Gibbs energy of a representative material; here we use values for Fe~3~O~4~ (magnetite) and V~3~O~5~ from @WEP_82.
+This correction is then applied to all of the aqueous species that have that metal.
+Finally, `mod.OBIGT()` is used to add the obtained energies to the OBIGT database.
+
+<button id="B-materials" onclick="ToggleDiv('materials')">Show code</button>
+<div id="D-materials" style="display: none">
+```{r materials, message = FALSE, results = "hide"}
+## Formation energies (eV / atom) for solids from MP website
+# mp-13, mp-1279742, mp-19306, mp-19770
+Fe.cr <- c(Fe = 0, FeO = -1.728, Fe3O4 = -1.859, Fe2O3 = -1.907)
+# mp-146, mp-18937, mp-1275946, mp-19094, mp-756395, mp-25279
+V.cr <- c(V = 0, V2O3 = -2.528, V3O5 = -2.526, VO2 = -2.485, V3O7 = -2.328, V2O5 = -2.295)
+
+# Convert formation energies from eV / atom to eV / molecule
+natom.Fe <- sapply(makeup(names(Fe.cr)), sum)
+Fe.cr <- Fe.cr * natom.Fe
+natom.V <- sapply(makeup(names(V.cr)), sum)
+V.cr <- V.cr * natom.V
+
+# Convert formation energies from eV / molecule to J / mol
+eVtoJ <- function(eV) eV * 1.602176634e-19 * 6.02214076e23
+Fe.cr <- eVtoJ(Fe.cr)
+V.cr <- eVtoJ(V.cr)
+
+# Gibbs energies of formation (J / mol) for aqueous species from Wagman et al.
+Fe.aq <- 1000 * c("Fe+2" = -78.90, "Fe+3" = -4.7, "FeO2-2" = -295.3,
+ "FeOH+" = -277.4, "FeOH+2" = -229.41, "HFeO2-" = -377.7,
+ "Fe(OH)2+" = -438.0, "Fe(OH)3" = -659.3,
+ "FeO2-" = -368.2, # SSWS97
+ "FeO4-2" = -322.2 # Mis73
+)
+Fe.aq <- 1000 * c("Fe+2" = -78.90, "FeOH+" = -277.4)
+
+V.aq <- 1000 * c("VO+2" = -446.4, "VO2+" = -587.0, "VO3-" = -783.6, "VO4-3" = -899.0,
+ "V2O7-4" = -1719, "HVO4" = -745.1, "HVO4-2" = -974.8,
+ "VOH2O2+3" = -523.4, "VO2H2O2+" = -746.3, "V2HO7-3" = 1792.2, "V2H3O7-" = -1863.8,
+ "HV10O28-5" = -7702, "H2V10O28-4" = -7723
+)
+V.aq <- 1000 * c("VO+2" = -446.4, "HVO4-2" = -974.8)
+
+# Gibbs energies of formation (J / mol) for solids from Wagman et al., 1982
+Fe3O4 <- 1000 * -1015.4 # magnetite
+V3O5 <- 1000 * -1803
+
+# Calculate correction for difference between reference and DFT energies (Persson et al., 2012)
+Fe.corr <- (Fe.cr["Fe3O4"] - Fe3O4) / 3
+V.corr <- (V.cr["V3O5"] - V3O5) / 2
+
+# Apply correction to Gibbs energies of aqueous species (Persson et al., 2012)
+nFe <- sapply(makeup(names(Fe.aq)), "[", "Fe")
+Fe.aq <- Fe.aq + nFe * Fe.corr
+nV <- sapply(makeup(names(V.aq)), "[", "V")
+V.aq <- V.aq + nV * V.corr
+
+# Add energies to OBIGT
+# Be sure to set E_units to Joules; the default in OBIGT is calories!
+# This function modifies OBIGT and returns the species indices of the affected species
+modfun <- function(x, state) sapply(seq_along(x), function(i) {
+ mod.OBIGT(names(x)[i], formula = names(x)[i], state = state, E_units = "J", G = x[i])
+})
+iFe.cr <- modfun(Fe.cr, "cr")
+iFe.aq <- modfun(Fe.aq, "aq")
+iV.cr <- modfun(V.cr, "cr")
+iV.aq <- modfun(V.aq, "aq")
+
+# Formation energies (eV / atom) for bimetallic solids from MP website
+# mp-1335, mp-1079399, mp-866134
+FeV.cr <- c(VFe = -0.129, V3Fe = -0.171, VFe3 = -0.129)
+# Convert energies and add to OBIGT
+natom.FeV <- sapply(makeup(names(FeV.cr)), sum)
+FeV.cr <- FeV.cr * natom.FeV
+FeV.cr <- eVtoJ(FeV.cr)
+iFeV.cr <- modfun(FeV.cr, "cr")
+```
+</div>
+
+This code produce species indices in the OBIGT database for Fe- and V-bearing aqueous species (`iFe.aq`, `iV.aq`), solids (`iFe.cr`, `iV.cr`), and bimetallic solids (`iFeV.cr`), which are used in the following diagrams.
+
+Now we set up the plot area and write a function to assign the colors for regions with two solids, one solid, and no solids; these are the web colors "burlywood", "antiquewhite" and "aliceblue" with some transparency added to show the underlying water stability region.
+The activities of aqueous species are set to the default values for the diagrams produced by the link "Generate Phase Diagram" -- "Aqueous Stability (Pourbaix)" on a material information page in the Materials Project (MP).
+The pH and Eh ranges are made relatively small in order to show just a part of the diagram.
+
+```{r mix, eval = FALSE, echo = 1:12}
+mat <- matrix(c(1, 2, 2, 3, 3, 0, 4, 4, 5, 5, 6, 6), byrow = TRUE, nrow = 2)
+layout(mat)
+plot.new()
+text(0.5, 0.5, lTP(25, 1), cex = 1.4)
+fill <- function(a) {
+ ifelse(grepl("cr,cr", a$species$state), "#DEB88788",
+ ifelse(grepl("cr", a$species$state), "#FAEBD788", "#F0F8FF88")
+)}
+loga.Fe <- -5
+loga.V <- -5
+pH <- c(4, 10)
+Eh <- c(-1.5, 0.5)
+
+# Fe-O-H diagram
+basis(c("VO+2", "Fe+2", "H2O", "e-", "H+"))
+species(c(iFe.aq, iFe.cr))
+species(1:length(iFe.aq), loga.Fe)
+aFe <- affinity(pH = pH, Eh = Eh)
+dFe <- diagram(aFe, fill = fill(aFe))
+water.lines(dFe, col = 4)
+title("Fe-O-H")
+# V-O-H diagram
+species(c(iV.aq, iV.cr))
+species(1:length(iV.aq), loga.V)
+aV <- affinity(aFe) # argument recall
+dV <- diagram(aV, fill = fill(aV))
+water.lines(dFe, col = 4)
+title("V-O-H")
+
+# Calculate affinities for bimetallic species
+species(iFeV.cr)
+aFeV <- affinity(aFe) # argument recall
+dFeV <- diagram(aFeV, plot.it = FALSE, bold = TRUE)
+# 1:1 mixture (Fe:V)
+aFeV11 <- mix(dFe, dV, dFeV, c(1, 1))
+diagram(aFeV11, fill = fill(aFeV11), min.area = 0.01)
+water.lines(dFe, col = 4)
+title("Fe:V = 1:1")
+# 1:3 mixture
+aFeV13 <- mix(dFe, dV, dFeV, c(1, 3))
+diagram(aFeV13, fill = fill(aFeV13), min.area = 0.01)
+water.lines(dFe, col = 4)
+title("Fe:V = 1:3")
+# 1:5 mixture
+aFeV15 <- mix(dFe, dV, dFeV, c(1, 5))
+diagram(aFeV15, fill = fill(aFeV15), min.area = 0.01)
+water.lines(dFe, col = 4)
+title("Fe:V = 1:5")
+```
+
+The next group of commands makes Eh-pH diagrams for the single-metal systems Fe-O-H and V-O-H.
+The output of `diagram()` is saved in `dFe` and `dV` for later use.
+
+```{r mix, eval = FALSE, echo = 14:28}
+```
+
+Then we calculate the affinities for the bimetallic species and save the output of `diagram()` in `dFeV` without making a plot, but formatting the names in bold.
+Finally, `mix()` is used to generate diagrams for three different compositions.
+The `min.area` argument in `diagram()` is used to remove labels for very small fields.
+
+```{r mix, echo = 27:47, message = FALSE, results = "hide", fig.width = 10, fig.height = 6.5, out.width = "100%", out.extra='class="full-width"', pngquant = pngquant}
+```
+
+In these diagrams, changing the Fe:V ratio affects the fully reduced metallic species.
+In the 1:1 mixture, the V~3~Fe + VFe~3~ assemblage is stable instead of VFe.
+This result is unlike Figure 1 of @SZS_17 but is consistent with the [MP page for VFe](https://doi.org/10.17188/1189535) where it shown to decompose to this assemblage.
+On the other hand, [V~3~Fe is stable](https://materialsproject.org/materials/mp-1079399/) in the 1:3 mixture.
+For an even higher proportion of V, the V + V~3~Fe assemblage is stable, which can also be seen in the Pourbaix diagram linked from the [MP page for V~5~FeO~12~](https://doi.org/10.17188/1305091).
+
+Because this example modified the thermodynamic data for some minerals that are used below, we should restore the default OBIGT database before proceeding to the next section.
+
+```{r reset, message = FALSE}
+reset()
+```
+
+## Mosaic Stacking 1
+
A mosaic diagram shows the effects of changing basis species on the stabilities of minerals.
The Fe-S-O-H system is a common example: the speciation of aqueous sulfur species affects the stabilities of iron oxides and sulfides.
Examples of mosaic diagrams with Fe or other single metals are given elsewhere.
-A mosaic series is when predominance fields for minerals calculated in one mosaic diagram are used as input to a second mosaic diagram, where the minerals are now themselves basis species.
+A mosaic stack is when predominance fields for minerals calculated in one mosaic diagram are used as input to a second mosaic diagram, where the minerals are now themselves basis species.
The example here shows the construction of a Cu-Fe-S-O-H diagram.
First we define the conditions and basis species.
@@ -123,7 +292,7 @@
It is important to put Cu^+^ first so that it will be used as the balance for the reactions with Cu-bearing minerals (which also have Fe).
Pyrite is chosen as the starting Fe-bearing basis species, which will be changed as indicated in `bases2`.
-```{r series1_1, results = "hide", message = FALSE}
+```{r stack1_1, results = "hide", message = FALSE}
logaH2S <- -2
T <- 200
pH <- c(0, 12, 500)
@@ -141,7 +310,7 @@
2. Predominance fields for the aqueous S species (italic blue text and dashed lines)
3. Stability areas for the Fe-bearing minerals (black text and solid lines)
-```{r series1_2, eval = FALSE, echo = 1:6}
+```{r stack1_2, eval = FALSE, echo = 1:6}
species(bases2)
mFe <- mosaic(bases1, pH = pH, O2 = O2, T = T)
diagram(mFe$A.bases, lty = 2, col = 4, col.names = 4,
@@ -158,7 +327,7 @@
```
Next we load the Cu-bearing minerals and calculate their affinities while changing *both* the aqueous sulfur species and the Fe-bearing minerals whose stability fields were just calculated.
-The latter step is the key to the mosaic series and is activated by supplying the calculated stabilities of the Fe-bearing minerals in the `predominant` argument.
+The latter step is the key to the mosaic stack and is activated by supplying the calculated stabilities of the Fe-bearing minerals in the `predominant` argument.
This is a list whose elements correspond to each group of changing basis species given in the first argument.
The NULL means that the abundances of S-bearing aqueous species are calculated according to the default in `mosaic()` (actually using `equilibrate()` to blend their transitions).
Because the Fe-bearing minerals are the second group of changing basis species (`bases2`), their stabilities are given in the second position of the `predominant` list.
@@ -168,18 +337,18 @@
After that we add the legend and title.
-```{r series1_2, echo=7:14, results = "hide", message = FALSE, fig.width = 6, fig.height = 5, out.width = "80%", fig.align = "center", pngquant = pngquant}
+```{r stack1_2, echo=7:14, results = "hide", message = FALSE, fig.width = 6, fig.height = 5, out.width = "80%", fig.align = "center", pngquant = pngquant}
```
This diagram has a distinctive chalcopyrite "hook" that is controlled on the low-pH side by the reaction with pyrite.
Only that reaction is shown in many published diagrams [e.g. @And75;@Gio02], but diagrams with a similar chalcopyrite wedge or hook can be seen in @BBR77 and @Bri80.
-## Mosaic Series 2
+## Mosaic Stacking 2
-The results of a mosaic series can also processed with `flatten()` to label each region with the minerals from both systems.
+The results of a mosaic stack can also processed with `flatten()` to label each region with the minerals from both systems.
For this example, the speciation of aqueous sulfur is not considered.
-```{r series2, results = "hide", message = FALSE, fig.width = 6, fig.height = 4, out.width = "100%", pngquant = pngquant}
+```{r stack2, results = "hide", message = FALSE, fig.width = 6, fig.height = 4, out.width = "100%", pngquant = pngquant}
layout(matrix(c(1, 2, 3, 3), nrow = 2), widths = c(1, 1.5))
# Fe-S-O-H diagram
basis(c("copper", "hematite", "S2", "oxygen", "H2O"))
@@ -188,7 +357,7 @@
abbrv <- info(species()$ispecies)$abbrv
dFe <- diagram(aFe, names = abbrv)
# Cu-S-O-H diagram based on reactions with the
-# stable Fe-bearing minerals (mosaic series)
+# stable Fe-bearing minerals (mosaic stack)
species(c("copper", "chalcocite", "covellite", "chalcopyrite", "bornite"))
mCu <- mosaic(bases, S2 = c(-34, -10), O2 = c(-55, -40),
T = 125, predominant = list(dFe$predominant))
@@ -335,7 +504,7 @@
This example was prompted by Figure 3 of @MH85; earlier versions of the diagram are in @HBL69 and @Hel70a.
-In some ways this is like the inverse of the **mosaic series** example.
+In some ways this is like the inverse of the **mosaic stacking** example.
There, reactions were balanced on Fe or Cu, and *f*~O<sub>2</sub>~ and pH were used as plotting variables.
Here, the reactions are balanced on O~2~ and implicitly on H^+^ through the activity ratios with *a*~Fe^+2^~ and *a*~Cu^+^~, which are the plotting variables.
@@ -344,7 +513,7 @@
### Mosaic Combo
-Instead of adding minerals with different metals by stacking mosaic diagrams (**mosaic series**), it may be possible to include two different metals in the basis species and formed species.
+Instead of adding minerals with different metals by stacking mosaic diagrams, it may be possible to include two different metals in the basis species and formed species.
The `mosaic()` and `equilibrate()` functions can be combined to balance on two different elements.
The example here is for N and C rather than metals.
Modified: pkg/CHNOSZ/vignettes/vig.bib
===================================================================
--- pkg/CHNOSZ/vignettes/vig.bib 2020-07-22 03:11:31 UTC (rev 571)
+++ pkg/CHNOSZ/vignettes/vig.bib 2020-07-22 17:20:44 UTC (rev 572)
@@ -517,6 +517,17 @@
doi = {10.1039/FT9928800803},
}
+ at Article{SSWS97,
+ author = {Shock, Everett L. and Sassani, David C. and Willis, Marc and Sverjensky, Dimitri A.},
+ journal = {Geochimica et Cosmochimica Acta},
+ title = {{I}norganic species in geologic fluids: {C}orrelations among standard molal thermodynamic properties of aqueous ions and hydroxide complexes},
+ year = {1997},
+ number = {5},
+ pages = {907--950},
+ volume = {61},
+ doi = {10.1016/S0016-7037(96)00339-0},
+}
+
@Article{SS98,
author = {Shock, Everett L. and Schulte, Mitchell D.},
journal = {Journal of Geophysical Research},
@@ -582,6 +593,17 @@
doi = {10.2475/ajs.288.1.19},
}
+ at Article{WEP_82,
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/chnosz -r 572
More information about the CHNOSZ-commits
mailing list