Applications of molecular modeling in heterogeneous catalysis research

Applied Catalysis A: General 200 (2000) 23–46

Review

Applications of molecular modeling in heterogeneous catalysis research Linda J. Broadbelt1 , Randall Q. Snurr∗ Department of Chemical Engineering and Institute for Environmental Catalysis, Northwestern University, Evanston, IL 60208, USA Received 16 February 2000; received in revised form 16 March 2000; accepted 16 March 2000

Abstract The application of molecular modeling in heterogeneous catalysis research as a complement to experimental studies has grown rapidly in recent years. This review summarizes methodologies for probing catalytic phenomena in terms of a hierarchical approach. The elements of the hierarchy are different computational methods at different time and length scales that may be linked together to answer questions spanning from the atomic to the macroscopic. At the most detailed level of description, quantum chemical calculations are used to predict the energies, electronic structures, and spectroscopic properties of small arrangements of atoms and even periodic structures. Atomistic simulations, using systems of hundreds or thousands of molecules, can be used to predict macroscopic thermodynamic and transport properties, as well as preferred molecular geometries. At the longest time and length scales, continuum engineering modeling approaches such as microkinetic modeling are used to calculate reaction rates, reactant conversion, and product yields and selectivities, using model parameters predicted by the other levels of the hierarchy. We highlight some interesting recent results for each of these approaches, stress the need for integrating modeling at widely varying time and length scales, and discuss current challenges and areas for future development. © 2000 Elsevier Science B.V. All rights reserved. Keywords: Quantum chemistry; Atomistic simulations; Microkinetic modeling; Molecular modeling; Hierarchical modeling; Adsorption; Diffusion; Mechanistic modeling

1. Introduction Traditionally, heterogeneous catalysis has largely been an experimental field. While this is still true today, it is increasingly recognized that computer mod∗ Corresponding author. Tel.: +1-847-467-2977; fax: +1-847-467-1018. E-mail addresses: [email protected] (L.J. Broadbelt), [email protected] (R.Q. Snurr). 1 Also corresponding author. Tel.: +1-847-491-5351; fax: +7-847-491-3728.

eling and simulation are important tools in the study and development of catalytic systems [1–6]. This is due to both recent theoretical developments and significant advances in computational power. Computer modeling has the potential to provide new insights into reaction pathways, to predict properties of catalysts that have not yet been synthesized, and to bring information for a given system from many different experimental techniques into a coherent picture. A close interaction of modeling and experiment is vital. Firstly, to achieve confidence in new modeling techniques, they must be carefully validated against experiment. In

0926-860X/00/$ – see front matter © 2000 Elsevier Science B.V. All rights reserved. PII: S 0 9 2 6 - 8 6 0 X ( 0 0 ) 0 0 6 4 8 - 7

24

L.J. Broadbelt, R.Q. Snurr / Applied Catalysis A: General 200 (2000) 23–46

doing this, it is useful to have comparisons with many different experiments, e.g. kinetic studies of reaction rates, thermodynamic information on adsorption, and spectroscopic data on molecular-level structure. Secondly, the most fruitful strategy for using modeling in catalysis research is likely to be a dual-feedback mode, where experiments are used to validate the modeling and modeling is used to explain experimental results, to suggest new experiments, and perhaps to substitute for experiments in the screening of different catalysts or reaction conditions. Computational modeling of catalysis is an exploding field. Widely available and easy-to-use quantum chemical packages [7–9] have made such calculations accessible to non-experts, and calculations that were done on supercomputers a few years ago are now done routinely with desktop personal computers. In recognition of computational quantum chemistry’s importance, Pople and Kohn received the 1998 Nobel Prize in Chemistry for their work in laying much of its groundwork. The 1999 Nobel Chemistry Prize to Ahmed Zewail for femtosecond spectroscopy of reactions also highlights that catalytic scientists increasingly have the tools to probe, understand, and control catalytic systems at the molecular level. The need for such molecular-level understanding is greater today than ever before, with industry under increasing economic and environmental pressure to develop more precise control of catalytic processes. Breakthrough advances in reaction selectivities and activities will only occur if we can understand and control catalytic systems at the molecular level. Catalysis depends on phenomena that span a wide range of time and length scales, and this makes modeling of catalytic systems a considerable challenge. Consider a full catalytic cycle, for example in a zeolite. First, reactant molecules must adsorb into the pores (physisorption). Then, they must diffuse through the pores and adsorb onto a reaction center (chemisorption). At this site, the molecule may react, perhaps requiring the presence of additional reactant molecules. Product molecules diffuse back out of the pores and desorb into the surrounding fluid phase to complete the cycle. These events occur on length scales ranging from the catalytic site (10−10 m) up to the size of the reactor (1 m). Relevant time scales may span from femtoseconds up to hours. In principle, one can predict the macroscopic phenomena of interest by

solving the governing microscopic equations, namely the Schrödinger equation of quantum mechanics. However, in practice, it will never be feasible to do this directly. Quantum mechanical calculations, even with today’s fastest computers, can only be carried out for hundreds of atoms, and the time scales of interest are many, many orders of magnitude beyond current capabilities. Spanning these gaps can only be accomplished by adopting a hierarchical approach. The strategy is to utilize different computational methods at different time and length scales and link them together to answer questions spanning from the atomic to the macroscopic. This kind of comprehensive approach is now emerging through the efforts of researchers around the world. The approach is summarized in Fig. 1. At the most detailed level of description, electronic structure calculations are used to predict the energetics of small arrangements of atoms. These results can provide energetic parameters that are fed into atomistic simulations. Monte Carlo and molecular dynamics (MD) simulations, using systems of hundreds or thousands of molecules, can be used to predict macroscopic thermodynamic and transport properties, as well as preferred molecular geometries. Preferred structures can then, in turn, be analyzed in more detail with the electronic structure methods. If the time and length scales are beyond those accessible to atomistic simulation, other statistical mechanical treatments are applied. For example, transition-state theory (TST) can be used to predict rates of chemical reactions or rates of slow physical processes. At the longest time and length scales, continuum engineering modeling approaches such as microkinetic modeling are used to calculate properties of interest, using model parameters predicted by other levels of the hierarchy. Note that different kinds of information, at different levels of detail, are obtained as the process proceeds. Feedback between different levels of the hierarchy is important in directing and refining the modeling performed at each level. The objective of this article is not to provide an exhaustive review of computational modeling in catalysis. The field is too large and growing too rapidly to make this practical. Instead, we highlight some interesting recent results, stress the need for modeling at widely varying time and length scales, and discuss areas for future development. Fig. 1 serves as an

L.J. Broadbelt, R.Q. Snurr / Applied Catalysis A: General 200 (2000) 23–46

25

Fig. 1. Hierarchical approach to modeling catalytic systems. Computational methods are shown in rectangles, with inputs and outputs in circles. General categories of computational methods are in bold (e.g. Quantum Mechanics), and specific examples are in regular lettering (e.g. Hartree–Fock). The hierarchical approach integrates together information from computational techniques at various length and time scales. Feedback between different levels of the hierarchy is important in directing and refining the modeling performed at each level.

outline of the topics to be discussed. We will begin at the smallest length scale and discuss quantum chemical calculations. We will next focus on obtaining information about adsorption and diffusion from atomistic simulations. We will then move to the macroscopic length scale and discuss microkinetic modeling. Finally, we will offer some concluding thoughts about future challenges. We hope to create more crossfertilization among researchers working in the different kinds of modeling shown in Fig. 1, so that, e.g. catalysis researchers using quantum chemical calculations may learn more about atomistic simulation and microkinetic modeling.

2. Quantum chemical calculations Quantum chemical calculations provide information about the smallest length scale of interest in the hierarchical approach to modeling heterogeneous catalytic processes. Based on solution of the Schrödinger equation, results from quantum mechanics may be used to predict molecular properties ranging from geometries to energetics to the distribution of electron density. In some cases, quantum chemical calculations offer the opportunity to probe details of catalytic chemistry that are difficult to obtain from experiment. For example, activation barriers for individual elementary steps

26

L.J. Broadbelt, R.Q. Snurr / Applied Catalysis A: General 200 (2000) 23–46

on surfaces can be calculated. Furthermore, detailed information about presumed active sites for catalysis may be obtained because the system is studied with explicit chemical detail. Various methods have been advanced to solve the Schrödinger equation, ranging from semiempirical methods to ab initio techniques. The basic principle that dictates the choice of a level of theory is the competition between speed and accuracy. Semiempirical methods are significantly less computationally demanding than techniques based on first principles. These methods introduce approximations to facilitate evaluation of terms introduced by electron–electron interactions and rely on parameterization against experimental data. The semiempirical method that has been used most widely for catalytic systems containing transition metals is ZINDO, which can also be useful for calculating electronic spectra [10]. As computational power has surged, techniques based on first principles, or ab initio methods, have been used more frequently. Examples of ab initio techniques include Hartree–Fock (HF), methods that more adequately represent electron correlation such as configuration interaction (CI) or Møller–Plesset perturbation methods (e.g. MP2), and density functional theory (DFT). A detailed description of the solution methodologies is beyond the scope of this review article, but complete details are given in several sources, e.g. [10,11]. In addition, excellent reviews of the application of quantum chemical methods to heterogeneous catalysis have been published [3–5]. Since the publication of these reviews, continued progress has been made in applying theoretical techniques to more realistic models of catalyst surfaces. We will therefore highlight what approaches are being advanced to model the actual size and composition of catalysts more accurately. In particular, we will focus on the use of cluster models, embedding schemes, and slab methods. As illustrated in Fig. 2, cluster models use a small number of atoms to represent the catalyst and often only capture the local environment around the active site. Embedding schemes treat a small region of the catalyst quantum mechanically and use classical methods to incorporate long-range interactions. Slab methods use periodic boundary conditions to treat the infinite solid. Many of the first ab initio calculations of interactions of adsorbates with surfaces used finite cluster

models because computational resources limited the number of atoms that could be included [5]. As significant gains in computational power have continued to become available at a modest expense, fewer cluster calculations are reported, particularly for metallic systems. However, in some cases, cluster calculations are still the appropriate choice. With improvements in synthetic capabilities, experimental studies of reactions over small, monodisperse metal clusters are increasingly being conducted, and these catalytic nanoclusters often have properties very different from those of the bulk metal. For example, nanometer-sized gold particles show catalytic activity for the oxidation of CO even though gold is conventionally thought to be an inert material [12]. Recent experiments by Heiz et al. [13] reveal the ability to tune the reaction rate of carbon monoxide oxidation on monodisperse platinum clusters by varying the cluster size from eight to 20 atoms, which is within the size range accessible to quantum chemical methods. The ruthenium clusters found in zeolites used for ammonia synthesis, on the order of 10 Å in diameter, can also be modeled directly [14]. To examine these systems using calculations, clusters are clearly the correct choice since they best mimic the experimental system. For modeling zeolite catalysis, clusters are also still widely used. For example, the unit cells of ZSM-5 and zeolite-␤ are comprised of 288 and 192 atoms, respectively, and therefore, only a representative portion of the framework that encompasses the local environment of the catalytic site is used in typical calculations. However, the size of the clusters extracted from the zeolite framework being studied is quickly growing. A recent paper by Sauer and coworkers reported the analysis of a 163-atom cluster using parallelized code [15]. One way to include longer-range effects without treating the whole system with ab initio quantum mechanics is to use embedded cluster calculations, in which a smaller number of atoms is treated quantum mechanically and the surrounding region is treated classically or at a lower level of theory. A large number of different implementations of these methods for different types of systems is available, and Gao [16] offers a recent review of the different approaches. One method that has been widely applied to investigate adsorption on metallic crystal faces is the many-electron embedding theory [17–19]. The local interactions between the adsorbate and a small

