manuscript.org

A Linear Response, DFT+U Study of Trends in the Oxygen Evolution Activity of Transition Metal Rutile Dioxides

1 Introduction

Density functional theory (DFT) is a first principles tool that can be used to understand catalytic processes and identify promising candidates through the calculation of kinetic and thermodynamic properties, which include formation energies, adsorption energies, and reaction barriers cite:noerskov-2009-towar,greeley-2006-comput,norskov-2002-univer-heter-catal,michaelides-2003-ident. Transition metal oxides (TMOs), a class of catalysts used in a wide variety of important chemical processes cite:weaver-2013-surfac-chemis,doyle-2013-redox,gouma-2011-nanos-polym, have thermodynamic and electronic properties that are difficult to capture accurately using standard exchange correlation functionals (LDA and GGA) cite:cohen-2008-insig-curren. The culprit of these inaccuracies is the self-interaction error produced by highly correlated electrons, such as the d-electrons in oxidized systems cite:wang-2006-oxidat,franchini-2007-groun. The Hubbard U (DFT+$U$) is the most feasible correction to account for the self-interaction error cite:anisimov-1997-first,anisimov-1991-band-mott, but its method of application is not trivial. The specific Hubbard U required for a given material can be empirically determined, but the experimental data required oftentimes is not available. For example, adsorption energies on well defined surfaces of oxides are typically difficult to measure cite:campbell-2013-enthal-entrop. Bulk oxidation energies can be used, but the Hubbard U values are typically reaction specific cite:wang-2006-oxidat,aykol-2014-local-gga. In contrast, the Hubbard U can also be calculated via a linear response method cite:cococcioni-2005-linear, but there have been few studies that use this method in for the calculation of catalytic properties.

One of the most studied reactions catalyzed by transition metal oxides is the oxygen evolution reaction cite:dau-2010,parent-2014-progr-base. The oxygen evolution reaction (OER) is the conversion of H₂O into protons, electrons and oxygen. The high energy of protons and electrons can be stored into the chemical bonds of hydrogen, alcohols, or hydrocarbons, while pure oxygen is a widely used oxidant in chemical industries and must be separated from N₂ if acquired from air. The observed trends in kinetics of OER on different catalysts can be related to calculated chemical and electronic properties transition metal oxides cite:rossmeisl-2007-elect,man-2011-univer,suntivich-2011,vojvodic-2011-optim. Key conclusions from these studies are that the adsorption energies of a few intermediates describe the activity trends, these adsorption energies scale with each other, and the scaling of adsorption energies produces an activity volcano with a theoretical activity limit. These conclusions were established without the Hubbard U. While a few studies have applied the Hubbard U to test cases cite:liao-2012-water,garcia-mota-2012-impor, it is still not clear whether the aforementioned conclusions still apply to with the application of the Hubbard U nor if the linear response U will lead to better agreement with experimental results.

In this study, we use DFT+$U$ coupled with the calculated linear response U to evaluate trends in activity of transition metal rutile dioxides for the oxygen evolution reaction. We apply an atomistic thermodynamic method that relates the activity of surfaces to differences in the adsorption energies OH, O, and OOH. We find that the application of any $U$ in almost all cases leads to more endothermic adsorption energies of all intermediates and these shifts in adsorption energies preserve the scaling relationships between OER intermediates calculated with $U=0$. The combination of both observations results in relatively small shifts of all systems to the weak binding side of the OER volcano. We find that these shifts leads to activity trends that are more consistent with experimental observations.

2 Methods

2.1 DFT calculation parameters

All DFT calculations were performed with \textsc{Quantum-ESPRESSO} cite:giannozzi-2009-quant-espres with the Perdew-Burke-Ernzerhof (PBE) exchange correlation functional cite:Perdew1996,Perdew1997a. The core electrons were described by the GBRV library of ultrasoft pseudopotentials cite:garrity-2014-pseud-dft. The kinetic energy cutoff for wavefunctions and the charge density were 40 and 500 Ry, respectively. For surface slabs, we used a 4 × 4 × 1 Monkhorst-Pack grid of k-points cite:monkhorst-1976-special-brill. All calculations were spin-polarized.

The general method for the calculation of the linear response U is described in a previous paper by Cococcioni and de Gironcoli cite:cococcioni-2005-linear. For calculation of linear response U values in the bulk, we applied perturbations up to $±$ 0.15 eV to both the metal and oxygen in 2 × 2 × 2 rutile supercells consisting of 48 atoms to ensure that interactions between the perturbations were and their periodic images were minimal.

2.2 Structural parameters

