A Reactive Molecular Dynamics Study of the Thermal Decomposition

Apr 14, 2009 - in sealed-for-life electric motors, in computer hard drives, and ..... neously delete the reactant molecule and insert the internally e...
0 downloads 0 Views 2MB Size
Download PDF
13670

J. Phys. Chem. B 2009, 113, 13670–13677

A Reactive Molecular Dynamics Study of the Thermal Decomposition of Perfluorodimethyl Ether† Bangwu Jiang, Myvizhi Esai Selvan, David J. Keffer,* and Brian J. Edwards Department of Chemical and Biomolecular Engineering, UniVersity of Tennessee, KnoxVille, Tennessee 37996-2200 ReceiVed: December 17, 2008; ReVised Manuscript ReceiVed: March 10, 2009

Classical reactive molecular dynamics (RMD) simulation is used to model the thermal decomposition of perfluorodimethyl ether (CF3OCF3), which is relevant as a simple molecule containing the necessary architectural elements to study the chemical stability of perfluoropolyether lubricants. The RMD algorithm employs nonreactive interaction potentials for the reactants and products. The reactivity is implemented through a coarse-grained simulation algorithm, incorporating elements from both the quantum and macroscopic descriptions of the reaction. The RMD scheme maps the quantum mechanically determined transition state onto a set of geometric triggers. When a configuration matching those triggers is found in the RMD simulation, the reaction instantaneously occurs. A brief, local equilibration process stabilizes the configuration, and the simulation continues. Using two geometric triggers, the RMD simulation can describe quantitatively the temperature dependence of the thermal decomposition of CF3OCF3, when compared to the quantum mechanical standard. 1. Introduction Perfluoropolyethers (PFPEs) form a class of lubricants that are broadly applied in oxygen service, in aircraft instrument bearings, in reactive chemical environments, in vacuum pumps, in sealed-for-life electric motors, in computer hard drives, and as high-temperature greases. They are particularly well-suited for these applications due to their numerous advantageous properties, including chemical inertness, thermal stability, nonflammability, radiation resistance, and shear stability.1-4 While these PFPEs are intrinsically stable, they decompose under certain conditions, for example, heating in the presence of metallic surfaces and Lewis acids.5 At low temperatures, the degradation of certain types of PFPEs has been induced by electron beam exposure.6 This degradation is primarily associated with scission of the ether linkage in the PFPE backbone.6 PFPEs are generally stable at room temperature, with thermal decomposition beginning at higher than 300 °C in an inert atmosphere.4 The thermal decomposition of perfluorodimethyl ether, CF3OCF3, the smallest molecule containing an ether linkage of the type found in a PFPE, was studied by Pacansky et al.7,8 through quantum mechanical (QM) calculations. They showed that the thermal decomposition of CF3OCF3 involves a “cyclic” transition state, rather than a direct C-O bond scission, to produce free radicals. In the present work, a method was developed and applied to study the thermal degradation of CF3OCF3 at high temperature using the classical, nonreactive force field within a reactive molecular dynamics (RMD) algorithm. The term reactive molecular dynamics has previously appeared in the literature.9-11 In this work, our algorithm falls within the broad categorization of RMD but uses a novel scheme for accounting for reaction. †

Part of the “H. Ted Davis Special Section”. * To whom correspondence should be addressed. E-mail: [email protected]. Phone: (865)-974-5322.

We present this RMD algorithm as a coarse-grained model on the basis that we incorporate the chemical reaction, which involves the redistribution of electronssa problem in the quantum mechanical domainsinto a classical RMD algorithm. With regards to the two paradigms of multiscale modeling, multiscale (LM) and closely-coupled multiscale (CCM), we classify this RMD algorithm in the LM approach, in which we map details obtained from the quantum mechanical analysis (performed first) into a RMD algorithm (performed last). Specifically, the coarse-graining in this work refers to the elimination of the quantum mechanical degrees of freedom associated with the change in electron distribution that occurs during chemical reaction. Rather than explicitly solve the Schro¨dinger equation to describe the reaction each time it takes place in the simulation, we perform the quantum mechanical analysis first. We then take information from the quantum mechanical analysis, including the activation energy, the heat of reaction, the reaction rate, and the geometry of the transition state, and incorporate it into a classical RMD algorithm, in which the key properties of the reaction, the activation energy, the heat of reaction, the reaction rate, are reproduced. The validation of a coarse-grained scheme is first and foremost one of internal consistency. Does the coarse-grained model deliver the same results as the more finely resolved model? In this work, we compare the results of the coarsegrained model (classical RMD) and the more finely resolved level (QM calculations) as validation of the reasonableness of the reduction in degrees of freedom. The most finely resolved description of chemical reactions is achieved through direct QM calculation. However, for reactive systems involving macromolecules or systems in which reaction must be coupled to a longer time scale phenomena such as transport, classical force field methods are one of the leading alternatives due to their much lower computational cost relative to QM calculations.12 There are a variety of existing coarse-grained approaches that incorporate chemical reactivity in a simulation without resorting