L.J. Broadbelt, R.Q. Snurr / Applied Catalysis A: General 200 (2000) 23–46

27

Fig. 2. Quantum chemical calculations of heterogeneous catalytic systems are typically carried out on one of three different representations of the catalyst: clusters, embedded clusters or slabs. (a) Nitrogen adsorption on ruthenium is studied using a small cluster of metal atoms to represent the extended metal surface [198]; (b) embedded clusters treat shaded central atoms quantum mechanically and atoms further away from the active site (shown as white atoms) classically (adapted from Fig. 1 of [20]); (c) slab calculations carried out by Klinke and coworkers [30–32] used a finite number of layers of metal atoms to represent the Ni(1 1 1) surface, placed adsorbates on the surface layer and replicated the unit cell in two dimensions.

28

L.J. Broadbelt, R.Q. Snurr / Applied Catalysis A: General 200 (2000) 23–46

cluster of surface and bulk atoms are rigorously optimized. The resulting system is then embedded within larger clusters to include bulk interactions. Recently, Yang and Whitten [20] investigated the chemisorption of OCN on Ni(1 0 0), using a three-layer, 64-atom cluster that embedding reduced to a 30-atom cluster. Although there were no experimental data to which the calculations could be compared, the results were in the range for other published experimental data for similar systems. van Santen and Neurock [5] caution that the approach is successful in modeling particular systems but does introduce some sense of arbitrariness that makes it difficult to adapt in a general sense. A paper published by Brändle et al. [21] is one example of a small number of cases in which it was possible to compare results for a large system using both an embedding scheme and the full quantum mechanical solution, using the same method and applying the same basis set. They examined proton siting and ammonia adsorption in chabazite, a zeolite with a small unit cell, employing a hybrid method, combined quantum mechanics/interatomic potential function (QM-Pot). They found that the combined QM-Pot relative stabilities and reaction energies deviated from the full QM results by only 4–9 kJ mol−1 . Furthermore, the electrostatic potentials calculated inside the zeolite pore by both methods were very similar, further justifying the use of the less computationally demanding approach. Embedded cluster calculations were a first step in being able to compare theoretical results with the rich database offered by surface science experiments. To enable even more direct comparison with experimental results, ab initio methods that treat the infinite solid with periodic boundary conditions are being used with increasing frequency. The methods are not simply a check on available experimental data but also offer new insight into unresolved questions. To highlight the versatility of these methods and their application to heterogeneous catalysis, several recent examples will be described. Many different systems have been examined by the Nørskov group, using periodic calculations, e.g. [22–25]. One system that they have examined in detail is the adsorption and dissociation of nitrogen on iron single-crystal surfaces [24,25]. DFT was used to examine differences of nitrogen adsorption on Fe(1 1 1), (1 0 0) and (1 1 0) surfaces. Periodic arrays of slabs of the different orientations were used, and self-consistent spin-polarized DFT calculations

with the generalized gradient approximation (GGA) were carried out. The pseudo-wavefunctions were represented using plane waves. The calculations revealed the stability of a c(2×2)-N/Fe(1 0 0) structure with the N atoms in the four-fold sites, consistent with the tendency for c(2×2) island formation found experimentally [26]. This study also reflects the ability of periodic calculations to examine the selectivity of different crystals faces for the same chemical process. To link these results to ammonia synthesis, the dissociation of nitrogen was studied in detail on Fe(1 1 1) [25]. Based on the calculations, a detailed picture of the dissociation process which is consistent with all experimental evidence was obtained. In particular, a new molecular precursor to dissociation accessible via two different routes was identified. The isolation of this species and the identification of two dissociation channels helped to resolve the question as to why experimental results disagreed about the presence of an activation barrier. In addition to the ability to examine different singlecrystal facets, periodic calculations may also be used to investigate bimetallic systems. Supported bimetallic catalysts are used by many processes in the petrochemical and fine chemicals industries [27]. Bimetallic systems have been investigated by Pallassana and coworkers who examined hydrogen chemisorption on Pd(1 1 1), Re(0 0 0 1), PdML /Re(0 0 0 1) and ReML / Pd(1 1 1), where the ML subscript denotes a single layer of one metal on top of the second as a bulk metal [28]. Gradient-corrected periodic DFT slab calculations were carried out using a planewave basis set. Calculations showed that hydrogen binding was strongest on the Re surfaces, and the weak binding energy observed for the Pd monolayer on Re(0 0 0 1) was consistent with experimental observations. Another advantage that periodic calculations offer is the ability to probe surface coverage effects and the extent of cluster edge effects. Pallassana et al. [29] examined the effect of surface coverage on the modes and energies of maleic anhydride adsorption. They showed that, at higher surface coverages, both the atop mode (less favorable at low surface coverages) and the di-␴ state were energetically similar. Klinke and coworkers [30–32] examined the dependence on surface coverage of the binding energy of small adsorbates relevant to Fischer–Tropsch synthesis on nickel and cobalt surfaces. Significant repulsion among carbon atoms was

L.J. Broadbelt, R.Q. Snurr / Applied Catalysis A: General 200 (2000) 23–46

29

Table 1 Effect of cluster size on the di-␴ adsorption energy of maleic anhydride on palladium (adapted from [29]) Total number of atoms in cluster

Number of atoms in first layer

Number of atoms in second layer

Binding energy (kJ mol−1 )

2 3 10 19 Pd(1 1 1) (θ=0.11 ML)

2 3 10 12 –

0 0 0 7 –

−76 −103 −87 −83 −84

observed as the coverage was increased to a monolayer, and as the surface coverage was increased further, stable graphite was formed. The dependence of the binding energy of hydrogen on these surfaces was minor as the coverage was varied. Pallassana et al. [29] also compared their results from periodic slab calculations with binding energies obtained using different cluster sizes. As summarized in Table 1, the dependence of the binding energy on Pd cluster size was significant for clusters with fewer than 10 atoms, but as the total cluster size was increased to 19 atoms, the binding energy was in excellent agreement with the results of the periodic slab calculations. These comparisons helped to verify that cluster edge effects were negligible for the 19-atom cluster and that the results are quantitatively reliable. Another method for computing ground-state electronic properties of large and/or disordered systems is the Car–Parrinello scheme, which combines MD with density functional calculations [10,33]. Boero et al. [34] recently applied this approach to simulate the polymerization of ethylene in a realistic Ziegler–Natta heterogeneous system. Both the deposition of the catalyst TiCl4 on the (1 1 0) active surface of a solid MgCl2 support and the polymer chain formation were simulated. The authors note that these simulations offered a way to study reaction mechanisms that were not currently accessible to experimental probes. This method of linking electronic structure calculations with the ensemble-based methods of statistical mechanics is an important advance. 3. Atomistic simulations In addition to the events at the active site, physical adsorption and diffusion can be important in a full catalytic cycle. These phenomena, however, occur on

longer time- and length-scales than those currently accessible to the detailed quantum chemical methods illustrated in the previous section. Also, they are generally governed by weak van der Waals forces, which are difficult to capture accurately with current quantum mechanical software. Thus, researchers turn to other computational tools; we focus here primarily on simulation methods using atomistic models. The properties of most interest are the adsorption isotherms, heats of sorption, diffusion coefficients, and activation energies for diffusion. In addition, the calculations should provide details of the microscopic structure and short-time motions that may be helpful in building our understanding of the systems of interest. Because adsorption thermodynamics and diffusion are especially important for catalysis in zeolites, much of the work in this area has focused on zeolites and other molecular sieves, and we will also focus on these systems for our discussion. The first ‘input’ to a molecular simulation is knowledge of the chemical constitution and structure of the zeolite and the adsorbate molecules. For many zeolites, the structures have been determined by X-ray diffraction studies [35], and the structures of small molecules are well known. 2 The zeolite framework is often assumed to be rigid. The other input required is a set of simple potential functions that describe the interaction energetics between adsorbates and the zeolite and among the adsorbate species themselves. Dispersion, repulsion, and electrostatic forces are typically accounted for, as well as intramolecular forces such as torsion angle energies for flexible sorbates [37]. Induced dipole and other forces may also be included if they are thought to be important. The total potential 2 For new zeolites where the structure is unknown, modeling can also be useful in elucidating the structure [36].

30

L.J. Broadbelt, R.Q. Snurr / Applied Catalysis A: General 200 (2000) 23–46

energy is usually approximated by a pairwise-additive potential summed over all pairs of atoms in the system [38]. Lennard–Jones parameters for the dispersion and repulsion energies can often be taken from previous studies in the literature or estimated from empirical expressions such as the Slater–Kirkwood formula [39], using physical properties of the atoms such as their polarizabilities and radii. Coulombic interactions are usually modeled by placing point charges on the atoms, with the partial charges usually obtained from quantum chemical calculations. The ‘outputs’ of a molecular simulation include thermodynamic properties, such as adsorption isotherms and heats of sorption; dynamic properties, such as the self-diffusion coefficient; and information on the molecular structure. Given the model, a molecular simulation relies on statistical mechanics to generate the desired information. In general, different simulation methods are required for predicting adsorption thermodynamics and diffusion. These are described in detail elsewhere [10,40,41], including several references that focus on microporous materials [42–44]. Brief descriptions are given below, along with examples from the literature related to catalysis. 3.1. Adsorption The simplest calculations of adsorption are in the low-loading Henry’s law region. The Henry’s constant, KH , can be calculated in milligrams sorbate per gram zeolite per atmosphere from the molecular model by evaluation of the following configurational integral [45,46]: KH =

Ma 2 8π Vs ρz RT R

×

exp[−βV z (r, Ψ, φ) − βV intra (φ)] dr dΨ dφ R exp[−βV intra (φ)] dφ (1)

where Ma is the sorbate molecular weight, Vs the zeolite volume, ρ z the bulk zeolite density, R the gas constant, T the temperature and where β is 1/kT with k being Boltzmann’s constant. The potential energy due to zeolite–sorbate interactions, Vz , is a function of the molecule’s position r, its orientation Ψ , and its internal degrees of freedom φ. The denominator is only integrated over the intramolecular degrees of

freedom. The integrals in Eq. (1) are easily evaluated using Monte Carlo integration or other numerical integration methods [47]. Complete isotherms are typically computed using grand canonical Monte Carlo (GCMC) simulations [48–52]. In a GCMC simulation, the volume, temperature, and chemical potential are fixed, and the number of molecules in the system fluctuates. For equilibrium between an adsorbed zeolite phase and a surrounding gas phase, the chemical potentials and temperatures must be equal in the two phases. By fixing the temperature and chemical potential in the simulation, the phase equilibrium condition is satisfied. For a single-component system, one obtains a point on the isotherm as follows. For a given gas-phase temperature and pressure, the chemical potential is calculated using the gas-phase equation of state. Then, the GCMC simulation is performed with the calculated chemical potential, the given temperature, and a zeolite volume of perhaps eight to 27 unit cells with periodic boundary conditions [40]. During the simulation, the average number of molecules in the system is calculated, yielding the amount adsorbed that corresponds to the gas-phase conditions (i.e. a point on the isotherm). Other properties such as the heat of sorption and preferred locations of sorbate molecules are also calculated. To sample from the grand canonical ensemble, three types of moves are performed. The first is a displacement of some set of molecular coordinates, such as a displacement of a molecule’s center of mass. The resulting change in the system energy is calculated and the move is accepted according to a Boltzmann weighting by the normal Metropolis method [40]. In the second type of move, a new molecule is inserted into the system at a randomly chosen position. The new configuration is accepted with a probability that depends on the new system energy and the chemical potential. In the third type of move, a molecule is randomly chosen to be removed, and again, the acceptance probability depends on the energy change and the chemical potential. Millions of these moves are performed to equilibrate the system and to obtain averages for properties of interest. Molecular simulations using atomistically-detailed models have been used to study adsorption in zeolites since the 1970’s [38,53–56], and recent years have produced rapid growth in the number of publications. An extensive tabular listing of molecular

