[CHNOSZ-commits] r586 - in pkg/CHNOSZ: . man vignettes
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sun Jul 26 16:54:54 CEST 2020
Author: jedick
Date: 2020-07-26 16:54:53 +0200 (Sun, 26 Jul 2020)
New Revision: 586
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/man/mosaic.Rd
pkg/CHNOSZ/vignettes/mklinks.sh
pkg/CHNOSZ/vignettes/multi-metal.Rmd
pkg/CHNOSZ/vignettes/vig.bib
Log:
Get formation energies from Materials API for multi-metal.Rmd
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2020-07-26 07:42:16 UTC (rev 585)
+++ pkg/CHNOSZ/DESCRIPTION 2020-07-26 14:54:53 UTC (rev 586)
@@ -1,6 +1,6 @@
Date: 2020-07-26
Package: CHNOSZ
-Version: 1.3.6-59
+Version: 1.3.6-60
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/man/mosaic.Rd
===================================================================
--- pkg/CHNOSZ/man/mosaic.Rd 2020-07-26 07:42:16 UTC (rev 585)
+++ pkg/CHNOSZ/man/mosaic.Rd 2020-07-26 14:54:53 UTC (rev 586)
@@ -14,7 +14,7 @@
\item{bases}{character, basis species to be changed in the calculation, or list, containing vectors for each group of changing basis species}
\item{bases2}{character, second set of changing basis species}
\item{blend}{logical, use relative abundances of basis species?}
- \item{stable}{list, stable species to use with blend = FALSE}
+ \item{stable}{list, previously determined stable species}
\item{...}{additional arguments to be passed to \code{\link{affinity}}}
}
Modified: pkg/CHNOSZ/vignettes/mklinks.sh
===================================================================
--- pkg/CHNOSZ/vignettes/mklinks.sh 2020-07-26 07:42:16 UTC (rev 585)
+++ pkg/CHNOSZ/vignettes/mklinks.sh 2020-07-26 14:54:53 UTC (rev 586)
@@ -120,3 +120,4 @@
sed -i 's/NaCl()/<a href="..\/html\/NaCl.html">NaCl()<\/a>/g' multi-metal.html
sed -i 's/solubility()/<a href="..\/html\/solubility.html">solubility()<\/a>/g' multi-metal.html
sed -i 's/subcrt()/<a href="..\/html\/subcrt.html">subcrt()<\/a>/g' multi-metal.html
+sed -i 's/convert()/<a href="..\/html\/util.units.html">convert()<\/a>/g' multi-metal.html
Modified: pkg/CHNOSZ/vignettes/multi-metal.Rmd
===================================================================
--- pkg/CHNOSZ/vignettes/multi-metal.Rmd 2020-07-26 07:42:16 UTC (rev 585)
+++ pkg/CHNOSZ/vignettes/multi-metal.Rmd 2020-07-26 14:54:53 UTC (rev 586)
@@ -122,15 +122,15 @@
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 two energies (among others): 1) the formation energy from the elements in their reference state, and 2) the energy above the convex hull defined by the stable materials.
-The energy above the hull is zero for stable materials, and greater than zero for metastable materials.
-Furthermore, the Pourbaix energy (Δ*G*~pbx~) refers to the energy above the hull in addition to potential terms associated with pH and Eh, and is used to predict the stable decomposition products.
-The parallel terminology used in CHNOSZ is that aqueous species or minerals have a 1) Gibbs energy of formation from the elements [Δ*G*° = f(*T*, *P*)], and 2) affinity of formation from the basis species [*A* = -Δ*G* = f(*T*, *P*, and activities of all species)].
+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 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.
+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 leads 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.
@@ -138,11 +138,14 @@
<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
+## 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.728, Fe3O4 = -1.859, Fe2O3 = -1.907)
+Fe.cr <- c(Fe = 0, FeO = -1.72803, Fe3O4 = -1.85868, Fe2O3 = -1.90736)
# 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)
+V.cr <- c(V = 0, V2O3 = -2.52849, V3O5 = -2.52574, VO2 = -2.48496, V3O7 = -2.32836, V2O5 = -2.29524)
# Convert formation energies from eV / atom to eV / molecule
natom.Fe <- sapply(makeup(names(Fe.cr)), sum)
@@ -177,7 +180,7 @@
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)
+# 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")
@@ -194,9 +197,9 @@
iV.cr <- modfun(V.cr, "cr")
iV.aq <- modfun(V.aq, "aq")
-# Formation energies (eV / atom) for bimetallic solids from MP website
+# Formation energies (eV / atom) for bimetallic solids from Materials API
# mp-1335, mp-1079399, mp-866134, mp-1200054, mp-504509 (triclinic FeVO4)
-FeV.cr <- c(FeV = -0.129, FeV3 = -0.171, Fe3V = -0.129, Fe2V4O13 = -2.236, FeVO4 = -1.756)
+FeV.cr <- c(FeV = -0.12928, FeV3 = -0.17128, Fe3V = -0.12854, Fe2V4O13 = -2.23622, FeVO4 = -1.75611)
# Convert energies and add to OBIGT
natom.FeV <- sapply(makeup(names(FeV.cr)), sum)
FeV.cr <- FeV.cr * natom.FeV
@@ -207,7 +210,7 @@
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 assign activities of aqueous species to match 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).
+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.
@@ -255,7 +258,7 @@
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.
We also 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.
-We now have all the ingredients needed to combine the Fe-bearing, V-bearing, and bimetallic species to generate a given composition.
+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].
@@ -312,11 +315,11 @@
Eh <- d11$vals$Eh[imax[2]]
points(pH, Eh, pch = 10, cex = 2, lwd = 2, col = 7)
stable <- d11$names[d11$predominant[imax]]
-text(pH, Eh, stable, adj = c(0.3, 2), cex = 1.2, col = 7)
+text(pH, Eh, stable, adj = c(0.5, -1), cex = 1.2, col = 7)
```
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 previous diagrams 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 (-Δ*G*~pbx~).
+Given the previous diagrams 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.
```{r FeVO4, echo = 22:34, message = FALSE, results = "hide", fig.width = 10, fig.height = 5, out.width = "100%", pngquant = pngquant}
@@ -325,7 +328,7 @@
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:7, message = FALSE, fig.keep = "none"}
+```{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]]
@@ -333,24 +336,30 @@
points(pH, Eh, pch = 10, cex = 2, lwd = 2, col = 7)
stable <- d11$names[d11$predominant[imax]]
text(pH, Eh, stable, adj = c(0.3, 2), cex = 1.2, col = 7)
+range(aFeVO4_vs_stable[d11$predominant == d11$predominant[imax]])
```
-To calculate the energy above the hull, let's define the stable species as the basis species.
-O~2~ is also needed to make a compositionally complete set, but it does not appear in the reaction to form or decompose FeVO~4~.
-`subcrt()` is used to automatically balance the reaction to form FeVO~4~ from the basis species and calculate the standard Gibbs energy of the reaction.
-The value of Δ*G*° in cal/mol (the default for `subcrt()`) is then converted to J/mol, then to eV/mol, and finally eV/atom.
+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.
+The `convert()` command divides this value by *-RT* (which yields `r logK` of the reaction), showing the same result as calculated above from the Pourbaix diagram.
+The value of Δ*G*° in cal/mol (the default for `subcrt()`) is converted to J/mol, then to eV/mol, and finally eV/atom.
+
```{r hull, echo = 1:6, message = FALSE}
b <- basis(c("Fe2O3", "Fe2V4O13", "O2"))
cal_mol <- subcrt("FeVO4", 1, T = 25)$out$G
+convert(cal_mol, "logK")
J_mol <- convert(cal_mol, "J")
eV_mol <- J_mol / 1.602176634e-19
eV_atom <- eV_mol / 6.02214076e23 / 6
round(eV_atom, 3)
-stopifnot(round(eV_atom, 3) == 0.411)
+stopifnot(round(eV_atom, 3) == 0.412)
```
-This is nearly equal to the value of 0.412 eV/atom for the energy above the hull for [triclinic FeVO~4~ on the MP website](https://materialsproject.org/materials/mp-504509/), showing that we successfully made a round trip starting with the input formation energies (eV/atom) from the MP website, to Gibbs energy (J/mol) in the OBIGT database, and back out to energy above the hull (eV/atom).
+This is equal to the value for the energy above the hull for [triclinic FeVO~4~ on the MP website](https://materialsproject.org/materials/mp-504509/), showing 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.
Modified: pkg/CHNOSZ/vignettes/vig.bib
===================================================================
--- pkg/CHNOSZ/vignettes/vig.bib 2020-07-26 07:42:16 UTC (rev 585)
+++ pkg/CHNOSZ/vignettes/vig.bib 2020-07-26 14:54:53 UTC (rev 586)
@@ -941,3 +941,13 @@
volume = {13},
doi = {10.1016/S0010-938X(73)80037-X},
}
+
+ at Article{OCJ_15,
+ author = {Shyue Ping Ong and Shreyas Cholia and Anubhav Jain and Miriam Brafman and Dan Gunter and Gerbrand Ceder and Kristin A. Persson},
+ journal = {Computational Materials Science},
+ title = {The {M}aterials {A}pplication {P}rogramming {I}nterface ({API}): {A} simple, flexible and efficient {API} for materials data based on {RE}presentational {S}tate {T}ransfer ({REST}) principles},
+ year = {2015},
+ pages = {209--215},
+ volume = {97},
+ doi = {10.1016/j.commatsci.2014.10.037},
+}
More information about the CHNOSZ-commits
mailing list