Surface-Specific DFT + U Approach Applied to Î±-Fe2O3(0001) - The

Article pubs.acs.org/JPCC

Surface-Specific DFT + U Approach Applied to α‑Fe2O3(0001) Xu Huang, Sai Kumar Ramadugu, and Sara E. Mason* Department of Chemistry, University of Iowa, Iowa City, Iowa 52242, United States S Supporting Information *

ABSTRACT: We report the bulk properties and ab initio thermodynamics surface free energies for α-Fe2O3(0001) using density functional theory (DFT) with calculated Hubbard U values for chemically distinct surface Fe atoms. There are strong electron correlation effects in hematite that are not welldescribed by standard DFT. A better description can be achieved by using a DFT + U approach in which U represents a Hubbard on-site Coulomb repulsion term. While DFT + U calculations result in improved predictions of the bulk hematite band gap, surface free energies using DFT + U total energies result in surface structure predictions that are at odds with most experimental results. Specifically, DFT + U predictions stabilize a ferryl termination relative to an oxygen termination that is widely reported under a range of experimental conditions. We explore whether treating chemically distinct surface Fe atoms with different U values can lead to improved bulk and surface predictions. We use a linear response technique to derive specific Ud values for all Fe atoms in several slab geometries. We go on to add a Coulomb correction, Up, to better describe the hybridization between the Fe d and oxygen p orbitals, accurately predicting the structural and electronic properties of bulk hematite. Our results show that the sitespecific Ud is a key factor in obtaining theoretical results for surface stability that are congruent with the experimental literature results of α-Fe2O3(0001) surface structure. Finally, we use a model surface reaction to trace how the various DFT + U methods affect the surface electronic structure and heterogeneous reactivity.

■

INTRODUCTION The surfaces of iron oxides are essential components in a wide range of environmental and technological processes such as contaminant adsorption1,2 and electrochemical water splitting.3−6 Fundamentally, the heterogeneous reactivity of iron oxides arises from the coupling of geometry and electronic structure, and therefore a reliable prediction of the surface phase diagram is essential to delineating structure−reactivity relationships. Hematite (α-Fe2O3) is the most stable and abundant iron oxide, and the (0001) plane is reported to dominate under natural conditions.7,8 Bulk α-Fe2O3 has a corundum-type structure, shown in Figure 1. The oxygen layers are parallel and packed along the surface normal (0001). Two Fe layers are between each of the two neighboring oxygen layers. Using R to denote the bulk stoichiometric stacking unit, −(Fe−O3−Fe)−, the three chemically distinct terminations that can result from cleaving the α-Fe2O3(0001) surface are R−Fe−O3−Fe (Feterm), R−Fe−O3 (O-term), and R−Fe−O3−Fe−Fe (Bi-term). A fourth surface structure discussed in the literature is the R− Fe−O3−FeO (ferryl) termination. Side views of these four terminations are shown in Figure 2. Different experimental methods have been applied to determine the surface structure using scanning tunneling microscopy (STM) and low-energy electron diffraction (LEED),7,9−20 the resonant tunneling model (RTM),12,21 crystal truncation rods (CTRs),7,22,23 resonant anomalous X-ray reflectivity (RAXR),24 X-ray photoelectron spectroscopy (XPS),16,22,25,26 and X-ray powder © XXXX American Chemical Society

Figure 1. (a) Rhombohedral and hexagonal bulk structure of the αFe2O3 unit cell. (b) Side view of the bulk α-Fe2O3 structure with the alternating spin states indicated by arrows along the c-axis.

diffraction (XRD)19,22 to investigate stable surface domains under different humidity, oxygen partial pressure, hydroxylation, and pH conditions. While the general consensus is that the Fe-term and O-term geometries (or their hydroxylated forms) are the dominant surface structures over the range of experimental conditions, there are also reports supporting a ferryl termination. Interestingly, the STM study of Lemire and co-workers15 suggests that ferryl groups may coexist with domains of the Fe-term surface, while Barbier and co-workers22 Received: December 11, 2015 Revised: February 12, 2016

A

DOI: 10.1021/acs.jpcc.5b12144 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

band gap.32,33,38 Convinced by the improvements seen for bulk hematite properties, hematite surface studies have also adapted DFT + U styles. For example, DFT + U studies have been carried out to study water splitting on hematite, showing that the Fe termination exhibits a relatively small energetic barrier to heterolytic water dissociation.39 However, when DFT + U calculations are used in conjunction with ab initio thermodynamics predictions of surface stability, previous studies report that the Fe-term surface is stable over most of the accessible range of oxygen chemical potential,8 with a small window of stability for the ferryl surface.33 Therefore, despite improved prediction of certain bulk properties, DFT + U methods do not appear to produce a reasonable α-Fe2O3(0001) phase diagram, and this puzzling result has been noted in the literature.26 In evaluating the performance of DFT + U calculations, the choice of the applied U value must be considered. One option is to empirically select a U value that is able to match experimental results for some physical property such as the lattice constant, magnetic moment, or the band gap.32,40−44 It has also been suggested to select a U value based on energetics between oxidized and reduced oxide forms.45 Alternatively, values of U can be derived from first-principles,46 and this approach has been applied to systems such as CeOx where the calculated U values for Ce2O3 bulk, for CeO2 bulk, and for the Ce(111) surface are reported to be 6.70, 5.09, and 5.20 eV, respectively. For Ce in cluster complexes such Ce3H3,6,7O7, the calculated U values vary from 2.1 to 2.4 eV.47 Several theoretical studies have addressed the structural dependence of U. Kulik and co-workers studied the potential energy surface of FeO+, suggesting that U is a function of Fe−O distance, denoted as DFT + U(R). This work went on to develop a self-consistent iteration process to calculate U.48,49 Using the self-consistent process, Hsu and co-workers studied low-spin LaCoO3 and found that the U derived from the iteration produces better agreement with experimental data.50 Another consideration in DFT + U applied to transition metal oxides is how to treat the oxygen p orbitals. Strong d- and p-orbital hybridization in transition metal oxides is supported by spectroscopic data.51,52 When the U correction is added to only the cation d orbital of an oxide it can disrupt the interaction with oxygen p states, as noted, for example, by Guo and co-workers.53 It has been put forward that it is also necessary to apply a Hubbard U on O p orbitals, denoted Up, to treat all the valence states with a more balanced Coulomb correction. Studies using Up on different transition metal oxide systems in the literature include bulk TiO2 with defects54,55 and bulk LiCoO2 and Co(II)(H2O)6 complex56 where lattice constants and band structures were improved. Also with the Up treatment, the observed spin state of a Ni-based magnet was confirmed as the ground state,57 and for a defect ZnO system,58 substitutional nitrogen (NO) was revealed to be a deep acceptor. In the present study, we aim to resolve the DFT + U methodology dependence in predicting the physical properties and surface stability of α-Fe2O3(0001). We systematically test and compare the bulk and surface properties predicted using standard DFT-GGA and various DFT + U approaches. Part of our approach is to derive surface-specific Ud values for the chemically distinct Fe atoms of the surface models.

Figure 2. Side views of the DFT-GGA-optimized geometries of the αFe2O3(0001) surface models included in this study.

observe the ferryl termination over a wide range of conditions but do not see the Fe-term structure in their experiments. By using an ab initio thermodynamics framework developed by Reuter and co-workers,27−29 density functional theory (DFT) calculations can be used to predict surface-phase diagrams as a function of oxygen chemical potential, μO. The αFe2O3(0001) phase diagram predicted by Wang and co-workers supports that the Fe-term and O-term surfaces are preferred for the (1 × 1) surface at low and high values of μO, respectively.20 In later work by Bergermayer and Schweiger, a ferryl termination was additionally considered and reported to have similar stability to the Fe-term surface.30 When the presence of H2O is considered, hydroxylated forms of the Fe-term and Oterm structures are predicted.7 However, all of the aforementioned theoretical predictions were based on DFT calculations using local exchange-correlation functionals. Owing to the 3d Fe electrons, hematite is a strongly correlated system in which the Fe d- and O p-orbitals hybridize significantly. DFT has a strong tendency to overdelocalize d electrons, resulting in the underestimation of the band gap in metal oxides.31 Specifically for hematite, DFT predicts only ≈0.3 eV band gap,32,33 while the observed band gap is 2.2 eV.34 The DFT + U method was developed35−37 to compensate for the on-site repulsion between d and f electrons. When applied to bulk hematite, DFT + U improves the prediction of certain bulk properties such as the magnetic moment and the

