Modeling the Oxygen Evolution Reaction on Metal Oxides: The

Article pubs.acs.org/JPCC

Modeling the Oxygen Evolution Reaction on Metal Oxides: The Infuence of Unrestricted DFT Calculations Rik V. Mom,† Jun Cheng,‡ Marc T. M. Koper,† and Michiel Sprik*,‡ †

Leiden Institute of Chemistry, Leiden University, P.O. Box 9502, 2300RA Leiden, The Netherlands Department of Chemistry, University of Cambridge, Cambridge CB2 1EW, United Kingdom

‡

S Supporting Information *

ABSTRACT: Due to the contraints imposed on the optimization of oxygen evolution reaction (OER) catalysts, the correlations between adsorption energies of OER intermediates have received considerable interest. Computation has made important contributions uncovering and elucidating these correlations. The calculations have predominantly been of the spin-restricted type. Using both the restricted and unrestricted formalism, we have performed a DFT−PBE study of the key OER intermediates on 10 nonmagnetic metal oxide surfaces. We show that, depending on the binding energy, unrestricted DFT calculations may yield considerably higher binding energies than restricted DFT calculations. At low binding energies, spin quenching of the adsorbates is inefficient, leading to a strong open shell character that is not well described by restricted DFT. At high binding energies, spin quenching is complete, so that the unrestricted DFT binding energy equals the restricted DFT binding energy. We conclude therefore that although strong binding systems can be properly described by restricted DFT, unrestricted DFT is a necessity for weak binding systems, to which photocatalysts typically belong. The use of unrestricted DFT has little effect on the shape of volcano plots when ΔGO − ΔGOH is chosen as a descriptor. The location of the data points on the plot does show a considerable shift for weak binding substrates.

■

INTRODUCTION Electrochemical and photochemical water splittings are promising methods for the creation of sustainable hydrogen fuel.1,2 Both processes are hampered by inefficiency, however, adding to their high costs. A major bottleneck is the overpotential associated with the oxygen evolution reaction (OER) at the (photo)anode. Modeling the OER provides atomistic information that can lead to design principles for catalytic anodes. For the oxygen reduction reaction on metals, this has proven very successful.3,4 Recently, focus has shifted to extend this approach to the OER on metal oxides.5−8,8−13 Although the mechanism of the OER on metal oxides can vary for different types of surfaces, the single site mechanism shown in Figure 1 is often stated as the predominant pathway.7,9 Here, we will only focus on the hydroxyl, oxyl, and peroxyl intermediates, as the second and third protoncoupled electron transfers (PCETs) are thought to be the most difficult steps.9 Henceforth, these intermediates shall be referred to as *OH, *O, and *OOH. The catalytic effect that the (photo)anode may have on the OER reaction arises from the tuning of the relative stabilities of the intermediates. Ideally, the free energy difference of each PCET step is 1.23 eV versus the standard hydrogen electrode, allowing for the OER to occur at its equilibrium potential, assuming barriers surmountable at room temperature. However, studies by several groups indicate that the ideal case cannot be achieved, because the binding energies of the intermediates are correlated.7,12,14 These scaling relationships © 2014 American Chemical Society

Figure 1. Oxygen evolution reaction mechanism. The horizontal lines denote the surface. The * sign marks an empty adsorption site.

impose restrictions on the possibilities to optimize the relative binding energies of the OER intermediates. Most notably, it has been reported for both metals and metal oxides that there is an essentially constant difference between the binding energies of hydroxyl and peroxyl.9,15 Because this difference deviates from the ideal value, the optimal catalyst still possesses an Received: September 19, 2013 Revised: February 1, 2014 Published: February 3, 2014 4095

dx.doi.org/10.1021/jp409373c | J. Phys. Chem. C 2014, 118, 4095−4102

The Journal of Physical Chemistry C

Article

parameters, shown in Table 1, were obtained from the literature.24−31 Where available, DFT−PBE lattice parameters

overpotential of approximately 0.4 V for the OER, with no reaction barriers assumed. Studies on the scaling relationships of the OER intermediates have been carried out using mostly spin-restricted calculations. Given that most OER adsorbates are paramagnetic without the presence of a surface, this choice is only justified in the case that complete spin quenching occurs when the adsorbate is placed on the surface. In this case, the spatial occupation of the alpha and beta spin orbitals is delocalized to such an extent that the difference in eigenvalues for the alpha and beta spin orbitals vanishes. For the metals this has been proven true, from the study of both adsorbate spin quenching and adsorption energies.16 For metal oxides, however, spin quenching was often found to be incomplete.17−23 In the present study, we investigate whether incomplete spin quenching generally occurs for the OER intermediates and what consequences this has for the estimated catalytic performance of metal oxides. To study the differences between spin-restricted and unrestricted DFT for the reaction energies of *OH, *O, and *OOH, we have performed calculations using both formalisms on 10 nonmagnetic metal oxides. Rutile, rock salt, and perovskite structures were studied. Using the localization of the spin density, here defined as the difference between the alpha and beta spin densities, we have investigated the extent of spin quenching. Consequences of differences in restricted and unrestricted DFT calculations are discussed in terms of volcano plots.