L.J. Broadbelt, R.Q. Snurr / Applied Catalysis A: General 200 (2000) 23–46

simulations addressing adsorption and diffusion in zeolites through early 1995 has been compiled by Bell et al. [44]. Bates and van Santen have also reviewed this topic recently [57]. Early work focused on single-component systems and simple spherical adsorbates such as xenon and methane [48–51,58–61]. This was due to limitations in computational resources and the need to test the methods on simple systems before attempting more complicated problems. Recent work has also been devoted to multicomponent systems [52,62–75] due to widespread use of zeolites in separation processes. We discuss here only a few representative examples, with special emphasis on systems of interest in zeolite catalysis and on recent methodological developments that have allowed more complicated systems to be examined. One naturally thinks of chemisorption on an active site as being important for catalysis, but the importance of physisorption in the zeolite catalysis of hydrocarbons has also been pointed out by many authors [76–81]. For example, several groups [78–80] have recently performed studies in H-ZSM-5 zeolite showing that apparent increases in cracking rate with increasing n-alkane size are mostly due to increases in the adsorption constant of the hydrocarbons. That is, the intrinsic reaction rate of n-alkane cracking is about the same for all chain lengths, and the increased activity of longer molecules is due to the higher concentration of physisorbed molecules within the ZSM-5 pore space. Similarly, in comparing the cracking activity of a particular molecule in different zeolites, one needs to know the adsorption isotherms and heats of adsorption to allow for meaningful comparison [82]. One problem in using sorption data in interpreting catalytic activities is that calorimetry and isotherm measurements cannot be performed under reaction conditions and extrapolations can be questionable. In addition, multicomponent adsorption measurements are difficult to perform. A combination of sorption experiments and molecular simulations, however, can provide a solution to this problem. The strategy is first to do simulations at low temperatures and compare with experiment to validate the modeling. Then, simulations are done at the reaction conditions of interest where experiments are difficult or impossible. Based on the discussion so far, however, such simulations would be doomed to failure; random

31

insertion moves required in conventional GCMC simulations have a vanishingly low acceptance probability for longer alkane species due to overlaps of the molecules with atoms of the zeolite framework [83]. Configurational-bias Monte Carlo (CB MC) methods [84–86], however, have been developed and applied to zeolite adsorption to overcome the problems associated with random insertions of larger molecules. In CB MC, a molecule is ‘grown’ atom by atom in the zeolite avoiding unfavorable overlaps with the zeolite atoms. This introduces a bias into the algorithm that must be removed by changing the acceptance rules for the move. CB MC moves can be incorporated into various types of Monte Carlo calculations such as canonical ensemble (NVT) Monte Carlo, GCMC [73,74], or Gibbs ensemble MC. Bates et al. [87,88] used CB MC to look at the energetics, siting, and conformations of n-alkanes in a variety of all-silica zeolites. They found that in the cages of zeolite A the molecules adsorb in highly coiled conformations, while in the more open supercages of faujasite the conformations are similar to those for gas-phase alkanes. They report heats of adsorption for n-alkanes of varying chain length as a function of zeolite pore diameter and point out that it is necessary to obtain information on the location of the molecules within the zeolites to fully understand the trends obtained. For example, one may quote a pore diameter of 4 Å for zeolite A, corresponding to the window diameter, but since the molecules preferentially adsorb in the alpha cages, the more relevant pore size here is the diameter of the alpha cage, which is approximately 11.5 Å. Maginn et al. [89] investigated low-occupancy adsorption thermodynamics of n-alkanes as long as C25 in silicalite, using CB MC to evaluate the configurational integral of Eq. (1). They found that at high temperature (650 K) the molecules probe both the straight and zigzag channels of silicalite, favoring high-entropy conformations that ‘bend’ to go through both kinds of channels. At lower temperatures (300 K), short chains continue to populate all regions of the zeolite, but chains longer than n-octane tend to align in the straight channels. This localization produces an increase around C8 in the slope of the isosteric heat versus carbon number as shown in Fig. 3. This illustrates the possible danger in simply extrapolating low-temperature adsorption measurements up to catalytic conditions — a problem that arises because of

32

L.J. Broadbelt, R.Q. Snurr / Applied Catalysis A: General 200 (2000) 23–46

Fig. 3. Isosteric heats of adsorption for n-alkanes in silicalite at 300 K calculated by Maginn et al. [89] from Monte Carlo integration. Localization in the straight channels for molecules larger than octane causes a change in the slope of the predicted isosteric heat versus carbon number. The lines shown in the figure are empirical correlations derived from experimental data. The correlations have been extrapolated beyond the range over which they were fit. Reprinted with permission from J. Phys. Chem. 99 (1995) 2057–2079. Copyright 1995 American Chemical Society.

the sometimes subtle interplay between energetic and entropic effects in zeolite adsorption. In their simulated isotherms for n-alkanes in silicalite, Smit and Maesen [90] found that short chains (C1 –C5 ) and long chains (longer than C10 ) display simple isotherms, but hexane and heptane have kinked isotherms. This is also observed experimentally. The simulations suggest a new kind of phase transition, arising from the interplay of chain length and the length of the zigzag channels. When these two are comparable, the molecules can ‘freeze’ into these channels, leaving the channel intersections free so that molecules can now also access the straight channels. The straight channels fill after the mild step in the isotherm. Vlugt et al. [73] predicted adsorption isotherms for linear and branched alkanes and their mixtures in silicalite. They also observed interesting siting behavior and stepped isotherms for branched 2-methylalkanes for all carbon numbers studied (C4 –C9 ). Macedonia and Maginn [91] used CB GCMC simulations to predict single-component and binary sorption of light hydrocarbons in silicalite

and used their simulated data to test the ability of ideal adsorbed solution theory to predict binary sorption from the single-component isotherms. Other kinds of biased GCMC insertion moves were developed to allow simulation of benzene and pxylene in silicalite [92,93]. The simulations allowed an explanation of experimentally observed steps in the isotherms. The steps for the aromatics are thought to be related to a phase transition in the silicalite zeolite itself that has been found experimentally; this is in contrast to the explanation given above for steps in the alkane isotherms in silicalite. The aromatic simulations consumed a large amount of computer time, even with the special biasing of insertions that was developed to make the calculations possible. The simulations, as well as other simulation and experimental work, indicate that aromatic molecules in silicalite are adsorbed in clearly identifiable adsorption sites, suggesting the lattice representation in Fig. 4. In order to exploit the computational efficiency of a lattice model while preserving the predictive capability of atomistic simulations, a hierarchical approach was subsequently developed in which short atomistic simulations are used to obtain parameters that are then fed into the more efficient, but more coarse-grained, lattice model [94]. The parameters describe the free energy in the different adsorption sites due to molecule–zeolite interactions, as well as interactions between neighboring molecules. The hierarchical approach reproduces atomistic simulation results quite well, but it is over an order of magnitude more efficient computationally. The hierarchical approach was recently extended to adsorbed mixtures, and the effects of next-nearest-neighbor interactions were also studied [72]. 3.2. Diffusion From the same atomistic models used to predict adsorption isotherms, one can calculate transport properties. Equilibrium MD simulations provide shorttime dynamics information, and for many systems, MD simulations can access the time scales of diffusive motion. While a Monte Carlo simulation provides averages by stochastically generating new system configurations, an MD simulation samples configuration space deterministically by integrating the classical Newtonian equations of motion. The forces

L.J. Broadbelt, R.Q. Snurr / Applied Catalysis A: General 200 (2000) 23–46

33

Fig. 4. Schematic drawing of the lattice of adsorption sites for benzene in silicalite. Sites in the channel intersections, straight channels, and zigzag channels are denoted by I, S, and Z, respectively. This coarse-grained picture of the system was suggested by results of atomistic GCMC simulations. Using this lattice model, adsorption isotherms were obtained using a predictive hierarchical approach. Reprinted with permission from J. Phys. Chem. 98 (1994) 5111–5119. Copyright 1994 American Chemical Society.

that atoms feel are calculated from the potential energy model, and finite difference schemes are used to integrate the equations of motion to generate trajectories in time [40,41,43,95]. The self-diffusivity can then be calculated from the molecules’ mean-squared displacement according to the Einstein equation. A large number of MD studies of diffusion in zeolites have appeared in recent years, and Demontis and Suffritti have written a comprehensive review [95]. Only a few examples are quoted here, one describing work used to interpret catalytic results and several describing recent developments in using MD to study multicomponent diffusion in zeolites. Webb and Grest [96] recently used MD simulations to calculate the self-diffusion coefficients of normal decane and n-methylnonanes (n=2, 3, 4, and 5) at catalytically relevant temperatures in seven 10-member-ring zeolites. Two general types of behavior can be observed in Fig. 5a and b for the change in self-diffusivity as the branch position is moved toward the center of the alkane chain. For zeolites MFI, MEL, and MTT, 3 the self-diffusivity 3 The zeolites are identified by their three-letter structure type codes [35].

decreases monotonically as expected from the bulkiness of the different isomers, but for zeolites FER, AEL, EUO, and TON, other trends are observed. The results in Fig. 5b are explained by Webb and Grest, using an analysis of structural features of the zeolites that provide unique hindrance to molecular motion along the main diffusion channels. For example, EUO has a one-dimensional 10-ring channel system that is interrupted periodically by alternating side pockets. The entrances to the side pockets are formed by 12-member rings, so molecules easily access the pockets, but diffusion is only possible in the main channel direction. The ability of the side pockets to trap C10 alkanes depends on the position of the branch point. Webb and Grest use their results to help understand catalytic behavior in these systems. They suggest that the order-of-magnitude differences in diffusion coefficient shown in Fig. 5a and b may play a role in determining the modified constraint index (CI∗ ), which is used to characterize 10-member-ring zeolites. CI∗ is defined [97] as the ratio of yields of 2-methylnonane to 5-methylnonane at 5% conversion of n-decane. Generally, a higher CI∗ indicates greater spatial constraint or smaller pore dimensions. While transition-state shape selectivity

34

L.J. Broadbelt, R.Q. Snurr / Applied Catalysis A: General 200 (2000) 23–46

Fig. 5. Molecular dynamics predictions of Webb and Grest [96] for the self-diffusivity of n-methylnonanes in 10-member-ring zeolites at 600 K. For the three zeolites shown in (a), the diffusivity decreases monotonically with higher branch position. For the zeolites shown in (b), the diffusivity is not a monotonic function of branch position. (Reprinted with permission from Catal. Lett. 56 (1998) 95–104. Copyright 1998 Baltzer Science Publishers B.V.)

has been quoted as the primary factor in determining methylnonane yields upon isomerizing n-decane, Webb and Grest argue that diffusion may also play a role. They find that, for five of the seven zeolites studied, ordering them according to their CI∗ is the same as ordering them according to the ratio of the diffusivities for 2-methylnonane to 5-methylnonane. Their results suggest that product shape selectivity, based on differences in product diffusion coefficient, cannot be fully discounted. Despite the success of MD for predicting singlecomponent self-diffusivities in zeolites, MD simulations of binary diffusion in zeolites have appeared only recently. In 1997, Snurr and Kärger [98] published results from MD simulations and pulsed field gradient (PFG) NMR diffusion measurements for

