HF: an ab Initio

The free energy profiles for the protonation of simple alkanes in SbF5/HF superacid solution have been studied by ab initio molecular dynamics simulat...
32 downloads 0 Views 234KB Size
11596

J. Phys. Chem. B 2002, 106, 11596-11605

Hydrocarbon Reactivity in the Superacid SbF5/HF: an ab Initio Molecular Dynamics Study Simone Raugei*,† and Michael L. Klein‡ International School for AdVanced Studies (SISSA/ISAS), Via Beirut 4, 34014 Trieste, Italy, and Center for Molecular Modeling and Department of Chemistry, UniVersity of PennsylVania, Philadelphia, PennsylVania 19104-6323 ReceiVed: June 27, 2002; In Final Form: August 23, 2002

The free energy profiles for the protonation of simple alkanes in SbF5/HF superacid solution have been studied by ab initio molecular dynamics simulations. It is shown that the presence of extremely reactive protons at high SbF5 concentrations open reaction channels for the C-H protonation with a considerably decreased barrier with respect to dilute solutions and the gas phase. It is proposed that this, along with the progressive presence of neutral SbF5, is responsible for the experimentally observed increased reactivity of these solutions as the concentration of SbF5 is increased. It has also been found that in solution an electrostatic stabilization of the transition state contributes to the lowering of the reaction barrier. In the case of ethane C-H protonation is in competition with C-C protonation. In fact, nearby, the transition state of C-C protonation, alternative reaction pathways that lead to the C-H protonation are accessible.

I. Introduction Since Olah’s pioneering work in the late 1960s, it has wellknown that saturated hydrocarbons react with liquid superacids under mild condition.1-4 The initial step is an oxidation process by the H+ ion to yield the carbenium ion, R+, through a pentacoordinated carbonium ion, RH2+,

RH + H+ h RH2+ h R+ + H2,

(1)

which is formed via a proton attack on a C-H or C-C bond. These nonclassical carbonium ions are characterized by threecenter, two-electron bonded (3c-2e) structures having pentacoordinated carbon atoms and bridging hydrogens:5

where R is an alkyl group and R′ is either an alkyl group or a H atom. On the basis of the reactivity of different alkanes in liquid superacids,2,6 Olah developed the concept of σ-bond basicity (reactivity) with the following qualitative order:7 tertiary-C-H > C-C > secondary-C-H > primary-C-H > CH4. Once formed, the carbenium cation R+ further reacts according to the usual carbocation chemistry (fragmentation, isomerization, and oligomerization). Due to the high reactivity of carbonium ions, most of the evidence supporting their existence relies on indirect observations, such as product selectivities and isotopic exchange, rather than on direct spectroscopy measurements. Only small alkonium ions, such as methonium (CH5+) and ethonium (C2H7+) cations, have been well characterized experimentally in the gas phase.8-10 * Corresponding author. E-mail: [email protected]. † International School for Advanced Studies (SISSA/ISAS). ‡ University of Pennsylvania.

The most striking observation about chemical reactivity in liquid superacid is that in SbF5/HSO3F 1:1 (magic acid) at 80150 °C, methane, which is the most stable alkane, undergoes hydrogen exchange and self-condensation to a trimethyl cation with liberation of H2.11 In SbF5/HF, i.e., the strongest known superacid, methane polycondensation takes place at only 60 °C.7 Much effort has been made to characterize the stationary points of protonated hydrocarbons using standard quantum chemical (gas-phase) methods, such as post-Hartree-Fock and density functional theory (DFT).12-21 These studies indicate that the interconversion among different isomeric structures of carbonium ions (bond-to-bond rearrangement22) is a facile process and should involve a low-energy barrier.23 On the other hand, there are few theoretical studies of alkane protonation using realistic models of liquid superacid media.17,24-26 Recently, Esteves et al.17,24,25 studied the protonation of C-H and C-C bonds of methane, ethane, propane, and isobutane in SbF5/HF clusters at the DFT/B3LYP level with a relativistic core potential for Sb. The superacid moiety was modeled as SbF5:HF and Sb2F11-:H2F+ complexes. They found that the transition states (TS’s) for the C-H and C-C protonations were quite similar to the respective carbonium ions interacting with the anion moiety, suggesting a late transition state. Moreover, in the case of C-H protonation, the proton was basically transferred to the alkane and another one is donated back to the SbF6- or Sb2F11ions. Studies on the protonation of methane in (HF)n-solvated H+ clusters (n ) 0-4), using DFT/B3LYP and MP2 levels of theory,26 showed that only the unsolvated H2F+ cation is strong enough to protonate methane, yielding the methonium ion, CH5+, solvated by HF. In contrast, (HF)1-4-solvated H2F+ chains do not seem to be sufficiently strong to yield solvated CH5+. However, such ions show up as parts of the transition state of the H/H exchange reactions. Many questions of primary importance remain open concerning the nature of the initial step, the true nature of the reaction intermediates, and the observed SbF5 concentration dependence of reactivity. Therefore, it is indispensable to employ reliable

10.1021/jp026395l CCC: $22.00 © 2002 American Chemical Society Published on Web 10/12/2002

Hydrocarbon Reactivity in Superacid SbF5/HF

J. Phys. Chem. B, Vol. 106, No. 44, 2002 11597

theoretical models that fully take into account the solvent at finite temperature. Herein, we report an ab initio (CarParrinello) molecular dynamics27 (CPMD) study of the early stages of C-H protonation of methane, ethane, and propane and C-C protonation of ethane in SbF5/HF solutions at 0 °C. CPMD simulations have already been applied to the study of superacids28-31 and the reactivity of liquid superacids.31 The results obtained show clearly how the technique is capable of revealing microscopic details of reaction mechanisms that are relevant for a deeper understanding of the experiments. This paper is organized as follows. In the next section we will briefly review some aspects of liquid SbF5/HF that are relevant for the ensuing discussion and further comment on our previous results of CPMD simulations on SbF5/HF solutions. In section III we will give the computational details of our simulations. Then, in section IV the results of our study will be given and analyzed. The main results obtained will be summarized in section V along with some additional remarks. II. SbF5/HF Solutions In HF-based superacids the superacidity is attributed to the existence of solvated H2F+ species. For instance, in dilute solutions of up to 10 mol % of SbF5 in HF, SbF5 is almost exclusively fluorinated to form the hexasolvated [SbF6(HF)6]anion and solvated H2F+ cation according to the following ionization reaction:

SbF5 + HF f SbF6- + H2F+

(2)

