Prediction of Stable Surfaces of Metal Oxides through the Unsaturated Coordination Index

This study proposes the unsaturated coordination index, σ, as a potential descriptor of the stability of metal-oxide surfaces cleaved from bulk. The value of σ, the number of missing bonds per unit area, can be obtained very quickly using only crystallographic data, namely, the bulk geometry. The surface energies of various binary oxides, with and without atom relaxation, were calculated. Their correlations with σ had good coefficients of determination (R2) values, particularly in high-symmetry crystals. The proposed descriptor is very useful for an initial evaluation of stable metal-oxide surfaces without conducting any surface model calculations.


INTRODUCTION
Chemical reactions take place on surfaces of solid materials in heterogeneous catalysis. Much effort has been devoted to gaining atomic-scale understandings of detailed surface structures and reaction mechanisms, which is necessary for improving catalytic activity or the rational design of novel heterogeneous catalysis. 1−14 Among solid materials, metal oxides are widely utilized in heterogeneous catalysis because of their high stabilities under experimental conditions in industrial processes, such as high temperature, pressure, and concentration of gaseous components. Slab-and-vacuum models (hereafter "slab models") have been used to model solid surface structures within the three-dimensional periodic boundary condition to explore plausible active sites and reaction mechanisms typically based on density functional theory. A large number of cleaved terminations can be generated within a relatively small supercell size limit in lowsymmetry crystals. Compared to metal surfaces, 13,14 the surface structures of metal oxides could be very complex. 15−18 However, most theoretical works focus on only the simplest termination, although slightly less stable surfaces are recognized to contain catalytically active sites in a limited number of individual cases. Major obstacles to the comprehensive screening of metal-oxide surfaces are the requirement of huge computational costs and the difficulty in modeling numerous accessible surfaces with different terminations. One of the authors developed algorithms to automatically generate nonpolar slab models using symmetry information. 19−21 These efforts enabled calculations of energetic and electronic properties of slab models, including surface energies (E surf ), surface defect energies (e.g., surface anion vacancy formation energies), small-molecule adsorption energies, ionic potentials and electron affinities. 20,22−32 This study proposes the unsaturated coordination index, σ, which is the number of missing bonds, or unsaturated coordination, per unit area. This value is derived from geometrical parameters of bulk and is designed as a descriptor of E surf . The validity of σ was verified using 12 terminations, including reconstructed surfaces, of rocksalt structure MgO and CaO, six terminations of anti-fluorite structure Li 2 O, Na 2 O, and K 2 O, 16 terminations of rutile structure TiO 2 , SnO 2 , and GeO, 18 terminations of anatase TiO 2 , 67 terminations of isostructural θ-Al 2 O 3 and β-Ga 2 O 3 , and 68 terminations of the monoclinic (m-) baddeleyite structure ZrO 2 . Good correlations between σ and E surf were found when atom positions were fixed to cleaved bulk, and the correlations were reasonable when atom positions were relaxed. The accessible surfaces of rutile and anatase structured TiO 2 , SnO 2 , GeO 2 , and m-ZrO 2 that do not spontaneously decompose into macroscopic facets of terminations with different orientations are also reported here.

