Benchmark ab initio and dynamical characterization of the stationary points of reactive atom + alkane and SN2 potential energy surfaces

ion means proton transfer from CH3I to F . PCCP Perspective


I. Introduction
Chemical thinking has been traditionally based on stationary structures of molecular systems. These stationary geometries (points) on a potential energy surface (PES) correspond to complexes (minima) and transition states (saddle points), which may play a key role in the dynamics and mechanisms of chemical reactions. Therefore, characterization of the structures and energies of these stationary points is essential to uncover reaction mechanisms, which is one of the main goals of chemistry. 1 The structure and energy of the transition state (TS) may determine the dynamics and outcome of a reaction, because reactant-like TSs (early barriers) are conveniently surmounted by faster collisions, whereas product-like late barriers prefer vibrational excitations to facilitate the reactivity. 2 Despite the key role of the stationary points of a reactive PES in chemistry, their experimental investigation is usually highly challenging or even impossible, because the complexes are often unstable and the transition states are not in equilibrium but sit at the top of an energy curve along the reaction coordinate. Some experimental insights can be obtained by advanced matrix isolation 3,4 or anion photo-electron spectroscopy (transitionstate spectroscopy) [5][6][7][8] and crossed-beam scattering, 9,10 but the complete characterization of the stationary points requires theoretical work, which may guide, explain, and complement experiments.

Balázs Olasz
Balázs Olasz received his MSc degree in pharmaceutics at the University of Szeged, Szeged, Hungary, in 2014. In 2015 he joined the Computational Reaction Dynamics Research Group and started his PhD work under the supervision of Gábor Czakó. He studied the dynamics of the F À + CH 3 I reaction using a new analytical potential energy surface. In 2018 he received a 6-month Richter Scholarship. In 2019 he defended his PhD dissertation based on 7 publications including high-profile papers in Chemical Science and Science Advances.

Dóra Papp
Dóra Papp received her PhD in theoretical chemistry at Eötvös University, Budapest, Hungary in 2017, then she joined Gábor Czakó's group as a postdoctoral researcher at the University of Szeged, Hungary. As an undergraduate student she did research in computational biochemistry to model protein aggregation and unfolding. She won the Scholarship of the Hungarian Republic, and completed a half-year-long research project at Chalmers University, Gothenburg, Sweden. During her PhD research she developed and applied a program which computes energies and lifetimes of ro-vibrational resonance states of polyatomic weaklybound molecules. Her current research interest involves PES development and chemical reaction dynamics.

István Szabó
István Szabó obtained his PhD in theoretical chemistry at Eötvös University, Budapest, Hungary in 2016. As a member of the Czakó Group he focused on PES developments and reaction dynamics simulations of S N 2 reactions.
Subsequently, he joined Edina Rosta's group at King's College London, UK to gain expertise in QM/MM and Markov model-based analysis and enhanced sampling techniques. Supported by the Cavendish Laboratory at the University of Cambridge, UK he revealed key binding properties of cucurbit[n]uril ''cages'' for selective drug detection. Currently, István is working as senior cheminformatician at ChemPass Ltd on the development of an AI-driven drug discovery platform.

Viktor Tajti
Viktor Tajti received his BSc in molecular bionics engineering and MSc in info-bionics engineering at the University of Szeged, Szeged, Hungary, in 2017 and 2019, respectively. In 2018 he took 2nd place at the University Competition of Research Students. He is currently a first-year PhD student at the University of Szeged in the group of Gábor Czakó. He has implemented several computer codes used in the group and continues his undergraduate research on the dynamics of the F À + CH 3 CH 2 Cl reaction.

