Predicting aqueous stability of solid with computed Pourbaix diagram using SCAN functional

In this work, using the SCAN functional, we develop a simple method on top of the Materials Project (MP) Pourbaix diagram framework to accurately predict the aqueous stability of solids. We extensively evaluate the SCAN functional’s performance in computed formation enthalpies for a broad range of oxides and develop Hubbard U corrections for transition-metal oxides where the standard SCAN functional exhibits large deviations. The performance of the calculated Pourbaix diagram using the SCAN functional is validated with comparison to the experimental and the MP PBE Pourbaix diagrams for representative examples. Benchmarks indicate the SCAN Pourbaix diagram systematically outperforms the MP PBE in aqueous stability prediction. We further show applications of this method in accurately predicting the dissolution potentials of the state-of-the-art catalysts for oxygen evolution reaction in acidic media.


INTRODUCTION
Predicting aqueous stability of solids is of central importance to materials development for electrochemical applications such as fuel cells, electrolyzers, and Li-air batteries [1][2][3][4] . The electrode potential-pH (aka. Pourbaix) diagram maps the most stable species in an aqueous environment and has become an invaluable tool in evaluating aqueous stability. Hansen et al. 1 developed the concept of surface Pourbaix diagram to identify the most stable surface of metal cathodes under the operating conditions of oxygen reduction reaction. Persson et al. 2 proposed an efficient scheme unifying the Gibbs free energies of ab initio calculated solids and the experimental aqueous ions for the Pourbaix diagram construction. The advent of large computational material databases (e.g., the Materials Project (MP)) has enabled highthroughput screening of materials of interest in terms of their aqueous stability [3][4][5][6] . Nevertheless, aqueous stability screening based on the current implementation of the MP Pourbaix diagram 2,5 leads to some unexpected inaccuracies. One example is selenium dioxide (SeO 2 ), which is incorrectly predicted to be stable in acid and water 7 . This qualitative inconsistency occurs in that the corrected formation energy of SeO 2 is an outlier in the MP data, exhibiting a large error (1.18 eV per formula unit). In the MP, a correction of 1.40 eV/O 2 is used for oxide compounds to calibrate the error in calculated formation energy, arising from the overestimated O 2 binding energy in the Perdew-Burke-Ernzerhof (PBE) functional 8 . Though this approach significantly improves the accuracy of computed formation energies and phase diagrams overall, large deviations can persist for many compounds like selenium (See Supplementary Fig. 1), which may affect predicted positions of equilibrium lines on the pourbaix diagram.
Since the spread of this error, even when systematically correcting formation energies, may have an impact on the positioning of aqueous phase equilibria, improving the overall accuracy of the physical methods to calculate formation energy of solids can significantly improve the reliability of aqueous stability prediction. Recently, the strongly constrained and appropriately normed (SCAN) functional has demonstrated such an improvement accuracy in both energetics and geometries 9-12 . Specifically, Zhang et al. 12 demonstrated that the formation energies of 196 binary compounds using the SCAN functional reduce the mean absolute error (MAE) to 0.1 eV/atom, halving that of the PBE functional. This remarkable improvement is attributed to its satisfaction of all known exact constraints applicable to a semilocal functional 9,10 . The SCAN functional still has difficulties in accurately predicting the formation energy for transition-metal compounds due to the self-interaction error 12 .
In this work, using the SCAN functional, we developed an approach that modifies the MP Pourbaix diagram construction to predict the aqueous stability of solids. The SCAN functional with Hubbard U corrections is used to improve performance in calculated formation energy for transition-metal oxides (TMOs). When compared to the MP PBE Pourbaix diagram, SCANcomputed Pourbaix diagrams systematically improve the accuracy of aqueous stability predictions, and therefore enable an avenue to more reliably screen materials by this figure of merit in highthroughput computations.

