## Abstract

The rate of ammonia synthesis over a nanoparticle ruthenium catalyst can be calculated directly on the basis of a quantum chemical treatment of the problem using density functional theory. We compared the results to measured rates over a ruthenium catalyst supported on magnesium aluminum spinel. When the size distribution of ruthenium particles measured by transmission electron microscopy was used as the link between the catalyst material and the theoretical treatment, the calculated rate was within a factor of 3 to 20 of the experimental rate. This offers hope for computer-based methods in the search for catalysts.

Detailed theoretical descriptions of the way in which solid surfaces act as catalysts for chemical reactions are now being obtained from density functional theory (DFT) calculations, which can be used to obtain the relevant activation energies. For example, Linic and Barteau have shown that a mean-field kinetic model of the selective oxidation of ethylene on an Ag catalyst, developed on the basis of DFT calculations, can describe experimental data (*1*). However, the mean-field description implicitly neglects the complexity associated with interactions between adsorbed surface species and the resulting existence of many different possible reaction paths. This problem was overcome recently by Reuter *et al*., in a DFT-based kinetic Monte Carlo description of the oxidation of carbon monoxide over a RuO_{2}(110) surface (*2*).

Here we take the further step of developing a kinetic description that includes the full complexity of interactions and reaction paths for a complete catalytic reaction under industrial conditions over a packed bed of a high-surface-area nanoparticle catalyst. Using the ammonia (NH_{3}) synthesis as the example, we show that DFT calculations can be used to directly predict a reaction rate for a supported nanoparticle Ru catalyst that is in good agreement with rate measurements performed over a wide range of industrially relevant synthesis conditions. The only experimental input included was the particle size distribution for the Ru catalyst, which was determined from transmission electron microscopy (TEM).

The synthesis of NH_{3} is probably the most studied reaction in heterogeneous catalysis, and it acts as the prototype reaction that has been used to develop many key concepts in the field (*3*). The best elementary metal catalysts (Ru and Fe) were discovered in large-scale screening experiments almost 100 years ago (*4*–*6*), and the nature of the rate-determining step for Fe-based catalysts, N_{2} dissociation, was pinpointed as early as 1934 (*7*, *8*). About 25 years ago, surface science studies became possible and revealed a detailed picture of the N_{2} dissociation process (*9*–*13*). It has been shown experimentally (*14*–*16*) and theoretically (*17*–*19*) that there is a direct link between the results of the ultra-low-pressure surface science results and NH_{3} synthesis data at elevated pressure and temperature. Most recently, DFT calculations were used to quantitatively outline the complete reaction mechanism with all elementary steps on Ru (*20*). It has also been shown that it is possible to predict and understand the trends in reactivity when DFT calculations of activation barriers and stabilities of intermediates for a series of different catalysts are combined with a simple microkinetic model (*21*). Our work here represents the final step in a complete first-principles description of the Ru-catalyzed NH_{3} synthesis.

The starting point for our calculation is the potential energy diagram for the full reaction (*20*) (Fig. 1). The calculations (*22*) show that step sites are much more reactive for N_{2} dissociation than the close-packed (001) surface, a finding that has been verified experimentally (*23*). By combining the results for step sites in Fig. 1 with harmonic transition state theory, we have shown that N_{2} dissociation is by far the slowest step under all realistic reaction conditions (*20*). We will exploit these findings by treating N_{2} dissociation as rate-limiting, and we only consider dissociation along step sites where the active (B_{5}) sites exist (*23*).

The results shown in Fig. 1 are for a low coverage of adsorbates. In reality, there will be other atoms and molecules adsorbed in the vicinity of the reacting molecule, and the activation energy for dissociation, *E*_{a,i}, and thus the rate constant, *k*_{i} = *ve*^{–Ea,i/kT} where *v* is the prefactor, *k* is the Boltzmann constant, and *T* is the temperature, will depend on the local environment (Fig. 2). Each local environment *i* will contribute to the total rate with a weight given by the probability *P*_{i} of finding this environment. The total rate *r* can be written as where *K*_{g} is the gas-phase equilibrium constant, and *p*_{N2}, *p*_{H2}, and *p*_{NH3} are the pressures of N_{2}, H_{2}, and NH_{3}, respectively. The term ensures that the gas phase equilibrium is established.

The probability *P*_{i} is given by the equilibrium between H_{2} and NH_{3} in the gas phase and adsorbed H, N, NH, NH_{2}, and NH_{3}; all of the steps after the rate-limiting N_{2} dissociation (Fig. 1) are in equilibrium. For this reason, we calculated the free energy of adsorption for all of the adsorbates. The entropy contributions were calculated on the basis of a complete normal mode analysis with the harmonic approximation (table S2) (*22*). We also calculated all possible nearest-neighbor interactions along the step (table S3). In principle, all of the interactions on the lower step should also be included. However, because we found that only H atoms appeared on the lower step, only the H-H interaction was needed. The use of grand canonical Monte Carlo simulations provides ensembles of the system in equilibrium and thus the probabilities *P*_{i}. Because N_{2} dissociates with one N atom at the top of the step and another at the bottom, an environment *i* is defined by having an empty site at the upper and lower step sites simultaneously and by the occupancy of the four neighboring sites.

In order to test the predictions of the rate for a real catalyst under industrially relevant conditions, we measured the activity for a Ru catalyst supported on magnesium aluminum spinel that contained 11.1 weight percent (wt %) Ru (*24*). The experimental results are in the form of the NH_{3} productivity from a plug flow reactor loaded with a catalyst containing 0.2 g of the 11.1 wt % Ru/MgAl_{2}O_{4} catalyst as a function of the input partial pressures, the flow, and the temperature. We integrated the NH_{3} synthesis rate down through the reactor in order to obtain the calculated NH_{3} productivity.