■

THEORETICAL METHODS Computational Details and Atomistic Thermodynamics. Spin-polarized DFT59,60 calculations were carried out using B


Article


Table 1. Comparison of the Calculated Bulk Properties of DFT-GGA and Different DFT + U Methods with Experimental Valuesa a experiment DFT-GGA Ud Ud+p

b

5.035 5.046(+0.22%) 5.107(+1.4%) 5.076(+0.81%)

c

gap character

b

c

13.747 13.809(+0.45%) 13.959(+1.5%) 13.822(+0.55%)

p−d d−d p−d p−d

μB d

4.0 3.15 3.44 3.67

B0 2.25−2.58c,e,f 179.5 176.3 207.2

Values in parentheses are percent error relative to experimental reference values. a and c refer to the α-Fe2O3 lattice constants, reported in Å. Gap character refers to the orbital character of the valence and conduction bands, with p referring to oxygen 2p and d referring to Fe 3d. Magnetic moments are reported in Bohr magneton units, μB, and bulk moduli, B0, are reported in units of GPa. bReference 76. cReference 72. dReference 73. e Reference 74. fReference 75. a

the QUANTUM-ESPRESSO (QE) package61 with the generalized gradient approximation (GGA) of Perdew, Burke, and Ernzerhof (PBE)62 to the exchange-correlation functional. Ultrasoft pseudopotentials63 were used with a plane-wave cutoff of 35 Ry for the wave function and 280 Ry for the charge density with a self-consistent field convergence criterion of 0.001 eV. The theoretical and experimental lattice constants of bulk α-Fe2O3, as well as percent error, can be found in Table 1. A 6 × 6 × 6 Monkhorst−Pack k-point mesh64 was used for bulk calculations, which was reduced to a 6 × 6 × 1 mesh for calculations of (0001) surfaces modeled in a (1 × 1) cell. Surfaces were modeled as symmetric slabs (with two equivalent exposed surfaces) with at least six oxygen atomic layers and 16 Å of vacuum space along the surface normal. All structures are optimized without constraint, with a force tolerance of 0.01 eV/ Å. The surface free energy, γ, is calculated within the firstprinciples thermodynamics framework of Reuter and coworkers27−29 and in the style detailed by Lo and co-workers for hematite.65 Briefly, at p = 0 and T = 0, the surface free energy is defined as 1 (Eslab − E bulk ) 2A ⎛3 ⎞ ⎤ 1 ⎡ 1 = ⎢⎣Eslab − NFeμFe 2 O3 + ⎝⎜ NFe − NO⎠⎟μO⎥⎦ 2 2A 2

Surface Models. Figure 1(a) shows the side view of a unit cell for bulk hematite, where a and b are equal and form an angle of 120° and c is along the (0001) direction perpendicular to the a−b plane. Figure 1(b) shows an expanded side view of the bulk structure. The spin states of Fe are shown in Figure 1(b), with arrows used to indicate the antiparallel spins of neighboring Fe bilayers. The overall spin state is antiferromagnetic (AF). We model the four surface terminations, Fe-term, O-term, Bi-term, and ferryl, as defined in the Introduction and shown in Figure 2. Additionally, two hydroxylated surfaces, Fe−OH-term (a hydroxylated form of the Fe-term surface) and OH-term (a hydroxylated form of the O-term surface) are included. Details of the hydroxylated surfaces have been recently reported.66 Side views of all of the optimized surface geometries are shown in Figure 2, with the common bulk structure shaded in gray. Figure 2 also defines atomic layer numbering that we use in discussing surface relaxations and derived Hubbard Ud results. The numbering scheme includes all distinct layers in the symmetric slab models. As reviewed in the Introduction, the prevalence of Fe-term and O-term surfaces (or their hydroxylated forms) is supported by a wealth of experimental studies, while domains of the ferryl termination may exist at high oxygen chemical potential conditions. The Bi-term surface is not reported experimentally and is included here as a theoretical point of comparison. Ud Calculations. Hubbard Ud values were derived using Cococcioni and de Gironcoli’s linear response theory46 method as implemented in QE, taking into account the reduced dimensionality of the slab models. A review of the linear response approach, as well as specific considerations made to use this method on slab geometries, is provided in the Supporting Information (SI). For the Fe Ud value in bulk αFe2O3, we use the hexagonal cell in Figure 1(a) with 30 atoms (12 Fe and 18 O). We obtained Ud = 3.81 eV from perturbation of Fe-d and O-p orbitals. The value of the perturbation parameter used in the linear response theory, α, was 0.06, as discussed in the SI and shown in Figure S.1. On the basis of results using the extrapolation scheme of Cococcioni and de Gironcoli, the value of Ud = 3.81 eV is converged within ±0.01 eV with respect to cell size. Further discussion of the variation in the calculated Ud values as a function of system size is presented in the SI. The reduced symmetry of the surface slabs leads to more chemically distinct Fe and O pairs. From the structural models in Figure 2, the terminal Fe and O atoms are obviously chemically distinct from the inner bulk-like species. There are also subtle differences in the layers of intermediate depth due to surface relaxations. For the Fe-term structure in Figure 2, for example, there are six types of Fe: Layer1, 3, 4, 6, 7, and 9and four types of oxygen: Layer

γ=

(1)

where Eslab and Ebulk are the total energies of the slab and bulk system with the same number of atoms as in the slab, and A is the surface area of the unit cell. Here the Gibbs free energy, G = Eel,slab + Evib + Eother,internal + PV − TS, is approximated by the dominant term, the electronic energy Eel, to give G ≅ Eel,slab.30 Similarly for the bulk system, there is Eel,bulk ≅ G = NFeμFe + NOμO. To connect the chemical potential of Fe, O, and bulk Fe2O3 and to set up the poor limit of μO, the equilibrium with the bulk is considered as μFebulkO = 2μFe + 3μO 2 3

(2)

which represents the limit as the decomposition of the material when μO is too low. The rich limit of μO is considered as a gasphase equilbrium 1 1 EO2 ≅ μO = μO 2 2 2

(3)

where EO2 is the calculated total energy of the O2 molecule. In addition to eq 2 for the lower limit of μO, we also consider the limit in terms of reducing hematite. Specifically, we also consider the limit at which α-Fe2O3 reduces to Fe3O4 or Fe1−xO,10,11,18,19 as will be discussed below. C


Article


Table 2. Percent Relaxations (Relative to the Corresponding Theoretical Bulk Spacings) for the Interlayer Spacings in αFe2O3(0001) Surface Modelsa layers 0−1 1−2 2−3 3−4 4−5 5−6 6−7 7−8 8−9 9−10

Fe-term Fe/O−Fe Fe−O O−Fe Fe−Fe Fe−O O−Fe Fe−Fe Fe−O O−Fe Fe−Fe

−57.91 6.95 −30.80 14.33 4.10 −2.45 0.29 −0.91 3.45

(−61.48) (6.95) (−40.05) (15.91) (6.04) (−4.66) (0.23) (−1.28) (3.51)

ferryl 6.92 −17.01 0.89 −32.60 15.67 −0.59 1.40 −0.46 −0.89 2.02

(72.11) (−60.63) (9.26) (−42.72) (16.79) (5.85) (−3.02) (−0.13) (−0.91) (2.97)

O-term

−4.16 −78.32 35.61 −5.98 14.96 −4.30 −0.96 0.89

Bi-term

(1.09) (−15.87) (5.31) (−1.59) (3.44) (2.83) (−3.41) (0.40)

19.64 −64.14 29.26 −1.12 −11.20 9.55 −3.67 −1.51 2.75 −2.71

(81.44) (−89.95) (24.04) (1.74) (−6.15) (6.87) (−12.10) (2.95) (−0.15) (3.04)

OH-Fe-term 12.48 15.48 −38.15 14.71 −0.22 5.54 −2.06 −0.31 1.15

OH-term

(9.64) (13.03) (−41.43) (15.76) (−0.78) (4.53) (−1.48) (−0.68) (0.71)

12.03 −24.50 8.78 0.87 3.22 −1.25 0.21 0.08

(5.73) (−26.64) (10.31) (−0.37) (5.04) (−0.24) (0.10) (−3.27)

a

Values are given for DFT-GGA calculations, with results for DFT + specific Ud in parentheses. The atomic layer numbering follows the scheme depicted in Figure 2.

Table 3. Surface-Specific Derived Ud Values for Fe Atoms (By Layer) in Different Surface Modelsa