The equilibrium volume, cell shape, atomic positions of all transition metal dioxide rutile structures were determined by constructing a polynomial equation of state and full relaxation of the shape and atomic coordinates. Ground state magnetic configurations were calculated for all materials, taking into consideration non-magnetic, ferromagnetic, and anti-ferromagnetic orderings.

All adsorption energies were performed on the (110) surface. Because of the large number of calculations we performed in this study, we chose to model the (110) surface as a two layer slab with terminating hydrogen atoms on the bottom layer. A similar two layer slab has been used in previous studies of oxygen evolution on MnO₂ and IrO₂ surfaces cite:steegstra-2013-revis-redox,busch-2012-water-oxidat. The validation of this smaller slab with respect to the typical four layer slab used in similar previous studies cite:rossmeisl-2007-elect,man-2011-univer,halck-2014-beyon is discussed in the results. Figure ref:surface-slabs shows the two layer slab and four layer slab used for validation along with the adsorption site used for all calculations, which is typically called the 5$cus$ site.

2.3 Atomistic thermodynamic framework for oxygen evolution

The atomistic thermodynamic framework we are using to study the oxygen evolution reaction has been used before cite:rossmeisl-2007-elect,man-2011-univer,mom-2014-model-oxygen,calle-vallejo-2013-oxygen, so we only briefly summarize it below. The mechanism of OER is assumed to proceed through four electron proton transfer steps and the OH, O, and OOH intermediates, shown below in acidic conditions.

\begin{align} \ce{H₂O + ∗} &→ \ce{^∗OH + H+ + e-}
\ce{^∗OH} &→ \ce{^∗O + H+ + e-} \ \ce{^∗O + H₂O} &→ \ce{^∗OOH + H+ + e-} \ \ce{^∗OOH} &→ ∗ + \ce{O2 + H+ + e-} \end{align}

At constant pH and with respect to the normal hydrogen electrode (NHE), the Gibbs free energy of each elementary step is shown below,

\begin{align} Δ G_1 &= Δ G_\ce{OH}
Δ G_2 &= Δ G_\ce{O} - Δ G_\ce{OH}\ Δ G_3 &= Δ G_\ce{OOH} - Δ G_\ce{O}\ Δ G_4 &= 4.92[\textrm{eV}] - Δ G_\ce{OOH} \end{align}

\noindent where the adsorption energy of OH, O, and OOH are as follows

\begin{align} Δ G_\ce{O} &= E_slab,\ce{O} - E_slab - (E_\ce{H2O} - E_\ce{H2})
Δ G_\ce{OH} &= E_slab,\ce{OH} - E_slab - (E_\ce{H2O} - \frac{1}{2}E_\ce{H2}) \ Δ G_\ce{OOH} &= E_slab,\ce{OH} - E_slab - (2E_\ce{H2O} - \frac{3}{2}E_\ce{H2}) \end{align}

\noindent where $E_slab,A$ is the total energy of slab with adsorbate $A$, $E_slab$ is the total energy of the bare slab, and $E_\ce{H2O}$ and $E_\ce{H2}$ is the total energy of H₂O and H₂ in an asymmetric box. All adsorbate and gas species included previously reported zero point energy corrections cite:man-2011-univer.

Because each reaction step involves the transfer of an electron to the electrode, applying a potential of U volts on the electrode with respect to NHE would result in a decrease of the $Δ G$ of each reaction step by U eV. When a potential is applied such that the $Δ G$ for all reaction steps is less than zero, all reaction steps are considered exothermic. The potential at which this happens minus 1.23 V is considered the theoretical overpotential, $η^OER$, and is the key metric we use to evaluate the activity of different catalysts for OER. The expression for $η^OER$ is shown below in Equation eqref:complete-eta,

\begin{equation} η^OER = \textrm{Max}[Δ G₁,Δ G₂,Δ G₃,Δ G₄]/e - 1.23V. \label{complete-eta} \end{equation}

The existence of scaling relationships between different reaction energies $Δ G$ gives rise to a descriptor and activity volcano where either $Δ G₂$ or $Δ G₃$ is the largest reaction energy and both of their magnitudes scale with the difference between the adsorption energies of O and OH ($Δ G_O - Δ G_OH$). The details of this analysis can be found in the seminal work that originally established this atomistic thermodynamics cite:rossmeisl-2007-elect,man-2011-univer, but we will be using this relationship to establish a similar volcano plot in our analysis.

3 Results and Discussion

3.1 Validation of the surface slab model