10.1021/jp811151m CCC: $40.75  2009 American Chemical Society Published on Web 04/14/2009

Thermal Decomposition of Perfluorodimethyl Ether to fully QM calculations. These methods include the surface hopping (SH)13 approach, the empirical valence bond (EVB)14 approach, and adiabatic reactive molecular dynamics (ARMD).15 Alternatively, there is a class of methods called empirical force field methods12 in which reactive potentials can be generated by inserting the detailed information regarding the breaking and forming of bonds into traditional interaction potentials.12 The reactive empirical bond order (REBO) potential16,17 and its revisions18 have been used to explore reactivity in covalent systems, such as silicon and hydrocarbons. The REBO-type potential can allow all possible reactions to occur; however, the potential has not been parametrized for the entire periodic table, including fluorine. In contrast to REBO, this RMD algorithm is implemented through a coarse-grained simulation algorithm rather than through manipulation of the interaction potential. One feature of a RMD algorithm is that one must know a priori the chemical reactions of interest. This has advantages and disadvantages. The obvious disadvantage is that the only reactions that will take place are those that were explicitly addressed in the algorithm. One advantage of the RMD algorithm is that the parametrization can be performed on a much shorter time scale than the parametrization of a generalized reactive potential since the goal of the RMD algorithm is to address only a single reaction or small set of reactions. An ideal RMD algorithm should incorporate elements from both the quantum and macroscopic descriptions of reaction. Stoliarov et al. applied this approach to simulate thermal decomposition reactions in poly(methyl methacrylate) and poly(isobutylene).11,19 The reaction was determined by the changing of bond order. A general mechanism for the thermal decomposition of vinyl polymers was formulated on the basis of the results of these simulations. Schmidt et al. investigated the proton transport in water using the coarse-grained RMD algorithm in a mixed molecular dynamics (MD) and Monte Carlo (MC) scheme, in which the migration of protons was treated as a two step process, the displacement and proton transfer.20 Lill et al. simulated proton transport through classical molecular dynamics simulation coupled with quantum mechanically derived proton hopping,21 which provides an attractive method to study proton transport in a large number of biological systems. Currently, we apply the RMD approach to investigate the thermal decomposition of perfluorodimethyl ether. The occurrence of reaction is based on the geometric change of the relative bonds, which is directly related to the transition-state structure of the reaction as obtained through quantum mechanics calculation. In this RMD algorithm, we use a published, nonreactive interaction potential for the reactants and products. We perform conventional nonreactive molecular dynamics simulations with an added RMD feature. At the end of each time step in the MD simulation, we evaluate the configuration of the reactant and, on the basis of a set of geometric triggers, determine whether a reaction will instantaneously occur. Reaction is followed by a local equilibrium to adjust the position of product molecules to satisfy the known value of the heat of reaction. After the local equilibration, the molecular dynamics simulation proceeds to the next time step. The remainder of this paper is organized as follows. A description of the coarse-grained RMD algorithm and the QM calculations are presented in section 2. Results and discussions from the RMD simulations are presented in section 3. The Conclusions are given in section 4.

J. Phys. Chem. B, Vol. 113, No. 42, 2009 13671 2. Quantum Mechanical Calculations and the Coarse-Grained Reactive Molecular Dynamics Algorithm In this section, we present the conceptual and practical considerations in the formulation of the RMD algorithm. Macroscopically, the general description of a chemical reaction requires the following parameters: the activation energy for the reaction, Ea, the heat of reaction, ∆Hr, and the rate constant prefactor, k0, which appears in the expression for the rate constant, k

( )

k ) k0 exp

-Ea RT

(1)

The activation energy, Ea, generally corresponds to the energy difference between the reactant and transition-state structures. The heat of reaction, ∆Hr, is the enthalpy difference between the isolated reactant and product molecules. Here, we consider the thermal decomposition of CF3OCF3 as one specific example of a chemical reaction. Specifically, our reaction is given by

CF3OCF3 f CF4 + COF2

(2)

The rate of this elementary unimolecular reaction (with units of mol/volume/time) is given by

r ) k[CF3OCF3]

(3)