2, 5, and 8. Therefore, in adapting the linear response method for determining Ud values in surface geometries, perturbations for each chemically distinct Fe and O pair are necessary to form the response matrices χIJ and χ0IJ (as defined in the SI). The surface-specific Fe Ud values are determined as the diagonal elements of the final U matrix. The perturbation parameter, α, was tested with different values from 0.01 to 0.15. The optimal one giving consistent and stable Ud values for each surface is chosen to be 0.06. Plots of calculated Ud values as a function of α are presented in Figure S.1. As noted in the Introduction, there are several reports in the literature that note the need to also use Up when modeling transition metal oxides. Previous studies have carefully considered the choice of Up value. It is possible to obtain Up values simultaneously to Ud values in the linear response approach. Mattioli and co-workers,54 in studying bulk TiO2, obtained a Up value of 10.65 eV from QE linear-response calculations. They compared the physical properties of the oxide using the derived value of Up to those obtained using Up = 5.90 eV, a value obtained from X-ray spectra.67 (There are also Auger spectra data used for obtaining Up values for copper oxides, which suggest Up values ranging from 4 to 7 eV.68) Mattioli et al. found that calculations using Up = 10.65 eV made the Ti−O bond too strong and therefore overestimated the bulk modulus and formation enthalpy. They concluded that calculations using a Up value of 5.90 eV gave better results than either calculations using no Up or the derived Up value. In other theoretical works, a range of Up values between 5.0 and 6.0 eV was found to achieve desirable accuracy in the prediction of lattice constants and band structure56,58 and the correct nature of band gap and the spin state of the ground state.57 For other transition metal oxides, bulk Up values can be obtained theoretically through Koopmans theorem in Hartree−Fock theory, and the resulting Up value for TiO2 is 5.2555 or 4.3 eV for ZnO.69 In the present study, we are interested in the variation of U as a function of varying chemical environments created by cleaving surfaces. We expect variations in Ud as a function of chemical environment due to the localized nature of the Fe d-band. (See the later discussion about Tables 2 and 3 and Figure S.3 in SI.) As oxygen only has electrons in delocalized sp states, changes in the bonding environment are not expected to affect the electronic structure as drastically as in iron. On the basis of the detailed previous investigations of Up and the expectation that iron will be more susceptible to surface-specific effects, we use the literature value of Up = 5.90

Fe layer

Fe-term

ferryl

0 1 3 4 6 7 9

5.90 4.20 4.26 4.30 4.22 4.38

8.52 4.42 4.26 4.27 4.19 4.23

O-term

Bi-term

OH-Fe-term

OH-term

5.39 5.56 4.17 4.00 4.04

2.80 3.54 3.87 3.66 3.80 3.89 3.74

4.63 4.48 4.25 4.24 4.22 4.22

4.50 4.32 4.15 4.19 4.17

a

All values are reported in eV. The Fe atomic layer numbering follows the scheme depicted in Figure 2.

eV for oxygen along with derived Fe Ud values and denote this method as “DFT + Ud+p.”

■

RESULTS AND DISCUSSION Bulk α-Fe2O3 Properties. Bulk properties calculated using different DFT methods are summarized in Table 1, with comparison to experimental data. The band structure of bulk αFe2O3 is presented as projected density of states (PDOS) plots shown in Figure 3. For standard DFT-GGA calculations, the theoretical lattice constants are in good agreement with experiment, with percent errors of +0.22%. The magnetic moment of Fe and the bulk modulus are significantly underestimated (by at least 20% for both). In terms of electronic structure, there is almost no band gap, consistent with previous DFT-GGA results for hematite.32,33,66,70 The PDOS also reveals the orbital character of the band gap and shows that the valence band maximum and conduction band minimum are both comprised of the Fe d-orbital, contradicting spectroscopic descriptions and calculations71,72 that assign the band gap character as O 2p-Fe 3d. For DFT + U calculations, our derived bulk Ud value is 3.81 eV (as described in Part C in Theoretical Methods), which is similar to the U values used in other theoretical studies.8,32,33,38,70 In Table 1, the calculated lattice constants exhibit 1.5% error from experimental values. The calculated bulk modulus was worsened in DFT + U relative to the standard DFT-GGA calculations, which can be interpreted to mean that adding Ud softens the structure. The calculated magnetic moment value was improved relative to standard DFT-GGA, and this is attributed to the localization Fe d-orbital and increased on-site Coulomb repulsion. The band gap in the D


Article


Figure 4. Change in 3-D charge density of using different U methods displaying in the rhombohedral α-Fe2O3 cell. Light (yellow) isosurfaces show regions of electron density gain, and dark (purple) isosurfaces show regions of electron density loss. (a) Charge density difference between DFT + Ud and DFT-GGA. (b) Charge density difference between DFT + Ud+p and DFT + Ud.

Fe, the electron density becomes more localized along the Fe− O bonding direction. Meanwhile, electron density around O is displaced away by the repulsion along the bonding direction from Fe, to occupy the nonbonding region. Therefore, the bonding between Fe and O is weakened, resulting in elongation of the Fe−O bond (corresponding to the expanded lattice constants), and the structure was softened (corresponding to the smaller bulk modulus), congruent with the calculated values in Table 1. Figure 4(b) shows the charge density difference between the DFT + Ud/Ud+p calculations, again with gain (loss) shown in yellow (purple). The shape and positions of the isosurfaces indicate that electron density in bonding states on oxygen increases, indicative of the improved overlap between Fe d- and O p-orbitals. This charge rearrangement is in line with the calculated physical properties presented in Table 1: The Fe−O bond length increases (resulting in smaller lattice constant values), and the structure was hardened (resulting in a larger bulk modulus value). From the analysis of the charge density and calculated properties of bulk hematite, we conclude that the Ud+p calculation balances the Coulomb repulsion between d- and p-orbitals to give an improved description of the band structure of hematite. On the basis of these results, we extend the Ud+p calculation approach to our surface calculations as well. Surface Structure and Surface-Specific Ud Values. The structural details of the optimized surface models are summarized in terms of the percent relaxations of each vertical layer relative to bulk positions, reported in Table 2. Both data for standard DFT-GGA and optimizations using DFT + Ud+p with our derived surface-specific Ud values and a constant Up = 5.90 eV are given. Results are given for all surface layers down to the inversion plane of the slab models, using the layer numbering defined in Figure 2. When no U is applied, for the Fe-term, ferryl, and Bi-term surfaces, there are significant contractions in the near surface layers, notably in the 1(Fe)−

Figure 3. PDOS of bulk α-Fe2O3 using different methods. The Fermi level is referenced to zero in all plots. (a) DFT-GGA. (b) DFT + Ud = 3.81 eV. (c) DFT + Ud+p, where Ud = 3.81 eV and Up = 5.90 eV.

PDOS is about 1 eV (Figure 3) and exhibits the correct nature of the band gap, indicating a p−d transition. All of our calculated results for bulk α-Fe2O3 using DFT and DFT + Ud are consistent with previous results32,33and show that some (but not all) physical property predictions are improved with the addition of Ud. However, when Up = 5.90 eV is also applied, the agreement with experiment for the physical properties unilaterally improves: The error in calculated lattice constants is less than 1%, and the magnetic moment, 3.67 μB, and the bulk modulus, 207.2 GPa, are also closer to experimental values of 4.0 μB73 and 225−258 GPa,72,74,75 respectively, as summarized in Table 1. The band gap also approaches the experimental value of 2.2 eV.76 Compared to DFT + Ud results, the DFT + Ud+p PDOS shows improved overlap between Fe d- and O p-orbitals in the valence band, especially in the energy region between 6 and 8 eV below the Fermi level. To visualize the effects of successively adding Ud and Up, three-dimensional plots of charge density differences are produced for the bulk, as shown in Figure 4. In Figure 4(a), the difference between the DFT + Ud and DFT charge density is shown, with gain (loss) shown in yellow (purple). Around E


Article