The TMOs we investigated are shown in Table ref:bulk along with their equilibrium lattice constants, magnetic structure, and calculated linear response U. The lattice coordinates and magnetic structure were then used to construct the two and four layer slabs, which are shown in Figure ref:surface-slabs (a), while the linear response U was used when assessing the OER activity trends of the different oxides.

				Magnetic	Linear
Compound	a	c	u	Structure	response U
TiO2	4.65	2.97	0.31	NM	4.95
CrO2	4.38	2.90	0.30	FM	7.15
MnO2	4.36	2.84	0.30	FM	6.63
NbO2	4.94	2.96	0.29	NM	3.32
MoO2	4.95	2.73	0.28	NM	4.83
RuO2	4.53	3.18	0.31	NM	6.73
RhO2	4.55	3.11	0.31	NM	5.97
ReO2	4.95	2.68	0.28	NM	5.27
IrO2	4.54	3.18	0.31	NM	5.91
PtO2	4.59	3.23	0.31	NM	6.25

We first validate the usage of the two layer surface model shown in Figure ref:surface-slabs (a). We motivate the usage of this slab because we are performing over 400 calculations using U values and would like to minimize the computational cost. We first calculate the adsorption energies of OH, O, and OOH on both the two layer and four layer slab at $U=0$ for all systems. Figure ref:surface-slabs (b) shows a parity plot between adsorption energies calculated on both slabs, and we see excellent agreement for OOH, good agreement for OH, and reasonable agreement for O. More importantly, Figure ref:surface-slabs (c) also shows that both sets of adsorption energies fall on the same scaling relationship. This suggests that a majority of the differences between the two adsorption energies are systematic, and that the underlying physics that results in the scaling relationships is the same for both the two and four layer slab. Following these results, we moved on to calculate adsorption energies on the two layer slab at $U>0$.

3.2 Variation of adsorption energies and scaling relationships with respect to U

For all materials, we calculated the adsorption energies by applying a $U=0$ eV to $U=8$ eV in 0.5 eV intervals. The starting geometry was taken from the relaxed structure of the calculation at $U=0$. The entirety of our results can be found in the supporting information cite:xu-suppor, but for brevity we discuss results for only NbO₂, IrO₂, TiO₂, and MnO₂ below. Observations for NbO₂ and IrO₂ were characteristic of early and late 4$d$ and 5$d$ transition metal dioxides, respectively. For 3$d$ systems, TiO₂ is a special case, and observations for MnO₂ and CrO₂ were similar.

We found that the application of $U>0$ had a number of systematic effects to the adsorption energies of OH, O, and OOH to the 4$d$ and 5$d$ TMO rutiles. These are summarized in Figure ref:4d-5d-adsorption. First, the application of U results in shifts to more endothermic adsorption energies of all species on all compounds (Figure ref:4d-5d-adsorption (a) and (c)). For low U values, these shifts are monotonic and smooth, but for high U values on early TMOs of MoO₂, NbO₂, and ReO₂, they deviate from the monotonic trend at low U values. It is likely that such high values of U are not appropriate for these early TMOs. Early TMOs have a smaller occupancy of $d$-electrons, and therefore one would expect a lower value of U is needed to correct the self-interaction error. This is supported by the lower linear response U for the early 4$d$ and 5$d$ TMOs (Table ref:bulk). We also observed that the calculated linear response U for all early 4$d$ and 5$d$ TMOs sits right at the point where the smooth, monotonic $Δ E_ads(U)$ behavior breaks down. This is shown for NbO₂ in Figure ref:4d-5d-adsorption (a) and the rest of the early TMOs in the supporting information. This is further evidence that high U values are not appropriate for early TMOs. For late TMOs of PtO₂, IrO₂, RuO₂, and RhO₂, the changes are smooth all the way up to a $U=8$, including their calculated linear response U values.

For 4$d$ and 5$d$ oxides, the U-induced endothermic changes of the adsorption energy preserve scaling relationships established at $U=0$. This is shown in Figure ref:4d-5d-adsorption (b) and (d). This is true for all U values tested on the 4$d$ and 5$d$ TMOs, including high U values on early TMOs. Our results further demonstrates the robustness of scaling relationships, showing that the additional physics via the Hubbard U does not lead to deviations of scaling relationships. This also demonstrates that correlations between the electronic structure and adsorption energies implied by the scaling relationships are also preserved with the addition of U. This conclusion is consistent with previous work that found similar electronic structure/activity correlations on doped TiO₂ with both DFT and DFT+$U$ results cite:garcia-mota-2012-elect-tio2.

