multi-metal.Rmd

---
title: "Diagrams with multiple metals"
author: "Jeffrey M. Dick"
output:
  html_vignette:
    mathjax: null
vignette: >
  %\VignetteIndexEntry{Diagrams with multiple metals}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
bibliography: [vig.bib, OBIGT.bib]
csl: elementa.csl
link-citations: true
---
<style>
/* https://gomakethings.com/how-to-break-an-image-out-of-its-parent-container-with-css/ */
@media (min-width: 700px) {
  .full-width {
    left: 50%;
    margin-left: -50vw;
    margin-right: -50vw;
    max-width: 100vw;
    position: relative;
    right: 50%;
    width: 100vw;
  }
}
@media (min-width: 1020px) {
  .full-width {
    left: 50vw; /* fallback if needed */
    left: calc(50vw - 160px);
    width: 1000px;
    position: relative;
    background-color: #9ecff7;
    padding:10px;
  }
}
/* Zero margin around pre blocks (looks more like R console output) */
pre {
  margin-top: 0;
  margin-bottom: 0;
}
</style>
<script>
function ToggleDiv(ID) {
  var D = document.getElementById("D-" + ID);
  var B = document.getElementById("B-" + ID);
  if (D.style.display === "none") {
    // open the div and change button text
    D.style.display = "block";
    B.innerText = "Hide code";
  } else {
    // close the div and change button text
    D.style.display = "none";
    B.innerText = "Show code";
  }
}
</script>

```{r setup, include=FALSE}
options(width = 80)
## Use pngquant to optimize PNG images
library(knitr)
knit_hooks$set(pngquant = hook_pngquant)
pngquant <- "--speed=1 --quality=0-25"
if (!nzchar(Sys.which("pngquant"))) pngquant <- NULL

## Resolution settings
# Change this to TRUE to make high-resolution plots
# (default is FALSE to save time in CRAN checks)
hires <- nzchar(Sys.getenv("CHNOSZ_BUILD_LARGE_VIGNETTES"))
res1.lo <- 150
res1.hi <- 256
res1 <- if(hires) res1.hi else res1.lo
res2.lo <- 200
res2.hi <- 400
res2 <- if(hires) res2.hi else res2.lo

## logK with a thin space 20200627
logK <- "log&thinsp;<i>K</i>"

## Set dpi 20231129
knitr::opts_chunk$set(
  dpi = if(nzchar(Sys.getenv("CHNOSZ_BUILD_LARGE_VIGNETTES"))) 100 else 72
)
```

```{r CHNOSZ_reset, include=FALSE}
library(CHNOSZ)
reset()
```

This vignette was compiled on `r Sys.Date()` with CHNOSZ version `r sessionInfo()$otherPkgs$CHNOSZ$Version`.

