Optimization of the facet structure of transition-metal catalysts applied to the oxygen reduction reaction

M. Núñez, J. L. Lansford and D. G. Vlachos*
Predicting the optimal structure for a catalytic material has been a long-standing goal, but typically an arbitrary active site on a uniform surface is modelled. Identification of the most-active facet structure for structure-sensitive chemistries, such as the oxygen reduction reaction, is lacking. Here we develop an approach to predict the optimal structure of a catalytic material by identifying the active site and identifying the density and spatial arrangement of such sites while minimizing the surface energy. We find that the theoretical peak performance predicted by linear scaling relations is unattainable because of the lack of suitable active sites on low-index planes, as well as geometric and stability constraints. A random array of vacancies results in a modest performance enhancement compared to ideal facets, whereas defect sites with a maximum density in disordered structures significantly increase the catalyst performance. We applied this methodology to the oxygen reduction reaction on defected Pt(111), Pt(100), Au(111) and Au(100) surfaces.ighly dispersed metal nanoparticles are used ubiquitously as catalysts. Substantial effort has been devoted over the past two decades to determine how reactions are affected by the
size and shape of nanoparticles and how best to maximize catalyst activity1–3. Techniques for the synthesis of shape-selecting nanopar- ticles have been developed4,5 and applied to various chemistries, which include formic acid oxidation6, acetylene hydrogenation7, NO reduction8, benzene hydrogenation9, trans-to-cis isomeriza- tion10 and 2-propenol oxidation11 to mention a few. Despite these advances, engineering the nanoparticle size and shape from first principles has remained elusive. In recent exciting studies on oxygen reduction reaction (ORR) on Pt(111) (ref. 12) and PtNi/C nanoparticles13, it was shown that surface defects can enhance activity. Similarly, a multiscale computational study demonstrated that patched bimetallic surfaces with vacancies on the top layer can deliver an unprecedented boost in activity for the ammonia decomposition reaction14; the size and shape of defects are obvi- ously important, but how to predict the optimal structures remains elusive. These early studies underline a rather novel and unexplored means to enhance catalyst activity whereby one purposely engineers defects on catalyst surfaces. Arbitrary defects, by contrast, may lead to a decrease in activity15.
The vast structure–activity design space, in which surfaces are purposely engineered with atomistic detail and their catalytic properties are computed on the fly, imposes a major computational challenge but promises to unlock unexplored structures with an improved performance. Due to the computational cost of first-prin- ciples methods, such as density functional theory (DFT), the direct evaluation of the activity of numerous structures is beyond the cur- rent supercomputer power. Instead, a surrogate microstructure– energy model is necessary. A recent breakthrough in correlating the binding energy with microstructure entails the work of Sautet and co-workers who related the binding energies of ORR interme- diates to the generalized coordination number (GCN) of a binding site12,16,17. An earlier study that described the adsorption of CO on Au

clusters with the coordination number of the binding site (a geomet- ric descriptor) provides further evidence that activity and microge- ometry can be correlated18. These structure-scaling relations enable the high throughput screening of microstructures by computing GCNs simply through counting the first- and second-nearest neigh- bours. By constructing volcano curves versus the GCN, Sautet and co-workers17 showed that an optimal site on a Pt(111) crystal exists and demonstrated experimentally that defects, indeed, increase the catalyst activity compared to that of the ideal surface. However, the experimental activity enhancement was modest compared to that of a structural model consisting of identical, single active sites that cor- responded to the volcano peak. This raises fundamental questions as to the predictive ability of models and whether microstructures that consisted of closely packed optimal active sites were experi- mentally made and stable, given that highly defected structures should possess higher surface free energies than ideal low Miller index planes. Despite these exciting developments, the maximum density of defects one can geometrically arrange in a structure, the effect of spatial organization of defects and the maximum activity achieved in an experiment remain unclear. Addressing these ques- tions requires further conceptual advances beyond the structure- based volcano curve. Hanselman and Gounaris19 proposed active sites given a descriptor, but did not consider stability requirements, which are critical, given that previous studies related active site sta- bility to activity and, in some cases, found a trade-off20,21.
Here we introduce a framework that can predict the microstruc- ture that results in the maximum activity of a catalyst surface along with its stability by using a multiobjective optimization approach. We designed the atomistic structure dynamically in terms of both optimal site activity and surface stability and addressed funda- mental questions that may close the gap between models and experiments. We demonstrated the approach on Pt(111), Pt(100), Au(111) and Au(100) facets for the ORR. The ORR at the cathode of proton membrane fuel cells is an important structure-sensitive chemistry22–25. Currently, the widespread commercialization of Pt/C

Catalysis Center for Energy Innovation (CCEI), Department of Chemical and Biomolecular Engineering, University of Delaware, Newark, DE, USA.
*e-mail: [email protected]

reduces an oxygen species, such as *OH and *OOH, with the addi- tion of an electron and hydrogen ion. It is well known that the binding energies of the *OH and *OOH intermediates vary among facets34–36. Based on the structure-scaling relations, the activity of a site depends on its coordination environment. Supplementary Section 1 describes how a site’s activity is computed based on structure-scaling relations using the GCN and experimental data. Structure-dependent volcano plots for Pt and Au versus the GCN are shown in Fig. 1. *OH removal is rate determining on the left leg of the volcano while *OOH formation is rate determining on the right leg. Au is a weak binding metal and has an optimal GCN of 5.75, which is smaller than the GCN of the (100) and (111) planes. In contrast, Pt binds oxygen species strongly and its activity peaks at a GCN of 8.29, which is higher than the value on either perfect facet. Clearly, ideal facets are suboptimal and to create defects can be beneficial. Figure 1 indicates that potentially different strategies of crystal imperfection may be necessary to maximize the activity
12 of weakly and strongly binding metals that lie on either side of the volcano curve. We elaborate on this point below.
The structures of the defected Pt(111), Pt(100), Au(111) and