Infrared32 and 1H NMR spectroscopy33,34 show that the cationic species in the solution are protonated HF chains, Hn+1Fn+, the length of which depends on the concentration of SbF5. With increasing concentration of SbF5, increasing amounts of oligomeric anions (Sb2F11-, Sb3F16-, etc.) are formed. These anions have been investigated experimentally with 19F NMR spectroscopy,33-38 but the actual extent of the oligomeric SbF5 chains is still unclear. However, solutions of SbF5 in HF are known to involve complex equilibria between various SbnF5n+1- anions.39 Recently, Culmann et al.33 found that when the concentration of SbF5 in HF is increased over 20 mol %, increasing amounts of neutral SbF5 are also present. Specifically, at a concentration of SbF5 of 50 mol %, free SbF5 is expected to be 1-2 mol % and to increase abruptly for higher SbF5 concentrations.33 Recent CPMD simulations of SbF5/HF 4, 50, and 100 mol % in SbF530,31,40 agreed well with these experimental findings. The low SbF5 concentration regime is characterized by a fully separated ion pair formed by the solvated excess proton and SbF6-. The 50 mol % SbF5 solution shows short protonated HF chains and SbnF5n+1- anions with n ) 1, 2, and 3, i.e., charged monomers, dimers, and trimers. The trimer was found in further simulations carried out after the publication of ref 31. In this second series of simulations, we also found the presence of Sb2F10 cyclic dimers (Figure 1), which are also very likely to occur in neat SbF5 along with SbnF5n cis-fluorinebridged linear and cyclic chains.40 Moreover, in the 50 mol % SbF5 solution, we also observed protons not solvated by HF molecules, which are almost freely exchanged between two adjacent SbnF5n+1- anions (Figure 1). The presence of this kind of proton, which we will refer to as a “hot proton” in the following, has also been proposed on the basis of 19F NMR experiments.33 We have previously interpreted the formation of the formyl cation in the 50 mol % SbF5 solution in terms of these extremely reactive protons.31 As will be discussed below, this type of proton may also play an important role in the

Figure 1. Snapshot of a configuration of the SbF5/HF solution with 50 mol % SbF5 showing the proton exchange between Sb2F11- and SbF6-. The “hot proton” is labeled as “H+”. A nearby neutral (SbF5)2 dimer is also evidenced. For the sake of clarity all the remaining atoms are not shown. Key: (big dark gray sphere) fluorine; (big bright gray sphere) antimony or carbon; (small bright sphere) hydrogen.

chemistry of hydrocarbons in SbF5/HF. Also, from CPMD studies it is evident that, in the high concentration regime, the H2F+ cations are associated in contact ion pairs SbnF5n+1-:H2F+, where H2F+ is directly bonded to one of the SbnF5n+1- anions via a Sb-F‚‚‚H hydrogen bond. There are few free HF molecules in the concentrated solution. Hence, the contact ion pair cannot evolve into a fully separated ion pair as found in the dilute solution. Instead, H2F+ cations stay close to the SbnF5n+1- anions being separated at most by one HF molecule. H2F+ cations shared between two SbnF5n+1- were also observed. This perfectly agrees with the proposal of the disruption of the protonated HF chains as the concentration of SbF5 increases. These findings are also in agreement with the recent study by Esteves et al.,41 who computed the gas-phase structure and the energies of a series of possible electrophilic species in SbF5/ HF solutions with 50 mol % SbF5. From these calculations, the species H+‚Sb2F11-, H2F+‚Sb2F11-, H3F2+‚Sb2F11-, and H4F3+‚ Sb2F11- were estimated to represent the 36.9%, 16.8%, 36.9%, and 9.4%, respectively, of the protonic species present in the 1:1 solution. III. Computational Details A. Calculation of the Free Energy Profile. CPMD simulations42 combine electronic structure calculation based on DFT with nuclei propagation via molecular dynamics algorithms. The advantage of the method lies essentially in the quality of the interatomic potential at all phase space points. Because the method allows an exact consideration of the anharmonic effects at finite temperature, a detailed study of the energy redistribution among the degrees of freedom as a function of time and of the polarization effects along the reaction pathway on the ground state potential energy surface is feasible. Because bond breaking and forming often require the crossing of a relatively large activation barrier, a method to sample the relevant reaction coordinate must be employed. A huge effort has been mounted to develop sampling techniques of reaction paths in the condensed phase, as is shown by the impressive number of papers devoted to this subject.43-47 To date, there is not any really efficient way to explore the phase space of a reactive system in the condensed phase, which can be used routinely in ab initio simulations, that does not rely at least on an initial guess of a reaction coordinate. From a computational point of view, one of the best techniques to be applied in CPMD simulations is the Blue Moon Ensemble constrained molecular dynamics,46,48 in which the system is forced to move along a

11598 J. Phys. Chem. B, Vol. 106, No. 44, 2002 suitable coordinate. Constrained CPMD simulations have already been successfully applied to the study of many chemical processes.31,49 The crucial technical details about the treatment of the constrained reaction coordinate, ζ, are given in the paper by Sprik and Ciccotti.46 Briefly, the holonomic constraint ζ ) ζ0 is included into the Car-Parrinello Lagrangian, L CP,

L ) L CP + λ(ζ - ζ′) where λ is the Lagrange multiplier associated with the chosen constraint. By straightforward thermodynamic integration, the method allows also the calculation of the free energy change, ∆A, along the reaction path

∆A ) -

∫ζζ

1

0

dζ′





∂H (ζ) ∂ζ

ζ′

where ζ0 and ζ1 are the initial and final values of the reaction coordinate ζ, and H is the Hamiltonian associated with L. The brakets 〈‚‚‚〉ζ′ indicate the ensemble average for the given value of the reaction coordinate, ζ′. The corresponding quantity is the so-called mean force. All the reactions discussed in the present paper were studied using a simple distance constraint ζ ) |RH - RX| between a (randomly) chosen H atom of the superacidic solution and a C atom or geometrical point of the molecule defined as a function of C atom positions as discussed below. This type of constraint, i.e., reaction coordinate, has been employed for all the hydrocarbons considered. In this case, the value of the mean force is simply given by the average of λ,fζ′ ) 〈λ〉ζ′. In our convention, a negative value of the mean force indicates that H is pulled in the direction of the C atom (or of the geometric point x), and a positive value indicates that H and C tend to separate. A zero value of the average of mean force means that we are exploring either (i) a “stationary point” ensemble (a basin of attraction of the potential energy surface or the transition state ensemble) or (ii) the lack of interaction between C and H. Clearly, for many chemical reactions in the liquid phase, the reaction path is a complicated function of several active coordinates of the reactants and solvent.50 Therefore, a simple bond distance between a given pair of atoms can at best give an approximate (overestimated) description of the free energy profile. To explore the reliability of the chosen coordinate, we have studied the protonation of methane in the 50 mol % SbF5 solution using the coordination constraint recently proposed by Sprik.51 With this constraint we focus on the number of hydrogens around the methane C atom, nH, rather than the relative distance between a hydrogen and the carbon. The coordination is gradually changed from nH ) 4 (ζ0 ) 4, i.e., methane) to nH ) 5 (ζ1 ) 5, i.e., methonium ion). We define a continuous proton coordination as NH

ζ ≡ nH )

∑ S(|rHk - rC|)

k)1

where NH is the total number of H atoms in solution. As discussed by Sprik,51 a convenient form for the weighting function, S(r), is the Fermi-like function

S(r) )

