Nucleophilic Substitution Dynamics: Comparing Wave Packet

Article pubs.acs.org/JPCA

Nucleophilic Substitution Dynamics: Comparing Wave Packet Calculations with Experiment M. Kowalewski,*,†,∥ J. Mikosch,‡ R. Wester,§ and R. de Vivie-Riedle† †

Department of Chemistry, Ludwig-Maximilians-Universität, D-81377 Munich, Germany Max-Born-Institute, Max-Born-Strasse 2A, D-12489 Berlin, Germany § Institut für Ionenphysik und Angewandte Physik, Universität Innsbruck, Technikerstraße 25/3, 6020 Innsbruck, Austria ∥ Division of Scientific Computing, Department of Information Technology, Uppsala University, SE-75105 Uppsala, Sweden ‡

S Supporting Information *

ABSTRACT: The reactive collision of chloride anions and methyl iodide molecules forming iodide anions and methyl chloride is a typical example of a concerted bimolecular nucleophilic substitution (SN2) reaction. We present wave packet dynamics calculations to investigate quantum effects in the collinear gas phase reaction. A new type of reduced coordinate system is introduced to allow for an efficient solution of the time-dependent Schrödinger equation on an ab initio potential energy surface. The reduced coordinates were designed to study the direct rebound mechanism under the Walden inversion. Especially the suppressed direct rebound mechanism at low collision energies, the quantum effects of the initial state preparation and the influence of the CH3 inversion mode are addressed. The internal energy distributions of the molecular product are evaluated from the wave packet calculations and compared to experimental results obtained with crossed-beam velocity map ion imaging. The observed reactivity is discussed in light of a dynamical barrier, a concept that is illustrated by the wave packet dynamics.

1. INTRODUCTION Nucleophilic substitution reactions in solution are an essential everyday tool in organic and inorganic synthesis. This reaction class is widely studied in the gas phase1−3 to discern the intrinsic reaction dynamics from the solvent effects. An interesting aspect of these studies is how the dynamics on the femtosecond time-scale relates to phase-space-based descriptions of the reaction outcome. The theoretical description of nucleophilic substitution reactions has been based on three main types of methods, statistical calculations of reaction rates,4 classical-trajectory chemical dynamics simulations,3 and quantum scattering calculations.5,6 Experimentally, most results have been obtained with the flow tube technique,7 with guided ion beam mass spectrometry,8 and with crossed-beam reactive scattering.9 Chemical dynamics simulations recently achieved good agreement with a crossed-beam experiment.10,11 Here we study the SN2-reaction Cl− + CH3I → ClCH3 + I−

of 0.4 eV a long-lived metastable complex was observed, whereas starting around 1 eV the classical direct rebound mechanism, which occurs along a collinear reaction path, was found to dominate. In addition an indirect roundabout reaction mechanism was discovered at 1.9 eV collision energy.9 In this work we study the direct rebound mechanism of the SN2 reaction by quantum dynamical wave packet calculations, allowing us to investigate the role of quantum effects in the initial state preparation and during the reaction. Bond breaking and formation occurs on the subpicosecond time scale and above. We inspect the dominant reaction mechanism: the socalled Walden inversion where the nucleophile attacks from the backside under the inversion of the CH3 group. We use specially designed reactive coordinates, which allow for an efficient description of this reaction, including the role of one of the spectator modes (CH3 umbrella motion). In addition the influence of C−I bond activation on the product outcome is tested. The associated energy redistribution among the internal molecular and translational degrees of freedom is followed in detail. It is known from experiments9 and previous theoretical studies13 that in this direct rebound mechanism the reactivity is suppressed even though the transition state (TS) is lower than

(1)

which is exoergic by about 0.5 eV. Its reaction probabiliy at room temperature has been derived to be about 9%, determined by the ratio of the measured versus the capturelimited rate coefficient.12 The dynamics of reaction 1 have been studied by crossed-beam ion imaging and on-the-fly molecular dynamics simulations.9 These studies revealed the appearance of three different reaction mechanisms: at a low collision energy © 2014 American Chemical Society

Received: April 23, 2014 Revised: June 3, 2014 Published: June 3, 2014 4661

dx.doi.org/10.1021/jp503974u | J. Phys. Chem. A 2014, 118, 4661−4669

The Journal of Physical Chemistry A

Article

the reactants. This question will be addressed in more detail by our wave packet calculations. The paper is organized as follows. In the next section a short summary of the experimental procedure is given. In section 3 the reduction scheme for the generation of the reactive coordinates and the setup for the quantum dynamical calculation is explained. In section 4 we present our calculations of the reaction probability, the internal energy distribution, which is compared to experimental data, and the influence of the CH3 umbrella motion. In section 5 we discuss the results and develop an intuitive picture of the dynamical barrier based on the data from the wave packet calculations.

Figure 1. Geometry for the gas phase collinear collision reaction of methyl iodide and chloride. The figure shows the choice of coordinates, which are of major importance for the reaction: The C−Cl and the C−I bonds are directly involved in the bond formation and bond cleavage process. The angle of the hydrogen atoms θ is called umbrella angle. It follows the reaction and undergoes an inversion.