Perspective PCCP
Our group has investigated the dynamics and mechanisms and/or characterized the stationary points of several polyatomic reactive systems such as reactions of atoms (F, O, Cl, Br, and/or I) with methane 11 and ethane 12 as well as ions (F À , Cl À , Br À , I À , OH À , SH À , CN À , NH 2 À , PH 2 À ) with methyl-, 13,14 ethyl-, 15 and/or amino-halides. 16 The story of our dynamics studies began with the atom + methane reactions, which have been investigated theoretically since the 90s by Espinosa-García and co-workers 17,18 based on semi-empirical PESs. In the 2000s a few ab initio-based analytical PESs were developed for the F and Cl + CH 4 reactions using different fitting strategies. [19][20][21][22] We reported our first analytical PES for the F + CH 4 39,40 and angular dependence of a TS barrier height. 10 In the case of S N 2 reactions one usually finds direct dynamics studies in the literature, where especially Hase and co-workers 41,42 have remarkable achievements. In 2013 we reported 43 the first high-level ab initio analytical PES for a S N 2 reaction (F À + CH 3 Cl) and investigated its dynamics with the quasiclassical trajectory (QCT) method. Later we also developed analytical PESs for the F À + CH 3 F and CH 3 I reactions and several other PESs are currently under development in our group. 44,45 These analytical PESs played a key role in the discovery of the double-inversion mechanism, 46,47 characterization of front-side complex formation, 48 quantum dynamics computations, 49 and comparisons with experiments. 50,51 The new findings motivated other theoretical groups; thus, Hase and co-workers 52 showed that double inversion is a non-intrinsic-reaction-coordinate pathway and Wang and coworkers 53,54 identified this mechanism in aqueous solutions. The theoretical study of chemical reactions begins with the ab initio characterization of the stationary points, which guides the full-dimensional analytical PES developments and the QCT 55 and/or quantum dynamics [56][57][58][59][60] simulations. We aim to provide the best technically feasible structures and relative energies of the stationary points utilizing sophisticated composite ab initio approaches. The composite electronic structure techniques combine different methods and basis sets to compute the most accurate results with affordable computational time. Our work uses the ideas of the focal-point analysis (FPA) approach, 61,62 which, unlike the black-box type ab initio thermochemistry protocols such as CBS-n, 63 Gn, 64 Wn, 65 HEAT, 66 etc., does not prescribe which specific methods and bases have to be used, but FPA suggests a 2-dimensional extrapolation scheme over methods (HF, MP2, CCSD, . . .) and basis sets (DZ, TZ, QZ, . . .) augmented with auxiliary corrections such as core correlation, relativistic effects, diagonal Born-Oppenheimer correction, 67 and zero-point vibrational energy. In our benchmark studies we combine the ideas of the FPA approach and the benefits of the novel explicitlycorrelated F12 correlation methods 68,69 to obtain the most accurate structures and relative energies of the stationary points characterizing PESs of chemical reactions. In Section II we describe the details of this benchmark composite ab initio approach highlighting examples from our own work. In Section III we briefly provide insight into the mechanisms of several atom/ion + molecule reactions revealed by the stationarypoint properties. In the following sections we address three questions which are rarely investigated. (1) The performance of ab initio methods and basis sets are usually tested at stationary points and/or along potential energy curves of diatomic molecules, whereas global PESs cover configurations far from the stationary geometries. Therefore, in Section IV we investigate the accuracy of different ab initio levels of theory at non-stationary geometries, thereby guiding PES developments. 70 (2) In Section V we show how the choice of the electronic structure theory affects the dynamics of a chemical reaction. 71 (3) Finally, in Section VI we review our numerical methods to uncover the role of the stationary points in the dynamics. 48,72 Our perspectives end with summary and conclusions in Section VII.

A. Structures
We usually compute the benchmark stationary-point structures using the explicitly-correlated CCSD(T)-F12b method 69 with the aug-cc-pVTZ basis set. 73 For Br and I small-core relativistic effective core potentials with the corresponding pseudopotential aug-cc-pVTZ-PP basis sets 74 are employed. Nowadays, CCSD(T)-F12b/aug-cc-pVTZ geometries can be obtained using the MOLPRO program package 75 for systems as large as F À + CH 3 CH 2 Cl as we reported 15   and post-reaction (POSTMIN) complexes and the Walden-inversion TS of the Cl À + CH 3 I reaction obtained by CCSD(T)-F12b as well as traditional CCSD(T) with the aug-cc-pVnZ (n = D, T, Q) basis sets. 76 As seen, in most cases, especially for the intramolecular distances, the CCSD(T)-F12b/aug-cc-pVTZ bond lengths agree with the QZ results within 0.001 Å, whereas the traditional CCSD(T) method gives much larger uncertainty. For example, the C-I distance at PREMIN is 2.

