Simulations to Cover the Waterfront for Iron Oxide Catalysis.

Hematite has been widely studied for catalytic water splitting, but the role of the interactions between catalytic sites is unknown. In this paper, we calculate the oxygen evolution reaction free energies and the surface adsorption distribution using a combination of density functional theory and Monte Carlo simulations to "cover the waterfront," or cover a wide range of properties with a simulation of the hematite surface under working conditions. First, we show that modeling noninteracting catalytic sites provides a poor explanation of hematite's slow reaction kinetics. The interactions between the catalytic site may hinder catalysis through the strong interactions of *OH2 and *OOH intermediates, which cause the reaction to revert back to the *O intermediate. Hence, neighboring interactions may be a possible reason for the abundant, experimentally observed *O intermediate on the surface. This study demonstrates how neighboring sites impact the energy required for catalytic steps, thus providing new avenues to improve catalysis by controlling neighboring site interactions.

Introduction

Hydrogen is an essential part of today's economy.[ , ] Current hydrogen production methods are highly polluting, thus encouraging the search for cleaner production methods. One method for producing hydrogen is water splitting in an electrochemical cell (PEC).[ , , , ] In a PEC, water is split into hydrogen and oxygen. The anode of the PEC produces oxygen by the oxygen evolution reaction (OER) or water oxidation. One of the biggest challenges in water splitting is finding an efficient catalyst for the OER. The oxygen evolution reaction requires the transfer of four electrons, making the ideal catalyst hard to find.

Water oxidation in a PEC requires a catalyst. One of the most studied catalysts is α‐Fe2O3 or hematite. Hematite is stable under alkaline conditions and has a favorable band gap for water splitting but suffers from slow OER kinetics. Though there have been many studies of hematite, there is still more to understand, such as the interactions of termination species.

The OER under alkaline conditions can be modeled in five reactions:

where * denotes a surface vacancy and *X denotes an adsorbed surface species.

Previous works concentrated on single‐site reactions[ , ] and two‐site species interactions that form new intermediates and different reactions. First, previous works calculated the free energy differences between two consecutive reactive species. Later, larger slabs were used to show that when terminations are far apart, they do not significantly affect the calculated free energy. In this paper, we investigate the effect of near‐neighbor sites on catalysis using a computational density functional theory method and Monte Carlo simulations. Since switching the location of two species does not affect the total energy, the number of slabs was reduced from 25 to 15. The slabs were created using MATLAB by combining four smaller slabs comprised of the five reaction intermediates. In each 2×2 slab, two intermediates were connected to a single iron atom. All other sites were taken to be hydroxyl terminated *OH sites, since the hydroxyl termination is the initial intermediate before the onset of catalysis. An example slab is portrayed in Figure 1. The lattice constants of each slab were a=b=10.20 Å, c=24.28 Å, α=β=90°, and γ=120°. We used the (0001) surface of hematite, which is one of the natural growth surfaces of hematite.

Figure 1

(top) Example slab of neighboring *OH2−*OOH sites. Surface *OH2 and *OOH sites circled. (bottom) An example cycle with *O as a neighbor and its free energy differences at pH=0.

Computational Details

Calculations were performed in VASP 5.4.4[ , , , ] using the Perdew‐Burke‐Ernzerhof (PBE) exchange‐correlation functional. The spin‐polarized DFT+U approach by Dudarev et al. was used with a U−J value of 4.3 eV for iron, as was derived by ab initio methods. Values for oxygen and hydrogen were set to zero. Core electrons were described by projector augmented wave (PAW) potentials.

The Kohn‐Sham equations were solved self‐consistently using a plane‐wave basis and three‐dimensional periodic boundary conditions. K‐space integration was performed using the tetrahedron method with Blöchl correction.[ , ] The plane‐wave energy cutoff was 700 eV. Gamma‐centered k‐grids of 2×2×1 were used for the slabs. These parameters converge the total energy to within 1 meV. Geometric relaxations were performed using the conjugate gradient (CG) method with a tolerance of 0.03 eV/Å for atomic forces.

For the Monte Carlo simulation, we created a simulated n‐by‐n grid of hexagonally shaped active sites. Each surface intermediate created a hexagonal grid of sites, as shown in Figure 2. The active site is at the center of the hexagons, where a water molecule would adsorb on an oxygen vacancy site. In this work, we chose a 30‐by‐30 grid for all simulations.

Figure 2

Hematite surface with a hexagonal grid. The blue site represents a chosen site, green sites are close neighbors that share an iron atom, and red sites are close neighbors that do not share an iron atom.