Table 1. Structural Data of Studied Surfaces material

orientation

a0 (Å)

SnO2 SnO2 red. TiO2 PbO2 PtO2 MgO CdO SrTiO3 BaTiO3 NaTaO3

(110) (110) (110) (110) (110) (100) (100) (100) (100) (100)

4.745 4.745 4.667 5.060 4.484 4.216 4.689 3.94 3.995 5.484

b0 (Å)

c0 (Å)

ref

4.539

3.185 3.185 2.975 3.455 3.136

7.795

4.034 5.521

21 21 22 23 24 25 26 27 28 29

were used. Except for reduced SnO2, all surfaces were stoichiometric. For reduced SnO2, the bridging oxygen atoms were removed from the top and bottom layers of the slab. Two adsorbates per unit cell were used in all calculations. For rutile and rock salt surfaces, the adsorbates were placed on both sides of the slabs, at the vacant sites on top of metal atoms. Exception to this was reduced SnO2, for which the bridge site is more favorable. For perovskite structures, both adsorbates were placed on the same side of the slab to allow for identical adsorption sites on top of the B atoms. Due to the choice of two adsorbates, all unit cells contain an even number of electrons, allowing for a closed shell calculation. In the literature, this is often accomplished by allowing for fractional occupation of orbitals. In principle, our method for creating a closed shell system could lead to one adsorbate being negatively charged and one positively charged. This unphysical situation, however, was never observed in the Mulliken charges of the adsorbates. In unrestricted calculations the high-spin state was assumed. Unrestricted test calculations using the low-spin state in all cases showed identical results to their high-spin analogue, because of the negligible spin coupling between the wellseparated adsorbates. The total energy of the clean surfaces was calculated using restricted DFT. Unrestricted test calculations yielded identical values, as expected. Restricted DFT calculations were always performed with a closed shell. Reaction energies were calculated from total energy calculations of relaxed surface slabs using eqs 1−3.

■

METHODS To describe the rutile structures of the SnO2, reduced SnO2, TiO2, PbO2, and PtO2, a 4 × 2 unit cell of the (110) surface was used. A slab thickness of five trilayers was chosen, as is shown in Figure 2a. For the rock salt MgO and CdO structures,

Figure 2. Unit cells used for the rutile, rock salt, and perovskite stuctures with hydroxide adsorbates. The top structures show the first layer of the slabs. Bottom structures show a side view of the complete slab. (a) Rutile, (b) rock salt, (c) perovskite.

ΔEOH = (E HO − S − OH + E H2 − ES − 2E H2O)/2

(1)

ΔEO = (EO − S − O + 2E H2 − ES − 2E H2O)/2

(2)

ΔEOOH = (E HOO − S − OOH + 3E H2 − ES − 4E H2O)/2

(3)

In the above equations, ES is the total energy of the clean surface, whereas EHO−S−OH, EO−S−O, and EHOO−S−OOH denote the energy of the surface with OH, O, and OOH adsorbed on it, respectively. EH2 and EH2O are the total energies of the hydrogen and water molecule in vacuum. The division by 2 is necessary because two adsorbates were used. Note that electrons and protons released in the reaction are taken into account as hydrogen. With this choice, all reaction energies are referenced to the standard hydrogen electrode (SHE). To convert the reaction energies to free energies, the zero-point energy and entropy changes as given by Valdés et al.10 were added. The values used are shown in Table 2. It is assumed that

a 3 × 3 unit cell of the (100) surface was used in combination with a slab thickness of seven layers, as is depicted in Figure 2b. The perovskite structures (general composition ABO3) SrTiO3, BaTiO3, and NaTaO3 were described using a 4 × 4 unit cell of the (100) surface with a slab thickness of eight layers containing the A and B cations in an alternating fashion; see Figure 2c. The topmost layer consisted of B cations. The bottom four layers were fixed in their initial positions. All slabs were separated by a vacuum layer of at least 20 Å . The lattice 4096



Article

In these cases, we switched to the CG method. To ensure as much consistency in the calculation method as possible, no dipole corrections were used for BaTiO3. Wave function convergence was defined as having a maximum electronic gradient of 3 × 10−7a−1 0 and an energy difference between selfconsistent field cycles smaller than 10−13Eh. For the geometry optimization, convergence criteria of 4.5 × 10−4Eh/a0 for the maximum force component and 3 × 10−3a0 for the maximum displacement were used. Because of their difficult convergence, NaTaO3 and MgO calculations had to be truncated at a maximum force component of 5.4 × 10−3Eh/a0 and 2.4 × 10 −2 E h /a 0, respectively. To ensure the same level of convergence, the structures were compared for similarity.

Table 2. Entropic and Zero-Point Energy Contributions to the Free Energies of Gaseous and Adsorbed OER Species As Calculated by Valdés et al.10 H2O *OH *O *OOH H2