Fig. 1 | Volcano map for ORR activities on Pt and Au. Current i is defined in equation (1). The GCNs of the optimal site are 5.75 and 8.29 for Au and Pt, respectively. The vertical dashed lines at GCN = 6.67 and 7.5 show the GCNs of the (111) and (100) planes for reference. *OH des and *OOH form refer to areas of the volcano map in which OH desorption and OOH formation, respectively, are the rate-controlling steps. Clearly, ideal facets
are suboptimal; for Au and Pt under- and overcoordinated sites are needed to enhance activity.

catalysts is hindered because of slow kinetics (a high overpotential), the high cost of Pt and catalyst degradation under harsh acidic con- ditions26–28. Consequently, there remains ample room for catalyst design to advance proton membrane fuel cell technology29,30. We report that the optimal active site changes among metals and facets of a single metal and with the density of active sites. Our approach systematically quantifies both how and why the active site changes as the natural coordination environment of the surface changes. This is an important method development for the catalysis community, as traditional methods of design31–33 compare only similar sites and surfaces to each other across different metals. Rather surprisingly, given the large number of distinct atomic configurations that lead to similar coordination environments, we found that only a small number of geometrically distinct sites contribute to the activity on defected low index planes. Disordered structures of such (highly packed) active sites deliver the best performance. Expectedly, the creation of defects increases the surface free energy. However, numerous metastable structures that possess lower activity than that of the peak of the volcano curve exist. Combined, the lack of sites on low index planes that match the peak of the volcano curve, packing constraints and surface stability rationalize that (geometric, volcano-based) peak activity is unattainable. Importantly, meta- stable structures are still an order of magnitude (or more) superior to ideal facets. In contrast, vacancies in the top layer at the same density as those of our optimal structures provide only a modest performance enhancement. To make the right geometric motif of the defects and to closely pack such defects are essential and key to materializing the benefits of the defect engineering of crystal sur- faces and provide a potential explanation of the modest activity seen experimentally. We introduce uncertainty quantification in the rate and optimal microstructures that arise from solvation, surface cov- erages and structure-scaling relations.
Results and discussion
Active site and optimal structures. The ORR mechanism is mod- elled as a series of four sequential elementary steps, each of which

Au(100) facets that maximize the current density and also minimize the surface energy are predicted using a two-level optimization pro- cedure. First, the simultaneous optimization of activity and stabil- ity is carried out (by weighing activity versus stability through a weighting factor, ω), as detailed in Methods. After this optimization, the surface is quenched to a local energy minimum on the coarse- grained potential energy surface (Methods). This quenching step accounts for the fact that very high surface energy sites will not exist under the reaction conditions if the surface can easily restructure to a lower energy state. Surface energies and current densities for Pt(111)-based metastable structures that result from the first opti- mization are shown in Fig. 2 (red circles). The other surfaces follow similar trends and are included in Supplementary Fig. 7. A nearly continuous set of solutions between the ideal facet (ω = 0) and the most-active (but not necessarily stable, ω = 1) structures is obtained. For ω > ~0.4, the activity and surface energy vary only slightly with ω; primarily, activity is optimized. In contrast, for low values of ω < ~0.2, the stability of the structure dominates optimization. The transition between the two regimes is fairly abrupt. Subsequent quenching to a local minimum energy lowers the surface energy at the expense of the activity, typically by a factor of 50–75% (Fig. 2, blue circles). The diffusion of surface atoms eliminates some active sites created during the multiobjective optimization.
The results indicate a profound effect of introducing defects to a structure whereby the catalyst activity can increase by nearly two orders of magnitude compared to the stable low-index plane. At the same time, they indicate that the stability and details of the micro- structure matter. We elaborate on these topics below.
The highest activity structures after the quenching step are shown in Fig. 3. Analysis of the structures identifies the active sites, shown in Fig. 4. For Pt(111), several types of cavities are active. Figure 4a shows a five-atom vacancy that exposes two adjacent active sites with a GCN of 8.25, which is near the Pt peak of 8.29. Structures at intermediate values of ω have this cavity as the predominant active site. A four-atom vacancy, shown in Fig. 4b, has a slightly higher GCN of 8.33 and is more active, but contributes less to the overall activity because there is only one active site in the cavity, that is, due to a lower density of such defects. Structures at a low ω consist of disjoint cavities (Fig. 4a,b). At higher activities, a denser packing of the active sites is required, which results in a predominant cav- ity, as shown in Fig. 4c. This includes vacancies, which branch out into other cavities, and allow for a dense network of active sites. The most-active structures consist of a disordered ensemble of these active sites, as can be seen in Fig. 3.
Calle-Vallejo et al.17 proposed that the active site on a defected Pt(111) was either a six-atom or five-atom vacancy that exposed
from low-activity/low-energy structures to high-activity/high- energy structures. All the Pareto fronts have linear portions, which are indicative of a single type of active site being added in a higher density as the activity increases. The constant slopes indicate pro- portional increases in surface energy and activity. The Pt(100) and Au(100) curves are entirely linear, which indicates that only one type of active site is present, consistent with the structural analysis presented above (Fig. 4). In contrast, the Pt(111) and Au(111) Pareto fronts have kinks that show the active site being added to the surface changed with an increasing activity, consistent with the structural analysis above. The Pareto front for Au(100) has a large gap. Only two structures are relevant: the ideal facet and the periodic struc- ture shown in Fig. 4g. Periodicity does not allow intermediate struc- tures with a lower density of active sites. Overall, only a relatively small number of distinct geometric motifs exist on each metal facet. Clearly, the creation of defects that give very active structures signif- icantly drives up the surface energy (Fig. 5). Yet, Fig. 2 indicates that even though the optimal structures possess a higher surface energy, numerous metastable structures of lower surface energy exhibit a nearly order of magnitude increase in activity (blue circles at large ω). This offers promise that a profound performance improvement