RESULTS
Formation enthalpy of oxides Though Zhang et al. 12 investigated the formation enthalpy of 196 binary compounds using the SCAN functional, only 46 compounds are oxides, and most importantly, some oxides of technologically relevance (e.g., IrO 2 and RuO 2 , both are popular oxygen evolution reaction (OER) catalysts) are not included in their benchmarking. In this work, we benchmark against a more comprehensive dataset of 114 binary oxides to demonstrate the performance of the SCAN functional in formation enthalpy calculations, which are particularly important in considerations of aqueous stability. However, the SCAN functional still performs poorly in calculated formation enthalpies for 3d TMOs (MAE = 0.205 eV/atom) due to self-interaction errors 12,13 . Using the MP methods, outlined in refs. 8,14 , we developed a set of Hubbard U corrections to aid the SCAN functional to accurately calculate formation enthalpies for TMOs (TM = V, Cr, Mn, Fe, Co, Ni, Mo) (Supplementary Table 1). These TMOs for Hubbard U corrections were selected by following the MP convention and TMOs studied in this work refer to these seven classes of TM oxides specifically 8,14 . We note that anion corrections from the standard suite of MP corrections 8 are omitted in our approach because the SCAN functional has a very small error in the calculated O 2 binding energy (5.27 eV vs. experimental 5.12-5.23 eV 15,16 ). All PBE functional data used in this paper were retrieved from the MP 5,17,18 . Figure 1 presents calculated formation enthalpies (ΔH calc ) of 114 binary oxides using the PBE(+U), SCAN, and SCAN(+U) functionals. The MAE of ΔH calc for PBE(+U), SCAN, and SCAN(+U) are 0.182, 0.094, and 0.072 eV/atom, respectively. The inset figure shows ΔH calc of TMOs using PBE+U, SCAN, and SCAN+U functionals. The non-U-corrected SCAN functional yields an MAE of 0.156 eV/atom in ΔH calc and larger than 0.084 eV/atom of PBE +U. Whereas, SCAN+U significantly reduces the MAE of ΔH calc from 0.156 to 0.025 eV/atom and is also smaller than PBE+U. In addition, SCAN+U using the U fitted to binary TMOs performs well on ternary systems, achieving an MAE of 0.053 eV/ atom on a benchmark of 47 ternary TMOs ( Supplementary Fig. 2).

