Benchmark Study on the Smallest Bimolecular Nucleophilic Substitution Reaction: H−+CH4 → CH4+H−

We report here a benchmark study on the bimolecular nucleophilic substitution (SN2) reaction between hydride and methane, for which we have obtained reference energies at the coupled cluster toward full configuration-interaction limit (CC-cf/CBS). Several wavefunction (HF, MP2, coupled cluster) and density functional methods are compared for their reliability regarding these reference data.


Introduction
Bimolecular substitution (S N 2) reactions play an important role in organic chemistry and in biochemistry (DNA replication mechanism). Interestingly, there is a profound solvent effect present which has a major effect on the reaction barriers and intermediates. For example, the prototypical S N 2 OPEN ACCESS reaction of chloride with methyl chloride shows in the gas phase a double-well potential (see Figure 1) with deep wells and a reduced barrier. On the other hand, in solution the energy profile turns basically into a unimodal reaction [1][2][3][4][5][6][7][8][9][10][11][12] (see Figure 1), accompanied by a significant increase of the reaction barrier. In previous studies  it was shown that coupled cluster methods in general give accurate results for the energy profile of S N 2 reactions, while density functional methods give qualitatively correct results but often underestimate barriers [15]. This has led to the design of new and improved density functionals (SSB-D [43], S12g [44] and S12h [44]), where in particular the latter hybrid functional (S12h) was shown to provide accurate results for the complete energy profile of S N 2 reactions. Here we have studied the smallest S N 2 reaction possible, between hydride and methane: For this reaction, we have been able to perform coupled cluster calculations [45] up to the level of CCSDT/aug-cc-pVTZ, and through extrapolation techniques we have obtained reference data at CC-cf/CBS (CC-cf=continued fraction [46] extrapolation toward full-CI limit; CBS = Complete Basis Set). Moreover, we have explored the energy profile for this reaction with 28 density functionals, (LDA, [47][48][49] [51,70] TPSSh, [71] SSB-D, [43] S12g, [44] S12h, [44] CAM-S12h, [44] CAM-B3LYP [72]), among which the most popular ones from the DFT2012 popularity poll [73] and the newly developed S12g/S12h functional [44].

Results and Discussion
The complete energy profile for the S N 2 reaction of H − +CH 4 was studied using both wavefunction and density functional methods. The reaction proceeds from the reactants (R, see Figure 1) towards a reactant complex (RC) and then crosses the central barrier (TS) to reach a product complex (PC) and finally products (P). The RC is reached early on, e.g., with a C-H(nucleophile) distance of some 3.84 Å; this is ca. 0.7 Å longer than the case of Cl − + CH 3 Cl (3.15 Å) even though the size of the nucleophile is probably much smaller here [note however that Pauling [74] reported a larger ionic radius for H − (2.08 Å) than for Cl − (1.81 Å), while Frecer [75,76] through Monte Carlo obtained values of 2.28 Å (H − ) and 2.30 Å (Cl − ) respectively]. This is consistent with a very weak ion-molecule interaction in the reactant complexes. Interestingly enough, at the highest level for which we could obtain the energies directly, CCSDT/atz (see Table 1), this leads to an energy profile of this gas-phase reaction that is more reminiscent of an S N 2 profile in solution [1,2,77] (see Figure 1). For example, the RC well is almost non-existent with a depth of 0.9 kcal·mol −1 and the barrier is quite steep with a value of 49.4 kcal·mol −1 .

Coupled Cluster Results
We extrapolated the coupled cluster energies to come as close as possible to the full-CI result, for which we use the continued-fraction approximant [46]. There are two possibilities for this extrapolation (Equations 1a,b), using either the CCSD(T) energy (as we do here in this first part), or the CCSDT energy (as discussed later on), for the ingredient. The resulting E CC-cf energies are also given in Table 1: The results of Table 1 make it clear that there is a clear basis set effect, where it is not really important to increase the basis set size but it is more important to include diffuse functions [42,78]. Given that we deal here with anionic species, for which diffuse functions are important [21,42,79], this comes as no surprise. There is also a significant electron-correlation effect, where energies obtained at the CCSDT level do not seem to have converged to the full-CI limit. For instance, the well depth increases with the atz basis from −0.  The results obtained with the different levels of coupled cluster method, i.e., CCSD, CCSD(T), CCSDT and CC-cf, together with RHF and MP2 (see Table 1) indicate that CCSD may be sufficient for getting good results. CCSD underestimates the RC well depth (e.g., by 0.3−0.6 kcal·mol −1 ), and overestimates the barrier by ca. 3−4 kcal·mol −1 . MP2 works in this respect even better with deviations from CC-cf that are twice as small, even though the computational effort is more or less the same as CCSD. As already noted often before, RHF cannot be trusted for these energy profiles as it gives barriers which are too large.