relatively small changes in geometry between DFT-GGA and DFT + Ud+p methods. There is about 10% and 3% further contraction in 3(Fe)−4(Fe) in Fe-term and OH-Fe-term compared with their no U results. For the OH−O−term surface, the 3(Fe)−4(Fe) contraction only increased by 2%. The ferryl and Bi-term surfaces show greater structural variation with Ud+p, where the interlayer contraction and expansion were both increased significantly. The FeO distance in 0(O)− 1(Fe) in ferryl increased by 72%, and the 0(Fe)−1(Fe) in Biterm increased by 81%. The 1(Fe)−2(O) distance contraction was also enhanced with 60% in ferryl and 90% in the Bi-term surface. The 3(Fe)−4(Fe) distance in ferryl also increased by 10%, while in the Bi-term surface it almost remained the same. The O-term is the only case that shows much less contraction and expansion in the interlayer distance compared with bulk when Ud+p is applied. The significant ≈80% contraction in 3(Fe)−4(Fe) distance now has reduced to only ≈16%, and the 35% 4(Fe)−5(O) expansion also reduced to only 5%. The calculated specific Ud values for distinct Fe atoms in the surface slab models are summarized in Table 3. The largest deviations from the bulk-derived Fe Ud value of 3.81 eV are seen in the layer 1(Fe) atoms in the Fe-term, ferryl, and O-term surfaces. The Ud values of the layer 1(Fe) atoms in the hydroxylated surfaces show less variation from the bulk-derived value. The Bi-term surface is the only incident of layer 1(Fe) Ud values that are smaller than the bulk value, with the topmost Fe Ud value of 2.8 eV of Ud. Compared with the Ud value derived from bulk α-Fe2O3 in this study, 3.81 eV, the calculated specific Ud values of the innermost Fe atoms in the slab models, namely, layer 6, 7, and 9 in Table 3, are generally 0.2−0.4 eV larger. The difference between the derived Ud values for the bulk and bulk-like slab Fe atoms can be attributed to the employed linear-response Udderivation process. As the linear-response approach depends on making small perturbations to atomic positions, it inherently depends on the cell size. However, going to larger and larger supercells becomes computationally prohibitive. As a compromise, an extrapolation scheme is used.46 However, owing to the reduced dimensionality of the surface slabs, the extrapolation along the surface normal direction is not possible. We further tested the effect of perturbation cell size by using Fe-term and O-term slabs of five, six, and seven stoichiometric units, as reported in detail in the SI and summarized in Figure S.2. For all three slab thicknesses, the innermost Fe atoms of the slab models exhibit derived Ud values that are ≈0.2−0.4 eV larger than the bulk-derived value. The largest variation in Ud value for any specific Fe atom as a function of slab layer thickness is on the order of 0.2 eV. As determined in a previous study by Rollmann and co-workers,32 the physical properties of hematite are not significantly affected by variations of this magnitude in the value of Ud. For example, over the range of Ud = 3.5−4.5 eV, the changes in bulk hematite volume, magnetic moment, and band gap are 0.025 Å3, 0.1 muB, and 0.3 eV, respectively. The derived specific Ud values therefore are determined to be robust within the dependence of calculated physical properties. We also consider the fact that comparing DFT + U total energies in which different values of U are applied can be problematic.77 We note that the calculation of γ using eq 1 depends on the difference between the surface and bulk Gibbs free energies, which allows for cancellation of error. The results shown in Tables 2 and 3 exhibit a correlation between surface layer spacings and calculated Ud values. In general, the outermost Fe layers that undergo greater vertical

Figure 5. PDOS comparison of the d-orbital of the layer 1(Fe) and the p-orbital of the topmost layer 0(O) in the ferryl surface. The Fermi level is referenced to zero in all plots. (a) DFT-GGA, (b) DFT + Ud, with Ud = 3.81 eV, (c) DFT + Ud using surface-specific values of Ud, and (d) DFT + Ud+p using surface-specific values of Ud and Up = 5.90 eV.

2(O) interlayer distance. For Fe-term and Bi-term surfaces, this contraction can be up to ≈60 and 75%, respectively. The 0(Fe)−1(Fe) distance on the Bi-term surface expands by up to ≈20%. For the hydroxylated surfaces, the OH-Fe-term does not exhibit large surface relaxations because the terminal Fe atom forms bonds with H2O and OH groups. The bulk O3−Fe−Fe− O3 units in general show ≈30% contraction between the 3(Fe)−4(Fe) layers for all the surfaces except the Bi-term structure where it was almost unchanged. For the O-term surface, this contraction is up to ≈80%. Meanwhile, the 2(O)− 3(Fe) and 4(Fe)−5(O) interlayer distance generally increased, except in the O-term surface in which the 2(O) layers contract inward by about 4%. Comparing with the O-term structure, the hydroxylated OH-O-term exhibits surface relaxations of smaller magnitudes, reflecting the greater extent of bond saturation of the hydroxylated oxygen atoms. Going deeper into the slabs, the interlayer spacings in all of the surfaces gradually return to bulk-like geometry, supporting that the slab thickness is sufficient for surface modeling. It has been noted in the literature that because U depends strongly on chemical environment it may be necessary to reoptimize the system geometry iteratively until convergence between U and the structure is achieved.47,48,50,77,78 To this end, in this work, we reoptimize the lattice constant at each level of theory (with the results detailed in Table 1), and the slab models are all reoptimized without constraint. Referring to Table 2, the Fe-term, OH-Fe-term, and OH-term surfaces show F


Article


Figure 6. Free energies of various surface models for α-Fe2O3 (0001). (a) DFT-GGA, (b) DFT + Ud with Ud = 3.81 eV, (c) DFT + Ud using surfacespecific values of Ud, and (d) DFT + Ud+p using surface-specific values of Ud and Up = 5.90 eV.

Ud is applied only to Fe, then the lower limit of μO is higher due to the increase in the total energies of the bulk α-Fe2O3 in eq 2, while the upper limit of μO is unaffected. To make the total energy results consistent, we use a derived Ud for Fe bulk metal (2.3 eV) to draw back the lower limit of μO in Figure 6(b) and (c). The inclusion of Ud on Fe bulk is motivated by the known error associated with mixing total energies from DFT-GGA and DFT-GGA + Ud calculations.77,79 When Up is applied to the Fe2O3, we also applied Up = 5.90 eV to determine the O2 total energy used in eq 3, which has the effect of shifting the upper limit of μO higher. Treating the O2 molecule with Up not only helps the cancellation of errors but also addresses the wellknown overbinding problem of GGA in the O2 molecule.62,80 As discussed previously,81 the DFT-GGA error in O2 does not cancel out in the energy sums and differences used in the ab initio thermodynamics framework. In the present work, we note that when Up is applied to the O2 the atomization energy error is reduced from 20% to 6%. Similarly, the error in the Fe2O3 formation energy is reduced from 28% in the DFT-GGA calculation to 7% in the DFT + Ud+p calculation. We therefore conclude that by using a Ud+p approach the O2 error largely cancels out without any further correction to the O2 total energy. We first discuss the phase diagram results for the surface models that do not include hydroxylation, namely, the Fe-term, O-term, Bi-term, and ferryl surfaces. Consistent with previous DFT-GGA results,20,30,33 the phase diagram in Figure 6(a) shows that as a function of increasing μO the Fe-term surface, ferryl surface, and O-term surface are favorable. As reviewed in the Introduction, the predicted stability of the ferryl surface is not broadly supported by experimental studies. The DFT + Ud phase diagram shown in Figure6(b) shows both the ferryl and O-term surfaces destabilized relative to the Fe-term surface, with the ferryl surface now lower in surface free energy than the O-term surface. This result does not provide an improved level of agreement with experimental observations. The phase diagram resulting from surface-specific Ud values shown in

contraction exhibit larger Ud values than those that do not deviate as much from the bulk structure. Likewise, terminal Fe atoms that exhibit vertical expansion have smaller Ud values. The correlation between surface structure and derived Ud values can be rationalized as follows: Surface layer contraction leads to shorter Fe−O bond lengths and greater Fe d- and O porbital hybridization, while expansion reduces interaction. Shorter Fe−O distances/greater hybridization between Fe and O indicates that the energy of the Fe d-bands will be more sensitive to changes in occupation and thus will result in larger Ud values in the linear response method.46 On the contrary, weaker Fe−O interaction (such as what is experienced in the terminal Fe layers in Bi-term surface) corresponds to a more metallic Fe d-band and thus a smaller Ud value. A similar trend had been revealed in previous work suggesting the strong correlation between U values and Fe−O bond length.48 The inverse relationship between surface relaxation and derived Ud values is demonstrated for the Fe-term and O-term surfaces in Figure S.3. However, the number of Fe−O bonds has little effect on the Ud values. For example, in the Fe-term surface, the Fe atoms in layer 1 have lost half of the Fe−O bonds relative to bulk Fe. However, the loss of bonds is compensated for by surface relaxations. Therefore, no simple prediction of the surface-specific Ud value can be made a priori based on atom coordination. Surface-Phase Diagrams and Heterogeneous Reactivity. The the various DFT methods are used to calculate γ, the results of which are presented as five α-Fe2O3(0001) surface-phase diagrams in Figure 6. The methods and their corresponding phase diagrams are DFT-GGA in Figure 6(a), GGA + Ud with the bulk-derived Ud value of 3.81 eV in Figure 6(b), GGA + Ud with surface-specific values of Ud from in Table 3(c), GGA + Ud+p using the bulk-derived Ud value and a constant Up value of 5.90 eV in Figure 6(d), and GGA + Ud+p using surface-specific values of Ud and a constant Up value of 5.90 eV in Figure 6(e). Immediately obvious from the plots in Figure 6 is the fact that the range of accessible μO varies. When G