2. EXPERIMENTAL SETUP The experimental data have been obtained with a crossed-beam ion imaging spectrometer that is designed to study ion− molecule reactive scattering.14 The experimental setup and the data acquisition have been described in previous publications.9,11,15,16 Here we only summarize the essential aspects of the experiment and the data analysis. Atomic chloride anions are produced in a pulsed electron impact ion source composed of a pulsed supersonic molecular beam that is irradiated by electrons with a few tens of electronvolt kinetic energy. Methyl tetrachloride (CCl4) in the supersonic beam is fragmented into Cl− and CCl3 by dissociative electron attachment. The chloride anions are accelerated toward the interaction region, where they cross a second pulsed supersonic jet, which provides the neutral CH3I reactants seeded at 10% in helium. The relative collision energy is tuned by changing the electric potential in the interaction region and thereby gently decelerating the ion packet prior to the collision with the neutral target beam. A velocity map imaging spectrometer is used to project the product ions emerging from the reactive collisions onto a position sensitive imaging detector with single-ion detection efficiency. The electric field of the spectrometer is pulsed on after the ion and neutral reactant packets have crossed each other to allow the collision to happen in an effectively field-free region. The imaging detector is operated in the slice imaging mode,17−19 which only detects product ions which emerge in the plane of the reactant beams. This yields the angle- and velocity-differential scattering distributions as a function of the relative collision energy. To monitor the collision energy, also the velocity distributions of the reactant beams are imaged. For this the neutral beam is ionized by electron impact. The product ion velocity distributions are transformed into the center-of-mass frame of reference. By integration over the scattering angle one obtains the distribution of translational energies. Using the known exoergicity of reaction 1 yields the distribution of the total internal energy for the reaction product CH3Cl.9

and the methyl iodide the most important coordinates are the iodine−carbon distance (rC−I) and the chlorine carbon distance (rC−Cl). During the reaction the angle of the H atoms in the CH3 group (θ) also changes, which is referred to as umbrella motion. The angle θ is defined as the angle between the H atoms and the plane perpendicular to the symmetry axis including the C atom. A negative sign of θ corresponds to the reactants and a positive sign corresponds to the products. The asymmetric CH3 deformations as well as the CH-stretching motion is neglected in this model. With the choice of the three internal coordinates rC−I, rC−Cl, and θ a quantum dynamical calculation on a regular spatial grid would be possible but not very efficient. The related potential energy surface (PES) would include high energy regions describing fragmentation processes, which are of no relevance at collision energies below 3.5 eV. Thus, large areas in the numerical propagation procedure would be included, which the wave packets will never explore. It is computationally more efficient to focus on a configuration space close to the minimum energy path of the bimolecular substitution reaction and to use the intrinsic reaction coordinate20 (IRC). The concept of an IRC related coordinate system was first described by Miller et al.21 in the reaction path Hamiltonian formalism and has been extended to fit different scenarios.22−24 In the following an IRC adapted coordinate system is described, which allows for an optimal numerical description in the wave packet simulations. The idea is to create a new set of coordinates (qi, i ∈ {1, 2, 3}), which includes the IRC as one coordinate. These coordinates are called further on IRC parallel-adapted coordinates (IRC-pad). For the transformation the internal coordinates rC−I, rC−Cl, and θ are expressed as functions of the IRC-pad coordinates. The mapping is unique and works as follows: The IRC coordinate is labeled as q1 and includes the TS and the close-by local minima formed due to the anion−dipole interaction. The TS is by definition the origin of the new coordinate system (q1, q2, q3). In q1 all three internal coordinates are contained. Vice versa, the internal coordinates can be regarded as a function of q1, e.g., rC−Cl ′ (q1). The coordinate q2 is defined as a parallel shift of the IRC path performed in the rC−I − rC−Cl (Figure 3) plane. The mapping is then obtained by using the functional dependence of r′C−Cl and rC−I ′ on q1 as auxiliary coordinate:

3. QUANTUM DYNAMICAL MODEL The full-dimensional description of the given reaction requires 15 degrees of freedom, which is still out of reach for reactive processes in a full quantum model. Thus, a reduced number of suitable reactive coordinates has to be chosen. We focus on the question of the essential coordinates needed to reproduce the features of the gas phase collision under the typical backside attack. A collinear geometry as depicted in Figure 1 is chosen, which implies a strict C3v symmetry. For the collision of the chloride 4662

rC″− Cl(q1 ,q2) = rC′ − Cl(q1) + q2

(2)

rC″− I(q1 ,q2) = rC′ − I(q1) + q2

(3)



Article

The third coordinate q3 is related to a parallel shift of the IRC path along θ: θ″(q1 ,q3) = θ′(q1) + q3

(4)

The super script ″ denotes that the internal coordinates are expressed in dependence of q1, q2, and q3. The choice of the parallel shifts creates an invertible bijective transformation. The IRC path and the PES is determined in a quantum chemical calculation (MP2/ECP46MWB/aug-cc-pVDZ25−27 MOLPRO28). Comparisons with stationary point calculations on the G2(+) level of theory29 show that the method chosen has sufficient accuracy. The obtained IRC path is then parametrized by a fit with the following functions to allow for a flexible mapping:

Figure 2. Plot of the coordinate mapping functions from eqs 5 and 6. ′ (q1) and rC−I ′ (q1) and (b) the sigmoid (a) shows the functions rC−Cl function for θ′(q1). For q1 → −∞ the functions describe an approach of the chloride anion to the methyl iodide. The umbrella angle θ and the C−I bond rC−I are in their equilibrium position, whereas the C−Cl distance decreases linearly. The opposite is described for the values q1 → ∞ for a methyl chloride with a departing iodide (the umbrella angle is adjusted to the equilibrium position of CH3Cl).

rX″(q1 ,q2) = rX′ (q1) + q2 = −aX (q1 − qX,0) + rX,0 1 ln(1 + eaX bX (q1− qX,0)) + q2 + bX (6)