1 exp[κ(r - rc)] + 1

Hydrogens in the interval rc - κ-1 < r < rc + κ-1 have a fractional weight, whereas atoms outside are counted as full or not. We have chosen the cutoff distance rc and width κ-1 to be

Raugei and Klein 1.561 and 0.236 Å, respectively. These values have been tuned by matching the C-H pair distribution function of CH5+ at 0 °C in liquid HF with that in the gas phase.52 As will be shown in section IVA, the distance and coordination constraints give comparable free energy barriers for the protonation of methane. We decided to carry out our study with a distance constraint because it does not require any additional parameter (i.e., position and width of the weighting function), which choice could make the use of the coordination constraint problematic for more complicated hydrocarbons. B. Ab Initio Molecular Dynamics. We used the empirical HCTH exchange-correlation functional proposed by Hamprecht et al.53 The validity of employing HCTH in studying SbF5 has already been discussed in previous publications.31,40 TroullierMartins pseudopotentials54 were adopted to describe the ionic cores along with the Kleiman-Bylander decomposition55 for F, and the Gauss-Hermit (GH)56 integration for Sb. For the latter, a good convergence is reached with a 10 point quadrature. The necessity of using the GH integration for Sb was discussed by Kim and Klein.30 The valence electron wave functions were expanded in a plane waves (PW) basis set up to the energy cutoff of 70 Ry. With this choice the molecular geometries converged to within about 0.3%.31,40 The equations of motion have been integrated with a time step of 0.09 fs using a fictitious electronic mass27 of 800 au. The temperature was controlled by coupling the system with a Nose´-Hoover chain of thermostats.57 Further computational details will be given in the following sections whenever they are deemed necessary. C. Simulation Systems. We have studied the reactivity of CH4 in SbF5/HF with 50 and 4 mol % SbF5 and the reactivity of C2H6 and C3H8 in SbF5/HF with 50 mol % SbF5. The starting configurations were obtained from our previous studies on the SbF5/HF solutions.30,31 The SbF5/HF with 4 mol % consists of 26 HF molecules and one SbF5 unit contained in a cubic cell of 9.8797 Å, whereas the solution with 50 mol % consits of 6 HF and 6 SbF5 molecules enclosed in a cubic cell of 10.8839 Å side. Both cells were periodically replicated in the three space directions. The hydrocarbon/SbF5/HF systems were prepared by inserting a hydrocarbon molecule in an equilibrated SbF5/HF solution near a proton chosen at random and increasing the volume of the system to allow for the van der Waals radius of the molecule. The volume was gradually increased in a finite number of steps (between 2 and 4 depending on the size of the molecule), keeping the geometry of the molecule fixed and reequilibrating (for about 0.5 ps) after each increase. When the desired volume was reached, the geometry-constrained molecule was released, the constraint on the chosen reaction coordinate was activated, and the system was further reequilibrated for 1.0 ps. The initial value of the reaction coordinate, i.e., constraint, was set to its (instantaneous) value at the end of the volume-adjusting procedure. The reaction coordinate was sampled forward and backward starting from this point. In each point of the reaction coordinated the system was equilibrated for about 1 ps, then the trajectories were collected for several picoseconds. Specifically, the length of the runs ranged between 2 and 5 ps, depending on the system. IV. Results and Discussion A. C-H Protonation. Gas Phase. We start by discussing C-H protonation of methane and ethane in SbF5/HF clusters, employing the same systems considered by Esteves et al.,24,25 namely the complexes SbF5:HF and Sb2F11-:H2F+.58 The computed structure of the TS’s at the HCTH/PW level are

Hydrocarbon Reactivity in Superacid SbF5/HF

J. Phys. Chem. B, Vol. 106, No. 44, 2002 11599

TABLE 1: Activation Energy (∆E‡) and Enthalpy (∆H‡) at 298 and 273 K (0 °C) for the H/H Exchange in Methane and Ethane Computed for Two Cluster Models of SbF5/HFa CH4 method/quantity, model B3LYP/6-31++G** HCTH/PW(70 Ry cutoff)

Sb2F11-:H2F+

SbF5:HF

Sb2F11-:H2F+

ref

70 77 75 76

90 88 94 91 92

66 61 60 60

82 81 86 84 85

25 25 this work this work this work

∆E‡ ∆H‡ (298 K) ∆E‡ ∆H‡ (298 K) ∆H‡ (273 K)

C2H6

SbF5:HF

a The energies and enthalpies at the HCTH/PW level are calculated using harmonic frequencies obtained by finite difference technique. Units in kJ mol-1.

Figure 2. Four different proton-transfer pathways, ζ. C-H protonation: (a) “hot proton” in the 50 mol % SbF5 solution; (b) proton directly bonded to a Sb-F unit in the 4% solution; (c) pendant proton of an H2F+ cation belonging to the contact ion pair SbnF5n+1-:H2F+ in the 4% solution. C-C protonation in C2H6: (d) distance constraint between a “hot proton” and the middle point, X, of the C-C bond.

similar to those found using B3LYP/6-31++G** by Esteves et al.24,25 In all the TS structures there are two protons forming an H2 moiety and thus a 3c-2e bond with the C atom. In particular, both methods agree on the value of the transition state energies and enthalpies with respect to the reactants (Table 1). We note that the activation free energies for the H/H exchange reaction of methane and ethane increase on going from H+‚SbF6- to H2F+‚SbF6-. For methane in (HF)n-solvated H+ clusters (n ) 1-4), Ahlberg et al.26 found that the activation energy for the H/H exchange with respect to (HF)nH+ + CH4 at B3LYP/6311++G(3df,2p) increases from 1.9 kJ mol-1 for n ) 1 to 98.8 kJ mol-1 for n ) 4. This indicates that the progressive solvation of the excess proton (or, more precisely, of H2F+) makes the proton itself less and less reactive. This may be due to the excess proton “delocalization” along the HF chain (see for example, ref 28). It is important to point out that, as observed by Marx and Parrinello,61 for fluxional molecules, such as carbonium cations (see below), harmonic analysis may be a severely limited tool. In particular, free energies computed from single point calcula-

tions using harmonic frequencies must be used with caution, because the entropic contribution could be affected by a considerable error. Liquid Phase. In the present CPMD study, we have followed the C-H protonation of methane for two different concentrations of SbF5, and of ethane and propane at 50 mol % SbF5. In principle, the protonation of a hydrocarbon in the SbF5/HF superacid might be carried out by (i) a “hot proton” (Figure 2a), (ii) the proton shared between a F atom of a SbF5 unit and HF chains (Figure 2b), (iii) the pendant proton of an H2F+ cation of a SbnF5n+1-:H2F+ contact ion pair (Figure 2c), or in the case of a dilute solution, (iv) an H2F+ cation well separated from the SbnF5n+1- anion in the bulk of the HF solution. The reader is referred to ref 31 for more details. We have investigated route i in SbF5/HF solutions with 50 mol % SbF5 and routes ii and iii for 4 mol % SbF5. As we have discussed in section II and more thoroughly in ref 31, route i for the protonation is very likely because there is a high probability of having SbnF5n+1-:H+ (or SbnF5n+1-:H2F+). We have not studied route iv because the reliability of the distance (or coordination) constraint as a reaction coordinate could break