B. Approaching the CCSD(T)/CBS limit
In order to determine the best technically feasible relative energies of the stationary points, we perform single-point energy computations at the most accurate, usually CCSD(T)-F12b/aug-cc-pVTZ, geometries. Our first goal is to approach the CCSD(T)/complete-basis-set (CBS) limit. The traditional FPA route 61,62 is to perform HF, 77 MP2, 78 CCSD, 79 and CCSD(T) 80 computations with the aug-cc-pVnZ, where n = D, T, Q,. . . basis sets, and extrapolate the HF energy and the correlation energy increments to the CBS limits. For extrapolation usually 2-parameter asymptotic formulae 81,82 are employed which give the best CBS estimates if the largest basis-set results are used. Note that for HF extrapolation traditionally a 3-parameter expression 83 was used, but recent benchmark studies 84 showed that a 2-parameter formula 81 provides slightly better CBS results. The energies obtained by different methods and basis sets are collected into a table, whose focal point is the ''best method''/CBS result. The lower-level computations involved in the FPA tables help to estimate the uncertainty of the final result, which is an important and useful feature of the FPA analysis. The convergence of the HF, MP2, CCSD, and CCSD(T) relative energies of the stationary points of Fig. 1 with respect to the n = D, T, Q, and 5 basis sets, and their extrapolated CBS limits are shown in Fig. 2. In all cases the electron correlation effects are significant, as the HF method provides errors of about 3-5 kcal mol À1 relative to the CCSD(T) energies. In most cases MP2 outperforms the CCSD method, but for the TS the MP2 results still have differences larger than 1 kcal mol À1 relative to CCSD(T). The HF method usually approaches its CBS limit fast, because HF converges exponentially with respect to n. 85 For PREMIN and WaldenTS the basis set dependence of the correlation methods is also not significant, whereas the depth of the POSTMIN well, relative to the reactants, increases by about 2 kcal mol À1 as n goes from 2 to 5. The CCSD(T)/5Z results approach the CCSD(T)/CBS limits within 0.5 kcal mol À1 in all cases. Fig . 3 compares the basis-set convergence of the standard CCSD(T) and the explicitly-correlated CCSD(T)-F12b methods for all the stationary points of the Cl À + CH 3 I system. As seen, the standard CCSD(T) method gives large errors of 2-3 kcal mol À1 for POSTMIN and I À + CH 3 Cl products when DZ and TZ bases are used, whereas CCSD(T)-F12b is converged within 1 kcal mol À1 . Furthermore, the CCSD(T)-F12b/aug-cc-pVQZ level usually agrees with the CBS limit within about 0.2 kcal mol À1 and the agreement is never worse than 0.5 kcal mol À1 , unlike in the CCSD(T) case. Thus, the explicitly-correlated CCSD(T)-F12b/aug-cc-pVQZ computations may replace the traditional CCSD(T) CBS extrapolations for chemically accurate (with uncertainty less than 1 kcal mol À1 ) benchmark energy determinations. More detailed and critical comparisons of the convergence of the standard and F12 correlation methods can be found in ref. [86][87][88]. We note that we usually use the F12b variant of the explicitly-correlated CCSD(T) method, because previous studies 89,90 showed that CCSD(T)-F12b has more monotonic basis set convergence behavior than the CCSD(T)-F12a method, though the two methods give very similar results. 70

C. Core and post-CCSD(T) correlations
The usual frozen-core electron correlation computations correlate the valence electrons only. The core-core and core-valence correlation effects can be taken into account by computing the difference between all-electron and frozen-core energies obtained by using the same basis set. We often determine the core-correlation effects using CCSD(T) or CCSD(T)-F12b with the aug-cc-pwCVTZ or cc-pCVTZ-F12 basis sets, 91,92 respectively.
Electron correlation contributions beyond the gold-standard CCSD(T) can be computed by the MRCC program 93 via the CCSDT, 94

D. Scalar and spin-orbit relativistic effects
Scalar relativistic effects are usually smaller than the core correlation corrections as shown in Section IV. For Br and I the scalar relativistic effects are approximated by the effective core potentials (ECPs). 74 For lighter atoms we usually neglect this small effect or perform all-electron relativistic computations using the second-order Douglas-Kroll (DK) 98 Hamiltonian. For the Br + CH 4 system we showed that the ECP and the DK computations provide similar results. 28 Spin-orbit (SO) coupling may be substantial for some openshell atoms and radicals, which can be computed using the Breit-Pauli operator in the interacting-states approach. 99 The different electronic states needed to set up the SO matrix can be obtained by multi-configurational self-consistent field (MCSCF) 100 or multi-reference configuration interaction (MRCI) 101 methods. In our studies we investigated the reactions of water, 102 methane, 11 and ethane 12 with halogen atoms, where SO effects can be significant in the entrance channel. We found that MRCI+Q/ aug-cc-pVDZ computations, where +Q denotes the Davidson correction, 103 with a minimal active space provide reasonably accurate SO corrections. 102 In a relativistic computation the ground electronic state of the halogen atom ( 2 P) is split into a SO ground ( 2 P 3/2 ) and a SO ( 2 P 1/2 ) excited state. As the halogen atom approaches a molecule the 2 P 3/2 state splits into a reactive