Validation
We first benchmark our method by validating the calculated Pourbaix diagram of three examples of Ti, Ta, and Se with respect to their experimental Pourbaix diagram. Since the construction of Pourbaix diagrams depends on the species included, we here use consistent species for direct comparison between experiment and computation. The primary difference between computed and experimental Pourbaix diagrams is that the chemical potentials of solids in the computed Pourbaix diagram are calculated using the SCAN(+U) functional instead of experimental data. For theoretical Pourbaix diagrams, we use experimental thermodynamic data of aqueous ions according to the integration scheme of Persson et al. 2 . To further validate the accuracy of our method in aqueous stability predictions, we also generate a SCAN calculated Pourbaix diagram by computing and including all materials in the given composition space from the MP (termed "complete Pourbaix diagram") to compare with the MP PBE results. This complete Pourbaix diagram serves as a valuable tool to predict the stability of a material with untabulated thermodynamic data. The experimental chemical potentials of solids and aqueous ions are from refs. [19][20][21] . Figure 2a shows the calculated Pourbaix diagram for Ti using the SCAN functional. Compared to the experimental diagram ( Fig.  2b), the solid and aqueous stability regions are extraordinarily well-represented. The Ti Pourbaix diagram presented in the literature 2,22 is generally constructed without TiH 2 , as shown in Supplementary Fig. 3, where an excellent agreement between our calculated and experimental results is achieved as well. The inclusion of TiH 2 will considerably reduce the stability region of  In the MP PBE Pourbaix diagram (Fig. 2d), the predicted stability region of Ti 2+ is much larger than that in the experimental diagram and TiO 2 dissolves to TiO 2þ 2 instead of being oxidized to Ti 2 O 7 at high potentials (E > 2.0 V) in acid. Figure 2e shows the calculated Ta Pourbaix diagram using the SCAN functional, which is consistent with the experimental results in Fig. 2f. Tantalum pentoxide, Ta 2 O 5 , is barely reactive in aqueous solution and often serves as a catalyst or catalyst stabilizer for oxygen evolution in both acid and base because of its high aqueous stability 25,26 . We note that when adding the aqueous ion TaO 2þ 2 in the Pourbaix diagram generation procedure (Supplementary Fig. 4), Ta 2 O 5 dissolves to TaO 2þ 2 in acid, contradicting experimental observations 25,26 . This is because Ta 2 O 5 is corroded in concentrated hydrofluoric or alkaline media 27 . Some care has to be taken when screening aqueous stable materials for electrocatalysis, where dilute concentrations are commonly used in experiment. In the complete Pourbaix diagram (Fig. 2g), highvalence tantalum oxides (e.g., Ta 2 O 7 ) and tantalum hydrides (e.g., TaH 2 ) are predicted to be stable under highly oxidative and reductive environments, respectively. Similar results are also observed in the MP PBE Pourbaix diagram (Fig. 2h).
In the above two examples, we have verified that our method makes predictions which have similar high qualitative accuracy as MP PBE 2,5 . Next, we present an example where our results are significantly different from MP PBE. Figure 2i shows the calculated selenium Pourbaix diagram using the SCAN functional, which closely reflects the experimental diagram, shown in Fig. 2j. The complete Pourbaix diagram using the SCAN functional is also the same as the experimental diagram (Fig. 2k). There is only one solid region (Se metal) in the SCAN calculated and experiment diagrams. In the MP PBE scheme, selenium oxides (SeO 2 , Se 2 O 5 , SeO 3 , SeO 4 ) are overstabilized in acid and appear prominently in the Pourbaix diagram, shown in Fig. 2l. This overstabilization results from the MP anion correction scheme, which overcorrects the calculated formation enthalpies of selenium oxides. For example, the MP anion correction shifts the PBE calculated formation enthalpy of SeO 2 from −2.111 to −3.515 eV/fu (ΔH exp = −2.334 eV/fu). 3d metal (Mn/Fe/Co/Ni) Pourbaix diagram By applying U corrections on top of O 2 correction 8,14 , the MP PBE scheme provided the state-of-the-art computational 3d metals (M = Mn/Fe/Co/Ni) Pourbaix diagrams (Fig. 3). Nevertheless, the scheme resulted in qualitative and quantitative inaccuracies when predicting the aqueous stability for metal (oxy)hydroxides and hydrides. This is because the affine corrections of the MP scheme may over-or undercorrect certain classes of materials or reference hydrogen chemical potentials in order to fix systematic errors in oxides. SCAN's accuracy in predicting the underlying thermochemistry makes affine corrections largely unnecessary, meaning that more accurate Pourbaix diagrams can likely be constructed using the simpler methodology herein going forward. Figure 3 presents the calculated and experimental 3d metal (M = Mn/Fe/Co/Ni) Pourbaix diagrams. In general, the predicted solid and aqueous ion regions in the SCAN calculated diagrams are consistent with the experimental results. A notable difference is that the relative regions of stable species in the SCAN calculated diagrams are slightly different from those in experimental diagrams. This difference can be attributed to the uncertainties of the SCAN+U functional in calculating redox reactions energies of M 2+ → M 3+ → M 4+ . Overall, the calculated Pourbaix diagrams using the SCAN+U functional show better agreement with the experimental diagrams than MP PBE results (Fig. 3). A detailed discussion of each 3d metal Pourbaix diagram is presented.
Mn. In Fig. 3a, the stability region of Mn(OH) 2 is slightly broader than that in the experimental diagram (Fig. 3b), resulting in the absence of Mn(OH) + . In the complete Pourbaix diagram (Fig. 3c), MnOOH appears in the stability regions of Mn 2 O 3 and MnO 2 , which is in line with the experimental observations 28 that MnOOH is synthesized at low temperature (<70°C) around pH = 11 and transforms to Mn 2 O 3 and MnO 2 at high temperature.
Fe. In Fig. 3e, the solid and aqueous ions regions are well predicted by our method, compared to the experimental diagram, shown in Fig. 3f. In the complete Pourbaix diagram (Fig. 3g), FeOOH shows up in the Fe 2 O 3 stability region. In experiment, FeOOH is found at low temperature and transforms to Fe 2 O 3 after thermal treatment 29,30 . Moreover, high-valence iron oxides (FeO 3 , Fe 4 O 13 ) and iron hydrides are found to be stable under extremely oxidative and reductive environments, respectively.
Co. In Fig. 3i, the stability regions of solids and aqueous ions are in good agreement with the experimental diagram (Fig. 3j). Similar to iron results, in the complete Pourbaix diagram (Fig. 3k), highvalence cobalt oxides and cobalt hydrides are found to be stable under extremely oxidative and reductive environments, respectively. The MP PBE results (Fig. 3l) shows large deviations from the experimental diagram, which is likely because of the inaccurate chemical potential of HCoO À 2 used (See Supplementary Table 2).
Ni. Compared to the experimental Pourbaix diagram (Fig. 3m), the relatively large stability region of Ni(OH) 2 results in the absence of Ni 3 O 4 in the calculated Pourbaix diagram (Fig. 3n). In the complete Pourbaix diagram (Fig. 3o), NiOOH appears in the regions of Ni 2 O 3 and Ni 3 O 4 , which is in line with the experimental observations that nickel oxyhydroxides often form instead of nickel oxides in solution 31,32 . The stability regions of solids and aqueous ions predicted by the MP PBE are in excellent agreement with the experiment diagram, except for the absence of nickel (oxy)hydroxides (Fig. 3p).