where X ∈ {C−Cl, C−I} (for the fitted parameters aθ, bθ, qθ,0, θ0, aX, bX, qX,0, rX,0 see Supporting Information). The mapping functions in eqs 5 and 6 are chosen to describe the asymptotic behavior correctly and to provide a good fit in the transition region. The angle θ″ of the umbrella motion is well described by a sigmoid function (eq 5), whereas for the halogen-carbon distance r″X the antiderivative of a sigmoid function describes the desired asymptotic properties (eq 6). The functions for θ″(q1,q3=0), rC−Cl ″ (q1,q2=0), and rC−I ″ (q1,q2=0) are shown in Figure 2. The parameters ai and bi in eqs 5 and 6 are used to obtain the correct scaling, whereas all parameters with a subscript “0” shift the functions to the equilibrium position of the respective internal coordinate. The above-described parametrization of the IRC path represents the typical motion of an SN2 backside attack including the umbrella motion of the CH3 group. To visualize the relation between the internal coordinates rC−Cl ″ , rC−I ″ , θ″ and the IRC-pad coordinates qi, the PES is depicted along rC−Cl and rC−I in Figure 3 together with the IRC-pad coordinates q1 and q2 (red arrows). The inset in Figure 3 shows the energy profile along the reaction path. In Figure 4 the PES is shown in dependence of the IRC-pad coordinates. The quantum dynamical calculations are performed with the coordinates q1, q2 (and q3 in the three-dimensional case). To carry out the wave packet calculations, the timedependent Schrödinger equation (TDSE) iℏ

∂ Ψ(q) = (T̂ + V̂ (q))Ψ(q) ∂t

Figure 3. Contour plot showing the relaxed potential energy surface with dependence on the halogen−carbon bond lengths. The red line depicts the parametrized reaction coordinate q1. For comparison, the energy profile along the reaction path is also shown in the inset, together with the corresponding chemical structures. The local minima are due to the anion dipole complexes, which is typical for the gas phase reaction.

Tq̂ ≃ −

(7)

ℏ2 2

∑∑ r

s

∂ ∂ Grs(q) ∂qr ∂qs

(8)

The first derivatives of the wave function with respect to qi are calculated by a fast Fourier transform. The G-Matrix elements Grs(q) combine the reduced mass with the coordinate transformation and are given by

is solved numerically on a spatial grid. The components of the vector q are the IRC-pad coordinates q1, q2, and q3. The time evolution is performed by a Chebychev propagation scheme.30 The PES V(q) is calculated with the same ab initio quantum chemistry methods as the IRC path. The kinetic operator T̂ is implemented by the G-Matrix method,31,32 which allows for flexible representation in any desired coordinate system:

3N

Grs(q) =

∑ i=1

4663

1 ∂qs ∂qr mi ∂xi ∂xi

(9)



Article

Ψ(q) = Φ0(q2) ϕ0(q3)

2 σr 2π

⎛ −(r″ (q ,q ) − r )2 ⎞ C − Cl 1 2 0 ⎟⎟e ikk 0q × exp⎜⎜ σ 2 r ⎝ ⎠

(10)

where Φ0 and ϕ0 are the uncoupled ground-state wave functions of the C−I bond and the umbrella mode obtained from the respective one-dimensional systems (note that in the asymptotic limit of separated reactants q2 and q3 correspond to rC−I and θ). The motion of the Cl− in the unbound direction is described by a Gaussian of width σr = a0 launched at a C−Cl distance r0 = 8 a0, where a0 is the Bohr radius. The specific collision energy is given through the momentum p = ℏkk0 characterized by the direction k0 and the magnitude k. For the two-dimensional model involving only q1 and q2 the product wave function in eq 10 is an exact eigenfunction to the Hamiltonian in the limit of separated reactants. In the threedimensional model involving q1, q2, and q3 the potential couplings between the C−I mode and the umbrella mode are neglected. A numerical test confirmed that this approximation has no significant influence on the results. The vector k0 has to contain the compensation for the kinetic couplings. For the given experimental conditions it is chosen to create exclusively a relative motion along q1. The direction of k0 can be achieved by solving the system of equations

Figure 4. PES represented in the IRC-pad coordinates q1 and q2 (q3 = 0). In this form the PES enters the quantum dynamics calculations. The red line and the red arrows indicate the linearized reaction path and the parallel shift as shown in Figure 3. The 2D wave packet dynamics have been carried out on this PES. The energy is limited to ≤5.44 eV.

where the sum over i involves all 3N degrees of freedom and xi are the Cartesian coordinates of the atoms in the center of mass (COM) frame. For practical reasons we first calculate the inverse of G because the derivatives ∂xi/∂qs can be evaluated from eqs 5 and 6 together with a Z-matrix describing the ″ , rC−Cl ″ , and θ″ transformation from the internal coordinates rC−I to the set of Cartesian coordinates {xi}. The Z-matrix format is the same as the one used in quantum chemistry programs to describe molecular geometries in terms of interatomic distances, interatomic angles, and dihedral angles (see ref 28, for example). To simulate the collision reaction in the IRC-pad coordinates, we also need to express the momentum in the initial wave packet accordingly. The challenge arises from the mixed derivatives in eq 8 inducing kinetic couplings among the IRC-pad coordinates. Thus, in the IRC-pad coordinates the direction of momentum and the direction of motion is sheared against each other. This fact has to be considered when momentum is added to the wave function at the beginning of a propagation. The initial conditions should be such that they represent the physical conditions of the collision experiment−a chloride approaching a methyl iodide in a well-defined vibrational state. For the head-on collision it is not enough to accelerate the chloride ion in the direction of q1 but it is necessary to correct for the virtual shearing forces caused by the mixed derivatives in the Hamilton, i.e., by the nondiagonal GMatrix. Those virtual forces arise from the transformation to the IRC-pad coordinate system. In a Cartesian coordinate system or in properly chosen skew coordinates33 these forces do not appear. The correction is implemented in the initial conditions in the following way. The relevant input parameter is the collision energy Ecoll of the colliding fragments Cl− and CH3I given in the COM frame. The physical constraint is that the umbrella mode and the C−I bond have to stay in a well-defined vibrational state during the approach as the intramolecular potential is flat (in the asymptotic region). This can be achieved by choosing an appropriate product wave function:

∂Tc(p) ∂q 0 = ∂t ∂p

with

k0 = 1

(11)

in the limit of separated reactants. Equation 11 represents the classical Hamilton equation of motion and Tc is the corresponding classical analogue of the kinetic energy T̂ q in the G-matrix representation from eq 8: Tc =

1 T p G(q1 = −∞,q2 = 0,q3 = 0)p 2

(12)

with the G-matrix G. The choice of the velocity vector (∂q0)/ (∂t) = (v1, 0, 0)T with v1 > 0 enforces the motion purely along the q1 axis. The procedure described above allows construction of an initial wave function that represents the given physical conditions within the IRC-pad coordinate model: a chloride anion approaching a methyl iodide with a defined relative collision energy Ecoll and its umbrella mode and C−I bond in a well-defined vibrational state. The wave function from eq 10 is an approximation to an eigenfunction of the reactants. Thus, the exact launch point r0 or the phase of the vibrational components of the wave packet do not affect the results. We then solve the TDSE from eq 7 numerically for different collision energies Ecoll to obtain the reaction probabilities and the energy redistribution in dependence of Ecoll. The quantities of interest are extracted from the time-evolved wave function. The separation into reactants and products is straightforward along the reaction coordinate q1. Moreover, using the ansatz from eq 10, it is also possible to separate the reaction product states into their vibrational and translational contributions.

4. RESULTS 4.1. Reaction Probabilities. In previous experiments a complex-mediated mechanism has been found to dominate for collision energies below 0.5 eV. This finding might be understood as a suppression of the direct rebound mechanism, 4664



Article

which has been observed9 for higher collision energies. It is therefore interesting to study the reaction probability for the direct rebound mechanism by means of wave packet simulations. 4.1.1. Two-Dimensional Model. We start with a twodimensional (2D) model that depends on the IRC-pad coordinates q1 and q2 and thus includes the two bond distances rC−Cl and rC−I as active coordinates. The umbrella angle θ is treated as a relaxed coordinate (q3 = 0) such that it minimizes the energy for a given set {q1, q2}. In this model it is possible to observe the direct collision dynamics. Processes which are due to rotation as well as the energy redistribution to internal degrees of freedom other than the C−I bond are not included. The latter processes typically lead to temporary trapping in the anion−dipole complex and an increase in the overall reaction probability but on a much longer time scale defined by the rotation of the complex. Furthermore, the geometry of the system is restricted to C3v neglecting all symmetry breaking modes. Thus, the presented reaction probabilities have to be considered as a relative measure for the short time dynamics. In Figure 5 the reaction

significantly increased. This effect cannot be explained by the additional energy deposited in the vibrational quanta, which is only 0.069 and 0.15 eV. Moreover, for higher collision energies an opposite behavior can be observed: the reaction probability is slightly suppressed above 0.9 eV. A similar effect was described for the system Cl− + CH3Br34 and will be discussed in section 5 in more detail. 4.1.2. Three-Dimensional Model: Influence of the Umbrella Mode. In a next step the umbrella mode, represented by the angle θ, is included as an active coordinate and the wave packet is propagated in the three-dimensional (3D) model. The reaction probabilities are compared for the 2D and 3D models in Figure 5 (black vs gray solid line). Apparently, the reaction probability is increased when the umbrella mode is included. The influence is strongest for values of Ecoll below 1 eV. Above 1.3 eV the reaction probabilities of the 2D and 3D models converge. Similar effects where the umbrella mode facilitates the reaction are described for the symmetric system Cl− + CH3Cl.5,35 On average the reaction probabilities are shifted by ≈−0.2 eV to smaller collision energies. Up to Ecoll = 0.75 eV it nearly matches the results from the full-dimensional classical dynamics simulations.13 The shift to lower collision energies can not be attributed to the difference in zero-point energy, which contributes only 0.09 eV. Direct excitation of the umbrella mode did not change the reaction probability. Thus, we can confirm that the low reduced mass of the CH3 inversion mode allows the group to adiabatically follow the product formation. This support is most sensitive for low collision energies; however, no fundamental difference in the mechanisms could be found in our simulations. 4.2. Internal Energy Distributions. Here we compare the ratios of internal and total available energy in the exit channel obtained in the experiment and in the quantum dynamics calculation. The available energy after the collision event is the sum of the collision energy Ecoll of the incoming fragments and an excess from the reaction enthalpy −ΔH. The energy transfer into rovibrational degrees of freedom is expressed by the fraction

Figure 5. Calculated reaction probabilities for the direct rebound mechanism versus collision energy. The black solid line shows the results from the 2D calculations. The comparison with the 3D model shows the influence of the umbrella mode (gray solid line), which slightly increases the reactivity (shifted by ≈−0.2 eV). The values from the reduced dimensional model have to be considered as relative measurements. The red diamonds indicate values from classical direct dynamics calculations for an impact parameter b = 0 (values taken from ref 13, enhanced by a factor of 30). A dramatic change in the probability happens if the C−I stretch mode is excited (dash blue line v = 1, dashed green line v = 2).