In contrast to our results on 4$d$ and 5$d$ TMOs, we found a mixture of results for 3d TMOs. Adsorption energies at $U>0$ on CrO₂ and MnO₂ gave similar results to each other, with adsorption on MnO₂ shown in Figure ref:3d-adsorption (a) and (b). With increasing U values, we observe a smooth monotonic increase in the adsorption energy, but at some intermediate U value adsorption of OOH on the surface is no longer stable for some species, shown by the lack of change in adsorption energy for $U>4$ eV for MnO₂. This is what gives rise to the breaking of the scaling relationship between OH and OOH at U ≈ 4 eV. We are unsure how to interpret the breaking of the surface-adsorbate bonds, but it is clear that even at low U values when adsorption was stable, the scaling relationships are preserved. This is consistent with our results on 4$d$ and 5$d$ TMO rutiles.

For TiO₂, application of U produces smooth, monotonic changes in the adsorption energy (Figure ref:3d-adsorption (c)), but interestingly the change in the OOH adsorption energy is exothermic upon increasing U. This was the only adsorption energy where the addition of U produced a more exothermic adsorption energy. Also unique to TiO₂ is that the scaling relationships are not preserved with the addition of U (Figure ref:3d-adsorption (d)). The relative change in the adsorption energy with respect to increasing U is also small. $Δ E_ads^OH$ changes by less than 0.1 eV by applying a U value of 8 eV.

There is still conflicting literature on how the Hubbard U should be implemented to capture accurate thermodynamic properties of Ti oxide systems cite:jain-2011,hu-2011-choic-u,yan-2013-calcul,aykol-2014-local-gga. Our results show this is still an open issue for adsorption on TiO₂. The Ti ion at the adsorption site of a stoichiometric TiO₂ has a $d⁰$ configuration and OH, O, or OOH primarily forms bonds with the 3$p$ electrons. Hence, adsorption induced changes to the electronic structure of the Ti $d$ electrons are subtle, which is reflected by the smaller change in the adsorption energies induced by adding a Hubbard U. This electronic structure phenomenon is typical for stoichiometric surfaces of closed shell materials, such as adsorption on stoichiometric alkaline-earth metal oxides cite:bajdich-2014-surfac-energ, and leads to deviations in trends of both adsorption and oxygen vacancy formation energies with respect to number of electrons cite:akhade-2012-effec-strain,calle-vallejo-2013-number,curnan-2014-effec-concen. This situation is not encountered in any other of the adsorption eneriges we studied.

Because of this unique change in the electronic structure caused by adsorption on TiO₂, we hypothesize that the application of the Hubbard U to the $d$-electrons of the TiO₂ after adsorption may require different treatment. The U we calculated for bulk TiO₂ likely does not describe the TiO₂ with a 3$p$ hole state. To resolve this special case, one might be required to calculate separate Hubbard U values of the Ti ion with and without an adsorbate and use the DFT+\mathit{U}(R) method to calculate an adsorption energy that takes changes in U into account cite:kulik-2011-accur-poten. Another possibility is the requirement of application of U to lattice oxygen 2$p$ states or Ti 3$p$ states. A relatively high U of 6 eV applied to the oxygen 2$p$ states was required to accurately capture hole states in SiO₂ doped with Al cite:nolan-2006-hole-al.

To summarize, we draw two main conclusions from our analysis of adsorption energies and scaling relationships with respect to increasing U values. With the exception of TiO₂, where the significance of the Hubbard U to calculate adsorption on TiO₂ remains unclear, the application of U produces more endothermic adsorption energies, and these changes in adsorption energy preserve the scaling relationships established at $U=0$. These conclusions further validate the scaling relationships and their usage for establishing models for catalytic reactions on TMOs cite:medford-2014-asses. The similar weakening of adsorption energies with respect to U also suggests that a majority of trend studies of adsorption on TMOs at $U=0$ are probably valid at $U>0$. Our results also provide researchers with useful estimates on the effect of U on adsorption energies. Having established some general rules between the Hubbard U, adsorption energies, and scaling relationships, we now move towards the specific application of OER and the usage of the linear response calculated U.

3.3 Activity trends with linear response U value