Fig. 2 | Current density and surface energy simulation results for numerous defected Pt(111)-based crystals that contain different active sites and/or density of active sites. a,b, Current densities (equation (2); a) and surface energies (equation (3); b) from the optimization are plotted against the activity weight ω. Red circles, values after the multiobjective optimization; blue circles, data after the energy-only minimization is
also applied.

a second-layer atom with a GCN of 8.00 or 8.17, respectively. The five-atom vacancy is similar to the structure in Fig. 4a and can be packed similarly on the surface. A maximal ordered packing of this site would result in activity 14 times higher than that of Pt(111). By comparison, our best structure is 46 times more active than Pt(111) (see below). We attribute the better performance of our best struc- ture to the slightly higher (about a factor of two) activity of our opti- mal active sites and the disordered packing of various types of active sites. Our exhaustive and systematic approach predicts more active sites than intuitive choices (that is, closer to the volcano peak) that can be packed densely on the surface. Additionally, our approach is general to other metals and structures.
The Pt(100)-based surfaces show more order. A four-atom vacancy that exposes a single active site in the second layer is the only active site found (Fig. 4d). The differing activities of Pt(100) surfaces are attributed only to the differing densities of this active site. The GCN of the active site is only 7.33, which falls short of the Pt peak. The (100) geometry is unable to possess a site with a higher GCN by introducing vacancies in the top layer.
Similar to Pt(111), Au(111) has multiple types of active sites. For the less-active structures, large cavities, as shown in Fig. 4e, are present. The active sites reside on the top layer edge sites along the cavity. All the atoms have coordination numbers of seven or higher and, therefore, surface stability is relatively maintained. The most-active structures contain a more active but less stable configuration of sites, as shown in Fig. 4f. Channels of vacan- cies expose edge sites on the top layer in a high density. Defected Au(100) has the periodic pattern of vacancies shown in Fig. 4g. As the optimal GCN is not very different from the ideal (100) sur- face, only a few vacancies are needed. For each of the active sites shown in Fig. 4, the predictive accuracy of our model is assessed with DFT, as detailed in Supplementary Section 2. Validation of the active sites allows for model refinement in an efficient hierar- chical refinement approach37–39.
In Fig. 5, current density and surface energy are plotted against each other for the Pareto optimal points. These structures span

is feasible via the defect engineering of experimental structures.
Catalyst activity and stability. The increase in current density due to the defects is summarized in Fig. 6. For Pt(111), our most-active structure consists of 49% vacancies on the top layer and is 46 times more active than an ideal Pt(111) surface. This is about an order of magnitude less active compared to the peak of the volcano curve due to the inability to create such an active site on Pt(111) and to geo- metric constraints, that is, the inability to have each site as the most active. In comparison, Calle-Vallejo et al.17 observed experimentally a 3.5-fold increase in activity due to defects. Our best structure is about an order of magnitude more active than the experimentally measured one. We attribute this difference either to an experi- mentally lower density of defects or to defects of lower activity; we cannot infer the precise reason without a detailed microscopic char- acterization of the experimental structures. However, we provided further insights by creating a random structure of defected Pt(111) with the same number of (49%) vacancies in the top layer as in our optimal structure. Then, we subjected the structure to annealing (by gradually decreasing the temperature in a Monte Carlo simulation) with a constant Pt loading of adatoms (Supplementary Fig. 8). This situation could relatively easily be realized in a real experiment in which a certain number of adatoms is added on Pt(111), for exam- ple, via atomic layer deposition, and the structure is then slightly annealed. The resulting average activity of the annealed structure is surprisingly comparable to (slightly lower than) the experimental value and double that of the perfect Pt(111) (Supplementary Fig. 8). Our results underline that the total density of vacancies or adatoms on the top layer is an important but not sufficient factor for improv- ing activity. Furthermore, we hypothesize that the experiments probably had a suboptimal distribution of defects. Microscopic engineering of the appropriate microgeometry of active sites and maximizing the density of such defects are crucial to materialize fully the benefit of crystal imperfection. Randomly created struc- tures with the correct density of vacancies or adatoms give only a modest boost in activity.
Other surfaces experience differing degrees of improvements when defects are introduced. Defects in Pt(100) increase the activ- ity by up to a factor of 15, which is threefold less than the enhance- ment on Pt(111) because the active site does not reach the volcano peak and therefore gains in current density are modest. Au(111)’s activity can be enhanced by almost three orders of magnitude due to how densely the active sites are packed along edges in the optimal structure. Au(100)’s activity, by contrast, benefits less from defects because it is already reasonably active relative to the volcano peak.

Fig. 3 | Structures of defected crystals. a–d, p(30 × 30) structure of Pt(111) (a), Pt(100) (b), Au(111) (c) and Au(100) (d) slabs with defects in the top layer, layer 1 (lighter shaded atoms). Atoms in layer 1 are added and removed during the optimization. Atoms in layers 2, 3 and 4 are held fixed. The structures shown are the highest activity structures for each surface achieved after the quenching step.