Article


Figure 7. Summary of theoretical and experimental surfaces predicted or observed as a function of μO. (a) and (b) are the most stable surfaces extracted from Figure 6 using DFT-GGA and DFT + Ud+p with derived values of Ud and Up = 5.90 eV, respectively. (c) to (e) are the experimental observations using different sample preparation methods as detailed in the figure. All of the x-axes are on the same scale. Black ranges in (a) and (b) and vertical bars in (c)−(e) indicate the Fe-term surface; red ones indicate the ferryl surface; and the blue indicates the O-term surface. The bars with half black and blue (D in (d); G, H, and J in (e)) indicate that both the Fe-term and O-term surfaces were observed, and the bar with black and red (E in (d)) indicates the presence of both Fe-term and ferryl surface was reported. Some surfaces were seen upon a wide range of μO,rel. They are shown in overlapping horizontal lines (F in (d) and I in e). The assignments of A−J to specific experimental references are as follows: A,9 B,25 C,17 D,19 E,15 F,22 G,20H,16 I,18 and J.14

Here we use the value from the JANAF thermochemical tables.82

Figure 6(c) also fails to produce results consistent with experimental reports, though the relative stabilities of the surfaces are reordered compared to the previous two methods. The phase diagram obtained using the bulk Ud value and Up = 5.90 eV in Figure 6(d) still suggests that the ferryl surface is more stable than the O-term surface. Finally, when both the surface-specific Ud values and Up = 5.90 eV are applied, for the resulting phase diagram (Figure 6(e)), the results show the ferryl surface to be less stable than the O-term surface, with the latter being preferred relative to the Fe-term surface at the upper limit of accessible μO. To further assess the calculated phase diagram, we consider the physical limits of μO. As can be seen in Figure 6(e), the Biterm surface is favorable at the low end of μO yet has never been reported by experiment. The oxygen-poor limit was defined as the starting point of metal oxide decomposition into bulk iron and molecular O2. However, an alternative lower limit can be determined by considering that α-Fe2O3 could first be reduced into Fe3O4 or FeO1−x. If we take into account the oxidation/reduction equilibrium between Fe3O4 in eq 4 below, then the oxygen-poor limit would be increased from −3.07 to −2.41 eV Fe2O3 ↔

2 1 Fe3O4 + O2 3 6

μO(T , pO ) = 2

⎛ p ⎞⎤ 1 ⎡ total ⎢EO + μÕ (T , p0 ) + kBT ln⎜ O02 ⎟⎥ 2 2 ⎢⎣ 2 ⎝ p ⎠⎥⎦ (5)

Finally, the relative oxygen chemical potential μrel,O, defined in eq 6 is determined using the experimental (T,p) conditions to use in comparison with the theoretical predictions μO,rel = μO(T , pO ) − 2

=

⎛ pO ⎞ 1 kBT ln⎜ 02 ⎟ 2 ⎝p ⎠

1 total [EO2 + μÕ (T , p0 )] 2 2

(6)

The calculated stable ranges of the surfaces at T = 0 K from Figure 6(a) and (e) were extracted in Figure 7(a) and (b), and the experimental observed surfaces were also sorted and converted in terms of μrel,O, as summarized in Figure 7. Figure 7 provides a convenient comparison of the DFT-GGA predictions (a) and the DFT + Ud+p results using the surfacespecific Ud values and the constant Up of 5.90 eV. When U is not used, the value of μO at which the O-term and Fe-term surface free energies cross is ≈−0.5 eV, while for the DFT + Ud+p results the crossover is at μO ≈ −0.2 eV. The DFT-GGA results predict the ferryl termination to be the most stable across μO values of −1.2 to 0.51 eV, while for the DFT + Ud+p results the ferryl surface is never predicted to be the most stable. Figure 7(c)−(e) summarizes the experimental surface studies of α-Fe2O3(0001) in terms of surface structure and μrel,O. The experimental results shown in Figure 7(c)−(e) show some significant disagreement. Figure 7(c) shows results from two studies using oxygen-plasma-assisted molecular beam epitaxy

(4)

Therefore, the theoretical ranges of μO over which the Bi-term surface is predicted to be the most stable are largely nonphysical. To aid the interpretation of the chemical potential axis and to quantitatively compare the stabilities in the phase diagram with experimental results, we translate a range of experimental conditions under which Fe2O3 surface structure was studied using the ideal gas approximation, as shown in eq 5, where the μ̃ O2(T,p0) indicates the contributions from vibrations and rotations of the O2 molecules and the gas entropy at p0 = 1 atm. H


Article


Table 4. Calculated Values of Ef Based on Equation 7 Using Different DFT Methods, Reported in eVa

(OPA-MBE), indicated as A9 and B,25 both with the result that only the Fe-term surface was observed. Figure 7(d) includes the results of four studies using surfaces prepared from natural crystal: C,17 D,19 E,15 and F.22 The results in Figure 7(d) for C−E are in reasonable agreement and either report only the Fe-term surface or both the Fe-term surface and the O-term surface (D) or ferryl sites existing in Fe-term domains (E). Figure 7(e) summarizes experimental results for samples grown on Pt(111) or Mo(111) substrate, G,20 H,16 I,18 and J,14 and all of these studies report either the Fe-term surface exclusively or with the O-term surface. A notable outlier among the experimental results can be seen in Figure 7(d)F, which reports the O-term and ferryl terminations over a range of μO, without ever observing the Fe-term surface. This is particularly interesting as the only other summarized experiment reporting the ferryl surface is that of Lemire and co-workers15 (Figure 7(d)E), which specifically notes that ferryl groups are seen in domains of the Fe-term surface. The differences in surface samples and surface preparations are all likely to contribute to the differences in observed surface structures. However, by comparing a relatively large number of experimental results to the various theoretical phase diagrams produced here, we conclude that the phase diagram obtained using surface-specific Ud values and Up = 5.90 eV shows the best qualitative agreement with the majority of experiments. Going on to consider the hydroxylated surfaces, the DFTGGA phase diagram in Figure 6(a) is in good agreement with previous results using comparable methods,7 with the Fe−OHterm and OH-term surfaces having the minimum surface free energies over most of the accessible range of μO. As discussed by Trainor and co-workers, crystal truncation rod diffraction results support that hydrated α-Fe2O3(0001) is dominated by domains consistent with our OH-term and Fe−OH-term models. A distinction between the phase diagram obtained using surface-specific Ud values and Up = 5.90 eV is applied (Figure 6(e)) compared to all of the other methods is that in the former the OH-term surface is predicted to be more stable than the Fe−OH-term surface. To study how the various DFT + U methods influence predictions of surface reactivity, we analyze the formation of the ferryl surface using a model heterogeneous reaction and electronic structure. The formation energy Ef of the ferryl surface from the Fe-term surface is calculated based on the oxidation reaction Ef =

1 1 ΔEtot = [Eferryl − (E Fe ‐ term + EO2)] 2 2

DFT-GGA DFT + Ud, bulk U(Fe) DFT + Ud+p, bulk U(Fe) DFT + Ud+p, specific U(Fe) a

Ef

FeO

−0.633 0.600 1.033 1.508

1.572 1.606 1.721 2.054

Also reported are the ferryl FO, in Å.