where the square brackets indicate molar concentration. There are several methods to describe this reaction, including experimental investigation, QM calculations, and MD simulation. We focus on the latter two methods in this work. 2.1. Quantum Mechanical Calculations. From the quantum mechanical point of view, the description of a chemical reaction is a matter of evaluating the structures and energies of the relevant atoms in the reactant, transition-state, and product configurations. The determination of the activation energy and heat of reaction at 0 K is therefore a matter of determining the energy of these three configurations. The reactant and product possess minima in the energy landscape, and the transition state is a saddle point. We applied density functional theory (DFT) to calculate the optimized geometries and harmonic vibrational frequencies using the B3LYP22-25 and Hartree-Fock methods with the 6-31G* basis set26,27 with Gaussian 98,28 all of which has been used earlier and found to be effective to compute the transition-state structure for CF3OCF3-AlF3.29 Transition-state theory30 was applied to estimate the reaction rate constants as functions of temperature. The Hartree-Fock method was used solely for the purpose of comparing the results of the present work with the published results of Pacansky et al.8 Although Pacansky et al. used the Mulliken computer code31 to carry out their ab initio calculations, the geometrical and energetic information that we obtained through Gaussian 9828 is completely consistent with their results. The B3LYP method, combined with the 6-31G* basis set, was applied to obtain vibrational frequencies and geometrical information to calculate the reaction constants via transition-state theory.30 No experimental thermokinetic data are available for the thermal decomposition of CF3OCF3. Therefore, the quantum methods and transition-state theory were used to estimate the rate constants of this reaction using the ideal polyatomic gas model of statistical mechanics.32 The same procedure was applied in previous work29 to estimate the reaction constants of the decomposition of CF3OCF3 in the presence of AlF3 according to the following formula30,33

13672

J. Phys. Chem. B, Vol. 113, No. 42, 2009

( )( ) ( )

k ) κ(T)

p0 RT0

-m

∆G0 kBT exp h RT

Jiang et al.

(4)

where κ(T) is the tunneling correction term, T is the absolute temperature, T0 and p0 are the temperature and pressure references, kB is the Boltzmann constant, m is the change in the number of molecules from the reactants to the transition state, h is Planck’s constant, R is the ideal gas constant, and ∆G0 is the change of the Gibbs free energy from the reactants to the transition state. The Gibbs free-energy change is expressed as

∆G0 ) ∆H0 - T∆S0

(5)

where ∆H0 and ∆S0 are the changes of enthalpy and entropy from the reactants to the transition state, which are estimated through vibrational frequencies obtained with B3LYP/6-31G*. The factor 0.9806 was applied to zero-point/thermal energies.34-36 The Wigner tunneling expression was applied to account for the tunneling contribution to the rate constant involving the transition state

κ(T) ) 1 +

( )

1 hνs 24 kBT

2

(6)

where νs is the imaginary frequency in the transition state. The vibrational mode with the lowest frequency in both the ground state and transition state was treated as internal rotation, calculated by assuming the torsional potential to have the simple form U(φ) ) V[1 - cos(σintφ)]2, using the tables in ref 37, and the remainder were treated harmonically. Eigenmovies of all of the normal modes of the ground state and transition state visually confirmed that the lowest frequency was indeed internal rotation.38,39 From the QM calculations, we therefore obtain the numerical values of Ea, ∆Hr, and k0. From the point of view of the RMD algorithm, the origin of these properties is unimportant; they could just have been equally generated from experiment. However, the QM calculations also provide the transition-state structure, which is used to define a set of geometric triggers that must be satisfied in the RMD algorithm in order for a reaction to take place. For this reason, even if experimental information had been used to provide Ea, ∆Hr, and k0, one would still need to know the ground-state structure of the reactants and the transition-state structure in order to identify meaningful geometric triggers in the coarse-grained model. 2.2. Interaction Potential in RMD. The basic interaction potential used in the simulations was the well-known universal force field (UFF) developed by Rappe et al.40 The potential energy includes contributions from bond stretching, ER, bond angle bending, Eθ, dihedral angle torsion, Eφ, and both intramolecular and intermolecular nonbonded van der Waals interactions, Evdw. Specifically, the bond-stretching energy assumes the form of a Hookean spring

1 ER ) kij(r - rij)2 2

(7)

where kij is the force constant for bond stretching and rij is the equilibrium bond length between atoms i and j. To describe the behavior of the CO bond stretching, the Morse function was applied solely to the CO bond stretch

ER ) D[1 - exp(-R(r - rij))]2

(8)

In this expression, D is the bond dissociation energy (D ) kij/ (2R2)), and R is 2.0 for the CO bond stretch.11 The expressions for bond angle, bond torsion, and nonbonded interactions can

Figure 1. Illustration of triggers used in the RMD algorithm for CF3OCF3.