TS (eV)

ZPE (eV)

0.67 (at 0.035 bar) 0 0 0 0.41

0.56 0.35 0.05 0.41 0.27

the zero-point energy and entropy changes are surface independent. For the zero-point energy change, this is justified by the very similar numbers other authors have found for different surfaces using both restricted and unrestricted DFT.7,9,10 The entropy change is strongly dominated by the gas phase reactants and products and can therefore also be assumed surface independent. The total energy of water was calculated in the gas phase. When the entropy at the vapor pressure is used, however, the free energies are effectively referenced to liquid water. Binding energies were calculated in a similar way as the reaction energies. Total energies of the surface slabs containing two adsorbates were referenced to the total energies of the clean surface and the free adsorbates in vacuum. The energies of the adsorbates in vacuum were calculated using unrestricted DFT, using the high-spin state in all cases. All calculations were performed using the freely available CP2K/Quickstep simulation package.32,33 The Perdew− Berke−Ernsthof (PBE) fuctional was used to describe exchange-correlation effects.34 Orbitals were described by an atom-centered, double-ζ Gaussian basis set, while an auxiliary plane wave basis set with a 280 Ry cutoff was used to re-expand the electron density.35,36 The total energy calculations of the reference molecules H2, H2, OH, O, and OOH were performed using the FULL MOLOPT Gaussian basis set36 to compensate for the lower flexibility of basis sets for small systems. In all cases, the Goedecker−Teter−Hutter (GTH) pseudopotentials were employed to describe to core electron distribution.37 Where possible, the Pulay mixing (DIIS) electronic minimizer was used.38 When this led to numerical instability, we continued the calculation using the more stable conjugate gradients (CG) minimizer.39 The geometry was optimized using the Broyden−Fletcher−Goldfarb−Shanno (BFGS) algorithm. For MgO and NaTaO3, this yielded numerical instability.

■

RESULTS AND DISCUSSION In order to test the agreement between our results and those found in the literature, we first compare reaction energies obtained using restricted DFT with results from Man et al.,9 who used a similar approach for their calculations on nonmagnetic surfaces. Figure 3a shows our results as markers with fitted solid lines and a fit to the restricted DFT results of Man et al. as dashed lines. Clearly, the trends in reaction energies show good agreement. A direct comparison of our results for TiO2, SnO2, reduced SnO2, PtO2, and SrTiO3 yielded an acceptable mean absolute difference of 0.25 eV, most of which originated from PtO2 and SrTiO3. For a tabulated version of our data, we refer to the Supporting Information. Given the good agreement found here, we can conclude that the most notable difference between our method and that of Man et al., the absence of fractional orbital occupation, does not have a strong influence on the calculated reaction energies. Having established the validity of our restricted DFT results, we shall use them as a reference for our unrestricted DFT reaction energies. In Figure 3b, the reaction energies obtained from unrestricted DFT calculations are shown as markers with fitted solid lines and restricted DFT fits as dotted lines. Note that the restricted DFT *O reaction energy is still used as the xaxis, to allow for easy comparison to Figure 3a. For the *OH, almost no change is observed. Indeed, when the difference between restricted and unrestricted calculations is plotted, as shown in Figure 3c, the results are scattered around zero. For *OOH a clear difference is observed, increasing at higher reaction energies (right-hand side of the graph). The energy

Figure 3. Reaction energies of *OH, *O, and *OOH intermediates calculated using (a) restricted DFT and (b) unrestricted DFT. In both cases the x-axis is the *O reaction energy calculated using restricted DFT. (c) Difference between the reaction energies calculated using restricted and unrestricted DFT. Triangles, *O data; dots, *OH data; squares, *OOH data; solid lines, fit to data; dashed lines, fit to data from Man et al.;9 dotted lines, fit to data from (a). 4097



Article

differences for *O show an even steeper dependence on the *O reaction energy, leading to differences as large as 1.3 eV. In all cases the difference between restricted and unrestricted DFT reaction energies vanishes at low reaction energies (left-hand side of the graph). Note at this point that the surfaces with the highest reaction energies bind the intermediates the most weakly. Therefore, it is clear that the unrestricted DFT results lead to stronger binding (as it should), an effect that is most prominent for weak binding systems. To corroborate our unrestricted DFT calculations, we have performed calculations on surface slabs more similar to those used by Man et al.,9 containing only one adsorbate. Using four trilayer slabs with the two bottom layers fixed, we studied SnO2, reduced SnO2, PbO2, and PtO2. A comparison of the results obtained for this type of calculation to that obtained for the calculations with two adsorbates is given in Table 3. Clearly, the agreement is rather good, thus showing the reliability of our method.

surface and assume its vacuum state. For all intermediates considered here, this is a paramagnetic state. However, restricted DFT calculations do not allow for this. Therefore, the adsorbate can only leave the surface as an anion or in an excited state. Thus, the restricted DFT adsorbate is still bound at positive binding energies (referenced to the paramagnetic intermediates in vacuum), when the unrestricted DFT adsorbate leaves the surface. In this limit, the scaling relationship of the unrestricted versus restricted DFT reaction energy is a flat line. To take the limit of zero binding energy and the limit of vanishing difference between restricted and unrestricted calculations into account, we replaced the linear fit functions used in Figure 3a by eq 4.

