On the development of a gold-standard potential energy surface for the OH⁻ + CH₃I reaction

Abstract

We report a story where CCSD(T) breaks down at certain geometries of the potential energy surface (PES) of the OH⁻ + CH₃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_N2 and proton-transfer cross sections.

---

The knowledge of the potential energy surface (PES) is essential to compute the dynamics and mechanisms of chemical reactions. For bimolecular nucleophilic substitution (S_N2) reactions, the use of the direct dynamics approach became popular due to its black-box nature and widespread applicability.^1–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. Therefore, in 2013, we started to develop analytical PESs for S_N2 reactions by fitting high-level ab initio energy points,⁷ thereby allowing efficient classical and/or quantum dynamics simulations.^7–13 We have already reported full-dimensional global PESs for the F⁻ + CH₃Y [Y = F, Cl, Br, I] reactions, which accurately describe both the S_N2 and proton-transfer channels.^8,14–16 The millions of trajectories computed on these PESs revealed a new retention mechanism for S_N2 reactions, called double inversion,⁸ front-side complex formation,¹² an unexpected leaving group effect,⁹ and unprecedented agreement with crossed-beam experiments.^9,17

Moving beyond the six-atomic systems, in the present study, we aim to develop a full-dimensional analytical ab initio PES for the OH⁻ + CH₃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¹⁸ experiments, the present PES will be the first analytical surface for the title reaction, and moreover the first PES for a 7-atomic S_N2 reaction. The efficient PES development is made possible by utilizing Robosurfer,¹⁵ 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.²³ Using 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,²⁴ often called the gold-standard of quantum chemistry, fail to provide an accurate PES for the OH⁻ + CH₃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.¹⁵ The 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⁻¹, 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⁻¹. 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,²⁷ BCCD/aug-cc-pVDZ,²⁸ CCSD-F12b/aug-cc-pVTZ,²³ and BCCD(T)/aug-cc-pVDZ²⁸ levels 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⁻¹ 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. This finding is not unprecedented as similar wrong (T) behavior was reported for potential energy curves of BH, HF, OH⁻, F₂, C₂, and BN.^29–31 To solve this problem, we investigate the performance of the Brueckner coupled cluster doubles (BCCD) and BCCD with perturbative triples (BCCD(T)) methods.²⁸

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.