Each active site has two bonds with neighboring iron atoms, as shown in Figure 2. Each iron atom has three bonds to active sites. In total, each active site (marked in blue) has six close neighbors, but only four of them (marked in green) share an iron atom with the site. Two sites (marked in red) do not share an iron atom with the active site, but were still considered first‐order neighbors in the calculation. Such neighboring relation is called an “interaction” in this paper, and we will show that such interactions affect the reaction energy. Under the “non‐interacting” sites model, reaction sites are unaffected by neighboring intermediates.

At each Monte Carlo iteration, a random site in the grid is chosen. The free energy change of the relevant catalysis reaction is calculated by taking the difference in the energy of the relaxed structure, gases, entropy, and zero‐point energy (ZPE) as used in other works[ , ] and shown in equation (6) for reaction (3) as an example. The effect of pH on the free energy is calculated from the Nernst equation. The external bias is modeled as a constant electric field across the slab.

In this equation, E is the relaxed‐structure energy as calculated by DFT, ZPE is the zero‐point energy, T is the absolute temperature, S is entropy, and U is the external bias relative to the RHE. The energy of hydrogen represents the proton that leaves *OH in the reaction. Using Equation (6), we calculated the free energy changes at pH=13.6 for all five reactions, as shown in Table 1.

Table 1

Free energy change of single‐site reactions at pH=13 and no bias with respect to SHE. The sum of free energy changes is not 1.7 eV, as theory suggests, since DFT calculations give a water oxidation potential of 1.11 V and not the theoretical 1.23 V.

Reaction Number	ΔG at pH=13.6 and no bias w.r.t SHE [eV]
(1)	0.09
(2)	−0.84
(3)	1.05
(4)	0.84
(5)	0.08

The slab energy difference was calculated with regard to near‐neighbors as described for single site reactions and averaged over the six neighbors. The ZPE and entropy values were assumed to be unaffected by neighbor interactions. Then, the transition probability was calculated as follows:

The simulation was run for a preset number of iterations. Each iteration result was saved in a MATLAB array.

Since the interaction between *OOH and *OH2 is a result of geometric relaxation, we assumed it has no kinetic barrier, so the transformation occurs automatically whenever *OOH and *OH2 sites are neighbors, regardless of the chosen site.

In addition to the effect of thermodynamics, kinetics is a key factor under working conditions. To measure the effect of incorporating kinetic barriers into the model, we considered the activation energy from Marcus, as shown in Equation (8) for Reactions (2)‐(5), which involve electron transfer.

where λ is the reorganization energy. We used a value of 1.5 eV from Graetzel's work on iron oxide, which involves the inner part of the reorganization energy. ΔG is calculated per reaction, as described earlier.

Since kinetic barriers involve an activation energy, we used Kinetic Monte Carlo (KMC) with the residence time algorithm (RTA) instead of the Metropolis algorithm. With RTA, we calculated reaction rates of all possible reactions on the grid by using Equation (9), where Ea,i is the activation energy of the reaction in the i′th site. Then we created a sequence of partial sums and define the sum of all N rates as R (Equation (10)). We assumed an identical k0 for all reactions, with a value of , as was done by Rajan and Carter.[ , ]

After calculating all sums, we draw a random number u from a uniform distribution and choose the reaction in site j, where j satisfies the term in Equation (12). The chosen reaction occurs with 100 % acceptance rate. For example, if the simulation has three sites with k1=2 Hz, k2=3 Hz, and k3=1 Hz, we get R1=2 Hz, R2=5 Hz, and R3=R=6 Hz. Then, to pick site 1, u needs to be between 0 and 1/3, a probability of . To pick site 2, u needs to be between 1/3 and 5/6, a probability of , and to pick site 3, u needs to be between 5/6 and 1, a probability of . In this way, the largest probability is to choose site 2, since it has the largest rate reflected in the largest range of u values between 1/3 and 5/6 Therefore, this method allows accounting for the relative rates of all sites.

Results and Discussion

We will present the results of hematite surface coverage by reaction intermediates in two parts. First, we will present the results of Monte Carlo simulations without the effects of near‐neighbor interactions between the species on the active sites. Then, we will present the results with near‐neighbor interactions, especially the effect of the interaction between the two intermediates *OOH−*OH2.

Without Neighbor Interactions

When sites do not interact with each other, the steady‐state configuration of the system can be calculated from transition probabilities:

where θi is the fraction of sites with the i'th species and is the probability that a catalysis reaction would occur should a site with the i'th species be chosen. The magnitude of probabilities determines the number of iterations required to achieve a steady state.