11600 J. Phys. Chem. B, Vol. 106, No. 44, 2002

Raugei and Klein

Figure 3. C-H protonation of methane in SbF5/HF: (a) mean force, fζ′, and (b) free energy, ∆A, profiles in SbF5/HF with 50 mol % SbF5 (circles) and 4 mol % SbF5 (triangles) at 0 °C. Key: (solid line) transfer of a “hot proton” in CH4 in SbF5/HF with 50 mol % SbF5; (dashed and dot dashed lines) transfer of a proton directly coordinated to a SbF5 moiety and the pendant proton of an H2F+ cation of a SbnF5n+1-: H2F+ contact ion pair in SbF5/HF with 4 mol % SbF5. See Figure 2 for the definition of the reaction coordinate.

TABLE 2: Activation Free Energy (kJ mol-1) Computed at 0 °C for the H/H Exchange Reaction in SbF5/HF Solutionsa hydrocarbon

50 mol % SbF5 route i

CH4 C2H6 C3H8 (C-H)I C3H8 (C-H)II

45 53 55 51

4 mol % SbF5 route i route ii 50

73

a (C-H)I and (C-H)II refer to the protonation of primary and secondary C-H bonds, respectively. See Figure 2 for the definition of the reaction paths.

down as a consequence of the complexity of the components of the reaction coordinate that involves a fast (coupled) excess proton/solvent reorganization, which is not known in advance.50 In Figure 3 the mean force, fζ′, and free energy, ∆A, profiles for these reaction paths computed for methane at 0 °C in liquid SbF5/HF with 50 and 4 mol % SbF5 are shown. With a linear constraint as a reaction coordinate, it was possible to study the free energy profile up to the TS. In fact, immediately after the TS, the CH5+ moiety becomes extremely reactive and promptly transfers back one of the unconstrained protons with the same mechanism proposed for the gas-phase reaction by Esteves et al.18 For the CH5+ cation the same remarks done by Trout and Parrinello31g and Sprik,51 regarding the H3O+ cation in liquid water, hold. The proton-transfer reaction that leads to the formation of RH2+ (or H3O+) has a relatively large activation barrier. On the other hand, the barrier for the reverse reaction is almost negligible, and this reaction can be termed as diffusion limited. Thus, it is a good approximation to consider the activation free energy of the direct reaction very close to reaction free energy itself. This assumption is also consistent with the idea of a late TS.18 The most remarkable feature of Figure 3 is the considerable lowering of the activation energy (Table 2) going from the gasphase clusters25,26 to the condensed-phase calculations when a

Figure 4. C-H protonation of methane in SbF5/HF 50 mol % SbF5 at 0 °C: (a) mean force, fζ′, and (b) free energy, ∆A, profiles computed by constraining the number of hydrogens around the C atom (see text). The incoming proton is a “hot proton”. In the inset of (b), the free energy as a function of the distance of the incoming proton from the C atom is also reported.

“hot proton” in the 50 mol % SbF5 solution is transferred. Similarly, the activation energy for the attack of a proton directly coordinated to a SbF5-HF moiety in the 4 mol % SbF5 solution is slightly higher. On the contrary, the attack of pendant proton (route ii) has an activation free energy (73 kJ mol-1) comparable to the values calculated by Ahlberg et al.26 in (HF)n-solvated H+ clusters with n ) 3 (81 kJ mol-1) and n ) 4 (99 kJ mol-1). This suggests for the methane H/H exchange in dilute solutions an asymptotic DFT value of the activation free energy around 100 kJ mol-1, in agreement with the experimental evidence.26 Before proceeding, we would like to point out that the statistical uncertainty on the height of the barrier computed from constrained dynamics is quite large. The running average of the mean force is converged within 5-10 kJ mol-1Å-1, as inferred from block averages along the trajectory. Thus, using simple error propagation principles, the uncertainty on the activation free energy can be estimated to be as large as 15 kJ mol-1. This is the typical accuracy allowed by ab initio simulations. Moreover, because of the presence of different equivalent protons, the precise specification of the reacting proton could be a limitation in terms of statistical mechanics. To support the results obtained from the distance constraint, we have also computed the free energy profile for the protonation of methane using the coordination constraint described in section IIIA. In Figure 4 the mean force and the free energy as a function of the C atom coordination number, nH, calculated for the transfer of a “hot proton” in SbF5/HF with 50 mol % is shown. Taking the zero of the free energy for nH ) 4.02 (see ref 52), the formation free energy for the CH5+ cation is about 49 kJ mol-1. Hence, assuming that the TS free energy is very close to formation free energy of CH5+, the free energy values obtained with distance and coordination constraints agree very well. This agreement makes us confident that a distance constraint represents, at least up to the TS, a good reaction coordinate. Also for ethane and propane we observe a decrease of the activation energy going from the gas-phase cluster models to liquid SbF5/HF with 50 mol % SbF5 (Tables 1 and 2). In C3H8

Hydrocarbon Reactivity in Superacid SbF5/HF

Figure 5. C-H protonation of ethane in SbF5/HF 50 mol % SbF5 at 0 °C: (a) mean force, fζ′, and (b) free energy, ∆A, profiles. The incoming proton is a “hot proton”.

the protonation of a primary C-H bond requires a slightly higher activation energy (+4 kJ mol-1) than a secondary C-H bond, in agreement with Olah’s σ-bond reactivity order. We also observe that quantum delocalization of protons, not present in our calculations, could lead to a further reduction of the reaction barrier because of quantum tunneling. The lower barrier for the “hot proton” attack may be a consequence of the higher mobility of this kind of proton (or of a proton directly coordinated to a SbF5). As we have previously discussed in section II and ref 31, a “hot proton” is scarcely bound to SbF5 units and, therefore, it can promptly react, as found for (HF)1, 2-H+ + CH4 by Ahlberg et al.26 The situation is different in the dilute solution. The excess proton is delocalized over one or more HF molecules, which strongly interact via a H-bond with other molecules. To protonate a hydrocarbon with such an excess proton, it is necessary to break the H-F covalent bond in the H2F+ and the F‚‚‚H hydrogen bond between the H2F+ cation and the nearby HF molecules. This is suggested by the mean force profiles for the protonation of CH4 shown in Figure 3a. At a large distance, the profiles for the transfer of a proton belonging to an HF chain in dilute SbF5/ HF are more repulsive than those for the transfer of an “hot proton” in the concentrated solution. This is particularly true for a pendant proton in the dilute solution. Of course, in the high concentration regime, there are different types of protons, other that the “hot protons”, that can attack a hydrocarbon. Hence, the actual free energy barrier (e.g., the barrier experimentally observable) will be between the value computed for the “hot proton” and a HF-solvated proton attack. From the simulations is evident a progressive intensification of the H‚‚‚F interactions between the forming carbonium cation and the F atoms of the solution, as can be qualitatively inferred from Figure 7 where representative configurations of the transition state ensemble for the protonation of CH4 and C2H6 are shown. More quantitative information is obtained from Figure 8 where the radial distribution functions, g(r)’s, between the H atoms of the CH5 or CH4 moieties and the nearby F atoms of SbF5 moieties calculated for different values of the reaction coordinate are shown. As the incoming proton approaches the carbon atom, the g(r) progressively shifts toward shorter