endothermic result is that obtained from the calculations using the surface-specific Ud values and Up = 5.90 eV, and the trends in Ef correlate directly with the FeO double bond length. The PDOS of the p-orbital of the doubly bound layer 0(O) and the d-orbital of the top layer 1(Fe) in the ferryl surface are plotted and presented in Figure 5. Figure 5(a) shows the PDOS of the ferryl surface modeled at the DFT-GGA level and shows that the electron density of Fe d- and O p-orbitals spans the conduction band, indicating strong hybridization that accounts for the favorable value of Ef and relatively short FeO distance. Figure 5(b) shows the PDOS obtained in calculations using bulk-derived Ud. The majority and minority Fe d-bands exhibit localization relative to the DFT-GGA results, indicating less interaction between Fe and O, but there is still an obvious bonding and antibonding overlap pattern at −8 to −6 eV and above the Fermi level, respectively. This is consistent with the slightly endothermic value of Ef = 0.600 eV for this level of theory and the elongated FeO distance. In Figure 5(c), the PDOS obtained using the bulk Ud value and Up = 5.90 eV is plotted. Here, both Fe d- and p-orbitals show strong localization with limited overlap, consistent with the more endothermic value of Ef = 1.033 eV and an increase in FeO bond distance. In Figure 5(d), the PDOS for calculations using the surface-specific Ud and Up = 5.90 eV is plotted. Here, the bonding and antibonding pattern has disappeared completely, consistent with the unfavorable value of Ef = 1.508 eV and the longest FeO bond distance of 2.054 Å. The analysis of the oxidation of the Fe-term surface using the various DFT approaches compared in this study provides an explanation of why the ferryl surface stability is overestimated in DFT-GGA and GGA + U approaches that do not treat the O p-states or treat the surface Fe uniquely: Not using U at all results in overly delocalized Fe d- and O p-states, which artificially enhances Fe−O interactions at the surface. The addition of Ud on Fe can localize the Fe d-orbital and reduce the Fe−O overlap. However, the addition of Up in Figure 5(c) and (d) further reduces the d- and p-orbital hybridization. These results are the opposite from what we have seen in Table 1 and Figures 3 and 4 that adding Up can restore Fe−O bonding, and this highlights that the effect of adding Up can differ from the bulk to the surface. In the bulk system, each Fe (or O) has six (or four) bonds with its counterions. On the surface, the layer 0(O) only binds to the layer 1(Fe). When Coulomb repulsion terms are added on both Fe and O, the top O can be repelled from the surface to the vacuum space, reducing the interaction with surface Fe. Therefore, the formation of the ferryl surface becomes increasingly less favorable.

(7)

where EO2, Eferryl, and EFe‑term are the total energies of molecule O2, the ferryl surface, and the Fe-term surface, respectively. The factor of 1/2 takes into account that our structural models are symmetric slabs with two equivalent exposed surfaces. We note that this model reaction is a simplification of the ferryl group formation mechanism put forward by Jarvis and Chaka.83 However, eq 7 provides a convenient means of tracking how the various DFT and DFT + U methods used here affect heterogeneous reactivity, which we go on to relate to substrate electronic structure. The calculated values of Ef, along with the ferryl FeO bond distance, are reported in Table 4 for the different DFT + U styles. The range in Ef as a function of method is greater than 2 eV, with the results changing from being exothermic (for DFTGGA) to endothermic for all other methods. The most

■

CONCLUSIONS The surface phase diagram of α-Fe2O3(0001) was revisited using different styles of DFT + U calculations. In general, previous DFT and DFT + U studies of the surface predict a I



■

stable ferryl termination in the (1 × 1) surface, while most experiments exclude this surface or suggest such a termination is only metastable and is not observed outside of domains.15,19,26 The fact that DFT + U better predicts the band gap of bulk hematite but does not improve the surface phase diagram motivated us to consider variations in the value of Ud for chemically distinct Fe atoms exposed at the surface. The notion that U can change with system geometry or chemical environment has been previously explored in molecular systems by Kulik and co-workers, in which U was determined to be a function of molecular configuration. Following other theoretical studies of transition metal oxides, we also include a Hubbard term for the oxygen 2p states, Up. We find that the DFT + Ud+p surface free energies using derived/surface-specific Ud values and a constant Up value of 5.90 eV result in a reasonable phase diagram when compared to a range of experimental reports. Specifically, we do not find the ferryl termination to be the most preferred surface structure in a (1 × 1) cell, though this result does not counter the possible stability of ferryl domains in more disordered systems. Further supporting our approach is the fact that using DFT + Ud+p leads to universal improvement in the prediction of bulk properties of hematite. Taking into account the variation of Ud increases the computational cost. However, we note that this approach is still more tractable then the application of hybrid functionals to periodic systems, as the long-range nature of Hartree exchange can incur two orders of magnitude more computational cost.53 Another loose comparison between the approach of using varying Ud values to the use of hybrid functionals is the fact that there is no universal amount of Hartree exchange that can achieve improved results for all systems. In the present study, we have used the well-studied α-Fe2O3(0001) surface to assess what level of theory best predicts the surface-phase diagram and bulk properties, and this will enable an informed choice of methods for future studies of materials for which less experimental information is available, including defect surfaces that may better represent real systems.

■

ACKNOWLEDGMENTS This work was supported by National Science Foundation (NSF) Grant CHE-1254127 and the University of Iowa College of Liberal Arts and Sciences. We acknowledge the Arctic Region Supercomputing Center and University of Iowa High Performance Computing for computational resources.

■

REFERENCES

(1) Brown, G. E., Jr.; Calas, G. Environmental MineralogyUnderstanding Element Behavior in Ecosystems. C. R. Geosci. 2011, 343, 90−112. (2) Brown, G. E., Jr.; Calas, G. Mineral-Aqueous Solution Interfaces and Their Impact on the Environment. Geochemical Perspectives 2012, 1, 483−742. (3) Sivula, K.; LeFormal, F.; Grätzel, M. Solar Water Splitting: Progress Using Hematite (α-Fe2O3) Photoelectrodes. ChemSusChem 2011, 4, 432−449. (4) Yatom, N.; Neufeld, O.; Toroker, M. C. Toward Setting the Debate on the Role of Fe2O3 Surface States for Water Splitting. J. Phys. Chem. C 2015, 119, 24789−24795. (5) Iandolo, B.; Hellman, A. The Role of Surface States in the Oxygen Evolotion Reaction in Hematite. Angew. Chem. 2014, 126, 13622−13626. (6) Hellman, A.; Pala, R. G. S. First-Principles Study of Photoinduced Water-Splitting on Fe2O3. J. Phys. Chem. C 2011, 115, 12901−12907. (7) Trainor, T. P.; Chaka, A. M.; Eng, P. J.; Newville, M.; Waychunas, G. A.; Catalano, J. G.; Brown, G. E., Jr. Structure and Reactivity of the Hydrated Hematite (0001) Surface. Surf. Sci. 2004, 573, 204−24. (8) Nguyen, M. T.; Seriani, N.; Gebauer, R. Water Adsorption and Dissociation on α-Fe2O3 (0001): PBE + U Calculations. J. Chem. Phys. 2013, 138, 194709−1−8. (9) Chambers, S. A.; Yi, S. I. Fe Termination for α-Fe2O3(0001) as Grown by Oxygen-Plasma-Assisted Molecular Beam Epitaxy. Surf. Sci. 1999, 439, L785−L791. (10) Condon, N. G.; Leibsle, F. M.; Lennie, A. R.; Murray, P. W.; Vaughan, D. J.; Thornton, G. Biphase Ordering of Iron Oxide Surfaces. Phys. Rev. Lett. 1995, 75, 1961−1964. (11) Condon, N. G.; Leibsle, F. M.; Lennie, A. R.; Murray, P. W.; Parker, T. M.; Vaughan, D. J.; Thornton, G. Scanning Tunnelling Microscopy Studies of α-Fe2O3(0001). Surf. Sci. 1998, 397, 278−287. (12) Eggleston, C. M. The Surface Structure of α-Fe2O3 (001) by Scanning Tunneling Microscopy: Implications for Interfacial Electron Transfer Reactions. Am. Mineral. 1999, 84, 1061−1070. (13) Eggleston, C. M.; Stack, A. G.; Rosso, K. M.; Higgins, S. R.; Bice, A. M.; Boese, S. W.; Pribyl, E. D.; Nichols, J. J. The Structure of Hematite (α−Fe2O3) (001) Surfaces in Aqueous Media: Scanning Tunneling Microscopy and Resonant Tunnerling Calculations of Coexisting O and Fe Terminations. Geochim. Cosmochim. Acta 2003, 67, 985−1000. (14) Ketteler, G.; Weiss, W.; Ranke, W. Surface Structures of αFe2O3(0001) Phases Determined by LEED Crystallography. Surf. Rev. Lett. 2001, 8, 661−683. (15) Lemire, C.; Bertarione, S.; Zecchina, A.; Scarano, D.; Chaka, A. M.; Shaikhutdinov, S.; Freund, H.-J. Ferryl (FeO) Termination in the Hematite α-Fe2O3(0001) Surface. Phys. Rev. Lett. 2005, 94, 166101. (16) Liu, S.; Wang, S.; Guo, J.; Guo, Q. Polarity and Surface Structural Evolution of Iron Oxide Films. RSC Adv. 2012, 2, 9938− 9943. (17) Lübbe, M.; Moritz, W. A LEED Analysis of the Clean Surface of α-Fe2O3(0001) and α-Cr2O3(0001) Bulk Single Crystals. J. Phys.: Condens. Matter 2009, 21, 134010. (18) Shaikhutdinov, S.; Weiss, W. Oxygen Pressure Dependence of the α-Fe2O3(0001) Surface Structure. Surf. Sci. 1999, 432, L627− L634. (19) Tang, Y.; Qin, H.; Wu, K.; Guo, Q.; Guo, J. The Reduction and Oxidation of α-Fe2O3 (0001) Surface Investigated by Scanning Tunneling Microscopy. Surf. Sci. 2013, 609, 67.