Fig. 4 | Top-down views of the active sites for each surface. Only the top and next layer, when its site is exposed, are shown. a–c, Pt(111) sites. d, Pt(100) site. e,f, Au(111) sites. g, Au(100) sites. Red atoms indicate the most-active sites, which occur in the second layer for Pt and the top layer for Au. Blue atoms are Pt and gold atoms are Au. Darker-shaded atoms indicate those in the bottom layers. The numbers in the bottom left of each panel indicate the GCN of the active site. Each surface has an optimal GCN that differs from the ideal one (peak of the volcano) and is determined from symmetry and the metal itself. The maximum activity requires a random packing of the sites to maximize the number of sites per unit area, and the optimal active site may change with an increased packing of the sites and thus with catalyst activity.

It is clear that the higher symmetry of the (111) plane enables a higher degree of engineering by introducing defects and therefore and activity improvement, compared to the (100) facet. For Pt, even the best structures fall short (by an order of magnitude) of reaching the peak activity of this catalyst. Pt cuboctahedral parti- cles with surface defects are highly desirable. In contrast, defected cubic particles of Au approach closely the optimal performance of this metal. Yet, Au is significantly less active than Pt even with the introduction of defects.

To assess the robustness of our predictions, we evaluated two potential sources of error. The first is the variability in the GCN correlations, that is, the imperfect fits to the DFT data. The uncer- tainty of the GCN model is propagated through the volcano model to produce a range of activity predictions for each GCN. The dif- ferences in activity between the optimal active sites and the ideal surfaces are beyond the statistical uncertainty of the GCN model, but the optimal coordination number (Supplementary Fig. 11) is relatively insensitive to this uncertainty. We found that the typical

Fig. 5 | Pareto front of current density and surface energy. Current density j (equation (2)) versus surface energy γ (equation (3)) of defected surfaces. a–d, The initial structures at the beginning of each simulation are Pt(111) (a), Pt(100) (c), Au(111) (b) and Au(100) (d). The black lines trace the data of the optimal structures via a linear spline of the plotted points and indicate the Pareto trade-off between stability and activity. Straight lines indicate that activity increases due to an increase in the density of the same active site; kinks in linearity are indicative of the creation of more than one active site. The most-active structures have the highest surface free energy and are thus less stable.

volcano curve without uncertainty overpredicts the maximum rate, which may explain why this is unattainable experimentally. Uncertainty and correlations reduce the expected maximum activ- ity by an order of magnitude. Supplementary Section 7 details the approach and results.
The second potential source of error is the effect of solvation and adsorbate interactions. To estimate these effects, a structure-depen- dent microkinetic model (MKM), which includes explicit solvation and coverage effects, was constructed (Supplementary Section 8). The structure optimization is repeated using the MKM. Although the predicted activity of the optimal structure changes quantita- tively, the catalyst structure in terms of defects remains essentially unchanged relative to the model without explicit solvation and cov- erage effects. Coverage effects have a much more pronounced effect on the optimum coordination number than differences that arise from the use of explicit or implicit solvation.
The proposed approach can be combined with ‘activity maps’33,40,41 to predict the structure of bimetallic catalysts42, and can easily be extended to multiple low Miller index planes of nanopar- ticles43–46 or entire nanoparticles via the Wulff construction34,35,47 and eventually discrete descriptions of the surface modelled with kinetic Monte Carlo simulations14,48–50 and selectivity problems51.
Conclusions
We have introduced a framework to enable atomistic level materials engineering. The framework predicts the optimal catalyst perfor- mance by simultaneously identifying both the optimal active site of a surface facet of a catalyst and the geometric arrangement of such sites that maximizes their density, while accounting for cata- lyst stability. Our approach offers insights unobtainable by intuition alone. We applied this framework to predict the optimal defects on Pt(111), Pt(100), Au(111) and Au(100) surfaces in catalysing the ORR chemistry of significance to fuel cells.
The results are rather revealing. Each material and facet has a dif- ferent optimal active site whose geometric structure is not intuitive. This is unlike volcano curves constructed in the past that assume

Fig. 6 | Activity of defected metastable structures compared to ideal crystals. The red bars show the activities of the most-active defected metastable surfaces (highest red in Fig. 5 and the structures shown in Fig. 3). The blue bars show the activities of the perfect, non-defected surfaces.
Numbers above the bars indicate the ratio between the defected and
non-defected surfaces. Green circles show the predicted activities for the optimal sites using DFT-computed binding energies on several platinum and gold sites for both the (111) and (100) facets. The platinum activities indicated by green circles are for sites with GCN values of 8.35, 8.33 and 7.33, whereas gold activities indicated by green circles are based on sites with GCN values of 5.67 and 5.75. Supplementary Section 2 gives the methods and a detailed discussion of the structures and corresponding energies outlined in Supplementary Fig. 3 and Supplementary Table 5, respectively. Black circles and error bars are the 90% confidence intervals around the median of the activities of the annealed structures (average of the blue circles in the high activity regimes of Fig. 2 and Supplementary Fig. 7), which are more realistic than the top activity of the bars. Error bars quantify the variability in how much activity is lost by destroying active sites due to processes such as diffusion. As a reference, Pt and Au peaks
correspond to the optimal site of the volcano curve assuming the density of such sites to be that of the (111) facet; such densities are not geometrically possible. Defect engineering can profoundly increase the activity over that of ideal crystal facets, but site geometry, density and stability should also be considered carefully.