Table 3. Comparison of Reaction Energies Obtained Using Unrestricted DFT Calculations with One Adsorbate and a Four Trilayer Slab and Calculations with Two Adsorbates and a Five Trilayer Slab (eV)

(4)

material SnO2 SnO2 red. PtO2 PbO2

ΔEOH 1 ad.

ΔEOH 2 ad.

ΔEO 1 ad.

ΔEO 2 ad.

ΔEOOH 1 ad.

ΔEOOH 2 ad.

1.55 −0.30

1.62 −0.42

4.33 0.83

4.44 0.48

3.74 3.39

4.39 3.20

0.23 2.03

0.26 2.14

2.39 4.76

2.53 4.71

3.40 4.69

3.48 4.74

⎛ ⎛ b − x ⎞⎞ ⎟⎟ f (x) = ax − ac ln⎜1 + exp⎜ ⎝ c ⎠⎠ ⎝ ⎛ ⎛ b − x ⎞⎞ ⎟⎟ + d + ac ln⎜exp⎜ ⎝ ⎝ c ⎠⎠

In this equation, a is the slope of the linear fit from the restricted DFT scaling relationship and d is the corresponding intercept. To constrain the function such that it reached the right maximum, reaction energy b was related to c such that limx →∞ f (x) = ΔEmax unrestricted

(5)

Here the maximum reaction energies are 5.24, 2.93, and 4.91 eV for *O, *OH, *OOH, respectively, the values of the reaction energy when it occurs in the gas phase. The resulting fit function contains only one fitting parameter. Even with this low flexibility, the fit agrees very well with the data, justifying our two-limit interpretation. However, other integrated sigmoid functions may give similar results, making further interpretation based on the fitting function difficult.

From Figure 3b it can be seen that the scaling relationships of *O and *OOH level off at high reaction energy, i.e., low binding energy. This can be explained if one considers the case of zero binding energy. In this case, the adsorbate will leave the

Figure 4. Relation between spin density localization and differences between restricted and unrestrited DFT calculations. The isosurfaces of the spin density of (from left to right) reduced SnO2, SrTiO3 , and MgO are shown at a density of 5 × 10−3, looking along the y-axis. Also shown is a plot of the total spin moment of the adsorbates, obtained from a Mulliken population analysis, versus its binding energy. Triangles, *O; dots, *OH; squares, *OOH. The solid lines represent Gaussian fits to the data, where only the width was used as a fitting parameter. A single fit was used to fit both *OH and *OOH data. 4098



Article

To give a more eleborate explanation for the differences in reaction energies observed in our comparison of unrestricted and restricted DFT, we shall interpret the data in terms of spin localization effects. In the limit of zero binding energy, the adsorbate will leave the surface. In this case, the spin density will be fully localized on the adsorbate. When the binding energy is going negative, the degree of hybridization of the adsorbate states with the substrate states increases. The hybridized states are more delocalized and will therefore possess a more diluted spin density; i.e., spin quenching will start to occur. Because the hybridization of adsorbate and surface states is not the only contribution to binding energy, the trend that follows when going to more and more negative binding may be somewhat scattered. Yet one can state that in the limit of very strong binding, strong hybridization is to be expected. Therefore, complete spin quenching will occur on very strong binding surfaces, thus erasing the difference between alpha and beta spin orbital eigenvalues. Hence there will be no difference between restricted and unrestricted DFT calculations in this limit. Alternatively, one can view this process from an ionic point of view.13 In this case, one can regard the transition from spin-localized to spin-quenched systems as a transition from covalent bonding to ionic bonding. The picture described above suggests that the spin density should delocalize with increasing binding energy. In Figure 4, the spin density of three surfaces is shown in order of decreasing binding strength. For the leftmost surface, reduced SnO2, the spin density distribution is very delocalized. At intermediate binding strength, represented here by SrTiO3, the spin density is still rather delocalized, but confined to the first surface layer. For the weak binding MgO, the spin density is almost completely localized on the adsorbate. These observations corroborate the connection between binding energy and spin localization. From the plot of adsorbate spin moment versus binding energy, shown in Figure 4, this is also clear. One has to note, though, that results are scattered somewhat. This could be expected based on the analysis in the last paragraph. The observation of various degrees of spin quenching, depending on the nature of the binding to the substrate, is in good agreement with literature reports on the adsorption of other adsorbates on metal oxides.17−23 Zhang et al. showed that SO2 can have a spin moment of 0 or 2, depending on the adsorption configuration. For metal atoms and clusters adsorbed on metal oxides, various others also reported a variety of spin-quenching behavior. To quantify the influence that the spin localization has on the energy difference in restricted and unrestricted DFT reaction energies, Figure 5 shows a plot of the adsorbate spin moment versus the observed energy difference. As expected, a clear correlation is observed. Note that all adsorbates seem to have a very similar response in energy difference to spin localization. Given this similar response, the origin of the much larger energy differences observed for *O relative to *OOH and *OH is the fact that oxyl possesses two unpaired spins in vacuum versus one for *OOH and *OH. The different restricted/unrestricted DFT response that is observed for *OH and *OOH can be explained by their difference in binding energies. From Figure 4 it is clear that hydroxyl and peroxyl show a similar degree of spin localization when bound equally strong. However, when bound to the same electrode, hydroxyl will always bind stronger. Therefore, it shows a lower degree of spin localization, and thus a smaller