J. Phys. Chem. B, Vol. 106, No. 44, 2002 11601

Figure 6. C-H protonation of propane in SbF5/HF 50 mol % SbF5 at 0 °C: (a) mean force, fζ′, and (b) free energy, ∆A, profiles. The incoming proton is a “hot proton”. Key: (continuous line) primary C-H bond; (dashed line) secondary C-H bond.

Figure 7. Selected configurations near the transition state ensemble for the protonation of CH4 and C2H6 in SbF5/HF 50 mol % SbF5 at 0°. The numbers indicate the depicted instantaneous F‚‚‚H distances. Key: (big dark gray sphere) fluorine; (big bright gray sphere) antimony or carbon; (small bright sphere) hydrogen.

distances, indicating a progressive intensification of the H‚‚‚F interactions. On average every hydrogen of a protonated carbon has a nearby fluorine atom of a SbF5 moiety, as can be seen from the instantaneous configurations shown in Figure 7. These interactions can be explained in terms of charge distribution around the forming carbonium cations. The atomic charges computed at B3LYP/6-311++G(d,p) for the species involved in the chemistry of CH4 and C2H6 in liquid SbF5/HF using both natural bond orbital and Mulliken population analyses are shown in Table 3. As can be seen, the excess charge on the carbonium cations is distributed over all the hydrogens of a protonated carbon with a slightly larger value on those of the H-H moiety. The same is true for the C-H protonation of propane (data not shown). As a result, there is an increase of the charge on the hydrogens going from the nonprotonated to the protonated alkane. As already mentioned, the C-H protonation is characterized by a late TS, which resembles very much the product. Therefore, near the TS the hydrogens of the forming carbenium cation can effectively interact with the F atoms of the surrounding SbF5 moieties or HF molecules. Of course, for bigger and strerically more hindered hydrocarbons this stabilization will be less important.

11602 J. Phys. Chem. B, Vol. 106, No. 44, 2002

Figure 8. Radial distribution functions between the hydrogens of a (C-H) protonated carbon and the surrounding fluorines calculated at different values of the reaction coordinate. To avoid a spurious peak, the contact with the fluorine to which the incoming H was coordinated has been removed. In the legend of each panel the numbers indicate the value of the reaction coordinate in Å. (C-H)I and (C-H)II stand for protonation of a primary and secondary C-H bond, respectively.

TABLE 3: Atomic Charges Computed at the B3LYP/ 6-311++G(d,p) Level for the Species Involved in the Chemistry of CH4 and C2H6 in Liquid SbF5/HFa

a The columns NBO and Mulliken report the charges derived from natural bond orbital and Mulliken population analysis, respectively.

A more careful analysis of the H‚‚‚F interactions in the 50 mol % SbF5 solution in terms of the lone pairs on the F atoms has also been carried out. It is evident that the H‚‚‚F axis points

Raugei and Klein toward the center of charge of a lone pair (computed as the center of charge of Wannier function associated with it, as briefly discussed below). This finding suggests a hydrogen bond interaction. However, the lifetime of this hypothetical bond is very short. On the time span by our simulations, only in the case of CH5+ does the longest H‚‚‚F “bond” last for about 1 ps, which is actually a value comparable with the H-bond lifetime in neat water. For the other hydrocarbons this time is much shorter (about 0.2-0.5 ps). In general, a hydrogen stays close to a certain F for few hundreds of femtoseconds and then “jumps” to another F atom. On the basis of these observations, we would rule out a genuine H-bond (i.e., accompanied by a charge transfer) in favor of an electrostatic H-bond-like interaction. C-H-protonated alkanes are fluxional molecules. The most theoretically studied carbonium ion, CH5+, reveals no equilibrium structure in the gas phase because of quantum effects and is characterized by a fast scrambling motion of the H atoms.61,62 Also in superacidic solution, near the TS, there is a considerable scrambling between the H atoms of the protonated moiety. The observed scrambling is faster for CH5+ than for the other hydrocarbons. The scrambling velocities have the following qualitative order: CH4 > C2H6 > C3H8 (primary-CH) > C3H8 (secondary-CH). With the idea of a late TS in mind, it is interesting to note that in the TS for the protonation of methane the scrambling motion is slower than in CH5+ in the gas phase.62 This is clearly due to the above-discussed solvent effect, as already discussed by Tse et al.62 for the clusters CH5+(H2)n, and partially to the presence of the constraint. To better characterize the chemical nature of the TSs, we have investigated the wave functions obtained in the simulation and compared them to the wave functions of the minimum structure of protonated hydrocarbons in the gas phase. The maximally localized Wannier functions technique has been employed.63 All the TSs are characterized by a H2 moiety strongly bound to a CHn moiety. In the forming CH5+, the H2 moiety breaks and forms on a time scale of few picoseconds. The center of charge of the 3c-2e bond is closer to the H2 moiety than to the C atom. In Figure 9 the radial distribution functions between the carbon atom of a CH5+ or CH4+ moiety and the center of charge of the Wannier orbitals associated with a 2c2e C-H bond or a 3c-2e C-[H2] bond are shown. The 3c-2e orbital gives rise to the farther peak. For the protonated methane the scrambling is faster and the splitting between the two peaks is smaller. In Figure 9 with arrows the distances between the C atom and the center of charge of a given orbital of optimized gas-phase cation structures are also shown. The comparison of gas-phase data with the finite temperature simulation, keeping in mind the fluxional nature of these ions, gives another indication of a late TS with a structure very close to that of the product. The results of the present and previous31,40 studies are fully consistent with several experimental findings. Sommer and coworkers6 observed for isobutane in SbF5/HF in the presence of CO at 0 °C an abrupt increase of the production of pivalic ester and a decrease of the H2 production as the concentration of SbF5 is increased above 20 mol %. They interpreted these findings in terms of the reduction of RH2+ by SbF5 that becomes more and more available for concentrations above 20%.33 In these conditions the protolytic reaction 1 is less favorable than the oxidative reaction

RH2+ + SbF5 h R+ + SbF3 + 2HF

(3)

Hydrocarbon Reactivity in Superacid SbF5/HF

Figure 9. Radial distribution functions, g(r), between the carbon atom of a CH5+ or CH4+ moiety and the center of charge of the Wannier orbitals associated with a 2c-2e C-H bond (e) or a 3c-2e C-[H2] bond (e′): (a) CH4 (continuous line, ζ ) 1.208 Å; dashed line, ζ ) 2.500 Å); (b) C2H6 (continuous line, ζ ) 1.275; dashed line, ζ ) 2.100 Å). The arrows indicate the corresponding distance in the minimum gasphase structure. The g(r)s are computed from 10 configurations equispaced in time. In the structures to the right, the small dark spheres represent the center of charge of Wannier orbitals.