methane and CF4 mixtures in silicalite. They found generally good agreement between the predicted and measured self-diffusivities. Jost et al. performed similar studies for methane and xenon mixtures in silicalite [99] and found that high xenon concentrations slowed the methane diffusivity in the mixture but that the diffusion of xenon was nearly unaffected by high methane concentrations. Sholl and Fichthorn [100] investigated a number of different mixtures in AlPO4 -5, using smart Monte Carlo, a method that combines short MD trajectories with a Monte Carlo sampling procedure. Gergidis and Theodorou [101] studied n-butane–methane mixtures in silicalite with MD simulations and used visualizations to reveal the jump-like character of intracrystalline motion for this system. Since catalytic applications of zeolites are not at equilibrium, it would be desirable to calculate the so-called ‘transport’ diffusivity as defined by Fick’s first law. The name transport diffusivity is often applied to the non-equilibrium Fick’s law diffusivity to distinguish it from the equilibrium self-diffusivity defined by the Einstein equation [42,43]. The problem of how to calculate a transport diffusivity for singlecomponent zeolite systems has been approached in several ways. Fritzsche et al. [102] designed non-equilibrium MD simulations that maintained a steady-state concentration gradient across the simulation cell to mimic a Wicke–Kallenbach permeation study. From the resulting flux, they obtained the transport diffusivity for methane in an A-type zeolite. They found that the concentration dependence was satisfactorily represented by the thermodynamic factor dln P/dln c as often assumed, where c is the concentration of adsorbate in equilibrium with pressure P. Maginn et al. [103] performed somewhat different simulations to mimic the transient conditions prevailing in an uptake experiment. In their technique, a concentration gradient was initially set up within the zeolite simulation cell and the microscopic equations of motion were integrated using a constant kinetic energy MD algorithm. By monitoring the relaxation of the concentration gradient as a function of time and fitting this relaxation to the continuum solution of Fick’s second law, the transport diffusivity for methane in silicalite was obtained. In a second method, Maginn et al. [103] applied a perturbing field to the system and observed the resulting steady-state

L.J. Broadbelt, R.Q. Snurr / Applied Catalysis A: General 200 (2000) 23–46

response of the system. This so-called ‘color field’ method was judged to be the better of the two approaches because of its computational efficiency and concerns that the gradient relaxation simulations were not in the linear response regime. The kinds of studies described so far (single-component and binary self-diffusivities and single-component transport diffusivities) represent an impressive demonstration of the ability of simulation to predict experimental diffusivities in zeolites. However, what the catalytic scientist really needs are multicomponent, counter-diffusivities at reaction temperatures [78]. Building on other work published to date, an approach has recently appeared for calculating binary transport diffusivities from MD simulations in zeolites [104]. Binary transport diffusivities were reported for 50:50 mixtures of methane and CF4 in a faujasite zeolite at several loadings. The cross-term diffusivities (methane-CF4 or CF4 -methane) were found to be an order of magnitude smaller than the main-term diffusivities (methane-methane or CF4 -CF4 ). Also reported were the phenomenological coefficients used in non-equilibrium thermodynamic treatments of transport [43]. The computation of the transport diffusivities requires considerably more computer time than simple MD simulations of equilibrium self-diffusivities, but with today’s computational power, such simulations will likely become common. The diffusion coefficients or phenomenological coefficients obtained from such simulations can be used directly in macroscopic engineering equations for transport under situations of co-diffusion, counter-diffusion, or other scenarios [105]. Then, one can evaluate the effect, if any, of the cross terms on real applications. The discussion of diffusion calculations so far has concentrated on MD simulations. There are many interesting systems, however, where the time scales of diffusive motion are beyond those accessible to current MD simulations. A rough guideline is that, with current computers, MD can simulate systems having diffusivities higher than about 10−11 m2 s−1 . To overcome the inherent limitations of MD, researchers have turned to hierarchical methods such as TST [106–111] and Brownian dynamics [112]. A discussion of such work is given in the review by Theodorou et al. [43]. TST is well suited to systems where diffusion occurs by infrequent hops between particular sites in

35

the zeolite. In typical TST calculations, an atomistic model like that used in MD is first used to identify the sorption sites and rate constants for the hops between them. Then, a more coarse-grained lattice model is used to study the long-time dynamics. Molecules hop between the lattice sites with probabilities governed by the TST rate constants. By incorporating reaction events into the lattice simulation, combined reaction and diffusion can be studied [113–118]. Such ‘kinetic Monte Carlo’ simulations can play a complementary role to microkinetic modeling, as described below.

4. Microkinetic modeling Microkinetic modeling is an ideal framework for assembling the microscopic information provided by atomistic simulations and electronic structure calculations to obtain macroscopic predictions of physical and chemical phenomena in systems involving chemical transformations [119]. The thrust of the approach is to write down the particular catalytic reaction mechanism of interest in terms of its most elementary steps. In contrast to more traditional approaches (e.g. power law kinetics, Langmuir– Hinshelwood–Hougen–Watson (LHHW) formulations), no rate-determining mechanistic step (RDS) is assumed. Microkinetic models are therefore much more widely applicable than traditional models which assume an RDS, since the RDS can change with reaction conditions. Because all postulated elementary steps are included explicitly, accurate rate parameters for all of the forward and reverse reactions are needed to solve the equations comprising the model. This requirement greatly increases the amount of information needed to create a microkinetic model, but this is also the power of the technique. Often, the information needed to create a microkinetic model cannot be determined from one set of experiments but requires the compilation of information from many different experiments and theoretical investigations at a variety of length scales. The interplay among the various levels in the modeling hierarchy of Fig. 1 is best illustrated through the details of the approach for developing a microkinetic model [119]. After the mechanistic steps of the catalytic reaction are proposed, balance equations for each of the species involved in the reaction, including

36

L.J. Broadbelt, R.Q. Snurr / Applied Catalysis A: General 200 (2000) 23–46

surface species and active sites, are assembled in the same manner routinely used to capture the chemistry of homogeneous reactions: Ωi =

n X νji rj

(2)

j =1

where Ω i is the turnover frequency (TOF) for species i in the units of molecules per site per second, n the total number of elementary reactions, ν j i the stoichiometric coefficient for species i in reaction j, and where rj is the rate of reaction j, in the same units as the TOF. The reaction rates require specification of rate constants. Although experimental values may be used when available and accurate, electronic structure calculations and TST can be used to afford rate constants for each step comprising the reaction mechanism. Furthermore, these computational techniques may be interfaced directly with the microkinetic model as subroutines to further integrate the levels of the modeling hierarchy. The Ω i values can then be coupled with the reactor design equation appropriate for the reactor configuration of interest. Example reactor design equations are provided in Eqs. (3) and (4) for a plug flow (PFR) and a transient continuously stirred tank (CSTR) reactor, respectively. PFR: Z 1 0 Ωi dSR (3) Fsi − Fsi = SR CSTR: Fsi dNs Fs Fsi dFs dFsi = (Fsi0 − Fsi + Ωi ) − + dt Ns Fs dt Ns dt (4) In the above equations, Fsi the molecular flow rate of species i per active site, SR the total number of active sites, Fs the total molecular flow rate per active site, Ns the total number of gas phase molecules per active site, and where superscript ‘0’ indicates an initial value. The system of equations formed by inserting the Ω i values into the reactor design equations can be solved by standard integration techniques [120], yielding both gas-phase concentrations and surface species’ coverages as a function of time (or distance down the reactor in the case of the PFR). By capturing the chemistry at this level of detail, valuable information can be obtained about the abundance of surface species,

the relative rates of reaction, and the extent to which each reversible reaction deviates from equilibrium. An early example of microkinetic analysis is the work of Bush and Dyer who recognized that a model that explicitly included surface species and did not assume an RDS provided the best opportunity for developing a model that could be used beyond the conditions for which it was parameterized [121]. In spite of this advantage, the widespread use of microkinetic models was retarded by the vast amount of information needed about interactions of chemical intermediates with complex, heterogeneous catalysts. With simultaneous advances in spectroscopic, isotopic tracing and other experimental methods, quantum chemical techniques and even solution methodologies, microkinetic models have become more prevalent. The pioneering work of Dumesic et al. and their publication of ‘The Microkinetics of Heterogeneous Catalysis’ [119] in 1993 greatly spurred the advance of microkinetic modeling. A recent volume of Catalysis Today [122] is testament to its position as an important component of the catalytic scientists’ toolbox. The microkinetic approach has been applied to numerous diverse chemistries including hydrocarbon cracking [123–133], isobutane dehydrogenation [134–136], ethane hydrogenolysis [137], ethylene hydrogenation [138–140], partial oxidation of methane [141], methane dimerization [142,143], methanation and Fischer–Tropsch synthesis [144–147], NO reduction [148], water gas shift chemistry [149–152], ammonia synthesis [153–167] and model systems [147,168,169]. To highlight the microkinetic approach and future directions for its development, we will discuss three specific chemistries: ammonia synthesis, hydrocarbon cracking and Fischer–Tropsch synthesis and related process chemistries. The chemistry that has been studied most widely using microkinetic modeling is ammonia synthesis, the production of ammonia from hydrogen and nitrogen, over both iron [153,154,161–167] and ruthenium [155–160] catalysts. The elementary steps of the models describing reaction over either metal are in general agreement: reaction proceeds through adsorption and dissociation of hydrogen and nitrogen and sequential hydrogenation of nitrogen to form ammonia. In some cases, the molecularly adsorbed precursor state of hydrogen [166,167] or nitrogen [158,160] is not included. Dumesic et al. [119] critically examined the

L.J. Broadbelt, R.Q. Snurr / Applied Catalysis A: General 200 (2000) 23–46

microkinetic models of Stoltze and Nørskov [166,167] and Bowker et al. [164,165] for ammonia synthesis over iron and were able to show that, with the correct parameter set, a model could extrapolate 11 orders of magnitude in pressure from adsorption and desorption studies under UHV conditions to reproduce the performance of an industrial iron catalyst. Sehested et al. [163] re-examined these models and that of Fastrup [161] to determine their ability to capture experimental data over a multiply-promoted iron catalyst. A binding energy derived for potassium-promoted Fe(1 1 1) at low coverage and variation of three other model parameters were both required to obtain an improved description of their data. One of the main differences among the models was the enthalpy of nitrogen chemisorption, which ranged from 44.4 to 106.9 kJ mol−1 . The values of 84.8 and 77.9 kJ mol−1 used by Stoltze and Nørskov [166,167] and Sehested et al. [163], respectively, are in better agreement with surface science studies over single crystals. However, adsorbate–adsorbate interactions are known to exist for dinitrogen adsorption [153]. Recent work using DFT by Mortensen et al. [25] may help to quantify the role of these interactions under working conditions of the catalyst. The recent interest in commercial ruthenium-based ammonia synthesis catalysts has spawned development of a microkinetic model by Hinrichsen and coworkers in parallel with experiments [156–160]. The model parameters were built up sequentially from temperature-programmed adsorption and desorption experiments, isotopic exchange reactions and temperature-programmed surface reaction experiments. The results of their microkinetic analysis suggested that the maximum rate of N2 dissociation would be several orders of magnitude lower than the observed rate of NH3 formation over the same catalyst if the low sticking coefficient obtained from Ru(0 0 0 1) single-crystal studies were used. Thus, the model was instrumental in pointing out the importance of interactions of ruthenium metal with the MgO support [156]. Similar to the success obtained for iron catalysts, the microkinetic model was able to bridge the ‘pressure gap’ effectively. Further, the work of Hinrichsen et al. provides an example of incorporating non-uniform behavior in the form of a coverage-dependent rate constant to improve the predictive capability of the model.