ASSOCIATED CONTENT

* Supporting Information S

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcc.5b12144. Review of the DFT + U method and linear response approach to calculate both bulk and surface Ud values is presented, along with considerations for applying the perturbation method to 2D slab models. Plots of calculated Ud values as a function of the perturbation parameter, α, on all six surface models, are reported in Figure S.1. Convergence test results for the calculated Ud values as a function of slab model layer thickness for the Fe-term and O-term surfaces are reported in Figure S.2. Plots of the correlation between calculated specific Ud values and the average change in the Fe−O bond lengths are presented in Figure S.3 (PDF)

■

Article

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest. J


Article

The Journal of Physical Chemistry C (20) Wang, X. G.; Weiss, W.; Shaikhutdinov, S. K.; Ritter, M.; Petersen, M.; Wagner, F.; Schlogl, R.; Scheffler, M. The Hematite (αFe2O3) (0001) Surface: Evidence for Domains of Distinct Chemistry. Phys. Rev. Lett. 1998, 81, 1038−41. (21) Eggleston, C. M.; Stack, A. G.; Rosso, K. M.; Bice, A. M. Adatom Fe(III) on the Hematite Surface: Observation of a Key Reactive Surface Species. Geochem. Trans. 2004, 5, 33−40. (22) Barbier, A.; Stierle, A.; Kasper, N.; Guittet, M.-J.; Jupille, J. Surface Termination of Hematite at Environmental Oxygen Pressures: Experimental Surface Phase Diagram. Phys. Rev. B: Condens. Matter Mater. Phys. 2007, 75, 233406. (23) Tanwar, K. S.; Petitto, S. C.; Ghose, S. K.; Eng, P. J.; Trainor, T. P. Fe (II) Adsorption on Hematite (0001). Geochim. Cosmochim. Acta 2009, 73, 4346−65. (24) Catalano, J. G.; Fenter, P.; Park, C.; Zhang, Z.; Rosso, K. M. Structure and Oxidation of Hematite Surfaces Reacted with Aqueous Fe(II) at Acidic and Neutral pH. Geochim. Cosmochim. Acta 2010, 74, 1498−12. (25) Thevuthasan, S.; Kim, Y. J.; Yi, S. I.; Chambers, S. A.; Morais, J.; Denecke, R.; Fadley, C.; Liu, P.; Kendelewicz, T.; Brown, G., Jr. Surface Structure of MBE-Grown α-Fe2O3 (0001) by IntermediateEnergy X-ray Photoelectron Diffraction. Surf. Sci. 1999, 425, 276−286. (26) Yamamoto, S.; Kendelewicz, T.; Newberg, J. T.; Ketteler, G.; Starr, D. E.; Mysak, E. R.; Andersson, K. J.; Ogasawara, H.; Bluhm, H.; Salmeron, M.; Brown, G. E., Jr.; Nilsson, A. Water Adsorption on αFe2O3(0001) at Near Ambient Conditions. J. Phys. Chem. C 2010, 114, 2256. (27) Reuter, K.; Scheffler, M. Composition, Structure, and Stability of RuO2(110) as a Function of Oxygen Pressure. Phys. Rev. B: Condens. Matter Mater. Phys. 2001, 65, 035406. (28) Sun, Q.; Reuter, K.; Scheffler, M. Effect of a Humid Environment on the Surface Structure of RuO2(110). Phys. Rev. B: Condens. Matter Mater. Phys. 2003, 67, 205424. (29) Reuter, K.; Scheffler, M. Composition and Structure of the RuO2(110) Surface in an O2 and CO Environment: Implications for the Catalytic Formation of CO2. Phys. Rev. B: Condens. Matter Mater. Phys. 2003, 68, 045407. (30) Bergermayer, W.; Schweiger, H.; Wimmer, E. Ab Initio Thermodynamics of Oxide Surfaces: O2 on Fe2O3(0001). Phys. Rev. B: Condens. Matter Mater. Phys. 2004, 69, 195409. (31) Cohen, A. J.; Mori-Sánchez, P.; Yang, W. Insights into Current Limitations of Density Functional Theory. Science 2008, 321, 792− 794. (32) Rollmann, G.; Rohrbach, A.; Entel, P.; Hafner, J. First-Principles Calculation of the Structure and Magnetic Phases of Hematite. Phys. Rev. B: Condens. Matter Mater. Phys. 2004, 69, 165107. (33) Rohrbach, A.; Hafner, J.; Kresse, G. Ab Initio Study of the (0001) Surfaces of Hematite and Chromia: Influence of Strong Electronic Correlations. Phys. Rev. B: Condens. Matter Mater. Phys. 2004, 70, 125426. (34) Mochizuki, S. Electrical-Conductivity of α-Fe2O3. Phys. Status Solidi A 1977, 41, 591−4. (35) Anisimov, V. I.; Zaanen, J.; Anderson, O. K. Bond Theory and Mott Insulators: Hubbard U Instead of Stoner T. Phys. Rev. B: Condens. Matter Mater. Phys. 1991, 44, 943−54. (36) Anisimov, V. I.; Gunnarsson, O. Density-Functional Calculation of Effective Coulomb Interactions in Metals. Phys. Rev. B: Condens. Matter Mater. Phys. 1991, 43, 7570−7574. (37) Anisimov, V. I.; Aryasetiawan, F.; Lichtenstein, A. I. FirstPrinciples Calculations of the Electronic Structure and Spectra of Stongly Correlated Systems: The LDA + U Method. J. Phys.: Condens. Matter 1997, 9, 767−808. (38) Liao, P.; Toroker, M. C.; Carter, E. A. Electron Transport in Pure and Doped Hematite. Nano Lett. 2011, 11, 1775−1781. (39) Pan, H.; Meng, X.; Qin, G. Hydrogen Generation by Water Splitting on Hematite (0001) Surfaces: First Principles Calculations. Phys. Chem. Chem. Phys. 2014, 16, 25442−25448. (40) Liao, P.; Carter, E. A. Ab Initio DFT + U Predictions of Tensile Properties of Iron Oxides. J. Mater. Chem. 2010, 20, 6703−6719.