Density Functional Energies
All density functionals show the correct energy profile with a shallow well for the RC (−0.4 to −2.7 kcal·mol −1 with the largest basis set aqz), and a substantial barrier (34.2 to 50.1 kcal·mol −1 with the aqz basis). Nevertheless, the results for the 28 density functionals show quite a diversity in the accuracy for the energy profile, even though some general trends are obvious: they tend to overestimate the RC well depth, and (generally) underestimate the reaction barrier. The least reliable functional is not surprisingly LDA, which for the RC predicts a well depth of −2.74 kcal·mol −1 , and places the TS at +34.16 kcal·mol −1 (e.g., a deviation of ca. 15 kcal·mol −1 ). Early GGA functionals (PBE, BP86 [80], BLYP) improve the barrier by ca. 5-7 kcal·mol −1 , and a further 3 kcal·mol −1 is obtained by the use of OPTX in OPBE/OLYP. Surprisingly, SSB−D predicts a barrier that is ca. 2.0 kcal·mol −1 lower than OPBE, even though prior studies showed them to behave similarly. This cannot be due to the inclusion of dispersion in SSB−D, since the effect of including Grimme's dispersion energy is limited (< 0.5 kcal·mol −1 , see Table 1). Even better results are obtained with hybrid functionals and reasonable results are obtained: the difference with the CC-cf results is now only 3−4 kcal·mol −1 (at a fraction of the computational cost) for the most often used hybrid functionals (B3LYP, PBE0, M06). The recently developed S12h and M06−2X bring the deviation from CC−cf down to ca. 2 kcal·mol −1 , while excellent results (deviation <0.5 kcal·mol −1 ) are obtained with four functionals: B2PLYP (48.53 kcal·mol −1 ), M06−L (49.45 kcal·mol −1 ), B97−3 (48.46 kcal·mol −1 ) and CAM−S12h (48.41 kcal·mol −1 ). Of these four, two give an excellent description of the RC well depth: M06−L (−1.18 kcal·mol −1 ) and CAM−S12h (−1.21 kcal·mol −1 ), while the other two underestimate it slightly (B97−3, −0.44 kcal·mol −1 ) or show a non-existent RC (B2PLYP, +0.04 kcal·mol −1 ). The non−existence of a RC happens also for OPBE and OLYP with the augmented basis set (adz, atz, aqz), where the optimization proceeds towards reactants.

Structural Parameters
All methods used here confirm the early character of the RC, with a distance between the nucleophile (Nu) and the central carbon observed within the range of 2.97 Å (B2PLYP) to 4.73 Å (RHF) (see Table 2). These two values are, together with LDA (3.06 Å), rather different from the other wavefunction and density functional methods that give values roughly between 3.4 and 4.0 Å. Moreover, this C−Nu distance is the one that distinguishes the several methods for the deviations with respect to the CCSDT/atz results. The variation for the other structural parameters is much smaller (see Table 2). Table 2. Structural parameters (Å, °) for stationary points, obtained with atz basis set.

Single−Point Calculations at CCSDT/atz Geometry
The comparison of the energy profiles for the different methods is of course influenced to some extent by the different geometries used with the different methods. Therefore, and in order to make an honest comparison between the different methods we also performed single-point energy calculations using the CCSDT/atz geometries. Given in Table 3 are the results for all wavefunction and density functional methods. By comparing the results from Table 3 with those from Table 1, it can be seen that the influence of the geometry is limited. The largest difference for the different methods observed is of the order of 0.4 kcal·mol −1 (e.g., for LDA, −2.74 kcal·mol −1 for the RC at the LDA geometry, and −2.35 kcal·mol −1 at the CCSDT/atz geometry). However, the typical deviation is less than 0.1 kcal·mol −1 (i.e., chemical accuracy); for instance, PBE−D 3 shows differences of 0.06 and 0.08 kcal·mol −1 for the RC and TS energies with the two different geometries. All of this indicates that the energy surface is quite flat, as was already obvious from Figure 1.

Competition with Proton Transfer Pathway
An alternative reaction is possible in which the hydride abstracts a proton from methane. This leads to the formation of a methyl anion and H 2 : This alternative process is however associated with an endothermic reaction energy of +21.35 kcal·mol −1 at S12h/aqz. This pathway is beyond the scope of the present investigation which focuses on the thermoneutral identity S N 2 reaction.

Experimental
All wavefunction based calculations (Hartree-Fock, Second-order Møller-Plesset Perturbation Theory, Coupled Cluster Theory) have been performed with the Coupled-Cluster techniques for Computational Chemistry (CFOUR) [81,82] program version 1.2, using a variety of Dunning's correlation-consistent basis sets [83,84]: cc−pVXZ (X=D, T, Q, abbreviated as dz, tz and qz respectively) and aug−cc−pVXZ (X=D, T, Q, abbreviated as adz, atz and aqz respectively). The continued-fraction approximant [46] for obtaining coupled−cluster energies toward the full configuration-interaction limit (CC−cf) has been used with equations 1a,b to reach completeness of electron-correlation energies. The NWChem program [85] version 6.1 was used for all density functional calculations.

Conclusions
We have performed a benchmark study on the smallest bimolecular nucleophilic substitution (S N 2) reaction possible: H − +CH 4 → CH 4 +H − , for which we obtained reference data at the near-full-CI coupled cluster limit using the continued−fraction approximant (CC-cf). Unlike typical S N 2 reactions in the gas phase, which usually show a double-well potential, the current reaction shows an energy profile that resembles more the unimodal profile of the S N 2 reaction in solution, with a relatively shallow reactant-complex well of only −1.19 kcal·mol −1 and a high barrier amounting to 48.87 kcal·mol −1 . All other computational methods also clearly show the steep reaction barrier that needs to be overcome (34−50 kcal·mol −1 ), and the very shallow wells of the (ion-dipole) reactant complex. All density functionals have the tendency to underestimate the reaction barrier, while for the RC the deviations compared to CC-cf are much smaller (< 0.5 kcal·mol −1 ). Excellent results have been obtained with B97-3, M06-L and the newly developed CAM−S12h functional.