PCCP Perspective
This journal is © the Owner Societies 2020 The other excited SO state (SO 3 ) correlating to 2 P 1/2 also approaches the non-reactive non-SO excited state. The SO coupling significantly affects the depths of the van der Waals wells in the entrance channel, but the effects are quenching or partially quenching at the other stationary points, depending on the SO value and the proximity of the stationary point to the reactants. As also shown in Fig. 5, the largest effects are found for the TSs of the I + C 2 H 6 reaction, usually 12-20%, because SO coupling is the strongest for I among halogens, and for the reactant-like hydrogen-abstraction TS of the F + C 2 H 6 system, because here the electronic structure of the F atom is only slightly perturbed. 12 For the products and productlike minima the SO effects are negligible.

E. Zero-point energy corrections
In order to compute relative energies which are comparable with experiment the zero-point energy (ZPE) corrections have to be determined. We usually perform harmonic frequency computations with the CCSD(T)-F12b method using double-or triple-zeta basis sets depending on the size of the system. Harmonic ZPE corrections are shown for the stationary points of the X + C 2 H 6 [X = F, Cl, Br, I] reactions in Fig. 6. As seen, the ZPE corrections are substantial, in the range of 2-6 kcal mol À1 for most cases, especially for non-reactant-like structures. Therefore, it is clear that the ZPE effects have to be considered if chemical accuracy is desired. Anharmonicity may cause about 5% uncertainty, which is usually less than 0.1-0.3 kcal mol À1 . If the experiments are performed at non-zero temperature, e.g. 298 K, thermal corrections also have to be calculated considering temperaturedependent electronic, translational, vibrational, and rotational enthalpy changes, as, for example, we did for the F + CH 4 -HF + CH 3 reaction in ref. 23.

F. Composite energy
Once the above-described energies and their auxiliary corrections are computed at the benchmark structures (see II. A), we can determine the high-accuracy relative energies of the stationary points as

G. Comparison with experiment
As mentioned in the Introduction, experimental determination of the relative energies of the stationary points is usually not feasible, except for the products (reaction enthalpies). Thus, the computed ZPE-corrected (adiabatic) relative energies can be compared with 0 K reaction enthalpies. In Table 1 we collected our computed benchmark adiabatic reaction energies and the available experimental results obtained from the 0 K heat of formation data of the Active Thermochemical Tables (ATcT). 104,105 As seen, for the 31 different atom/ion + molecule reactions, theory agrees with experiment with a root-mean-square deviation of only 0.34 kcal mol À1 , and the largest discrepancy (0.86 kcal mol À1 ), where experiment has a substantial uncertainty of AE0.48 kcal mol À1 , is still below 1 kcal mol À1 . Thus, this comparison demonstrates that modern ab initio theory is capable to reproduce experiment within chemical accuracy, thereby confirming the accuracy of the theoretical predictions of the experimentally not available chemical properties. Note, that if static electron correlation is significant for a stationary point, which is not the case for the systems considered in the present study, multi-reference methods should be used to achieve high accuracy.

