- Email: [email protected]
Chapter 14 Quantum chemical topology and reactivity: A comparative static and dynamic study on a SN2 reaction
Elsevier AMS
Ch14-N52719
10-10-2006
8:04 a.m.
Page:287
Trimsize:165×240 MM
Theoretical Aspects of Chemical Reactivity A. Toro-Labbé (Editor) © 2007 Published by Elsevier B.V.
Chapter 14
Quantum chemical topology and reactivity: A comparative static and dynamic study on a SN 2 reaction Laurent Joubert, Ilaria Ciofini, and Carlo Adamo Laboratoire d’Electrochimie et de Chimie Analytique, UMR CNRS 7575, Ecole Nationale Supérieure de Chimie de Paris, 11, rue Pierre et Marie Curie, 75231 Paris Cedex 05, France
Abstract Ab initio molecular dynamic simulations, using density functional theory (DFT) and the recent atom-centered density-matrix propagation (ADMP) method, were employed to study the bond breaking and formation for a case-study SN 2 reaction. Using the real space partition scheme of Bader’s quantum chemical topology (QCT), we performed a population analysis to examine intra- and intermolecular electronic charge transfer along the ADMP trajectory. These results were compared to a static approach, which is performed along the intrinsic reaction coordinate (IRC) path. Although similar features are found for both static and dynamic approaches, the QCT analysis allows to rationalize the differences observed during the formation of the ion–molecule complex. In particular, the dynamic approach suggests a stronger electron exchange tending to spontaneously maximize both covalent and noncovalent (i.e. electrostatic) interactions.
1. Introduction The interpretation of quantum chemical calculations in terms of classical chemical concepts is not an easy task. In particular, electron transfer is strongly related to the definition of atoms and bonds in a molecular system. Beyond the simplistic picture of Lewis1 based on intuitive concepts, two different ways, trying to give quantitative information on the nature of chemical bond, cohabit. On the one hand, classical localization procedures of the wave function give different partitions of the Hilbert’s space 287
Elsevier AMS
Ch14-N52719
288
10-10-2006
8:04 a.m.
Page:288
Trimsize:165×240 MM
Quantum chemical topology and reactivity
leading to atomic charges, such as Mulliken or Löwdin ones,23 or to hybrid orbitals.4−6 More recently, topological analyses based on the electron density, such as Bader’s quantum chemical topology (QCT)78 or on derived functions, such as the electron localization function (ELF),9 have proven their reliability as robust methods, crafted to put in evidence the subtle properties of bonding and reaction mechanisms.10−13 Nowadays, a wide choice of tools are available to chemists for the theoretical analysis of any kind of chemical reactions. These tools are even more powerful when coupled with efficient algorithms for the exploration of the potential energy surface, such as the intrinsic reaction coordinate (IRC) algorithm.1415 Calculation of any topological or electronic variable along the reaction path allows for a clear-cut vision of the reaction mechanism and for rationalizing the structure/energies properties of the surface extrema (reactants, transition states, and products).1617 It must be noticed, anyway, that all these approaches have been developed in the framework of a time-independent (i.e. static) view. Only recently Gross and co-workers have formulated a suitable, but rather complex approach including the time evolution of ELF.18 At the same time, ab intio dynamics, such as the Carr–Parrinello approach,19 have recently reached a mature state, and a blooming of reactivity studies can be found in the literature.20−22 Based on an extended Lagrangian molecular dynamics (MD) scheme, this Born–Oppenheimer approach evaluates the potential energy surface at the density functional theory (DFT) level, and both the electronic and nuclear degrees of freedom are propagated as dynamical variables.23 The calculation of the time variation of any chemicophysical properties (observable or not) is then easily computed using structure snapshots extracted from the trajectory and a statistical average.23 A tentative to explore the difference between static and dynamics ab initio approaches has been done,24 but is has been focused only on the thermochemistry, since the choice of unique and meaningful parameters for comparison is not trivial. In the present article, we report a study concerning the reaction mechanism of a prototype reaction using both static and dynamic approaches to explore a DFT potential surface. The static approach is the standard IRC model, while the dynamic one is based on a Carr–Parrinello method performed with localized (Gaussian) orbitals, the so-called atom-centered density matrix propagation (ADMP) model.25 Our aim is to elucidate the differences, and the common aspects, between the two approaches in the analysis of bond breaking/formation. To this end, we have chosen topological quantities as probe molecular descriptors. The so-called Walden inversion, an SN 2 reaction shown in Figure 1, has been chosen as the model, since it is a well-studied reaction, from both a static and a dynamical point of view.26−36 In the present study, we limit our analysis to the gas-phase reaction, and we do not consider solvent effects that are known to strongly modify the potential energy surface (PES).29
2. Computational details All the calculations were carried out with the Gaussian03 program,37 using the hybrid PBE0 functional3038 and the 6-31+G(d,p) basis set.39 Starting from the transition state, we have calculated reaction pathways for both static and dynamic approaches. On the one hand, we followed the minimum energy path between the transition state and
Elsevier AMS
Ch14-N52719
10-10-2006
C. Adamo et al.
8:04 a.m.
Page:289
Trimsize:165×240 MM
289
the stable ion–molecule complex that is, using mass-weighted coordinates, the internal reaction coordinate (IRC).15 On the other hand, the dynamic simulation was performed at 298 K in the canonical ensemble, for a total simulation time of 1.5 ps. The starting point (transition state) was previously optimized with the same basis set and functional used in the simulation. The well-known Velvet algorithm40 was employed for the integration of the equations of motions using a time step of 0.25 fs. The fictitious mass of the electron was set to 0.20 amu, with a scaling for both core and valence electrons as described in Schlegel et al.25 The velocities of the nuclei were scaled each five time steps to ensure a constant temperature within a T = 5 K of tolerance. The stability of the simulations was monitored by checking at each step the idempotency of the density matrix (within a 10−12 threshold) and the so-called adiabaticity index (within a 10−4 threshold, see Shlegel et al. (2001)25 for more details). In practice, we restricted our trajectory analysis to the time frame of interest, scrutinizing the path only during the formation of the stable ion–molecule complex. All along this truncated trajectory (less than 100 fs in time scale), we extracted a number of snapshots, taken each 2 fs, for the subsequent structural and electronic analyses of the reaction. The charge transfer processes were investigated using a QCT approach: the “Atoms in Molecule” (AIM) partition scheme of the electron charge density.78 According to this theory, topological atoms are defined as regions in real space consisting of a bundle of electron density gradient paths attracted to a nucleus. This partition allows evaluating atomic properties, defined as volume integrals over non-overlapping atomic basins. In particular, the population associated with an atom is simply the volume integral of the electronic charge density over the basin. In this work, we calculated and compared the variations of atomic basin populations all along the selected static and dynamic paths. Moreover, we monitored the variations of typical electron density properties calculated at a bond critical point (BCP) to characterize a chemical bond (see for instance Popelier (2000)8 for more details on BCPs). These calculations were performed using the Topmod package41 and a locally developed code.
3. Results The prototypical SN 2 reaction studied in this work (see figure 1) presents in the gas phase a double-well PES with two equivalent local minima corresponding to the formation of a pre- and a post-reaction ion–molecule complex Cl− · · · CH3 Cl and a transition state (TS) of D3h symmetry Cl · · · CH3 · · · Cl− . Owing to the overall symmetry of the reaction, we will restrain our study to a reaction path from the transition state to one of the equivalent minima.
3.1. Tuning the DFT approach The present study relies on the parameter-free exchange-correlation PBE0 functional. Our choice of the functional was guided by its good performance in predicting thermodynamics data for the current SN 2 reaction, using a large basis set (6-311+(3df,3pd)30 ). Since a large number of calculations are needed in order to evaluate the variations of the topological quantities along the trajectories, such large basis is not suitable. For
Elsevier AMS
Ch14-N52719
10-10-2006
8:04 a.m.
290
Page:290
Trimsize:165×240 MM
Quantum chemical topology and reactivity
[Cl…CH3…Cl]– ΔEovr ClCH3+ Cl–
Energy
Cl– + CH3Cl
ΔE#
ΔEcomp Cl– … CH3Cl
ClCH3…Cl–
Reaction coordinate
Figure 1 Sketch of energy profile for the SN 2 reaction under study
this reason, we performed some preliminary tests to assess the validity of the PBE0 functional with a medium-sized basis set, 6-31+G(d,p). We first examined the geometrical features of the key structures, i.e. the geometries of the ion–molecule complex and of the transition state. Structural and energetic results are summarized in Tables 1 and 2, respectively. From a structural point of view, our data are very close to the MP2 results32 for the transition state structure, with a deviation of about 0.01 Å for the carbon–chlorine distance. The same trends are observed for the C-Cl2 bonding distance in the ion–molecule complex. In contrast, a larger error (0.20 Å) is observed for the C Cl1 distance, this effect being strictly related to the small basis set considered. In fact, when a larger basis, 6-311+G(2d,p), is considered this distance is significantly augmented (3.23 Å) toward the MP2 value. The most significant thermodynamic quantities are the complexation energy of the ion–molecule complex Ecomp , the activation energy, i.e. the relative energy of the D3 h saddle point with respect to the ion–molecule complex E # , and the overall barrier Eovr , defined as the difference between these two energies. These data are reported in Table 2 and compared with previous theoretical29−35 or experimental results.36 Whereas the computation of initial closed-shell reagents does not generate particular difficulties, the determination of the energy of the charged transition state energy Cl · · · CH3 · · · Cl− by DFT approaches is more involved. In fact, the majority of standard functionals, and in particular those resting on the Generalized Gradient Approximation (GGA),
Table 1 Computed geometry parameters (in Å) for the equilibrium ion–dipole complex and the transition state Complex
HFa /TZ3P + R + 2fd MP2a /TZ3P + R + 2fd DFT-LDAa /TZ3P + R + 2fd DFTb -BP/TZ + 2P DFTc -BP/PAW DFTd -mPW1PW/6-31 + Gdp DFT-PBE0/6-31+G(d,p)
Transition state
RC-Cl1
RC-Cl2
RC-C11 = RC-C12
3 37 3 27 2 98 3 10 3 15 3 01 3 07
1 82 1 81 1 81 1 84 1 89 1 83 1 83
2 39 2 32 2 28 2 34 2 37 2 33 2 33
a: Ref [31], b: Ref [32], c: Ref [34], d: Ref. [29].
Elsevier AMS
Ch14-N52719
10-10-2006
8:04 a.m.
Page:291
C. Adamo et al.
Trimsize:165×240 MM
291
Table 2 Complexation energies, overall energies, and barrier heights (in kcal/mol) for the studied SN 2 reaction Ecomp HFa /TZ3P + R + 2fd; ZPE corrected MP2a /TZ3P + R + 2fd; ZPE corrected MP4b TZ + 2P; ZPE corrected G2c ; ZPE corrected G3c ; ZPE corrected CBS-QB3(+)c ; ZPE corrected DFTb -BP/TZ + 2P; ZPE corrected DFTd -BP/PAW DFTe -B3LYP/6-31G(d); ZPE corrected DFTf -mPW1PW/6-31+G(d,p); ZPE corrected DFT-PBE0/6-31+G(d,p); ZPE corrected experimentg
−8 1 −10 5 −10 6 −10 8 −11 2 −10 7 −10 3 −8 3 −9 5 −9 8 −10 1 −12 ± 2
Eovr
E #
7 6 3 5 1 8 3 1 1 8 2 4 −5 7 −3 0 −0 9 0 7 0 5 3/1 ± 1
15 7 14 0 12 4 13 8 13 0 13 1 4 6 5 4 8 7 10 5 10 6 13 ± 2
a: Ref [31], b: Ref [32], c: Ref. [33], d: Ref [34], e: Ref [35], f: Ref. [29], g: Ref. [36].
fail in determining the energy barriers. This effect has been related to the the selfinteraction error in the exchange part as well as in the correlation part, which implies an exaggerated delocalization of the electron density that overstabilize the transition state (for a discussion on this point see Toulouse et al. (2002)42 ). In contrast, PBE0 calculations provide accurate results, close to the post-HF (Hartree-Fock) calculations and in the range of the experimental estimations. Nevertheless, this result can be due to a simple error compensation between the GGA and HF contributions.
3.2. Static and dynamic topological analyses The energy variations along the intrinsic reaction path (IRP) and the ADMP trajectory are reported in Figures 2 and 3, respectively. We want to underline that the IRP corresponds to a unique and constrained minimum energy path and does not contain any time information, whereas the ADMP profile is computed along a trajectory issued by a dynamics simulation. Therefore, any direct (i.e. point-to-point) comparison is misleading. Nevertheless, their global differences can be discussed. The minimum energy constraint imposes a unique and distinctive IRP that continuously decreases from the transition state to the minimum. In contrast, the energy globally decreases along the ADMP trajectory, but it appears clearly that the reaction does not follow this minimum energy path. Actually, the minimum is rapidly reached after 96 fs, but the curve exhibits some energy risings corresponding to punctual destabilizations of the whole system. In order to rationalize the differences observed between the static approach and dynamic process, we investigated the intra- and intermolecular charge transfers along the different reaction paths. As a first step, we looked at the AIM atomic populations for the two states of interest, i.e. the transition state and ion–molecule complex. The results are gathered in Table 3. At the transition state, each chlorine atom bears a negative charge of −0 70e− . This electron excess corresponds to the sharing of the unit
Elsevier AMS
Ch14-N52719
10-10-2006
292
8:04 a.m.
Page:292
Trimsize:165×240 MM
Quantum chemical topology and reactivity
Potential energy (a.u.)
–959.990
–959.995
–960.000
–960.005
–960.010 10
0
20
30
50
40 1/2
MWC step (amu
60
.bohr)
Figure 2 Variations of the potential energy along the IRP from the transition state to the ion–molecule complex
–959.985
Potential energy (a.u.)
–959.990
–959.995
–960.000
–960.005
–960.010 0
12
24
36
48
60
72
84
96
Time (fs)
Figure 3 Variations of the potential energy along an ADMP trajectory from the transition state to ion–molecule complex
charge of the complex increased by a supplementary electron transfer of approximately 0 20 e− from the hydrogens to each of the chlorine atom. Along the IRP, the formation of the ion–dipole molecule corresponds to a global charge transfer of 0 25 e− from the forming chloromethane moiety to the leaving chlorine atom. This charge transfer is slightly increased by 0 01 e− when considering the dynamic process. This global molecule–ion interaction is accompanied by an intramolecular charge transfer in the forming chloromethane molecule. This phenomenon is evidenced by a slight increase of the methyl moiety electronic population, i.e. 0.05 and 0 06 e− at the end of the static
Elsevier AMS
Ch14-N52719
10-10-2006
8:04 a.m.
Page:293
C. Adamo et al.
Trimsize:165×240 MM
293 Table 3 Computed AIM basin populations e− for the transition state and the ion–molecule complex Atom
Transition state
Cl1 H1 H2 H3 C Cl2
Ion–molecule complex
17 70 0 83 0 83 0 83 6 11 17 70
17 95 0 84 0 88 0 86 6 10 17 38
and dynamic process, respectively. A deeper insight into the charge transfer processes can be reached by examining the variations of the atomic electronic populations along the IRP and ADMP trajectory. Figures 4 and 5 represent the relative variations of selected AIM basin populations with respect to the initial populations in the transition state configuration for both the static and dynamic approaches, respectively. The global molecule-to-ion charge transfer is materialized by two essential curves. The first one (squares) corresponds to the monotonically increasing atomic population of the leaving chlorine atom Cl1, tending to the population of a chloride anion. In parallel, a second curve (solid), grouping the atomic populations of the entire molecule CH3 Cl, shows the corresponding and decreasing electronic population of the forming molecule. When the ion–dipole complex is formed, the global electronic transfer between the two moieties corresponds to the aforementioned net charge of 0 25e− or 0 26e− , depending on the approach chosen to study the reactivity of the system, i.e. a static or a dynamic
Variations of AIM basin populations (e–)
0.4 0.3 0.2 0.1 0.0 –0.1 –0.2 –0.3 –0.4 0
10
20
30
40
50
60
MWC step (amu1/2.bohr)
Figure 4 Variations of AIM populations along the IRP. The superbasin corresponding to the forming CH3 Cl molecule (5 basins) is represented by the solid curve. Other curves with markers: carbon basin (lozenges), Cl1 basin (squares), Cl2 basin (circles), and sum of hydrogen basins (triangles) populations
Elsevier AMS
Ch14-N52719
10-10-2006
8:04 a.m.
294
Page:294
Trimsize:165×240 MM
Quantum chemical topology and reactivity
Variations of AIM basin populations (e–)
0.4 0.3 0.2 0.1 0.0 –0.1 –0.2 –0.3 –0.4 0
12
24
36
48
60
72
84
96
Time (fs)
Figure 5 Variations of AIM populations along the ADMP trajectory. All symbols are defined in Figure 4
one. If we now examine the variations of the individual atomic populations, substantial differences are observed between the static and dynamic models of reactivity. Three supplementary curves are presented on the abovementioned Figures 4 and 5. They correspond to the variations of the atomic populations of (i) the chlorine atom Cl1 of the forming chloromethane molecule, (ii) the three gathered hydrogen atoms, and (iii) the carbon atom. In Figure 4, i.e. along the IRP, the detailed analysis of these variations of the atomic populations allows us to propose a three-step charge transfer mechanism. In the first part of the path, corresponding approximately to 1 0 amu1/2 bohr, the electron excess on the chlorine atom Cl2 is transferred almost totally to the other chlorine through the atomic basins of the hydrogens. This particular point is evidenced by examining the variations of the atomic populations of the two chlorine atoms that almost compensate each other. In other words, the electronic flux coming from the Cl2 chlorine basin and entering the three hydrogen basins corresponds to the exiting flux that penetrates the Cl1 chlorine basin. Besides this main electron transfer, a second one is observed between the carbon and the three hydrogen atoms. In fact, along the path, the carbon–hydrogen bond lengths systematically increase when the methyl group adopts a C3v geometry in the forming molecule and points to the leaving chlorine atom. This phenomenon is here increased by a favorable increasing electrostatic interaction with this chlorine atom. Therefore, the lengthening of these bond lengths explains the local charge transfer observed from carbon to hydrogens. Finally, we note substantial fluctuations in the population variations of the linked carbon and chlorine atoms. Again, these fluctuations compensate each other. They result from two concerted effects. First, a small fraction of the electron flux between chlorine atoms penetrates the carbon basin, resulting in a small increase of the corresponding population that tends to be slightly in excess. Second, the formation of the carbon–chlorine bond counters this flux by inducing a small charge transfer from the carbon to the chlorine atom. In the second step of the charge transfer mechanism, i.e. up to 3 0 amu1/2 bohr, the electronic loss on the chlorine atom Cl2
Elsevier AMS
C. Adamo et al.
Ch14-N52719
10-10-2006
8:04 a.m.
Page:295
Trimsize:165×240 MM
295
becomes stronger that the electronic gain on the leaving chlorine Cl1. Therefore, the electron flux that enters the hydrogen basins becomes larger than the exiting flux. As a result, the electronic population in the hydrogen basins increases. At a lesser extent, the electronic charge transfer from carbon to hydrogens increases as well, reinforcing the electronic population of these basins. Furthermore, we note that the fluctuations observed along the carbon and chlorine (Cl2) curves are amplified when approaching the equilibrium structure of the chloromethane. The last step of the charge transfer mechanism corresponds approximately to the second half of the IRP. Structural changes are minimal in the chloromethane molecule, resulting in equilibrium for carbon and hydrogen basin populations. The chlorine atom Cl1 continues to move away from the molecule implying a strong decrease in the chlorine to chlorine charge transfer, the population of the Cl1 atom approaching those of a chloride anion. In Figure 5, i.e. along the ADMP trajectory, the analysis of atomic population variations suggests a charge transfer mechanism at variance with the static IRC mechanism. We can decompose it in five steps. We first note that the system is destabilized during the first 8 fs of the simulation (see Figure 3a), certainly due to an excess of initial kinetic energy. During that period, population variations are negligible, except for the cyclic fluctuations of the carbon and hydrogen populations corresponding to the C-H stretching vibrational mode. These fluctuations are present all along the trajectory. Between 8 and 16 fs, the charge transfer begins between chlorine atoms. This second dynamic step is very similar to the first part of the IRC path, where an electron excess on the chlorine atom Cl2 is transferred almost totally to the other chlorine through the atomic basins of the hydrogens. The third step corresponds roughly to a frame between 16 and 30 fs and can be compared to the second step of the IRC charge transfer mechanism. In fact, the curve corresponding to the population variations of the Cl2 chlorine atom strongly deviates from the one corresponding to the population variations of the whole chloromethane molecule. Meanwhile, both electronic populations of carbon and hydrogen atoms substantially increase, indicating an intramolecular charge transfer from the chlorine Cl2 basin to the carbon and hydrogen basins. In contrast with the static approach, no fluctuations are observed for the carbon and chlorine Cl2 curves, except for the cyclic variations on the carbon curve corresponding to the C-H stretching vibration mode. At this point, we want to remark that all the obtained results are preserved when different initial conditions (kinetic energy and temperature) are considered, since the excess energy mainly changes the population of the CH vibrational states. This is not surprising since the Cl-CH3 Cl starting ion–complex exhibits a rather poor energy transfer between the inter- and intravibrational modes.43 The third step corresponds to a time frame between 30 and 50 fs. We observe that the electronic population of the hydrogen basins strongly increases due to a substantial electron transfer from the carbon basin. Moreover, the corresponding population variations almost compensate each other. All these subtle variations are related to structural changes during the dynamics. In Figure 6, we isolated the population variations of the carbon and hydrogen basins together with the sum of the variations of C-H distances. During the time frame of interest, we note a substantial lengthening of these distances. This phenomenon was also observed along the static IRP, but with a lower amplitude. Here, the C-H stretching vibrational mode elongates the corresponding bond lengths, strengthening the electrostatic attraction between the leaving chlorine atom and hydrogen atoms and favoring the electron flux
Elsevier AMS
Ch14-N52719
10-10-2006
8:04 a.m.
Trimsize:165×240 MM
Quantum chemical topology and reactivity 0.4
0.08
0.3
0.06
0.2
0.04
0.1
0.02 0
0
–0.1
–0.02
–0.2
–0.04
–0.3
–0.06
–0.4
0
12
24
36
48
60
72
84
Variations of the sum of C-H distances (Å)
Variations of AIM basin populations (e–)
296
Page:296
–0.08 96
Time (fs)
Figure 6 Comparison between the variations of selected basin populations (C and H) and the sum of C-H distances (solid line) along a dynamic pathway from the transition state to the ion–molecule complex. Curves with markers: carbon basin (lozenges) and sum of hydrogen basins (triangles) populations
between them. This noncovalent effect should favor the stabilization of the system, but we know that the complex is strongly destabilized within this time frame, as evidenced by a sudden potential energy rising in Figure 3. This surprising behavior can also be rationalized by examining the structural evolution of the system. In Figure 7, we represented the variations of the potential energy together with the evolution of the H-C-H valence angles, i.e. we plotted the variations of the sum of these three angles along the trajectory. In the time frame of interest, we can see that the sum of angles decreases
–959.985
370
Potential energy (a.u.)
350 –959.995 340 –960.000 330 –960.005
–960.010
320
0
12
24
36
48
60
72
84
Sum of H-C-H valence angles (°)
360
–959.990
310 96
Time (fs)
Figure 7 Comparison between the potential energy variations (in bold) and the sum of the three H-C-H valence angles (solid) in the ion–molecule complex
Elsevier AMS
Ch14-N52719
10-10-2006
8:04 a.m.
Page:297
C. Adamo et al.
297 0.5
2.4 2.3
0.4
2.2
0.3
2.1
0.2
2.0 1.9
0.1
1.8
0
1.7 1.6
Rho(rc) et L(rc) (a.u.)
C–Cl2 distance (Å)
Trimsize:165×240 MM
0
12
24
36
48
60
72
84
–0.1 96
Time (fs)
Figure 8 Left Y axis: variations along the ADMP trajectory of the C-Cl2 distance (in bold). Right Y axis: variations of the electron density (solid) and the function Lr (dashed) calculated at the C-Cl2 bond critical point r = rc
to only 310 , a global minimum along the whole trajectory. In other words, the strong electrostatic attraction between the Cl1 chlorine atom and hydrogen atoms induces not only a lengthening of the C-H bond lengths but also a substantial closing of the H-C-H valence angles. The latter effect may explain partly the destabilization of the system by increasing electrostatic and Pauli repulsions between the hydrogen atoms. Furthermore, an examination of electron density properties at the C-Cl2 BCP reveals interesting information on the bond creation process that helps to explain this surprising potential energy rising. In Figure 8, we plotted the variations of the C-Cl2 distance along the trajectory, together with the variations of the density r and function Lr = − r calculated at the corresponding BCP. In the time frame of interest, we note a systematic decrease of the C-Cl2 distance down to the global minimum value of 1.682 Å at = 50fs, about 0.15 Å shorter than the equilibrium value. Meanwhile, the BCP electron density strongly increases to finally reach a value of 0 22 e− bohr −3 , characteristic of a strong covalent bond. The Lr function confirms these trends, with an increasing positive value at BCP that indicates a strong charge concentration. In short, the electron exchange is maximized during the dynamic simulation, favoring the penetration of the chlorine atom electron cloud into the one of the methyl moiety. This effect destabilizes the system by strongly enhancing the Pauli repulsion between these atoms. Finally, the two last steps of the dynamic charge transfer mechanism correspond to the second half of the trajectory. First, between 50 and 78 fs, the system, strongly destabilized, tends to reduce quickly its potential energy. Thus, we note a sudden lengthening of the C-Cl2 bond and an opening of the H-C-H valence angles (see Figure 7) coupled with a shortening of the C-H bond lengths (see Figure 6). This last effect induces a strong diminution of the charge density transferred from the carbon to the hydrogen atoms. Moreover, it is important to point out an inverted electron flux from the ion to the chloromethane molecule tempting to equilibrate the destabilized system. Finally, the last part of the dynamic trajectory, beyond 78 fs, corresponds to the effective formation of a
Elsevier AMS
Ch14-N52719
298
10-10-2006
8:04 a.m.
Page:298
Trimsize:165×240 MM
Quantum chemical topology and reactivity
stable ion–molecule complex and the end of the charge transfer between the molecule and the chlorine atom stabilized by the completion of its external valence shell.
4. Some comments on the static vs. dynamic descriptions Our results show how both ADMP and IRP approaches suggest that the driving force for bonding formation is the charge transfer process. This force, clearly present in the static approach, is so enhanced in the dynamics simulation that it allows for the penetration of the incoming chlorine into electron cloud of the methyl moiety. This results in a very short C-Cl2 distance and in high energy. The subsequent switching on of the Pauli repulsion permits to reach the final state, through a large C-Cl2 amplitude motion. Before this reactive state, the reactants are in a kind of “harmonic reversible” regime where each part still preserves its own electronic identity. This picture is consistent with previous analyses, based on different theoretical models. Among others, we want to recall that similar but maybe more qualitative conclusions were drawn several years ago by Shaik on the basis of molecular orbitals arguments resting on valence bond (VB) calculations.44 In the same philosophy, but more recently, it has been suggested that, within a VB approach, the overlap between the active orbitals of the incoming Cl− and CH3 Cl rules the overall reactive process.43 Moreover, the maximization of the charge transfer strongly reminds the “maximum localization hybrid orbitals” overlap principle suggested more than 30 years ago by Del Re. This criterion has been used as a driving force in a chemical reactivity study based on semi-empirical qualitative approaches.45 Analogous results were also reached by Toro–Labbé using the principle of maximum hardness.46 What is new in our analysis is that the dynamics simulations, based on a ab initio approach using Gaussian orbitals, well underlines this charge transfer phenomenon and its amplitude as driving force. At the same time, we find in a modern and accurate simulation qualitative descriptions and working hypothesis carried out in more approximate schemes. Reassuringly, analyses based on different theoretical backgrounds will converge toward the same description of the same phenomenon, when correctly applied.
5. Conclusion We have analyzed the mechanism of bond breaking and formation in a prototype SN 2 reaction using a static and dynamic approach, based on both DFT and localized (Gaussian) basis sets. Using the real space partition scheme of Bader’s QCT, we performed a population analysis to examine intra- and intermolecular electronic charge transfers along the dynamics trajectory and the static reaction profile. QCT analysis allows the use of well-defined molecular descriptors for the same overall phenomenon (the reaction) in two different frameworks. Our results show that similar features are present in the static and dynamic approaches, but some differences can be observed during the formation of the ion–molecule complex. In particular, the dynamic approach suggests a stronger electron exchange between the different moieties tending to spontaneously maximize both covalent and noncovalent interactions.
Elsevier AMS
Ch14-N52719
C. Adamo et al.
10-10-2006
8:04 a.m.
Page:299
Trimsize:165×240 MM
299
Acknowledgments The authors thank Prof. V. Barone (Naples, Italie) for fruitful discussions and Dr. Michele Pavone (Princeton, USA) for helpful suggestions in the ADMP simulations.
References 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11. 12. 13. 14. 15. 16. 17. 18. 19. 20. 21. 22. 23. 24. 25.
26. 27. 28. 29. 30. 31. 32. 33. 34.
G. N. Lewis, J. Am. Chem. Soc. 38 (1916) 762. R. S. Mulliken, J. Chem. Phys. 23 (1955) 1833. P. O. Löwdin, J. Chem. Phys. 18 (1950) 365. G. Del Re, Theor. Chim. Acta 1 (1963) 188. J. Pipek, P. G. Mezey, J. Chem. Phys. 90 (1989) 4916. A. E. Reed, L. A. Curtiss, F. Weinhold, Chem. Rev. 88 (1988) 899–926. R. F. W. Bader, Atoms in Molecules: A Quantum Theory, Oxford University Press, Oxford, 1990. P. L. A. Popelier, Atoms in Molecules: An Introduction, Pearson Education, London, 2000. A. D. Becke, K. E. Edgecombe, J. Chem. Phys. 92 (1990) 5397. B. Silvi, A. Savin, Nature 371 (1994) 683. S. Berski, J. Andrés, B. Silvi, L. R. Domingo, J. Phys. Chem. A 107 (2003) 6014. P. Geerlings, F. De Proft, W. Langenaeker, Chem. Rev. 103 (2003) 1793. P. L. A. Popelier, L. Joubert, J. Am. Chem. Soc. 124 (2002) 8725. K. Fukui, Acc. Chem. Res. 14 (1981) 363. C. Gonzalez, H. B. Schlegel, J. Chem. Phys. 90 (1989) 2154. M. Solà, A. Toro-Labbé, J. Phys. Chem. A 103 (1999) 8847. V. Polo, J. Andrés, J. Comput. Chem. 26 (2005) 1427. T. Burns, M. A. L. Marques, E. K. U. Gross, Phys. Rev. A 71 (2005) 10501. R. Car, M. Parrinello, Phys. Rev. Lett. 55 (1985) 2471. T. Ziegler, J. Autschbach, Chem. Rev. 105 (2005) 2695. D. Gleich, J. Hutter Chem. Eur. J. 10 (2004) 2435–2444. Y. Tateyama, J. Blumberger, M. Sprik, I. Tavernelli, J. Chem. Phys. 122 (2005) 234505. D. Marx, J. Hutter, in “Modern Methods and Algorithms of Quantum Chemistry Proceedings”, J. Grotendorst (Ed.), Forschungszentrum Jülich, NIC Series, Vol. 3, page 339 (2000). S. C. Ammal, H. Yamataka, M. Aida, M. Dupuis, Science 299 (2003) 1555. H. B. Shlegel, J. M. Millam, S. S. Iyengar, G. A. Voth, A. D. Daniels, G. E. Scuseria, M. J. Frish, J. Chem. Phys. 114 (2001) 9758; S. S. Iyengar, H. B. Shlegel, J. M. Millam, G. A. Voth, G. E. Scuseria, M. J. Frish, J. Chem. Phys. 115 (2001) 10291; H. B. Shlegel, S. S. Iyengar, X. Li, J. M. Millam, G. A. Voth, G. E. Scuseria, M. J. Frish, J. Chem. Phys. 117 (2002) 8694. W. P. Ho, D. G. Truhlar, J. Am. Chem. Soc. 117 (1995) 10726. J. M. Gonzales, R. S. Cox III, S. T. Brown, W. D. Allen, H. F. Schaefer III, J. Phys. Chem. A 105 (2001) 11327. B. Ensing, E. J. Meijer, P. E. Blöchl, and E. J. Baerends, J. Phys. Chem. A 105 (2001) 3300. M. Cossi, C. Adamo, V. Barone, Chem. Phys. Lett. 297 (1998) 1. C. Adamo, V. Barone, J. Chem. Phys. 110 (1999) 6158. B. D. Wladkowski, K. F. Lim, W. D. Allen, J. I. Brauman, J. Am. Chem. Soc. 114 (1992) 9136 and references therein. L. Deng, V. Branchadell, T. Ziegler, J. Am. Chem. Soc. 116 (1994) 10645. S. Parthiban, G. D. Oliveira, J. M. L. Martin, J. Phys. Chem. A 105 (2001) 895. S.-Y. Yang, P. Fleurat-Lessard, I. Hristov, T. Ziegler, J. Phys. Chem A 108 (2004) 9461.
Elsevier AMS
300
Ch14-N52719
10-10-2006
8:04 a.m.
Page:300
Trimsize:165×240 MM
Quantum chemical topology and reactivity
35. A. Streitwieser, G. S. Choy, F. J. Abu-Hasanayn, J. Am. Chem. Soc. 110 (1997) 5013. 36. S. E. Barlow, J. M. Van Doren, V. M. Bierbaum, J. Am. Chem. Soc. 110 (1988) 7240; J. W. Larson, T. B. McMahon, J. Am. Chem. Soc. 107 (1985) 766; C. Li, P. Ross, J. E. Szulejko, T. B. McMahon, J. Am. Chem. Soc. 118 (1996) 9360; B. D. Wladkowski, J. I. Brauman, J. Phys. Chem. 97 (1993) 13158. 37. M. J. Frisch, G. W. Trucks, H. B. Schlegel et al., GAUSSIAN 03, Revision B.04, Gaussian, Inc., Pittsburgh, PA, 2003. 38. M. Ernzerhof, G. E. Scuseria, J. Chem. Phys. 109 (1999) 911. 39. M. M. Francl, W. J. Petro, W. J. Hehre, J. S. Binkley, M.-H. Gordon, D. J. DeFree, J. A. Pople, J. Chem. Phys. 77, (1982) 3654. 40. W. C. Swope, H. C. Andersen, P. H. Berend, K. R. Wilson, J. Chem. Phys. 76 (1982) 637. 41. S. Noury, X. Krokidis, F. Fuster, B. Silvi, Topmod Package, 1997. 42. J. Toulouse, A. Savin, C. Adamo, J. Chem. Phys. 117 (2002) 10465. 43. J. J. Blavins, D. L. Cooper, P. B. Karadakov, J. Phys. Chem. A 914 (2004) 108. 44. S. S. Shaik, J. Am. Chem. Soc. 103 (1981) 3692. 45. A. Rastelli, A. S. Pozzoli, G. Del Re, J. Chem. Soc. Perkin Trans 2 (1972) 1571. 46. A. Toro-Labbé, J. Phys. Chem. 103 (1999) 4398.
Ch14-N52719
10-10-2006
8:04 a.m.
Page:287
Trimsize:165×240 MM
Theoretical Aspects of Chemical Reactivity A. Toro-Labbé (Editor) © 2007 Published by Elsevier B.V.
Chapter 14
Quantum chemical topology and reactivity: A comparative static and dynamic study on a SN 2 reaction Laurent Joubert, Ilaria Ciofini, and Carlo Adamo Laboratoire d’Electrochimie et de Chimie Analytique, UMR CNRS 7575, Ecole Nationale Supérieure de Chimie de Paris, 11, rue Pierre et Marie Curie, 75231 Paris Cedex 05, France
Abstract Ab initio molecular dynamic simulations, using density functional theory (DFT) and the recent atom-centered density-matrix propagation (ADMP) method, were employed to study the bond breaking and formation for a case-study SN 2 reaction. Using the real space partition scheme of Bader’s quantum chemical topology (QCT), we performed a population analysis to examine intra- and intermolecular electronic charge transfer along the ADMP trajectory. These results were compared to a static approach, which is performed along the intrinsic reaction coordinate (IRC) path. Although similar features are found for both static and dynamic approaches, the QCT analysis allows to rationalize the differences observed during the formation of the ion–molecule complex. In particular, the dynamic approach suggests a stronger electron exchange tending to spontaneously maximize both covalent and noncovalent (i.e. electrostatic) interactions.
1. Introduction The interpretation of quantum chemical calculations in terms of classical chemical concepts is not an easy task. In particular, electron transfer is strongly related to the definition of atoms and bonds in a molecular system. Beyond the simplistic picture of Lewis1 based on intuitive concepts, two different ways, trying to give quantitative information on the nature of chemical bond, cohabit. On the one hand, classical localization procedures of the wave function give different partitions of the Hilbert’s space 287
Elsevier AMS
Ch14-N52719
288
10-10-2006
8:04 a.m.
Page:288
Trimsize:165×240 MM
Quantum chemical topology and reactivity
leading to atomic charges, such as Mulliken or Löwdin ones,23 or to hybrid orbitals.4−6 More recently, topological analyses based on the electron density, such as Bader’s quantum chemical topology (QCT)78 or on derived functions, such as the electron localization function (ELF),9 have proven their reliability as robust methods, crafted to put in evidence the subtle properties of bonding and reaction mechanisms.10−13 Nowadays, a wide choice of tools are available to chemists for the theoretical analysis of any kind of chemical reactions. These tools are even more powerful when coupled with efficient algorithms for the exploration of the potential energy surface, such as the intrinsic reaction coordinate (IRC) algorithm.1415 Calculation of any topological or electronic variable along the reaction path allows for a clear-cut vision of the reaction mechanism and for rationalizing the structure/energies properties of the surface extrema (reactants, transition states, and products).1617 It must be noticed, anyway, that all these approaches have been developed in the framework of a time-independent (i.e. static) view. Only recently Gross and co-workers have formulated a suitable, but rather complex approach including the time evolution of ELF.18 At the same time, ab intio dynamics, such as the Carr–Parrinello approach,19 have recently reached a mature state, and a blooming of reactivity studies can be found in the literature.20−22 Based on an extended Lagrangian molecular dynamics (MD) scheme, this Born–Oppenheimer approach evaluates the potential energy surface at the density functional theory (DFT) level, and both the electronic and nuclear degrees of freedom are propagated as dynamical variables.23 The calculation of the time variation of any chemicophysical properties (observable or not) is then easily computed using structure snapshots extracted from the trajectory and a statistical average.23 A tentative to explore the difference between static and dynamics ab initio approaches has been done,24 but is has been focused only on the thermochemistry, since the choice of unique and meaningful parameters for comparison is not trivial. In the present article, we report a study concerning the reaction mechanism of a prototype reaction using both static and dynamic approaches to explore a DFT potential surface. The static approach is the standard IRC model, while the dynamic one is based on a Carr–Parrinello method performed with localized (Gaussian) orbitals, the so-called atom-centered density matrix propagation (ADMP) model.25 Our aim is to elucidate the differences, and the common aspects, between the two approaches in the analysis of bond breaking/formation. To this end, we have chosen topological quantities as probe molecular descriptors. The so-called Walden inversion, an SN 2 reaction shown in Figure 1, has been chosen as the model, since it is a well-studied reaction, from both a static and a dynamical point of view.26−36 In the present study, we limit our analysis to the gas-phase reaction, and we do not consider solvent effects that are known to strongly modify the potential energy surface (PES).29
2. Computational details All the calculations were carried out with the Gaussian03 program,37 using the hybrid PBE0 functional3038 and the 6-31+G(d,p) basis set.39 Starting from the transition state, we have calculated reaction pathways for both static and dynamic approaches. On the one hand, we followed the minimum energy path between the transition state and
Elsevier AMS
Ch14-N52719
10-10-2006
C. Adamo et al.
8:04 a.m.
Page:289
Trimsize:165×240 MM
289
the stable ion–molecule complex that is, using mass-weighted coordinates, the internal reaction coordinate (IRC).15 On the other hand, the dynamic simulation was performed at 298 K in the canonical ensemble, for a total simulation time of 1.5 ps. The starting point (transition state) was previously optimized with the same basis set and functional used in the simulation. The well-known Velvet algorithm40 was employed for the integration of the equations of motions using a time step of 0.25 fs. The fictitious mass of the electron was set to 0.20 amu, with a scaling for both core and valence electrons as described in Schlegel et al.25 The velocities of the nuclei were scaled each five time steps to ensure a constant temperature within a T = 5 K of tolerance. The stability of the simulations was monitored by checking at each step the idempotency of the density matrix (within a 10−12 threshold) and the so-called adiabaticity index (within a 10−4 threshold, see Shlegel et al. (2001)25 for more details). In practice, we restricted our trajectory analysis to the time frame of interest, scrutinizing the path only during the formation of the stable ion–molecule complex. All along this truncated trajectory (less than 100 fs in time scale), we extracted a number of snapshots, taken each 2 fs, for the subsequent structural and electronic analyses of the reaction. The charge transfer processes were investigated using a QCT approach: the “Atoms in Molecule” (AIM) partition scheme of the electron charge density.78 According to this theory, topological atoms are defined as regions in real space consisting of a bundle of electron density gradient paths attracted to a nucleus. This partition allows evaluating atomic properties, defined as volume integrals over non-overlapping atomic basins. In particular, the population associated with an atom is simply the volume integral of the electronic charge density over the basin. In this work, we calculated and compared the variations of atomic basin populations all along the selected static and dynamic paths. Moreover, we monitored the variations of typical electron density properties calculated at a bond critical point (BCP) to characterize a chemical bond (see for instance Popelier (2000)8 for more details on BCPs). These calculations were performed using the Topmod package41 and a locally developed code.
3. Results The prototypical SN 2 reaction studied in this work (see figure 1) presents in the gas phase a double-well PES with two equivalent local minima corresponding to the formation of a pre- and a post-reaction ion–molecule complex Cl− · · · CH3 Cl and a transition state (TS) of D3h symmetry Cl · · · CH3 · · · Cl− . Owing to the overall symmetry of the reaction, we will restrain our study to a reaction path from the transition state to one of the equivalent minima.
3.1. Tuning the DFT approach The present study relies on the parameter-free exchange-correlation PBE0 functional. Our choice of the functional was guided by its good performance in predicting thermodynamics data for the current SN 2 reaction, using a large basis set (6-311+(3df,3pd)30 ). Since a large number of calculations are needed in order to evaluate the variations of the topological quantities along the trajectories, such large basis is not suitable. For
Elsevier AMS
Ch14-N52719
10-10-2006
8:04 a.m.
290
Page:290
Trimsize:165×240 MM
Quantum chemical topology and reactivity
[Cl…CH3…Cl]– ΔEovr ClCH3+ Cl–
Energy
Cl– + CH3Cl
ΔE#
ΔEcomp Cl– … CH3Cl
ClCH3…Cl–
Reaction coordinate
Figure 1 Sketch of energy profile for the SN 2 reaction under study
this reason, we performed some preliminary tests to assess the validity of the PBE0 functional with a medium-sized basis set, 6-31+G(d,p). We first examined the geometrical features of the key structures, i.e. the geometries of the ion–molecule complex and of the transition state. Structural and energetic results are summarized in Tables 1 and 2, respectively. From a structural point of view, our data are very close to the MP2 results32 for the transition state structure, with a deviation of about 0.01 Å for the carbon–chlorine distance. The same trends are observed for the C-Cl2 bonding distance in the ion–molecule complex. In contrast, a larger error (0.20 Å) is observed for the C Cl1 distance, this effect being strictly related to the small basis set considered. In fact, when a larger basis, 6-311+G(2d,p), is considered this distance is significantly augmented (3.23 Å) toward the MP2 value. The most significant thermodynamic quantities are the complexation energy of the ion–molecule complex Ecomp , the activation energy, i.e. the relative energy of the D3 h saddle point with respect to the ion–molecule complex E # , and the overall barrier Eovr , defined as the difference between these two energies. These data are reported in Table 2 and compared with previous theoretical29−35 or experimental results.36 Whereas the computation of initial closed-shell reagents does not generate particular difficulties, the determination of the energy of the charged transition state energy Cl · · · CH3 · · · Cl− by DFT approaches is more involved. In fact, the majority of standard functionals, and in particular those resting on the Generalized Gradient Approximation (GGA),
Table 1 Computed geometry parameters (in Å) for the equilibrium ion–dipole complex and the transition state Complex
HFa /TZ3P + R + 2fd MP2a /TZ3P + R + 2fd DFT-LDAa /TZ3P + R + 2fd DFTb -BP/TZ + 2P DFTc -BP/PAW DFTd -mPW1PW/6-31 + Gdp DFT-PBE0/6-31+G(d,p)
Transition state
RC-Cl1
RC-Cl2
RC-C11 = RC-C12
3 37 3 27 2 98 3 10 3 15 3 01 3 07
1 82 1 81 1 81 1 84 1 89 1 83 1 83
2 39 2 32 2 28 2 34 2 37 2 33 2 33
a: Ref [31], b: Ref [32], c: Ref [34], d: Ref. [29].
Elsevier AMS
Ch14-N52719
10-10-2006
8:04 a.m.
Page:291
C. Adamo et al.
Trimsize:165×240 MM
291
Table 2 Complexation energies, overall energies, and barrier heights (in kcal/mol) for the studied SN 2 reaction Ecomp HFa /TZ3P + R + 2fd; ZPE corrected MP2a /TZ3P + R + 2fd; ZPE corrected MP4b TZ + 2P; ZPE corrected G2c ; ZPE corrected G3c ; ZPE corrected CBS-QB3(+)c ; ZPE corrected DFTb -BP/TZ + 2P; ZPE corrected DFTd -BP/PAW DFTe -B3LYP/6-31G(d); ZPE corrected DFTf -mPW1PW/6-31+G(d,p); ZPE corrected DFT-PBE0/6-31+G(d,p); ZPE corrected experimentg
−8 1 −10 5 −10 6 −10 8 −11 2 −10 7 −10 3 −8 3 −9 5 −9 8 −10 1 −12 ± 2
Eovr
E #
7 6 3 5 1 8 3 1 1 8 2 4 −5 7 −3 0 −0 9 0 7 0 5 3/1 ± 1
15 7 14 0 12 4 13 8 13 0 13 1 4 6 5 4 8 7 10 5 10 6 13 ± 2
a: Ref [31], b: Ref [32], c: Ref. [33], d: Ref [34], e: Ref [35], f: Ref. [29], g: Ref. [36].
fail in determining the energy barriers. This effect has been related to the the selfinteraction error in the exchange part as well as in the correlation part, which implies an exaggerated delocalization of the electron density that overstabilize the transition state (for a discussion on this point see Toulouse et al. (2002)42 ). In contrast, PBE0 calculations provide accurate results, close to the post-HF (Hartree-Fock) calculations and in the range of the experimental estimations. Nevertheless, this result can be due to a simple error compensation between the GGA and HF contributions.
3.2. Static and dynamic topological analyses The energy variations along the intrinsic reaction path (IRP) and the ADMP trajectory are reported in Figures 2 and 3, respectively. We want to underline that the IRP corresponds to a unique and constrained minimum energy path and does not contain any time information, whereas the ADMP profile is computed along a trajectory issued by a dynamics simulation. Therefore, any direct (i.e. point-to-point) comparison is misleading. Nevertheless, their global differences can be discussed. The minimum energy constraint imposes a unique and distinctive IRP that continuously decreases from the transition state to the minimum. In contrast, the energy globally decreases along the ADMP trajectory, but it appears clearly that the reaction does not follow this minimum energy path. Actually, the minimum is rapidly reached after 96 fs, but the curve exhibits some energy risings corresponding to punctual destabilizations of the whole system. In order to rationalize the differences observed between the static approach and dynamic process, we investigated the intra- and intermolecular charge transfers along the different reaction paths. As a first step, we looked at the AIM atomic populations for the two states of interest, i.e. the transition state and ion–molecule complex. The results are gathered in Table 3. At the transition state, each chlorine atom bears a negative charge of −0 70e− . This electron excess corresponds to the sharing of the unit
Elsevier AMS
Ch14-N52719
10-10-2006
292
8:04 a.m.
Page:292
Trimsize:165×240 MM
Quantum chemical topology and reactivity
Potential energy (a.u.)
–959.990
–959.995
–960.000
–960.005
–960.010 10
0
20
30
50
40 1/2
MWC step (amu
60
.bohr)
Figure 2 Variations of the potential energy along the IRP from the transition state to the ion–molecule complex
–959.985
Potential energy (a.u.)
–959.990
–959.995
–960.000
–960.005
–960.010 0
12
24
36
48
60
72
84
96
Time (fs)
Figure 3 Variations of the potential energy along an ADMP trajectory from the transition state to ion–molecule complex
charge of the complex increased by a supplementary electron transfer of approximately 0 20 e− from the hydrogens to each of the chlorine atom. Along the IRP, the formation of the ion–dipole molecule corresponds to a global charge transfer of 0 25 e− from the forming chloromethane moiety to the leaving chlorine atom. This charge transfer is slightly increased by 0 01 e− when considering the dynamic process. This global molecule–ion interaction is accompanied by an intramolecular charge transfer in the forming chloromethane molecule. This phenomenon is evidenced by a slight increase of the methyl moiety electronic population, i.e. 0.05 and 0 06 e− at the end of the static
Elsevier AMS
Ch14-N52719
10-10-2006
8:04 a.m.
Page:293
C. Adamo et al.
Trimsize:165×240 MM
293 Table 3 Computed AIM basin populations e− for the transition state and the ion–molecule complex Atom
Transition state
Cl1 H1 H2 H3 C Cl2
Ion–molecule complex
17 70 0 83 0 83 0 83 6 11 17 70
17 95 0 84 0 88 0 86 6 10 17 38
and dynamic process, respectively. A deeper insight into the charge transfer processes can be reached by examining the variations of the atomic electronic populations along the IRP and ADMP trajectory. Figures 4 and 5 represent the relative variations of selected AIM basin populations with respect to the initial populations in the transition state configuration for both the static and dynamic approaches, respectively. The global molecule-to-ion charge transfer is materialized by two essential curves. The first one (squares) corresponds to the monotonically increasing atomic population of the leaving chlorine atom Cl1, tending to the population of a chloride anion. In parallel, a second curve (solid), grouping the atomic populations of the entire molecule CH3 Cl, shows the corresponding and decreasing electronic population of the forming molecule. When the ion–dipole complex is formed, the global electronic transfer between the two moieties corresponds to the aforementioned net charge of 0 25e− or 0 26e− , depending on the approach chosen to study the reactivity of the system, i.e. a static or a dynamic
Variations of AIM basin populations (e–)
0.4 0.3 0.2 0.1 0.0 –0.1 –0.2 –0.3 –0.4 0
10
20
30
40
50
60
MWC step (amu1/2.bohr)
Figure 4 Variations of AIM populations along the IRP. The superbasin corresponding to the forming CH3 Cl molecule (5 basins) is represented by the solid curve. Other curves with markers: carbon basin (lozenges), Cl1 basin (squares), Cl2 basin (circles), and sum of hydrogen basins (triangles) populations
Elsevier AMS
Ch14-N52719
10-10-2006
8:04 a.m.
294
Page:294
Trimsize:165×240 MM
Quantum chemical topology and reactivity
Variations of AIM basin populations (e–)
0.4 0.3 0.2 0.1 0.0 –0.1 –0.2 –0.3 –0.4 0
12
24
36
48
60
72
84
96
Time (fs)
Figure 5 Variations of AIM populations along the ADMP trajectory. All symbols are defined in Figure 4
one. If we now examine the variations of the individual atomic populations, substantial differences are observed between the static and dynamic models of reactivity. Three supplementary curves are presented on the abovementioned Figures 4 and 5. They correspond to the variations of the atomic populations of (i) the chlorine atom Cl1 of the forming chloromethane molecule, (ii) the three gathered hydrogen atoms, and (iii) the carbon atom. In Figure 4, i.e. along the IRP, the detailed analysis of these variations of the atomic populations allows us to propose a three-step charge transfer mechanism. In the first part of the path, corresponding approximately to 1 0 amu1/2 bohr, the electron excess on the chlorine atom Cl2 is transferred almost totally to the other chlorine through the atomic basins of the hydrogens. This particular point is evidenced by examining the variations of the atomic populations of the two chlorine atoms that almost compensate each other. In other words, the electronic flux coming from the Cl2 chlorine basin and entering the three hydrogen basins corresponds to the exiting flux that penetrates the Cl1 chlorine basin. Besides this main electron transfer, a second one is observed between the carbon and the three hydrogen atoms. In fact, along the path, the carbon–hydrogen bond lengths systematically increase when the methyl group adopts a C3v geometry in the forming molecule and points to the leaving chlorine atom. This phenomenon is here increased by a favorable increasing electrostatic interaction with this chlorine atom. Therefore, the lengthening of these bond lengths explains the local charge transfer observed from carbon to hydrogens. Finally, we note substantial fluctuations in the population variations of the linked carbon and chlorine atoms. Again, these fluctuations compensate each other. They result from two concerted effects. First, a small fraction of the electron flux between chlorine atoms penetrates the carbon basin, resulting in a small increase of the corresponding population that tends to be slightly in excess. Second, the formation of the carbon–chlorine bond counters this flux by inducing a small charge transfer from the carbon to the chlorine atom. In the second step of the charge transfer mechanism, i.e. up to 3 0 amu1/2 bohr, the electronic loss on the chlorine atom Cl2
Elsevier AMS
C. Adamo et al.
Ch14-N52719
10-10-2006
8:04 a.m.
Page:295
Trimsize:165×240 MM
295
becomes stronger that the electronic gain on the leaving chlorine Cl1. Therefore, the electron flux that enters the hydrogen basins becomes larger than the exiting flux. As a result, the electronic population in the hydrogen basins increases. At a lesser extent, the electronic charge transfer from carbon to hydrogens increases as well, reinforcing the electronic population of these basins. Furthermore, we note that the fluctuations observed along the carbon and chlorine (Cl2) curves are amplified when approaching the equilibrium structure of the chloromethane. The last step of the charge transfer mechanism corresponds approximately to the second half of the IRP. Structural changes are minimal in the chloromethane molecule, resulting in equilibrium for carbon and hydrogen basin populations. The chlorine atom Cl1 continues to move away from the molecule implying a strong decrease in the chlorine to chlorine charge transfer, the population of the Cl1 atom approaching those of a chloride anion. In Figure 5, i.e. along the ADMP trajectory, the analysis of atomic population variations suggests a charge transfer mechanism at variance with the static IRC mechanism. We can decompose it in five steps. We first note that the system is destabilized during the first 8 fs of the simulation (see Figure 3a), certainly due to an excess of initial kinetic energy. During that period, population variations are negligible, except for the cyclic fluctuations of the carbon and hydrogen populations corresponding to the C-H stretching vibrational mode. These fluctuations are present all along the trajectory. Between 8 and 16 fs, the charge transfer begins between chlorine atoms. This second dynamic step is very similar to the first part of the IRC path, where an electron excess on the chlorine atom Cl2 is transferred almost totally to the other chlorine through the atomic basins of the hydrogens. The third step corresponds roughly to a frame between 16 and 30 fs and can be compared to the second step of the IRC charge transfer mechanism. In fact, the curve corresponding to the population variations of the Cl2 chlorine atom strongly deviates from the one corresponding to the population variations of the whole chloromethane molecule. Meanwhile, both electronic populations of carbon and hydrogen atoms substantially increase, indicating an intramolecular charge transfer from the chlorine Cl2 basin to the carbon and hydrogen basins. In contrast with the static approach, no fluctuations are observed for the carbon and chlorine Cl2 curves, except for the cyclic variations on the carbon curve corresponding to the C-H stretching vibration mode. At this point, we want to remark that all the obtained results are preserved when different initial conditions (kinetic energy and temperature) are considered, since the excess energy mainly changes the population of the CH vibrational states. This is not surprising since the Cl-CH3 Cl starting ion–complex exhibits a rather poor energy transfer between the inter- and intravibrational modes.43 The third step corresponds to a time frame between 30 and 50 fs. We observe that the electronic population of the hydrogen basins strongly increases due to a substantial electron transfer from the carbon basin. Moreover, the corresponding population variations almost compensate each other. All these subtle variations are related to structural changes during the dynamics. In Figure 6, we isolated the population variations of the carbon and hydrogen basins together with the sum of the variations of C-H distances. During the time frame of interest, we note a substantial lengthening of these distances. This phenomenon was also observed along the static IRP, but with a lower amplitude. Here, the C-H stretching vibrational mode elongates the corresponding bond lengths, strengthening the electrostatic attraction between the leaving chlorine atom and hydrogen atoms and favoring the electron flux
Elsevier AMS
Ch14-N52719
10-10-2006
8:04 a.m.
Trimsize:165×240 MM
Quantum chemical topology and reactivity 0.4
0.08
0.3
0.06
0.2
0.04
0.1
0.02 0
0
–0.1
–0.02
–0.2
–0.04
–0.3
–0.06
–0.4
0
12
24
36
48
60
72
84
Variations of the sum of C-H distances (Å)
Variations of AIM basin populations (e–)
296
Page:296
–0.08 96
Time (fs)
Figure 6 Comparison between the variations of selected basin populations (C and H) and the sum of C-H distances (solid line) along a dynamic pathway from the transition state to the ion–molecule complex. Curves with markers: carbon basin (lozenges) and sum of hydrogen basins (triangles) populations
between them. This noncovalent effect should favor the stabilization of the system, but we know that the complex is strongly destabilized within this time frame, as evidenced by a sudden potential energy rising in Figure 3. This surprising behavior can also be rationalized by examining the structural evolution of the system. In Figure 7, we represented the variations of the potential energy together with the evolution of the H-C-H valence angles, i.e. we plotted the variations of the sum of these three angles along the trajectory. In the time frame of interest, we can see that the sum of angles decreases
–959.985
370
Potential energy (a.u.)
350 –959.995 340 –960.000 330 –960.005
–960.010
320
0
12
24
36
48
60
72
84
Sum of H-C-H valence angles (°)
360
–959.990
310 96
Time (fs)
Figure 7 Comparison between the potential energy variations (in bold) and the sum of the three H-C-H valence angles (solid) in the ion–molecule complex
Elsevier AMS
Ch14-N52719
10-10-2006
8:04 a.m.
Page:297
C. Adamo et al.
297 0.5
2.4 2.3
0.4
2.2
0.3
2.1
0.2
2.0 1.9
0.1
1.8
0
1.7 1.6
Rho(rc) et L(rc) (a.u.)
C–Cl2 distance (Å)
Trimsize:165×240 MM
0
12
24
36
48
60
72
84
–0.1 96
Time (fs)
Figure 8 Left Y axis: variations along the ADMP trajectory of the C-Cl2 distance (in bold). Right Y axis: variations of the electron density (solid) and the function Lr (dashed) calculated at the C-Cl2 bond critical point r = rc
to only 310 , a global minimum along the whole trajectory. In other words, the strong electrostatic attraction between the Cl1 chlorine atom and hydrogen atoms induces not only a lengthening of the C-H bond lengths but also a substantial closing of the H-C-H valence angles. The latter effect may explain partly the destabilization of the system by increasing electrostatic and Pauli repulsions between the hydrogen atoms. Furthermore, an examination of electron density properties at the C-Cl2 BCP reveals interesting information on the bond creation process that helps to explain this surprising potential energy rising. In Figure 8, we plotted the variations of the C-Cl2 distance along the trajectory, together with the variations of the density r and function Lr = − r calculated at the corresponding BCP. In the time frame of interest, we note a systematic decrease of the C-Cl2 distance down to the global minimum value of 1.682 Å at = 50fs, about 0.15 Å shorter than the equilibrium value. Meanwhile, the BCP electron density strongly increases to finally reach a value of 0 22 e− bohr −3 , characteristic of a strong covalent bond. The Lr function confirms these trends, with an increasing positive value at BCP that indicates a strong charge concentration. In short, the electron exchange is maximized during the dynamic simulation, favoring the penetration of the chlorine atom electron cloud into the one of the methyl moiety. This effect destabilizes the system by strongly enhancing the Pauli repulsion between these atoms. Finally, the two last steps of the dynamic charge transfer mechanism correspond to the second half of the trajectory. First, between 50 and 78 fs, the system, strongly destabilized, tends to reduce quickly its potential energy. Thus, we note a sudden lengthening of the C-Cl2 bond and an opening of the H-C-H valence angles (see Figure 7) coupled with a shortening of the C-H bond lengths (see Figure 6). This last effect induces a strong diminution of the charge density transferred from the carbon to the hydrogen atoms. Moreover, it is important to point out an inverted electron flux from the ion to the chloromethane molecule tempting to equilibrate the destabilized system. Finally, the last part of the dynamic trajectory, beyond 78 fs, corresponds to the effective formation of a
Elsevier AMS
Ch14-N52719
298
10-10-2006
8:04 a.m.
Page:298
Trimsize:165×240 MM
Quantum chemical topology and reactivity
stable ion–molecule complex and the end of the charge transfer between the molecule and the chlorine atom stabilized by the completion of its external valence shell.
4. Some comments on the static vs. dynamic descriptions Our results show how both ADMP and IRP approaches suggest that the driving force for bonding formation is the charge transfer process. This force, clearly present in the static approach, is so enhanced in the dynamics simulation that it allows for the penetration of the incoming chlorine into electron cloud of the methyl moiety. This results in a very short C-Cl2 distance and in high energy. The subsequent switching on of the Pauli repulsion permits to reach the final state, through a large C-Cl2 amplitude motion. Before this reactive state, the reactants are in a kind of “harmonic reversible” regime where each part still preserves its own electronic identity. This picture is consistent with previous analyses, based on different theoretical models. Among others, we want to recall that similar but maybe more qualitative conclusions were drawn several years ago by Shaik on the basis of molecular orbitals arguments resting on valence bond (VB) calculations.44 In the same philosophy, but more recently, it has been suggested that, within a VB approach, the overlap between the active orbitals of the incoming Cl− and CH3 Cl rules the overall reactive process.43 Moreover, the maximization of the charge transfer strongly reminds the “maximum localization hybrid orbitals” overlap principle suggested more than 30 years ago by Del Re. This criterion has been used as a driving force in a chemical reactivity study based on semi-empirical qualitative approaches.45 Analogous results were also reached by Toro–Labbé using the principle of maximum hardness.46 What is new in our analysis is that the dynamics simulations, based on a ab initio approach using Gaussian orbitals, well underlines this charge transfer phenomenon and its amplitude as driving force. At the same time, we find in a modern and accurate simulation qualitative descriptions and working hypothesis carried out in more approximate schemes. Reassuringly, analyses based on different theoretical backgrounds will converge toward the same description of the same phenomenon, when correctly applied.
5. Conclusion We have analyzed the mechanism of bond breaking and formation in a prototype SN 2 reaction using a static and dynamic approach, based on both DFT and localized (Gaussian) basis sets. Using the real space partition scheme of Bader’s QCT, we performed a population analysis to examine intra- and intermolecular electronic charge transfers along the dynamics trajectory and the static reaction profile. QCT analysis allows the use of well-defined molecular descriptors for the same overall phenomenon (the reaction) in two different frameworks. Our results show that similar features are present in the static and dynamic approaches, but some differences can be observed during the formation of the ion–molecule complex. In particular, the dynamic approach suggests a stronger electron exchange between the different moieties tending to spontaneously maximize both covalent and noncovalent interactions.
Elsevier AMS
Ch14-N52719
C. Adamo et al.
10-10-2006
8:04 a.m.
Page:299
Trimsize:165×240 MM
299
Acknowledgments The authors thank Prof. V. Barone (Naples, Italie) for fruitful discussions and Dr. Michele Pavone (Princeton, USA) for helpful suggestions in the ADMP simulations.
References 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11. 12. 13. 14. 15. 16. 17. 18. 19. 20. 21. 22. 23. 24. 25.
26. 27. 28. 29. 30. 31. 32. 33. 34.
G. N. Lewis, J. Am. Chem. Soc. 38 (1916) 762. R. S. Mulliken, J. Chem. Phys. 23 (1955) 1833. P. O. Löwdin, J. Chem. Phys. 18 (1950) 365. G. Del Re, Theor. Chim. Acta 1 (1963) 188. J. Pipek, P. G. Mezey, J. Chem. Phys. 90 (1989) 4916. A. E. Reed, L. A. Curtiss, F. Weinhold, Chem. Rev. 88 (1988) 899–926. R. F. W. Bader, Atoms in Molecules: A Quantum Theory, Oxford University Press, Oxford, 1990. P. L. A. Popelier, Atoms in Molecules: An Introduction, Pearson Education, London, 2000. A. D. Becke, K. E. Edgecombe, J. Chem. Phys. 92 (1990) 5397. B. Silvi, A. Savin, Nature 371 (1994) 683. S. Berski, J. Andrés, B. Silvi, L. R. Domingo, J. Phys. Chem. A 107 (2003) 6014. P. Geerlings, F. De Proft, W. Langenaeker, Chem. Rev. 103 (2003) 1793. P. L. A. Popelier, L. Joubert, J. Am. Chem. Soc. 124 (2002) 8725. K. Fukui, Acc. Chem. Res. 14 (1981) 363. C. Gonzalez, H. B. Schlegel, J. Chem. Phys. 90 (1989) 2154. M. Solà, A. Toro-Labbé, J. Phys. Chem. A 103 (1999) 8847. V. Polo, J. Andrés, J. Comput. Chem. 26 (2005) 1427. T. Burns, M. A. L. Marques, E. K. U. Gross, Phys. Rev. A 71 (2005) 10501. R. Car, M. Parrinello, Phys. Rev. Lett. 55 (1985) 2471. T. Ziegler, J. Autschbach, Chem. Rev. 105 (2005) 2695. D. Gleich, J. Hutter Chem. Eur. J. 10 (2004) 2435–2444. Y. Tateyama, J. Blumberger, M. Sprik, I. Tavernelli, J. Chem. Phys. 122 (2005) 234505. D. Marx, J. Hutter, in “Modern Methods and Algorithms of Quantum Chemistry Proceedings”, J. Grotendorst (Ed.), Forschungszentrum Jülich, NIC Series, Vol. 3, page 339 (2000). S. C. Ammal, H. Yamataka, M. Aida, M. Dupuis, Science 299 (2003) 1555. H. B. Shlegel, J. M. Millam, S. S. Iyengar, G. A. Voth, A. D. Daniels, G. E. Scuseria, M. J. Frish, J. Chem. Phys. 114 (2001) 9758; S. S. Iyengar, H. B. Shlegel, J. M. Millam, G. A. Voth, G. E. Scuseria, M. J. Frish, J. Chem. Phys. 115 (2001) 10291; H. B. Shlegel, S. S. Iyengar, X. Li, J. M. Millam, G. A. Voth, G. E. Scuseria, M. J. Frish, J. Chem. Phys. 117 (2002) 8694. W. P. Ho, D. G. Truhlar, J. Am. Chem. Soc. 117 (1995) 10726. J. M. Gonzales, R. S. Cox III, S. T. Brown, W. D. Allen, H. F. Schaefer III, J. Phys. Chem. A 105 (2001) 11327. B. Ensing, E. J. Meijer, P. E. Blöchl, and E. J. Baerends, J. Phys. Chem. A 105 (2001) 3300. M. Cossi, C. Adamo, V. Barone, Chem. Phys. Lett. 297 (1998) 1. C. Adamo, V. Barone, J. Chem. Phys. 110 (1999) 6158. B. D. Wladkowski, K. F. Lim, W. D. Allen, J. I. Brauman, J. Am. Chem. Soc. 114 (1992) 9136 and references therein. L. Deng, V. Branchadell, T. Ziegler, J. Am. Chem. Soc. 116 (1994) 10645. S. Parthiban, G. D. Oliveira, J. M. L. Martin, J. Phys. Chem. A 105 (2001) 895. S.-Y. Yang, P. Fleurat-Lessard, I. Hristov, T. Ziegler, J. Phys. Chem A 108 (2004) 9461.
Elsevier AMS
300
Ch14-N52719
10-10-2006
8:04 a.m.
Page:300
Trimsize:165×240 MM
Quantum chemical topology and reactivity
35. A. Streitwieser, G. S. Choy, F. J. Abu-Hasanayn, J. Am. Chem. Soc. 110 (1997) 5013. 36. S. E. Barlow, J. M. Van Doren, V. M. Bierbaum, J. Am. Chem. Soc. 110 (1988) 7240; J. W. Larson, T. B. McMahon, J. Am. Chem. Soc. 107 (1985) 766; C. Li, P. Ross, J. E. Szulejko, T. B. McMahon, J. Am. Chem. Soc. 118 (1996) 9360; B. D. Wladkowski, J. I. Brauman, J. Phys. Chem. 97 (1993) 13158. 37. M. J. Frisch, G. W. Trucks, H. B. Schlegel et al., GAUSSIAN 03, Revision B.04, Gaussian, Inc., Pittsburgh, PA, 2003. 38. M. Ernzerhof, G. E. Scuseria, J. Chem. Phys. 109 (1999) 911. 39. M. M. Francl, W. J. Petro, W. J. Hehre, J. S. Binkley, M.-H. Gordon, D. J. DeFree, J. A. Pople, J. Chem. Phys. 77, (1982) 3654. 40. W. C. Swope, H. C. Andersen, P. H. Berend, K. R. Wilson, J. Chem. Phys. 76 (1982) 637. 41. S. Noury, X. Krokidis, F. Fuster, B. Silvi, Topmod Package, 1997. 42. J. Toulouse, A. Savin, C. Adamo, J. Chem. Phys. 117 (2002) 10465. 43. J. J. Blavins, D. L. Cooper, P. B. Karadakov, J. Phys. Chem. A 914 (2004) 108. 44. S. S. Shaik, J. Am. Chem. Soc. 103 (1981) 3692. 45. A. Rastelli, A. S. Pozzoli, G. Del Re, J. Chem. Soc. Perkin Trans 2 (1972) 1571. 46. A. Toro-Labbé, J. Phys. Chem. 103 (1999) 4398.