E surf Calculations.
First-principles calculations were conducted using the projector augmented-wave method 33 as implemented in the VASP code. 34,35 The plane-wave basis set cutoff was 400 eV. As in our previous studies, 26,29 the Perdew− Burke−Ernzerhof functional tuned for solids (PBEsol) 36 within the generalized gradient approximation (GGA) was used because it provides reasonable bulk energetics and crystal structures, for instance, compared to the standard PBE-GGA functional 37 as shown in a previous systematic study of polymorphs of groups I−VI binary oxides. 38 Furthermore, the Hubbard U based on Dudarev's formulation 39 was additionally considered on 4d states of Ti and Zr. The effective U value of U−J (denoted as U eff ) was set at 3 eV, which was used to obtain fitted elemental-phase reference energies (FERE) by Stevanovićet al. 40 The +U correction was applied to Ti and Zr even though Ti and Zr in TiO 2 and ZrO 2 have formally d 0 electronic configurations because corrections would be necessary when treating defects where electrons could locally enter the d states. Internal coordinates and lattice parameters were relaxed in bulk calculations, and lattice parameters were fixed in slab model calculations. Spin polarization was considered in slab calculations with an initial ferromagnetic spin ordering.
Slab models under the three-dimensional periodic boundary condition were used to analyze surfaces, where slabs infinitely extending parallel to the ab-plane alternate with vacuum layers along the c-axis. The slab thickness in the models was larger than 14.5 Å and three repeat units, and the minimum vacuum thickness was 12 Å. In addition to cubic and tetragonal crystals, monoclinic β-Ga 2 O 3 , θ-Al 2 O 3 , and ZrO 2 were considered in this study because their low symmetry and existence of inversion centers result in a large number of nonpolar surfaces within a reasonable range of the basal area of the slab model, S. The E surf is defined as where E slab and E bulk are the energy of the slab without defects and the energy of the constituents of the slab when in a perfect bulk, respectively. The coefficient of 2 accounts for both sides of the slab. Note that all used slab models are stoichiometric, and each surface (top and bottom) is symmetric.

Definition of the Unsaturated Coordination Index (σ).
The unsaturated coordination index, σ, is defined as follows. A stoichiometric slab-and-vacuum model (slab model) is assumed, and bonds form between pairs of atoms with a distance shorter than an arbitrary cutoff, r cut . The coordination number of atom i, CN i , is the number of bonds of atom i, or in other words, the number of atoms within a sphere of radius r cut around atom i. The average CN of bulk, ⟨CN⟩ bulk , is given by where M is the number of atoms in the bulk crystal. The number of missing bonds in a slab model, ΔCN, is the number of missing CN compared to perfect bulk. When there are N atoms in the slab model The unsaturated coordination number is ΔCN normalized by the surface area of the slab model, which is = S CN/2 (4) where the factor 2 accounts for both sides of the slab.
The value of r cut used to obtain CN i can be any value larger than the longest "NN" bond and smaller than the shortest "second NN" bond length. Practically, the "nearest neighbor (NN)" bond length is not a fixed value but is within a certain range, except in the simplest crystals ( Table 1). The calculated rutile structure TiO 2 has two "NN" bond lengths of 1.974 and 1.998 Å. Similarly, the "second NN (2NN)" may have a range. The "2NN" in rutile TiO 2 is a O−O bond with length 2.565 Å. There is also a 2.809 Å O−O bond, which might be considered as part of the extended "2NN" bond or a "third NN" bond, but in any case, r cut should be smaller than 2.565 Å and bond lengths larger than this are irrelevant. Values of r cut in this study can be any value within R cut1 and R cut2 , which are arbitrary lengths larger and smaller than NN and 2NN lengths, respectively. Here, R cut1 and R cut2 were obtained with 0.1 Å intervals.
A robust method to identify missing bonds for any cleaved termination was proposed by Mackenzie et al., 41,42 where bonds are defined as vectors and their spatial relations are evaluated against a "dividing plane", but our procedure is more straightforward because the initial coordination of each atom does not need to be tracked and information on the direction of the bonds is not necessary. Figure 1a is a schematic showing how σ is derived. When a surface structure is obtained by cleaving the bulk, the surface atoms should miss bonds from the bulk, generating surface  atoms with unsaturated coordination. At this point, the surface terminations, which miss more bonds, tend to be unstable; σ reflects the density of missing bonds among the surface atoms per surface area. Figure 1b,c shows (110) and (102) surfaces of rutile structure TiO 2 as the representatives of high and low σ terminations. A large σ indicates a greater extent of unsaturated coordination, or a higher density of missing bonds, which is expected to result in a higher E surf .