The only information about the catalyst we needed in order to compare calculations to experiments is the number of active sites per gram of catalyst. Thus, we characterized the activated catalyst by TEM, and we obtained the particle size distribution by analyzing ∼10^{3} particles. The active-site density is not directly given by the particle size distribution, because only the B_{5} sites are active for N_{2} dissociation. Steps can be clearly observed in the TEM image of the Ru particle (Fig. 3), but it is not possible to count them. To estimate the fraction of B_{5} sites as a function of particle size, we calculated the surface energy of all the low-energy facets of Ru and used them in a Wulff construction to give the basic particle shape. From DFT calculations, we found that the (001)/(101) edges lower their energy by a reconstruction in which the edge row of atoms is removed. This process gives rise to the steps that contain B_{5} sites along the edge. Our estimate of the particle morphology is quite crude, and it neglects the interaction with the support and changes in the surface energies caused by adsorption. However, the shape is similar to that found in the TEM image, and the fraction of B_{5} (step) sites was insensitive to such details. The lack of active sites on the smallest Ru particles is in full agreement with an experimental observation that the activity of a Ru catalyst can increase as a result of sintering the smallest Ru particles in a supported catalyst (*25*). The active-site density could now be calculated, and by folding it with the measured size distribution, we obtained the number of active sites per Ru mass. The largest uncertainty in this estimate is the implicit assumption that the Wulff polyhedra are complete, i.e., that there are no partly filled layers. Our estimate is likely a lower bound to the number of active sites.

A direct comparison of the calculated and measured NH_{3} productivity at a number of different conditions is shown in Fig. 4A. Given that there are no fitted parameters, the agreement is excellent. The overall rate is too small by a factor of ∼3 to 20, and there is a tendency for the calculated inhibition by NH_{3} to be slightly too weak.

The calculations provide insight into the exact nature of the active sites for the NH_{3} synthesis reaction over Ru. In Fig. 4B, we give results that show the local environments contributing to the synthesis in one specific slice of the reactor. The largest contribution to the total rate does not come from the local environment with the lowest activation energy (Fig. 2). The state of the surface during synthesis conditions is such that the configuration with adsorbed H on the upper step dominates the total rate.

The agreement between theory and experiment in Fig. 4A might seem surprising, considering that the inherent accuracy of the DFT calculations is of the order 0.2 to 0.3 eV (*26*). An error in an activation energy of 0.25 eV, for instance, gives rise to an error in the rate of an elementary reaction of a factor of 148 at 600 K. The important point is that the total reaction rate is considerably less sensitive to the absolute error than the rate of the individual steps. We tested the sensitivity of the results to the main approximation in the DFT calculations, the exchange correlation energy functional. All of the calculations were done with the RPBE functional (*26*), but we also calculated all of the parameters using another functional, PW91 (*27*). We then interpolated all of the calculated energies between the two: *E*(*x*) = *xE*_{RPBE} + (1 – *x*)*E*_{PW91}. The plot of the rate as a function of the mixing parameter *x* (Fig. 4C) shows a very weak dependence, because there is a compensation effect between the different steps in the total reaction (*28*). Changing from RPBE to PW91 decreased the barrier for N_{2} dissociation substantially (by 0.6 eV) and made the dissociation of N_{2} much faster. The dissociative chemisorption energy of H_{2} changed simultaneously from –0.36 (RPBE) to –0.52 (PW91) eV, and the N_{2} chemisorption energy changed from –0.8 (RPBE) to –1.4 (PW91) eV. These changes increased the coverage (through the equilibrium with H_{2} and NH_{3} in the gas phase) and decreased the number of free sites for dissociation. Because the barrier for dissociation and the stability of the intermediates on the surface vary together, a change in the overall interaction strength has only a small effect on the net rate. It is this compensation between two opposite trends that gives rise to a volcano in the catalytic activity as a function of the bond strength (*28*). The compensation effect is largest near to the maximum of the volcano (that is, for the best catalysts).

The insensitivity of the total catalytic rate to errors in the overall adsorption energies raises the question of the origin of the remaining discrepancy between the measured and calculated results (Fig. 4A). One possibility is that the number of active sites is underestimated. Another possible source of discrepancy could be systematic errors in the description of the difference in bonding of different adsorbates or configurations. We tested the sensitivity of the results to such relative errors by decreasing the stability of adsorbed H relative to the NH_{x} species by 0.06 eV. It turns out that this is enough to make the calculated points match the measured ones completely at all temperatures and flows. Small relative errors are therefore a likely source of the underestimate of the rate in Fig. 4A. By the same token, the ability of the calculations to predict the absolute rate to within a factor of 3 to 20 suggests that the relative errors are only of the order 0.06 eV (∼*kT* at reaction conditions) for the exchange correlation functional used here.

Colinear variations in activation energies and bond strengths are found quite generally (*21*), so we would expect compensation effects to be found for other reactions. This built-in insensitivity to absolute errors, together with the higher accuracy of the DFT methods for relative energies, offers hope that DFT calculations can give a good overall description of the catalytic activity of other reactions as well.

**Supporting Online Material**

www.sciencemag.org/cgi/content/full/307/5709/555/DC1

Materials and Methods

Figs. S1 to S6

Tables S1 to S4

References and Notes