37

A second chemistry to which microkinetic modeling has been broadly applied is catalytic cracking. The importance of fluidized catalytic cracking to the refining industry has led to a range of models being developed over the years, and most recently, microkinetic models have been advanced to capture the changes in product selectivities and conversion with different catalyst formulations and process conditions for various model feeds [123–133]. The first effort to put the kinetics of catalytic cracking on firm mechanistic ground using microkinetic modeling was by Dumesic and coworkers [123–128]. In early work, they developed a microkinetic model for isobutane cracking over calcined and steamed Y-zeolite catalysts based on carbocation surface chemistry. Using the experimental literature as a guide, a serviceable reaction mechanism was formulated that consisted of 21 reaction steps, including initiation, oligomerization, ␤-scission, olefin desorption, isomerization and hydride transfer. Kinetic parameters were estimated using a combination of TST, the Evans–Polanyi correlation [170], and gas-phase thermodynamic data. To account for the catalyst surface, gas-phase heats of formation were corrected by a catalyst-specific parameter, 1H+ , which was the heat of stabilization of a proton in the zeolite framework. Alternatively, heats of formation of surface species could be specified, using heats of adsorption obtained experimentally or from atomistic simulations as described above. The microkinetic model was able to capture the experimental data very well, particularly at low conversions. Drawing upon the advantages of the microkinetic approach, they quantified the coverages of various carbenium ion intermediates as a function of reactant conversion. A different combination of surface coverages and fitted parameters could give results in equal agreement with the experimental data, and thus, the surface coverages reported are not absolute. However, the model results did allow comparison among different catalysts. More recently, they were able to quantify the links between steaming, catalyst cracking activity and acid site strength [128]. A series of papers by Watson and coworkers [129–133] extended the microkinetic modeling approach of Dumesic and coworkers [123–128] to other model feeds, including longer chain alkanes, alkylbenzenes and their mixtures, with equal success. The most notable achievement in the work of Watson et al. was the ability

38

L.J. Broadbelt, R.Q. Snurr / Applied Catalysis A: General 200 (2000) 23–46

to predict binary mixture data, using model parameters previously determined from single-component data without adjustment. This highlights the power of the microkinetic modeling approach. Since the chemistry is described at a fundamental level, interactions between components can be described without redevelopment or reparameterization of the model. Another complex chemistry that involves transformations of many different surface species to which microkinetic modeling has been applied is the hydrogenation of carbon monoxide. Methanation over nickel-based catalysts where molecular weight growth is limited and high yields of C1 species result has been most extensively examined [144,145,147]. Alstrup [144] developed a microkinetic model for methanation based on CO dissociation and stepwise hydrogenation of surface carbon. Excellent agreement between the model results and the experimental data of Goodman et al. [171] over a Ni(1 0 0) surface was obtained. However, in order to obtain this agreement, the coverage of reactive carbon was treated as a parameter and a rate-controlling step was assumed; both assumptions are counter to the true microkinetic modeling approach. Alstrup did demonstrate that traditional models previously published did not agree well either with the single-crystal data or with methanation data over a nickel foil, a second data set that he was able to capture with his model. A microkinetic model developed by Aparicio [145] was applied to both methanation and steam reforming, chemistries which have common elementary steps. His model was able to capture the steam reforming data of Xu and Froment [172] through adjustment of eight parameters. As a test of the model’s predictive capability, the model was used without adjustment of any of the parameters to describe published data for the kinetics of CO methanation over several Ni/Al2 O3 catalysts (see Table 2 in [145]). Although the absolute prediction of the reaction rate was an order of magnitude too low, all of the major trends were correctly predicted. However, when the model was used to model the single-crystal data of Goodman et al., the carbon coverage had to be restricted to a maximum of one-half of a monolayer. To improve upon this model, Dooling et al. [147] incorporated elements of non-uniformity into a microkinetic model of methanation. Although surface non-uniformity (e.g., structure sensitivity, adsorbate-induced restructuring, adsorbate–adsorbate

interactions) has long been known to affect kinetics on heterogeneous catalytic surfaces [173–175], its incorporation into all types of kinetic models is difficult. To introduce known elements of surface non-uniformity, Dooling et al. included the change in heat of adsorption of carbon monoxide with coverage, the dependence of binding energy of carbon [30] and CH [31] with coverage and the mutual enhancement of coadsorbed carbon monoxide and hydrogen adsorption. The non-uniform microkinetic model was able to capture the transient data of Aparicio more fully than previous microkinetic models through the adjustment of only three parameters. In addition, the model was better able to capture the single-crystal data of Goodman et al. with no adjustment of the parameters or artificial capping of the maximum carbon coverage as shown in Fig. 6. The collapse of the variable pressure data to a single line at low temperature was reproduced, and the point at which the three pressure curves begin to deviate as a function of temperature was also captured well. The microkinetic models of methanation, a subset of Fischer–Tropsch (FT) synthesis, are relatively compact, consisting of nine to 15 reversible elementary steps. In order to describe molecular weight growth chemistry involved in higher hydrocarbon synthesis, carbon–carbon bond formation reactions must be included, and species with multiple carbon atoms must be tracked. Writing down all of the species leading to compounds with 15 or more carbon atoms would be tedious, time-consuming and prone to error. To facilitate construction of microkinetic models of large systems like FT chemistry, Broadbelt and coworkers have developed algorithms for computer-aided construction of reaction mechanisms [146,176–178]. Although the numbers of reactions and species are large, the number of reaction types involved in FT synthesis is small, as suggested by mechanistic postulates for initiation, propagation and termination put forth in the literature (e.g. [179,180]). With computer generation of reaction mechanisms, the modeler can focus on the information the model provides, rather than the tedious manual construction of the model. Although the algorithms and their implementation vary, many common features are found. The user provides a representation of reactants’ structure and information about the chemical reactions or products of interest, i.e. the mechanistic postulates based on experimental evidence. The structure of the reactants is transformed according to

L.J. Broadbelt, R.Q. Snurr / Applied Catalysis A: General 200 (2000) 23–46

39

Fig. 6. Prediction of the data of Goodman et al. [171] for methanation over a Ni(1 0 0) surface using a non-uniform microkinetic model with deactivation [147]. Reproduced with permission from Langmuir 15(18) (1999) 5846–5856. Copyright 1999 American Chemical Society.

the chemical reactions to form reaction intermediates and products, and the connectivity and uniqueness of the products are determined. Finally, the relationships among the species are compiled into a reaction mechanism and the controlling rate equations that describe it. In applying this approach to FT synthesis [146], it was necessary to quantify adsorption energies of each surface species that was generated. Because the majority of these values are not available, quantum chemistry was used to quantify interactions of FT-derived species on nickel and cobalt surfaces [30–32]. Because it was prohibitive to obtain adsorption energies of large intermediates from quantum chemical calculations, a phenomenological method [181–185] was used to link adsorption energies of high molecular weight fragments to atomic adsorption energies. A model with over 50 species and 200 reactions was able to model single-crystal data over a Co(0 0 0 1) surface [186] and capture the vast difference in molecular weight growth potential of nickel and cobalt surfaces. As alluded to above, an important factor in guiding the kinetics over heterogeneous catalysts is non-uniform behavior, i.e. crystal face selectivity, adsorbate-induced changes in the catalyst or adsorbate–adsorbate interactions. It is possible to use several different approaches to incorporate non-uniformity into microkinetic models. The most straightforward is to introduce coverage-dependent

activation energies for desorption. Using information from temperature-programmed desorption studies, microcalorimetric data or even quantum chemical calculations, much can be learned and quantified about adsorbate–adsorbate interactions. Examples in which this approach was adopted include ammonia synthesis over ruthenium [160] and hydrogen oxidation [187]. To include the effects of adsorbate–adsorbate interactions on the energetics of surface reactions, the concept of linear free energy relationships, which relate activation energies to adsorption energies, can be used successfully to incorporate non-uniformity into microkinetic models [147]. To take more complex interactions between multiple components into account, the functional dependence of multicomponent adsorption energies on coverage must be determined. However, ad hoc approaches may also be adopted. For example, extra sites were included to increase the capacity for hydrogen adsorption in the presence of carbon monoxide in the case of methanation chemistry [147]. It is also possible to extend microkinetic models further by including different types of sites with distinct adsorption energies that communicate through diffusion with no specific spatial relationship or explicitly including a spatial variable into the formulation of the model equations. However, developing the equations to include all of these complex interactions is more difficult and significantly

40

L.J. Broadbelt, R.Q. Snurr / Applied Catalysis A: General 200 (2000) 23–46

more information is required. Thus, the majority of the microkinetic models to date in the literature have used a mean field approach even when known non-uniformity is present. As an alternative to microkinetic models, kinetic Monte Carlo simulations maintain information about explicit spatial relationships among species on the catalyst surface [188,189] but are currently less widely adopted in practice. Such schemes provide a natural way to link atomistic modeling and reaction engineering. It is interesting to note that microkinetic modeling and atomistic simulations have traditionally been practiced by different communities of researchers, but lattice simulations have emerged as ‘common ground’ that facilitate description of reaction and diffusion in heterogeneous catalytic systems. Although microkinetic modeling has been shown to be a powerful tool to interpret phenomena in catalytic reaction systems, examples using microkinetic modeling to date have focused on incorporating chemical detail alone in systems where mass transfer effects are not important. Many industrial systems that may benefit from detailed microkinetic modeling are dominated by mass transfer limitations, and therefore, the microscopic description of the kinetics at the catalyst surface must be linked correctly with descriptions of external and internal mass transfer. Specifically, the concentrations used to calculate the TOFs on the catalyst surface must be related to characteristics of the bulk fluid as illustrated in Fig. 7. An important parameter characterizing internal mass transfer resistance

Fig. 7. Transport limitations both external to the catalyst surface and within the catalyst pores can affect the concentration at the active site (adapted from [199]).

is a multicomponent diffusion coefficient, which for microporous materials such as zeolites is generally not available. However, molecular simulations offer a method to obtain this information as described above.

5. Concluding thoughts 5.1. State-of-the-art The wide range of time and length scales at play in heterogeneous catalysis have led to the development and application of various modeling techniques suited to answering different questions at different scales. The examples described in this review give some impression of the current state of molecular modeling for catalytic systems. Quantum chemical methods are used routinely now to calculate the energies, electronic structures, and spectroscopic properties of small arrangements of atoms and even periodic structures. One can search for minimum-energy configurations of the atoms, and with somewhat more difficulty, saddle points can be found, corresponding to reaction transition states. Infrared spectra can be calculated with standard software. Energies of adsorption can be computed. Atomistic modeling, using empiricallyderived potential interactions, can provide adsorption isotherms, heats of adsorption, and diffusion coefficients from Monte Carlo and MD simulations. For systems where the diffusion times are beyond those accessible to MD, TST can be used to calculate rates of hopping between preferred sites, and these rates can be fed into more coarse-grained lattice models to predict the diffusion coefficient. Such lattice models can also be used to model combined reaction and diffusion in porous solids or on surfaces. Microkinetic modeling is used to link molecular-level information about reactants, products and reactive intermediates on heterogeneous surfaces to macroscopic kinetic observations. A mechanism is postulated based on experimental and theoretical evidence, rate constants for each elementary step are specified from experiment, correlations or theory, and the mechanism is combined with the appropriate reactor design equations. The full set of equations is solved, i.e. no rate-determining step (RDS) is assumed, and relative reaction rates, coverages of surface intermediates, reactant conversion, and product yields and selectivities are obtained.