III. Applications to atom/ion-molecule reactions
We applied the above-described benchmark ab initio composite techniques to characterize the stationary points of several atom + methane/ethane 11,12 and ion + methyl/ethyl-halide [13][14][15] reactions. As examples, Fig. 7 (Fig. 5), whereas in the case of S N 2 reactions the iondipole complexes can be below the reactants by 18 kcal mol À1 , as shown in Fig. 8. For X + C 2 H 6 substitution is a higher-energy channel via barriers between 20 and 80 kcal mol À1 depending on X and the leaving group. Besides H substitution, found for X + CH 4 reactions as well, 106,107 for X + C 2 H 6 CH 3 substitution can also occur. 12 Both product channels can be obtained via the usual Walden-inversion mechanism or through a front-side attack retention pathway as shown in Fig. 7. Front-side attack TS is also found for S N 2 reactions, with a classical barrier of 30.0 kcal mol À1 for F À + CH 3 CH 2 Cl, whereas the Walden-inversion TS is submerged by 11.3 kcal mol À1 . For S N 2 reactions our dynamics simulations revealed a new double-inversion retention pathway, 46 initiated by a proton-abstraction induced inversion followed by a second inversion via the Walden-inversion TS. For the first step of the double-inversion process we found a TS, which, in the case of F À , OH À , and NH 2 À nucleophiles, is below the front-side attack TS, 14,97 thereby opening the lowest energy retention pathway for several S N 2 reactions. (For F À + CH 3 CH 2 Cl the double-inversion   104 of the Active Thermochemical Tables (ATcT). 105 Uncertainties are derived from the uncertainties of each 0 K enthalpy of formation given in ATcT using the Gaussian error-propagation law. c Absolute energy differences (in kcal mol À1 ) between theory and experiment. d Obtained from the non-SO benchmark classical energy of 5.32 kcal mol À1 (ref. 27), D ZPE (CCSD(T)-F12b/aug-cc-pVTZ) of À4.08 kcal mol À1 , and D SO of +0.02 kcal mol À1 . e Obtained from the benchmark classical energy of À18.07 kcal mol À1 (ref. 15) and a corrected D ZPE of À4.06 kcal mol À1 .

PCCP Perspective
This journal is © the Owner Societies 2020 Phys. Chem. Chem. Phys., 2020, 22, 4298--4312 | 4305 classical barrier height is 20.7 kcal mol À1 , as seen in Fig. 8.) Unlike for the S N 2 reactions of methyl-halides, for ethyl-halide systems, bimolecular elimination (E2) leading to, for example, Cl À + HF + C 2 H 4 (Fig. 8) can also occur via syn and anti pathways in competition with the S N 2 channel. 15,108 As seen, the determination of the stationary points of reactive PESs provides a good picture about the possible reaction channels and pathways and their energetic requirements. However, to reveal the importance of the different mechanisms reaction dynamics simulations are necessary on full-dimensional PESs. We have developed such PESs and investigated the dynamics of several atom/ion + molecule reactions. 11,13 Detailed review of these dynamics studies can be found in ref. 11 and 13, here we just highlight the most important features of our work in context of the literature.
The key of our dynamics studies is that we represent the PESs by analytical functions obtained by fitting high-level ab initio energy points. 11,13,25 We can construct these PESs using a few tens of thousands of energies instead of billions of on-the-fly gradients needed for a direct dynamics study. Thus, the analytical PESs allow efficient and accurate dynamical investigations using either the QCT or quantum methods. 11,13 For atom + methane reactions such PESs were developed 17-22 before our benchmark work 23 published first for the F + CH 4 system in 2009. The unique feature of our atom + methane PESs was that we proposed 23,26-28 several composite ab initio methods to compute accurate energy points within affordable computational time based on the ideas described in Sec. II. Furthermore, we reported SO-corrected fully ab initio PESs for the first time for the F/Cl/Br + methane reactions. 26,28,109 Since then several groups have followed our ideas and developed SOcorrected PESs for reactions of halogen atoms with different molecules. 36,[110][111][112][113] In the 2010s the use of the explicitly-correlated CCSD(T)-F12 methods has become widespread for PES developments, 36,44,45,[114][115][116] which may diminish the significance of the traditional (non-F12-based) composite methods.
For S N 2 reactions our group has played a pioneering role in developing analytical PESs. 13 As mentioned in the Introduction we developed 43 the first full-dimensional high-level ab initio analytical PES for the F À + CH 3 Cl S N 2 reaction in 2013 and later we reported 44,45 analytical PESs for other S N 2 reactions as well. Unlike the traditional direct dynamics studies, 9,41,42,52,108,[122][123][124][125][126] the analytical PESs made the computations of millions of trajectories possible, allowing the discovery of low-probability reaction channels and determination of statistically accurate differential cross sections. Therefore, the analytical PESs played a key role in revealing a new reaction mechanism, called double inversion, 46 for S N 2 reactions and achieving unprecedented agreement between theory and detailed crossed-beam experiments. 50,51 As mentioned above more details about the dynamics can be found in ref. 11 and 13, here we discuss the effects of the choice of the electronic structure theory on the PES development 70 and dynamics 71 in Sections IV and V, respectively. Furthermore, we show the role of the stationary points in the dynamics of a S N 2 reaction 72 in Section VI.

IV. On the choice of the ab initio level for PES developments
In Section II we discussed the accuracy of the different ab initio levels of theory for stationary-point properties. However, global  (Note that the same finding is found for the other systems as well. 70 ) Increasing the basis to aug-cc-pVTZ the RMS drops to 1 kcal mol À1 . If we use either CCSD(T)-F12a or CCSD(T)-F12b, which virtually give the same results, the RMS becomes less than 1 kcal mol À1 even with the aug-cc-pVDZ basis. Therefore, the CCSD(T)-F12 methods are highly recommended for PES developments. Obviously, higherlevel electronic structure calculations are directly related with higher computational cost. If tens of thousands of points are necessary to describe the PES in polyatomic systems, the computational effort could be very expensive, thus the careful choice of the ab initio method and basis is necessary. Fig. 10 shows the total correlation, the core correlation, and the scalar relativistic effects for the above-mentioned six reactions. As seen, electron correlation results in energy effects as large as 5 to 20 kcal mol À1 in relative energies, showing that the Hartree-Fock method is an unreasonable choice for PES developments and direct dynamics studies. Core correlation and scalar relativity affect the PESs by 0.2-0.4 and B0.1 kcal mol À1 , respectively; therefore, these effects may be considered in spectroscopic studies, but may be negligible in PES developments for reaction dynamics computations.

