1174
Ind. Eng. Chem. Res. 2003, 42, 1174-1183
Parameter Optimization of Molecular Models: Application to Surface Kinetics S. Raimondeau,†,‡ P. Aghalayam,†,§ A. B. Mhadeshwar,| and D. G. Vlachos*,| Department of Chemical Engineering, University of Massachusetts-Amherst, Amherst, Massachusetts 01003, and Department of Chemical Engineering and Center for Catalytic Science and Technology, University of Delaware, Newark, Delaware 19716
A multistep methodology is explored for the optimization of the parameters of molecular models. This methodology is based on response surface methods and circumvents the practically impossible rigorous parameter fitting of computationally intensive molecular models. The approach is demonstrated for the parameters of the catalytic oxidation reaction of CO on platinum modeled by kinetic Monte Carlo simulations. After implementation of an initial reaction mechanism, sensitivity analysis is carried out to determine the key (active) reaction parameters for selected experimental points. Next, solution mapping is used to parametrize the model responses as low-degree polynomials of the active parameters. Finally, optimization of the active parameters is performed using simulated annealing. The optimized parameters are contrasted to those of a continuum-type mean-field model. The effect of surface diffusion in determining intrinsic kinetic parameters is also addressed, and the possibility of bridging the pressure gap via multiscale simulations is discussed. Introduction Molecular models, such as molecular dynamics and Monte Carlo (MC) methods, are widely used to gain insight into complex thermodynamic and dynamic behavior in a number of applications ranging from vaporliquid equilibrium, to adsorption and diffusion in microporous materials and on surfaces, to crystal growth, to surface reactions.1-6 Molecular models often involve large numbers of input parameters, for example, in the intermolecular potential of hydrocarbons in zeolites or in a large kinetic mechanism.1,7 There is currently a need to incorporate realistic parameters into these models so further quantitative model predictions and guidance of experiments can be realized. These different model parameters are usually fitted to quantum mechanical calculations or experimental data by a trialand-error procedure until an adequate fit is achieved. Because of the significant computational cost of molecular models, it is impractical to rigorously optimize their parameters and obtain confidence intervals because of the large number (which are often many thousands) of function evaluations, i.e., molecular simulations, needed in the optimization process. A further complication, specific to surface reactions, in optimizing model parameters against experimental data entails the inability of MC simulations to treat very fast surface diffusion processes. Because surface diffusion is not rate-determining under all conditions, it is desirable to * Corresponding author. Tel.: (302) 831-2830. Fax: (302) 831-2085. E-mail:
[email protected]. † University of Massachusetts-Amherst. ‡ Current address: ExxonMobil Research and Engineering Company, 3225 Gallows Rd., Room 7A0671, Fairfax, VA 22037-0001. § Current address: Department of Chemical Engineering, Indian Institute of Technology-Bombay, Powai, Mumbai 400076, India. | University of Delaware.
systematically identify conditions under which surface diffusion is important in extracting intrinsic surface kinetics. To overcome these difficulties, here, we explore a mathematically systematic multistep approach, which has been previously employed in continuum scale models only,8,9 to optimize parameters of molecular models. This approach is illustrated in the catalytic oxidation of CO on a platinum surface modeled using kinetic MC simulations. The choice of this model reaction is driven in part by its importance in the automotive catalytic converter and the fact that it is a subset of more complex reaction systems, such as the oxidation of natural gas to syngas. However, the ideas discussed here can be applied to more complex surface reaction mechanisms and any other molecular-scale model. A comparison of optimized parameters of conventional mean-field (MF) models and MC simulations is also presented to assess the extent to which optimized reaction rate constants of MF models can differ from those of kinetic MC simulations. Finally, the role of microscopic inhomogeneities in extracting intrinsic kinetics and the possibility of computer-assisted bridging of the so-called pressure gap are discussed. Models and Experimental Data The steps in the catalytic oxidation of CO on platinum are listed in Table 1. They include the reversible adsorption of CO, the reversible dissociative adsorption of O2, and the bimolecular surface reaction between CO* and O*. The adsorption rate constants are computed using the kinetic theory of ideal gases.8 The desorption and reaction rate constants are computed using an Arrhenius expression.8 Preexponentials and sticking coefficients are hereafter termed prefactors. The optimization is illustrated for two models, namely, a continuum-type, deterministic MF theory and stochastic, kinetic MC simulations. Deterministic Mean-Field Model. In MF models, species are assumed to be distributed homogeneously
10.1021/ie0202470 CCC: $25.00 © 2003 American Chemical Society Published on Web 02/26/2003
Ind. Eng. Chem. Res., Vol. 42, No. 6, 2003 1175 Table 1. Steps in the Catalytic Oxidation of CO along with Transition Probabilities Used in the MC Simulations, Initial Kinetic Parameters, MF-Optimized Kinetic Prefactors, and MC-Optimized Prefactors in the Absence and Presence of Surface Diffusiona
reaction
transition probabilities
Ao (s-1) or so14
Ea14 (kcal/mol)
MF-optimized prefactor
MC-optimized prefactor without diffusion
MC-optimized prefactor with diffusion
1. CO + * f CO*
Γˆ aCO ) kaCOPyCOP*
0.89
2. CO* f CO + *
Γˆ dCO ) kdCOPCO*
9.40 × 1016
3. O2 + 2* f 2O*
a a Γˆ O ) kO PyO2P*P*/* 2 2
0.10
4. 2O* f O2 + 2*
d d ) kO PO*PO*/O* Γˆ O 2 2 Γˆ r ) kr(PCO*PO*/CO* + PO*PCO*/O*)
1.00 × 1013
51.0
-
0.15c -
0.89 [0.81, 0.95]b 0.92c 5.82 × 1017 d [4.9, 7.2] × 1017 b 4.50 × 1017 c 0.09 [0.088, 0.097]b 0.11c -
4.90 × 109
11.0
-
-
-
5. CO* + O* f CO2
0.0
44.0
0.0
0.84 [0.82, 0.87]b 3.46 × 1017 d [3.0, 4.0] × 1017 b 0.09 [0.085, 0.093]b
0.91 0.87c 2.01 × 1017 d 1.96 × 1017 c 0.15
a P denotes the probability of picking species x and P x x/y is the conditional probability of picking species y knowing that species x has been chosen (for modeling details, see ref 17). The prefactors of steps 4 and 5 are left unoptimized. The adsorbate-adsorbate interaction b in CO desorption is taken to be 5 kcal/mol. 85% confidence interval for the optimized prefactors for the case of complete factorial design. c Prefactors using the SAB method. d Preexponentials A are in s-1. o
through space. Such models hold when the diffusion of all species is Fickian and infinitely fast. The MF equations are10 (note that the desorption of O* is included and, for simplicity, the saturation values are taken to be 1)
dθCO* ) kaCOyCO[(1 - (θCO*)2)] - kdCOθCO* - krθCO*θO* dt (1) dθO* a d 2 ) 2kO y θ 2 - 2kO θ - krθCO*θO* 2 O * 2 O* dt
(2)
θCO* + θO* + θ* ) 1
(3)
where yi is the mole fraction of gaseous species i; k is the rate constant for each step, θi is the surface coverage of species i; and the superscripts a, d, and r denote adsorption, desorption, and reaction, respectively. The dependence of the CO sticking coefficient on the CO* surface coverage results from a precursor-type adsorption process11,12 and has been employed in recent simulations in different ways. Lateral adsorbate-adsorbate interactions lead to the multiplying factor eωθCO*/RT of kd, where ω is the repulsive interaction energy (for its value, see Table 1). A coverage dependence of the activation energy of surface reaction step can also be included in simulations, for example, by employing the bond-order conservation (BOC) theory,13 as has been done in continuum models.14 However, because surface reaction is not rate-determining under the conditions we study, we ignore such a dependence, as it plays no role. At steady state, the MF equations are solved using Newton’s technique. It has been recognized for a long time that the layer of adsorbates can be spatially nonuniform. Microscopicscale spatial inhomogeneities can arise whenever a nonlinear term appears in the corresponding MF equations, i.e., whenever the rates of microscopic processes depend on the local environment in a nonlinear manner. MF models are then inadequate to quantitatively describe the system response such as reaction rates. Discrete particle methods, such as kinetic MC simulations, are ideally suited to handle such inhomogeneities.
The work of Ziff and co-workers has catalyzed developments in this area.15 Such simulations are discussed next. Kinetic Monte Carlo Simulations. There are various types of MC simulations. When species are relatively strongly adsorbed on a surface, i.e., the various activated processes are relatively slow compared to the vibrational time, species reside in local minima of the potential energy surface that can be approximated as points on an appropriate lattice (lattice MC simulations). Under these circumstances, activated processes happen infrequently (rare-event dynamics) in comparison to the vibrational time, with certain transition probabilities per unit time that, in general, vary with time and in space. MC simulations are stochastic in nature and follow the history of individual species in a Markov sense by comparing normalized transition probabilities to a random number. Because of their discrete nature, MC simulations are very intensive, and as a result, they are limited to short length and time scales. The specific lattice MC algorithms used in the modeling of surface reactions and growth processes belong to the general class of dynamic MC simulations. In particular, in the continuous time or kinetic MC method,16 first, the probabilities of all microscopic processes (adsorptions, desorptions, surface reactions, etc.) are computed. Next, a microscopic process is randomly chosen on the basis of its transition probability. Finally, a site or a pair of sites is randomly picked from a list of sites or pairs associated with the selected process, and the event is executed. The efficiency and accuracy of kinetic MC algorithms for surface reactions is discussed in ref 17. For the specific problem addressed here, the transition probabilities per unit time for each process are listed in the second column of Table 1 (more general cases are described in ref 17). For CO-CO lateral firstnearest-neighbor interactions, kd should be multiplied by ∑surfaceewCO-COnco/RT normalized by the number of occupied sites, where nco is the number of CO first nearest neighbors of each CO species and wCO-CO is the repulsion interaction energy per pair. More details in treating interactions are presented in refs 17 and 18. Each MC simulation is split into windows, i.e., intervals consisting of a certain number of MC events.
1176
Ind. Eng. Chem. Res., Vol. 42, No. 6, 2003
During each window, bookkeeping is done so that spatiotemporal averages of reaction rates and coverages can be computed at the end of a window, as proposed in ref 17. The number of MC events in each window controls the level of noise and the apparent output discretization accuracy in the time domain. This number should not be large to avoid significant coarse graining of the dynamics, but the window should contain a certain number of events from each process so that reasonable statistics can be collected. Windows consisting of 102-104 MC events typically provide good statistics. Of course, this generalization depends on the disparity of time scales between the microscopic processes simulated. A simulation is run until steady state is reached. Two criteria are used to assess whether steady state has been reached as described in refs 17 and 18. The first criterion examines mass conservation for each species, i.e., whether a normalized residual of each material balance is smaller than some low tolerance
|Ran - Rdn - Rrn|/Ran <
(4)
where R represents the rate; n is the window number; and the subscripts a, d, and r refer to adsorption, desorption, and reaction, respectively. Because the actual reaction rates can be large or small, we normalize the rates to the rate of adsorption of the corresponding species. Regarding the steady-state tolerance, for a lattice of 120 × 120 and the above window size, steadystate residuals on the order of 10-2-10-3 are typical. The second criterion examines the existence of an inflection point in spatiotemporal average coverages and rates between successive windows. For the spatiotemporal average X of a given quantity, one examines whether
(Xn - Xn+1)(Xn+1 - Xn+2) < 0
(5)
Once both steady-state criteria have been met, a longer window of MC events is used to compute spatially averaged steady-state reaction rates and coverages. Most kinetic MC simulations dealing with surface reactions have been conducted in the absence of surface diffusion. In systems of interest, surface diffusion exhibits a significantly lower activation energy than desorption, causing a disparity in time scales that renders MC simulations with realistic parameters impractical.19 We return to this point in the section Role of Surface Diffusion in Selecting Target Experimental Data Points. A significant fraction of simulations that have included surface diffusion (albeit slow compared to real systems) have addressed the case of Fickian diffusion, i.e., the case where diffusion is unaffected by the presence of intermolecular forces and causes mixing of adsorbates. Here, we treat the case of non-Fickian surface diffusion, which is probably more realistic for interacting species. Surface diffusion is modeled following Arrhenius kinetics, i.e., the activation energy depends on the binding energy of the departing site alone, as detailed in ref 20. The activation energy consists of a zerocoverage term plus a term that accounts for lateral interactions using the typical pairwise-additivity assumption. Lateral interactions are taken equal to those of desorption. Finally, the vibrational frequency of surface diffusion jumps is taken as equal to that of desorption.
In this model, a square, two-dimensional 120 × 120 lattice of active sites is considered (unless otherwise specified). Even though this facet does not precisely represent the experimental crystallographic plane discussed below, it is sufficient to demonstrate the proposed computational approach. Similar model simplification was followed in ref 21. Experimental Data. The kinetic parameters are refined here using the experiments of ref 22 carried out under ultrahigh-vacuum conditions on a Pt(111) surface with a fixed total flux of 1.65 × 1014 molecules/(cm2 s) at various surface temperatures and a CO/O2 ratio varying from 0 to 1. These data are shown in Figure 1a. Initial Estimation of Parameters and Model Predictions For the optimization, initial values of model parameters are needed. Although intensive research has been conducted on CO oxidation on Pt and parameters have been published (albeit quite different in some cases), e.g., refs 23-25, we prefer a bottom-up multiscale approach. This multiscale approach enables one to construct reaction mechanisms for reaction systems more complex than CO for which limited information is available. We have recently demonstrated the power of this approach using mean-field models for hydrogen, carbon monoxide, and methane oxidation on platinum.7,8,14 Variations of the multiscale methodology have been reviewed in refs 19 and 26. For large reaction mechanisms, first-principles simulations become too intensive, so we rely on semiempirical methods with minimal input from experiments or quantum mechanics. The activation energies along with their dependence on coverage can be obtained relatively accurately using semiempirical methods, such as the BOC technique.13,27 The ability of the model to describe experiments over a wide range of temperature (because of the large variation in Boltzmann factors), without optimization, provides further credibility for our approach. Therefore, these parameters are not subject to optimization. Prefactors, such as sticking coefficients, are more difficult to obtain from first-principles simulations. Furthermore, less is known about the role of interactions in prefactors. For this reason, we choose some of the prefactors as the target parameters for optimization. Judging from our experience from the optimization of MF models,8,14 order-of-magnitude estimates of prefactors, based on the type of reaction as determined from transition state theory,28 can be sufficient. The initial parameters developed from this procedure optimized against different experimental data in ref 14 are listed in Table 1 and are used for optimization of parameters in MC simulations. Note that direct chemisorption of CO was assumed in ref 14 instead of the precursor type of adsorption used here. This modification was found to be essential for correctly predicting the slope (linearity of reaction rate with varying composition for fuellean mixtures) in Figure 1a. Figure 1a shows the reaction rates (in monolayers per second) predicted using both the MF model (solid line) and the MC model (filled circles connected with dashed line) versus the CO/O2 ratio using the initial prefactors (Table 1), along with the experimental data (open circles) at a surface temperature of 500 K. No surface diffusion is included in these MC simulations. The error
Ind. Eng. Chem. Res., Vol. 42, No. 6, 2003 1177
Figure 1. Reaction rate in monolayers per second versus the CO/O2 ratio for the MC model (filled circles connected with dashed line) and the MF model (solid curve) along with the experimental data (open circles) of ref 22 at a surface temperature of 500 K. Panels b-d show three snapshots (at conditions indicated by arrows) of the surface, which is mainly covered by O* under fuel-lean conditions and CO* at higher CO/O2 ratios.
bars for selected compositions indicate two standard deviations associated with the inherent noise of MC simulations. Although the qualitative trend in the experimental data is well captured by the initial set of prefactors (note that this set of experiments was not used in the optimization in ref 14), refinement of the kinetic parameters is needed to quantitatively describe these data. The peak in the rate often corresponds to a stoichiometric surface coverage of O* and CO*, marking the transition between a fuel-lean mixture with an O*covered surface and a fuel-rich mixture with a partially CO*-covered surface as discussed below. Figure 1b-d shows snapshots from MC simulations at the conditions indicated by the arrows and illustrate the inhomogeneous distribution of adsorbates. Note that, in MC simulations, the maximum does not correspond to equal coverages of the two adsorbates. Identification of Important Parameters An important step of the multistep optimization methodology is to determine the important (active) parameters affecting the model responses, i.e., the reaction rate in our case. This is an important step in developing polynomials for parameter optimization, described in the next section. This task can be achieved through a brute-force sensitivity analysis (SA). SA results also provide physical insight into the ratelimiting step(s) of a process. This step entails perturbation (e.g., reduction) of all parameters by a sufficiently large fraction and recalculation of model responses such
as reaction rates, coverages, etc., that are experimentally available.29 The model response chosen here is the rate of production of CO2. Although large perturbations cause measurable changes in model responses that are beyond model uncertainties, they could also shift the system into a different regime of parameter space, e.g., from fuel-lean to fuel-rich conditions. This can lead to wrong conclusions regarding the important parameters in each regime. A complication intrinsic to molecular models is related to the inherent noise present in MC simulations. For meaningful conclusions, the standard deviation in the response (here the reaction rate) should be smaller than the response change induced by perturbing the prefactors. Figure 2a shows the standard deviation in the reaction rate versus the inverse of the square root of the lattice size for the same number of MC steps (corresponding to the same real time) for a CO/O2 ratio of 0.3 at a surface temperature of 500 K. The standard deviation scales well with the inverse of the square root of the lattice size, and thus, the resulting SA performed on larger lattices is even more reliable for smaller perturbations. As an example, Figure 2b shows corresponding steady-state reaction rates versus time for the nominal (unperturbed) case and for a reduction of 30% in the CO desorption preexponential. Under these conditions, it can be seen that the mean values are distinct from each other. Such simple scaling results could indicate appropriate sizes for the perturbations in model parameters.
1178
Ind. Eng. Chem. Res., Vol. 42, No. 6, 2003
Figure 2. (a) Standard deviation (open circles) vs the inverse of the square root of the lattice size. The straight line is a leastsquares fit to the data. (b) Reaction rate for a CO/O2 ratio of 0.3 at a surface temperature of 500 K for the nominal case (solid line) and for a reduction of 30% in the CO desorption preexponential (dashed line). For panel b, the lattice size is 120 × 120.
The brute-force normalized sensitivity coefficient is defined as (Rperturbed - Rnominal)/0.3Rnominal, where Rperturbed is the value of the reaction rate recomputed after a 30% reduction in the prefactor and Rnominal is the nominal value of the reaction rate. Figure 3 summarizes the SA results for both (a) the MF and (b) the MC models. First it is seen that SA analysis, conventionally applied to continuum models, works very well for stochastic simulations as well. To compare the changes in model response and uncertainty in reaction rate, the vertical dotted area in Figure 3b indicates the standard deviation normalized by the nominal reaction rate and by the perturbation factor of 0.3. Any change within this area is within the uncertainty range. Physically, it is found that the reaction rate is practically independent of the O2 desorption and the bimolecular surface reaction between CO* and O*. Under fuel-lean conditions (left of the maximum reaction rate), the reaction rate is mainly controlled by the adsorption of CO because of the blocking of the surface by O*. For stoichiometric and fuel-rich mixtures, three kinetic parameters (L ) 3) influence the reaction rate in both models, namely, CO adsorption, CO desorption, and O2 adsorption. SA conducted at 600 K indicates that the reaction rate depends only on the adsorption of CO and O2. As CO desorption represents an important step under some conditions (lower temperatures), a temperature of 500 K was chosen to carry out the optimization. Numerically, an important observation is that the results of MF and MC simulations are comparable. Thus, SA in MF models of surface reactions, which can
Figure 3. Results of SA in reaction rate for a reduction of 30% in the prefactors at three different CO/O2 ratios (0.05, 0.3, and 0.6) and a surface temperature of 500 K for (a) the MF model and (b) the MC model. For fuel-lean mixtures, the adsorption of CO is the key step controlling the reaction rate. For higher CO contents, the adsorption of CO and O2 and the desorption of CO are important. The MF results are comparable to the MC ones. The shaded bar in panel b indicate the maximum uncertainty range in reaction rate among the three compositions (also see text for an exact definition).
be obtained efficiently by brute force or even more efficiently by linear algebra, can be sufficient for the optimization of MC parameters. Furthermore, the size of the perturbations that need to be introduced in kinetic MC simulations for sufficient changes in model responses compared to the uncertainty discussed above can be estimated using MF models, which are computationally much less expensive. Role of Surface Diffusion in Selecting Target Experimental Data Points. Target points from the available experimental data must be chosen, and SA must be performed at each point to identify the active parameters. From an optimization point of view, the number of points chosen must be larger than the number of parameters fitted. However, which points should be chosen is an open question. Here, we discuss how points can be chosen so that intrinsic kinetics (without surface transport effects and microscopic inhomogeneities) can be determined. Realistic surface diffusion cannot be rigorously modeled because of the disparity in time scales introduced by this microprocess. Indeed, the inclusion of realistic diffusion rates leads to insufficient MC sampling of the rest of the microprocesses. A number of methods have been proposed to overcome this deficiency, including partial equilibration of the surface by allowing a certain number of surface diffusion events after each successful microprocess, quasi-equilibrium models, a hybrid MF/ MC approach, and mesoscopic models (see ref 19 for a
Ind. Eng. Chem. Res., Vol. 42, No. 6, 2003 1179
Figure 4. (a) Surface coverages for MC and (b) ratio of adsorption to adsorption plus desorption for both models versus the CO/O2 ratio using the optimized parameters at surface temperatures of 500 K (lower curves) and 600 K (upper curves). For sufficiently fuel-lean mixtures, the surface is mainly covered with O*, so the reaction probability of each absorbing CO molecule is high. CO* is in partial equilibrium for sufficiently fuel-rich mixtures at 500 K. The shading indicates an approximate range where surface diffusion can be the rate-determining step. The MF model can also indicate partial equilibrium conditions. The inset in panel b shows the relative change in reaction rate upon inclusion of surface diffusion at the three CO/O2 ratios indicated at a surface temperature of 500 K. The zero-coverage activation energy for surface diffusion has been taken to be 38 kcal/mol.
review). We chose an approach closely related to the partial equilibration approach30,31 that is outlined below. To reduce computational cost, we propose that most (or even all) of the target points could be chosen at conditions where surface diffusion is not the ratedetermining step. One method of identifying the importance of surface diffusion is via SA, similar to the one discussed above. As an example, the inset of Figure 4b shows the relative change in reaction rate upon inclusion of surface diffusion for three compositions. This sensitivity analysis indicates that surface diffusion is very important only near the maximum reaction rate. Another economical method of determining the importance of surface diffusion is discussed next. Figure 4 shows the surface coverage (panel a) and the ratio of adsorption rate to the sum of the adsorption and desorption rates (panel b) versus the CO/O2 ratio for MC simulations carried out at a surface temperature of 500 K. For sufficiently fuel-lean mixtures, CO* surface diffusion should have no influence on the reaction rate, because the surface is mostly covered by O* and the probability of reaction of every adsorbing CO is high. For sufficiently fuel-rich mixtures, on the other hand, CO can be assumed to be in partial equilibrium, as the ratio of adsorption to adsorption plus desorption ap-
proaches 0.5 (Figure 4b). In this regime, the adsorption and desorption of CO are significantly faster than the surface reaction, as evidenced by the material balance of CO* (see eq 4) and the nearly equal adsorption and desorption rates. The fast adsorption and desorption processes give rise to a partial equilibrium in the adsorbed layer. As a result, the surface diffusion of CO*, an inherently equilibrium process, does not affect the overall reaction rate. At intermediate compositions, corresponding to the peak in the reaction rate, the partial equilibrium condition breaks down, and the CO* surface diffusion cannot be neglected. The shaded area of Figure 4b provides a rough guide as to the composition regime where surface diffusion is important. The results of the inset of Figure 4b, which are computationally more expensive, lead to the same conclusions, that is, surface diffusion cannot be neglected near the peak in reaction rate. Consequently, to obtain intrinsic kinetic parameters, one could choose the target points away from the maximum reaction rate. Figure 4b indicates that the MF model can again be a rough guide for determining partial equilibrium and potentially the role of surface diffusion (this might appear contradictory to the infinite Fickian diffusion assumption of MF models). Note that, at higher temperatures, e.g., 600 K, the partial equilibrium assumption is less well valid, even for fuel-rich conditions, as shown in Figure 4b. Therefore, at higher temperatures, surface diffusion plays a more dominant role over a wider composition range. In our simulations, to accurately capture experiments performed at higher temperatures where the system does not reach partial equilibrium and explore how one can extract intrinsic kinetics even when surface diffusion is important, one target point near the maximum rate is included in the optimization. MC simulations with surface diffusion are conducted by gradually reducing the zero-coverage activation energy (starting from that of desorption) for the CO* surface diffusion step (O* diffusion is neglected). In the absence of surface patterns, we have found that a plateau in the reaction rate is reached, when surface diffusion is sufficiently faster (but still slow compared to experimental values) than the rest of the processes. An example is depicted in Figure 5a. Two snapshots of the surface are also shown without (panel b) and with (panel c) surface diffusion. In the presence of sufficiently fast surface diffusion, CO* diffuses and reacts with O*, leading to higher reaction rates and more vacant sites as observed in panel c. Optimization Procedure and Validation Here, we discuss the optimization procedure that makes estimation of parameters of molecular models a feasible task. Then, we compare model predictions to experimental data, and finally, we discuss differences between MC and MF parameters, i.e., the role of heterogeneities in extracting intrinsic kinetics, as well as the effect of surface diffusion on parameter estimation. Solution Mapping Methodology. To overcome the significant computational cost of direct optimization of parameters of molecular models, we propose to use a solution mapping or response surface technique. Response surface methods are used to approximate multidimensional surfaces with low-degree polynomials, which are subsequently used for parameter optimization.32 The basic idea is to develop an approximate but
1180
Ind. Eng. Chem. Res., Vol. 42, No. 6, 2003
where Y is the model response (rate of CO2 production in the considered case), L is the number of active parameters (L ) 2 or 3 in the considered case, depending on the CO/O2 ratio), the a’s are the coefficients in the polynomial, and the x’s are the active parameters (rate constants k) normalized by the range of sampling f using the following relations
Figure 5. (a) Reaction rate versus the zero-coverage activation energy of CO* surface diffusion at a CO/O2 ratio of 0.2. The plateau in reaction rate indicates sufficiently fast non-Fickian surface diffusion. Snapshots of the surface corresponding to a MC simulation (b) with no surface diffusion and (c) with non-Fickian surface diffusion.
relatively accurate model (a response surface), which represents the solution of the molecular model, but using a limited number of simulations. This approximate model is then subject to rigorous optimization techniques. Because of its simplicity, the approximate model overcomes the problem of expensive function evaluations encountered in the direct use of molecular models. Furthermore, various simulation programs or packages can be used, making analysis of very different data (e.g., NMR data for self-diffusivity and adsorption isotherms) tractable by creating a joint objective function for optimization. This approach has proved to be robust and very efficient in the optimization of parameters of large reaction mechanisms using continuum models.8,9 In particular, we parametrize the reaction rates as low-degree polynomials of the active parameters determined through the SA. The polynomials are functions of the active parameters and are given as L
ln Y ) a0 +
L
L
aixi + ∑∑aijxixj + ... ∑ i)1 i)1 jgi
(6)
xi )
ln(ki/k0i) ln fi
(7)
fi )
x
(8)
k0i ) xkmax,ikmin,i
(9)
kmax,i kmin,i
Here, kmax,i, kmin,i, and k0i represent the maximum, minimum, and nominal values, respectively, of the ith active parameter. The upper and lower bounds of the uncertainty correspond to x ) +1 and x ) -1, respectively, and the nominal value corresponds to x ) 0. One such polynomial is developed for each target experimental point. In particular, five different CO/O2 ratios (0.05, 0.2, 0.3, 0.6, and 1.0) are considered. Subsequently, factorial design is used to generate, in an optimal manner, sampling points in a window of parameter space around the initial values listed in Table 1.32 As a general rule, an order of magnitude around the nominal value is used as the range of sampling for preexponentials and sticking coefficients, while ensuring that the maximum sticking coefficient does not exceed 1. At these points, a small number of computations are performed to determine the model responses (vector Y). A second-degree polynomial is used here to represent the model responses derived from a 3L complete factorial design template, i.e., the reaction rate as a function of active parameters, because it yields satisfactory results. This template yields a factorial design matrix Z discussed in more detail in ref 32. Using this factorial design template, the associated coefficients (vector B) for a second-degree polynomial are determined using B ) (ZTZ)-1ZTY. The computational time needed to derive these polynomials in the absence of surface diffusion is about 0.5 h per run (for good statistics) times 27 runs, i.e., about 15 h of CPU time on a single processor (Alpha 533 MHz) are required. On the other hand, in the presence of nonFickian surface diffusion, using a 33 template requires about 10 days of CPU time on a single processor. These times should be compared to a few years needed to perform a direct optimization using the MC model for function evaluation, assuming 105 function evaluations. Once the reaction rate has been accurately parametrized by a polynomial, simulated annealing (a stochastic optimization technique33) is used to minimize the distance of the response surfaces from the experimental data, according to N
f)
Wi(Yiexp - Yimodel)2 ∑ i)1
(10)
i are the experimental points and where Yiexp and Ymodel response surface predictions, respectively; Wi is the weight for each set (here taken as 1); and N is the number of experimental points considered. The simulated annealing MC simulation should not be confused
Ind. Eng. Chem. Res., Vol. 42, No. 6, 2003 1181
ai ) si ln fi )
ln Y|R - ln Y|-R 2R
(12)
where R indicates the size of the perturbation in xi (0 e R e 1). The coefficients of the second-order terms are given as
aii ) aij )
Figure 6. Optimized reaction rates versus the CO/O2 ratio for the MF model (solid lines) and the MC model (dotted lines), along with experimental data for the three different surface temperatures indicated. Optimization was performed at 500 K only. The predictions are in good agreement with the experimental results at 400 and 600 K. The number of MC simulations is comparable to that in Figures 1-3. Surface diffusion is included in all MC simulations at 600 K, where partial equilibrium breaks down, but only near the maximum at lower temperatures. The dashed line at 500 K shows MF predictions obtained using the optimized parameters of the MC model.
with the kinetic MC simulations used for reaction modeling. Note that the uncertainty in the experimental points and/or model predictions (due to the inherent MC noise) can be taken into account by using nonequal weights. Finally, as mentioned above, multiple experimental data can be analyzed by simply adding a term for each set of experiments to eq 10. The optimized parameters along with confidence intervals are listed in Table 1. Figure 6 illustrates the reaction rate predicted using the optimized set of kinetic parameters along with the experimental points (symbols) versus the CO/O2 ratio at the three different temperatures indicated. At 500 K (targeted temperature), the maximum rate is overpredicted by the MF model and less so by the MC model, but overall, the experimental data are well captured by both models. Model validation at 400 and 600 K, depicted in the same graph, shows good agreement between the experimental data and the predictions of the optimized models. Alternative Response Surface Methods. In addition to complete factorial design used above with three major levels, other surface response surface techniques can be used, including two major levels, fractional factorial design, etc. Here, we demonstrate a recently introduced sensitivity-analysis-based (SAB) method in developing low-degree polynomials. In the SAB method, the coefficients of the polynomials are determined from the normalized sensitivity analysis coefficients (see ref 34 for details). In particular
a0 ) ln Y(0)
(11)
The first-order terms are obtained from first-order sensitivity coefficients, commonly defined as si ) ∂ ln Y/∂ ln ki. Sensitivity coefficients are calculated either locally or by finite difference. In the considered case, we follow a central finite difference approach, which gives
(
( )
) (
1 si,i - s-i,i ln fi 2 2R
(13)
)
1 si,j - s-i,j 1 sj,i - s-j,i ln fj + ln fi 2 2R 2 2R
(14)
where s(i,j ) ∂ ln Y(xi)(R)/∂ ln kj gives the second-order sensitivity coefficients obtained by finite differencing the first-order sensitivity coefficients s. Similarly, the third-order coefficients are calculated (if desired) as
( (
) )
aiii )
1 si,i - 2s0,i + s-i,i ln fi 6 R2
(15)
aiij )
1 si,j - 2s0,j + s-i,j ln fj 2 R2
(16)
The SAB method was shown to lead to third-degree polynomials with reduced computational cost compared to complete factorial design when local sensitivity analysis can be used.34 The active parameters (in the range from -1 to +1) are perturbed by (R (R ) 0.1) to calculate the sensitivity analysis matrix. Aside from the possibility of reducing the computational cost, what makes this approach interesting is that the sensitivity matrix computed using MC simulations is noisy. Therefore, the success of the SAB technique is not a priori clear. As an illustrating example, results obtained in the absence of surface diffusion using both the complete factorial design and the SAB method are shown in Figure 7, and the optimized parameters in the absence and presence of surface diffusion are summarized in Table 1. Comparable fits between the two surface response methods are depicted for all cases. Because a brute-force sensitivity analysis is employed in the SAB method, the number of numerical simulations is (2L + 1)2 compared to 2L + 1 described in ref 34 for systems where local sensitivity analysis is feasible. This is because, at all 2L + 1 points in the hypercube, an extra (2L) × (2L + 1) simulations are needed to generate SA coefficients, thus, resulting in a total of (2L + 1)2 simulations for L parameters. This number can be halved if forward finite difference is used instead. In contrast, for a complete factorial design, depending on the number of major levels, either 3L (three major levels, viz., -1, 0, and +1) or 2L + 2L + 1 (two major levels, viz., -1 and +1) simulations are required. In fractional factorial design, 2L-M + 2L + 1 simulations are needed. In the particular case of CO oxidation, where only three parameters (L ) 3) are important for fuel-rich conditions, 49 numerical simulations are required for the SAB method for each target point, whereas for the complete factorial design method, 27 numerical simulations (3L) or 15 numerical simulations (2L + 2L + 1) are sufficient to yield second-degree polynomials. For this specific case, the SAB method is computationally more expensive than the factorial design methods. However, a crossover occurs at L ) 9,
1182
Ind. Eng. Chem. Res., Vol. 42, No. 6, 2003
scale simulations along with the approach described herein can be used to extract intrinsic kinetics that, in turn, can be used in very different operating conditions. In this regard, the approach proposed here might help one to bridge the pressure gap between surface science studies and industrial conditions of catalytic processes. Conclusions
Figure 7. Optimized reaction rates versus the CO/O2 ratio using the MC model, along with experimental data for a surface temperature of 500 K. Predictions using the complete factorial design method and the SAB method are compared in the absence of surface diffusion.
where the complete factorial design two-level method becomes computationally more expensive than the SAB method with central finite difference. Therefore, the SAB method is certainly promising when molecular models with large numbers of parameters are to be optimized. Role of Inhomogeneities and Surface Diffusion in Extracting Intrinsic Kinetics The optimized parameters of the MF and MC models are close enough, and only the preexponential for CO desorption appears to be statistically different. The implications of this finding are that an initial optimization of active parameters using an MF model can reduce computational cost by guiding MC simulations, for example, in picking the factorial design window. Despite the small difference in optimized parameters, noticeable differences in predictions are found when the refined reaction parameters derived using the MF model are used in the MC model (and vice versa) at 500 K. As an example, the dashed line in Figure 6 at 500 K shows model predictions obtained with the MF model using the optimized parameters of the MC simulator. However, at 600 K, where CO desorption is not a key step, the two optimized mechanisms yield practically the same curve (not shown). Simulations away from the maximum reaction rate have confirmed that surface diffusion plays a minor role when the surface is in partial equilibrium. Therefore, one should expect that MF-optimized parameters will be only a reasonable starting point of MC simulators, but will lack quantitative predictive capability. Finally, in comparing kinetic parameters when surface diffusion is left out (taken as zero) to those where surface diffusion is included, it is clear that surface diffusion plays a role in extracting kinetic parameters. The better prediction of experimental data near the maximum, seen in Figure 7 compared to Figure 6, is not a failure of the optimization. Instead, we think that it reflects experimental errors and/or deficiencies of the model. The latter might include the lack of oxygen mobility, inaccurate representation of the crystallographic plane, incompleteness of microscopic processes, etc. In closing, it is clear that molecular and/or multi-
In summary, we have demonstrated the feasibility of optimizing intrinsic kinetic parameters in molecular models. At present, it is not clear whether the lack of significant differences between the optimized parameters of MF and MC models is generic for surface reactions. Moderate differences between models offer an advantage in optimization, as mean-field theory can be viewed as a lower-level model in the entire hierarchy of surface models. As such, MF models can serve in fast mapping of operation diagrams, sensitivity analysis, and identification of conditions where surface diffusion is not rate-limiting, aspects that can reduce the total number of MC simulations. Extensions of the multistep optimization approach to consider more complex phenomena, such as surface reconstruction, are also possible. Finally, the proposed approach provides a computational route for closing the pressure gap between surface science experiments and actual industrial conditions. Acknowledgment Acknowledgment for partial support of this work is made to the National Science Foundation (CAREER Award CTS-9702615) and to the U.S. Department of Energy (Award DE-FC26-00NT41027). However, any opinions, findings, conclusions, or recommendations expressed herein are those of the authors and do not necessarily reflect the views of the DOE. Literature Cited (1) Auerbach, S. M. Theory and simulation of jump dynamics, diffusion and phase equilibrium in nanopores. Int. Rev. Phys. Chem. 2000, 19, 155. (2) Gomer, R. Diffusion of adsorbates on metal surfaces. Rep. Prog. Phys. 1990, 53, 917. (3) Nieminen, R.; Jansen, A. Monte Carlo simulations of surface reactions. Appl. Catal. A: Gen. 1997, 160, 99. (4) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Oxford Science Publications: Oxford, U.K., 1989. (5) Gilmer, G. Computer models of crystal growth. Science 1980, 208, 355. (6) Keil, F. J.; Krishna, R.; Coppens, M. O. Modeling of diffusion in zeolites. Rev. Chem. Eng. 2000, 16, 71. (7) Aghalayam, P.; Park, Y. K.; Fernandes, N. E.; Papavassiliou, V.; Mhadeshwar, A. B.; Vlachos, D. G. A C1 mechanism for methane oxidation on platinum. J. Catal. 2003, 213, 23. (8) Aghalayam, P.; Park, Y. K.; Vlachos, D. G. Construction and optimization of detailed surface reaction mechanisms. AIChE J. 2000, 46, 2017. (9) Frenklach, M.; Wang, H.; Rabinowitz, M. J. Optimization and analysis of large chemical kinetic mechanisms using the solution mapping methodsCombustion of methane. Prog. Energy Combust. Sci. 1992, 18, 47. (10) Zhdanov, V. P.; Kasemo, B. Kinetic phase transitions in simple reactions on solid surfaces. Surf. Sci. Rep. 1994, 20, 111. (11) Kisliuk, P. The sticking probabilities of gases chemisorbed on the surfaces of solids. J. Phys. Chem. Solids 1957, 3, 95. (12) Gasser, R.; Smith, E. A surface mobility parameter for chemisorption. Chem. Phys. Lett. 1967, 1, 457. (13) Shustorovich, E.; Sellers, H. The UBI-QEP method: A practical theoretical approach to understanding chemistry on transition metal surfaces. Surf. Sci. Rep. 1998, 31, 1.
Ind. Eng. Chem. Res., Vol. 42, No. 6, 2003 1183 (14) Aghalayam, P.; Park, Y. K.; Vlachos, D. G. A detailed surface reaction mechanism for CO oxidation on Pt. Symp. (Int.) Combust. 2000, 28, 1331. (15) Ziff, R. M.; Gulari, E.; Barshad, Y. Kinetic phase transitions in an irreversible surface-reaction model. Phys. Rev. Lett. 1986, 56, 2553. (16) Bortz, A. B.; Kalos, M. H.; Lebowitz, J. L. A new algorithm for Monte Carlo simulations of Ising spin systems. J. Comput. Phys. 1975, 17, 10. (17) Reese, J. S.; Raimondeau, S.; Vlachos, D. G. Monte Carlo algorithms for complex surface reaction mechanisms: Efficiency and accuracy. J. Comput. Phys. 2001, 173, 302. (18) Raimondeau, S.; Vlachos, D. G. The role of adsorbate-layer nonuniformities in catalytic reactor design: Multiscale simulations for CO oxidation on Pt. Comput. Chem. Eng. 2002, 26, 965. (19) Raimondeau, S.; Vlachos, D. G. Recent developments on multiscale, hierarchical modeling of chemical reactors. Chem. Eng. J. 2002, 90, 3. (20) Vlachos, D. G.; Schmidt, L. D.; Aris, R. Kinetics of faceting of crystals in growth, etching, and equilibrium. Phys. Rev. B 1993, 47, 4896. (21) Zhdanov, V.; Kasemo, B. Simulations of the reaction kinetics on nanometer supported catalyst particles. Surf. Sci. Rep. 2000, 39, 29. (22) Liu, J.; Xu, M.; Zaera, F. Determination of the rate-limiting step in the oxidation of CO on Pt(111) surfaces. Catal. Lett. 1996, 37, 9. (23) Hickman, D. A.; Schmidt, L. D. Steps in CH4 oxidation on Pt and Rh surfaces: High-temperature reactor simulations. AIChE J. 1993, 39, 1164. (24) Deutschmann, O.; Schmidt, R.; Behrendt, F.; Warnatz, J. Numerical modeling of catalytic ignition. In Twenty-Sixth Symposium (International) on Combustion; The Combustion Institute: Pittsburgh, 1996; p 1747. (25) Rinnemo, M.; Kulginov, D.; Johansson, S.; Wong, K. L.; Zhdanov, V. P.; Kasemo, B. Catalytic ignition in the CO-O2 reaction on platinum: Experiment and simulations. Surf. Sci. 1997, 376, 297.
(26) Raimondeau, S.; Aghalayam, P.; Vlachos, D. G.; Katsoulakis, M. Bridging the gap of multiple scales: From microscopic, to mesoscopic, to macroscopic models. In Foundations of Molecular Modeling and Simulation; AIChE Symposium Series No. 325; American Institute of Chemical Engineers (AIChE): New York, 2001; Vol. 97, pp 155-158. (27) Shustorovich, E. Bond-order conservation approach to chemisorption and heterogeneous catalysis: Applications and implications. Adv. Catal. 1990, 37, 101. (28) Dumesic, I. A.; Rud, D. F.; Aparicio, L. M.; Rekoske, J. E.; Revino, A. A. The Microkinetics of Heterogeneous Catalysis; American Chemical Society: Washington, DC, 1993. (29) Tomlin, A. S.; Turanyi, T.; Pilling, M. J. Mathematical tools for the construction, investigation, and reduction of combustion mechanisms. Elsevier Sci. J. 1997, 35, 293. (30) Goodman, R. H.; Graff, D. S.; Sander, L. M.; Leroux-Hugon, P.; Clement, E. Trigger waves in a model for catalysis. Phys. Rev. E 1995, 52, 5904. (31) Sales, J. L.; Unac, R. O.; Gargiulo, M. V.; Bustos, V.; Zgrablich, G. Monte Carlo simulation of temperature programmed desorption spectra: A guide through the forest for monomolecular adsorption on a square lattice. Langmuir 1996, 12, 95. (32) Box, G. E. P.; Draper, N. R. Empirical Model-Building and Response Surfaces; Wiley: New York, 1987. (33) Kirkpatrick, S.; Gelatt, C. D., Jr.; Vecchi, M. P. Optimization by simulated annealing. Science 1983, 220, 671. (34) Davis, S. G.; Mhadeshwar, A. B.; Vlachos, D. G.; Wang, H. A new approach to response surface development for detailed gas-phase and surface reaction kinetic model development and optimization. Int. J. Chem. Kinet., manuscript submitted.
Received for review April 3, 2002 Revised manuscript received December 19, 2002 Accepted January 20, 2003 IE0202470