L.J. Broadbelt, R.Q. Snurr / Applied Catalysis A: General 200 (2000) 23–46

5.2. Challenges Much exciting work is being done in molecular-level modeling of heterogeneous catalysis, but many challenges remain if these methods are to live up to their full promise. Problems remain at each of the length scales shown in Fig. 1, and the even larger challenge is to better link these qualitatively different types of calculations into a comprehensive approach. High-level quantum mechanical calculations remain prohibitively expensive for large systems, and one is always faced with trade-offs in terms of system size, quality of basis sets, and level of theory. The question of accuracy becomes especially important if one considers that an error of 8 kJ mol−1 in an activation energy leads to an error of approximately a factor of 10 (at 420 K) in predicted rate because of the Arrhenius dependence. Predictions with even higher accuracy remain a long-term goal. Another goal is to model ever-bigger systems with good accuracy. Slab calculations may be used for catalysts with regular repeat structures but are not yet well suited for systems with large unit cells or irregular defects. Periodic calculations do permit surface coverage effects to be explored, and thus, information about adsorbate–adsorbate interactions may be fed to other levels of the hierarchy. Embedding schemes also allow larger systems to be studied with true-order N performance. For example, Ellis and coworkers carry out quantitative first principles calculations on fragments of an extended system and couple them together through mutual boundary conditions and constraints on the electron density [190]. Modeling bigger systems will be essential for addressing many questions, e.g. the role of metal–support interactions in supported catalysts. Improvements in computing speed will continue to increase the amount of detail that can be included in atomistic modeling. For simulations of adsorption and diffusion relevant to zeolite catalysis, researchers can now include framework flexibility, polarization energies, and zeolitic cations. These features bring a new level of realism that was neglected in earlier work, but they make the simulations quite time-consuming. The force fields remain the Achilles’ heel of these methods. Two avenues for force field improvement should be simultaneously pursued: development of better empirical force fields and development of force fields based on quantum chemical calculations. Force fields for

41

polar molecules interacting with cations or catalytic sites are especially problematic due to the difficulties in dealing with electron transfer and polarization. Simulation of longer time scales remains a problem. MD simulations have inherent time-scale limitations that cannot be solved simply by faster computers. Multiple time step methods [191] are an improvement, but some problems will always have time scales orders of magnitude beyond what can be accessed with MD. For systems where molecules diffuse by hopping between preferred sites, TST can be used as described above. Such methods, however, do not yet appear in commercial molecular modeling software. It is not a trivial task to ‘generalize’ TST algorithms such that they can be written into software that could be readily applied to a variety of problems. Applications of TST to date have tended to require specially written software for particular purposes. Another method for simulating longer time scales was suggested recently by Voter [192,193]. This ‘hyper-MD’ method is fundamentally based on TST but can handle correlated events, and its further development appears promising. Microkinetic modeling assembles molecular-level information obtained from quantum chemical calculations, atomistic simulations and experiment to quantify the kinetic behavior at given reaction conditions on a particular catalyst surface. To formulate a microkinetic model, a reaction mechanism must be postulated. Although experimental evidence, chemical intuition and theory are all used together to propose a likely mechanism, full reaction mechanisms on heterogeneous surfaces are rarely known, particularly for complex systems in which many reaction intermediates may be involved. However, even in the absence of knowledge of the ‘true’ mechanism, microkinetic modeling provides a framework for testing various mechanistic postulates and may offer suggestions about the likelihood of experimentally-detected surface intermediates playing a key role in the observed kinetics. Specifying the rate parameters for each reaction in the mechanism presents particular challenges to developing microkinetic models. Much success has already been achieved by Dumesic and coworkers and other groups of researchers by using a combination of experimental values, kinetics correlations such as the Evans–Polanyi relationship, and TST to estimate activation energies and frequency factors. Often parameters must still be fitted to experimental data, and this

42

L.J. Broadbelt, R.Q. Snurr / Applied Catalysis A: General 200 (2000) 23–46

limits the predictive capability that microkinetic modeling inherently offers. Continued improvement in the accuracy and speed of calculating rate constants using ab initio quantum chemistry is one promising avenue for improving the predictions of microkinetic models. The majority of microkinetic models to date have adopted a mean field approach for describing the interactions of reactive intermediates with the catalyst surface. However, catalysts are known to be non-uniform and have also been shown to be dynamic under reaction conditions. Incorporating more realistic representations of catalyst surfaces into microkinetic models requires detailed catalyst characterization and implies that a larger set of parameters must be known. Furthermore, a spatial variable would need to be added to traditional microkinetic models to take into account local differences in reactivity. Explicitly modeling deactivation phenomena presents additional challenges for microkinetic modeling. For example, the chemical steps leading to coke formation on acid catalysts are not well resolved, and therefore, it is difficult to postulate a serviceable reaction mechanism for deactivation. Finally, examples using microkinetic modeling to date have focused primarily on incorporating only chemical detail in systems where mass transfer effects are not important. To include mass transfer effects, multicomponent diffusion coefficients must be specified, and a more explicit description of the catalyst (e.g. pore architecture, location of active sites etc.) may be required. Finally, there is the ‘grand’ challenge of linking these methods together. Reliable, easy-to-use commercial software is available for quantum chemistry. In some sense, this code is completely general and independent of the physical system of interest — the job is always to solve the Schrödinger equation. MD codes are similarly general — they solve Newton’s equations, given the potential interactions among atoms. Hierarchical modeling, on the other hand, is often system-specific, and generalized software has yet to appear. Some methods that may be helpful in creating better links between existing methods include the following. The electronegativity equalization method (EEM) may be useful in incorporating electronic structure information into atomistic simulations. Using a DFT-based formulation, EEM treats individual atoms as centers having an equilibrium electronegativity and hardness [194,195]. During a simulation, these constants and the external point charges on other atoms

determine the charge on individual atoms, thus incorporating charge fluctuations into the simulation. It may also be possible to use the Fukui function [196] and other reactivity indexes to link quantum chemical methods with atomistic simulations, in a scheme that uses a classical simulation to search for configurations that are most likely to react. These could then be further analyzed by detailed DFT calculations. Better linking of atomistic simulations with reaction engineering models may come from developments in TST, kinetic Monte Carlo, and field theories [197]. 5.3. The future Progress continues to be made in overcoming the challenges outlined above, and molecular modeling is rapidly becoming an essential tool for the design of catalysts and other materials. Molecular modeling can provide trends and new insights that may be difficult to determine experimentally, and it can be used to screen potential catalysts for particular applications — complementing and, in some cases, even replacing expensive experiments. Instead of synthesizing and testing a large number of potential catalysts, researchers in the future may use molecular modeling to ‘test’ the original candidates and subsequently perform experimental tests only on a small number of the most promising candidates. Modeling may also suggest radically different candidates that were not previously considered or help in optimization. It seems clear that molecular modeling is poised to play a key role in the development of future catalytic systems, especially in conjunction with new synthesis and characterization methods that allow ever-increasing control over the nature of the catalyst.

Acknowledgements We thank Profs. Alex Bell, Donald Ellis, Mike Klein, Harold Kung, and Doros Theodorou and Dr. Jim Rekoske for discussions that have helped form our thoughts on the topics of this article. Louis Clark is thanked for assistance with Fig. 1. We also thank the EMSI program of the National Science Foundation and the US Department of Energy Office of Science (CHE-9810378) at the Northwestern University Institute for Environmental Catalysis.

L.J. Broadbelt, R.Q. Snurr / Applied Catalysis A: General 200 (2000) 23–46

References [1] R.A. van Santen, Theoretical Heterogeneous Catalysis, World Scientific, Singapore, 1991. [2] R.A. van Santen, Chem. Eng. Sci. 45 (1990) 2001. [3] M. Witko, J. Mol. Catal. 70 (1991) 277. [4] F. Ruette (Ed.), Quantum Chemistry Approaches to Chemisorption and Heterogeneous Catalysis, Kluwer Academic Publishers, Boston, 1992. [5] R.A. van Santen, M. Neurock, Catal. Rev.-Sci. Eng. 37 (1995) 557. [6] J.W. Andzelm, A.E. Alvarado-Swaisgood, F.U. Axe, M.W. Doyle, G. Fitzgerald, C.M. Freeman, A.M. Gorman, J.R. Hill, C.M. Kolmel, S.M. Levine, P.W. Saxe, K. Stark, L. Subramanian, M.W. van Daelen, E. Wimmer, J.M. Newsam, Catal. Today 50 (1999) 451. [7] C. Pisani, J. Mol. Struct. (Theochem) 463 (1999) 125. [8] P. Blaha, K. Schwarz, P. Dufek, R. Augustyn, WIEN95: A Full Potential Linearized Augmented Plane Wave Package for Calculating Crystal Properties, Technical University Vienna, 1995. [9] M.J. Frisch, G.W. Trucks, H.B. Schlegel, G.E. Scuseria, M.A. Robb, J.R. Cheeseman, V.G. Zakrzewski, J.A. Montgomery Jr., R.E. Stratmann, J.C. Burant, S. Dapprich, J.M. Millam, A.D. Daniels, K.N. Kudin, M.C. Strain, O. Farkas, J. Tomasi, V. Barone, M. Cossi, R. Cammi, B. Mennucci, C. Pomelli, C. Adamo, S. Clifford, J. Ochterski, G.A. Petersson, P.Y. Ayala, Q. Cui, K. Morokuma, D.K. Malick, A.D. Rabuck, K. Raghavachari, J.B. Foresman, J. Cioslowski, J.V. Ortiz, A.G. Baboul, B.B. Stefanov, G. Liu, A. Liashenko, P. Piskorz, I. Komaromi, R. Gomperts, R.L. Martin, D.J. Fox, T. Keith, M.A. Al-Laham, C.Y. Peng, A. Nanayakkara, C. Gonzalez, M. Challacombe, P.M.W. Gill, B. Johnson, W. Chen, M.W. Wong, J.L. Andres, C. Gonzalez, M. Head-Gordon, E.S. Replogle, J.A. Pople, GAUSSIAN 98, Gaussian, Pittsburgh, PA, 1998. [10] A.R. Leach, Molecular Modelling: Principles and Applications, Addison-Wesley Longman, Harlow, England, 1996. [11] F. Jensen, Introduction to Computational Chemistry, Wiley, New York, 1999. [12] M. Haruta, Catal. Today 36 (1997) 153. [13] U. Heiz, A. Sanchez, S. Abbet, W.-D. Schneider, J. Am. Chem. Soc. 121 (1999) 3214. [14] D.J. Dooling, L.J. Broadbelt, in: G.F. Froment, K.C. Waugh (Eds.), Dynamics of Surfaces and Reaction Kinetics in Heterogeneous Catalysis, Vol. 109, Elsevier, Amsterdam, 1997, p. 251. [15] J. Sauer, U. Eichler, U. Meier, A. Schafer, M. von Arnim, R. Ahlrichs, Chem. Phys. Lett. 308 (1999) 147. [16] J. Gao, in: K.B. Lipkowitz, D.B. Boyd (Eds.), Reviews in Computational Chemistry, Vol. 7, VCH, New York, 1996, p. 119. [17] P. Cremaschi, J.L. Whitten, Theor. Chem. Acta 72 (1987) 485. [18] J.L. Whitten, H. Yang, Int. J. Quantum Chem. 29 (1995) 41. [19] J.L. Whitten, Chem. Phys. 177 (1993) 387.