Figure 5. Relation between adsorbate spin moment difference in reaction energy for restricted and unrestricted DFT calculations and total spin moment of the adsorbates as obtained from a Mulliken population analysis versus the observed . Triangles, *O; dots, *OH; squares, *OOH.

difference between restricted and unrestricted DFT calculations. Having explained the trends in the differences between the reaction energies calculated using restricted and unrestricted DFT, we shall now discuss their consequences for the estimated catalytic activity. In Figure 6a, the free energy scaling relationships calculated using unrestricted DFT are shown for all adsorbates. Assuming that either the transition from *OH to *O or the transition from *O to *OOH is the potentialdetermining step, the thermodynamic overpotential is defined as9 η = max(ΔGO − ΔGOH , ΔGOOH − ΔGO) − 1.23

(6)

With this definition, the thermodynamic overpotential gives the largest deviation from the ideal energy cost per electron transfer. This gives a lower limit to the kinetic overpotential measured experimentally. A good catalyst therefore possesses a low thermodynamic overpotential. Given eq 6 and Figure 6a, it follows that the thermodynamic overpotential has a pyramidal dependence on all of the adsorbate free energies. Therefore, they could all serve as the descriptor, the x-axis of a volcano plot. A more convenient choice, however, is ΔGO − ΔGOH, as this choice allows for cancellation of the possible systematic errors made in calculating ΔGO and ΔGOH. Furthermore, the shape of the volcano plot is now only determined by ΔGOOH − ΔGO. The reason for this is that at the right side of the pyramid η = ΔGO − ΔGOH − 1.23, so the trend is fixed. Because of the differences in reaction energies calculated using restricted and unrestricted DFT, one would also suspect differences in the estimated catalytic activity. A volcano plot comparing these estimations, using ΔGO − ΔGOH as a descriptor, is shown in Figure 6b. Clearly, the shapes of the restricted and unrestricted DFT volcano plots are very similar. The reason for this can found in Figure 6c, where the differences in the thermodynamic overpotential are shown. On the left-hand side and around the center of the volcano, the differences are scattered around zero. On the right-hand side, there are marked differences. However, on this side of the volcano the trend is fixed by definition. Note that when one is interested in the activity of a specific material rather than the trend, using unrestricted DFT can of course make a significant difference. 4099



Article

Figure 6. (a) Free energy scaling relationships calculated using unrestricted DFT. Triangles, *O data; dots, *OH data; squares, *OOH data; solid lines, fit to data. (b) Volcano plots constructed from restricted DFT (squares) and unrestricted DFT free energies (triangles). The x-axis uses restricted DFT data for the restricted DFT volcano plot and unrestricted DFT data for the unrestricted DFT volcano plot. (c) Difference in thermodynamic overpotential calculated using restricted and unrestricted DFT.

adsorbate. The added electron density on the adsorbate will cause (partial) spin quenching, because the adsorbate tends toward its negative ion (OH−, O2−, OOH−), which has a closed shell. Note that this is a different spin-quenching mechanism than that observed here. For our substrates, the oxygen atoms of a large part of the substrate cause spin quenching rather than a single metal atom. This difference may affect the binding energy at which the transition from a spin-localized to a spinquenched adsorbate occurs. We should note at this point that although unrestricted calculations may capture more of the essential chemistry of the metal oxide−OER adsorbate system, it is not our aim in the present work to give accurate predictions for the catalytic activities of the systems under study. Apart from the fact that our model is crude, for instance, due to the absence of a solvent, the PBE functional is known to have severe deficiencies.41 Most notably, DFT−PBE calculations overestimate the degree of delocalization of states. As we have stressed above, the delocalization of states affects the spin localization. Therefore, PBE may underestimate the effects of spin localization on the volcano plot. It may well be that the transition from spin-localized to spin-quenched systems is located further left on the volcano plot. In this case, the shape of the volcano plot could be affected. More importantly, it could also affect the prediction of the catalytic activity of materials at the top of the volcano, which are the best electrocatalysts. The use of a better quality method may at least partially resolve the overdelocalization problem and therefore give more trustworthy results. Recently, PBE+U is gaining popularity as a possible solution for overdelocalization errors. However, the U is typically only applied to the metal atoms and therefore has little effect on the localization behavior of the states that are almost completely localized on the oxygen atoms of the system. Because the spin localization phenomena are usually dominated by the oxygen atoms, PBE+U may not improve the spin localization much. Alternatively, a hybrid exchange-correlation functional can be used. Although expensive, a hybrid functional improves the localization of states in a general and relatively unambiguous way. Future work will have to show how this affects the results found here.