be found in the original UFF force field40 and previous work.41-43 The CO bond stretch functional form was changed to aid in the local equilibration step of the algorithm. In the simulations, the CO bond does not break without the RMD algorithm. 2.3. Conceptual Implementation and Simulation Details. The basis of any classical molecular dynamics algorithm is solving the equations of motion for the atomic trajectories via numerical integration with respect to time. The RMD algorithm modifies a conventional nonreactive MD algorithm by inserting the possibility of chemical reaction at the end of each MD time step. The RMD component of the simulation can be broken into three steps, (1) identification of molecules which satisfy the triggers for reaction, (2) instantaneous reaction, and (3) local equilibration. These three steps of reaction satisfy the QM description of reaction in a coarse-grained manner. The RMD triggers dictate that the molecule has the proper starting point for the reaction. The instantaneous reaction in RMD coarse-grains out the details of the reaction pathway. The local equilibration dictates that the products have a proper ending point for the reaction. With proper starting and ending points, coupled with parametrization of triggers to fit k0 and Ea and local equilibration to fit ∆Hr, the RMD algorithm is consistent with both quantum mechanical and macroscopic descriptions of reaction. These three steps are described in more detail below. The RMD algorithm should be consistent with both the QM description of the reaction as well as the macroscopic description of reaction. One of the challenges of a RMD algorithm is finding a faithful mapping of the QM description of reaction onto the classical simulation. This mapping is performed by examining the QM ground-state structure of the reactant and comparing it with the QM transition-state structure. In the RMD simulation, the molecule must move toward the QM transition-state structure. Of course, it will never reach the transition state due to the nonreactive potential. Therefore, the geometrical triggers for the reaction must capture the initial movement of the molecule toward the transition state rather than the transition state itself. At the same time, guided by the macroscopic model, we wish to include as little information as possible about the transition state in order to achieve an efficient, coarse-grained model which captures the essential, defining features of the reaction. The geometrical triggers used in RMD are illustrated in Figure 1. Initially, we investigated the reaction under the trigger scheme in which a single geometrical trigger, Rmax, was applied to monitor the CO bond length. Whenever a CO bond length was longer than Rmax, the reaction instantly took place. The second scheme was to monitor the variation of bond lengths for both CO bonds. In addition to Rmax, another geometrical trigger, Rmin, was used to monitor simultaneously another bond length. In a CF3OCF3 molecule, when one CO bond length is larger than Rmax and the other CO bond length is smaller than Rmin, the

Thermal Decomposition of Perfluorodimethyl Ether

J. Phys. Chem. B, Vol. 113, No. 42, 2009 13673

Figure 2. The optimized geometries for CF3OCF3 (right) and its transition-state structure (left) obtained using B3LYP with the 6-31G* basis set.

CF3OCF3 molecule is instantly decomposed into CF4 and COF2. We investigated different sets of numerical values for Rmax and Rmin in order to fit k0 and Ea obtained from the QM calculations. While there are other possible geometric triggers, like the CF bond distance, there was no need to define them since we were able to capture the reaction kinetics with the above two triggers. In the reaction step of this RMD algorithm, we instantaneously delete the reactant molecule and insert the internally equilibrated structures of the product molecules with the same centers of mass as the corresponding components of the reactant molecule. In all likelihood (especially in dense systems), this molecule exchange results in an increase in the potential energy due to intermolecular overlap of atoms. Therefore, we require a local relaxation with the intention of disturbing the trajectories being generated in the RMD simulation as little as possible and simultaneously satisfying the heat of reaction, ∆Hr. The local relaxation, which matches the change in the total energy of the system during the reaction to the ∆Hr, can be achieved by a number of means. The increase in energy can be either distributed between the kinetic and potential energies of the product molecules or involve the equilibration of the neighboring atoms. In the current scenario, the product atoms have the same velocity as that of the reactant atoms. In other words, the instantaneous reaction has caused a modification only in the potential energy, while the kinetic energy of the system has been maintained the same. We then moved only the two product molecules directly related to this specific reaction to reduce the potential energy such that the difference between the energy of the system before and after the reaction was within a specified tolerance of the target, ∆Hr. The equipartition of energy between the kinetic and potential energy is automatically achieved during the subsequent MD steps. It is acknowledged that reactions with larger ∆Hr may require a local equilibration scheme involving a change in both kinetic and potential energies, but for this example, it sufficed to simply adjust the potential energy. The three steps in the RMD algorithm were incorporated in an existing in-house-built classical MD code where the equations of motion are integrated using the two time scale reversible reference system propagator algorithm44 (r-RESPA) and where periodic boundary conditions are applied. The time step for the intermolecular interaction was 2.0 and 0.2 fs for intramolecular interaction. We used 512 (or 216 in the single-trigger case) CF3OCF3 molecules at the density of 0.051 g/mL at 2500-3500 K, which is a dense gas system. We simulated at high temperatures in order to observe a statistically significant number of reactions within the nanosecond duration of the NVT simulations. The temperature was controlled by a Nose´-Hoover thermostat.45,46 The nonbonded interactions were cut off at a radii of 12.5 Å, and the long-range interactions were accounted

TABLE 1: Optimized Bond Distance (Å) and Angle (degree) of Compounds, Calculated Using B3LYP/6-31G* parameters

ground state

transition state