Under standard operating conditions, T=298.15 K, pH=13.6, and no bias, all transition probabilities are very close to zero since the free energy difference of Reaction (3) is approximately 1 eV and . An increased bias is required to lower the free energy difference to initiate catalysis in the simulation and in the experiment.

Increasing the bias to 1 V lowers all free energy differences below zero except for Reaction (3), which has with no bias, and Reaction (1), which is independent of bias and pH. As a result, vacancies and *OH sites dominate the surface, while *OH2, *O, and *OOH sites have approximately equal coverages, as shown in Figure 3.

Figure 3

Number of intermediate species distributions without neighboring interactions averaged over 1000 Monte Carlo steps (top) and spatial distribution of intermediates (bottom). Simulations were run at pH=13.6 and a simulated voltage of 1.5 V.

When the bias is increased beyond 1.06 V, all free energy differences except for reaction (1) become negative, so the ratio of vacancies to other species becomes approximately 30 : 1 : 1 : 1 : 1 since and . A further increase in bias has no effect since the probabilities are unaffected when the free energy differences remain below zero. Such a high coverage of surface vacancies does not agree with the experimental results and may mean that neighboring interactions need to be considered.

With Neighbor Interactions

When DFT was used to simulate the *OOH and *OH2 sites as near‐neighbors, they immediately transformed into *O and *OH sites, respectively. The transformation resulted from geometric relaxation of a slab with neighboring *OOH and *OH2 sites, which was unstable and resulted in a state with *O and *OH sites and a desorbed water molecule.

This interaction of species occurs due to two weak bonds: the O−H bond of *OH2 and the O−O bond of *OOH. The free energy change of the transition of *OH2 to *OH is close to zero with most neighbors when pH=0 and negative when pH=13.6, as shown in Table 2, which means it is almost energetically equivalent to that of *OH. The free energy change of the conversion of *O to *OOH is over 1 eV and the highest of all the reactions with some neighbors, which means that removing the extra OH is energetically favorable. Since *OH2 and *OOH are the largest species with three adsorbed atoms each, they can interact spatially more easily.

Table 2

Free energy differences (in eV) of all reactions with a near‐neighbor constant species. Numbers were calculated with T=298.15 K, pH=13.6, and no applied bias. (a) Reactions involving the *OOH−*OH2 neighbors that automatically result in *OH−*O neighbors after geometry optimization. (b) Reactions involving *OOH−*OH2 that will not occur since the interaction occurs before the reaction.

Constant species Reaction	*	*OH2	*OH	*O	*OOH
*→OH2	−0.072	−0.082	−0.007	−0.253	0.153a
*OH2→*OH	−0.690	−0.615	−0.982	0.339	−0.867b
*OH→*O	−0.183	−0.062	1.259	0.881	1.312
*O→*OOH	1.558	1.964a	0.758	1.189	0.544
*OOH→*	0.241	0.017b	0.194	−0.936	0.078

*

To show that neighboring species are possible, we calculated the formation energies of all 15 slabs. The formation energies of the five slabs with *OH are zero, since we take the *OH intermediate to be the “neutral” intermediate. All other formation energies are negative, except for the formation of *O−*OOH intermediates sharing an iron atom. However, the formation energy of the *O−*OOH slab is 0.106 eV, the lowest, in absolute values, of all formation energies, so the creation of *O and *OOH intermediates that share an iron atom is not impossible. All formation energies and calculation method are detailed in the Supporting Information and Table S3.

The addition of neighboring active site adsorbate interactions changes the simulation in two ways. The first is that the transition probabilities change according to the free energy changes detailed in Table 2, which may change the steady‐state coverage of the iron oxide surface. The second is the addition of the strong *OOH−*OH2 interaction that inhibits further reactivity for the participating species. The stability of the *OOH−*OH2 interaction becomes the only relevant reaction that affects the progress of water oxidation, since at increases simulated external bias, the free‐energy differences decrease or negative, making all transition probabilities closer and eventually equal to one, except for Reaction (1), which has bias‐independent free energy.

The free energy differences in Table 2 reveal possible reaction paths in which neighboring species enable reactions at lower voltages. For example, when going from *O−*OH neighbors to *V−*OOH neighbors, there are several possible reaction paths, as shown in Figure 4. However, some paths undergo reactions with higher free‐energy differences, which increases the overpotential.

Figure 4

Different reaction paths from *OH−*O to *‐*OOH with pH=13.6. The four reaction paths are created by choosing one path between (a) and (b) and one between (c) and (d). The intermediates in each step are shown in labels. * alone denotes a vacant site.