The observed differences in the estimation of the catalytic activity are particularly important for photocatalysts, which are typically weak binding systems that appear on the right side of the volcano plot. As an example, the overpotential for CdO is 0.73 eV lower when calculated using unrestricted rather than restricted DFT. In studies on photocatalysts that used restricted GGA−DFT, calculations on the bare, H2O-covered, and OHcovered surfaces remain valid. Yet data on *O and *OOH should be reconsidered. Even if the reaction mechanism is not the one assumed here, the intermediates will be of a similar type. Therefore, they will also be prone to spin localization, giving rise to a necessity to use unrestricted DFT. From the similarity in shape of the restricted and unrestricted DFT volcano plots, we can understand why other authors that have worked on the OER or oxygen reduction reaction with unrestricted DFT did not report deviations from restricted DFT results. To the best of our knowledge, no unrestricted calculations have been performed on metals or metal oxides that fall on the far right side of the pyramid,4−7,9 where differences are largest. For metals, it would seem that there generally is no possibility for obtaining very small binding energies. The reason for this, as the Hammer−Nørskov model explains,40 is that all metals have coupling of the s-band to the adsorbate, which yields a considerable minimum for the binding energy. This minimum seems to be sufficient to cause spin quenching.16 Metal oxides, which do not have an s-band, are able to show much weaker binding. The binding energy found in our study for *O on CdO, for instance (unrestricted DFT), is some 1.8 eV smaller than the smallest binding energy found by Nørskov et al.4 (for Au). Another important factor that has obscured the importance of unrestricted DFT in the calculations on OER intermediates on metals oxides is that it is almost exclusively applied to magnetic surfaces. For these, no direct comparison between restricted and unrestricted DFT data is possible, so that only the trend in the volcano plot can be used for comparison. Similar arguments hold for the scaling relationships, where again none of the studied substrates fall in the weak binding region. Indeed from the data of Carter et al.7 it can be deduced that spin quenching of the adsorbates occurred on their substrates. It was shown that the spin moment on the dopants in doped hematite changes significantly upon adsorption of OH, O, or OOH (but not H2O), with corresponding oxidation of the dopant by the 4100



■

Article

CONCLUSION Using both restricted and unrestricted DFT, we have calculated the reaction energies of three key OER intermediates on 10 nonmagnetic metal oxides. We show that significant differences exist between the restricted and unrestricted DFT results. The differences are most pronounced for the *O intermediate, for which the unrestricted DFT result was maximally 1.3 eV lower in energy. For *OOH, this maximum is 0.5 eV. For *OH, almost no changes are observed. The difference between the restricted and unrestricted DFT reaction energies decreases with increasing binding strength, vanishing for strong binding systems. The observations can be explained in terms of the spin localization on the adsorbate. For weak binding systems, the adsorbates retain the paramagnetic character they have in vacuum. This situation is not well described by restricted DFT, yielding higher energies than its unrestricted counterpart. For strong binding systems, however, spin quenching occurs. In this case the difference between restricted and unrestricted DFT calculations vanishes. The different response of *OH, *O, and *OOH toward the use of unrestricted DFT originates from their different numbers of unpaired electrons and their differences in binding energies. Because *O has two unpaired spins, compared to one for *OH and *OOH, it benefits most of the use of unrestricted DFT. The difference between *OH and *OOH originates from the higher binding energy of *OH. Because *OH binds stronger, it shows more spin quenching. Therefore, *OH shows practically no difference between restricted and unrestricted DFT calculations, whereas a clear difference is observed for *OOH. Volcano plots, which estimate the catalytic activity of surfaces, show a very similar shape when constructed from restricted and unrestricted DFT results with ΔGO − ΔGOH as a descriptor. However, the location of data points on the curve does change, particularly for weak binding substrates. Because photocatalysts typically belong to this category, the use of unrestricted DFT is vital for the study of this class of materials. Because of the well-known delocalization error of the PBE functional,41 the transition from spin-localized to spinquenched systems may in reality occur at higher binding energies, possibly also affecting the best electrocatalysts at the top of the volcano plot. To show whether this is indeed the case, studies using more accurate exchange-correlation potentials like hybrid functionals will be necessary.

■

Emmanuel College at Cambridge for a research fellowship. The calculations in this work have been performed using an allocation of computer time on HECToR, the U.K.’s high-end computing resource funded by the Research Councils, as part of a grant to the UKCP consortium.

■