V. Effects of the level of electronic structure theory on the reaction dynamics
Many studies investigated the accuracy of the various electronic structure methods for energy computations; however, little is known about their effects on the dynamics of chemical reactions. In 2018 we developed 71 20 different PESs for the F À + CH 3 I reaction using several ab initio (HF, MP2, MP2-F12, CCSD, CCSD-F12b, CCSD(T), CCSD(T)-F12b, OQVCCD(T) 117 ) and density functional theory (DFT) (B97-1, 118 PBE0, 119 M06-2X, 120 B2PLYP 121 ) methods with double-and/or triple-zeta basis sets. Then, quasiclassical trajectory computations were performed on these PESs and the effects of the level of electronic structure theory on the cross sections, reaction probabilities, angular and product internal energy distributions were revealed.
The cross sections of the S N 2 (I À + CH 3 F) and protonabstraction (HF + CH 2 I À ) channels obtained at the different levels of theory are shown in Fig. 11. As seen, the reactivity of both channels significantly depends on the level of electronic structure theory. The MP2 S N 2 cross sections are about 50-80% of the HF value, depending on the basis set, whereas for the Fig. 9 Potential energy diagram and RMS errors of different standard and explicitly-correlated (F12) frozen-core (FC) and all-electron (AE) ab initio levels of theory for the F À + CH 3 Cl S N 2 reaction. 70 The RMS errors are based on 15 energy points, obtained by varying R CX , R CY , R CH , and y HCY covering energies as indicated along the relative energy axis, and are relative to all-electron CCSD(T)-F12b/cc-pCVQZ-F12 reference data. Fig. 10 Electron correlation, core electron correlation, and scalar relativistic effects obtained as RMS deviations of 15 Hartree-Fock/aug-cc-pCVQZ, frozen-core CCSD(T)/aug-cc-pCVQZ, and Douglas-Kroll all-electron CCSD(T)/aug-cc-pCVQZ energies, respectively, relative to all-electron CCSD(T)/ aug-cc-pCVQZ energy points for six different benchmark reactions. 70

PCCP Perspective
This journal is © the Owner Societies 2020 Phys. Chem. Chem. Phys., 2020, 22, 4298--4312 | 4307 abstraction channel MP2 gives 2-3 times larger reactivity than HF without dispersion correction. In the case of HF and standard MP2 methods, increasing the basis size from DZ to TZ decreases the S N 2 reactivity by about 10-30% and increases the abstraction probability by about 100%. The MP2-F12 method gives TZ quality results with a DZ basis and in the case of the explicitly-correlated methods the DZ basis provides basis-set-converged cross sections as the MP2/TZ, MP2-F12/DZ, and MP2-F12/TZ results are virtually the same for both channels (Fig. 11). The CCSD(T) and CCSD(T)-F12b methods increase the corresponding MP2 cross sections by 20-40% and about 100% for the S N 2 and abstraction channels, respectively. In the case of the S N 2 channel the DFT functionals significantly overestimate the most accurate CCSD(T)-F12b cross sections. The B97-1 functional gives twice as large reactivity, whereas M06-2X overestimates CCSD(T)-F12b by about 30%. However, for the abstraction channel, B97-1 agrees well with the CCSD(T)-F12b result, but M06-2X gives twice larger reactivity. Product internal energy distributions for the S N 2 channel obtained on the different PESs are shown in Fig. 12, allowing comparison with the experimental results 122 of the Wester group. At low collision energy (7.4 kcal mol À1 ) the reaction produces internally hot CH 3 F molecules and the distributions peak at the highest available energy (indicating complex-forming indirect dynamics), whereas at high collision energy (35.3 kcal mol À1 ) the distributions become much broader. The HF method significantly overestimates the product internal energy by about 20 kcal mol À1 , in accord with the overestimated exothermicity. The correlation methods capture the main experimental features; the agreement is very good at low collision energy, but at high collision energy significant differences can be observed. MP2 produces too cold, whereas DFT gives too hot internal energy distributions at collision energy of 35.3 kcal mol À1 . The best agreement between theory and experiment is seen for the OQVCCD(T) method, 117 which may perform better than CCSD(T) at multi-reference configurations. More work toward this direction would be desired in the near future.
As the above findings show one should be aware of the fact that the different choices of the electronic structure theory can significantly affect the quantitative outcomes of the reaction dynamics simulations. This conclusion is in agreement with that of Hase and co-workers 123 who compared MP2 and B97-1 direct dynamics results for the F À + CH 3 I S N 2 reaction and found that B97-1 significantly overestimates the MP2 reactivity in accord with Fig. 11. If one is to develop the first PES for a system, we recommend using at least CCSD(T)-F12a/b with a DZ basis, otherwise some of the quantitative results may have large uncertainties.