R(C1-O9) R(C2-O9) R(C1-F6) R(C1-F7) R(C1-F8) R(C2-F3) R(C2-F4) R(C2-F5) R(F8 · · · C2) ∠O9C1F6 ∠O9C1F7 ∠O9C1F8 ∠O9C2F3 ∠O9C2F4 ∠O9C2F5 ∠C1O9C2 ∠F6C1F8 ∠F6C1F7 ∠F7C1F8 ∠F3C2F4 ∠F3C2F5 ∠F5C2F4 ∠C1F8C2 ξ-F8C1O9C2 ξ-F5C2O9C1 ξ-F7C1O9C2

1.402 1.393 1.292 1.292 1.301 1.294 1.297 1.297 2.571 108.3 108.3 110.1 106.0 110.4 110.3 121.3 109.7 110.6 108.3 110.1 110.1 110.0 69.2 0.0 -60.9 120.0

1.250 2.171 1.346 1.346 1.556 1.301 1.282 1.288 1.989 119.6 119.6 103.4 145.8 87.5 82.3 97.3 102.2 106.8 102.1 112.5 109.2 118.2 95.4 20.7 133.8 113.8

for. The simulations were run long enough (50-60 ns) to observe the decomposition of a significant fraction of reactant molecules. 3. Results and Discussion 3.1. Reaction Constants through Quantum Calculations. The reaction mechanics of the thermal decomposition of CF3OCF3 have been elaborated by Pacansky et al.8 with different levels of theory. The results from B3LYP/6-31G* are similar to those in ref 8. Here, we provide the key structural parameters relative to this work. The optimized geometries for CF3OCF3 and its transition structure computed through the B3LYP/6-31G* level of theory are illustrated in Figure 2, which demonstrates that in the transition state, one CO bond length increases while the other becomes shorter, relative to the equilibrium CO bond length. The main parameters of the two structures are listed in Table 1. Figure 3 demonstrates that the energy variation along the reaction coordinate validates the transition-state structure; the energy gradually decreases on both sides of the identified transition-state structure. The reaction heat, entropy, and Gibbs free energy of the rate-limiting step, the formation of the transition-state structure, were calculated using statistical mechanics32 and transition-state theory30 at temperatures between

13674

J. Phys. Chem. B, Vol. 113, No. 42, 2009

Jiang et al.

Figure 3. The energy difference along the reaction coordinate computed using B3LYP/6-31G*. All of the data are relative to the energy of the transition state (TS).

Figure 5. Snapshots from RMD simulation (Rmax ) 2.1 Å, T ) 2473 K); blue represents the reactant (CF3OCF3) and red the products (CF2O and CF4).

Figure 4. Arrhenius plot of thermal decomposition rates of CF3OCF3 calculated using B3LYP/6-31G*.

1000 and 3500 K. The corresponding results are shown in Figure 4. The solid line is obtained by fitting the data calculated from quantum mechanics using eq 1. The corresponding values of the parameters are Ea ) 345 kJ/mol and k0 ) 2.57 × 1014 s-1. The heat of reaction, ∆Hr, is 14.1 kJ/mol. 3.2. RMD with a Single Trigger. In the RMD algorithm, one of the most important points is to determine what kinds of triggers are needed to describe the reaction at the atomistic level. For thermal decomposition of CF3OCF3, our initial attempt was to use as few parameters as possible to describe the reaction. Thus, we attempted to model the reaction with a single trigger, Rmax, as mentioned above, which is the requirement for the CO bond length and corresponds to the distance between C2 and O9 in the transition-state structure shown in Figure 2. During the RMD simulation, the number of reaction events is recorded as a function of simulation time at the different temperatures for several values of Rmax. Snapshots at the different reaction times are shown in Figure 5 for Rmax ) 2.1 Å. First- and second-order macroscopic reaction rate equations were applied to describe the change of reactant molecules along the reaction course. The first-order reaction equation is N0 - N ) N0 exp(-kτ), derived from eq 3, where N0 is the initial number of CF3OCF3, N is the number of CF3OCF3, which changes with the simulation time, k is the reaction rate constant, and τ is the reaction time. The first-

Figure 6. The change in the number of reactant molecules (N) during the RMD simulation (Rmax ) 2.1 Å, T ) 2473 K) as a function of time. The solid curve represents the fit to a first-order reaction equation.

order reaction equation consistently describes the reaction better than the second-order equation for the different cases; therefore, we used the first-order reaction equation to fit the data to obtain the reaction rate constant. Moreover, for a unimolecular reaction, the first-order reaction rate is the obvious choice. We investigated the possibility of the secondorder reaction to allow for the intermolecular collisions to possibly provide sufficient energy to force the CF3OCF3 into a configuration satisfying the trigger. Figure 6 displays the change in the number of reactant molecules during the RMD for the case Rmax ) 2.1 Å at the temperature of 2473 K. The solid curve is the fit calculated using the first-order reaction equation. The reaction rate constant is obtained from the fitting. For each value of the trigger, Rmax, we investigated the reaction at different temperatures. Then, the reaction rates at these temperatures were fitted to eq 1 to estimate the prefactor and activation energy. The resulting Arrhenius plot with the trigger Rmax ) 2.1 Å is shown in Figure 7. One can see that the simulation yields reaction rates that are well described by eq 1.