the same active site for all materials. For materials on the left of the structure-based volcano curve (activity versus GCN), overcoordi- nated sites (vacancies) are needed, whereas on the right of the vol- cano curve, undercoordinated sites along the edges are needed. The active site on Pt(111) and Au(111) changes with increasing defect density due to packing constraints. Interestingly, the optimal facet structures contain a disordered ensemble of a small number of dis- tinct sites. Our results indicate that the optimal structures are sig- nificantly inferior to the peak of the volcano due to the lack of ideal active sites, packing constraints on low index planes and higher sur- face energies compared to the ideal facets. Yet, metastable defected structures can still be significantly better than perfect crystals. In contrast, randomly put structures give marginal activity improve- ments. Defect engineering, whereby active sites of the right geomet- ric motif are created and stabilized, is thus imperative to materialize significant performance enhancements.
For the studied ORR catalysts, both the maximum activity and median of activities of the annealed structures is ranked Pt(100) < Au(111) < Au(100) < Pt(111). Geometrically, all the surfaces other

than Pt(100) are able to produce nearly optimal sites, but the density of these sites limits the activity. The activity of Au improves much more than that of Pt because the low-coordinated active sites on Au are easier to produce in abundance than the high-coordinated active sites on Pt. Similarly, the (111) facets exhibit a higher degree of engineering than the (100) facets because the hexagonal sym- metry offers greater flexibility compared to the square symmetry. These specific results are system dependent. Our studies pertain to defected surfaces in both vacuum and in aqueous environments at 0.9 V relative to the standard hydrogen electrode. Both the elec- trolyte and voltage can alter the results. A nanoparticle study in an HClO4 environment52 agrees with our results that Pt(111) is more active for ORR than Pt(100), and studies using H2SO4 as an electro- lyte53 suggest that sulfates result in Pt(100) being more active for the reaction than Pt(111). Adding chloride ions to the H2SO4 solution makes Pt(111) a more favourable surface54. Although the methods provided are universal, parts of the model must be parameterized for specific reaction conditions.
The proposed structure optimization allows the identification of the active sites and of the optimal arrangement of such sites, as well as a comparison of different materials and facets. Uncertainty in the structure-scaling relation and correlations in the param- eters lower the expected activity by an order of magnitude but, interestingly, only close to the optimal coordination number. In contrast to solvation-based interactions and the error in the struc- ture-scaling relation, lateral interactions have a pronounced effect the Pt(110) surface61. Our approach allows one to account for the adsorbate effects on surface stability62 by considering the relative adsorption energies on different sites. The stability of the catalyst in the presence of adsorbates will be worth exploring further in dedicated future work.
Catalyst structure and optimization. The catalyst surface is modelled as a four- layer (111) or (100) slab, shown in Fig. 3, of a p(30 × 30) unit cell size and with periodicity in the x–y directions. The layers in the slab are indexed from top to bottom. The top layer (1) may include atoms or vacancies, that is, it is a defected layer, whereas deeper layers contain no defects (vacancies in deeper layers can easily be considered in our model, but their geometric accessibility should also be taken into account). Each microstructure is uniquely described by a vector c which lists the occupancies (1 if present and 0 if absent) of the metal atoms at each lattice position. Each exposed binding site contributes to the overall activity dependent on the atom’s GCN.
For each microstructure, the GCNs for all the atoms in the slab are computed rapidly using graph theory. Details are described in Supplementary Section 3.
Given the GCN of a site, the binding energies of *OH and *OOH are computed according to Supplementary equations (5) and (6), respectively, and are used to estimate the site-specific activity (equation (1)). The activity of the structure is computed by summing over each active site (equation (2)). The structure’s surface energy and total activity are non-dimensionalized to E and I, respectively, for numerical convenience, as detailed in Supplementary Section 3.
We performed a two-level optimization. First, we simultaneously minimized the surface energy and maximized the activity using a weighted sum for the objective function, as is standard in multiobjective optimization63,64. A weighted sum of objectives is capable of identifying any Pareto optimal point that does not lie below the line interpolating between the extrema of the Pareto front:
F = × I + (1− ) × E (5)
0 ≤ ω ≤1 is a weighting factor that emphasizes the surface energy minimization

on the optimum coordination number. Although the absolute catalyst activity is sensitive to such model uncertainties, the con- cepts of defect engineering are robust outcomes of the model. This capability can be extended to any chemistry to enable significant improvements of catalysts.

Methods
Coordination numbers as a descriptor of activity and stability. We compute the current i of an individual site as:

(ω = 0) or activity maximization (ω = 1). A simulated annealing Metropolis Monte Carlo algorithm was applied, at each event a position in the top layer of the slab was chosen and a move that changed the occupancy was attempted (an occupied site becomes empty and vice versa) according to the Metropolis algorithm using the change in the objective function ΔF (details in Supplementary Section 3).
Inverting the sign for the activity in the objective function caused our algorithm to create structures with an activity lower than that of an ideal surface, which demonstrates how specific types of defects, rather than arbitrary ones, are needed to achieve a higher activity. An example is shown in Supplementary Section 6.
Given that the multiobjective optimization strikes a balance between competing factors, the resulting structures can be unstable to small perturbations, especially when there is a large weight on activity. For this reason, after the multiobjective optimization, relaxation was performed to minimize the surface energy locally, whereby atoms in the top layer were allowed to diffuse using a downhill energy minimization approach. The structure converges to the nearest