This vignette is associated with a paper that has been published in *Applied Computing and Geosciences* ([Dick, 2021](https://doi.org/10.1016/j.acags.2021.100059 "Diagrams with multiple metals in CHNOSZ")).
Please consider citing that paper if you use the functions or diagrams described here.

The plots in this vignette were made using the following resolution settings.
Low resolutions are used for submitting the package to CRAN.
High resolutions are used if the `CHNOSZ_BUILD_LARGE_VIGNETTES` environment variable is set.

```{r res, results = "asis", echo = FALSE}
cat("```")
cat("\n")
cat(paste0(if(hires) "# " else "", "res1 <- ", res1.lo))
cat("\n")
cat(paste0(if(hires) "" else "# ", "res1 <- ", res1.hi))
cat("\n")
cat(paste0(if(hires) "# " else "", "res2 <- ", res2.lo))
cat("\n")
cat(paste0(if(hires) "" else "# ", "res2 <- ", res2.hi))
cat("\n")
cat("```")
```

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 **mashing**, **mixing**, **mosaic stacking**, and **secondary balancing**.

## Mashing

Mashing or simple overlay refers to independent calculations for two different systems that are displayed on the same diagram.

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 after the first are taken from the previous command.
This allows calculations to be run at the same conditions for a different system.
This feature is also used in other examples in this vignette.

```{r mash, echo = 1:8, eval = FALSE}
par(mfrow = c(1, 2))
basis("CHNOS+")
species(c("CH4", "CO2", "HCO3-", "CO3-2"))
aC <- affinity(pH = c(0, 14), O2 = c(-75, -60))
dC <- diagram(aC, dx = c(0, 1, 0, 0), dy = c(0, 1, 0, 0))
species(c("H2S", "HS-", "HSO4-", "SO4-2"))
aS <- affinity(aC)  # argument recall
dS <- diagram(aS, add = TRUE, col = 4, col.names = 4, dx = c(0, -0.5, 0, 0))
aCS <- mash(dC, dS)
srt <- c(0, 0, 90, 0, 0, 0, 90, 0, 0, 0)
cex.names <- c(1, 1, 0.8, 1, 1, 1, 1, 1, 1, 1)
dy <- c(0, 0, 0, -0.2, 0, 0, 0, 0, 0, 0)
diagram(aCS, srt = srt, cex.names = cex.names, dy = dy)
legend("topright", legend = lTP(25, 1), bty = "n")
```

The second diagram is just like the first, except the function `mash()` is used to label the fields with names of species from both systems, and a legend is added to indicate the temperature and pressure.

```{r mash, echo = 9:14,  results = "hide", message = FALSE, fig.width = 10, fig.height = 5, out.width = "100%"}
```

Note that these are predominance diagrams, so they show only the species with highest activity; there is in fact a distribution of activities of aqueous species that is not visible here.

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.

## Mixing 1

As shown above, mashing two diagrams is essentially a simple combination of the two 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.
Before getting started, it may be helpful to clarify some terminology.
In the materials science community, materials are characterized by several energies (among others): 1) the formation energy from the elements in their reference state, 2) the energy above the convex hull, which is zero for stable materials, and greater than zero for metastable materials, and 3) the Pourbaix energy difference (Δ*G*~pbx~), which refers to the energy of a given material with respect to the stable solids and aqueous species as a function of pH and Eh.
The parallel terminology used in CHNOSZ is that aqueous species or minerals have a 1) standard Gibbs energy of formation from the elements, Δ*G*° = *f*(*T*, *P*), which is available from the [OBIGT database](OBIGT.html), 2) standard Gibbs energy of reaction from the stable species, and 3) affinity of formation from the basis species, *A* = -Δ*G* = *f*(*T*, *P*, and activities of all species).
As used in CHNOSZ, "formation reaction" refers to formation from the basis species, not from the elements.
The basis species **are not** in general the stable species, so we begin by identifying the stable species in the system; the difference between *their* affinities and the affinity of any other species corresponds to -Δ*G*~pbx~.

First we need to assemble the standard Gibbs energies of the solids and aqueous species.
For solids, values of formation energy from the elements (in eV/atom) computed using density functional theory (DFT) were retrieved from the [Materials API](https://github.com/materialsproject/mapidoc) [@OCJ_15] and are converted to units of J/mol.
The Materials Project (MP) website also provides these values, but with fewer decimal places, which would lead to a small rounding error in the comparison of energy above the hull at the end of this example.
For aqueous species, values of standard Gibbs energy of formation from the elements 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 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 in CHNOSZ.

<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 Materials API, e.g.
# from pymatgen import MPRester
# m = MPRester("USER_API_KEY")
# m.query(criteria={"task_id": "mp-1279742"}, properties=["formation_energy_per_atom"])
# mp-13, mp-1279742, mp-19306, mp-19770
#Fe.cr <- c(Fe = 0, FeO = -1.72803, Fe3O4 = -1.85868, Fe2O3 = -1.90736)  # 20201109
Fe.cr <- c(Fe = 0, FeO = -1.72768, Fe3O4 = -1.85838, Fe2O3 = -1.90708)  # 20210219
# mp-146, mp-18937, mp-1275946, mp-19094, mp-756395, mp-25279
#V.cr <- c(V = 0, V2O3 = -2.52849, V3O5 = -2.52574, VO2 = -2.48496, V3O7 = -2.32836, V2O5 = -2.29524)  # 20201109
V.cr <- c(V = 0, V2O3 = -2.52787, V3O5 = -2.52516, VO2 = -2.48447, V3O7 = -2.32789, V2O5 = -2.29480)  # 20210219

# 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
# Most are from Wagman et al., 1982
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
)

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
)

# 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 standard 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
# This function modifies OBIGT and returns the species indices of the affected species
modfun <- function(x, state, model = NULL) sapply(seq_along(x), function(i) {
  # We explicitly set the units to Joules (this is the default as of CHNOSZ 2.0.0)
  if(is.null(model)) mod.OBIGT(names(x)[i], formula = names(x)[i], state = state, E_units = "J", G = x[i])
  else mod.OBIGT(names(x)[i], formula = names(x)[i], state = state, model = model, E_units = "J", G = x[i])
})
# We need model = "CGL" to override the Berman model for some minerals 20220919
iFe.cr <- modfun(Fe.cr, "cr", model = "CGL")
iFe.aq <- modfun(Fe.aq, "aq")
iV.cr <- modfun(V.cr, "cr", model = "CGL")
iV.aq <- modfun(V.aq, "aq")

# Formation energies (eV / atom) for bimetallic solids from Materials API
# mp-1335, mp-1079399, mp-866134, mp-558525, mp-504509 (triclinic FeVO4)
#FeV.cr <- c(FeV = -0.12928, FeV3 = -0.17128, Fe3V = -0.12854, Fe2V4O13 = -2.23833, FeVO4 = -1.75611)  # 20201109
FeV.cr <- c(FeV = -0.12815, FeV3 = -0.17114, Fe3V = -0.12832, Fe2V4O13 = -2.23967, FeVO4 = -1.75573)  # 20210219
# 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 produces 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 assign activities of aqueous species to 10^-5^, which is the default value for diagrams on the MP website (from the page for a material: "Generate Phase Diagram" -- "Aqueous Stability (Pourbaix)").
The following commands compute Eh-pH diagrams for the single-metal systems Fe-O-H and V-O-H.
The pH and Eh ranges are made relatively small in order to show just a part of the diagram.
The diagrams are not plotted, but the output of `diagram()` is saved in `dFe` and `dV` for later use.

```{r mixing1, eval = FALSE, echo = 1:14}
par(mfrow = c(1, 3))
loga.Fe <- -5
loga.V <- -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 = c(4, 10, res1), Eh = c(-1.5, 0, res1))
dFe <- diagram(aFe, plot.it = FALSE)
# 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, plot.it = FALSE)

# 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)
a11 <- mix(dFe, dV, dFeV, c(1, 1))
# Adjust labels 20210219
iV2O3 <- info("V2O3")
iFeO <- info("FeO", "cr")
iFe3V <- info("Fe3V")
srt <- rep(0, nrow(a11$species))
srt[a11$species$ispecies == paste(iFeO, iV2O3, sep = ",")] <- 90
srt[a11$species$ispecies == paste(iV2O3, iFe3V, sep = ",")] <- -13
diagram(a11, min.area = 0.01, srt = srt)
title("Fe:V = 1:1")
label.figure(lTP(25, 1), xfrac = 0.12)
# 1:3 mixture
a13 <- mix(dFe, dV, dFeV, c(1, 3))
srt <- rep(0, nrow(a13$species))
srt[a13$species$ispecies == paste(iFeO, iV2O3, sep = ",")] <- 90
srt[a13$species$ispecies == paste(iV2O3, iFe3V, sep = ",")] <- -13
diagram(a13, min.area = 0.01, srt = srt)
title("Fe:V = 1:3")
# 1:5 mixture
a15 <- mix(dFe, dV, dFeV, c(1, 5))
iFeV3 <- info("FeV3")
srt <- rep(0, nrow(a15$species))
srt[a15$species$ispecies == paste(iFeO, iV2O3, sep = ",")] <- 90
srt[a15$species$ispecies == paste(iV2O3, iFe3V, sep = ",")] <- -13
srt[a15$species$ispecies == paste(iV2O3, iFeV3, sep = ",")] <- -13
diagram(a15, min.area = 0.01, srt = srt)
title("Fe:V = 1:5")
```

Then we calculate the affinities for the bimetallic species and save the output of `diagram()` in `dFeV`, again without making a plot, but formatting the names in bold.
Note that `diagram()` uses different colors for regions with two solids, one solid, and no solids, including some transparency to show the underlying water stability region that is plotted first.

Now we have all the ingredients needed to combine the Fe-bearing, V-bearing, and bimetallic species to generate a given composition.
The `mix()` function is used to calculate the affinities of formation from basis species for all combinations of aqueous species and minerals that satisfy each of three different compositions.
Finally, the `diagram()`s are plotted; the `min.area` argument is used to remove labels for very small fields.
Regarding the legend, it should be noted that although the DFT calculations for solids are made for zero temperature and zero pressure [@SZS_17], the standard Gibbs energies of aqueous species [e.g. @WEP_82] are modified by a correction term so that they can be combined with DFT energies to reproduce the experimental energy for dissolution of a representative material for each metal at 25 °C and 1 bar [@PWLC12].

```{r mixing1, echo = 16:47, message = FALSE, results = "hide", fig.width = 9, fig.height = 3, 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 FeV~3~ + Fe~3~V assemblage is predicted to be stable instead of FeV.
This result is unlike Figure 1 of @SZS_17 but is consistent with the [MP page for FeV](https://doi.org/10.17188/1189535) where it is shown to decompose to this assemblage.
On the other hand, [FeV~3~ is stable](https://next-gen.materialsproject.org/materials/mp-1079399/) in the 1:3 mixture.
For an even higher proportion of V, the V + FeV~3~ assemblage is stable, which can be seen for instance in the Pourbaix diagram linked from the [MP page for FeV~5~O~12~](https://doi.org/10.17188/1305091).

Let's make another diagram for the 1:1 Fe:V composition over a broader range of Eh and pH.
The diagram shows a stable assemblage of Fe~2~O~3~ with an oxidized bimetallic material, [Fe~2~V~4~O~13~](https://next-gen.materialsproject.org/materials/mp-1200054/).
```{r FeVO4, eval = FALSE, echo = 1:29}
layout(t(matrix(1:3)), widths = c(1, 1, 0.2))
par(cex = 1)
# Fe-bearing species
basis(c("VO+2", "Fe+2", "H2O", "e-", "H+"))
species(c(iFe.aq, iFe.cr))$name
species(1:length(iFe.aq), loga.Fe)
aFe <- affinity(pH = c(0, 14, res2), Eh = c(-1.5, 2, res2))
dFe <- diagram(aFe, plot.it = FALSE)
# V-bearing species
species(c(iV.aq, iV.cr))$name
species(1:length(iV.aq), loga.V)
aV <- affinity(aFe)  # argument recall
dV <- diagram(aV, plot.it = FALSE)
# Bimetallic species
species(iFeV.cr)
aFeV <- affinity(aFe)  # argument recall
dFeV <- diagram(aFeV, plot.it = FALSE, bold = TRUE)
# 1:1 mixture (Fe:V)
a11 <- mix(dFe, dV, dFeV, c(1, 1))
# Adjust labels 20210219
iV2O3 <- info("V2O3")
iFe3V <- info("Fe3V")
iVO4m3 <- info("VO4-3")
iFe2O3 <- info("Fe2O3")
srt <- rep(0, nrow(a11$species))
srt[a11$species$ispecies == paste(iV2O3, iFe3V, sep = ",")] <- -13
srt[a11$species$ispecies == paste(iFe2O3, iVO4m3, sep = ",")] <- 90
d11 <- diagram(a11, min.area = 0.01, srt = srt)
water.lines(d11, col = "orangered")

# Calculate affinity of FeVO4
species("FeVO4")
aFeVO4 <- affinity(aFe)  # argument recall
# Calculate difference from stable species
aFeVO4_vs_stable <- aFeVO4$values[[1]] - d11$predominant.values
# Overlay lines from diagram on color map
diagram(a11, fill = NA, names = FALSE, limit.water = FALSE)
opar <- par(usr = c(0, 1, 0, 1))
col <- rev(topo.colors(128)) # No hcl.colors() in R < 3.6.0
if(getRversion() >= "3.6.0") col <- rev(hcl.colors(128, palette = "YlGnBu", alpha = 0.8))
image(aFeVO4_vs_stable, col = col, add = TRUE)
par(opar)
diagram(a11, fill = NA, add = TRUE, names = FALSE)
water.lines(d11, col = "orangered")
thermo.axis()

imax <- arrayInd(which.max(aFeVO4_vs_stable), dim(aFeVO4_vs_stable))
pH <- d11$vals$pH[imax[1]]
Eh <- d11$vals$Eh[imax[2]]
points(pH, Eh, pch = 10, cex = 2, lwd = 2, col = "gold")
stable <- d11$names[d11$predominant[imax]]
text(pH, Eh, stable, adj = c(0.5, -1), cex = 1.2, col = "gold")

# Make color scale 20210228
par(mar = c(3, 0, 2.5, 2.7))
plot.new()
levels <- 1:length(col)
plot.window(xlim = c(0, 1), ylim = range(levels), xaxs = "i", yaxs = "i")
rect(0, levels[-length(levels)], 1, levels[-1L], col = rev(col), border = NA)
box()
# To get the limits, convert range of affinities to eV/atom
arange <- rev(range(aFeVO4_vs_stable))
# This gets us to J/mol
Jrange <- convert(arange, "G")
# And to eV/atom
eVrange <- Jrange / 1.602176634e-19 / 6.02214076e23 / 6
ylim <- formatC(eVrange, digits = 3, format = "f")
axis(4, at = range(levels), labels = ylim)
mtext(quote(Delta*italic(G)[pbx]*", eV/atom"), side = 4, las = 0, line = 1)
```

We then compute the affinity for formation of a metastable material, in this case triclinic FeVO~4~, from the same basis species used to make the previous diagrams.
Given the diagram for the stable Fe-, V- and bimetallic materials *mixed with the same stoichiometry* as FeVO~4~ (1:1 Fe:V), the difference between their affinities of formation and that of FeVO~4~ corresponds to the Pourbaix energy difference (-Δ*G*~pbx~).
This is plotted as a color map in the second diagram.
(See the source of this vignette for the code used to make the scale bar.)

```{r FeVO4, echo = 31:44, message = FALSE, results = "hide", fig.width = 11, fig.height = 5, out.width = "100%", pngquant = pngquant}
```

Now we locate the pH and Eh that maximize the affinity (that is, minimize Δ*G*~pbx~) of FeVO~4~ compared to the stable species.
In agreement with @SZS_17, this is in the stability field of Fe~2~O~3~ + Fe~2~V~4~O~13~.

```{r Gpbx_min, echo = 2:8, message = FALSE, fig.keep = "none"}
plot(1:10) # so we can run "points" in this chunk
imax <- arrayInd(which.max(aFeVO4_vs_stable), dim(aFeVO4_vs_stable))
pH <- d11$vals$pH[imax[1]]
Eh <- d11$vals$Eh[imax[2]]
points(pH, Eh, pch = 10, cex = 2, lwd = 2, col = "gold")
stable <- d11$names[d11$predominant[imax]]
text(pH, Eh, stable, adj = c(0.3, 2), cex = 1.2, col = "gold")
(Apbx <- range(aFeVO4_vs_stable[d11$predominant == d11$predominant[imax]]))
```

Although one point is drawn on the diagram, FeVO~4~ has the same Pourbaix energy difference with respect to the entire Fe~2~O~3~ + Fe~2~V~4~O~13~ field, as shown by the `range()` command (the values are dimensionless values of affinity, *A*/(*RT*) = -Δ*G*~pbx~/(*RT*)).
This can occur only if the decomposition reaction has no free O~2~ or H~2~, and means that in this case Δ*G*~pbx~ in the Fe~2~O~3~ + Fe~2~V~4~O~13~ field is equal to the energy above the hull.

To calculate the energy above the hull "by hand", let's set up the basis species to be the stable decomposition products we just found.
O~2~ is also needed to make a square stoichiometric matrix (i.e. same number of elements and basis species), but it does not appear in the reaction to form FeVO~4~ from the basis species.
`subcrt()` is used to automatically balance the formation reaction for 1 mole of FeVO~4~ and calculate the standard Gibbs energy of the reaction.
We first test that `r logK` of the reaction (calculated with `convert()`, which divides Δ*G*° by *-RT*) is the same as the dimensionless affinity for FeVO~4~ calculated above.
Then, the value of Δ*G*° in J/mol is converted to eV/mol, and finally eV/atom.

```{r hull, message = FALSE}
b <- basis(c("Fe2O3", "Fe2V4O13", "O2"))
J_mol <- subcrt("FeVO4", 1, T = 25)$out$G
stopifnot(all.equal(rep(convert(J_mol, "logK"), 2), Apbx))
eV_mol <- J_mol / 1.602176634e-19
eV_atom <- eV_mol / 6.02214076e23 / 6
round(eV_atom, 3)
stopifnot(round(eV_atom, 3) == 0.415)
```

This is equal to the value for the energy above the hull / atom for [triclinic FeVO~4~ on the MP website](https://next-gen.materialsproject.org/materials/mp-504509/) (0.415 eV, accessed on 2020-11-09 and 2021-02-19).
This shows that we successfully made a round trip starting with the input formation energies (eV/atom) from the Materials API, to standard Gibbs energy (J/mol) in the OBIGT database, and back out to energy above the hull (eV/atom).

The concept of using the stable minerals and aqueous species to calculate reaction energetics is formalized in the `mosaic()` function, which is described next.
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

*For a function that implements the workflow described below, see `stack_mosaic()` (added in CHNOSZ 2.0.0).*

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 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.
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 `Fe.cr`.

```{r stack1_1, results = "hide", message = FALSE}
logaH2S <- -2
T <- 200
pH <- c(0, 14, res2)
O2 <- c(-48, -33, res2)
basis(c("Cu+", "pyrite", "H2S", "oxygen", "H2O", "H+"))
basis("H2S", logaH2S)
S.aq <- c("H2S", "HS-", "HSO4-", "SO4-2")
Fe.cr <- c("pyrite", "pyrrhotite", "magnetite", "hematite")
Fe.abbrv <- c("Py", "Po", "Mag", "Hem")
```

Now we calculate affinities for minerals in the Fe-S-O-H system that take account of the changing aqueous sulfur species in `S.aq`.
The result is used to make different layers of the diagram (1 and 2 are both made by the first call to `diagram()`):

1. Water stability region (gray shading)
2. Predominance fields for the aqueous S species (blue text and dashed lines)
3. Stability areas for the Fe-bearing minerals (black text and lines)

```{r stack1_2, eval = FALSE, echo = 1:4}
species(Fe.cr)
mFe <- mosaic(S.aq, pH = pH, O2 = O2, T = T)
diagram(mFe$A.bases, lty = 2, col = 4, col.names = 4, italic = TRUE, dx = c(0, 1, 0, 0), dy = c(-1.5, 0, 1, 0))
dFe <- diagram(mFe$A.species, add = TRUE, lwd = 2, names = Fe.abbrv, dx = c(0, 0.5, 0, 0), dy = c(-1, 0, 0.5, 0))
FeCu.cr <- c("chalcopyrite", "bornite")
Cu.cr <- c("copper", "cuprite", "tenorite", "chalcocite", "covellite")
FeCu.abbrv <- c("Ccp", "Bn", "Cu", "Cpr", "Tnr", "Cct", "Cv")
species(c(FeCu.cr, Cu.cr))
mFeCu <- mosaic(list(S.aq, Fe.cr), pH = pH, O2 = O2,
              T = T, stable = list(NULL, dFe$predominant))
col <- c(rep("#FF8C00", 2), rep(2, 5))
lwd <- c(2, 1, 1, 1, 1, 1, 1)
dy = c(0, 0, 0, 0, 0, 1, 0)
diagram(mFeCu$A.species, add = TRUE, col = col, lwd = lwd, col.names = col, bold = TRUE, names = FeCu.abbrv, dy = dy)
TPS <- c(describe.property(c("T", "P"), c(T, "Psat")), expression(sum(S) == 0.01*m))
legend("topright", TPS, bty = "n")
title("Cu-Fe-S-O-H (minerals only)", font.main = 1)
```

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 stack and is activated by supplying the calculated stabilities of the Fe-bearing minerals in the `stable` 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()`, which uses `equilibrate()` to compute the continuous transition between them ("blending").
Because the Fe-bearing minerals are the second group of changing basis species (`Fe.cr`), their stabilities are given in the second position of the `stable` list.
The result is used to plot the last layer of the diagram:

4. Stability areas for Cu-bearing minerals (red text and lines; orange for chalcopyrite and bornite)

After that we add the legend and title.

```{r stack1_2, echo=5:17, results = "hide", message = FALSE, fig.width = 6, fig.height = 5, out.width = "75%", fig.align = "center", pngquant = pngquant}
```

This diagram has a distinctive  chalcopyrite "hook" surrounded by a thin bornite field.
Only the chalcopyrite-bornite reaction in the pyrite field is shown in some published diagrams [e.g. @And75;@Gio02], but diagrams with a similar chalcopyrite wedge or hook geometry can be seen in @BBR77 and @Bri80.

## Mosaic Stacking 2

The previous diagram shows the relative stabilities of minerals only.
The next diagram adds aqueous species to the system.
The position of the boundaries between the stability fields for minerals and aqueous species are calculated for a given activity for the latter, in this case 10^-6^.

<button id="B-stack2" onclick="ToggleDiv('stack2')">Show code</button>
<div id="D-stack2" style="display: none">
```{r stack2, eval = FALSE}
# Define system
pH <- c(0, 14, res2)
O2 <- c(-48, -33, res2)
T <- 200
logmS <- -2
m_NaCl <- 0.1
logm_aq <- -6 # for both Fe- and Cu-bearing aq species
# Basis species
S.aq <- c("H2S", "HS-", "HSO4-", "SO4-2")
# Minerals
Fe.cr <- c("pyrite", "pyrrhotite", "magnetite", "hematite")
Fe.abbrv <- c("Py", "Po", "Mag", "Hem")
FeCu.cr <- c("chalcopyrite", "bornite")
Cu.cr <- c("copper", "cuprite", "tenorite", "chalcocite", "covellite")
FeCu.abbrv <- c("Ccp", "Bn", "Cu", "Cpr", "Tnr", "Cct", "Cv")
# Aqueous species
iFe.aq <- retrieve("Fe", c("S", "O", "H", "Cl"), "aq")
Fe.aq <- info(iFe.aq)$name
iCu.aq <- retrieve("Cu", c("S", "O", "H", "Cl"), "aq")
Cu.aq <- info(iCu.aq)$name
# Expressions for making the legend
TPexpr <- describe.property(c("T", "P"), c(T, "Psat"))
Sexpr <- as.expression(bquote(sum(S) == .(10^logmS)*m))
NaClexpr <- as.expression(bquote(NaCl == .(m_NaCl)*m))
aqexpr <- as.expression(bquote("("*aq*")"[italic(i)] == 10^.(logm_aq)*m))

# Setup basis species
basis(c("Cu+", "pyrite", "H2S", "oxygen", "H2O", "H+", "Cl-"))
basis("H2S", logmS)
nacl <- NaCl(m_tot = m_NaCl, T = T, P = "Psat")
basis("Cl-", log10(nacl$m_Cl))
# Fe-bearing minerals
species(Fe.cr)
# Add aqueous species 20210220
species(iFe.aq, logm_aq, add = TRUE)
mFe <- mosaic(S.aq, pH = pH, O2 = O2, T = T, IS = nacl$IS)
# Start plot with just the fields for transparency effect
dFe <- diagram(mFe$A.species, lwd = 0, names = FALSE)

# Cu-bearing minerals
species(c(FeCu.cr, Cu.cr))
# Add aqueous species 20210220
species(iCu.aq, logm_aq, add = TRUE)

## Mosaic with all Fe species as basis species
#mFeCu <- mosaic(list(S.aq, c(Fe.cr, Fe.aq)), pH = pH, O2 = O2, T = T, stable = list(NULL, dFe$predominant))
# Use only predominant Fe species as basis species (to speed up calculation) 20210224
predom <- dFe$predominant
ipredom <- sort(unique(as.numeric(predom)))
for(i in seq_along(ipredom)) predom[dFe$predominant == ipredom[i]] <- i
Fe.predom <- c(Fe.cr, Fe.aq)[ipredom]
# Use loga_aq argument to control the activity of aqueous species in mosaic calculation 20220722
# c(NA, logm_aq) means to use:
#   basis()'s value for logact of aqueous S species
#   logm_aq for logact of aqueous Fe species
mFeCu <- mosaic(list(S.aq, Fe.predom), pH = pH, O2 = O2, T = T, stable = list(NULL, predom), IS = nacl$IS, loga_aq = c(NA, logm_aq))

# Adjust labels
bold <- c(rep(TRUE, length(FeCu.abbrv)), rep(FALSE, length(Cu.aq)))
names <- c(FeCu.abbrv, Cu.aq)
srt <- dx <- dy <- rep(0, length(names))
cex <- rep(1, length(names))
dx[names == "Cu"] <- -1.5
dx[names == "Bn"] <- 1.4
dx[names == "CuHS"] <- 1
dx[names == "Cu+"] <- -0.5
dy[names == "Cu"] <- 3
dx[names == "Cct"] <- -2
dy[names == "Cct"] <- 4
dy[names == "CuHS"] <- 1
dy[names == "Bn"] <- -0.9
dx[names == "CuCl2-"] <- -1
dy[names == "CuCl2-"] <- 2
cex[names == "Bn"] <- 0.8
srt[names == "Bn"] <- 85
col.names <- col <- rep(2, nrow(mFeCu$A.species$species))
# Highlight Ccp and Bn in orange
col[1:2] <- "#FF8C00"
col.names[1:2] <- "#FF8C00"
# Thick line around Ccp field
lwd <- rep(1, nrow(mFeCu$A.species$species))
lwd[1] <- 2
diagram(mFeCu$A.species, add = TRUE, lwd = lwd, col = col, col.names = col.names, names = names, bold = bold, dx = dx, dy = dy, cex.names = cex, srt = srt)
# Add second Cu label
text(12.3, -47, "Cu", col = 2, font = 2)

# Plot the Fe-system lines and names "on top" so they are not covered by fill colors
diagram(mFe$A.bases, add = TRUE, lty = 2, col = 4, names = FALSE, fill = NA)
bold <- c(rep(TRUE, length(Fe.abbrv)), rep(FALSE, length(Fe.aq)))
names <- c(Fe.abbrv, Fe.aq)
srt <- dx <- dy <- rep(0, length(names))
cex <- rep(1, length(names))
dy[names == "Hem"] <- 0.5
dy[names == "Mag"] <- -0.8
dx[names == "Hem"] <- -0.5
dx[names == "Mag"] <- 0.25
dx[names == "Fe+2"] <- 0.5
dx[names == "FeO2-"] <- 1
dy[names == "FeO2-"] <- -3
dx[names == "HFeO2-"] <- -0.5
dy[names == "HFeO2-"] <- 1
srt[names == "Mag"] <- 83
srt[names == "FeSO4"] <- 90
diagram(mFe$A.species, add = TRUE, lwd = 2, names = names, bold = bold, dx = dx, dy = dy, cex.names = cex, srt = srt, fill = NA)
legend("topright", c(TPexpr, Sexpr, NaClexpr, aqexpr), bty = "n")
title("Cu-Fe-S-O-H-Cl (minerals and aqueous species)", font.main = 1)

# Restore default OBIGT database
OBIGT()
```
</div>

```{r stack2, echo = FALSE, results = "hide", message = FALSE, fig.width = 6, fig.height = 5, out.width = "75%", fig.align = "center", pngquant = pngquant}
```

After running the code to make this diagram, we can list the reference keys for the minerals and aqueous species.

```{r stack2.refs, message = FALSE}
minerals <- list(Fe.cr = Fe.cr, Cu.cr = Cu.cr, FeCu.cr = FeCu.cr)
aqueous <- list(S.aq = S.aq, Fe.aq = Fe.aq, Cu.aq = Cu.aq)
allspecies <- c(minerals, aqueous)
iall <- lapply(allspecies, info)
allkeys <- lapply(iall, function(x) thermo.refs(x)$key)
allkeys
```

The next code chunk prepends `@` to the reference keys and uses the chunk option [`results = 'asis'`](https://bookdown.org/yihui/rmarkdown-cookbook/results-asis.html) to insert the citations into the R Markdown document in chronological order.

```{r stack2.cite, results = "asis"}
allyears <- lapply(iall, function(x) thermo.refs(x)$year)
o <- order(unlist(allyears))
keys <- gsub("\\..*", "", unique(unlist(allkeys)[o]))
cat(paste(paste0("@", keys), collapse = "; "))
```

## Mixing 2

The previous diagram shows a stability boundary between chalcopyrite and bornite but does not identify the stable *assemblages* that contain these minerals.
This is where `mix()` can help.
Following the workflow described in **Mixing 1**, we first calculate individual diagrams for Fe-S-O-H and Cu-S-O-H, which are overlaid on the first plot and saved in `dFe` and `dCu`.
We then calculate the affinities for the bimetallic Cu and Fe minerals and run them through `diagram()` without actually making a plot, but save the result in `dFeCu`.
Then, we combine the results using `mix()` to define different proportions of Fe and Cu.

<button id="B-mixing2" onclick="ToggleDiv('mixing2')">Show code</button>
<div id="D-mixing2" style="display: none">
```{r mixing2, eval = FALSE}
par(mfrow = c(2, 2))

logaH2S <- -2
T <- 200
pH <- c(0, 14)
O2 <- c(-60, -25)
basis(c("Cu+", "Fe+2", "H2S", "oxygen", "H2O", "H+"))
basis("H2S", logaH2S)

S.aq <- c("H2S", "HS-", "HSO4-", "SO4-2")
Fe.cr <- c("pyrite", "pyrrhotite", "magnetite", "hematite")
Cu.cr <- c("copper", "cuprite", "tenorite", "chalcocite", "covellite")
FeCu.cr <- c("chalcopyrite", "bornite")

species(Fe.cr)
mFe <- mosaic(S.aq, pH = pH, O2 = O2, T = T)
diagram(mFe$A.bases, lty = 2, col = 4, col.names = 4, italic = TRUE, dx = c(0, 1, 0, 0), dy = c(-1, 0, -2, 0))
names <- info(info(Fe.cr))$abbrv
dFe <- diagram(mFe$A.species, add = TRUE, names = names)

species(Cu.cr)
mCu <- mosaic(S.aq, pH = pH, O2 = O2, T = T)
names <- info(info(Cu.cr))$abbrv
col.names <- rep(2, length(names))
col.names[1] <- 0
dCu <- diagram(mCu$A.species, add = TRUE, col = 2, col.names = col.names, names = names)
text(12, -55, "Cu", col = 2)
legend("topright", legend = lTP(T, "Psat"), bty = "n")
title(paste("Fe-S-O-H and Cu-S-O-H; Total S =", 10^logaH2S, "m"))

species(FeCu.cr)
mFeCu <- mosaic(S.aq, pH = pH, O2 = O2, T = T)
names <- info(info(FeCu.cr))$abbrv
dFeCu <- diagram(mFeCu$A.species, plot.it = FALSE, names = names)

fill <- function(a) {
  ifelse(grepl("Ccp", a$species$name), "#FF8C0088",
    ifelse(grepl("Bn", a$species$name), "#DC143C88", NA)
)}
srt <- function(a) ifelse(a$species$name %in% c("Mag+Bn", "Mag+Ccp"), 80, 0)
cex.names <- function(a) ifelse(a$species$name %in% c("Mag+Bn", "Mag+Ccp", "Mag+Cct"), 0.8, 1)
dx <- function(a) sapply(a$species$name, switch, "Mag+Bn" = 0.15, "Mag+Cct" = 0.5, 0)
dy <- function(a) sapply(a$species$name, switch, "Py+Bn" = -1, "Po+Bn" = -0.8, "Po+Cu" = -0.8, "Mag+Ccp" = -1, 0)

a11 <- mix(dFe, dCu, dFeCu)
diagram(a11, fill = fill(a11), srt = srt(a11), min.area = 0.01, dx = dx(a11), dy = dy(a11), cex.names = cex.names(a11))
title("Fe:Cu = 1:1")

a21 <- mix(dFe, dCu, dFeCu, c(2, 1))
diagram(a21, fill = fill(a21), srt = srt(a21), min.area = 0.01, dx = dx(a21), dy = dy(a21), cex.names = cex.names(a21))
title("Fe:Cu = 2:1")

a12 <- mix(dFe, dCu, dFeCu, c(1, 2))
diagram(a12, fill = fill(a12), srt = srt(a12), min.area = 0.01, dx = dx(a12), dy = dy(a12), cex.names = cex.names(a12))
title("Fe:Cu = 1:2")
```
</div>

```{r mixing2, echo = FALSE, results = "hide", message = FALSE, fig.width = 8, fig.height = 6.5, out.width = "100%", pngquant = pngquant}
```

These diagrams show that changing the amounts of the metals affects the stability of minerals involved in reactions with chalcopyrite.
At a 1:1 ratio of Fe:Cu, chalcopyrite is a stable single-mineral assemblage.
At a 2:1 ratio, pyrite, pyrrhotite, or magnetite can coexist in a two-phase assemblage with chalcopyrite.
At a 1:2 ratio, an assemblage consisting of the two bimetallic minerals (chalcopyrite and bornite) is stable.

## Mosaic Stacking 3

The results of a mosaic stack can also be processed with `mash()` to label each region with the minerals from both systems.
For this example, the speciation of aqueous sulfur is not considered; instead, the fugacity of S~2~ is a plotting variable.
The stable Fe-bearing minerals (Fe-S-O-H) are used as the changing basis species to make the diagram for Cu-bearing minerals (Cu-Fe-S-O-H).
Then, the two diagrams are mashed to show all minerals in a single diagram.
Greener colors are used to indicate minerals with less S~2~ and more O~2~ in their formation reactions.

<button id="B-stack3" onclick="ToggleDiv('stack3')">Show code</button>
<div id="D-stack3" style="display: none">
```{r stack3, eval = FALSE}
T <- 125
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", "H+", "Cl-"))
bFe <- species(c("hematite", "magnetite", "pyrite"))$name
aFe <- affinity(S2 = c(-34, -10, res1), O2 = c(-55, -40, res1), T = T)
# Order species by a function of composition to make colors
oFe <- order(aFe$species$S2 - aFe$species$O2)
fill <- terrain.colors(length(oFe), alpha = 0.2)[oFe]
abbrv <- info(aFe$species$ispecies)$abbrv
dFe <- diagram(aFe, names = abbrv, fill = fill)
title("Fe-S-O-H")

# Cu-Fe-S-O-H diagram based on reactions with the
# stable Fe-bearing minerals (mosaic stack)
bCu <- species(c("copper", "chalcocite", "covellite", "chalcopyrite", "bornite"))$name
mCu <- mosaic(bFe, S2 = c(-34, -10, res1), O2 = c(-55, -40, res1),
              T = T, stable = list(dFe$predominant))
oCu <- order(mCu$A.species$species$S2 - mCu$A.species$species$O2)
fill <- terrain.colors(length(oCu), alpha = 0.2)[oCu]
abbrv <- info(mCu$A.species$species$ispecies)$abbrv
dCu <- diagram(mCu$A.species, names = abbrv, fill = fill, dx = c(0, 0, 0, 0, 1.8))
title("Cu-Fe-S-O-H")

# Mash the diagrams and adjust labels
aFeCu <- mash(dFe, dCu)
names <- aFeCu$species$name
dx <- dy <- srt <- rep(0, length(names))
cex <- rep(1, length(names))
cex[names %in% c("Hem+Ccp", "Hem+Cv")] <- 0.8
srt[names %in% c("Mag+Cu", "Hem+Cu")] <- 90
srt[names %in% c("Mag+Bn", "Hem+Bn")] <- 63
srt[names %in% c("Mag+Ccp", "Hem+Ccp")] <- 68
srt[names %in% c("Py+Bn", "Py+Cv")] <- 90
dx[names == "Hem+Ccp"] <- -0.4
dy[names == "Hem+Ccp"] <- -0.5
oFeCu <- order(aFeCu$species$S2 - aFeCu$species$O2)
fill <- terrain.colors(length(oFeCu), alpha = 0.2)[oFeCu]
diagram(aFeCu, cex.names = cex, srt = srt, dx = dx, dy = dy, fill = fill)
legend("topleft", legend = lTP(T, "Psat"), bg = "white")
title("Cu-Fe-S-O-H")
```
</div>

```{r stack3, echo = FALSE, results = "hide", message = FALSE, fig.width = 6, fig.height = 4, out.width = "100%", pngquant = pngquant}
```

The resulting diagram is similar to Figure 2 of @Sve87; that diagram also shows calculations of the solubility of Cu and concentration of SO~4~^-2^ in model Cu ore-forming fluids.
The `solubility()` function can be used to calculate the total concentration of Cu in different complexes in solution (listed in the `iaq` argument).
The `bases` argument triggers a `mosaic()` calculation, so that the solubility corresponds that that of stable minerals at each point on the diagram.
The pH for these calculations is set to 6, and the molality of free Cl^-^, which affects the formation of the Cu chloride complexes, is estimated based on the composition of fluids from Table 2 of @Sve87 (ca. 80000 mg Cl / kg H~2~O) and the `NaCl()` function in CHNOSZ.
This also gives an estimated ionic strength, which is used in the following `mosaic()` and `affinity()` calls to calculate activity coefficients.

<button id="B-solubility" onclick="ToggleDiv('solubility')">Show code</button>
<div id="D-solubility" style="display: none">
```{r solubility, eval = FALSE}
par(mfrow = c(1, 3))
basis("pH", 6)
# Estimate the molality of Cl for ca. 80,000 mg/kg solution (Table 2 of Sverjensky, 1987)
m_tot <- 80000 / mass("Cl") / 1000
calc <- NaCl(T = T, m_tot = m_tot)
# Use log molality here, not log activity, because
# activity coefficients are calculated by setting IS below
basis("Cl-", log10(calc$m_Cl))
# Initial setup: dissolve a single mineral to form aqueous Cu complexes
species("copper")
iaq <- retrieve("Cu", c("O", "H", "Cl", "S"), "aq")

# Function to calculate solubility of Cu for stable assemblages of Fe and Cu minerals
# (i.e. equilibrium is imposed with all of these minerals, not only Cu(s))
mfun <- function() {
  s <- solubility(iaq, bases = list(bFe, bCu), S2 = c(-34, -10, res1), O2 = c(-55, -40, res1),
    T = T, IS = calc$IS, stable = list(dFe$predominant, dCu$predominant))
  s <- convert(s, "ppm")
  diagram(aFeCu, names = NA, col = "gray", fill = fill)
  diagram(s, type = "loga.balance", levels = 10^(-3:3), add = TRUE)
  diagram(s, type = "loga.balance", levels = 35, add = TRUE, lwd = 3, col = 6, contour.method = NA)
}
# DIAGRAM 1
mfun()
title("Cu (ppm)")

# Calculate logK for CuCl2- dissociation at 125 °C
logK <- subcrt(c("CuCl2-", "Cu+", "Cl-"), c(-1, 1, 2), T = T)$out$logK
# Sverjensky (1987) used Helgeson (1969) value, which is ca. -5.2
dlogK <- logK - -5.2
# Calculate the difference in ΔG° corresponding to this logK difference
dG_J <- convert(dlogK, "G", T = convert(T, "K"))
# We should use calories here because the database values are in calories 20220604
stopifnot(info(info("CuCl2-"))$E_units == "cal")
dG_cal <- convert(dG_J, "cal")
# Apply this difference to the CuCl2- entry in OBIGT
newG <- info(info("CuCl2-"))$G + dG_cal
mod.OBIGT("CuCl2-", G = newG)

# Do the same thing for CuCl3-2
logK <- subcrt(c("CuCl3-2", "Cu+", "Cl-"), c(-1, 1, 3), T = T)$out$logK
dlogK <- logK - -5.6
dG_J <- convert(dlogK, "G", T = convert(T, "K"))
stopifnot(info(info("CuCl3-2"))$E_units == "cal")
dG_cal <- convert(dG_J, "cal")
newG <- info(info("CuCl3-2"))$G + dG_cal
mod.OBIGT("CuCl3-2", G = newG)

# DIAGRAM 2
mfun()
title("Cu (ppm)", line = 1.7)
CuCl2 <- expr.species("CuCl2-")
CuCl3 <- expr.species("CuCl3-2")
title(bquote("Helgeson (1969)"~.(CuCl2)~and~.(CuCl3)), line = 0.9)

# Set up system to dissolve S2(gas)
basis(c("S2", "copper", "hematite", "oxygen", "H2O", "H+", "Cl-"))
basis("pH", 6)
species("S2")
# Calculate concentration of SO4-2
iaq <- info("SO4-2")
s <- solubility(iaq, S2 = c(-34, -10, res1), O2 = c(-55, -40, res1), T = T, IS = calc$IS, in.terms.of = "SO4-2")
s <- convert(s, "ppm")
# DIAGRAM 3
diagram(aFeCu, names = NA, col = "gray", fill = fill)
diagram(s, type = "loga.balance", levels = 10^(-3:3), add = TRUE)
diagram(s, type = "loga.balance", levels = 35, add = TRUE, lwd = 3, col = 6, contour.method = NA)
title(bquote(bold(.(expr.species("SO4-2"))~"(ppm)")))
```
</div>

```{r solubility, echo = FALSE, results = "hide", message = FALSE, fig.width = 7, fig.height = 3, out.width = "100%", fig.align = "center", pngquant = pngquant}
```

After running the code above, we can inspect the value of `calc` to show the estimated ionic strength and activity of Cl^-^; the latter is very close to unity.
```{r NaCl}
# Ionic strength
calc$IS
# Logarithm of activity of Cl-
log10(calc$m_Cl * calc$gam_Cl)
```

The thick magenta lines indicate the 35 ppm contour for Cu and SO~4~^-2^.
The first plot shows a lower Cu solubility in this region compared to Figure 2 of @Sve87.
The difference is likely due to lower stabilities of Cu(I) chloride complexes in the default OBIGT database, compared to those available at the time [@Hel69].
For the second plot, the standard Gibbs energies of CuCl~2~^-^ and CuCl~3~^-2^ are adjusted so that the `logK` for their dissociation reactions at 125 °C matches values interpolated from Table 5 of @Hel69.
Here are the `logK` values after the adjustment, followed a `reset()` call to compare the values with the default database, which is also used for later examples in this vignette.
(`T` was set to 125 above.)
```{r logK, message = FALSE}
# logK values interpolated from Table 5 of Helgeson (1969)
subcrt(c("CuCl2-", "Cu+", "Cl-"), c(-1, 1, 2), T = T)$out$logK
subcrt(c("CuCl3-2", "Cu+", "Cl-"), c(-1, 1, 3), T = T)$out$logK
# Default OBIGT database
reset()
subcrt(c("CuCl2-", "Cu+", "Cl-"), c(-1, 1, 2), T = T)$out$logK
subcrt(c("CuCl3-2", "Cu+", "Cl-"), c(-1, 1, 3), T = T)$out$logK
```

The higher stability of these complexes from @Hel69 causes the 35 ppm contour to move closer to the position shown by @Sve87.

Interestingly, the calculation here also predict substantial increases of Cu concentration of Cu at high *f*~S<sub>2</sub>~ and low *f*~O<sub>2</sub>~, due to the formation of bisulfide complexes with Cu.
The aqueous species considered in the calculation can be seen like this:
```{r iCu.aq}
names(iCu.aq)
```

CuHS and Cu(HS)~2~^-^ can be excluded by removing S from the `retrieve()` call above (i.e. only `c("O", "H", "Cl")` as the elements in possible ligands); doing so precludes a high concentration of aqueous Cu in the highly reduced, sulfidic region.

The third plot for the concentation of SO~4~^-2^ is simply made by using `affinity()` to calculate the affinity of its formation reaction as a function of *f*~S<sub>2</sub>~ and *f*~O<sub>2</sub>~ at pH 6 and 125 °C, then using `solubility()` to calculate the solubility of S~2~(gas), expressed in terms of moles of SO~4~^-2^ in order to calculate parts per million (ppm) by weight.

## Secondary Balancing

Predominance diagrams in CHNOSZ are made using the maximum affinity method, where the affinities of formation reactions of species are divided by the *balancing coefficients* [@Dic19].
Usually, these balancing coefficients are taken from the formation reactions themselves; for example, if they are the coefficients on the Fe-bearing basis species, then the reactions are said to be "balanced on Fe".

Some diagrams in the literature are made with secondary balancing constraints in addition to the primary ones.
For example, reactions of Fe-bearing minerals are balanced on Fe, and reactions of Cu-bearing minerals are balanced on Cu; these are both primary balancing coefficients.
Then, reactions between all minerals are balanced on H^+^ as the secondary balancing coefficients.
The core concept is to apply the secondary balance while also maintaining the primary balance; a method to do this has been implemented in the `rebalance()` function.

Different parts of the script to make the diagrams are described below; press the button to show the entire script.

<button id="B-rebalance" onclick="ToggleDiv('rebalance')">Show code</button>
<div id="D-rebalance" style="display: none">
```{r rebalance, eval = FALSE}
mat <- matrix(c(1, 1, 2, 2, 3, 3, 4, 4, 4, 5, 5, 5), byrow = TRUE, nrow = 2)
layout(mat)
par(font.main = 1)

basis(c("Fe+2", "Cu+", "hydrogen sulfide", "oxygen", "H2O", "H+"))
xlab <- ratlab("Fe+2", "Cu+")

### PRIMARY balancing

# Only Fe-bearing minerals
species(c("pyrite", "pyrrhotite", "magnetite", "hematite"))
aFe <- affinity("Fe+2" = c(0, 12), O2 = c(-40, -16), T = 400, P = 2000)
names <- info(aFe$species$ispecies)$abbrv
dFe <- diagram(aFe, xlab = xlab, names = names)
title(bquote("Only Fe; 1° balance:" ~ .(expr.species(dFe$balance))))
label.plot("A")

# Only Cu-bearing minerals
species(c("covellite", "chalcocite", "tenorite", "cuprite"))
aCu <- affinity(aFe)  # argument recall
names <- info(aCu$species$ispecies)$abbrv
dCu <- diagram(aCu, xlab = xlab, names = names)
title(bquote("Only Cu; 1° balance:" ~ .(expr.species(dCu$balance))))
label.plot("B")

# Only Fe- AND Cu-bearing minerals
species(c("chalcopyrite", "bornite"))
aFeCu <- affinity(aFe)
names <- info(aFeCu$species$ispecies)$abbrv
dFeCu <- diagram(aFeCu, xlab = xlab, balance = "H+", names = names)
title(bquote("Only Fe+Cu; 1° balance:" ~ .(expr.species(dFeCu$balance))))
label.plot("C")

### SECONDARY balancing

# Fe- or Cu-bearing minerals
ad1 <- rebalance(dFe, dCu, balance = "H+")
names <- info(ad1$species$ispecies)$abbrv
d1 <- diagram(ad1, xlab = xlab, balance = 1, names = names)
title(bquote("Only Fe or Cu; 2° balance:" ~ .(expr.species("H+"))))
label.plot("D")

# All minerals
d1$values <- c(dFe$values, dCu$values)
ad2 <- rebalance(d1, dFeCu, balance = "H+")
names <- info(ad2$species$ispecies)$abbrv
diagram(ad2, xlab = xlab, balance = 1, names = names)
title(bquote("Fe and/or Cu; 2° balance:" ~ .(expr.species("H+"))))
label.plot("E")

db <- describe.basis(3)
leg <- lex(lTP(400, 2000), db)
legend("bottomleft", legend = leg, bty = "n")
```
</div>

We first define basis species to contain both Cu- and Fe-bearing species.
The \emph{x} axis is the ratio of activities of Fe^+2^ and Cu^+^; the label is made with `ratlab()`.
```{r rebalance, eval = FALSE, echo = 5:6}
```

We then calculate the diagrams for the primary balancing coefficients, for the groups of only Fe-, only Cu-, and only Fe+Cu-bearing minerals.
It is obvious that the first two systems are balanced on Fe and Cu, respectively, but the third has a somewhat unusual balance: H^+^.
See Reaction 4 of @MH85 for an example.
```{r rebalance, eval = FALSE, echo = 10:33}
```

Now comes the secondary balancing, where all reactions, not only that between bornite and chalcopyrite, are balanced on H^+^.
We first rebalance the diagrams for the Fe- or Cu-bearing minerals to make diagram D.
Note that after secondary balancing with `rebalance()`, the argument `balance = 1` should be used in `diagram()` to prevent further balancing.
This is because `rebalance()` preserves the primary balancing for Fe- and Cu-bearing minerals (internally the "plotvals" components of `dFe` and `dCu`).

Then we rebalance diagrams D and C to make the final diagram in E.
The fields in this diagram are labeled with mineral abbreviations from the OBIGT database.
```{r rebalance, eval = FALSE, echo = 36:49}
```

```{r rebalance, echo = FALSE, results = "hide", message = FALSE, fig.width = 6.5, fig.height = 5, out.width = "100%", fig.align = "center", pngquant = pngquant}
```

The final diagram is like one shown in Figure 5 of @Bri80 and Figure 5 of @MH85.

*Challenge*: Although the diagram here is drawn only for H~2~S in the basis species, take it a step further and make a mosaic diagram to account for the stability of HSO~4~^-^ at high oxygen fugacity.

## Other Possibilities

Conceptually, the methods described above treat different metal-bearing elements as parts of distinct chemical systems that are then joined together.
Other methods may be more suitable for considering multiple metals (or other elements) in one system.

### Balancing on a Non-Metal

As shown in the **secondary balancing** example, there is no requirement that the balancing coefficients come from a metal-bearing species.
It is possible to make diagrams for minerals with different metallic elements simply by using a non-metallic element as the primary balance.
Here is an example for the Cu-Fe-S-O-H system.
The reactions are balanced on O~2~, which means that no O~2~ appears in the reaction between any two minerals, but Fe^+2^ and/or Cu^+^ can be present, depending on the chemical composition.
Saturation limits are shown for species that have no O~2~ in their formation reactions.

```{r non-metal, results = "hide", message = FALSE, fig.width = 6, fig.height = 5, out.width = "80%", fig.align = "center", pngquant = pngquant}
basis(c("Fe+2", "Cu+", "hydrogen sulfide", "oxygen", "H2O", "H+"))
basis("H2S", 2)
species(c("pyrite", "magnetite", "hematite", "covellite", "tenorite",
          "chalcopyrite", "bornite"))
a <- affinity("Cu+" = c(-8, 2, 500), "Fe+2" = c(-4, 12, 500), T = 400, P = 2000)
names <- info(a$species$ispecies)$abbrv
d <- diagram(a, xlab = ratlab("Cu+"), ylab = ratlab("Fe+2"), balance = "O2", names = names)
title(bquote("Cu-Fe-S-O-H; 1° balance:" ~ .(expr.species(d$balance))))
# Add saturation lines
species(c("pyrrhotite", "ferrous-oxide", "chalcocite", "cuprite"))
asat <- affinity(a)  # argument recall
names <- asat$species$name
names[2] <- "ferrous oxide"
diagram(asat, type = "saturation", add = TRUE, lty = 2, col = 4, names = names)
legend("topleft", legend = lTP(400, 2000), bty = "n")
```

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 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.

More common diagrams of this type are balanced on Si or Al.
See `demo(saturation)` for an example in the H~2~O-CO~2~-CaO-MgO-SiO~2~ system.

### Mosaic Combo

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 compares two methods applied to N-, C-, and N+C-bearing species because bimetallic aqueous species are not currently available in the OBIGT database.
The total activities used here are modified from the example for sedimentary basin brines described by @Sho93, which is also the source of the thermodynamic parameters for acetamide.

1. The "effective equilibrium constant" (*K*^eff^) method [@RBG_21] is used to calculate the activity of acetamide for a total activities of neutral and ionized species, i.e. ∑(ammonia and ammonium) and ∑(acetic acid and acetate).

2. Using the mosaic combo method, the `mosaic()` command calculates equilibrium activities of NH~3~ and NH~4~^+^ for a given total activity of N in the basis species, and calculates the corresponding affinities of the formed species.
Then, the `equilibrate()` command calculates equilibrium activities of the formed species for given total activity of C, and combines them with the activities of the changing basis species (NH~3~ and NH~4~^+^).

The mosaic combo method (solid black line) produces results equivalent to those of the *K*^eff^ method (dashed blue line).

<button id="B-mosaic-combo" onclick="ToggleDiv('mosaic-combo')">Show code</button>
<div id="D-mosaic-combo" style="display: none">
```{r mosaic-combo, eval = FALSE}
## METHOD 1: Keff equation (Robinson et al., 2021)

# A function to calculate Keff for any combination of T and pH
Keff <- function(T = 25, pH = 7) {
  # Make T and pH the same length
  len <- max(length(T), length(pH))
  T <- rep(T, length.out = len)
  pH <- rep(pH, length.out = len)
  # Calculate activity of H+
  aH <- 10^(-pH)
  # Calculate logKs
  logK1 <- subcrt(c("acetic acid", "NH3", "acetamide", "H2O"), c(-1, -1, 1, 1), T = T)$out$logK
  logK2 <- subcrt(c("acetic acid", "acetate", "H+"), c(-1, 1, 1), T = T)$out$logK
  logK3 <- subcrt(c("NH3", "H+", "NH4+"), c(-1, -1, 1), T = T)$out$logK
  # Calculate Ks
  K1 <- 10^logK1
  K2 <- 10^logK2
  K3 <- 10^logK3
  # Calculate Keff (Eq. 7)
  Keff <- K1 * (1 + K2 / aH) ^ -1 * (1 + K3 * aH) ^ -1
  Keff
}

# Calculate logKeff as a function of pH at 100 °C
res <- 128
pH <- seq(0, 14, length.out = res)
T <- 100
logKeff <- log10(Keff(pH = pH, T = T))

# Calculate activity of acetamide for
# acetic acid + acetate = 0.01 m
# ammonia + ammonium = 0.001 m
logAc <- log10(0.01)
logAm <- log10(0.001)
logAcAm <- logKeff + logAc + logAm

## METHOD 2: Mosaic combo

# Define total activities
a_N <- 0.001
# This is 2 * 0.01 because acetic acid has 2 carbons
a_C <- 2 * 0.01
loga_N <- log10(a_N)
loga_C <- log10(a_C)
# Setup basis species
basis(c("CO2", "NH3", "O2", "H2O", "H+"))
basis("NH3", loga_N)
# Load all C-bearing species (including acetamide)
species(c("acetamide", "acetic acid", "acetate"))
# Calculate distribution of C-bearing species accounting for ammonia/ammonium speciation
m <- mosaic(c("NH3", "NH4+"), pH = c(0, 14, res), T = T)
e <- equilibrate(m, loga.balance = loga_C)

# Plot and label diagram
# Start with empty diagram
diagram(e, ylim = c(-8, -0), lty = 0, names = FALSE)
# Add pH = 6 line
abline(v = 6, col = "gray60", lty = 5)
# Add line for acetamide activity calculated with Keff
lines(pH, logAcAm, col = 4, lwd = 6, lty = 2)
# Add lines from CHNOSZ calculations
diagram(e, add = TRUE, lty = c(2, 3, 1, 2, 3), lwd = c(1, 1, 2, 1, 1), dx = c(-0.2, 0.2, -2.5, 0, 0), dy = c(0.1, 0.1, -1, 0.1, 0.1), srt = c(0, 0, 52, 0, 0))
tN <- paste("Total N in basis species =", format(a_N, scientific = FALSE), "m")
tC <- paste("Total C in formed species =", format(a_C, scientific = FALSE), "m")
title(main = paste(tN, tC, sep = "\n"), font.main = 1)
legend("topright", legend = lTP(T, "Psat"), bty = "n")
# Check that we got equal values
stopifnot(all.equal(as.numeric(e$loga.equil[[3]]), logAcAm, tol = 1e-3, scale = 1))
```
</div>

```{r mosaic-combo, echo = FALSE, results = "hide", message = FALSE, fig.width = 6, fig.height = 5, out.width = "80%", fig.align = "center", pngquant = pngquant}
```

The diagram shows the ionization of acetic acid and NH~3~ at different pHs.
The predicted appearance of acetamide (CH~3~CONH~2~) is a consequence of the interaction between the N-bearing and C-bearing species, and is analogous to the formation of a bimetallic complex.

*Thanks to Kirt Robinson for the feature request and test case used in this example.*

## Document History

* 2020-07-15 First version.
* 2021-03-01 Improve mineral abbreviations and placement of labels; use updated DFT energies from Materials Project; add Mosaic Stacking 2 (minerals and aqueous species); add *K*~eff~ calculation; add Δ*G*~pbx~ color scale.
* 2024-03-25 Use orange for both chalcopyrite and bornite in Mosaic Stacking 1 and 2

## References