43

[20] H. Yang, J.L. Whitten, Surf. Sci. 401 (1998) 312. [21] M. Brändle, J. Sauer, R. Dovesi, N.M. Harrison, J. Chem. Phys. 109 (1998) 10379. [22] B. Hammer, J.K. Nørskov, Phys. Rev. Lett. 79 (1997) 4441. [23] J.J. Mortensen, Y. Morikawa, B. Hammer, J. Nørskov, J. Catal. 169 (1997) 85. [24] J.J. Mortensen, M.V. Ganduglia-Pirovano, L.B. Hansen, B. Hammer, P. Stoltze, J.K. Nørskov, Surf. Sci. 422 (1999) 8. [25] J.J. Mortensen, L.B. Hansen, B. Hammer, J.K. Nørskov, J. Catal. 182 (1999) 479. [26] G. Ertl, M. Grunze, M. Weiss, J. Vac. Sci. Technol. 13 (1976) 314. [27] J.H. Sinfelt, Bimetallic Catalysts: Discoveries, Concepts and Applications, Wiley, New York, 1983. [28] V. Pallassana, M. Neurock, L.B. Hansen, B. Hammer, J.K. Nørskov, Phys. Rev. B 60 (1999) 6146. [29] V. Pallassana, M. Neurock, G.W. Coulston, J. Phys. Chem. B 103 (1999) 8973. [30] D.J. Klinke, S. Wilke, L.J. Broadbelt, J. Catal. 178 (1998) 540. [31] D.J. Klinke, D.J. Dooling, L.J. Broadbelt, Surf. Sci. 425 (1999) 334. [32] D.J. Klinke, L.J. Broadbelt, Surf. Sci. 429 (1999) 169. [33] R. Car, M. Parrinello, Phys. Rev. Lett. 55 (1985) 2471. [34] M. Boero, M. Parrinello, K. Terakura, Surf. Sci. 438 (1999) 1. [35] W.M. Meier, D.H. Olson, Atlas of Zeolite Structure Types, Butterworth-Heinemann, London, 1992. [36] M. Falcioni, M.W. Deem, J. Chem. Phys. 110 (1999) 1754. [37] G.C. Maitland, M. Rigby, E.B. Smith, W.A. Wakeham, Intemolecular Forces: Their Origin and Determination, Clarendon Press, Oxford, 1981. [38] A.G. Bezus, A.V. Kiselev, A.A. Lopatkin, Q.D. Pham, J. Chem. Soc., Faraday Trans. II 74 (1978) 367. [39] A. Bondi, Physical Properties of Molecular Crystals, Liquids, and Glasses, Wiley, New York, 1968. [40] M.P. Allen, D.J. Tildesley, Computer Simulation of Liquids, Clarendon Press, Oxford, 1987. [41] D. Frenkel, B. Smit, Understanding Molecular Simulation: From Algorithms to Applications, Academic Press, San Diego, 1996. [42] J. Kärger, D.M. Ruthven, Diffusion in Zeolites and Other Microporous Solids, Wiley, New York, 1992. [43] D.N. Theodorou, R.Q. Snurr, A.T. Bell, in: G. Alberti, T. Bein (Eds.), Comprehensive Supramolecular Chemistry: Solid-state Supramolecular Chemistry: Twoand Three-Dimensional Networks, Vol. 7, Pergamon Press, Oxford, 1996, p. 507. [44] A.T. Bell, E.J. Maginn, D.N. Theodorou, in: G. Ertl, H. Knözinger, J. Weitkamp (Eds.), Handbook of Heterogeneous Catalysis, VCH, Weinheim, 1997. [45] A.V. Kiselev, A.A. Lopatkin, A.A. Shulga, Zeolites 5 (1985) 261. [46] R.L. June, A.T. Bell, D.N. Theodorou, J. Phys. Chem. 94 (1990) 1508. [47] W.H. Press, B.P. Flannery, S.A. Teukolsky, W.T. Vetterling, Numerical Recipes: The Art of Scientific Computing, Cambridge University Press, Cambridge, 1986.

44

L.J. Broadbelt, R.Q. Snurr / Applied Catalysis A: General 200 (2000) 23–46

[48] J.L. Soto, A.L. Myers, Mol. Phys. 42 (1981) 971. [49] G.B. Woods, J.S. Rowlinson, J. Chem. Soc., Faraday Trans. 85 (1989) 433. [50] F. Karavias, A.L. Myers, Mol. Simulation 8 (1991) 23. [51] R.Q. Snurr, R.L. June, A.T. Bell, D.N. Theodorou, Mol. Simulation 8 (1991) 73. [52] D.M. Razmus, C.K. Hall, AIChE J. 37 (1991) 769. [53] R.I. Derrah, D.M. Ruthven, Can. J. Chem. 53 (1975) 996. [54] A.V. Kiselev, Adv. Chem. 102 (1971) 37. [55] R.-G. Kretschmer, K. Fiedler, Z. Phys. Chemie (Leipzig) 258 (1977) 1045. [56] H.J.F. Stroud, E. Richards, P. Limcharoen, N.G. Parsonage, J. Chem. Soc., Faraday Trans. I 72 (1976) 942. [57] S.P. Bates, R.A. van Santen, Adv. Catal. 42 (1998) 1. [58] R.F. Cracknell, K.E. Gubbins, Langmuir 9 (1993) 824. [59] V. Lachet, A. Boutin, R.J.-M. Pellenq, D. Nicholson, A.H. Fuchs, J. Phys. Chem. 100 (1996) 9006. [60] C.J. Jameson, A.K. Jameson, B.I. Baello, H.-M. Lim, J. Chem. Phys. 100 (1994) 5965. [61] P.R. Van Tassel, H.T. Davis, A.V. McCormick, J. Chem. Phys. 98 (1993) 8919. [62] F. Karavias, A.L. Myers, Mol. Simulation 8 (1991) 51. [63] M.W. Maddox, J.S. Rowlinson, J. Chem. Soc., Faraday Trans. 89 (1993) 3619. [64] P.R. Van Tassel, H.T. Davis, A.V. McCormick, Langmuir 10 (1994) 1257. [65] P.R. Van Tassel, H.T. Davis, A.V. McCormick, Mol. Simulation 17 (1996) 239. [66] J.A. Dunne, A.L. Myers, D.A. Kofke, Adsorption 2 (1996) 41. [67] C.J. Jameson, A.K. Jameson, H.M. Lim, J. Chem. Phys. 104 (1996) 1709. [68] M. Heuchel, R.Q. Snurr, E. Buss, Langmuir 13 (1997) 6795. [69] V. Lachet, A. Boutin, B. Tavitian, A.H. Fuchs, Faraday Discuss. 106 (1997) 307. [70] V. Lachet, A. Boutin, B. Tavitian, A.H. Fuchs, J. Phys. Chem. B 102 (1998) 9224. [71] L.A. Clark, A. Gupta, R.Q. Snurr, J. Phys. Chem. B 102 (1998) 6720. [72] K.F. Czaplewski, R.Q. Snurr, AIChE J. 45 (1999) 2223. [73] T.J.H. Vlugt, R. Krishna, B. Smit, J. Phys. Chem. B 103 (1999) 1102. [74] M.D. Macedonia, E.J. Maginn, Mol. Phys. 96 (1999) 1375. [75] A. Gupta, L.A. Clark, R.Q. Snurr, Langmuir, 16 (2000) 3910. [76] J.F. Denayer, G.V. Baron, W. Souverijns, J.A. Martens, P.A. Jacobs, Ind. Eng. Chem. Res. 36 (1997) 3242. [77] E.G. Derouane, J. Mol. Catal. A 134 (1998) 29. [78] W.O. Haag, Stud. Surf. Sci. Catal. 84 (1994) 1375. [79] T.F. Narbeshuber, H. Vinek, J.A. Lercher, J. Catal. 157 (1995) 388. [80] J. Wei, Chem. Eng. Sci. 51 (1996) 2995. [81] B.A. Williams, S.M. Babitz, J.T. Miller, R.Q. Snurr, H.H. Kung, Appl. Catal. A 177 (1999) 161. [82] S.M. Babitz, B.A. Williams, J.T. Miller, R.Q. Snurr, W.O. Haag, H.H. Kung, Appl. Catal. A 179 (1999) 71. [83] B. Smit, J.I. Siepmann, Science 264 (1994) 1118. [84] J.I. Siepmann, D. Frenkel, Mol. Phys. 75 (1992) 59.

[85] D. Frenkel, G.C.A.M. Mooij, B. Smit, J. Phys. Condens. Matter 4 (1992) 3053. [86] J.J. dePablo, M. Laso, U.W. Suter, J. Chem. Phys. 96 (1992) 3295. [87] S.P. Bates, W.J.M. van Well, R.A. van Santen, B. Smit, J. Phys. Chem. 100 (1996) 17573. [88] S.P. Bates, W.J.M. van Well, R.A. van Santen, B. Smit, J. Am. Chem. Soc. 118 (1996) 6753. [89] E.J. Maginn, A.T. Bell, D.N. Theodorou, J. Phys. Chem. 99 (1995) 2057. [90] B. Smit, T.L.M. Maesen, Nature 374 (1995) 42. [91] M.D. Macedonia, E.J. Maginn, Fluid Phase Equilibria 158–160 (1999) 19. [92] R.Q. Snurr, A.T. Bell, D.N. Theodorou, J. Phys. Chem. 97 (1993) 13742. [93] L.A. Clark, R.Q. Snurr, Chem. Phys. Lett. 308 (1999) 155. [94] R.Q. Snurr, A.T. Bell, D.N. Theodorou, J. Phys. Chem. 98 (1994) 5111. [95] P. Demontis, G.B. Suffritti, Chem. Rev. 97 (1997) 2845. [96] E.B. Webb, G.S. Grest, Catal. Lett. 56 (1998) 95. [97] J.A. Martens, M. Tielen, P.A. Jacobs, J. Weitkamp, Zeolites 4 (1984) 98. [98] R.Q. Snurr, J. Kärger, J. Phys. Chem. B 101 (1997) 6469. [99] S. Jost, N.K. Bär, S. Fritzsche, R. Haberlandt, J. Kärger, J. Phys. Chem. B 102 (1998) 6375. [100] D.S. Sholl, K.A. Fichthorn, J. Chem. Phys. 107 (1997) 4384. [101] L.N. Gergidis, D.N. Theodorou, J. Phys. Chem. B 103 (1999) 3380. [102] S. Fritzsche, R. Haberlandt, J. Kärger, Z. Phys. Chemie 189 (1995) 211. [103] E.J. Maginn, A.T. Bell, D.N. Theodorou, J. Phys. Chem. 97 (1993) 4173. [104] M.J. Sanborn, R.Q. Snurr, Sep. Purif. Technol. 20 (2000) 1. [105] M.J. Sanborn, R.Q. Snurr, in preparation. [106] R.L. June, A.T. Bell, D.N. Theodorou, J. Phys. Chem. 95 (1991) 8866. [107] R.Q. Snurr, A.T. Bell, D.N. Theodorou, J. Phys. Chem. 98 (1994) 11948. [108] F. Jousse, S.M. Auerbach, J. Chem. Phys. 107 (1997) 9629. [109] C. Tunca, D.M. Ford, J. Chem. Phys. 111 (1999) 2751. [110] D.S. Sholl, C.K. Lee, J. Chem. Phys. 112 (2000) 817. [111] A. Gupta, R.Q. Snurr, in preparation. [112] E.J. Maginn, A.T. Bell, D.N. Theodorou, J. Phys. Chem. 100 (1996) 7155. [113] J. Kärger, M. Petzold, H. Pfeifer, S. Ernst, J. Weitkamp, J. Catal. 136 (1992) 283. [114] B. Frank, K. Dahlke, G. Emig, Chem. Ing. Tech. 12 (1992) 1104. [115] C. Rödenbeck, J. Kärger, K. Hahn, J. Catal. 157 (1995) 656. [116] B.L. Trout, A.K. Chakraborty, A.T. Bell, Chem. Eng. Sci. 52 (1997) 2265. [117] E. Klemm, J. Wang, G. Emig, Chem. Eng. Sci. 52 (1997) 3173. [118] M.S. Okino, R.Q. Snurr, H.H. Kung, J.E. Ochs, M.L. Mavrovouniotis, J. Chem. Phys. 111 (1999) 2210. [119] J.A. Dumesic, D.F. Rudd, L.M. Aparicio, J.E. Rekoske, A.A. Trevino, The Microkinetics of Heterogeneous Catalysis, American Chemical Society, Washington, DC, 1993.