fint =

Ecoll − ΔH − ⟨E P⟩ Ecoll − ΔH

(13)

where ΔH = −0.53 eV and ⟨EP⟩ is the mean translational energy of the products. Thus, the numerator describes the energy deposited in the internal degrees of freedom, called the internal energy. In the simulation ⟨EP⟩ is obtained by analyzing the wave function after the reaction products have reached the asymptotic region in the exit channel. This wave function is then projected onto the vibrational and free particle wave functions analogue to eq 10 to qualify and quantify the energy stored in the vibrational degrees of freedom. In the 2D model only the C−Cl bond can store internal energy, whereas in the 3D model the energy can also be transferred into the umbrella mode of CH3Cl. In the experiment, the average internal energy and the fraction f int of internal versus total available energy have been determined for collision energies in the COM frame between 0.4 and 1.5 eV (see section 2). The results are shown in Figure 6 together with the calculated values. The latter ones are scaled by a factor of 2 to allow for a qualitative comparison, which also implies that the three reactive coordinates already allow us to reproduce approximately 50% of the total internal energy. The overall trend of the relative energy transfer is well

probabilities for different collision energies are shown. The result of the 2D model is shown as a black solid line. Even though there is no net barrier for the direct rebound channel a significant collision energy of about 0.55 eV is required to obtain a reaction probability of 50%. For lower energies the major part of the wave packet is reflected back. The velocity map images from ref 9 show a similar qualitative trend for a suppressed direct rebound mechanism at low energies. Fulldimensional classical trajectory calculations (red diamonds)13 shows the same trend up to 0.75 eV, besides a scaling factor. An explanation for this gradual increase in the reaction probability with increasing collision energy is discussed in section 5. Another interesting question that can be addressed in the quantum dynamical model is the influence of the initial vibrational state on the reaction outcome. For this purpose we replace the vibrational ground-state wave function Φ0 in eq 10 by a vibrational excited-state wave function of the C−I mode. The results are shown in Figure 5 (dotted and dashed-dotted line) for a vibrational quantum number of v = 1 and v = 2. For collision energies below 0.6 eV the reaction probability is now 4665



Article

most of the transferred collision energy is stored in the C−Cl stretching mode. Next we address the central question why the reaction probability increases only gradually with collision energy even though the TS for the direct rebound mechanism lies below the energy of the reactants. On the short time scale at first contact of the reaction partners the mechanism can not be solely described by the stationary points of the PES. The main finding for the low collision energy regime is that the energy transfer from the C−Cl to the C−I bond is quite inefficient. This is in agreement with results obtained by classical trajectory simulations.3 From the quantum dynamics simulations, we attribute that to the strong curvature in the reaction path (Figure 3). The accelerated wave packet does not follow the IRC path. Rather it is back-reflected and dynamical effects dominate the reaction. To classify these effects, the path of the wave packet is analyzed for different collision energies. The components of the wave function with a momentum toward the products are extracted and treated separately. The time-dependent expectation values of position ⟨q1⟩(t), and ⟨q2⟩(t) are then used to draw a simplified trajectory. An exemplary path is shown in Figure 7 for a collision energy of 0.84 eV (black line). The

Figure 6. Partitioning of the energy f int into the internal degrees of freedom of the product after the collision event. The calculated values from the 2D (black line) and the 3D (gray line) models are scaled by a factor of 2 to compare with the experimental values (blue diamonds).

reproduced; with increasing collision energy, a decreasing fraction of the total energy is partitioned into internal energy.

5. DISCUSSION In the previous sections we studied various aspects of the direct rebound mechanism such as the dependence on the collision energy and the influence of vibrational excitations. Only minor influence was found for the umbrella mode. The reaction probabilities obtained from the 2D and 3D models show that the inclusion of the umbrella bending coordinate facilitates the reaction also for low collision energies but does not change the reaction probabilities fundamentally. Less pronounced is the influence on the internal energy distribution. Moreover, an initial vibrational excitation in the umbrella mode with one quantum (0.16 eV) does not lead to a significant change in the reaction probability. From the wave packet calculations it can be concluded that the inversion of the CH3 group has the role of a spectator mode supporting the reaction adiabatically. These findings are in qualitative agreement with theoretical studies36,37 of the symmetric reaction Cl− + CH3Cl. Also in this case an excitation in the C−Cl bond leads to an increase in reactivity whereas the excitation of the umbrella mode shows no significant influence on the reaction probability. The investigation of the collision induced energy transfer into the vibrational modes of the products shows that with increasing collision energy the fraction of internal excitation decreases. This trend can be explained by the observation that the wave packet spends more time at the TS for low collision energies. The energy transfer to the internal degrees of freedom occurs via the potential in the region around the TS and thus depends on q1. In this region the Hamiltonian is no longer separable into vibrational and translational contributions, as it is in the asymptotic limit (eq 10), and collision energy can be transferred to vibrational energy. An incoming wave packet with lower collision energy passes the TS more slowly and thus the time for energy redistribution is longer. When the wave packet approaches the exit channel, the coupling between translational and vibrational energy disappears again and the newly formed CH3Cl remains in a fixed distribution of vibrational states. The experimental value at Ecoll = 0.39 eV in Figure 6 does not agree well with the theoretical trend. The discrepancy can be explained by the fact that the regime below 0.5 eV is dominated by the indirect, anion−dipole complex mediated mechanism, which is observed in the crossed-beam imaging data9 and not included in the theoretical model. From the comparison of the 2D and 3D models, we can conclude that