We next evaluate the effect of applying a calculated linear response Hubbard U to the activity trends for OER. We focus our analysis on the IrO₂, PtO₂, RuO₂, and RhO₂ oxides in our study. We choose only these materials for a number of reasons. First, from Pourbaix diagrams, one can easily see that CrO₂, MoO₂, NbO₂ and ReO₂ are not stable in either acidic or alkaline OER conditions cite:pourbaix-1974-atlas. In contrast, IrO₂, PtO₂, RuO₂, and RhO₂ are predicted to be stable at acidic OER conditions and in some cases have been observed in situ in experimental work cite:pourbaix-1974-atlas,sanchez-2014-in-situ,silva-2000-in-ruo2. MnO₂ was not used in this comparison for two reasons. First, it is still unclear whether MnO₂ is the active species at OER conditions. Recent studies have identified that the Mn³⁺ as the active species in OER cite:ramirez-2014-evaluat-mnox,gorlin-2010. Second, our results point towards OOH desorption at the linear response, calculated U value. TiO₂ was not used in our comparison due to our conclusion that our DFT+\textit{U} method did not seem appropriate for an accurate calculation of adsorption energies and it is not a good OER catalyst.

Figure ref:OER-volcano shows the changes in the activity of the selected oxides as one applies the linear response, calculated Hubbard U. As expected from the observed preservation of scaling relationships, the changes in the adsorption energy produced by applying the linear response U for all species results in movement along the weak binding and strong binding legs of the volcano, but not changes in the activity volcano itself. Furthermore, all species are moved towards the weaker binding leg of the volcano, which is explained by the universal weakening of adsorption energies caused by applying the Hubbard U.

The combination of these two observations leads to changes in the relative ordering of activity. With DFT, we predict the activity trend to be RhO₂ $>$ IrO₂ $>$ PtO₂ $>$ RuO₂. With the addition of the calculated Hubbard U, we now predict the activity to be IrO₂ $>$ RhO₂ $>$ RuO₂ $>$ PtO₂. The ordering with the addition of the Hubbard U shows better agreement with experiments, which has been observed as RuO₂ ≈ IrO₂ $>$ RhO₂ $>$ PtO₂ cite:miles-1976-period,cherevko-2014-dissol-noble. We still obtain discrepancy with regard to the activity of RuO₂ with respect to RhO₂, but the addition of U improves agreement with experimentally observed trends. IrO₂ and RuO₂ move towards the top of the volcano from the strong binding side, while RhO₂ and PtO₂ move away from the top of the volcano on the weak binding side. The combination of these two effects corrects the incorrect ordering of RhO₂ $>$ IrO₂ and PtO₂ $>$ RuO₂. We note that previous results observed a different ordering between these compounds, found to be RuO₂ $>$ PtO₂ ≈ RhO₂ $>$ IrO₂, at $U=0$. We associate these slight differences with differences in pseudopotentials, calculation parameters, and the implementation of different surface models. However, both set of results saw IrO₂ and RuO₂ on the strong binding (left) side of the volcano and RhO₂ and PtO₂ on the weak binding (right) side of the volcano. Hence, it is likely the application of the linear response U to those results should give similar improvements to those seen here, with IrO₂ predicted to be more active than PtO₂ and RhO₂.

We also comment that though changes in ordering are observed, the absolute changes in reaction energies are relatively small. The changes in reaction energy with the application of the calculated U value was on the order of 0.2 $∼$ 0.4 eV, which in no case was enough to move a species from the strong binding to weak binding side of the volcano. Hence, we propose that large scale screening studies based on correlations between adsorption energies done without and with the Hubbard U should produce similar conclusions, except perhaps near the top of the volcano.

4 Conclusions

To summarize, we have performed a DFT+\textit{U} study on the adsorption of OER intermediates on the (110) surface of rutile transition metal dioxides. Our analysis focused on changes in the adsorption energy, scaling relationships, and activity trends by applying a range of Hubbard U values in addition to the linear response, calculated U value. We find that with the exception of TiO₂, the application of a large range of Hubbard U values produces more endothermic adsorption energies and preserves scaling relationships established at $U=0$. We also find that when linear response U values applied, the relative ordering of the activity of IrO₂, PtO₂, RuO₂, and RhO₂ oxides improves with respect to experimental observations. Our work reveals a number of universal relationships between the Hubbard U and catalytic processes on transition metal oxides.

\begin{acknowledgement} This work was partially supported by the IMI Program of the National Science Foundation under Award No. DMR 08-43934. We gratefully acknowledge support from the DOE Office of Science Early Career Research program (DE-SC0004031). \end{acknowledgement}

\begin{suppinfo} Full details of the computational setup and analysis of the computations is available. All input and output files can be found online at http://dx.doi.org/10.5281/zenodo.12635, and the supporting information contains scripts used for generating and analyzing these files. \end{suppinfo}

bibliography:references.bib

\begin{tocentry}

\includegraphics{./figures/TOC.png}

\end{tocentry}