VI. Role of the stationary points in the reaction dynamics
Reaction pathways are traditionally identified by visually inspecting several classical trajectory animations. One may   12 Normalized product internal energy distributions for the F À + CH 3 I S N 2 reaction obtained on various ab initio-and DFT-based analytical PESs (see Fig. 11 for notations) at collision energies of 7.4 and 35.3 kcal mol À1 . 71 The experimental data are taken from ref. 122. observe that trajectories approach several stationary points before reaching the product region. One may also find interesting cases where reaction pathways avoid an energetically favorable minimum, as Hase and co-workers 124 found in the case of the OH À + CH 3 F S N 2 reaction. Thus, trajectory animations are usually very useful to provide a qualitative picture about the role of the stationary points in the dynamics and mechanisms of a chemical reaction. However, one cannot watch millions of trajectories; therefore, a quantitative analysis technique is needed. Recently we have developed such techniques as described briefly below.
Let us start the story with our joint experimental-theoretical study on the dynamics of the F À + CH 3 Cl S N 2 reaction. 50 Experiment revealed that the F À + CH 3 Cl reaction is more direct than F À + CH 3 I, in agreement with theory. We speculated, on one hand, that the deep F À Á Á ÁICH 3 front-side minimum plays a key role in the dynamics of the F À + CH 3 I reaction, steering the reactants into a non-reactive orientation, thereby making the reaction indirect. On the other hand, the shallow F À Á Á ÁClCH 3 complex does not divert the reactants away from the reactive F À Á Á ÁH 3 CCl minimum. To quantify this prediction, we developed a trajectory orthogonal projection (TOP) method, 48 which projects the position of F À onto one-or two-dimensional subspaces of the entrance channel. In the F À + CH 3 I case we orthogonally project the Cartesian coordinates of F À along trajectories onto the C-I axis or one of the I-C-H planes and compute the distribution of the projected positions averaged over trajectories and time. TOP revealed that F À spends significant time in the front-side complex region of the F À + CH 3 I reaction, whereas front-side complex formation is negligible in the F À + CH 3 Cl reaction. Following some pioneering work, 125 our study 48 provided the first quantitative dynamical characterization of frontside complex formation in S N 2 reactions.
For the F À + CH 3 I reaction about 15 stationary points have been found 45,126 by different electronic structure theories as shown in Fig. 13. In order to quantitatively characterize their role in the dynamics we developed a method which assigns every trajectory structure to a stationary point based on the best overlap of the geometries. 72 In our implementation the best overlap is determined using an exact Eckart-transformation method, 127,128 which has been successfully used in our group for mode-specific quasiclassical polyatomic product analysis. 43,128,129 In short, we move both the stationary-point and the actual trajectory structures into the center of mass frame and construct a pseudo-rotational matrix, 127,128 which transforms the actual configuration into the Eckart frame 127,128 corresponding to the stationary-point geometry. If we take the permutational symmetry properly into account, we obtain the best overlap between the two structures. We perform this transformation for all the stationary points, and the assignment is made by minimizing the root-mean-square distances of the Cartesian coordinates of the actual trajectory geometry and the stationary points with respect to the different stationary-point structures. Note that similar automated reaction mechanism assignment technique was also reported by Taketsugu and co-workers, 130 where the trajectory geometries are assigned to structures along intrinsic reaction coordinates using a minimum-distance method via the Kabsch algorithm. 131 Fig. 13 shows the stationary-point probability distributions for the different mechanisms of the F À + CH 3 I reaction. In all cases the trajectories spend significant time in the front-side (FSMIN) and hydrogen-bonded (HMIN) minimum wells and near the hydrogen-bonded transition state (HTS). Interestingly, the formation of the traditional ion-dipole complex (PREMIN) is negligible for this reaction. This finding does not mean that the trajectory-point probability density is not high near