ΔGORR is the change in Gibbs energy of the limiting reaction and kB is the Boltzmann constant. The constant ic = 3.68 × 1011 kA mol–1 is used to match the experimental value of the Pt(111) current density at a cell potential of U = 0.9 V and temperature T = 298 K reported by Calle-Vallejo et al.12. The current density j of a microstructure composed of many sites is then computed from the activities of all the surface sites normalized by the exposed surface area Asurf:

local minimum on the coarse-grained potential energy surface, typically at the expense of some loss in activity. The entire algorithm is summarized in a pseudocode in Supplementary Section 3.
For a given metal facet, simulations were run for ω values using a random initial structure and random seed for each one. The inherent randomness of the optimization resulted in some structures having superior properties to others (fluctuations in the blue circles in Fig. 2). Structures of high current density and low surface energy are the most interesting. Supplementary Section 4 describes

Asurf k
Supplementary equations (5) and (6) in conjunction with equation (1) describe the structure sensitivity of the current density.
The stability of an extended surface is approximated by its surface energy (γ), which is computed by dividing the formation energy (Eform) of a slab Asurf: = Eform

how a Pareto optimality criterion was used to select these structures.

Data availability
All code and data that went into generating the results in this paper can be found at https://github.com/VlachosGroup/ORR-Optimization. The GitHub repository includes code for the following: structure generation, MKM, DFT data used in parameterizing the Hamiltonian used in the MKM and optimization algorithms.

Code availability

The formation energy is rapidly calculated using a surrogate model derived from tight-binding theory55,56:

All code that went into generating the results in this paper can be found at https:// github.com/VlachosGroup/ORR-Optimization. All the code is implemented in Python and requires the atomic simulation environment to represent molecular
structures and for file input and output (IO). The directory structure is described in the readme.rst file.
The summation is taken over all atoms k, each with coordination number CNi. CNmax = 12 is the maximum coordination number in a face-centred cubic lattice. Ecoh is the cohesive energy of the metal. To validate equation (4), DFT calculations were performed as detailed in Supplementary Section 2. Although we are interested in the surface energy for catalyst design purposes, experiments show that adsorbates can restructure a surface57. High-coverage CO on the Pt(100) and Pt(110) surfaces can even lift reconstruction entirely58–60. Low-energy electron diffraction studies showed that oxygen, however, does not result in restructuring of

References
1. Van Santen, R. A. Complementary structure sensitive and insensitive catalytic relationships. Acc. Chem. Res. 42, 57–66 (2009).
2. Koper, M. T. M. Structure sensitivity and nanoscale effects in electrocatalysis.
Nanoscale 3, 2054–2073 (2011).