Thermal Decomposition of Perfluorodimethyl Ether

Figure 7. Arrhenius plot of thermal decomposition rates of CF3OCF3 from the RMD simulations for the single-trigger case; Rmax ) 2.1 Å.

J. Phys. Chem. B, Vol. 113, No. 42, 2009 13675

Figure 9. Comparison of the reaction rate constants obtained from the quantum calculation and RMD simulation with different sets of Rmax and Rmin as a function of temperature. In the legend, the first number is Rmax and the second is Rmin.

TABLE 2: The Prefactor, k0, and Activation Energy, Ea, for the Different Triggers, Rmax, Obtained through RMD Simulations with a Single Trigger Rmax/Å

Ea/kJ mol-1

k0/1012 s-1

1.8 1.9 2.0 2.1 2.3

114 169 226 286 376

1.02 2.28 4.87 10.86 24.69

The prefactor and activation energy for the different triggers, Rmax, are listed in Table 2 and are plotted in Figure 8. This figure shows that both the activation energy and the prefactor increase with the increasing trigger size. Clearly, the activation energy, Ea, increases linearly with Rmax, but the prefactor, k0, increases nonlinearly. Note that the left ordinate (for Ea) is linear, and the right one (for k0) is a log scale. In terms of the activation energy, for the simulation to match the QM result of 345 kJ/ mol, Rmax should be around 2.23 Å, which is reasonably close to the value of the transition-state structure, 2.17 Å. However, with this trigger value, the corresponding k0 from RMD is 1.8 × 1013 s-1, which is much lower than the value from the QM calculations of 2.57 × 1014 s-1. On the basis of this evidence, we conclude that it is not possible to model accurately the thermal decomposition of CF3OCF3 with a single trigger. Therefore, the next logical step of the analysis was to introduce a second trigger.

Figure 8. The prefactor, k0, and activation energy, Ea, obtained from RMD simulation for the different triggers, Rmax. The dashed line shows the activation energy from the QM calculations.

Figure 10. Arrhenius plot of thermal decomposition rates of CF3OCF3 from the two-trigger reactive molecular simulation with Rmax ) 1.65 Å and Rmin ) 1.19 Å.

3.3. RMD with Two Triggers. Since the single-trigger RMD algorithm was incapable of simultaneously generating the target values of Ea and k0, we next used a two-trigger RMD algorithm. While there are several geometrical parameters important to the formation of the transition-state structure in the thermal decomposition of CF3OCF3 since the reaction involves the breakage of one CO bond and the transfer of another CO bond from a single bond to double bond, it was natural to choose geometrical parameters related to the change of CO bond lengths, Rmin and Rmax. Their definitions were given above. The determination of the values of Rmin and Rmax that generate the target values of Ea and k0 is a nonlinear twovariable optimization problem. Each iteration in this optimization procedure requires a set of RMD simulations across a range of temperatures in order to extract Ea and k0 from the RMD simulations. The geometrical information from the QM calculations can be used to guide this search. In Figure 9, we plot the reaction rate as a function of temperature for various combinations of Rmax and Rmin. The reaction rate from the QM calculations is also shown. From the plot, we observe that at the same Rmax, a decrease in Rmin will decrease the reaction rate constants, while at the same Rmin, an increase of Rmax will decrease the reaction rate constants. This is not surprising since choosing triggers farther from equilibrium

13676

J. Phys. Chem. B, Vol. 113, No. 42, 2009