RESULTS AND DISCUSSION
3.1. Details of Slab Models. 3.1.1. MgO and CaO. The space group of these rocksalt structure compounds is Fm3̅ m. Slab models with the minimum choice of S less than four times that of the smallest S, which is for the (111) orientation, were analyzed. There is a total of six terminations of type 1 in Tasker's definition 43 and nonpolar type A in the definition by Hinuma et al. 19 Additionally, the (111) surface is polar (type 3) according to Tasker but can be made nonpolar by surface reconstruction of a supercell. Three types of reconstruction patterns were considered. One is the "row" (R) pattern where one atom wide rows of one-half of surface atoms are removed. In the "zigzag" (Z) pattern, the surface atoms form a zigzag pattern when viewed from the direction normal to the surface. From another point of view, two atoms wide rows of one-half of surface atoms are removed. The "octopolar" (O) pattern has 2 × 2 ordering, where 75 and 25% of atoms in the outermost layer and second layer from the surface are removed, respectively. 44 The outermost layer can be anions (A) or cations (C), and thus six reconstructions (combinations of A or C and R, Z, or O) were studied. The geometries of slab models and values of E surf with fixed and relaxed atom positions are given in Tables S1 and S2, and the terminations are shown in Figure S1.
The lowest E surf surface is (100), followed by (310), (210), and (110). As discussed afterward, the (h10) surfaces may be viewed as a vicinal surface of steps and terraces, with the step vector parallel to [010] and edge vector parallel to [001], on the (100) terrace plane. A smaller h (>0) results in a higher density of steps and, therefore, a higher E surf . The octopolar pattern has the lowest E surf among the (111) reconstructions, as expected. 45 The choice of the outermost element (anion or cation) does not change σ in these terminations, and the difference in E surf was relatively small at less than 10% in the terminations considered in this study. The σ of the (111) terminations is high and does not affect the prediction of low E surf terminations.