L.J. Broadbelt, R.Q. Snurr / Applied Catalysis A: General 200 (2000) 23–46 [120] L.R. Petzold, DASSL Differential/Algebraic System Solver, Sandia National Laboratories, Livermore, CA, 1983. [121] S.F. Bush, P. Dyer, Proc. R. Soc. London A 351 (1976) 33. [122] A.S. McLeod, Catal. Today 53 (1999) 159. [123] J.E. Rekoske, R.J. Madon, L.M. Aparicio, J.A. Dumesic, Stud. Surf. Sci. Catal. 75 (1993) 1653. [124] G. Yaluris, J.E. Rekoske, L.M. Aparicio, R.J. Madon, J.A. Dumesic, J. Catal. 153 (1995) 54. [125] G. Yaluris, J.E. Rekoske, L.M. Aparicio, R.J. Madon, J.A. Dumesic, J. Catal. 153 (1995) 65. [126] G. Yaluris, R.J. Madon, D.F. Rudd, J.A. Dumesic, Ind. Eng. Chem. Res. 33 (1994) 2913. [127] G. Yaluris, R.J. Madon, J.A. Dumesic, J. Catal. 165 (1997) 205. [128] G. Yaluris, R.J. Madon, J.A. Dumesic, J. Catal. 186 (1999) 134. [129] B.A. Watson, M.T. Klein, Ind. Eng. Chem. Res. 35 (1996) 1506. [130] B.A. Watson, M.T. Klein, R.H. Harding, Int. J. Chem. Kinet. 29 (1997) 545. [131] B.A. Watson, M.T. Klein, R.H. Harding, Ind. Eng. Chem. Res. 36 (1997) 2954. [132] B.A. Watson, M.T. Klein, R.H. Harding, Energy & Fuels 11 (1997) 354. [133] B.A. Watson, M.T. Klein, R.H. Harding, Appl. Catal. A: Gen. 160 (1997) 13. [134] R.D. Cortright, J.A. Dumesic, J. Catal. 148 (1994) 771. [135] R.D. Cortright, J.A. Dumesic, J. Catal. 157 (1995) 576. [136] R.D. Cortright, E. Bergene, P. Levin, M. Natal-Santiago, J.A. Dumesic, Stud. Surf. Sci. Catal. 101 (1996) 1185. [137] S.A. Goddard, M.D. Amiridis, J.E. Rekoske, N. CardonaMartinez, J.A. Dumesic, J. Catal. 117 (1989) 155. [138] R.D. Cortright, S.A. Goddard, J.E. Rekoske, J.A. Dumesic, J. Catal. 127 (1991) 342. [139] J.E. Rekoske, R.D. Cortright, S.A. Goddard, S.B. Sharma, J.A. Dumesic, J. Phys. Chem. 96 (1992) 1880. [140] S.A. Goddard, R.D. Cortright, J.A. Dumesic, J. Catal. 137 (1992) 186. [141] M.D. Amiridis, J.E. Rekoske, J.A. Dumesic, D.F. Rudd, N.D. Spencer, C.J. Pereira, AIChE J. 37 (1991) 87. [142] L.M. Aparicio, S.A. Rossini, D.G. Sanfilippo, J.E. Rekoske, A.A. Trevino, J.A. Dumesic, Ind. Eng. Chem. Res. 30 (1991) 2114. [143] D. Wolf, Catal. Lett. 27 (1994) 207. [144] I. Alstrup, J. Catal. 151 (1995) 216. [145] L.M. Aparicio, J. Catal. 165 (1997) 262. [146] D.J. Klinke II, L.J. Broadbelt, Chem. Eng. Sci. 54 (1999) 3379. [147] D.J. Dooling, J.E. Rekoske, L.J. Broadbelt, Langmuir 15 (1999) 5846. [148] J.A. Dumesic, N.Y. Topsoe, T. Slabiak, P. Morsing, B.S. Clausen, E. Tornqvist, H. Topsoe, M. Sinev, G.V. Vanommen, M. Deboer, L. Holmid, G.C. Bond, L. Vonhippel, P. Forzatti, J. Armor, I.S. Nam, A.T. Bell, J. Blancon, Stud. Surf. Sci. Catal. 75 (1993) 1325. [149] C.V. Ovesen, B.S. Clausen, B.S. Hammershoi, G. Steffensen, T. Askgaard, I. Chorkendorff, J.K. Norskov, P.B. Rasmussen, P. Stolze, P. Taylor, J. Catal. 158 (1996) 170.

45

[150] C.R.F. Lund, Ind. Eng. Chem. Res. 35 (1996) 3067. [151] C.R.F. Lund, Ind. Eng. Chem. Res. 35 (1996) 2531. [152] E. Tserpe, K.C. Waugh, Stud. Surf. Sci. Catal. 109 (1997) 401. [153] L.M. Aparicio, J.A. Dumesic, Top. Catal. 1 (1994) 233. [154] M. Muhler, F. Rosowski, G. Ertl, Catal. Lett. 24 (1994) 317. [155] M. Muhler, R. Rosowski, O. Hinrichsen, A. Hornung, Stud. Surf. Sci. Catal. 101 (1996) 317. [156] F. Rosowski, O. Hinrichsen, M. Muhler, G. Ertl, Catal. Lett. 36 (1996) 229. [157] O. Hinrichsen, F. Rosowski, M. Muhler, G. Ertl, Chem. Eng. Sci. 51 (1996) 1683. [158] O. Hinrichsen, F. Rosowski, A. Hornung, M. Muhler, G. Ertl, J. Catal. 165 (1997) 33. [159] O. Hinrichsen, F. Rosowski, M. Muhler, G. Ertl, Stud. Surf. Sci. Catal. 109 (1997) 389. [160] O. Hinrichsen, Catal. Today 53 (1999) 177. [161] B. Fastrup, J. Catal. 168 (1997) 235. [162] K.C. Waugh, Catal. Today 53 (1999) 161. [163] J. Sehested, C.J.H. Jacobsen, E. Tornqvist, S. Rokni, P. Stoltze, J. Catal. 188 (1999) 83. [164] M. Bowker, I. Parker, K.C. Waugh, Surf. Sci. 197 (1988) L223. [165] M. Bowker, I.B. Parker, K.C. Waugh, Appl. Catal. 14 (1985) 101. [166] P. Stoltze, J.K. Nørskov, J. Catal. 110 (1988) 1. [167] P. Stoltze, J.K. Nørskov, Phys. Rev. Lett. 55 (1985) 2502. [168] L.J. Broadbelt, J.E. Rekoske, Chem. Eng. Sci. 51 (1996) 3337. [169] L.J. Broadbelt, J.E. Rekoske, in: G.F. Froment, K.C. Waugh (Eds.), Dynamics of Surfaces and Reaction Kinetics in Heterogeneous Catalysis, Vol. 109, Elsevier, Amsterdam, 1997, p. 341. [170] M.G. Evans, M. Polanyi, Trans. Faraday Soc. 34 (1938) 11. [171] D.W. Goodman, R.D. Kelley, T.E. Madey, J.T. Yates, J. Catal. 63 (1980) 226. [172] J. Xu, G.F. Froment, AIChE J. 35 (1989) 88. [173] M. Boudart, G. Djega-Mariadassou, Kinetics of Heterogeneous Catalytic Reactions, Princeton University Press, Princeton, NJ, 1984. [174] M. Boudart, Ind. Eng. Chem. Res. 28 (1989) 379. [175] M. Boudart, AIChE J. 2 (1956) 62. [176] L.J. Broadbelt, S.M. Stark, M.T. Klein, Chem. Eng. Sci. 49 (1994) 4991. [177] L.J. Broadbelt, S.M. Stark, M.T. Klein, Ind. Eng. Chem. Res. 33 (1994) 790. [178] L.J. Broadbelt, S.M. Stark, M.T. Klein, Comput. Chem. Eng. 20 (1996) 113. [179] J.P. Hindermann, G.J. Hutchings, A. Kiennemann, Catal. Rev.-Sci. Eng. 35 (1993) 1. [180] G.P. van der Laan, A.A.C.M. Beenackers, Catal. Rev.-Sci. Eng. 41 (1999) 255. [181] E. Shustorovich, Surf. Sci. 163 (1985) L730. [182] E. Shustorovich, Surf. Sci. Rep. 6 (1986) 1. [183] E. Shustorovich, Surf. Sci. 181 (1987) L205. [184] E. Shustorovich, J. Mol. Catal. 54 (1989) 301. [185] E. Shustorovich, A.T. Bell, Surf. Sci. 248 (1991) 359.

46

L.J. Broadbelt, R.Q. Snurr / Applied Catalysis A: General 200 (2000) 23–46

[186] J.J.C. Geerlings, M.C. Zonnevylle, C.P.M. de Groot, Surf. Sci. 241 (1991) 302. [187] Y.K. Park, P. Aghalayam, D.G. Vlachos, J. Phys. Chem. A 103 (1999) 8101. [188] A.S. McLeod, Catal. Today 53 (1999) 289. [189] H. Persson, P. Thormahlen, V.P. Zhdanov, B. Kasemo, Catal. Today 53 (1999) 273. [190] D.E. Ellis, D. Guenzburger, Adv. Quantum Chem. 34 (1999) 51. [191] M.E. Tuckerman, B.J. Berne, G.J. Martyna, J. Chem. Phys. 97 (1992) 1990. [192] A. Voter, Phys. Rev. Lett. 78 (1997) 3908.

[193] A. Voter, J. Chem. Phys. 106 (1997) 4665. [194] W. Mortier, W. Ghosh, S. Shankar, J. Am. Chem. Soc. 108 (1986) 4315. [195] K. Van Genechten, W. Mortier, P. Geerlings, J. Chem. Phys. 86 (1987) 5063. [196] R. Parr, W. Yang, J. Am. Chem. Soc. 106 (1984) 4049. [197] M.W. Deem, AIChE J. 44 (1998) 2569. [198] D.J. Dooling, R.J. Nielsen, L.J. Broadbelt, Chem. Eng. Sci. 54 (1999) 3399. [199] G.F. Froment, K.B. Bischoff, Chemical Reactor Analysis and Design, Wiley, New York, 1990.