When the simulated applied voltage is higher than all free energy differences (exemplified by the extreme value of 3 V), all reaction probabilities go to one except for that of Reaction (1) next to *OH2. After several thousand cycles of Monte Carlo simulations, such as the one shown in Figure 5 at a steady state, there is a majority of *O intermediates.

Figure 5

Average simulated distribution of species by number over 1000 frames (top). Spatial distribution of termination species with a simulated voltage of 3 V and pH=13.6 (bottom). Both simulations were run with 30×30 sites for 10,000 steps.

The number of *OH and *OOH intermediates is almost equal and amounts to approximately half of the number of *O intermediates. According to the simulation results, *OH2 accounts for a small portion of the surface, and vacancies account for approximately 60 % of the number of *O intermediates. Overall, vacancies that are generated by the reaction cover approximately 23 % of the surface, *OH2 species cover approximately 2 %, *OH and *OOH species cover approximately 18 % each, and *O intermediates cover approximately 38 %. The surface coverage percentage fluctuates by approximately 1 % for each intermediate when averaged over 100 or more iterations at times, e. g., iterations 5000–5100 vs. iterations 5100–5200.

Simulations with an activation energy transition probability require a high input voltage to observe some reactions. For example, when the voltage lowers ΔG to zero, as shown in Equation (6), the activation energy changes according to Equations (8) to . At room temperature, the exponent according to Equation (9) is approximately 10−6. Such value, compared to other rates with an exponent closer to unity, will create a rate of about six orders of magnitude smaller than other rates, so the reaction with such an activation energy will rarely happen.

Raising the input voltage to 2 V lowers the barriers even more, so the only reactions with ΔG+λ>0 are the transitions from *OH and *O sites, which have approximately equal probabilities. Because of the *OOH−*OH2 interaction, *O becomes the dominant species on the surface, as shown in Figure 6.

Figure 6

Average simulated distribution of species by number over 10000 frames (top). Spatial distribution of termination species with a simulated voltage of 2 V and pH=13.6 (bottom). Both simulations were run with 30×30 sites for 100,000 steps.

When activation energy is used for transition probability, *O becomes even more dominant than when only thermodynamics is used. *OH has the next highest abundance, covering 13 % of the surface, and other species cover 3 %–5 % each. Applying an even higher voltage, one that reduces the sum of free energy change and reorganization energy below zero, produces the same results as the thermodynamic simulation.

In conclusion, we have shown a new possible explanation for the abundance of *O intermediates during OER catalysis by hematite. Previous experimental and theoretical work associated the dominance of a *O intermediate with a spectroscopic peak that remained apparent during the operation due to comparable surface absorption wavelengths and ab initio calculations for this specific *O intermediate.[ , ] To do so, they used a single‐site model to show that the *O intermediate is dominant on the surface. In this work, we show that the single‐site approach may not be enough for a full description of catalysis since near‐neighbor interactions affect the surface coverage and that the source of the *O dominance is, in fact, the near‐neighbor interactions.

The mechanism of *O generation can be explained by the stability of the intermediates in the model accounting for near‐neighbor interactions. If there is proximity between two active sites holding *OH2 and *OOH interactions, they break down into *OH and *O intermediates and release a water molecule because both intermediates are unstable when they are near each other. The *O to *OOH reaction has a high free energy change, which makes the reverse reaction more favorable. In contrast, the *OH2 to *OH reaction has a negative free energy change, making this reaction favorable. Together, the two neighboring species are inclined to break down into new intermediates. Attractive forces between OH and H create a water molecule from the broken‐down species. The resulting intermediates (*OH and *O) have the highest free energy change for proceeding further in their reaction; thus, catalysis is hindered. The transition of *OH2 and *OOH intermediates to *OH and *O is found to be a reason for the domination of *O on the surface, and this result is in agreement with experimental measurements.[ , ]

With this work, we demonstrated how neighboring sites may use each other's presence as neighbors to reduce the energy required for catalytic steps. To the best of our knowledge, this work is the first to show the effects of neighbor‐site species interactions that do not create a new intermediate. For water splitting with hematite, the interaction of neighbors appears to be an inherent problem, so the thermodynamic limit of the OER might be tighter than previously thought.

Conflict of interest

The authors declare no conflict of interest.

As a service to our authors and readers, this journal provides supporting information supplied by the authors. Such materials are peer reviewed and may be re‐organized for online delivery, but are not copy‐edited or typeset. Technical support issues arising from supporting information (other than missing files) should be addressed to the authors.

Supporting Information

Click here for additional data file.

Supporting Information

Click here for additional data file.