operated by SbF5 with production of SbF3 and HF, as earlier suggested by Olah et al.65 Our results support this proposal. In fact, at high SbF5 concentration, we have found the presence of (SbF5)2 near a “hot proton” (see Figure 1). Clearly, the limited size of our simulation cell inhibits the formation of larger quantities of neutral SbF5 aggregates and we cannot quantify their role. In this regard, it is worth recalling that our simulation on pure liquid SbF5 indicated a large presence of SbF5 monomers and dimers along with more extended linear and cyclic chains. Most importantly, the results of the present study suggest that the abrupt increase in the isobutane reactivity at high SbF5 concentration may also be a consequence of the increased number of both “free” SbF5 and “hot protons”. The increased quantity of carbonium cations and of neutral SbF5 could drastically increase the reactivity of hydrocarbons at high SbF5 concentrations. B. C-C Protonation. In the gas phase, we studied the C-C protonation of C2H6 only for the H+‚SbF6- cluster model. Also in this case, the most relevant geometric parameters of the TS structure agree fairly well with those of Esteves et al.24 The activation enthalphy at the HCTH/PW level computed at 298 and 273 K is 105 and 103 kJ mol-1, respectively, and has to be compared with 120 kJ mol-1 calculated by Esteves et al.24 (these authors found an activation enthalphy of 146 kJ mol-1 for the H2F+‚Sb2F11- cluster model). To study the free energy profile in the condensed phase, we have started constraining the distance between a selected “hot proton” and the middle point of the C-C bond, as schematically shown in Figure 2d. With this constraint the C-C protonation was not observed. A C-H protonation takes place instead. In fact, as the system approaches ζ ) 1.275 Å (approximately around the extremum of the mean force profile, Figure 10), a stronger and stronger H‚‚‚H interaction characterizes the C2H6:H+ moiety. This is evident from a direct inspection of the distance between the incoming proton and the hydrogens of the ethane methyl groups (data not shown).

J. Phys. Chem. B, Vol. 106, No. 44, 2002 11603

Figure 10. C-C protonation of ethane in SbF5/HF 50 mol % SbF5 at 0 °C. The chosen reaction coordinate, ζ, is the distance between the incoming proton and the middle point of the C-C bond. (a) Mean force, fζ′, profiles from distance constraint without (empty triangles) and with (empty circles) the H-X-C angle constrained to 90°; with filled circles the average value of the Lagrange multiplier associated with the angle constraint is also reported (units kJ mol-1deg-1). (b) Free energy, ∆A, profile obtained from the dynamics with the H-X-C angle constrained. The incoming proton is a “hot proton”.

To observe the C-C protonation, we had to constraint the H-X-C angle to 90° as well. The free energy profile shown in Figure 10 is characterized by a maximum at ζ ) 1.20 Å and a minimum corresponding to the solvated [CH3-H-CH3]+ cation. Because a second constraint was employed, it is deemed necessary to verify the nature of the maximum. Several free trajectories were monitored starting from configurations equispaced in time taken from the runs with ζ ) 1.15, 1.20, and 1.25 Å, corresponding to the point immediately after (A), very near (N), and before (B) the position of the maximum, respectively. All the trajectories starting from A and B (10 for each value of ζ) relaxed to the [CH3-H-CH3]+ cation and to C2H6 plus solvated H+, respectively. The trajectories starting from N (for a total of 16) relaxed mostly to separated C2H6 + H+ (8 out of 16) or the C-C protonated ethane (6 out of 16), and in a minor amount to the C-H (2 out of 16) protonated ethane. Although the number of free trajectories we ran is relatively small, these results indicate that the maximum along the free energy profile with the constrained angle can be considered as belonging to a transition state ensemble that connects the free to the C-C protonated ethane. We would like also to point out that the value of the mean force associated with the angle constraint is very small and oscillates around zero (Figure 10). These findings lead one to conclude that a proton attack perpendicular to the C-C bond represents one possible minimum reaction path, or at least a good approximation to it. It is also important to point out that the constrained molecular dynamics implies a quasi-equilibrium transformation. During the exploration of the constrained phase space it is possible to access another viable reaction channel, which requires an activation energy of the order of a few kT. Thus, it is evident that near the maximum of the mean force there is a viable reaction channel that leads to the C-H protonated isomer. Hence, as the proton approaches the C-C bond, small fluctuations in the C-X-H angle can result in the C-H protonated alkane.

11604 J. Phys. Chem. B, Vol. 106, No. 44, 2002

Figure 11. Bond to bond rearrangement of the C-C protonated ethane to C-H protonated ethane: (a) mean force, fζ′, and free energy, ∆A, profiles as a function of the distance between the bridging proton and the middle point of the C-C bond.

The computed activation free energy is ∆A‡ ) 78 kJ mol-1 (Figure 10). The value of the activation energy of ∆E‡ ) 100 kJ mol-1, computed in the gas phase for the H+‚SbF6- cluster model, indicates that also for the C-C protonation the attack of a “hot proton” has a low activation barrier, even though not as pronounced as for C-H protonation. Whereas the C-H protonation gives rise to unstable species that immediately decompose, the C-C protonated ethane is stable in superacid solution on the time scale of our simulation. In fact, the reverse reaction presents a relatively high activation free energy. A free (without any constraint) 15 ps long CPMD simulation of [CH3-H-CH3]+ in SbF5/HF has not shown other easily accessible reaction channels with activation energies within kT. Also, the increase of the temperature to 500 K has not disclosed any thermally accessible reaction channel. To study one of the most likely reaction channels of [CH3H-CH3]+, namely the bond to bond rearrangement with the formation of the C-H protonated ethane, the X‚‚‚H distance, ζ, was once again used as a reaction coordinate. Starting from [CH3-H-CH3]+, we have increased ζ and the resulting free energy profile is shown in Figure 11. At ζ ) 1.1 Å a TS is observed. The symmetry of this TS is similar to that computed for the rearrangement in the gas phase.13 As observed above, the H‚‚‚H interaction seems to drive the system toward the C-H protonated isomer. Immediately after the TS the protons of the CH4+ moiety become so reactive that the forming C-H protonated isomer decomposes into C2H6 + “H+”. The computed activation free energy for the isomerization process of ∆A‡ ) 12 kJ mol-1 must be contrasted with the activation energy ∆E‡ ) 19 kJ mol-1 computed in the gas phase at the same level of theory (∆E‡ ) 30 kJ mol-1 at MP2/6-311G** level13). V. Conclusions The protonation of simple alkanes in the liquid superacid SbF5/HF has been investigated by using ab initio molecular dynamics studies. The model systems we considered allowed us to rationalize some important aspects of the chemical