simply reduces the probability of finding a molecule with that configuration. The combination of Rmin ) 1.19 Å and Rmax ) 1.65 Å gave the best match between the RMD simulation results and QM calculations. The fitting of eq 1 to this case is shown in Figure 10. The fitting specifies an activation energy of 341 kJ/mol and a prefactor of 2.26 × 1014 s-1, which match the QM values of 345 kJ/mol and 2.57 × 1014 s-1 to within 1.2 and 12%, respectively. At this point, we conclude that a two-trigger RMD algorithm is sufficient to reproduce the reaction rate and its temperature dependence for the thermal decomposition of CF3OCF3. Finally, we would like to point out that the numerical values of Rmin and Rmax do not correspond precisely to any dimension in the QM transition state because we are performing a coarse-grained mapping from the QM description to this model with a classical and nonreactive interaction potential. Therefore, we expect the numerical values of the trigger to correspond only semiquantitatively. The values that match quantitatively to the output of the coarse-grained simulation are Ea, k0, and ∆Hr. Since the parametrization is an optimization problem, it is also possible for another choice of the triggers to satisfy the same Ea and k0. 4. Conclusions A conventional, well-tested, nonreactive interaction potential has been integrated into a RMD algorithm to model the thermal decomposition of CF3OCF3. The reactivity of the RMD algorithm is implemented via a three-step process involving satisfaction of reaction triggers, instantaneous reaction, and local equilibration. The RMD algorithm is firmly connected to quantum mechanical calculations through the definition of geometrical triggers based on the transitionstate structures obtained from QM calculations. The RMD algorithm is firmly connected to macroscopic descriptions of the chemical reaction through the direct evaluation of activation energy, rate constant prefactor, and heat of reaction from the simulation, which can be used in macroscopic reaction rate expressions. For the specific case of the thermal decomposition of CF3OCF3, we found that a single trigger was insufficient to satisfy simultaneously both the activation energy and the reaction rate constant. However, a two-trigger scheme, in which decomposition takes place when one CO bond is greater than a certain distance while the other CO bond is shorter than a certain distance, did provide a satisfactorily quantitative description of the decomposition reaction. This RMD algorithm has general applicability to a variety of reactions. The algorithm also has practical advantages. The RMD scheme does not add significant computational time to a conventional nonreactive MD simulation. Also, there is a relatively small investment of time in the parametrization of the triggers, compared to the development of a generalized reactive potential, because the number of parameters is much fewer and because we have a good initial guess for the trigger values that can be obtained based on the QM transition state. Moreover, the mapping of the geometric triggers gives the RMD algorithm a built-in adaptability to its local environment, such as temperature, which has been demonstrated in the results. The RMD algorithm also gives rise to the possibility of the reaction rate responding to changes in other environmental factors like density or confinement. Finally, it allows us to simulate systems in which transport properties are coupled to the reaction kinetics by bridging the disparate time scales.