Li 2 O, Na 2 O, and K 2 O.
The space group of these antifluorite structure compounds is also Fm3̅ m. Slab models with the minimum choice of S less than four times that of the smallest S, which is for the (111) orientation, were analyzed. There are three type 1 and nonpolar type A terminations, (110), (211), and (310), and three Tasker type 2 and nonpolar type B terminations, (111), (311), and (331). The geometries of slab models and values of E surf with fixed and relaxed atom positions (fixed and relaxed E surf , respectively) are given in Tables S3−S5, and the terminations are shown in Figure S2. All missing bonds at the lowest E surf termination, (111), are perpendicular to the surface and there are no bonds parallel to the surface. In contrast, all missing bonds at the third lowest E surf termination, (110), are not perpendicular to the surface and there are bonds parallel to the surface. The second lowest E surf termination, (331), has (111) termination-type steps on the (110) terrace surfaces.
3.1.3. Rutile Structure TiO 2 , SnO 2 , and GeO 2 . The space group of these rutile structure compounds is P4 2 /mnm. Slab models with the minimum choice of S less than four times that of the smallest S, which is for the (001) orientation, were analyzed. 14 orientations and 16 terminations were obtained, which are all Tasker type 2 and nonpolar type B except for the (001) termination that is Tasker type 1 and nonpolar type A. The geometries of slab models and values of fixed and relaxed E surf are given in Tables S6−S8, and the terminations are shown in Figure S3.
Surfaces subject to facet decomposition were investigated using the procedure by Hinuma et al. 46 Surface optimization using the USPEX code 47 Tables S6−S8. 3.1.4. Anatase Structure TiO 2 . The space group of anatase TiO 2 is I4 1 /amd. Slab models with the minimum choice of S less than four times that of the smallest S, which is for the (001) orientation, were analyzed. 13 orientations and 18 terminations were obtained. There are three Tasker type 1 and nonpolar type A terminations, (100), (110), and (310), and the rest are Tasker type 2 and nonpolar type B terminations. The geometries of slab models and values of fixed and relaxed E surf are given in Table S9, and the terminations are visualized in Figure S4. Surfaces subject to facet decomposition were investigated as in Section 3.1.3. The terminations that do not decompose into facets are in the order of increasing relaxed E surf : (101) and (001).
3.1.5. β-Ga 2 O 3 and θ-Al 2 O 3 . Compounds β-Ga 2 O 3 and θ-Al 2 O 3 share the same crystal structure with space group C2/m. Slab models with the minimum choice of S less than four times that of the smallest S, which is for the (100) orientation, were considered. There are 34 orientations and 67 terminations each for both Ga 2 O 3 and Al 2 O 3 . The obtained surfaces are type 2, or nonpolar type B except for the (010) surface that was Tasker type 1 and nonpolar type A, respectively. The geometries of slab models are given in Tables S10 and S11, and the terminations are illustrated in Figures S1−S9 of ref 29. 3.1.6. m-ZrO 2 . The space group of m-ZrO 2 is P21/c. Slab models with the minimum choice of S less than three times that of the smallest S, which is for the (001) orientation, were analyzed. 28 orientations and 68 terminations were obtained, which are all Tasker type 2 and nonpolar type B. The geometries of slab models and values of fixed E surf are given in Table S12, and the terminations are shown in Figures S5−S9.
Surfaces subject to facet decomposition were investigated according to the procedure in Hinuma et al. 46 Slab calculations with atom positions relaxed could not be calculated for a number of slabs in Table S12. Slabs with different thicknesses were available to obtain E surf in some of those cases, and the geometries of slabs for relaxation calculations are given in Table S13. The terminations that did not decompose into facets are in the order of increasing relaxed E surf : (111̅ )A,  Table S12. Here, surface optimization using the USPEX code 47−50 was conducted for the third−sixth lowest E surf terminations, (211̅ ), (102), (110), and (112̅ ). To reduce the computational cost, the General Utility Lattice Program (GULP) code 51,52 was used as the energy calculator in USPEX calculations together with interatomic potentials by Woodley et al. 53 Considering the four USPEX-relaxed terminations results in facet decomposition of (102̅ ) to (112̅ ) and (11̅ 2̅ ). Figure 2 plots the fixed E surf against σ, which would reflect the effect of σ more than the relaxed ones. Cubic crystals are discussed first. The coefficients of determination (R 2 ) in rocksalt structure MgO and CaO are good at 0.85 and 0.80, respectively (Figure 2a,b). The nonpolar type C (111) terminations, which appear as three distinct dots at the same σ, are below the trend of the remaining nonpolar type A terminations. The R 2 for the nonpolar type A terminations only are 0.92 and 0.91 for MgO and CaO, respectively. The reconstructions for the (111) terminations cannot be distinguished using σ, but the highest E surf termination among the six (111) terminations is 27 and 38% larger than the lowest E surf termination in MgO and CaO, respectively. Further estimation of the E surf requires additional descriptors from alternative approaches. The  (Figure 3a−e, respectively) except for the two highest relaxed E surf terminations in MgO. The R 2 of relaxed E surf is lower than that of the fixed E surf , but the order of relaxed E surf can still be estimated reliably using bulk crystal structure information only.