(1) Dau, H.; Limberg, C.; Reier, T.; Risch, M.; Roggan, S.; Strasser, P. The Mechanism of Water Oxidation: From Electrolysis via Homogeneous to Biological Catalysis. ChemCatChem 2010, 2, 724− 761. (2) Maeda, K.; Domen, K. Photocatalytic Water Splitting: Recent Progress and Future Challenges. J. Phys. Chem. Lett. 2010, 1, 2655− 2661. (3) Stamenkovic, V. R.; Fowler, B.; Mun, B. S.; Wang, G.; Ross, P. N.; Lucas, C. A.; Markovic̀, N. M. Improved Oxygen Reduction Activity on Pt3Ni(111) via Increased Surface Site Availability. Science 2007, 315, 493−497. (4) No̷ rskov, J. K.; Rossmeisl, J.; Logadottir, A.; Lindqvist, L. Origin of the Overpotential for Oxygen Reduction at a Fuel-Cell Cathode. J. Phys. Chem. B 2004, 108, 17886−17892. (5) Fang, Y.-H.; Liu, Z.-P. Mechanism and Tafel Lines of ElectroOxidation of Water to Oxygen on RuO2(110). J. Am. Chem. Soc. 2010, 132, 18214−22. (6) Chen, J.; Selloni, A. Water Adsorption and Oxidation at the Co3O4(110) Surface. J. Phys. Chem. Lett. 2012, 3, 2808−2814. (7) Liao, P.; Keith, J. A.; Carter, E. A. Water Oxidation on Pure and Doped Hematite (0001) Surfaces: Prediction of Co and Ni as Effective Dopants for Electrocatalysis. J. Am. Chem. Soc. 2012, 134, 13296−309. (8) Valdés, A.; Brillet, J.; Grätzel, M.; Gudmundsdóttir, H.; Hansen, H. A.; Jónsson, H.; Klüpfel, P.; Kroes, G.-J.; Le Formal, F.; Man, I. C. Solar Hydrogen Production with Semiconductor Metal Oxides: New Directions in Experiment and Theory. Phys. Chem. Chem. Phys. 2012, 14, 49−70. (9) Man, I. C.; Su, H.-Y.; Calle-Vallejo, F.; Hansen, H. A.; Martìnez, J. I.; Inoglu, N. G.; Kitchin, J.; Jaramillo, T. F.; No̷ rskov, J. K.; Rossmeisl, J. Universality in Oxygen Evolution Electrocatalysis on Oxide Surfaces. ChemCatChem 2011, 3, 1159−1165. (10) Valdés, A.; Qu, Z.-W.; Kroes, G.-J.; Rossmeisl, J.; No̷ rskov, J. K. Oxidation and Photo-Oxidation of Water on TiO2 Surface. J. Phys. Chem. C 2008, 2, 9872−9879. (11) Rossmeisl, J.; Qu, Z.-W.; Zhu, H.; Kroes, G.-J.; No̷ rskov, J. K. Electrolysis of Water on Oxide Surfaces. J. Electroanal. Chem. 2007, 607, 83−89. (12) Valdès, A.; Kroes, G.-J. First Principles Study of the PhotoOxidation of Water on Tungsten Trioxide (WO3). J. Chem. Phys. 2009, 130, 114701. (13) Cheng, J.; Sulpizi, M.; Vandevondele, J.; Sprik, M. Hole Localization and Thermochemistry of Oxidative Dehydrogenation of Aqueous Rutile TiO2(110). ChemCatChem 2012, 2, 636−640. (14) Rossmeisl, J.; Logadottir, A.; No̷ rskov, J. K. Electrolysis of Water on (Oxidized) Metal Surfaces. Chem. Phys. 2005, 319, 178−184. (15) Koper, M. T. M. Thermodynamic Theory of Multi-Electron Transfer Reactions: Implications for Electrocatalysis. J. Electroanal. Chem. 2011, 660, 254−260. (16) Fajìn, J. L. C.; Cordeiro, D. S.; Natàlia, M.; Gomes, J. R. B.; Illas, F. On the Need for Spin Polarization in Heterogeneously Catalyzed Reactions on Nonmagnetic Metallic Surfaces. J. Chem. Theory Comput. 2012, 8, 1737−1743. (17) Yudanov, I.; Pacchioni, G.; Neyman, K.; Rösch, N. Systematic Density Functional Study of the Adsorption of Transition Metal Atoms on the MgO(001) Surface. J. Phys. Chem. B 1997, 5647, 2786− 2792. (18) Cruz Hernandez, N.; Màrquez, A.; Sanz, J. F.; Gomes, J. R. B.; Illas, F. Density Functional Theory Study of Co, Rh, and Ir Atoms Deposited on the r-Al2O3(0001) Surface. J. Phys. Chem. B 2004, 3, 15671−15678.

ASSOCIATED CONTENT

S Supporting Information *

Tabulation of all calculated reaction energies and adsorbate spin moments; description of the conversion from reaction energy to free energy change. This material is available free of charge via the Internet at http://pubs.acs.org/.

■

REFERENCES

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Phone: +44(0)1223 336314. Fax: +44(0)1223 336422. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS R.V.M. thanks Leiden University, the Leiden University Fund. and the Erasmus program for financial support. J.C. thanks 4101



Article