Raugei and Klein reactivity as a function of the SbF5 concentration. Highly concentrated SbF5/HF solutions are characterized by the presence of extremely reactive “hot protons” that are almost freely exchanged between SbF6 units. In the high concentration regime it is also very likely to have the presence of nonionized SbF5, which can be involved in the oxidative transformation of a newly formed carbonium cation RH2+ to a carbenium cation R+ with formation of H2 and SbF3.2,6 For the C-H protonation, it has been found that the activation free energy for the transfer of a scarcely bound “hot proton” is considerably lower than that of an excess proton in a HF chain. Although it is very likely to have such protons for high concentrations of SbF5, in solution there are many other topologically different protons that can attack a hydrocarbon molecule. In particular, for methane our calculations predict a rapid increase of the activation free energy as the transferred proton is more and more solvated by HF molecules. Therefore, the experimentally measured protonation barrier will be between the lower limit value represented by the “hot proton” attack and the upper limit value represented by the attack conducted by a protonated HF chain. We also remark that the computed activation free energies for “hot proton” attack on a C-H bond in methane, ethane, and propane have comparable values, indicating the high reactivity of these protons. We would like also to stress that the importance of the “hot protons” diminishes as the SbF5 is reduced. CPMD simulations have revealed a highly dynamical nature of the transition state for the C-H protonation. The transition state ensemble is characterized by H-bond-like interactions between the hydrogens of a protonated alkane moiety and the nearby F atoms. These very weak bonds continuously break and form (the average lifetime is shorter than 1 ps). This type of interaction can be explained by observing that, in the transition state, the excess charge on the forming carbonium cation is distributed over all the H atoms belonging to the C atom involved in the reaction, which therefore can favorably interact with the nearby electronegative fluorines. The study of the C-C protonation of ethane has revealed that from the minimum reaction path, which corresponds to an approach of the proton nearly perpendicular to the C-C bond, the system can easily gain access to reaction channels that lead to the C-H protonate isomer. We are aware that the small size of the systems studied will impose limitations in describing structural features of the high concentration solution. However, it is reasonable to assume that the local structure around an SbF6 unit is, at least qualitatively, well reproduced and, hence, that the main features of the protonation mechanism are well described. The statistical uncertainty represents another source of error that makes the free energy profiles reported qualitative rather than quantitative. Nevertheless, the results obtained give a vivid description of the first stages of chemistry of small hydrocarbons in a real superacid solution. Acknowledgment. The present study was supported in part by the National Science Foundation (NSF) under the grant NSF CHE0205146. Some of the calculations were carried out using Lemieux, the Terascale Computing System installed at Pittsburgh Supercomputing Center. We thank NPACI for generous allocation of CPU time. References and Notes (1) Pines, H. The chemistry of catalytic hydrocarbon conVersions; Academic Press: New York, 1981. (2) Olah, G. A.; Prakash, G. K.; Sommers, J. Superacids; John Wiley and Sons: New York, 1985.

Hydrocarbon Reactivity in Superacid SbF5/HF (3) Olah, G. A. Angew. Chem., Int. Ed. Engl. 1993, 32, 767. (4) Olah, G. A. J. Org. Chem. 2001, 66, 5943. (5) McMurry, J. E.; Lectka, T. Acc. Chem. Res. 1992, 25, 47. (6) Sommer, J.; Bukala, J.; Hachoumy, M.; Jost, R. J. Am. Chem. Soc. 1997, 119, 3274. (7) Olah, G. A.; Prakash, G. K. S.; Williams, R. E.; Field, L. D.; Wade, K. Hypercarbon Chemistry; Wiley: New York, 1987. (8) Field, F. H.; Munson, M. S. B. J. Am. Chem. Soc. 1965, 87, 3289. (9) French, M.; Kebarle, P. Can. J. Chem. 1975, 53, 2268. (10) Yeh, L. I.; Rice, J. M.; Lee, Y. T. J. Am. Chem. Soc. 1989, 111, 5597. (11) Olah, G. A.; Schlosberg, R. H. J. Am. Chem. Soc. 1968, 90, 2726. Olah, G. A.; Klopman, G.; Schlosberg, R. H. J. Am. Chem. Soc. 1969, 91, 3261. (12) Sieber, S.; Buzek, P.; Schleyer, P. v. R.; Koch, W.; de M. Carneiro, J. W. J. Am. Chem. Soc. 1993, 115, 259. (13) de M. Carneiro, J. W.; von Schleyer, P. R.; Saunders, M.; Remington, R.; H. F. S., III; Rauk, A.; Sorensen, T. S. J. Am. Chem. Soc. 1994, 116, 3483. (14) Collins, J.; O’Malley, P. J. J. Chem. Soc., Faraday Trans. 1996, 92, 4347. (15) Mota, C. J. A.; Esteves, P. M.; Ramı´rez-Solı´s, A.; Herna´ndezLamoneda, R. J. Am. Chem. Soc. 1997, 119, 5193. (16) Esteves, P. M.; Mota, C. J. A.; Ramı´rez-Solı´s, A.; Herna´ndezLamoneda, R. J. Am. Chem. Soc. 1998, 120, 3213. (17) Esteves, P. M.; Alberto, G. G. P.; Ramı´rez-Solı´s, A.; Mota, C. J. A. J. Am. Chem. Soc. 1999, 121, 7345. (18) Esteves, P. M.; Alberto, G. G. P.; Ramı´rez-Solı´s, A.; Mota, C. J. A. J. Phys. Chem. A 2000, 104, 6233. (19) Okulik, N.; Peruchena, N.; Esteves, P. M.; Mota, C. J. A.; Jubert, A. H. J. Phys. Chem. A 2000, 104, 7586. (20) Esteves, P. M.; Alberto, G. G. P.; Ramı´rez-Solı´s, A.; Mota, C. J. A. J. Phys. Chem. A 2001, 105, 4308. (21) Okulik, N. B.; Sosa, L. G.; Esteves, P. M.; Mota, C. J.; Jubert, A. H.; Peruchena, N. M. J. Phys. Chem. A 2002, 106, 1584. (22) Olah, G. A.; Rasul, G. Acc. Chem. Res. 1996, 30, 245. (23) Goeppert, A.; Sassi, A.; Sommer, J.; Esteves, P. M.; Mota, C. J. A.; Karlsson, A.; Ahlberg, P. J. Am. Chem. Soc. 1999, 121, 10628. (24) Esteves, P. M.; Ramı´rez-Solı´s, A.; Mota, C. J. A. J. Braz. Chem. Soc. 2000, 11, 345. (25) Esteves, P. M.; Ramı´rez-Solı´s, A., Mota, C. J. A. J. Phys. Chem. B 2001, 105, 4331. (26) Ahlberg, P.; Karlsson, A.; Goeppert, A.; Lill, S. O. N.; Dine´r, P.; Sommer, J. Chem. Eur. J. 2001, 7, 1936. Ahlberg, P.; Karlsson, A.; Goeppert, A.; Lill, S. O. N.; Dine´r, P.; Sommer, J. Chem. Eur. J. 2001, 7, 2501 (erratum). (27) Car, R.; Parrinello, M. Phys. ReV. Lett. 1985, 55, 2471. (28) Kim, D.; Klein, M. L. J. Am. Chem. Soc. 1999, 121, 11251. (29) Kim, D.; Klein, M. L. Chem. Phys. Lett. 1999, 308, 235. (30) Kim, D.; Klein, M. L. J. Phys. Chem. B 2000, 104, 10074. (31) Raugei, S.; Klein, M. L. J. Phys. Chem. B 2001, 105, 8212. (32) Bonnet, B.; Mascherpa, G. Inorg. Chem. 1980, 19, 785. (33) Culmann, J.-C.; Fauconet, M.; Jost, R.; Sommer, J. New. J. Chem. 1999, 23, 863. (34) Gillespie, R. J.; Moss, K. C. J. Chem. Soc. A 1996, 1966, 1770. (35) Dean, P. A. W.; Gillespie, R. J. J. Am. Chem. Soc. 1969, 91, 7260. (36) Dean, P. A. W.; Gillespie, R. J. J. Am. Chem. Soc. 1969, 91, 7264. (37) Bacon, J.; Dean, P. A. W.; Gillespie, R. J. Can. J. Chem. 1969, 47, 1655. (38) Bacon, J.; Dean, P. A. W.; Gillespie, R. J. Can. J. Chem. 1970, 48, 3413. (39) Donnell, T. A. Superacids and acids melts as inorganic chemical reaction media; VCH Publisher: New York, 1993. (40) Raugei, S.; Klein, M. L. J. Chem. Phys. 2002, 116, 7087. (41) Esteves, P. M.; Ramirez-rez Solı´s, A.; Mota, C. J. A. J. Am. Chem. Soc. 2002, 124, 2672.

