Figure 5.3. The relative free energy as a function of the mid-plane reaction coordinate for the chain termination process (Scheme 5.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 33 is depicted.
Figure 5.3 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 experimental 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.133 These values are also in reasonable agreement with previously calculated energy barriers calculated from "static" simulations.29,59
Figure 5.3 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 ~ 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 profiles was also observed with static calculations first presented in Chapter 3.
Figure 5.4. Selected structural and energetic quantities as a function of the reaction coordinate from the combined QM/MM simulation. (a) and (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 33.
Selected geometric and energetic quantities are plotted in Figure 5.4 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 5.5. Figure 5.4a traces the C-H, C2-H and Ni-H bond distances to portray the transfer process (see Scheme 5.1 for labeling convention). As the reaction proceeds from the -complex through to the transition state, the C-H bond slow 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.
33
Figure 5.5. 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 33 is RC=0.526. Bond distances reported are in Ångstroms.
The distances plotted in Figure 5.4b shows 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 occupied the axial positions as shown in Figure 5.5, the expansion results in an increased steric interaction with the iso-propyl 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 5.4c 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 5.4c, therefore, expresses the increase in the steric interactions between the aryl 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 about ∆Evdw = 0.5 kcal/mol.
The applicability and efficiency of the combined QM/MM ab initio molecular dynamics approach to study transition metal catalysis has been demonstrated. The PAW QM/MM method will allow for ab initio level molecular dynamics simulations of extended systems to be performed. For our research in olefin polymerization catalysis, this will allow us to incorporate realistic ligand effects into our molecular dynamics simulations. We feel this will be an important advance in computational modeling of single-site catalysts since there is a trend towards larger more complex ligand systems, particularly in the area of stereospecific -olefin polymerization.88,123,167
5.5 QM/MM Multiple Time Step Dynamics
Interest in the linear scaling of electronic structure calculations has surged with the recent developments from the groups of Yang,6,168 Head-Gordon5 and Scuseria.4 Although these techniques allow the time of an electronic structure calculation to scale linearly with the size of the system, they do not address the problem of non-linear scaling of the geometry optimization or sampling of configuration space. In this section, we address the issue within the framework of the QM/MM ab initio molecular dynamics methodology. Our goal is to increase the configurational sampling of the MM region without significantly increasing the computational effort put into the QM region.
A QM/MM simulation will most often involve a small QM region embedded within a much larger MM domain, such as with the simulation of the active site chemistry within an enzyme. Associated with the larger size of the MM region is a potential energy surface that is more complex and one which is likely to posses a higher degree of configurational variability. This necessitates an increased degree of sampling in the MM region in order to obtain meaningful ensemble averages from a simulation. In principle, increased sampling of the MM region can be achieved by the technique of mass rescaling. In classical molecular dynamics, configurational averages do not depend on the masses of the nuclei,◊ and therefore the true nuclear masses can be replaced with more convenient values.150,169 Rescaling of nuclear masses to smaller values, which allows the particles to move faster, increases the rate of configurational sampling which is proportional to the inverse square root of the mass, m-1/2, of the particles. Therefore, by decreasing the masses of the MM nuclei relative to those in the QM region, one can differentially sample the configuration space of the two regions. However, if the same time step is used for both regions a limit to the practical over sampling is quickly reached. This is because as we rescale the masses and increase the speed of the MM nuclei, a smaller time step is required in order to accurately integrate the equations of motion of the fast moving system. Thus, when the QM and MM regions are propagated simultaneously, the "slow" QM region will be propagated with unnecessarily small time steps. This is undesirable because the QM derived forces will be changing slowly and much time will be wasted recalculating these forces at each of the small time steps.
Figure 5.6. Representation of a simple multiple time step scheme within the Car-Parrinello framework. The wave function and the QM subsystems are propagated a time step ∆t, while the MM subsystem is over sampled with a smaller time step ∆t/n.
One way to overcome this problem is to propagate the two regions asynchronously with what has been termed multiple time-step molecular dynamics.170 The multiple time-step method in the Car-Parrinello QM/MM framework is represented in Figure 5.6 where the wave function and the QM nuclei are propagated with a large time step, ∆t, while the less computationally demanding MM region is n times over sampled with a time step of ∆t/n. Thus, the multiple time step method in combination with mass rescaling can be applied to the QM/MM molecular dynamics methodology as to differentially sample the QM and MM regions, thereby increasing the sampling of the MM domain without increasing the computational expenditure in the QM region.
The multiple time step techniques have been developed since 1978170 to more efficiently treat systems with high and low-frequency motion and/or short and long-range forces. However, the simple algorithm often used (as depicted in Figure 5.6) has no rigorous basis in theory.171 Recently, Tuckerman and co-workers158,172,173 have pioneered the development of multiple time-step methods where there is a rigorous separation of time scales. We have adopted the reversible multiple time step algorithm of Tuckerman et al.158 which allows for numerically stable, energy conserving multiple time step molecular dynamics to be performed. In our implementation the original formalism which was based on the velocity Verlet propagation algorithm has been modified to accommodate the standard Verlet algorithm used in the PAW program.
Share with your friends: |