Figure 7. Contour plot of the PES highlights the region close to the TS. The gray dashed line indicates the IRC path. The black solid line depicts a trajectory evaluated for the expectation value of a wave packet accelerated by a collision energy of 0.84 eV. Shown in the inset are the corresponding energy profiles. They reveal that the actual path of the wave packet has to overcome an effective barrier due to the dynamics.

comparison with the IRC (dashed line) shows a large deviation especially in the region around the TS and discloses the reason for suppressed reactivity in the direct collision process. A large fraction of the wave packet is back-reflected at low collision energies without transfer of kinetic energy into the C−I bond. From this picture it also becomes obvious that the energy transfer is mediated by the shape of the dynamically reached potential around the TS. It determines whether the wave packet is back-reflected or deflected toward the products. The higher the initial momentum of the wave packet the less its motion follows the IRC path and the higher values on the PES are reached in the TS region. The inset in Figure 7 shows the effective potential under the trajectory (V(⟨q1⟩(t), ⟨q2⟩(t))) for a collision energy of 0.84 eV. Its maximum describes the 4666



Article

because with increasing quantum numbers the delocalization of the wave packet increases accordingly and back-reflection beyond the barrier also contributes. The observation of a strongly promoting effect of vibrational pre-excitation is consistent with the classical expectation for a pseudotriatomic reaction with a slightly late dynamic barrier.43 This is furthered by the fact that the attacking and the leaving ions are heavy as compared to the transferred CH3 group, which results in a strong bend of the PES in mass-skewed coordinates. The latter also leads to the expectation of a significant vibrational excitation of the product, in line with the presented results.

average potential barrier, which has to be overcome by the wave packet to form the products. In contrast to the IRC path the maximum of the dynamic path is well above the TS by 0.75 eV. From this picture a more intuitive concept of a dynamic barrier can be invoked. Similar concepts involving effective barriers are also mentioned in refs 38−41. We define the dynamic barrier as the potential energy maximum experienced along the path of a wave packet for a given collision energy: Edyn(Ecoll) = max{V(⟨q1⟩(t), ⟨q2⟩(t))}. In Figure 8 the height of

6. CONCLUSION We have shown wave packet calculations for a fundamental SN2 reaction in a reduced dimensionality model that resolves the direct rebound mechanism. An intuitive concept for the representation of reactive coordinates based on the reaction path has been introduced. The usage of the reaction path and its parallel shifts allows for an efficient description on a spatial grid covering specifically the potential energy range of interest. The presented IRC-pad coordinates may also be used for other reactions, which are well described by an unique IRC path. A 2D model involving the halogen carbon distances and a 3D model including the umbrella motion have been implemented. A comparison of the 2D and the 3D models demonstrates that the 2D model is sufficient to obtain a qualitative description of the direct collision reaction in the energy range 0.5−1.4 eV. Moreover, the CH3 umbrella motion has been identified to support the reaction. However, it is not essential for a qualitative description of the direct rebound mechanism. For the internal energy distribution of the products a comparison between theoretical calculations and velocity map imaging measurement has been made. The theoretical model is in good agreement with the amount of energy that is retained in the CH3Cl product. Only a scaling of the calculated data is needed to reproduce the experimental results. The wave packet simulations also lead to qualitatively new insight into the direct rebound mechanism of gas phase nucleophilic substitution reactions. The analysis of the wave packet path shows how the probability for the reaction is determined by the shape of the PES around the TS. This is described using the dynamical barrier concept, derived from the wave packet traces. This barrier is found to be substantially higher than the transition state and increases with collision energy. Thus, the reduction to a one-dimensional picture gives an intuitive explanation for the increased reactivity with collision energy. Overall, our quantum dynamics calculations show that the reaction is not sufficiently characterized by the stationary points on the PES but the regions at high energies become important due to the dynamical properties.

Figure 8. Height of the dynamical barrier Edyn versus the collision energy Ecoll. The black dots are the maximum potential energy values extracted from the wave packet simulations. The gray line is a linear fit to these points (Edyn ≈ 0.72Ecoll − 0.05 eV).

the dynamic barrier is plotted against the collision energy. The curve shows a linear behavior with a slope of 0.72, which means that the dynamic TS is ≈72% of the collision energy. The dynamical barrier can be thought of being associated with the dynamical reaction path and adapts the familiar concepts of stationary points on the PES to its dynamical variant. In contrast to the classical transition-state theory,42 which describes statistical processes, the dynamical barrier describes the nonstatistical properties of the reaction. Thus, it cannot be connected to a reaction rate. The concept of the dynamical barrier can also explain the observed change in reactivity when the C−I bond is excited (Figure 5). Unlike an excitation of the umbrella mode this leads to a significant change in the reaction probability. One or more quanta in the C−I stretch mode increase the spatial extend along this mode in the incoming wave packet. This spreading modifies the reflection behavior in the region around the TS. The nodal structure in the wave function of the vibrational excited C−I bond is of special importance to the reaction. In case of a one-quantum excitation (v = 1, Figure 5, blue dashed line) the local maxima of the probability density do not coincide with the equilibrium distance of the C−I bond. Instead one of the lobes of the wave packet resembles a shortened C−I bond and thus experiences a higher dynamical barrier. The other lobe of the wave packet corresponds to a slightly elongated C−I bond facilitating the reaction already at low collision energies. This geometrical interpretation provides an alternative explanation for the reactivity under vibrational excitation.The nearly constant reaction probability between 0.4 and 1.2 eV reflects the competition between the elongated C−I bond and the shortened C−I bond, where the latter one gains influence with increasing collision energy. Beyond 1.2 eV the reaction is dominated by the part of the wave packet corresponding to the shortened C−I bond, leading to an overall reduction of the reaction probability compared to a C−I vibration in its ground state. The concept of the dynamical barrier is still visible in the results for v = 2 (Figure 5, green dashed line). However, the analysis becomes more complex