3. Carchini, G. et al. How theoretical simulations can address the structure and activity of nanoparticles. Top. Catal. 56, 1262–1272 (2013).
4. Roldan Cuenya, B. Metal nanoparticle catalysts beginning to shape-up.
Acc. Chem. Res. 46, 1682–1691 (2012).
5. Somorjai, G. A. & Park, J. Y. Colloid science of metal nanoparticle catalysts in 2D and 3D structures. Challenges of nucleation, growth, composition, particle shape, size control and their influence on activity and selectivity. Top. Catal. 49, 126–135 (2008).
6. Xia, X. et al. Facile synthesis of palladium right bipyramids and their use as seeds for overgrowth and as catalysts for formic acid oxidation. J. Am. Chem. Soc. 135, 15706–15709 (2013).
7. Kim, S. K. et al. Performance of shape-controlled Pd nanoparticles in the selective hydrogenation of acetylene. J. Catal. 306, 146–154 (2013).
8. Renzas, J. R., Zhang, Y., Huang, W. & Somorjai, G. A. Rhodium nanoparticle shape dependence in the reduction of NO by CO. Catal. Lett. 132,
317–322 (2009).
9. Bratlie, K. M., Lee, H., Komvopoulos, K., Yang, P. & Somorjai, G. A. Platinum nanoparticle shape effects on benzene hydrogenation selectivity. Nano Lett. 7, 3097–3101 (2007).
10. Lee, I., Delbecq, F., Morales, R., Albiter, M. A. & Zaera, F. Tuning selectivity in catalysis by controlling particle shape. Nat. Mater. 8, 132–138 (2009).
11. Mostafa, S. et al. Shape-dependent catalytic properties of Pt nanoparticles.
J. Am. Chem. Soc. 132, 15714–15719 (2010).
12. Calle-Vallejo, F. et al. Finding optimal surface sites on heterogeneous catalysts by counting nearest neighbors. Science 350, 185–189 (2015).
13. Dubau, L. et al. Defects do dcatalysis: CO monolayer oxidation and oxygen reduction reaction on hollow PtNi/C nanoparticles. ACS Catal. 6, 4673–4684 (2016).
14. Guo, W. & Vlachos, D. G. Patched bimetallic surfaces are active catalysts for ammonia decomposition. Nat. Commun. 6, 8619–8625 (2015).
15. Calle-Vallejo, F. et al. Why conclusions from platinum model surfaces do not necessarily lead to enhanced nanoparticle catalysts for the oxygen reduction reaction. Chem. Sci. 8, 2283–2289 (2017).
16. Calle-Vallejo, F., Martínez, J. I., García-Lastra, J. M., Sautet, P. & Loffreda, D.
Fast prediction of adsorption properties for platinum nanocatalysts with generalized coordination numbers. Angew. Chem. Int. Ed. 53, 8316–8319 (2014).
17. Calle-Vallejo, F., Loffreda, D., Koper, M. T. M. & Sautet, P. Introducing structural sensitivity into adsorption–energy scaling relations by means of coordination numbers. Nat. Chem. 7, 403–410 (2015).
18. Mpourmpakis, G., Andriotis, A. N. & Vlachos, D. G. Identification of descriptors for the CO interaction with metal nanoparticles. Nano Lett. 10, 1041–1045 (2010).
19. Hanselman, C. L. & Gounaris, C. E. A mathematical optimization framework for the design of nanopatterned surfaces. AIChE J. 62, 3250–3263 (2016).
20. Zhuang, H., Tkalych, A. J. & Carter, E. A. Surface energy as a descriptor of catalytic activity. J. Phys. Chem. C 120, 23698–23706 (2016).
21. Rong, X., Parolin, J. & Kolpak, A. M. A fundamental relationship between reaction mechanism and stability in metal oxide catalysts for oxygen evolution. ACS Catal. 6, 1153–1158 (2016).
22. Perez-Alonso, F. J. et al. The effect of size on the oxygen electroreduction activity of mass-selected platinum nanoparticles. Angew. Chem. Int. Ed. 51, 4641–4643 (2012).
23. Markovic, N., Gasteiger, H. & Ross, P. N. Kinetics of oxygen reduction on Pt(hkl) electrodes: implications for the crystallite size effect with supported Pt electrocatalysts. J. Electrochem. Soc. 144, 1591–1597 (1997).
24. Lee, S. W. et al. Role of surface steps of Pt nanoparticles on the electrochemical activity for oxygen reduction. J. Phys. Chem. Lett. 1, 1316–1320 (2010).
25. Guerin, S., Hayden, B. E., Pletcher, D., Rendall, M. E. & Suchsland, J.-P. A combinatorial approach to the study of particle size effects on supported electrocatalysts: oxygen reduction on gold. J. Comb. Chem. 8, 679–686 (2006).
26. Gasteiger, H. A., Kocha, S. S., Sompalli, B. & Wagner, F. T. Activity benchmarks and requirements for Pt, Pt-alloy, and non-Pt oxygen reduction catalysts for PEMFCs. Appl. Catal. B 56, 9–35 (2005).
27. Wang, B. Recent development of non-platinum catalysts for oxygen reduction reaction. J. Power Sources 152, 1–15 (2005).
28. Jinnouchi, R., Suzuki, K. K. T. & Morimoto, Y. DFT calculations on electro-oxidations and dissolutions of Pt and Pt–Au nanoparticles. Catal. Today 262, 100–109 (2016).
29. Stephens, I. E. L., Bondarenko, A. S., Grønbjerg, U., Rossmeisl, J. & Chorkendorff, I. Understanding the electrocatalysis of oxygen reduction on platinum and its alloys. Energy Environ. Sci. 5, 6744–6762 (2012).
30. Yamamoto, K. et al. Size-specific catalytic activity of platinum clusters enhances oxygen reduction reactions. Nat. Chem. 1, 397–402 (2009).
31. Nørskov, J. K. et al. Universality in heterogeneous catalysis. J. Catal. 209, 275–278 (2002).
32. Nørskov, J. K. et al. Origin of the overpotential for oxygen reduction at a fuel-cell cathode. J. Phys. Chem. B 108, 17886–17892 (2004).