Jiang et al. Acknowledgment. This work has been supported by Air Force Office of Scientific Research through Contract # FA 9550-05-1-0342. The authors wish to acknowledge the resources of the Center for Computational Sciences at Oak Ridge National Laboratory, which is supported by the Office of Science of the DOE under Contract DE-AC0500OR22725. References and Notes (1) Helmick, L. S.; Gschwender, L. J.; Sharma, S. K.; Snyder, C. E.; Liang, J. C.; Fultz, G. W. Tribol. Trans. 1997, 40, 393. (2) Jones, W. R.; Bierschenk, T. R.; Juhlke, T. J.; Kawa, H.; Lagow, R. J. Ind. Eng. Chem. Res. 1988, 27, 1497. (3) Snyder, C. E.; Gschwender, L. J.; Scott, O. L. Tribol. Trans. 1995, 38, 733. (4) Bell, G. A.; Howell, J.; Del Pesco, T. W. Synthetic Lubricants and High-Performance Functional Fluids; Marcel Dekker, Inc.: New York, 1999. (5) Kasai, P. H.; Wheeler, P. Appl. Surf. Sci. 1991, 52, 91. (6) Pacansky, J.; Waltman, R. J.; Jebens, D.; Heery, O. J. Fluorine Chem. 1997, 82, 85. (7) Pacansky, J.; Waltman, R. J. J. Fluorine Chem. 1997, 83, 41. (8) Pacansky, J.; Waltman, R. J. J. Fluorine Chem. 1997, 82, 79. (9) Nyden, M. R.; Forney, G. P.; Brown, J. E. Macromolecules 1992, 25, 1658. (10) Nyden, M. R.; Coley, T. R.; Mumby, S. Polym. Eng. Sci. 1997, 37, 1496. (11) Stoliarov, S. I.; Westmoreland, P. R.; Nyden, M. R.; Forney, G. P. Polymer 2003, 44, 883. (12) Santiso, E. E.; Gubbins, K. E. Mol. Simul. 2004, 30, 699. (13) Tully, J. C.; Preston, R. K. J. Chem. Phys. 1971, 55, 562. (14) Warshel, A.; Weiss, R. M. J. Am. Chem. Soc. 1980, 102, 6218. (15) Danielsson, J.; Meuwly, M. J. Chem. Theory Comput. 2008, 4, 1083. (16) Tersoff, J. Phys. ReV. Lett. 1986, 56, 632. (17) Tersoff, J. Phys. ReV. Lett. 1988, 61, 2879. (18) Stuart, S. J.; Tutein, A. B.; Harrison, J. A. J. Chem. Phys. 2000, 112, 6472. (19) Stoliarov, S. I.; Lyon, R. E.; Nyden, M. R. Polymer 2004, 45, 8613. (20) Schmidt, R. G.; Brickmann, J. Ber. Bunsen-Ges. Phys. Chem. Chem. Phys. 1997, 101, 1816. (21) Lill, M. A.; Helms, V. J. Chem. Phys. 2001, 115, 7993. (22) Lee, C. T.; Yang, W. T.; Parr, R. G. Phys. ReV. B 1988, 37, 785. (23) Becke, A. D. J. Chem. Phys. 1996, 104, 1040. (24) Miehlich, B.; Savin, A.; Stoll, H.; Preuss, H. Chem. Phys. Lett. 1989, 157, 200. (25) Parr, R. G.; Yang, W. Density-functional Theory of Atoms and Molecules; Oxford University Press: New York, 1989. (26) Petersson, G. A.; Allaham, M. A. J. Chem. Phys. 1991, 94, 6081. (27) Petersson, G. A.; Bennett, A.; Tensfeldt, T. G.; Allaham, M. A.; Shirley, W. A.; Mantzaris, J. J. Chem. Phys. 1988, 89, 2193. (28) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scusera, G. E.; Robb, M. A.; Cheeseman, J. R.; Zakrzewski, V. G.; Montgomery, J. A.; Stratmann, R. E., Jr.; Burant, J. C.; Dapprich, S.; Millam, J. M.; Daniels, A. D.; Kudin, K. N.; Strain, M. C.; Farkas, O.; Tomasi, J.; Barone, V.; Cossi, M.; Cammi, R.; Mennucci, B.; Pomelli, C.; Adamo, C.; Clifford, S.; Ochterski, J.; Petersson, G. A.; Ayala, P. Y.; Cui, Q.; Morokuma, K.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Cioslowski, J.; Ortiz, J. V.; Baboul, A. G.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Gomperts, R.; Martin, R. L.; Fox, D. J.; Keith, T.; Al-Laham, M. A.; Peng, C. Y.; Nanayakkara, A.; Gonzalez, C.; Challacombe, M.; Gill, P. M. W.; Johnson, B.; Chen, W.; Wong, M. W.; Andres, J. L.; Gonzalez, C.; Head-Gordon, M.; Replogle, E. S.; Pople, J. A. Gaussian 98; Gaussian, Inc.: Pittsburgh, PA, 1998. (29) Jiang, B. W.; Keffer, D. J.; Edwards, B. J. J. Phys. Chem. A 2008, 112, 2604. (30) Glasstone, S.; Laidler, K. J.; Eyring, H. The Theory of Rate Process; McGraw-Hill: New York, 1941. (31) Rice, J. E.; Horn, H.; Lengsfield, B. H.; McLean, A. D.; Carter, J. T.; Replogle, E. S.; Barnes, L. A.; Maluendes, S. A.; Lie, G. C.; Gutowski, M.; Rudge, W. E.; Sauer, S. P. A.; Lindh, R.; Andersson, K.; Chevalier, T. S.; Widmark, P. O.; Bouzida, D.; Pacansky, G.; Singh, K.; Gillan, C. J.; Carnevali, P.; Swope, W. C.; Liu, B. Mulliken: A Computational Quantum Chemistry Program; Almaden Research Center, IBM Research Division: San Jose, CA, 1995. (32) McQuarrie, D. A. Statistical Mechanics; University Science Books: Sausalito, CA, 2000. (33) Leininger, J. P.; Minot, C.; Lorant, F.; Behar, F. J. Phys. Chem. A 2007, 111, 3082. (34) Foresman, J. B.; Frisch, A. J. Exploring Chemistry with Electronic Structure Methods; Gaussian, Inc: Pittsburgh, PA, 1993. (35) Scott, A. P.; Radom, L. J. Phys. Chem. 1996, 100, 16502.

Thermal Decomposition of Perfluorodimethyl Ether (36) Bauschlicher, C. W.; Partridge, H. J. Chem. Phys. 1995, 103, 1788. (37) Pitzer, K. S.; Gwinn, W. D. J. Chem. Phys. 1942, 10, 428. (38) Kassaee, M. H.; Keffer, D. J.; Steele, W. V. J. Mol. Struct.: THEOCHEM 2007, 802, 23. (39) Kassaee, M. H.; Keffer, D. J.; Steele, W. V. J. Chem. Eng. Data 2007, 52, 1843. (40) Rappe, A. K.; Casewit, C. J.; Colwell, K. S.; Goddard, W. A.; Skiff, W. M. J. Am. Chem. Soc. 1992, 114, 10024. (41) Jiang, B. W.; Keffer, D. J.; Edwards, B. J. J. Fluorine Chem. 2006, 127, 787.

J. Phys. Chem. B, Vol. 113, No. 42, 2009 13677 (42) Jiang, B.; Crawford, N. J.; Keffer, D. J.; Edwards, B. J.; Adcock, J. L. Mol. Simul. 2007, 33, 871. (43) Jiang, B.; Kim, J. M.; Keffer, D. J.; Edwards, B. J. Mol. Simul. 2008, 34, 231. (44) Tuckerman, M.; Berne, B. J.; Martyna, G. J. J. Chem. Phys. 1992, 97, 1990. (45) Hoover, W. G. Phys. ReV. A 1985, 31, 1695. (46) Nose, S. Mol. Phys. 1984, 52, 255.

JP811151M