■

ASSOCIATED CONTENT

S Supporting Information *

Table with parameters for eqs 5 and 6. Plot of the PES in coordinates q1 and q3. Plot of the G-Matrix elements along q1. This material is available free of charge via the Internet at http://pubs.acs.org/.

■

AUTHOR INFORMATION

Corresponding Author

*M. Kowalewski: e-mail, [email protected]. 4667



Article

Notes

(17) Gebhardt, C. R.; Rakitzis, T. P.; Samartzis, P. C.; Ladopoulos, V.; Kitsopoulos, T. N. Slice imaging: A new approach to ion imaging and velocity mapping. Rev. Sci. Instrum. 2001, 72, 3848−3853. (18) Li, W.; Huang, C. S.; Patel, M.; Wilson, D.; Suits, A. Stateresolved reactive scattering by slice imaging: A new view of the Cl +C2H6 reaction. J. Chem. Phys. 2006, 124, 011102. (19) Lin, J. J.; Zhou, J.; Shiu, W.; Liu, K. Application of time-sliced ion velocity imaging to crossed molecular beam experiments. Rev. Sci. Instrum. 2003, 74, 2495−2500. (20) Fukui, K. The path of chemical reactions - the IRC approach. Acc. Chem. Res. 1981, 14, 363−368. (21) Miller, W. H.; Handy, N. C.; Adams, J. E. Reaction path Hamiltonian for polyatomic molecules. J. Chem. Phys. 1980, 72, 99− 112. (22) Ruf, B. A.; Miller, W. H. A new (cartesian) reaction-path model for dynamics in polyatomic systems, with application to H-atom transfer in malonaldehyde. J. Chem. Soc., Faraday Trans. 2 1988, 84, 1523−1534. (23) Giese, K.; Kü hn, O. The all-Cartesian reaction plane Hamiltonian: Formulation and application to the H-atom transfer in tropolone. J. Chem. Phys. 2005, 123, 054315. (24) Kraka, E. Reaction path Hamiltonian and the unified reaction valley approach. Comp. Mol. Sci. 2011, 1, 531−556 and references therein.. (25) Bergner, A.; Dolg, M.; Küchle, W.; Stoll, H.; Preuß, H. Ab initio energy-adjusted pseudopotentials for elements of groups 13−17. Mol. Phys. 1993, 80, 1431−1441. (26) Dunning, T. H. Gaussian basis sets for use in correlated molecular calculations. I. The atoms boron through neon and hydrogen. J. Chem. Phys. 1989, 90, 1007−1023. (27) Woon, D. E.; Dunning, T. H. Gaussian basis sets for use in correlated molecular calculations. III. The atoms aluminum through argon. J. Chem. Phys. 1993, 98, 1358−1371. (28) Werner, H. J.; Knowles, P. J.; Lindh, R.; Manby, F. R.; Schütz, M.; Others, MOLPRO, version 2006.1, a package of ab initio programs. 2006. (29) Glukhovtsev, M. N.; Pross, A.; Radom, L. Gas-Phase NonIdentity SN2 Reactions of Halide Anions with Methyl Halides: A High-Level Computational Study. J. Am. Chem. Soc. 1996, 118, 6273− 6284. (30) Ezer, H. T.; Kosloff, R. An accurate and efficient scheme for propagating the time dependent Schrödinger equation. J. Chem. Phys. 1984, 81, 3967−3971. (31) Wilson, E. B.; Decius, J. C.; Cross, P. C. Molecular Vibrations: The Theory of Infrared and Raman Vibrational Spectra; Dover Publications: New York, 1980. (32) Stare, J.; Balint-Kurti, G. G. Fourier Grid Hamiltonian Method for Solving the Vibrational Schrodinger Equation in Internal Coordinates: Theory and Test Applications. J. Phys. Chem. A 2003, 107, 7204−7214. (33) Tannor, D. J. Introduction to Quantum Mechanics: A TimeDependent Perspective; University Science Books: Sausolito, CA, 2006. (34) Schmatz, S.; Clary, D. C. Quantum scattering calculations on the SN2 reaction Cl− + CH3Br → ClCH3 + Br−. J. Chem. Phys. 1999, 110, 9483−9491. (35) Clary, D. C.; Palma, J. Quantum dynamics of the Walden inversion reaction Cl− + CH3Cl → ClCH3 + Cl−. J. Chem. Phys. 1997, 106, 575−583. (36) Vande Linde, S. R.; Hase, W. L. A direct mechanism for SN2 nucleophilic substitution enhanced by mode-selective vibrational excitation. J. Am. Chem. Soc. 1989, 111, 2349−2351. (37) Vande Linde, S. R.; Hase, W. L. Trajectory studies of SN2 nucleophilic substitution. I. Dynamics of Cl−+CH3Cl reactive collisions. J. Chem. Phys. 1990, 93, 7962−7980. (38) Hase, W. L. Simulations of Gas-Phase Chemical Reactions: Applications to SN2 Nucleophilic Substitution. Science 1994, 266, 998−1002.

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS Financial support by the Deutsche Forschungsgemeinschaft through the excellence cluster Munich Centre for Advanced Photonics and the Sonderforschungsbereich 749 and the Centre for Interdisciplinary Mathematics, Uppsala, are gratefully acknowledged. M.K. and R.d.V. thank Regina Wyrwich for carrying out some of the QD calculations. We thank William L. Hase for valuable discussions. We dedicate this publication to the 75th birthday of Dr. Klaus Römer.

