On the development of a gold-standard potential energy surface for the OH (cid:2) + CH 3 I reaction †

We report a story where CCSD(T) breaks down at certain geometries of the potential energy surface (PES) of the OH (cid:2) + CH 3 I reaction. To solve this problem, we combine CCSD-F12b and Brueckner-type BCCD(T) methods to develop a full-dimensional analytical PES providing method-and basis-converged statistically-accurate S N 2 and proton-transfer cross sections.

The knowledge of the potential energy surface (PES) is essential to compute the dynamics and mechanisms of chemical reactions.2][3][4][5][6] Direct dynamics simulations determine the potential gradients on-the-fly along the trajectories using an electronic structure program package, which is very time consuming, and thus only a low level of theory can be afforded.5][16] The millions of trajectories computed on these PESs revealed a new retention mechanism for S N 2 reactions, called double inversion, 8 front-side complex formation, 12 an unexpected leaving group effect, 9 and unprecedented agreement with crossed-beam experiments. 9,17oving beyond the six-atomic systems, in the present study, we aim to develop a full-dimensional analytical ab initio PES for the OH À + CH 3 I reaction.Following the pioneering direct dynamics studies 3,5,18,19 as well as crossed-beam 3,5,20 and selected ion flow tube 18 experiments, the present PES will be the first analytical surface for the title reaction, and moreover the first PES for a 7-atomic S N 2 reaction.The efficient PES development is made possible by utilizing Robosurfer, 15 our very recently developed program package enabling automatic construction of analytical PESs for reactive systems.Our plan is to build a PES at the MP2/ aug-cc-pVDZ level of theory 21,22 and then to re-compute the energy points using the more accurate and more time-consuming CCSD(T)-F12b/aug-cc-pVTZ level of theory. 23Using this strategy, we develop a good MP2/aug-cc-pVDZ surface, which performs very well in quasiclassical trajectory (QCT) simulations.However, unexpectedly, the CCSD(T)-F12b PES gives many unphysical trajectories, for example, at zero impact parameter (b = 0) about 10% of the trajectories break up to many energetically forbidden fragments.Does the CCSD(T) method, 24 often called the gold-standard of quantum chemistry, fail to provide an accurate PES for the OH À + CH 3 I reaction?Is there a solution to this problem?Can we develop an accurate analytical PES for this system?The answer for the third question is a yes.If one is interested in knowing how, we describe the details below along with the answers to the other two questions.
We develop a full-dimensional ab initio PES for the title reaction by fitting roughly 37 000 MP2/aug-cc-pVDZ energy points using the permutationally invariant polynomial approach, 25,26 while the fitting set is generated via the Robosurfer program system. 15The computational details of the automated PES development are given in the ESI.† The PES points are recomputed using the explicitly-correlated CCSD(T)-F12b/aug-cc-pVTZ level of theory and the distributions of the differences between the CCSD(T)-F12b and MP2 relative energies are shown in Fig. 1.The root-mean-square (RMS) energy difference is 4.5 kcal mol À1 , which is typical when one compares MP2 and CCSD(T) relative energies.However, we find 17 CCSD(T)-F12b points with potential energies below the corresponding MP2 data by more than 50 kcal mol À1 .Such energy points may cause deep holes on the CCSD(T)-F12b PES, which may be responsible for the large fraction of unphysical trajectories.Indeed, upon the removal of the aforementioned 17 configurations from the dataset, the probability of the unphysical b = 0 trajectories drops from 6-13% to 2-3%.To get a deeper insight into the origin of this issue, we recompute the energy points using the DF-MP2-F12/aug-cc-pVTZ, 27 BCCD/aug-cc-pVDZ, 28 CCSD-F12b/aug-cc-pVTZ, 23 and BCCD(T)/aug-cc-pVDZ 28  of theory.Using the CCSD-F12b method, the 17 low-energy points disappear, and we do not find any configuration whose relative energy is at least 20 kcal mol À1 below the MP2 results (Fig. 1).Thus, we find that the perturbative (T) correction is too negative in some cases, and this is the reason why the CCSD(T)-F12b PES breaks down in certain regions.0][31] To solve this problem, we investigate the performance of the Brueckner coupled cluster doubles (BCCD) and BCCD with perturbative triples (BCCD(T)) methods. 28s shown in Fig. 1, unlike CCSD(T)-F12b, the BCCD(T) method does not give relative energies below the MP2 data by more than 50 kcal mol À1 .Thus, BCCD(T) seems to be a promising method to incorporate the (T) correction without making ''holes'' on the PES.Since an explicitly-correlated version of the BCCD(T) method is not widely available in quantum chemistry software and the 2.9 kcal mol À1 RMS difference between DF-MP2-F12/aug-cc-pVTZ and MP2/aug-cc-pVDZ relative energies indicates significant basis set effects, we propose a composite ab initio method defined as CCSD-F12b/aug-cc-pVTZ + BCCD(T)/aug-cc-pVDZ where the good basis-set convergence is ensured by the first term and the (T) correction is described by the Brueckner coupled cluster theory.This composite method solves the issue of the 17 low-energy outliers as seen in Fig. 1.
To further investigate the accuracy of the different levels of electronic structure theory, we select several representative geometries and compute their relative energies using various ab initio methods and basis sets.Two representative examples are shown in Fig. 2 and additional results are presented in the ESI.† On the one hand, the left panel of Fig. 2 shows a typical case where all the ab initio levels provide relative energies within a 4 kcal mol À1 window, BCCD(T) agrees with CCSD(T) within 0.1 kcal mol À1 , and CCSD(T) differs from CCSDT by only 0.01 kcal mol À1 showing the excellent performance of the perturbative (T) approximation.On the other hand, the right panel of Fig. 2 shows an extreme configuration with strong correlation effects.Using the aug-cc-pVDZ basis set, the HF method 32 provides an energy of 94.2 kcal mol À1 relative to the reactants, the corresponding MP2 result is only 61.8 kcal mol À1 , which is relatively close to the {CCSD, BCCD, OQVCCD} values of {56.5, 51.2, 51.0} kcal mol À1 , where OQVCCD denotes optimized-orbital quasi-variational coupled cluster doubles. 33However, considering the (T) corrections, striking differences are found.The {CCSD(T), BCCD(T), OQVCCD(T)} methods give relative energies of {À61.0,19.2, 14.7} kcal mol À1 , showing the breakdown of the CCSD(T) method.If we do not use the pertubative (T) approximation and perform the full CCSDT computations with single, double, and triple excitations, 34 the above relative energy becomes 29.7 kcal mol À1 , which clearly shows the issue with Fig. 1 Distributions of the number of structures as a function of relative potential energy differences from the MP2/aug-cc-pVDZ reference data.The composite energies are defined as CCSD-F12b/aug-cc-pVTZ + BCCD(T)/aug-cc-pVDZ À BCCD/aug-cc-pVDZ. the CCSD(T) method.The best agreement between CCSDT and the (T) methods is found for BCCD(T) as seen from the above data and also supported by several additional test computations reported in the ESI † (Fig. S1).This motivated us to define the composite energy as given in eqn (1), which provides a relative energy of 31.1 kcal mol À1 in good agreement with the CCSDT value (29.7 kcal mol À1 ) even in this problematic case.
To get into the physical insight, we can conclude that the failure of the traditional perturbative (T) approximation occurs at homolytic bond (e.g., C-I) dissociation, where HF orbital quasi-degeneracy and/or wrong t-amplitudes cause the breakdown of the perturbation theory. 29This problematic behavior is also related to significant multi-reference character predicted by large T 1 -diagnostic 35 values (40.03) as shown in Table S2 (ESI †).In these strongly correlated cases, the use of the non-HF Brueckner reference, which provides the best overlap with the exact wave function, may be more robust than the traditional HF-based approaches.
After the above-described ab initio investigations, we develop eight different analytical PESs by fitting MP2/aug-cc-pVDZ, DF-MP2-F12/aug-cc-pVTZ, BCCD/aug-cc-pVDZ, CCSD-F12b/ aug-cc-pVTZ, BCCD(T)/aug-cc-pVDZ, CCSD(T)-F12b/aug-cc-pVTZ, CCSD(T)-F12b/aug-cc-pVTZ 0 , and composite [eqn (1)] ab initio data, where 0 denotes the removal of the 17 outlier configurations.The classical relative energies of the stationary points corresponding to the S N 2 pathways of the OH À + CH 3 I reaction obtained on the different PESs are shown in Fig. 3. Our all-electron CCSDT(Q)/ complete-basis-set-quality benchmark data are also available for comparison. 36Electron correlation effects are significant for the product region, and both MP2 and CCSD methods differ from the benchmark data by 2-3 kcal mol À1 in the case of the CH 3 OHÁ Á ÁI À minimum and the CH 3 OH + I À products.(T) correction is needed to achieve accuracy within 0.7 kcal mol À1 .Basis set effects are important for the front-side (FS) attack and double-inversion (DI) TSs as, without explicit correlation, the DZ results differ from the TZ ones by about 1.5 (FSTS) and 2.5 (DITS) kcal mol À1 , respectively.Thus, as seen, the (T) method with a TZ-quality basis and/or explicit correlation is needed to obtain a chemically accurate PES.Indeed, the CCSD(T)-F12b/aug-cc-pVTZ 0 PES provides stationary-point relative energies with only 0.43 kcal mol À1 RMS and 0.66 kcal mol À1 maximum deviations from the benchmark data.However, as we discussed above, CCSD(T)-F12b fails at certain non-stationary regions of the PES, where only the composite method defined in eqn (1) performs satisfactorily.Therefore, it is comforting to find that the composite PES also provides excellent stationary-point energies, in agreement with the benchmark data with only 0.47 kcal mol À1 RMS and 0.62 kcal mol À1 maximum deviations.We can conclude that when CCSD(T)-F12b performs well, the present composite method provides the same energies as CCSD(T)-F12b/ aug-cc-pVTZ within 0.1 or 0.2 kcal mol À1 (Fig. 3), whereas in certain regions where CCSD(T)-F12b suffers from serious breakdown, the glory of the composite method is obvious (Fig. 2, right panel).
In total, more than 1 million trajectories are computed for the OH À + CH 3 I reaction on the eight different PESs at collision energies (E coll ) of 5, 20, and 50 kcal mol À1 .Computational details and opacity functions (Fig. S2, ESI †) are given in the ESI, † and cross sections for the S N 2 and proton-abstraction channels are shown in Fig. 4. Furthermore, Fig. 4 also shows the cross sections of the rejected unphysical trajectories (resulting in energetically nonavailable products), which allows the assessment of the quality of the PESs.The original CCSD(T)-F12b PES gives large cross sections of 13.0, 7.0, and 4.7 bohr 2 for these unphysical products at collision energies of 5, 20, and 50 kcal mol À1 , respectively.At the highest E coll , the cross-section of the unphysical trajectories is comparable to those of the S N 2 (9.2 bohr 2 ) and abstraction products (22.5 bohr 2 ).Removing 17 low-energy outliers from the CCSD(T)-F12 dataset reduces the unphysical cross sections to 1.8, 0.5, and 1.3 bohr 2 , in the above order.Fortunately, using the composite PES, these troublesome values become negligible such as 0.3, 0.1, and 0.2 bohr 2 , respectively.Despite the different unphysical cross sections, it is important and comforting to find that the reactive S N 2 and proton-abstraction cross sections are almost the same on the CCSD(T)-F12b, CCSD(T)-F12b 0 , and composite PESs as seen in Fig. 4. Considering the lower-level PESs, which can be important to assess the accuracy of direct dynamics studies, significant basis-set effects are found for the S N 2 channel, because the DZ basis overestimates the TZ cross sections by 10-100%, depending on the ab initio method and E coll .Electron correlation effects are

Fig. 2
Fig. 2 Energies of two representative structures relative to the reactants obtained by different ab initio levels of theory.

Fig. 3
Fig. 3 Classical relative energies (kcal mol À1 ) of the stationary points characterizing the OH À + CH 3 I S N 2 reaction obtained on the different analytical PESs.The composite PES energies are defined in eqn (1) and the benchmark data are taken from ref. 36.

Fig. 4
Fig. 4 Cross sections of the S N 2 (I À + CH 3 OH), proton-abstraction (H 2 O + CH 2 I À ), and rejected (unphysical) channels of the OH À + CH 3 I reaction obtained on the different analytical ab initio PESs at collision energies of 5, 20, and 50 kcal mol À1 .