33. Norskov, J. K., Bligaard, T., Rossmeisl, J. & Christensen, C. H. Towards the computational design of solid catalysts. Nat. Chem. 1, 37–46 (2009).
34. Tripković, V., Cerri, I., Bligaard, T. & Rossmeisl, J. The influence of
particle shape and size on the activity of platinum nanoparticles for oxygen reduction reaction: a density functional theory study. Catal. Lett. 144, 380–388 (2014).
35. Tritsaris, G. A., Greeley, J., Rossmeisl, J. & Nørskov, J. K. Atomic-scale modeling of particle size effects for the oxygen reduction reaction on Pt. Catal. Lett. 141, 909–913 (2011).
36. Li, H., Zhao, M. & Jiang, Q. Cohesive energy of clusters referenced by Wulff construction. J. Phys. Chem. C 113, 7594–7597 (2009).
37. Salciccioli, M., Stamatakis, M., Caratzoulas, S. & Vlachos, D. G. A review of multiscale modeling of metal-catalyzed reactions: mechanism development for complexity and emergent behavior. Chem. Eng. Sci. 66, 4319–4355 (2011).
38. Sutton, J. E., Guo, W., Katsoulakis, M. A. & Vlachos, D. G. Effects of correlated parameters and uncertainty in electronic-structure-based chemical kinetic modelling. Nat. Chem. 8, 331–337 (2016).
39. Ulissi, Z. W., Medford, A. J., Bligaard, T., Nørskov, J. K. & Nørskov, J. K. To address surface reaction network complexity using scaling relations machine learning and DFT calculations. Nat. Commun. 8, 14621 (2017).
40. Salciccioli, M. & Vlachos, D. G. Kinetic modeling of Pt catalyzed and computation-driven catalyst discovery for ethylene glycol decomposition. ACS Catal. 1, 1246–1256 (2011).
41. Hansgen, D. A., Vlachos, D. G. & Chen, J. G. Using first principles to predict bimetallic catalysts for the ammonia decomposition reaction. Nat. Chem. 2, 484–489 (2010).
42. Stephens, I. E. L. et al. Tuning the activity of Pt(111) for oxygen electroreduction by subsurface alloying. J. Am. Chem. Soc. 133, 5485–5491 (2011).
43. Herron, J. A., Scaranto, J., Ferrin, P., Li, S. & Mavrikakis, M. Trends in formic acid decomposition on model transition metal surfaces: a density functional theory study. ACS Catal. 4, 4434–4445 (2014).
44. Ferrin, P. & Mavrikakis, M. Structure sensitivity of methanol electrooxidation on transition metals. J. Am. Chem. Soc. 131, 14381–14389 (2009).
45. Falsig, H. et al. Trends in the catalytic CO oxidation activity of nanoparticles.
Angew. Chem. 120, 4913–4917 (2008).
46. Jiang, T. et al. Trends in CO oxidation rates for metal nanoparticles and close-packed, stepped, and kinked surfaces. J. Phys. Chem. C 113, 10548–10553 (2009).
47. Blaylock, D. W., Zhu, Y.-A. & Green, W. H. Computational investigation of the thermochemistry and kinetics of steam methane reforming over a multi-faceted nickel catalyst. Top. Catal. 54, 828–844 (2011).
48. Lin, S. et al. Influence of step defects on methanol decomposition: periodic density functional studies on Pd (211) and kinetic Monte Carlo simulations. Phys. Chem. C 117, 451–459 (2013).
49. Yang, L., Karim, A. & Muckerman, J. T. Density functional kinetic Monte Carlo simulation of water–gas shift reaction on Cu/ZnO. J. Phys. Chem. C 117, 3414–3425 (2013).
50. Stamatakis, M., Chen, Y. & Vlachos, D. G. First-principles-based kinetic Monte Carlo simulation of the structure sensitivity of the water gas shift reaction on platinum surfaces. J. Phys. Chem. C 115, 24750–24762 (2011).
51. Christopher, P. & Linic, S. Engineering selectivity in heterogeneous catalysis: Ag nanowires as selective ethylene epoxidation catalysts. J. Am. Chem. Soc. 130, 11264–11265 (2008).
52. Shao, M., Peles, A. & Shoemaker, K. Electrocatalysis on platinum nanoparticles: particle size effect on oxygen reduction reaction activity. Nano Lett. 11, 3714–3719 (2011).
53. Wang, C., Daimon, H., Onodera, T., Koda, T. & Sun, S. A general approach to the size- and shape-controlled synthesis of platinum
nanoparticles and their catalytic reduction of oxygen. Angew. Chem. Int. Ed.
47, 3588–3591 (2008).
54. Stamenkovic, V., M. Markovic, N. & Ross, P. N. Structure-relationships in electrocatalysis: oxygen reduction and hydrogen oxidation reactions on Pt(111) and Pt(100) in solutions containing chloride ions. J. Electroanal. Chem. 500, 44–51 (2001).
55. Tománek, D. et al. Simple theory for the electronic and atomic structure of small clusters. Phys. Rev. B 28, 665–673 (1983).
56. Eichler, A. et al. Structural and electronic properties of rhodium surfaces: an ab initio approach. Surf. Sci. 346, 300–321 (1996).
57. Somorjai, G. A. & Van Hove, M. A. Adsorbate-induced restructuring of surfaces. Prog. Surf. Sci. 30, 201–231 (1989).
58. Martin, R., Gardner, P. & Bradshaw, A. M. The adsorbate-induced removal of the Pt{100} surface reconstruction Part II: CO. Surf. Sci. 342, 69–84 (1995).
59. Schwegmann, S., Tappe, W. & Korte, U. Quantitative structure analysis of a disordered system: RHEED study of the CO induced (1 × 2) → (1 × 1) structure transition of Pt(110). Surf. Sci. 334, 55–76 (1995).
60. Karakatsani, S., Ge, Q., Gladys, M. J., Held, G. & King, D. A. Coverage- dependent molecular tilt of carbon monoxide chemisorbed on Pt{110}: a combined LEED and DFT structural analysis. Surf. Sci. 606, 383–393 (2012).
61. Walker, A. V., Klötzer, B. & King, D. A. Dynamics and kinetics of oxygen dissociative adsorption on Pt{110}(1×2). J. Chem. Phys. 109, 6879–6888 (1998).
62. Tománek, D. Simple criterion for the reconstruction of clean and adsorbate- covered metal surfaces. Phys. Lett. A 113, 445–448 (1986).
63. Czyzżak, P. & Jaszkiewicz, A. Pareto PT-100 simulated annealing—a metaheuristic technique for multiple‐objective combinatorial optimization. J. Multi‐Criteria Decis. Anal. 7, 34–47 (1998).
64. Marler, R. T. & Arora, J. S. Survey of multi-objective optimization methods for engineering. Struct. Multidiscip. Optim. 26, 369–395 (2004).

Acknowledgements
Material developed by Núñez and Vlachos is based on work supported by the US Department of Energy Office of Science, Office of Advanced Scientific Computing Research and Applied Mathematics programme under award no. DE-SC0010549. Funding for Lansford was provided by the Defense Advanced Research Project Agency under grant W911NF-15-2-0122. We thank F. Calle-Vallejo for providing the original data of his group’s work. We thank S. Giles for useful discussion on electrochemistry.

Author contributions
M.N. implemented the optimization methodology, performed the calculations and analysed the results. J.L.L. performed the DFT calculations and MKM. D.G.V. conceived the problem and provided supervision. All the authors contributed to writing the paper.

Competing interests
The authors declare no competing interests.

Additional information
Supplementary information is available for this paper at https://doi.org/10.1038/ s41557-019-0247-4.
Reprints and permissions information is available at www.nature.com/reprints.
Correspondence and requests for materials should be addressed to D.G.V.
Publisher’s note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
© The Author(s), under exclusive licence to Springer Nature Limited 2019