(19) Markovits, A.; Paniagua, J.; Lòpez, N.; Minot, C.; Illas, F. Adsorption Energy and Spin State of First-Row Transition Metals Adsorbed on MgO(100). Phys. Rev. B 2003, 67, 115417. (20) Di Valentin, C.; Pacchioni, G.; Bredow, T.; Dominguez-Ariza, D.; Illas, F. Bonding of NO to NiO(100) and NixMg1−xO(100) Surfaces: A Challenge for Theory. J. Chem. Phys. 2002, 117, 2299. (21) Zhang, C.; Lindan, P. J. D. A Density Functional Theory Study of Sulphur Dioxide Adsorption on Rutile TiO2(110). Chem. Phys. Lett. 2003, 373, 15−21. (22) Helali, Z.; Markovits, A.; Minot, C.; Abderrabba, M. First-Row Transition Metal Atoms Adsorption on Rutile TiO2(110) Surface. Struct. Chem. 2012, 23, 1309−1321. (23) Sousa, C.; Tosoni, S.; Francesc, I. Theoretical Approaches to Excited-State-Related Phenomena in Oxide Surfaces. Chem. Rev. 2013, 113, 4456−4495. (24) Evarestov, R. A.; Bandura, A. V.; Proskurov, E. V. Plain DFT and Hybrid HF-DFT LCAO Calculations of SnO2 (110) and (100) Bare and Hydroxylated Surfaces. Phys. Status Solidi B 2006, 243, 1823−1834. (25) Cheng, J.; Sprik, M. Acidity of the Aqueous Rutile TiO2(110) Surface from Density Functional Theory Based Molecular Dynamics. J. Chem. Theory Comput. 2010, 2, 880−889. (26) Scanlon, D.; Kehoe, A.; Watson, G.; Jones, M.; David, W.; Payne, D.; Egdell, R.; Edwards, P.; Walsh, A. Nature of the Band Gap and Origin of the Conductivity of PbO2 Revealed by Theory and Experiment. Phys. Rev. Lett. 2011, 107, 246402. (27) Range, K.-J.; Rau, F.; Klement, U.; Heyns, A. M. High Pressure Synthesis of Single Crystals. Mater. Res. Bull. 1987, 22, 1541−1547. (28) Zhu, Y.; Chen, G.; Ye, H.; Walsh, A.; Moon, C.; Wei, S.-H. Electronic Structure and Phase Stability of MgO, ZnO, CdO, and Related Ternary Alloys. Phys. Rev. B 2008, 77, 245209. (29) Piskunov, S.; Heifets, E.; Eglitis, R. I.; Borstel, G. Bulk Properties and Electronic Structure of SrTiO3, BaTiO3, PbTiO3 Perovskites: An Ab Initio HF/DFT Study. Comput. Mater. Sci. 2004, 29, 165−178. (30) Salehi, H.; Shahtahmasebi, N.; Hosseini, S. M. Band Structure of Tetragonal BaTiO3. Eur. Phys. J. B 2003, 32, 177−180. (31) Ahtee, M.; Darlington, C. Structures of NaTaO3 by Neutron Powder Diffraction. Act. Crystallogr. 1980, 36, 1007−1014. (32) The CP2K developers group, http://www.cp2k.org, accessed November 25th, 2012. (33) VandeVondele, J.; Krack, M.; Mohamed, F.; Parrinello, M.; Chassaing, T.; Hutter, J. QUICKSTEP: Fast and Accurate Density Functional Calculations Using a Mixed Gaussian and Plane Waves Approach. J. Comput. Phys. Commun. 2005, 167, 103−128. (34) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865−3868. (35) Lippert, G.; Hutter, J.; Parrinello, M. A Hybrid Gaussian and Plane Wave Density Functional Scheme. Mol. Phys. 1997, 92, 477− 487. (36) VandeVondele, J.; Hutter, J. Gaussian Basis Sets for Accurate Calculations on Molecular Systems in Gas and Condensed Phases. J. Chem. Phys. 2007, 127, 114105. (37) Goedecker, S.; Teter, M.; Hutter, J. Separable Dual-Space Gaussian Pseudopotentials. Phys. Rev. B 1996, 54, 17031710. (38) Hamilton, T. P.; Pulay, P. Direct Inversion in the Iterative Subspace (DIIS) Optimization of Open-Shell, Excited-State, and Small Multiconfiguration SCF Wave-Functions. J. Chem. Phys. 1986, 184, 5728−5735. (39) VandeVondele, J.; Hutter, J. An Efficient Orbital Transformation Method for Electronic Structure Calculations. J. Chem. Phys. 2003, 118, 4365−4369. (40) Hammer, B. Special Sites at Noble and Late Transition Metal Catalysts. Top. Catal. 2006, 37, 3−16. (41) Cohen, A. J.; Mori-Sànchez, P.; Yang, W. Challenges for Density Functional Theory. Chem. Rev. 2012, 112, 289−320.

4102


Modeling the Oxygen Evolution Reaction on Metal Oxides: The

Recommend Documents