© Copyright 1997 by the American Chemical Society
VOLUME 101, NUMBER 40, OCTOBER 2, 1997
LETTERS A Combined Car-Parrinello QM/MM Implementation for ab Initio Molecular Dynamics Simulations of Extended Systems: Application to Transition Metal Catalysis Tom K. Woo,† Peter M. Margl,† Peter E. Blo1 chl,‡ and Tom Ziegler*,† Department of Chemistry, UniVersity of Calgary, 2500 UniVersity DriVe, N.W., T2N 1N4 Calgary, Alberta, Canada, and IBM Research DiVision, Zu¨ rich Research Laboratory, Sa¨ umerstrasse 4, CH-8803 Ru¨ schlikon, Switzerland ReceiVed: May 29, 1997; In Final Form: August 4, 1997X
We test a new implementation of the combined quantum mechanics and molecular mechanics (QM/MM) methodology applied within the Car-Parrinello framework to perform ab initio molecular dynamics simulations of extended systems. The novel method is applied to determine the free energy barrier of the chain termination process in a nickel diimine based ethylene polymerization catalyst of the type (ArNdC(R)-C(R)dNAr)Ni(II)-R′+, where R ) Me and Ar ) 2,6-C6H3(i-Pr)2. In this combined QM/MM ab initio molecular dynamics simulation, the Ni diimine core was treated at the Becke88-Perdew86 DFT level while the large substituted aryl rings were treated by the AMBER molecular mechanics force field. A 39 000 time step slow growth simulation of the termination process at 300 K has been performed providing a free energy barrier of ∆Fq ) 14.8 kcal/mol by thermodynamic integration. This is in excellent agreement with the experimental termination barrier of ∆Gq ≈ 16 kcal/mol. Without the bulky ligands, the analogous pure QM simulation provided a free energy barrier of ∆Fq ) 9.8 kcal/mol.
I. Introduction There have been many approaches introduced, such as multiple time scale techniques1 and linear scaling methods,2,3 to increase the efficiency of Car-Parrinello4 based ab initio molecular dynamics simulations to treat extended systems. Taking a different route, we have implemented the combined quantum mechanics and molecular mechanics (QM/MM) method5 of Singh6 and Field7 into the Car-Parrinello ab initio molecular dynamics framework (technical details to be published elsewhere8). The combined QM/MM method has recently gained significant attention6,7,9-13 because of its potential to simulate * Author to whom correspondence should be addressed. E-mail: ziegler@ zinc.chem.ucalgary.ca. † University of Calgary. ‡ Zu ¨ rich Research Laboratory. X Abstract published in AdVance ACS Abstracts, September 15, 1997.
S1089-5647(97)01729-X CCC: $14.00
large systems in a detailed, yet efficient, manner. Despite the attention the methodology has received, there are still many issues to be resolved in the combined QM/MM field, and this is particularly true when the QM/MM boundary crosses covalent bonds. Here it is our intent to present the first calculations using a combination of QM/MM coupling with dynamical free energy calculations using the Car-Parrinello molecular dynamics approach. Recently, we have applied a “static” QM/MM approach to study a Ni(II) diimine based olefin polymerization catalyst11 and found that the QM/MM approach provided results that were in excellent agreement with experiment. Specifically, the method was applied to a Ni(II) diimine catalyst of the type (ArNdC(R)-C(R)dNAr)Ni(II)-R′+, 1, where R ) Me and Ar ) 2,6-C6H3(i-Pr)2, which was discovered by Brookhart and co-workers14,15 and is currently being developed for com© 1997 American Chemical Society
7878 J. Phys. Chem. B, Vol. 101, No. 40, 1997 SCHEME 1
mercialization by DuPont.16 In this polymerization system the bulky aryl groups play a crucial role, since without the bulky substituents the catalyst acts only as a dimerization catalyst due to the favorability of the β-elimination chain termination process. From the structure of the catalyst, it is evident that the bulky aryl substituents partially block the axial coordination sites of the Ni center, and it was proposed by Johnson et al.14 that it is likely this steric feature which impedes the termination relative to the insertion process. Our “static” QM/MM calculations11 confirmed this notion. It was found that the termination transition state has both axial coordination sites occupied whereas the insertion transition state has only the equatorial sites occupied. In that same work,11 the enthalpic termination barrier was calculated to be ∆Eq ) 18.6 kcal/mol, which is in satisfactory agreement with the experimental17 free energy barrier of ∆Gq ≈ 16 kcal/mol.18 In this study we intend to apply and validate our combined ab initio molecular dynamics Car-Parrinello QM/MM methodology by determining the free energy barrier of the hydrogen transfer to the monomer chain termination process (Scheme 1). The bulky R ) Me and Ar ) 2,6-C6H3(i-Pr)2 will be treated by a molecular mechanics potential, whereas the Ni diimine core including the growing chain and monomer will be treated by a density functional potential.
Letters interactions were neglected. Torsional and angle parameters not available in the AMBER95 force field were replaced with similar parameters from the same force field.22 The details of the ab initio molecular dynamics methodology are essentially the same as in previous works23-25 with minor modifications described below. The accuracy of the CP-PAW method for geometries and energetics has been established by Blo¨chl.19 Periodic boundary conditions were used with a unit cell spanned by the lattice vectors ([0.0 8.5 8.5][8.5 0.0 8.5][8.5 8.5 0.0]). All simulations were performed with the local density approximation in the parametrization of Perdew and Zunger26 with gradient corrections due to Becke27 and Perdew.28,29 The frozen core approximation was utilized with an Ar core used for Ni and a He core used for the first-row elements. A Nose´ thermostat30 was used to maintain an average temperature of 300 K. Separate thermostats were used for the QM and MM regions. Masses of the nuclei were rescaled to 50.0 amu for Ni, 2.0 amu for N and C, and 1.5 amu for H. Since we do not discuss time-dependent properties and since configurational ensemble averages remain unchanged under a rescaling of the masses, this technique is appropriate.31 The free energy barriers were calculated using the “slow growth” technique, which has been detailed and demonstrated on other olefin polymerization catalysts.23-25 With the slow growth technique, a reaction coordinate (RC) is constrained during the dynamics and slowly varied from a value characteristic of the initial state to a value characteristic of the final state of interest. The Helmholtz free energy difference ∆F between two arbitrary points λ ) 0 and λ ) 1 along the RC is determined as
∆F )
II. Computational Details All reported molecular dynamics simulations were carried out with the Car-Parrinello projector augmented wave (CPPAW) code developed by Blo¨chl19 and extended8 to include the AMBER9520 molecular mechanics force field in a manner similar to the method introduced by Singh and Kollman.6 We have used dummy hydrogen atoms to cap the exposed valences due to the covalent bonds that cross the QM-MM boundary. The dummy atoms are replaced by the so-called link atoms, which are denoted with an asterisk in structure 1. The quantum mechanical calculations were performed on the 26-atom Ni diimine molecule [(HNdC(H)-C(H)dNH)Ni-propyl + ethene]+ for both the combined QM/MM simulation and the pure QM simulation. An augmented AMBER95 force field20 was utilized to describe the molecular mechanics potential. Employing the AMBER atom type labels as described in ref 20, the diimine carbon was assigned with atom type “CM” parameters, the diimine N with “N2”, aryl ring carbon atoms with “CA”, aryl ring hydrogen atoms with “HA” and the remaining carbon and hydrogen atoms of the MM region with “CT” and “HC”, respectively. The reacting ethene monomer was assigned with sp2 “C” van der Waals parameters. Alkyl carbon and hydrogen atoms of the active site were assigned “CT” and “HC” van der Waals parameters, respectively. Ni was assigned the "Ni4+2" van der Waals parameters of Rappe´’s UFF.21 Nonbonded interactions between the QM and MM regions were treated with a molecular mechanics van der Waals potential. Electrostatic
∫01〈fλ〉|λ,T dλ
(1)
where fλ is the force on the reaction coordinate and the number of samples at each point λ equals 1 in the slow growth limit. The “midplane” slow growth reaction coordinate32 was utilized. This reaction coordinate is defined as the ratio r/R, where R is the length of the vector between C2 and Cβ (see Scheme 1 for labeling) and r is the length of projection of the Cβ-Hβ bond vector onto the vector R. This reaction coordinate constrains the transferring H atom to lie on the plane perpendicular to and passing through the end point of the vector r. The midplane reaction coordinate has been previously demonstrated.23 In the simulations that are presented here, the total scan time chosen was about 39 000 time steps. All calculations were performed on an IBM RS/6000 model 3CT workstation with 64 MB of memory. Each of the simulations required 2 weeks of CPU time. III. Results and Discussion We have performed a unidirectional scan of the chain termination process33 (Scheme 1) commencing from the resting state π-complex to the approximate transition state. The computational effort has been focused on simulating the first half of the hydrogen transfer process in order to estimate the free energy barrier for the process. Two simulations have been performed, one that models the bulky substituents with a combined QM/MM potential and one pure QM simulation that neglects the bulky substituents altogether. In these simulations a midplane reaction coordinate which is defined in the Computational Details section, is utilized. In the initial π-complex, the value of the RC is initially small (approximately 0.18), indicating that the H atom to be transferred is still bound to the Cβ of the alkyl chain. As the reaction progresses the midplane value increases. When the value of the reaction coordinate is
Letters
J. Phys. Chem. B, Vol. 101, No. 40, 1997 7879
Figure 1. Relative free energy as a function of the midplane reaction coordinate for the chain termination process (Scheme 1). The solid profile refers to the combined QM/MM simulation whereby the bulky alkyl and aryl substituents are treated by a molecular mechanics force field. The dotted profile refers to the pure QM simulation on the model system without the bulky substituents. The position of the snapshot structure 2 is depicted.
0.5, this indicates that the transferred Hβ atom lies on a plane midway between the two atoms Cβ and C2. Figure 1 displays the free energy as a function of the midplane reaction coordinate for both simulations. The plots reveal that the bulky ligands significantly increase the free energy barrier for chain termination. The pure QM simulation, which does not account for the bulky ligands, provides a free energy barrier of ∆F ) 9.8 kcal/mol, whereas the combined QM/MM simulation provides a barrier of 14.8 kcal/mol. The combined QM/MM free energy barrier of 14.8 kcal/mol is in excellent agreement with the experimental18 free energy barrier of ≈16 kcal/mol. We note that this simulation does not include any quantum effects, namely the correction from the zero-point energy, which we expect to be small.23 These values are also in reasonable agreement with previously calculated energy barriers calculated from “static” simulations.11,34 Figure 1 reveals that there is a plateau in the free energy profile for the simulation without the bulky substituents, indicating the presence of a transient double olefin hydride complex where the RC is approximately 0.5. When the bulky ligands are added with the molecular mechanics potential, the plateau vanishes, leaving a single distinct termination transition state. This interesting feature of the energy profiles35 was also observed with static calculations.11,34 Selected geometric and energetic quantities are plotted in Figure 2 as a function of the slow growth reaction coordinate for the QM/MM simulation. A snapshot structure of the QM/ MM simulation at the approximate transition state is illustrated in Figure 3. Figure 3a traces the Cβ-Hβ, C2-Hβ, and Ni-Hβ bond distances to portray the transfer process (see Scheme 1 for labeling convention). As the reaction proceeds from the π-complex to the transition state, the Cβ-Hβ bond slowly increases as the C2-Hβ distance correspondingly decreases. The transition state occurs at a midplane value of approximately 0.52, where the hydrogen is roughly midway between the two carbon atoms. As the hydrogen is transferred, it is pulled into the Ni atom as demonstrated by the Ni-Hβ distance, which varies from roughly 1.75 Å in the resting state to roughly 1.5 Å in the transition state region. The distances plotted in Figure 2b show the expansion of the active site during the hydrogen transfer process, which peaks at the transition state region. Since the ethene and propyl moieties of the active site occupy the axial positions as shown in Figure 3, the expansion results in an increased steric
Figure 2. Selected structural and energetic quantities as a function of the reaction coordinate from the combined QM/MM simulation: (a), (b) selected distances; (c) the molecular mechanics van der Waals interaction energies between the two bulky aryl groups and the active site ethene and propyl groups relative to the resting state value. The energies plotted have been “smoothed” and are the running average over 500 time steps. The vertical line running through all three plots represents the position of the snapshot structure 2.
Figure 3. Snapshot structure from the transition state region of the combined QM/MM molecular dynamics simulation. The molecular mechanics atoms are ghosted for clarity, while the dummy hydrogen atoms are omitted. The value of the midplane reaction coordinate in 2 is RC ) 0.526. Bond distances reported are in angstroms.
interaction with the isopropyl substituents of the aryl rings. In the present QM/MM model, the steric interaction between these groups is accounted for by the molecular mechanics van der Waals potential. Plotted in Figure 2c is the sum of all the van der Waals interaction energies between the aryl rings and the active site moieties as a function of the reaction coordinate (relative to the initial resting state value). Figure 2c, therefore, expresses the increase in the steric interactions between the aryl
7880 J. Phys. Chem. B, Vol. 101, No. 40, 1997 rings and the active site ethene and propyl fragments as the reaction progresses. The interaction involving the propyl group peaks at the transition state region and amounts to roughly 1 kcal/mol. On the other hand, the interaction involving the ethene molecule peaks somewhat prior to the transition state region at a reaction coordinate value of roughly 0.35, and at the transition state the steric interaction has diminished somewhat, amounting to only ≈∆Evdw ) 0.5 kcal/mol. A more detailed analysis of the role of the bulky substituents in this polymerization system has been presented in our “static” QM/MM study.11 IV. Conclusions A new implementation to carry out Car-Parrinello ab initio molecular dynamics simulations of extended systems using a combined quantum mechanics and molecular mechanics potential is presented. The applicability and efficiency of the approach to study transition metal catalysis is demonstrated. In particular, we have examined the chain termination process in the Brookhart-type Ni(II) diimine olefin polymerization catalyst and found that the calculated free energy barrier of ∆Fq ) 14.8 kcal/mol at 25 °C is in excellent agreement with the experimental result of ∆Gq ≈ 16 kcal/mol at 0 °C. In this simulation the QM/MM boundary crosses several covalent bonds, such that the Ni diimine core is treated by the DFT potential while the large bulky substituents are treated with the AMBER molecular mechanics force field. We are currently applying the methodology to problems more amenable to the dynamics scheme such as the monomer capture barrier in the Brookhart catalyst system,36 where we are utilizing a QM/MM multiple time step methodology.8 Acknowledgment. This work has been supported by the National Sciences and Engineering Research Council (NSERC) of Canada, as well as by the donors of the Petroleum Research Fund, administered by the American Chemical Society (ACSPRF No. 31205-AC3) and by Novacor Research and Technology Corporation (NRTC) of Calgary. The authors acknowledge Professor M. Brookhart and Dr. L. Deng for the valuable discussions. T.K.W. thanks NSERC, the Alberta Heritage Scholarship Fund, and the Izaak Walton Killam memorial foundation for graduate fellowships. References and Notes (1) Tuckerman, M. E.; Parrinello, M. J. Chem. Phys. 1994, 101, 1316. (2) Galli, G.; Parrinello, M. Phys. ReV. Lett. 1992, 69, 3547. (3) Mauri, F.; Galli, G.; Car, R. Phys. ReV. B 1993, 47, 9973. (4) Car, R.; Parrinello, M. Phys. ReV. Lett. 1985, 55, 2471. (5) Aqvist, J.; Warshel, A. Chem. ReV. (Washington, D.C.) 1993, 93, 2523. (6) Singh, U. C.; Kollman, P. A. J. Comput. Chem. 1986, 7, 718. (7) Field, M. J.; Bash, P. A.; Karplus, M. J. Comput. Chem. 1990, 11, 700. (8) Woo, T. K.; Margl, P. M.; Blo¨chl, P. E.; Ziegler, T., in preparation. (9) Bash, P. A.; Field, M. J.; Davenport, R.; Ringe, D.; Petsko, G.; Karplus, M. Biochemistry 1991, 30, 5826.
Letters (10) Tun˜o´n, I.; Martins-Costa, M. T. C.; Millot, C.; Ruis-Lo´pes, M. F. J. Chem. Phys 1997, 106, 3633. (11) Deng, L.; Woo, T. K.; Cavallo, L.; Margl, P. M.; Ziegler, T. J. Am. Chem. Soc. 1997, 119, 6177. (12) Maseras, F.; Morokuma, K. J. Comput. Chem. 1995, 16, 1170. (13) Gao, J. In ReViews in Computational Chemistry; Lipkowitz, K. B., Boyd, D. B., Eds.; VCH: New York, 1996; Vol. 7. (14) Johnson, L. K.; Killian, C. M.; Brookhart, M. J. Am. Chem. Soc. 1995, 117, 6414. (15) Johnson, L. K.; Mecking, S.; Brookhart, M. J. Am. Chem. Soc. 1996, 118, 267. (16) Haggin, J. Chem. Eng. News 1996, 74 (Feb 5), 6. (17) Reference 14. Polymerization of 1.6 × 10-6 mol of [(ArNdC(R)C(R)dNAr)Ni(CH3)(OEt2)]+[B(3,5-C6H3(CF3)2)4]-, where Ar ) 2,6C6H3(i-Pr)2 and R ) Me in 100 mL of toluene at 0 °C for 15 min. (18) The free energy barrier of olefin insertion is experimentally estimated to be ∆Gq ) 10-11 kcal/mol (Professor Maurice Brookhart, Department of Chemistry, University of North Carolina at Chapel Hill, a private communication.) The weight-average molecular weight, Mw, of 8.1 × 105 g/mol provides an estimate for the ratio of termination events to insertion events of 1:28900. Applying Boltzmann statistics to this ratio gives a ∆∆Gq of 5.6 kcal/mol and an estimate for the termination barrier of 15.516.5 kcal/mol. In calculating ∆∆Gq we have assumed that every hydrogen transfer event leads to the loss of the chain and, consequently, chain termination. (19) Blo¨chl, P. E. Phys. ReV. B 1994, 50, 17953. (20) Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Gould, I. R., Jr.; K. M. M.; Ferguson, D. M.; Spellmeyer, D. C.; Fox, T.; Caldwell, J. W.; Kollman, P. A. J. Am. Chem. Soc. 1995, 117, 5179. (21) Rappe´, A. K.; Casewit, C. J.; Colwell, K. S.; Goddard, W. A., III; Skiff, W. M. J. Am. Chem. Soc. 1992, 114, 10024. (22) CM-N2-CA, Ni-N2-CA, and N2-CM-CT angle parameters were replaced by CM-N*-CT, H-N2-CA, and CM-CM-CT angle parameters of the AMBER95 force field, respectively. The X-N2-CM-X torsional parameters was replaced by the X-N2-CA-X parameter of the AMBER95 force field. The X-Ni-N2-X torsional parameter was assigned a barrier of zero. (23) Woo, T. K.; Margl, P. M.; Blo¨chl, P. E.; Ziegler, T. Organometallics 1997, 16, 3454. (24) Woo, T. K.; Margl, P. M.; Lohrenz, J. C. W.; Blo¨chl, P. E.; Ziegler, T. J. Am. Chem. Soc. 1996, 118, 13021. (25) Margl, P.; Lohrenz, J. C. W.; Blo¨chl, P.; Ziegler, T. J. Am. Chem. Soc. 1996, 118, 4434. (26) Perdew, J. P.; Zunger, A. Phys. ReV. B 1981, 23, 5048. (27) Becke, A. Phys. ReV. A 1988, 38, 3098. (28) Perdew, J. P. Phys. ReV. B 1986, 34, 7406. (29) Perdew, J. P. Phys. ReV. B 1986, 33, 8822. (30) Nose´, S. Mol. Phys. 1984, 52, 255. (31) De Raedt, B.; Sprik, M.; Klein, M. L. J. Chem. Phys. 1984, 80, 5719. (32) Blo¨chl, P. E. Unpublished work. (33) Animations of the simulations presented here are available at our research group’s web site at . (34) Deng, L.; Margl, P. M.; Ziegler, T. J. Am. Chem. Soc. 1997, 119, 1094. (35) In both simulations depicted in Figure 1, the reverse scan from the transition state back to the resting state has not been performed. The difference between forward and reverse barriers (hysteresis) is a measure of the statistical error incurred with the slow growth method. It is our experience from other small compounds and a reaction coordinate where the force of constraint pushing the system across the barrier is nearly parallel to the free energy gradient that the hysteresis effect is small (1-2 kcal/ mol).23 In the present case, one might expect a larger hysteresis due to the low-frequency modes of the ligands in the QM/MM simulation. In this case, however, the good agreement in the free energy barriers with static QM/ MM calculations indicates that the effect is small as well. (36) Woo, T. K.; Margl, P. M.; Deng, L.; Ziegler, T., in preparation.
