J. Phys. Chem. B 1999, 103, 3691-3698
3691
A QM/MM Direct Dynamics Trajectory Investigation of Trimethylene Decomposition in an Argon Bath Kim Bolton† and William L. Hase* Department of Chemistry, Wayne State UniVersity, Detroit, Michigan 48202
Charles Doubleday, Jr.* Department of Chemistry, Columbia UniVersity, New York, New York 10027 ReceiVed: July 13, 1998; In Final Form: February 22, 1999
The decomposition of trimethylene embedded in an argon environment is studied using classical trajectory simulations. The trimethylene intramolecular forces are calculated directly from semiempirical electronic structure theory at each trajectory step, and the intermolecular forces are determined from Lennard-Jones 12-6 potential energy functions. When the argon is in a low to moderate density fluid phase, it does not affect the rate of trimethylene decomposition or product branching ratios, but at high densities, the H-transfer reaction leading to propene is favored over cyclization.
Introduction The intramolecular and unimolecular dynamics of trimethylene and its role in the isomerization of isotopically substituted cyclopropanes have been the subject of extensive experimental1-10 and theoretical11-28 investigation. Recently, classical trajectory simulations based on two different potential energy surfaces (PESs) have provided insight into the microscopic details of these reactions.24-28 One set of studies24-27 employed the direct dynamics technique,29 where the energy and forces are obtained directly from electronic structure theory calculations at each trajectory step. These studies were based on a semiempirical PES, with AM1 parameters modified to reproduce CASSCF stationary point properties and experimental activation energies.24 The resulting PES was called AM1-SRP (AM1 with specific reaction parameters30,31). The other study employed an analytic PES where three “critical” degrees of freedom, the central trimethylene CCC angle and the two terminal methylene torsions, were described by functions fit to CASPT2N/6-31G*// GVB-PP(1)/6-31G* energies, and the remaining degrees of freedom were described by molecular mechanics functions.28 Simulations of the thermal isomerization of cyclopropane1,2-d2 based on these PESs yielded similar ratios of the rate constants for double methylene rotation, k12, to single methylene rotation, k1 (k12/k1 ) 2.3-3.5 on the AM1-SRP PES26,27 and 4.7 on the analytic PES28). These studies also gave clear evidence that the dynamics of trimethylene is nonstatistical at 695 K.24-28 The simulated k12/k1 ratios lie between the experimentally deduced value of 1.0 obtained by Baldwin and collaborators6 and 5-42 obtained by Berson et al.3,4 The similarity of the k12/ k1 values obtained from the two PESs and their approximate agreement with experiment supports the use of these PESs to gain insight into trimethylene dynamics at 695 K. Additional studies, based on the AM1-SRP PES, examined trimethylene decomposition at fixed energies.24,25 For energies less than 77 kcal mol-1 (including zero point energy), the
trimethylene dynamics is nonstatistical, in agreement with the data obtained under thermal conditions. At higher energiess between 107 and 146.4 kcal mol-1sthe dynamics is statistical and the decomposition single exponential, in agreement with the molecular beam data of Zewail and co-workers.10 At 146.4 kcal mol-1 the simulated unimolecular rate constant of 11 ps-1 is similar to the 8.3 ( 1.4 ps-1 measured experimentally. This indicates that the AM1-SRP PES also yields valid trimethylene dynamics and decomposition at elevated as well as thermal energies. In the above simulations the dynamics and decomposition of isolated trimethylene molecules were studied, thereby assuming that collisional effects in the molecular beam and thermal experiments are insignificant. Due to the short trimethylene lifetime (tenths of picoseconds10,24,27), this is expected to be a good assumption.32 However, most organic reactions that are performed in the gas, liquid, and matrix environments are expected to be sensitive to collisional and frictional33,34 effects. An example is the trans f cis isomerization of stilbene, where increasing the collision frequency (buffer gas pressure) leads to an initial increase in the isomerization rate, followed by a decrease at higher collision frequencies (solvent densities).34,35 In the present study, collisional and frictional effects on the intra- and unimolecular dynamics of trimethylene are examined. This molecule was chosen since previous studies have provided a detailed microscopic understanding of the unimolecular dynamics of isolated trimethylene, i.e., in a collision free environment. A mixed QM/MM approach36-39 is used to simulate the decomposition of trimethylene embedded in an Ar environment. Trimethylene intramolecular forces are obtained directly from the AM1-SRP PES, and intermolecular forces are obtained from pairwise Lennard Jones (LJ) 12-6 energy functions. The sensitivity of trimethylene decomposition to collisional and frictional effects is examined over a range of Ar temperatures and densities. QM/MM Potential Energy Surface
†
Present address: Go¨teborg University, Department of Chemistry, Physical Chemistry, S-412 96, Go¨teborg, Sweden.
The QM/MM Hamiltonian has the form
10.1021/jp982988d CCC: $18.00 © 1999 American Chemical Society Published on Web 04/16/1999
3692 J. Phys. Chem. B, Vol. 103, No. 18, 1999
H ) HQM + HMM + HQM/MM
Bolton et al.
(1)
where HQM is the quantum mechanical Hamiltonian, HMM is the molecular mechanics Hamiltonian (which could, in general, include any type of analytic potential energy function), and HQM/MM describes the coupling between the QM and MM subregions. In the laboratory fixed frame the kinetic energy T(p) is a function of the atomic momenta only,40 and eq 1 can be expressed as
H(q, p) ) VQM(q) + VMM(q) + VQM/MM(q) + T(p) (2) For the Ar-trimethylene system, VQM(q) is the AM1-SRP PES describing the trimethylene intramolecular energy, VMM(q) is the Ar-Ar intermolecular potential, and VQM/MM(q) describes the Ar-trimethylene interactions. Trimethylene Intramolecular PES. Features of the AM1SRP PES, such as stationary point structures, energies, and harmonic vibrational frequencies, have been detailed elsewhere.24 Using the half-electron method41 available in MOPAC 7.0,42 configuration interaction was included in the AM143 wavefunction and the resulting PES was fine-tuned to reproduce CASSCF ab initio stationary point features23 and experimental activation energies.1,2,5 This was done by scaling the two-center one-electron (resonance) integrals associated with the terminal C-C interactions (to fit cyclization transition state properties) and those associated with C-H interactions between the central hydrogens and the terminal carbons (to fit features of the H-transfer barrier). Inverted Morse-type potential energy terms were appended to the AM1 Hamiltonian to obtain the correct cyclopropane and propene exothermicities. The resulting surface, which bears a close resemblance to the ab initio surface, is illustrated in Figure 1. The curves that link the saddle point structures merely serve to show the connectivity. The disrotatory saddle point, which has two imaginary frequencies and is 2.2 kcal mol-1 above the conrotatory saddle point, and the saddle point for single terminal isomerization, which is 1.4 kcal mol-1 above the conrotatory saddle point, are not shown in Figure 1. The disrotatory saddle point structure is Cs and has a CCC angle of 109.6°. At the saddle point for single terminal rotation, the methylene terminals are staggered by 90° and the CCC angle is 111.9°. Whereas this saddle point connects the cis and trans isomers of the (deuterated) cyclopropane on the ab initio surface,23 it connects the two enantiomeric structures of trimethylene on the AM1-SRP PES. Two stationary points particularly important for the trimethylene dynamics reported here are the trimethylene minimum and the H-atom transfer transition state. The AM1-SRP harmonic vibrational frequencies for these stationary points, listed in ref 24, are in good agreement with ab initio values.23,44 The transition states for geometric isomerization have frequencies similar to those for trimethylene. Intermolecular PES. The intermolecular potential between the Ar atoms, VMM(q), and between the Ar and trimethylene, VQM/MM(q), is described by a sum of pairwise LJ 12-6 functions
[(σr ) - (σr ) ]
V ) 4
12
6
(3)
for Ar-Ar, Ar-C, and Ar-H interactions. The parameters are taken from the PES developed by Lim and Gilbert,45 where collisional energy transfer between inert gases and large hydrocarbons was studied, and are Ar-Ar ) 0.23 kcal mol-1, σAr-Ar ) 3.47 Å, Ar-C ) 0.10 kcal mol-1, σAr-C ) 3.39 Å, Ar-H ) 0.05 kcal mol-1, and σAr-H ) 3.23 Å. The cutoff for
Figure 1. . Minimum potential energy along the cyclopropane T trimethylene T propene reaction path.
the intermolecular interactions is 7.5 Å. Using larger cutoffs did not change the simulation results. For this investigation, highly accurate intermolecular forces are not necessary since qualitative aspects of the Ar environment are investigated. Interactions between atomic partial charges, which partially account for polarization effects and are included in the study of Lim and Gilbert, are not included in this work. By excluding polarization effects, which would also have to be incorporated into the AM1-SRP Hamiltonian,36 the focus is on collisional and frictional effects only, not on electronic changes of the trimethylene PES due to the presence of the environment. This is referred to as mechanical embedding of trimethylene in an Ar bath.38 Trajectory Methods Trajectory Initialization. The simulations employed periodic boundary conditions,46 with the trimethylene initialized at the center of a 20 Å × 20 Å × 20 Å box (or 16 Å × 16 Å × 16 Å at higher densities) and surrounded by a thermalized Ar bath. The trimethylene had an initial energy (including the 47.3 kcal mol-1 zero point energy) of 54.6 or 146.4 kcal mol-1 and was rotationally cold. Under collision free conditions, trimethylene decomposition is nonstatistical at the lower energy and statistical at the higher energy.23 The efficient microcanonical47 and quasiclassical normal mode48 sampling procedures, which have been detailed elsewhere,23 were used to select the initial trimethylene Cartesian coordinates and momenta. By using two sampling schemes, the effects of the initial conditions on the ensuing intra- and unimolecular dynamics were ascertained. The efficient microcanonical sampling (EMS) scheme, which takes anharmonic effects in account, prepares microcanonical ensembles of the total trimethylene classical phase space. The dividing surface between trimethylene and cyclopropane was defined in terms of the trimethylene CCC angle. The transition states (TS) are at CCC ) 94° for an energy of 54.6 kcal mol-1 and CCC ) 89° for 146.4 kcal mol-1. EMS sampling over a range of CCC angles showed23 that these are the classical anharmonic variational49-51 transition states for these energies. The dividing surface between trimethylene and propene is assumed to be at the H-transfer barrier (see Figure 1) and is thus energy independent. For quasiclassical normal mode sampling, the initial conditions are selected from a statistical sampling of quasiclassical, harmonic vibrational states at the H-transfer barrier.23 Each state
QM/MM Investigation of Trimethylene with an energy less than the total trajectory energy has an equal probability of being selected. The difference in the total energy and the energy of the selected state is placed in the reaction coordinate, giving it an energy distribution in accord with RRKM theory.52,53 The initial reaction coordinate momentum is directed toward trimethylene. Random phases are selected for the normal mode coordinates and momenta before the Cartesian coordinates are determined from the standard linear transformation.54 Since this transformation is not exact at finite energies, the Cartesian coordinates and momenta are scaled to obtain the desired trimethylene energy. This scaling is not severe, with the average difference between the unscaled and scaled energies being ≈10% of the desired energy. The Ar bath was initialized using a thermal Monte Carlo (MC) sampling procedure.46 To initiate the sampling, the Ar atoms were placed on a grid surrounding the trimethylene, such that they were uniformly distributed in the periodic box. By fixing the trimethylene coordinates and momenta during the Ar thermalization, the trimethylene energy distribution, chosen by EMS or quasiclassical sampling, was not affected. Although the thermal conditions obtained from the MC sampling do not depend on the initial Ar positions, the number of MC steps required to achieve thermalization is sensitive to their initial location. For moderate Ar densities, 6 × 10-3 e F e2.5 × 10-2 atoms Å-3, thermalization was obtained in 105 MC steps (for all temperatures 100 K e T e 1000 K). Thermal equilibrium conditions were identified by a constant average potential energy.46 At higher densities, F ) 6-7 × 10-2 atoms Å-3, about 5 times as many MC steps were required for thermalization. Due to the computational expense (of MC and MD), the box length was decreased from 20 to 16 Å for these high densities. All Ar atoms were moved at each MC step, and the step size was chosen so that ≈50-60% of trial configurations were accepted in the Markov chain. At the lowest density (F ) 6.25 × 10-3 atoms Å-3) there were 50 Ar atoms in the periodic box, and at the highest density (F ) 7.32 × 10-2 atoms Å-3) there were 300 atoms. The initial conditions for the first MD trajectory in each ensemble were obtained after the above MC “meltdown” period. The initial conditions of subsequent trajectories were obtained after 5000 MC steps, where the sampling was started from the initial Ar positions of the previous trajectory. Using these positions, instead of reinitializing the Ar atoms on their grid positions, allowed for thermalization after the smaller number of MC steps, thereby improving the simulation efficiency. A total of 5000 steps were required to eliminate correlation between the initial conditions of successive trajectories and to ensure rethermalization of the Ar bath surrounding new trimethylene geometries (chosen for each trajectory from EMS or quasiclassical sampling). To summarize, there are two steps in the sampling of initial conditions. First, 54.6 or 146.0 kcal mol-1 of vibrational energy is added to trimethylene with EMS or quasiclassical normal mode sampling. Trimethylene is then placed in the center of a box, with periodic boundary conditions, and surrounded by an argon bath with an equilibrium temperature and density. Initially, trimethylene is in a nonequilibrium state with respect to the bath, since its coordinates and momenta are held fixed when the bath is equilibrated. Trajectory Integration. Cartesian coordinates and momenta were integrated by solving Hamilton’s equations with a combined fourth order Runge-Kutta and sixth order AdamsMoulton predictor-corrector integration algorithm55,56 as implemented in VENUS-MOPAC.57 At each time step Schro¨dinger’s
J. Phys. Chem. B, Vol. 103, No. 18, 1999 3693 equation is solved for the trimethylene intramolecular energies and forces using the appropriate MOPAC 7.042 subroutines. Once a converged SCF wavefunction has been obtained and configuration interaction included, the trimethylene atomic forces are determined analytically. By starting the SCF convergence using the optimized molecular orbital coefficients from the previous trajectory step, only 6-10 SCF iterations are required at each step. The intermolecular forces are the first derivatives of the LJ 12-6 energy functions with respect to the Cartesian coordinates and are added to the trimethylene intramolecular forces to give the total forces. The velocity scaling procedure of Berendsen et al.58 was used to maintain a constant Ar temperature. At each time step the Ar velocities are scaled by
[
λ) 1+
)]
(
∆t T -1 τ Tins
1/2
(4)
where ∆t is the trajectory time step, τ is a smoothing factor, T is the desired temperature and Tins is the instantaneous temperature
Tins )
2 (3N - 3)kb
N
1
mVi2 ∑ 2 i)1
(5)
where N is the number of Ar atoms, kb is the Boltzmann factor, m is the Ar mass, and Vi are the atomic velocities. If τ ) ∆t, the Ar temperature is set exactly to T at each time step, which may induce artificial noise or chaos into the simulation. For τ greater than ∆t, the Ar temperature approaches, but does not equal, T and the stochastic perturbations to the atomic velocities are not as severe. A larger value of τ also allows for local heating of the Ar bath in the vicinity of the trimethylene. Different temperature scaling rates were investigated by varying τ around 10 fs, and since the simulation results were insensitive to these variations, a τ of 10 fs (∆t ) 0.25 fs) was used for the simulations reported here. The insensitivity to variations in τ indicates that similar simulation results would be obtained with other algorithms, e.g., the Nose´-Hoover thermostat,59 for maintaining the temperature of the Ar bath. Trajectories were integrated with a step size of 0.25 fs, which in the absence of temperature scaling (τ f ∞), conserved the energy to four significant figures over the length of each trajectory. The trajectories were propagated until either cyclopropane or propene was formed. The criterion for cyclopropane formation was that the central CCC angle decreased to 75°. This is 19° smaller than at the variational cyclization transition state for a trimethylene energy of 54.6 kcal mol-1 and 14° smaller than that at 146.4 kcal mol-1 and facilitated the monitoring of TS recrossings. The criterion for propene formation was that a central CH bond was at least 2.3 Å with the H-atom being closer than 1.3 Å to a terminal carbon. Since these configurations lie on the propene side of the H-transfer barrier, this criterion also allows for evaluation of TS recrossings. Analysis of Trajectory Data. Collisional energy transfer was studied by monitoring trimethylene intramolecular and rotational energies over time. In particular, collisional energy transfer can be observed as a broadening of the trimethylene energy distributions, where each distribution is obtained from energies monitored over the trimethylene lifetime. Distributions were obtained by calculating energies every tenth trajectory step (i.e., every 2.5 fs) for all 150 trajectories in the ensemble and binning these data to yield histograms. The resulting distributions are averages over the lifetimes of the trimethylene molecules.
3694 J. Phys. Chem. B, Vol. 103, No. 18, 1999
Bolton et al.
In accordance with the no recrossing assumption of RRKM theory,52,53 cyclopropane:propene branching ratios, Ncyc:Nprop, were obtained from the first crossing of the cyclization or H-transfer barrier, i.e., Ncyc is the number of trajectories crossing the cyclization barrier and Nprop the number crossing the H-transfer barrier. Effects of recrossings were determined by comparing Ncyc:Nprop with the ratio of the number of trajectories forming cyclopropane and propene products. The branching ratios are normalized, i.e.,
BR )
(
)(
)
Ncyc Nprop : Ncyc + Nprop Ncyc + Nprop
(6)
Decay rate constants were also obtained from the first time that trajectories crossed the cyclization or H-transfer barrier. In accordance with RRKM theory,52,53 they are determined by fitting the simulated lifetime distribution, P(t), with a single exponential function
1 dN(t) ) ke-kt N(0) dt
P(t) ) -
(7)
where N(t) is the number of trajectories that have not reacted at time t and k is the rate constant. If the decay is biexponential, it is fit with the sum of two exponentials, i.e.,
S(t) )
N(t) ) Ae-k1t + Be-k2t N(0)
(8)
where A and B are the fraction of trajectories decomposing with rate constants k1 and k2, respectively. Effects of barrier recrossings are evaluated using eqs 7 and 8, but with t being the final time of barrier crossing before product formation, and comparing rate constants determined in this manner with those obtained from the first time of barrier crossing. Three properties that provide structural and dynamic information of the Ar environment are the reduced density, F*, the diffusion constant, D, and the radial distribution function, g(r).46,60 The reduced density is
F* )
Nσ3 ) Fσ3 V
(9)
where N is the number of Ar atoms, V is the volume of the periodic box, and σ3 is the excluded volume of an Ar atom, calculated from the LJ σ. The presence of the trimethylene, which slightly reduces the volume V, is ignored. The diffusion rate constant is46
D)
MSD 6t
(10)
where the mean square displacement (MSD) over time t ) tf ti is N
(rj(tf) - rj(ti))2> ∑ j)1
< MSD )
N
(11)
The sum in eq 11 is over all Ar atoms and rj(ti) and rj(tf) are the positions of atom j at the initial and final times ti and tf, respectively. The averaging is over all time origins, ti, along the trajectory with t ) tf - ti held fixed. Since relatively long trajectories are required to obtain the diffusion constants and calculating AM1-SRP forces is computationally demanding,
TABLE 1: Cyclopropane:Propene Branching Ratios for G* ) 0.52 E ) 54.6 kcal mol-1
E ) 146.4 kcal mol-1
T/K
EMSa
barrierb
EMS
barrier
100 200 300 500 800 1000
0.86:0.14 0.82:0.18 0.88:0.12 0.83:0.17 0.83:0.17 0.85:0.15
0.77:0.23 0.83:0.17 0.81:0.19 0.88:0.12 0.81:0.19 0.79:0.21
0.45:0.55 0.47:0.53 0.50:0.50 0.52:0.48 0.57:0.43 0.54:0.46
0.43:0.57 0.47:0.53 0.47:0.53 0.40:0.60 0.47:0.53 0.44:0.56
a Trajectories were initialized with EMS. sampling. b Trajectories were initialized with normal mode sampling at the H-transfer barrier.
these simulations were done in the absence of trimethylene. The presence of trimethylene is not expected to significantly affect the diffusion rates since the radial distribution function, which measures the structure of the Ar environment, is the same in the presence or absence of trimethylene. To determine the diffusion rate constant, a 100 ps trajectory was propagated with a step size of 1 fs. The radial distribution function is46
g(r) )
V
δ(r - rij)> ∑i ∑ j*i