Correlation of σ and E surf .
An excellent R 2 is not necessarily good news for σ. If every missing bond carries the same energy penalty, the energy penalty per unit surface area when atoms are fixed to cleaved bulk, E surf , should be proportional to the density of missing bonds per unit surface area, which is σ. In other words, E surf = kσ should hold, where k is a crystal-dependent coefficient. This is clearly not the case in cubic crystals. For instance, the relation between E surf and σ of MgO is E surf = 2908σ − 239. The range of E surf and σ lies between roughly 0 and 300 meV/ Å 2 and 0.1 and 0.2 Å −2 , respectively, while the intercept of the regression line at σ = 0 is extremely low at −239 meV/Å 2 .
The E surf of (100), (310), (210), and (110) can be described very well using a steps-on-terraces model, which is informative in understanding the possible cause of the nonzero intercept. Figure 4a shows the clean (100) surface, and Figure 4b  of the lattice parameter a, d, the number of atoms with one and two missing bonds per terrace area of a 2 , b 1 , and b 2 , and surface energy per terrace area, surf , are shown for MgO, with atom positions fixed to cleaved bulk, in Table 2. Note that although we refer to E′ surf as the "surface energy", it is defined as the energy penalty per terrace area upon surface formation, to be exact, not the surface energy of the terrace (100) and vicinal (310), (210), and (110) surfaces. Atoms with one missing bond are on the terraces, whereas those with two missing bonds are at the step edges, as illustrated in Figure 4a,c. The ratio of step edge atoms increases, and therefore b 2 increases, as h in (h10) decreases. E′ surf vs b 2 is shown for MgO and CaO in Figure 4e. There is a clear linear relationship, which is consistent with the assumption that E′ surf is determined by a sum of energy penalties in the linear form E′ surf = E 1 b 1 + 2E 2 b 2 . Here, E 1 and E 2 are energies required to break a bond in atoms that end up with one and two missing bonds, respectively. Using b 1 + b 2 = 4 based on the number of atoms per terrace area of a 2 , the above equation can be rewritten as E′ surf = 4E 1 + {E 1 + 2(E 2 − E 1 )}b 2 . The value of E 1 is obtained from the intercept at b 2 = 0.
In a hypothetical world where E 2 = E 1 , the surface energy per terrace area of a 2 is E′ surf_hyp = 4E 1 + E 1 b 2 using E 1 , which is shown with empty symbols in Figure 4e. The slope for the actual surface energy, E 1 + 2(E 2 − E 1 ) is substantially larger than that for the hypothetical surface energy, E 1 . This result suggests that breaking the second bond has a much larger energy penalty than breaking the first bond (E 2 ≫ E 1 ), which is consistent with chemical intuition. The hypothetical surface energy per surface area, surf hyp , for MgO and CaO is plotted against σ in Figure 4f using empty symbols. The linear fit of empty symbols in Figure 4f passes the origin by definition. These results demonstrate that the differences in bond breaking energies between various terminations are reflected as a linear relation between E surf and σ, suggesting that σ is an excellent descriptor of E surf in these rocksalt structure crystals despite the fact that the original underlying assumption that the same energy is required to break any bond does not hold. A more substantial intercept means that the assumption is less valid. The steps-onterrace analysis clarified the origin of the high R 2 linear fit with a very low intercept, but such discussion can be applied only to