As 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⁻¹. 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⁻¹ 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 − BCCD/aug-cc-pVDZ, (1)

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⁻¹ window, BCCD(T) agrees with CCSD(T) within 0.1 kcal mol⁻¹, and CCSD(T) differs from CCSDT by only 0.01 kcal mol⁻¹ 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³² provides an energy of 94.2 kcal mol⁻¹ relative to the reactants, the corresponding MP2 result is only 61.8 kcal mol⁻¹, which is relatively close to the {CCSD, BCCD, OQVCCD} values of {56.5, 51.2, 51.0} kcal mol⁻¹, where OQVCCD denotes optimized-orbital quasi-variational coupled cluster doubles.³³ However, 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⁻¹, 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,³⁴ the above relative energy becomes 29.7 kcal mol⁻¹, which clearly shows the issue with 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⁻¹ in good agreement with the CCSDT value (29.7 kcal mol⁻¹) even in this problematic case.

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

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.²⁹ This problematic behavior is also related to significant multi-reference character predicted by large T₁-diagnostic³⁵ values (>0.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′, and composite [eqn (1)] ab initio data, where ′ denotes the removal of the 17 outlier configurations. The classical relative energies of the stationary points corresponding to the S_N2 pathways of the OH⁻ + CH₃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.³⁶ Electron correlation effects are significant for the product region, and both MP2 and CCSD methods differ from the benchmark data by 2–3 kcal mol⁻¹ in the case of the CH₃OH⋯I⁻ minimum and the CH₃OH + I⁻ products. (T) correction is needed to achieve accuracy within 0.7 kcal mol⁻¹. 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⁻¹, 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′ PES provides stationary-point relative energies with only 0.43 kcal mol⁻¹ RMS and 0.66 kcal mol⁻¹ 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⁻¹ RMS and 0.62 kcal mol⁻¹ 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⁻¹ (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).

Fig. 3 Classical relative energies (kcal mol⁻¹) of the stationary points characterizing the OH⁻ + CH₃I S_N2 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.

In total, more than 1 million trajectories are computed for the OH⁻ + CH₃I reaction on the eight different PESs at collision energies (E_coll) of 5, 20, and 50 kcal mol⁻¹. Computational details and opacity functions (Fig. S2, ESI†) are given in the ESI,† and cross sections for the S_N2 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 non-available 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² for these unphysical products at collision energies of 5, 20, and 50 kcal mol⁻¹, respectively. At the highest E_coll, the cross-section of the unphysical trajectories is comparable to those of the S_N2 (9.2 bohr²) and abstraction products (22.5 bohr²). 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², in the above order. Fortunately, using the composite PES, these troublesome values become negligible such as 0.3, 0.1, and 0.2 bohr², respectively. Despite the different unphysical cross sections, it is important and comforting to find that the reactive S_N2 and proton-abstraction cross sections are almost the same on the CCSD(T)-F12b, CCSD(T)-F12b′, 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_N2 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 found to be more pronounced at higher collision energies, and the (T) correction increases the S_N2 cross sections by 50–100%. Due to the opposite effects of the basis and correlation, the error cancelation improves the accuracy of the MP2/aug-cc-pVDZ cross-section predictions. For the proton-abstraction channel, the reactivity shows less significant method- and basis-set dependence (Fig. 4). Both the S_N2 and proton-abstraction cross sections decrease with increasing E_coll, as expected in the case of exothermic, barrier-less reactions. Note that unlike in the case of the X⁻ + CH₃Y reactions,³⁷ where X and Y are halogens, for OH⁻ + CH₃I, the proton-abstraction channel is slightly exothermic.^3,18 At E_coll = 5 kcal mol⁻¹ the S_N2 and abstraction reactivity is similar, whereas at higher collision energies, the energetically less favored proton-abstraction channel clearly dominates over the S_N2 reaction, as the CH₂I⁻/I⁻ product cross section ratios, on the most accurate composite PES, are 1.4, 2.7, and 2.2 at E_coll of 5, 20, and 50 kcal mol⁻¹, respectively. This finding is in qualitative agreement with the B97-1/ECP/d direct dynamics results of the Hase group, as they reported CH₂I⁻/I⁻ ratios of 0.6, 2.3, 1.7, and 2.0, at E_coll of 1, 12, 23, and 46 kcal mol⁻¹, respectively.³ Considering the absolute cross sections, direct dynamics found S_N2 cross sections of 38.9 ± 11.1 and 15.4 ± 3.6 bohr² at E_coll = 23 and 46 kcal mol⁻¹, respectively, which may be compared to our composite results of 28.7 and 9.3 bohr² corresponding to E_coll = 20 and 50 kcal mol⁻¹, respectively. This comparison indicates that B97-1 overestimates the S_N2 reactivity. For the proton-abstraction channel the trend is less obvious, because direct dynamics gives 64.6 ± 11.4 and 31.4 ± 5.7 bohr²,³ whereas the composite PES provides 76.2 and 20.9 bohr², in order at the above E_coll values. One may notice the large error bars of the direct dynamics cross sections resulting from the relatively low total number of trajectories (3736),³ whereas the present analytical PESs make the computation of more than a million trajectories possible, thereby providing results without any significant statistical uncertainty.

Fig. 4 Cross sections of the S_N2 (I⁻ + CH₃OH), proton-abstraction (H₂O + CH₂I⁻), and rejected (unphysical) channels of the OH⁻ + CH₃I reaction obtained on the different analytical ab initio PESs at collision energies of 5, 20, and 50 kcal mol⁻¹.

The present methodology tests may provide guidance for future PES developments for similar challenging systems and the new high-level ab initio PES opens the door for detailed classical and/or quantum dynamics investigations allowing critical comparisons with experiments.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

We thank the NKFIH, K-125317, the Ministry of Human Capacities, Hungary grant 20391-3/2018/FEKUSTRAT, the Momentum (Lendület) Program of the Hungarian Academy of Sciences, and the National Young Talent Scholarship (NTP-NFTÖ-18-B-0399 for D. A. T.) for financial support.

References

1. J. Xie and W. L. Hase, Science, 2016, 352, 32.
2. J. Xie, R. Otto, J. Mikosch, J. Zhang, R. Wester and W. L. Hase, Acc. Chem. Res., 2014, 47, 2960.
3. J. Xie, R. Sun, M. R. Siebert, R. Otto, R. Wester and W. L. Hase, J. Phys. Chem. A, 2013, 117, 7162.
4. Y.-T. Ma, X. Ma, A. Li, H. Guo, L. Yang, J. Zhang and W. L. Hase, Phys. Chem. Chem. Phys., 2017, 19, 20127.
5. J. Xie, J. Zhang, R. Sun, R. Wester and W. L. Hase, Int. J. Mass Spectrom., 2019, 438, 115.
6. H. Tachikawa and M. Igarashi, Chem. Phys., 2006, 324, 639.
7. I. Szabó, A. G. Császár and G. Czakó, Chem. Sci., 2013, 4, 4362.
8. I. Szabó and G. Czakó, Nat. Commun., 2015, 6, 5972.
9. M. Stei, E. Carrascosa, M. A. Kainz, A. H. Kelkar, J. Meyer, I. Szabó, G. Czakó and R. Wester, Nat. Chem., 2016, 8, 151.
10. Y. Wang, H. Song, I. Szabó, G. Czakó, H. Guo and M. Yang, J. Phys. Chem. Lett., 2016, 7, 3322.
11. Y. Li, Y. Wang and D. Y. Wang, J. Phys. Chem. A, 2017, 121, 2773.
12. I. Szabó, B. Olasz and G. Czakó, J. Phys. Chem. Lett., 2017, 8, 2917.
13. B. Olasz and G. Czakó, Phys. Chem. Chem. Phys., 2019, 21, 1578.
14. I. Szabó, H. Telekes and G. Czakó, J. Chem. Phys., 2015, 142, 244301.
15. T. Győri and G. Czakó, J. Chem. Theory Comput., 2020, 16, 51.
16. B. Olasz, I. Szabó and G. Czakó, Chem. Sci., 2017, 8, 3164.
17. M. Stei, E. Carrascosa, A. Dörfler, J. Meyer, B. Olasz, G. Czakó, A. Li, H. Guo and R. Wester, Sci. Adv., 2018, 4, eaas9544.
18. J. Xie, S. C. Kohale, W. L. Hase, S. G. Ard, J. J. Melko, N. S. Shuman and A. A. Viggiano, J. Phys. Chem. A, 2013, 117, 14019.
19. J. Xie, J. Zhang and W. L. Hase, Int. J. Mass Spectrom., 2015, 378, 14.
20. R. Otto, J. Brox, S. Trippel, M. Stei, T. Best and R. Wester, Nat. Chem., 2012, 4, 534.
21. C. Møller and M. S. Plesset, Phys. Rev., 1934, 46, 618.
22. T. H. Dunning, Jr., J. Chem. Phys., 1989, 90, 1007.
23. T. B. Adler, G. Knizia and H.-J. Werner, J. Chem. Phys., 2007, 127, 221106.
24. K. Raghavachari, G. W. Trucks, J. A. Pople and M. Head-Gordon, Chem. Phys. Lett., 1989, 157, 479.
25. B. J. Braams and J. M. Bowman, Int. Rev. Phys. Chem., 2009, 28, 577.
26. J. M. Bowman, G. Czakó and B. Fu, Phys. Chem. Chem. Phys., 2011, 13, 8094.
27. H.-J. Werner, T. B. Adler and F. R. Manby, J. Chem. Phys., 2007, 126, 164102.
28. K. A. Brueckner, Phys. Rev., 1954, 96, 508.
29. J. B. Robinson and P. J. Knowles, J. Chem. Phys., 2013, 138, 074104.
30. D. A. Mazziotti, Phys. Rev. Lett., 2008, 101, 253002.
31. D. Kats, D. Kreplin, H.-J. Werner and F. R. Manby, J. Chem. Phys., 2015, 142, 064111.
32. W. J. Hehre, L. Radom, P. V. R. Schleyer and J. A. Pople, Molecular Orbital Theory, Wiley, New York, 1986.
33. J. B. Robinson and P. J. Knowles, Phys. Chem. Chem. Phys., 2012, 14, 6729.
34. J. Noga and R. J. Bartlett, J. Chem. Phys., 1987, 86, 7041.
35. T. J. Lee and P. R. Taylor, Int. J. Quantum Chem., 1989, S23, 199.
36. D. A. Tasi, Z. Fábián and G. Czakó, J. Phys. Chem. A, 2018, 122, 5773.
37. I. Szabó and G. Czakó, J. Phys. Chem. A, 2015, 119, 3134.

---

Footnote

† Electronic supplementary information (ESI) available: Additional computational details, discussion, and ab initio test results as well as opacity functions. See DOI: 10.1039/c9cp07007a

---