Applications
In this section, we showcase applications of our method in accurately predicting the aqueous stability of catalysts for electrocatalysis. IrO 2 and RuO 2 are the state-of-the-art electrocatalysts for OER in acidic media. The typical OER working conditions in acid are pH = 0 or 1 and potentials of 1.23-2.0 V. Singh et al. 33 proposed that the aqueous stability can be quantitatively evaluated by computing the material's Gibbs free energy difference (ΔG pbx ) with respect to the stable domains in the Pourbaix diagram as a function of pH and potential. We will use this concept to evaluate the aqueous stability of IrO 2 and RuO 2 at OER working conditions. Figure 4a, b shows the calculated Ir Pourbaix diagram and the Pourbaix decomposition free energy of IrO 2 from the potential 0 to 2.0 V at pH = 0, respectively. We find that at pH = 0, IrO 2 transforms to IrO 3 at the potential of 1.55 V and IrO 3 dissolves to IrO À 4 at the potential of 1.60 V, which is consistent with the experimental observations that IrO 2 dissolves at potentials of 1.50-1.60 V in OER [34][35][36] . In experiment, the formation of highvalence iridium oxides, e.g., IrO 3 is commonly regarded as the origin of instability of IrO 2 during OER 35,37 . Based on the density functional theory (DFT) calculated phase diagram ( Supplementary  Fig. 5), IrO 3 is predicted to be highly unstable in the SCAN functional with a decomposition energy of 159 meV/atom, whereas it is predicted be stable in the PBE functional. Considering that IrO 3 is only detected as an intermediate species during OER in experiment 35,37 , we argue that the SCAN functional may give a reasonable phase stability for IrO 3 . As shown in Supplementary   Figure 4a shows the calculated Ru Pourbaix diagram. There are only two solid regions, Ru and RuO 2 . In acid, the ruthenium metal is severely corroded, forming Ru 3+ , and RuO 2 has also a narrow stable region. In experiment, it is generally believed that RuO 2 exists in the form of hydrated oxide, RuO 2 ⋅mH 2 O (1 ≤ m ≤ 2) and could be oxidized or reduced into amorphous Ru 2 O 5 and Ru (OH) 3 ⋅H 2 O, respectively 21 . Figure 4b shows the calculated Pourbaix decomposition free energy of RuO 2 as a function of potentials at pH = 1. We observe that RuO 2 dissolves to H 2 RuO 5 at 1.34 V, which is in line with the experimental results that RuO 2 has a dissolution onset potential of 1.37 V in 0.1 M H 2 SO 4 34 . In the MP PBE Pourbaix diagram ( Supplementary Fig. 7), RuO 4 is predicted to be stable from 0.92 to 2.0 V. Experiments indicate that RuO 4 is unstable and readily dissolves in aqueous solution by forming H 2 RuO 5 21 .

DISCUSSION
Distinct from the MP method 2 , the chemical potential of liquid water in this work is evaluated by using the gas-phase H 2 O at 0.035 bar, since at this pressure the gas-phase H 2 O is in equilibrium with liquid water at 298.15 K 38 . Using SCAN, this results in a very accurate ΔG f = −2.461 eV/fu without any affine corrections. In the MP framework of Pourbaix diagram construction, the chemical potential of hydrogen is corrected such that experimental Gibbs free formation energy of water is exactly replicated by construction. We use a chemical potential of −3.66 eV/atom for hydrogen, equal to half of the Gibbs free energy of hydrogen gas computed using the SCAN functional. We therefore regard the errors in the Gibbs free energy of formation for (oxy)hydroxides and hydrides as those of the SCAN functional itself, instead of from applied hydride corrections (e.g., 0.7-0.8 eV/ H 2 ) in the standard PBE functional 2,39 .
In the reaction energy calculation, we neglect the zero-point energy (E ZPE ) and the integrated heat capacity from 0 to 298.15 K (δH) for both solids and gases by assuming that the differences in these quantities between reactants and products are negligible at room temperature. Previous studies in the literature show that E ZPE and δH of solids and gases are comparable [39][40][41] , but we note that estimation of E ZPE and δH is possible, albeit expensive, with phonon or molecular dynamics calculations, and a systematic approach towards these that can further improve our approach merits a study in its own right.
In conclusion, by building on the MP Pourbaix diagram framework, we have developed a simple method to accurately predict the aqueous stability of solids. Benchmark results show that the calculated Pourbaix diagram using the SCAN functional improves the accuracy of predicted aqueous stability in a large number of cases of interest compared to the MP PBE results. We have also demonstrated the applications of this method in quantitatively predicting the dissolution potential of the state-ofthe-art OER catalysts in acid, for which both the experimental and/ or the MP PBE predictions show large deviations. Our approach enables reliably searching materials in terms of their aqueous stability for electrochemical applications.