ACS Omega
http://pubs.acs.org/journal/acsodf Article a subset of terminations with steps in the same direction on the same terrace termination. An infinite combination of step directions and terrace orientations is necessary to describe the full set of orientations. Therefore, using additional structural information, such as the coordination of atoms at the surface, was not considered in our search for E surf descriptors.
To evaluate the performance of this new explanatory variable, we tested it for more complex oxide surfaces in addition to the cubic systems with high symmetry. In tetragonal crystals, the two lowest fixed E surf terminations, (110) and (100), in rutile structure TiO 2 , SnO 2 , and GeO 2 are the terminations with the smallest and second smallest σ (Figure 2f−h, respectively), and the lowest fixed E surf termination in anatase TiO 2 , (101)A, has the smallest σ ( Figure 2f). The range of points in rutile and anatase structure TiO 2 are shown together in Figure 2f and these overlap, suggesting that the bond strengths are similar in the two polymorphs. The R 2 is between 0.64 and 0.72, which is lower than the cubic crystals but still reasonable. The results using relaxed E surf (Figure 3f−h) show similar trends but with slightly lower R 2 .
In addition, monoclinic crystals were also studied. The lowest fixed E surf termination in Al 2 O 3 and Ga 2 O 3 , (100)A, has the fourth smallest σ (Figure 2i,j, respectively). The lowest fixed E surf termination of ZrO 2 , (111)A, has the smallest σ (Figure 2k). These statements also hold for relaxed E surf (Figure 3i−k). The R 2 for fixed E surf is relatively high at between 0.58 and 0.75, and relaxing atoms reduces R 2 to In panels (e) and (f), filled symbols indicate actual surface energies, and the empty symbols represent hypothetical surface energies where the surface energy is proportional to the missing bond density and the energy required to break a bond is always the same as that in the (100) termination. Table 2. Information Used for the Analysis of MgO in Figure 4 termination The E surf vs σ plots discussed above tend to correlate linearly for both fixed and relaxed E surf . Therefore, the coefficient of determination (R 2 ) was discussed primarily as a measure of error. The descriptor σ is useful when the order of terminations sorted by E surf is close to that sorted by σ, and thus a linear correlation is not necessary for this purpose. However, a linear correlation helps the estimation of the absolute value of E surf for a specific surface with a certain σ when the linear trend between E surf and σ was derived from some trial calculations. Other statistical methods, such as the (absolute) mean error and the standard deviation of the error from the linear fit, may be scientifically interesting but are not considered because practical usage is very limited.
The five terminations, in ascending order, with smallest σ, five lowest fixed E surf , and five lowest relaxed E surf for each oxide, are summarized in Tables S14−S25. There is excellent consistency in the rocksalt and fluorite structure oxides, and the order is reproduced somewhat in rutile and anatase structure oxides. However, the prediction performance is limited in Ga 2 O 3 , Al 2 O 3 , and ZrO 2 with monoclinic symmetry because there are more terminations with close E surf compared to high-symmetry crystals. The lowest fixed E surf have very low relaxed E surf in the studied oxides, and thus once the former is correctly identified, the latter is identified too. Figure 5 plots σ against the ratio of E surf between relaxed and fixed atom positions. The cubic crystals (Figure 5a−e) have high R 2 values of more than 0.8, whereas those of the others are poor at about 0.3 at most. Intuitively, a larger unsaturation of coordination, or larger σ, implies more room for relaxation, thereby reducing the ratio. This is clearly observed in the cubic crystals but not obvious in other crystals. One reason may be that the high symmetry of cubic crystals limits the variation of relaxation routes; existence of diverse relaxation opportunities would result in different extents of relaxation for surfaces with similar E surf or σ and thereby worsening the correlation between the ratio and σ.
Finally, we discuss another limitation of σ. The values of fixed and relaxed E surf are plotted together for all considered surfaces in Figures S10 and S11, respectively. Points for the rocksalt and anti-fluorite structures are consistently below points for the other structures, and thus the absolute value of E surf cannot be estimated from the value of σ.

SUMMARY
The unsaturated coordination index, σ, was defined as the number of missing bonds, or unsaturated coordination, per unit area. The validity of σ as a descriptor of E surf was verified using 12 terminations, including reconstructed surfaces, of rocksalt structure MgO and CaO, six terminations of antifluorite structure Li 2 O, Na 2 O, and K 2 O, 16 terminations of rutile structure TiO 2 , SnO 2 , and GeO, 18 terminations of anatase TiO 2 , 67 terminations of isostructural θ-Al 2 O 3 and β-Ga 2 O 3 , and 68 terminations of monoclinic baddeleyite structure ZrO 2 . Very good correlations of 0.80 or higher between fixed E surf and σ were found for the cubic crystals, and the correlation was still good at 0.58 or higher for other crystals. The lowest E surf termination, for both fixed and relaxed E surf , consistently had small σ. The underlying assumption is that σ is proportional to fixed E surf when breaking of any bond carries the same energy penalty, which is not obviously true, but still σ turns out to be a reasonable descriptor of both fixed and relaxed E surf . The developed σ contributes to the computational design of solid materials, such as heterogeneous catalysts, from unexplored surfaces of metal oxides. While the user-defined parameters R cut1 and R cut2 continue to introduce a level of arbitrariness in the missing bonds, the development of a simplified descriptor is necessary to describe crystal structures effectively. The proposed descriptor σ is intrinsically more effective when the nature of bonds is not diverse. Alternate descriptors need to be designed to predict low E surf with further confidence, especially in crystals with diverse types of bonds. ■ ASSOCIATED CONTENT
Information on the used slab models; structure of the considered surface; and plots of fixed and relaxed E surf for all considered surfaces against σ (PDF) ■