J. Phys. Chem. B, Vol. 106, No. 44, 2002 11605 (42) Marx, D.; Hutter, J. In Modern Methods and algorithms of quantum chemistry; Grotendorst, J., Ed.; NIC series; John von Neumann Institute for Computing: Ju¨lich, 2000; Vol. 1, pp 301-349. (43) Passerone, D.; Parrinello, M. Phys. ReV. Lett. 2001, 87, 108302. (44) Elber, R. Reaction path studies of biological molecules. In Recent deVelopments in theoretical studies of proteins; Elber, R., Ed.; World Scientific: Singapore, 1996. (45) Bolhuis, P. G.; Chandler, D.; Dellago, C.; Geissler, P. L. Annu. ReV. Phys. Chem. 2002, 53, 291. (46) Sprik, M.; Ciccotti, G. J. Chem. Phys. 1998, 109, 7737. (47) Jarzynski, C. Phys. ReV. Lett. 1997, 78, 2690. (48) Blo¨ch, P. E.; Senn, H. M.; Togni, A. In Transition state modelling for catalysis; Truhlar, D. G., Morokuma, K., Eds.; ACS Symposium Series, American Chemical Society, Washington, DC, 1999; pp 89-99. (49) (a) Raugei, S.; Cardini, G.; Schettino, V. J. Chem. Phys. 2001, 114, 4089. (b) Ensing, B.; Meijer, P. E. B. E. J.; Baerends, E. J. J. Phys. Chem. A 2001, 105, 3300. (c) Boero, M.; Parrinello, M.; Hu¨ffer, S.; Weiss, H. J. Am. Chem. Soc. 2000, 122, 501. (d) Boero, M.; Morikawa, Y.; Terakura, K.; Ozeki, M. J. Chem. Phys. 2000, 112, 9549. (e) Mundy, C. J.; Hutter, J.; Parrinello, M. J. Am. Chem. Soc. 2000, 122, 4837. (f) Mortensen, J. J.; Parrinello, M. J. Phys. Chem. A 2000, 104, 2901. (g) Trout, B. L.; Parrinello, M. J. Phys. Chem. B 1999, 103, 7340. (h) Meijer, E.; Sprik, M. J. Am. Chem. Soc. 1998, 120, 6345. (i) Meijer, E.; Sprik, M. J. Phys. Chem. A 1998, 102, 2893. (j) Boero, M.; Parrinello, M.; Terakura, K. J. Am. Chem. Soc. 1998, 120, 2746. (k) Raugei, S.; Cardini, G.; Schettino, V. J. Chem. Phys. 1999, 111, 10887. (50) For a detailed discussion, see for example: Geissler, P. L.; Dellago, C.; Chandler, D. J. Phys. Chem. B 1999, 103, 3706. See also ref 45. (51) Sprik, M. Chem. Phys. 2000, 258, 139. (52) As noted by Sprik,51 the finite width of the weigth function S(r) gives rise to huge restoring forces. This makes it impossible to study the system for nH ) 5 (or nH ) 4), and these preliminary simulations were performed with nH ) 4.99. This also represents the maximum value of the coordination in the study of the protonation of CH4 in SbF5/HF. The other extremum for the study of this reaction was nH ) 4.02. (53) Hamprecht, F. A.; Cohen, A. J.; Tozer, D. J.; Handy, N. C. J. Chem. Phys. 1998, 109, 6264. (54) Troullier, N.; Martins, J. L. Phys. ReV. B 1991, 43, 1993. (55) Kleinman, L.; Bylander, D. M. Phys. ReV. Lett. 1982, 48, 1425. (56) Stich, I.; Car, R.; Parrinello, M.; Baroni, S. Phys. ReV. B 1989, 39, 4997. (57) Martyna, G. J.; Klein, M. L.; Tuckerman, M. J. Chem. Phys. 1992, 97, 2635. (58) The calculations were performed with the HCTH exchange and correlation functional in a PW basis set (70 Ry cutoff) in a cubic box of 15 Å side. The interaction with periodic images was filtered out using the Barnett-Landman implementation of the Hockney method.59,60 (59) Hockney, R. W. Methods Comput. Phys. 1970, 9, 136. (60) Barnett, R. N.; Landman, U. Phys. ReV. B 1993 48, 2081. (61) Marx, D.; Parrinello, M. Science 1999, 286, 1051. Marx, D.; Parrinello, M. Z. Phys. D 1997, 41, 253. Marx, D.; Parrinello, M. Nature 1995, 375, 216. (62) Tse, J. S.; Klung, D. D.; Laasonen, K. Phys. ReV. Lett. 1995, 74, 876. (63) The wave functions in the PW calculations are completely delocalized and therefore the maximally localized Wannier functions technique64 has been used to facilitate the interpretation of the electronic structure in terms of doubly occupied localized orbitals. However, because the full information contained in the Wannier function is difficult to represent, we have drastically contracted this information and monitored only the centers of charge of the Wannier functions, the so-called Wannier Functions Centers. (64) Silvestrelli, P. L.; Marzari, N.; Vanderbilt, D.; Parrinello, M. Solid State Commun. 1998, 107, 7. (65) Olah, G. A.; Halpern, Y.; Mo, Y. K. J. Am. Chem. Soc. 1973, 95, 4960.