Density functional theory
The spin-polarized DFT calculations were performed using the Vienna ab initio simulation package (VASP) within the projector-augmented wave method 42,43 . The SCAN functional 9 was used for all structural relaxations and energy calculations. A plane wave energy cutoff of 520 eV was used, and the electronic energy and atomic forces were converged to within 10 −5 eV and 0.02 eV/Å, respectively. The Brillouin zone was integrated with a k-point density of at least 1000 per reciprocal atom 5 . We followed the MP selection of TMOs (M = V, Cr, Mn, Fe, Co, Ni, Mo) for Hubbard U corrections 8,14,44 . Note that the U correction for tungsten (W) oxides was not considered in this work due to the poor convergence of SCAN calculations. Different magnetic configurations were only considered for Mn, Fe, Co, and Ni TMOs in the Pourbaix diagram construction. All crystal structure manipulations and data analysis were carried out using the Pymatgen software package 17 .

Pourbaix diagram
All the Pourbaix diagrams were constructed using the MP method developed by Persson 2 as implemented in Pymatgen 17 . In the MP method, the stable domains on the Pourbaix diagram are determined based on the knowledge of all possible equilibrium redox reactions in the chemical composition of interest. In an aqueous medium under a given pH (−log [H + ]) and potential (E), the following redox reaction occurs: At equilibrium, the Gibbs free energy change (ΔG rxn ) of this reaction can be related to E using the Nernst equation where ΔG o rxn is the Gibbs free energy change of the reaction at standard conditions, F is the Faraday constant, R is the gas constant, and T is the temperature. a is the activity coefficient. The most stable species in aqueous solutions can be therefore determined by minimizing (ΔG rxn + nFE) across all possible reactions under certain pH and applied potential.
The chemical potential of a material under standard conditions can be expressed as where E DFT is the DFT calculated total energy, E ZPE is the zero-point energy, δH is the integrated heat capacity from 0 to 298.15 K, and TS 0 is the entropy contribution at the standard conditions. For gas-phase species, the DFT total energy was evaluated by a molecule in a 15 × 15 × 15 Å box; S 0 data were collected from the experimental thermodynamic database 45 . For solid phase species, the entropy contribution was neglected since it was very small compared to the gas species. The E ZPE and δH for both solid and gas species were neglected by assuming that these energies have similar contributions between species and therefore cancel out in reaction energies. The chemical potentials of aqueous ionic species were obtained from experimental data with a correction of the referenced solid phase energy difference between DFT calculations and experiments (μ solid-solid ) μ 0 ion ¼ μ 0 ion;exp þ Δμ solidÀsolid : The experimental formation Gibbs free energy of referenced solid and aqueous ions was collected from refs. [19][20][21] .

Energy correction
In the MP PBE Pourbaix diagram construction, there are four primary corrections (anion correction, U correction for TMOs, hydroxide/peroxide correction, and water correction) applied to calibrate the calculated chemical potential of solids. Specifically, an energy correction of 1.40 eV/O 2 was used for oxide compounds (See Supplementary Fig. 1). The Hubbard U correction was adopted for selected TMOs 5,8 . The hydroxides correction was used to calibrate the errors in calculated formation energy of hydroxides. The water correction was to adjust the total energy of water in order to reproduce the chemical potential of water 2 . In the SCAN Pourbaix diagram construction, only the Hubbard U correction was used for selected TMOs.

DATA AVAILABILITY
All data necessary to support the findings of this study are available in the Supplementary Information. Further data and codes can be made available from Z.W.
Z. Wang et al.