PCCP Perspective
PREMIN, but the configuration space of the PREMIN-like structures is much more localized than that of the HMIN-like geometries and this is the reason of the significantly higher probability of HMIN formation/assignment. The S N 2 trajectories are usually trapped in the post-reaction ion-dipole complex (POSTMIN) region, whereas for the induced-inversion (which provides an inverted reactant) and proton-abstraction channels POSTMIN formation is not significant (the small non-zero probabilities indicate recrossing dynamics). The indirect dynamics of the F À + CH 3 I reaction is underpinned by the fact that most of the proton-abstraction stationary points participate in the S N 2 channels as well. In order to provide more insights into the dynamics of the reaction, we computed stationarypoint transition-probability matrices, 72 as shown for the S N 2 inversion channel in Fig. 13. The matrix is nearly symmetric showing many forward-backward transitions between stationary points, confirming again the indirect nature of the F À + CH 3 I reaction.
The trajectories usually enter into the HMIN and HTS regions, WALDENTS is most likely approached from HTS and PREMIN, and the S N 2 products are usually formed via POSTMIN. We also proposed to apply various distance and energy constraints into the analysis, 72 which may provide additional insights into the mechanisms of chemical reactions.

VII. Summary and conclusions
Benchmark ab initio characterization of the stationary points of reactive PESs is the first step toward understanding the dynamics and mechanisms of chemical reactions. We determine the best technically feasible relative energies of the stationary points at CCSD(T)-F12b/triple-zeta geometries as given in eqn (1). The CCSD(T)/CBS limit may be obtained by extrapolation of traditional CCSD(T)/aug-cc-pVnZ [n = 4 and 5] energies or using the explicitlycorrelated CCSD(T)-F12b method with a quadruple-zeta basis. The core (D core ) and post-CCSD(T) correlation effects are usually a few tenths of kcal mol À1 and sometimes, but not always, cancel each other. Scalar relativistic effects (D rel ) are described by effective core potentials for heavy atoms, e.g., Br and I, and often neglected for first-and second-row elements as D rel , which can be obtained by Douglas-Kroll computations, is usually less than D core . Spin-orbit corrections (D SO ) can be determined by MRCI computations using the Breit-Pauli operator in the interacting-states approach 99 and can be significant for heavy open-shell species, such as Br and I atoms, and their weakly-bound complexes. Zero-point-energy corrections (D ZPE ) can be as large as a few kcal mol À1 , thus, cannot be neglected to achieve good agreement with experiment. The CCSD(T)-F12 methods with double-or triple-zeta basis sets are recommended for D ZPE computations. Considering all the above energy terms, quantum chemistry can provide definitive relative energies with uncertainties well below 1 kcal mol À1 , as confirmed by comparisons to measured 0 K reaction enthalpies.
In the present paper we focus on single-reference coupledcluster computations, which usually give accurate stationarypoint properties as demonstrated here. However, it should be noted that CCSD(T) and CCSD(T)-F12 methods may fail to provide a good description of certain regions of the PES, especially where several coupled configurations come into play. In this case multi-reference methods such as MRCI 101 or MRCC 132-134 should be used.
The above-described benchmark composite ab initio methods were applied to several atom + alkane and ion + molecule reactions. [11][12][13][14][15] The stationary-point structures and energies guide full-dimensional analytical PES developments, which allow efficient dynamical investigations. We developed and have been developing such PESs for several reactions, 11,13 which revealed a new doubleinversion mechanism, 46 unexpected leaving-group effect, 50 and front-side complex formation for S N 2 reactions, 48 as well as in the case of atom + alkane reactions extended the validity of the Polanyi rules 26,38 and mapped the angle dependence of a transition-state barrier. 10 We tested the performance of several ab initio methods and basis sets at non-stationary geometries and concluded that the explicitly-correlated CCSD(T)-F12 methods are strongly recommended for PES developments. 70 We also showed that the results of dynamics simulations such as cross sections, reaction probabilities, etc., may depend significantly, for example, by factors of 2, on the level of electronic structure theory. 71 Converged cross sections may be obtained by using a CCSD(T)-F12 method with a double-zeta basis set.
We developed numerical analysis techniques to uncover the role of the stationary points in the dynamics of chemical reactions. Trajectory orthogonal projections, 48 Eckart-transformation-based stationary-point assignments, 72 and stationary-point transition probability matrices 72 reveal the probability distributions of the trajectory geometries in different stationary-point regions and transition probabilities between them, thereby uncovering frontside complex formation in S N 2 reactions and various reaction pathways. These numerical analysis methods complement traditional trajectory animations, and provide quantitative links between static stationary-point properties and reaction dynamics simulations. Therefore, we hope that our perspectives strengthen the connections between the fields of clamped-nuclei electronic structure theory and reaction dynamics.

Conflicts of interest
There are no conflicts of interest to declare.