■

REFERENCES

(1) Laerdahl, J. K.; Uggerud, E. Gas phase nucleophilic substitution. Int. J. Mass. Spectrom. 2002, 214, 277−314 and references therein. (2) Mikosch, J.; Weidemüller, M.; Wester, R. On the dynamics of chemical reactions of negative ions. Int. Rev. Phys. Chem. 2010, 29, 589−617. (3) Manikandan, P.; Zhang, J.; Hase, W. L. Chemical Dynamics Simulations of X− + CH3Y → XCH3 + Y− Gas-Phase SN2 Nucleophilic Substitution Reactions. Nonstatistical Dynamics and Nontraditional Reaction Mechanisms. J. Phys. Chem. A 2012, 116, 3061−3080. (4) Craig, S. L.; Brauman, J. I. Perturbed equilibria and statistical energy redistribution in a gas-phase SN2 reaction. Science 1997, 276, 1536−1538. (5) Schmatz, S. Quantum Dynamics of Gas-Phase SN2 Reactions. Chem. Eur. J. Chem. Phys. 2004, 5, 600−617. (6) Hennig, C.; Schmatz, S. Four-dimensional quantum study on exothermic complex-forming reactions: Cl−+CH3Br → ClCH3+Br−. J. Chem. Phys. 2005, 122, 234307. (7) Viggiano, A. A.; Morris, R. A.; Paschkewitz, J. S.; Paulson, J. F. Kinetics of the gas-phase reactions of chloride anion, Cl− with CH3Br and CD3Br: experimental evidence for nonstatistical behavior? J. Am. Chem. Soc. 1992, 114, 10477−10482. (8) Angel, L. A.; Ervin, K. M. Dynamics of the Gas-Phase Reactions of Fluoride Ions with Chloromethane. J. Phys. Chem. A 2001, 105, 4042−4051. (9) Mikosch, J.; Trippel, S.; Eichhorn, C.; Otto, R.; Lourderaj, U.; Zhang, J. X.; Hase, W. L.; Weidemüller, M.; Wester, R. Imaging Nucleophilic Substitution Dynamics. Science 2008, 319, 183−186. (10) Zhang, J.; Mikosch, J.; Trippel, S.; Otto, R.; Weidemueller, M.; Wester, R.; Hase, W. L. F−+CH3I → FCH3+I− Reaction Dynamics. Nontraditional Atomistic Mechanisms and Formation of a HydrogenBonded Complex. J. Phys. Chem. Lett. 2010, 1, 2747−2752. (11) Mikosch, J.; Zhang, J.; Trippel, S.; Eichhorn, C.; Otto, R.; Sun, R.; DeJong, W.; Weidemüller, M.; Hase, W. L.; Wester, R. Indirect dynamics in a highly exoergic substitution reaction. J. Am. Chem. Soc. 2013, 135, 4250−4259. (12) Gronert, S.; DePuy, C. H.; Bierbaum, V. M. Deuterium isotope effects in gas-phase reactions of alkyl halides: distinguishing E2 and SN2 pathways. J. Am. Chem. Soc. 1991, 113, 4009−4010. (13) Zhang, J.; Lourderaj, U.; Sun, R.; Mikosch, J.; Wester, R.; Hase, W. L. Simulation studies of the Cl− + CH3I SN2 nucleophilic substitution reaction: Comparison with ion imaging experiments. J. Chem. Phys. 2013, 138, 114309. (14) Wester, R. Velocity map imaging of ion−molecule reactions. Phys. Chem. Chem. Phys. 2014, 16, 396−405. (15) Mikosch, J.; Trippel, S.; Otto, R.; Eichhorn, C.; Hlavenka, P.; Weidemüller, M.; Wester, R. Kinematically complete reaction dynamics of slow ions. J. Phys. Conf. Series 2007, 88, 012025. (16) Trippel, S.; Stei, M.; Otto, R.; Hlavenka, P.; Mikosch, J.; Eichhorn, C.; Lourderarj, U.; Zhang, J. X.; Hase, W. L.; Weidemüller, M.; Wester, R. Kinematically complete chemical reaction dynamics. J. Phys. Conf. Ser. 2009, 194, 012046. 4668



Article

(39) Lara, M.; Jambrina, P. G.; Varandas, A. J. C.; Launay, J. M.; Aoiz, F. J. On the role of dynamical barriers in barrierless reactions at low energies: S(1D) + H2. J. Chem. Phys. 2011, 135, 134313. (40) Pijper, E.; Kroes, G. J.; Olsen, R. A.; Baerends, E. J. The effect of corrugation on the quantum dynamics of dissociative and diffractive scattering of H2 from Pt(111). J. Chem. Phys. 2000, 113, 8300−8312. (41) Angel, L. A.; Ervin, K. M. Gas-Phase SN2 and Bromine Abstraction Reactions of Chloride Ion with Bromomethane: Reaction Cross Sections and Energy Disposal into Products. J. Am. Chem. Soc. 2003, 125, 1014−127. (42) Levine, R. D. Molecular Reaction Dynamics; Cambridge University Press: Cambridge, U.K., 2005. (43) Polanyi, J. C.; Wong, W. H. Location of Energy Barriers. I. Effect on the Dynamics of Reactions A+BC. J. Chem. Phys. 1969, 51, 1439−1450.

4669


Nucleophilic Substitution Dynamics: Comparing Wave Packet

Recommend Documents