(41) Liao, P.; Carter, E. A. Ab Initio Density Functional Theory + U Predictions of the Shear Response of Iron Oxides. Acta Mater. 2010, 58, 5912−5925. (42) Canepa, P.; Schofield, E.; Chadwick, A. V.; Alfredsson, M. Phys. Chem. Chem. Phys. 2011, 13, 12826−12834. (43) Plata, J. J.; Marquez, A. M.; Sanz, J. F. Communication: Improving the Density Functional Theory+U Description of CeO2 by Including the Contribution of the O 2p Electrons. J. Chem. Phys. 2012, 136, 041101. (44) Huda, M. N.; Walsh, A.; Yan, Y.; Wei, S. H.; Al-Jassim, M. M. Electronic, Structural, and Magnetic Effects of 3d Transition Metals in Hematite. J. Appl. Phys. 2010, 107, 123712. (45) Hu, Z. P.; Metiu, H. Choice of U for DFT + U calculations for titaniun dioxides. J. Phys. Chem. C 2011, 115, 5841−45. (46) Cococcioni, M.; de Gironcoli, S. Linear Response Approach to the Calculation of the Effective Interaction Parameters in the LDA + U Method. Phys. Rev. B: Condens. Matter Mater. Phys. 2005, 71, 035105. (47) Lu, D.; Liu, P. Rationalization of the Hubbard U Parameter in CeOx from First Principles: Unveiling the Role of Local Structure in Screening. J. Chem. Phys. 2014, 140, 084101. (48) Kulik, H. J.; Marzari, N. Accurate Potential Energy Surfaces with a DFT + U(R) Approach. J. Chem. Phys. 2011, 135, 194105. (49) Kulik, H. J.; Marzari, N. Transition-Metal Dioxides: A Case for the Intersite Term in Hubbard-Model Functionals. J. Chem. Phys. 2011, 134, 094103. (50) Hsu, H.; Umemoto, K.; Cococcioni, M.; Wentzcovitch, R. FirstPrinciples Study for Low-Spin LaCoO3 with a Structurally Consistent Hubbard U. Phys. Rev. B: Condens. Matter Mater. Phys. 2009, 79, 125124. (51) Perdew, J. P.; Zunger, A. Self-interaction Correction to Densityfunctional Approximations for Many-electron Systems. Phys. Rev. B: Condens. Matter Mater. Phys. 1981, 23, 5048−79. (52) Svane, A.; Gunnarsson, O. Transition-Metal Oxides in the SelfInteraction-Corrected Density-Functional Formalism. Phys. Rev. Lett. 1990, 65, 1148−1151. (53) Gou, G. Y.; Bennett, J. W.; Takenaka, H.; Rappe, A. M. Post Density Functional Theoretical Studies of Highly Polar Semiconductive Pb(Ti1−xNix)O3−x Solid Solutions: Effects of Cation Arrangement on Band Gap. Phys. Rev. B: Condens. Matter Mater. Phys. 2011, 83, 205115. (54) Mattioli, G.; Alippi, P.; Filippone, F.; Caminiti, R.; Bonapasta, A. A. Deep versus Shallow Behavior of Intrinsic Defects in Rutile and Anatase TiO2 Polymorphs. J. Phys. Chem. C 2010, 114, 21694−21704. (55) Morgan, B. J.; Watson, G. W. Polaronic Trapping of Electrons and Holes by Native Defects in Anatase TiO2. Phys. Rev. B: Condens. Matter Mater. Phys. 2009, 80, 233102. (56) Mattioli, G.; Giannozzi, P.; Bonapasta, A. A.; Guidoni, L. Reaction Pathways for Oxygen Evolution Promoted by Cobalt Catalyst. J. Am. Chem. Soc. 2013, 135, 15353−15363. (57) Cao, C.; Hill, S.; Cheng, H. P. Strongly Correlated Electrons in the [Ni(hmp) (ROH)X]4 Single Molecule Magnet: A DFT+U Study. Phys. Rev. Lett. 2008, 100, 167206. (58) Lany, S.; Zunger, A. Generalized Koopmans Density Functional Calculations Reveal the Deep Acceptor State of NO in ZnO. Phys. Rev. B: Condens. Matter Mater. Phys. 2010, 81, 205209. (59) Hohenberg, P.; Kohn, W. Inhomogeneous Electron Gas. Phys. Rev. 1964, 136, B864−71. (60) Kohn, W.; Sham, L. J. Self-Consistent Equations Including Exchange and Correlation Effects. Phys. Rev. 1965, 140, A1133−8. (61) Giannozzi, P. e. a.; et al. QUANTUM-ESPRESSO: A Modular and Open-Source Software Project for Quantum Simulation of Materials. J. Phys.: Condens. Matter 2009, 21, 395502. (62) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865−8. (63) Vanderbilt, D. Soft Self-Consistent Pseudopotentials in a Generalized Eigenvalue Formalism. Phys. Rev. B: Condens. Matter Mater. Phys. 1990, 41, 7892−5. (64) Monkhorst, H. J.; Pack, J. D. Special Points for Brillouin-Zone Integrations. Phys. Rev. B 1976, 13, 5188−5192. K


Article

The Journal of Physical Chemistry C (65) Lo, C. S.; Tanwar, K. S.; Chaka, A. M.; Trainor, T. P. Density functional theory study of the clean and hydrated hematite (11¯02) surfaces. Phys. Rev. B: Condens. Matter Mater. Phys. 2007, 75, 075425. (66) Goffinet, C. J.; Mason, S. E. Comparative DFT Study of InnerSphere As(III) Complexes on Hydrated α-Fe2O3(0001) Surface Models. J. Environ. Monit. 2012, 14, 1860−71. (67) Knotek, M. L.; Feibelman, P. J. Ion Desorption by Core-Hole Auger Decay. Phys. Rev. Lett. 1978, 40, 964−967. (68) Rao, C. N. R. Chemistry of High Temperature Superconductors; World Scientific Publishing Co. Pte. Ltd.: Singapore, 1991. (69) Lany, S.; Zunger, A. Polaronic Hole Localization and Multiple Hole Binding of Acceptors in Oxide Wide-Gap Semiconductors. Phys. Rev. B: Condens. Matter Mater. Phys. 2009, 80, 085202. (70) Kiejna, A.; Pabisiak, T. Surface Properties of Clean and Au or Pd Covered Hematite (α-Fe2O3) (0001). J. Phys.: Condens. Matter 2012, 24, 095003. (71) Zaanen, J.; Sawatzky, G. A.; Allen, J. W. Band Gaps and Electronic Structure of Transition-Metal Compounds. Phys. Rev. Lett. 1985, 55, 418−21. (72) Catti, M.; Valerio, G.; Dovesi, R. Theoretical Study of Electronic, Magnetic, and Structural Properties of α-Fe2O3 (Hematite). Phys. Rev. B: Condens. Matter Mater. Phys. 1995, 51, 7441−7450. (73) Coey, J. M. D.; Sawatzky, G. A. A Study of Hyperfine Interactions in the System (Fe1−xRhx)2O3 Using the Mössbauer Effect. J. Phys. C: Solid State Phys. 1971, 4, 2386−2407. (74) Rozenberg, G. K.; Dubrovinsky, L. S.; Pasternak, M. P.; Naaman, O.; Bihan, T. L.; Ahuja, R. High-Pressure Structural Studies of Hematite Fe2O3. Phys. Rev. B: Condens. Matter Mater. Phys. 2002, 65, 064112. (75) Sato, Y.; Akimoto, S. Hydrostatic Compression of Four Corundum-Type Compounds: α−Al2O3, V2O3, Cr2O3 and α−Fe2O3. J. Appl. Phys. 1979, 50, 5285−5291. (76) Finger, L. W.; Hazen, R. M. Crystal-Structure and Isothermal Compression of Fe2O3,Cr2O3, and V2O3 to 50 Kbars. J. Appl. Phys. 1980, 51, 5362−5367. (77) Xu, Z.; Joshi, Y. V.; Raman, S.; Kitchin, J. R. Accurate Electronic and Chemical Properties of 3d Transition Metal Oxides Using a Calculated Linear Response U and DFT + U(V) Method. J. Chem. Phys. 2015, 142, 144701. (78) Sorkin, A.; Iron, M. A.; Truhlar, D. G. Density Functional Theory in Transition-Metal Chemistry: Relative Energies of LowLying States of Iron Compounds and the Effect of Spatial Symmetry Breaking. J. Chem. Theory Comput. 2008, 4, 307−315. (79) Jain, A.; Hautier, G.; Ong, S. P.; Moore, C. J.; Fischer, C. C. Formation Enthalpies by Mixing GGA and GGA + U Calculations. Phys. Rev. B: Condens. Matter Mater. Phys. 2011, 84, 045115. (80) Wang, L.; Maxisch, T.; Ceder, G. Oxidation Energies of Transition Metal Oxides Within the GGA+U Framework. Phys. Rev. B: Condens. Matter Mater. Phys. 2006, 73, 195107. (81) Mason, S. E.; Iceman, C. R.; Trainor, T. P.; Chaka, A. M. Density Functional Theory Study of Clean, Hydrated, and Defective Alumina (11¯02) Surfaces. Phys. Rev. B: Condens. Matter Mater. Phys. 2010, 81, 125423. (82) Stull, D. R.; Prophet, H. JANAF Thermochemical Tables, 2nd ed.; U.S. National Bureau of Standards (U.S. EPO, Washington, D.C.), 1971. (83) Jarvis, E. A. A.; Chaka, A. M. Oxidation Mechanism and Ferryl Domain Formation on the α-Fe2O3(0001) Surface. Surf. Sci. 2007, 601, 1909−14.

L


Surface-Specific DFT + U Approach Applied to Î±-Fe2O3(0001) - The

Recommend Documents