Pericyclic reaction benchmarks: hierarchical computations targeting CCSDT(Q)/CBS and analysis of DFT performance

Hierarchical, convergent ab initio benchmark computations were performed followed by a systematic analysis of DFT performance for five pericyclic reactions comprising Diels-Alder, 1,3-dipolar cycloaddition, electrocyclic rearrangement, sigmatropic rearrangement, and double group transfer prototypes. Focal point analyses (FPA) extrapolating to the ab initio limit were executed via explicit quantum chemical computations with electron correlation treatments through CCSDT(Q) and correlation-consistent Gaussian basis sets up to aug′-cc-pV5Z. Optimized geometric structures and vibrational frequencies of all stationary points were obtained at the CCSD(T)/cc-pVTZ level of theory. The FPA reaction barriers and energies exhibit convergence to within a few tenths of a kcal mol−1. The FPA benchmarks were used to evaluate the performance of 60 density functionals (eight dispersion-corrected), covering the local-density approximation (LDA), generalized gradient approximations (GGAs), meta-GGAs, hybrids, meta-hybrids, double-hybrids, and range-separated hybrids. The meta-hybrid M06-2X functional provided the best overall performance [mean absolute error (MAE) of 1.1 kcal mol−1] followed closely by the double-hybrids B2K-PLYP, mPW2K-PLYP, and revDSD-PBEP86 [MAE of 1.4–1.5 kcal mol−1]. The regularly used GGA functional BP86 gave a higher MAE of 5.8 kcal mol−1, but it qualitatively described the trends in reaction barriers and energies. Importantly, we established that accurate yet efficient meta-hybrid or double-hybrid DFT potential energy surfaces can be acquired based on geometries from the computationally efficient and robust BP86/DZP level.


Introduction
Owing to their robust and versatile nature, pericyclic reactions are highly useful chemical transformations in organic and organometallic chemistry. 1 A particular advantage of this class of reactions is the ability to increase molecular complexity in a single reaction step, often with high stereoselectivity. For this reason, these reactions have been widely utilized in syntheses of a broad range of target compounds, including complex natural products, as well as molecular species with applications in material sciences and medicinal chemistry. 2 An important feature of pericyclic reactions, including the herein studied Diels-Alder (DA) reaction, 1,3-dipolar cycloaddition (1,3-DC), electrocyclic rearrangement (ER), sigmatropic rearrangement (SR), and double group transfer (DGT) (Scheme 1), is that they proceed via a fully conjugated cyclic transition state. 1 Density functional theory computations have been paramount for elucidating the nature of such concerted chemical transformations. 3,4 However, to obtain accurate results an appropriate density functional must be chosen among myriad possibilities, a decision that is hardly transparent and one that requires highly accurate reference data for calibration.
In recent decades, various studies have been dedicated to finding the best density functional approximation for either a large set of various pericyclic reactions 5 or only a subset of them, such as Diels-Alder 6 or 1,3-dipolar cycloaddition reactions. 7 Most of these studies benchmarked density functionals by means of a 'single-shot' energy calculation with a composite method of some complete basis set (CBS) type; however, the true quality of the reference data is unknown. In fact, Karton and Goerigk recently discovered that the oftenused CBS-QB3 method fails to give accurate reaction barriers for four classes of pericyclic reactions (26 reactions in total), 5c exhibiting maximum and RMS errors of À4.7 and 2.5 kcal mol À1 , respectively, relative to the W2-F12 protocol. This demonstrates the necessity of having indubitable reference data for the potential energy surfaces (PESs) of pericyclic reactions. Another study by Karton in 2019 8 computed the barrier heights of a diverse set of 28 reactions that proceed via ring-forming transition states using the W3lite-F12 protocol. As detailed below, part of their investigation overlaps with the DA, ER, and SR reactions studied here; however, our results are more extensive and rigorous, pushing further toward the limits of electronic structure theory.
The first goal of this study is to employ high-level ab initio methods to compute definitive structures and energetics for the paradigmatic pericyclic reactions depicted in Scheme 1. We employ the highest-level correlated wave functions [CCSD(T)/ cc-pVTZ] used to date to optimize geometric structures for these reactions. A majority of previous studies on pericyclic reactions have utilized DFT geometries. 5a,c,6,7a,b Nevertheless, a few ab initio structures can be found, computed with, for example, CCSD(T)/ 6-31G (d,p) for DGT, 9a CCSD/6-311++G(df,p) for ER, 9b and MP2/ 6-31G(d,p) for 1,3-DC, 9c but these are of significantly lower quality than the geometrical structures presented here. To pinpoint reaction energetics, hierarchical focal point analyses (FPA) 10 are executed to converge toward both the one-and n-particle limits of ab initio quantum chemistry by means of correlation-consistent Gaussian basis sets [aug 0 -cc-pVXZ, X = D, T, Q, 5] of systematically increasing flexibility and polarization conjoined with the highorder coupled-cluster electron correlation series [HF, MP2, CCSD, CCSD(T), CCSDT, CCSDT(Q)].
The second goal of this study is to employ the FPA benchmarks to evaluate the performance of 60 diverse density functionals, including eight with dispersion corrections, implemented with Slater-type TZ2P and QZ4P basis sets. The functionals range in quality from the local-density approximation (LDA) to generalized gradient approximations (GGAs), meta-GGAs, hybrids, metahybrids, double-hybrids, and range-separated hybrids. An efficient and accurate protocol is established for the study of pericyclic reactions that relies on the computation of BP86/DZP optimum geometries for final TZ2P single-point energy computations with the meta-hybrid functional M06-2X, or the alternatives B2K-PLYP, mPW2K-PLYP, and revDSD-PBEP86.

Benchmark geometric structures, vibrational frequencies, and relative energies
Fully optimized geometric structures and harmonic vibrational frequencies of the reactants, products, and transition states on the potential energy surfaces of the pericyclic reactions were computed using coupled-cluster electronic wave functions with single, double, and perturbative triple excitations [CCSD(T)] 11 in conjunction with the correlation-consistent polarized valence basis set cc-pVTZ. 12 The character of each stationary point was verified by its vibrational frequencies, and the normal mode associated with the single imaginary frequency of the transition states was inspected to ensure that it pertained to the reaction of interest. The CCSD(T)/cc-pVTZ harmonic vibrational frequencies also provide high-quality zero-point vibrational energies (ZPVEs) and thus are provided in Table S1 (ESI †).
To obtain final energetics, focal point analyses 11 were executed on the CCSD(T)/cc-pVTZ geometric structures using the following hierarchy of quantum chemical methods: HF, MP2, 13 CCSD, 14 CCSD(T), 11 CCSDT, 15 and CCSDT(Q) 16 in conjunction with the cc-pVDZ 9 and aug 0 -cc-pVXZ (X = D, T, Q, 5) basis sets to Scheme 1 Paradigmatic pericyclic reactions studied in this work. systematically proceed toward the CBS limit. These computations were carried out using the Molpro 2010 package, 17 except that Kállay's MRCC 18 program was employed in the CCSDT and CCSDT(Q) cases. The aug 0 -cc-pVXZ basis sets are composed of aug-cc-pVXZ 12b functions on the heavy atoms and standard Dunning sets (cc-pVXZ) 12a on the hydrogens. The following extrapolation equations were utilized for the Hartree-Fock 19 (E HF ) and correlation energies 20 (e): where X is the cardinal number of the aug 0 -cc-pVXZ basis sets [X = (3,4,5) for HF; (4,5) for MP2; and (3,4) for CCSD and CCSD(T)]. Ancillary FPA reaction energies and barrier heights obtained with the non-augmented cc-pVXZ series are presented in Table S2 (ESI †); these results differ from the more tightly converged aug 0 -cc-pVXZ predictions by no more than 0.17 kcal mol À1 , in accord with previous findings. 21 Increments for the higher-order CCSDT and CCSDT(Q) levels of theory were obtained using the following additivity schemes: Numerous studies have demonstrated the effectiveness of such CCSDT and CCSDT(Q) corrections computed with the double-z basis set, a consequence of the general basis set insensitivity of high-order correlation increments. 22 Core-electron correlation effects were quantified by subtracting all-electron (AE) and frozen-core (FC) CCSD(T) energies computed with the cc-pCVTZ basis set: 12c A first-order relativistic correction, 23 D(rel), was obtained from the one-electron mass-velocity and Darwin terms at the CCSD(T)/cc-pVTZ level of theory. The diagonal Born-Oppenheimer correction (DBOC) 24 was computed at the HF/aug-cc-pVTZ level of theory. The D(rel) and D(DBOC) calculations were performed with the CFOUR 2.0 package. 25 In all cases the auxiliary D corrections are less than 0.25 kcal mol À1 in magnitude, confirming the use of frozen-core, nonrelativistic, clampednucleus electronic wave functions in the primary FPA determinations. In particular, the minuscule DBOC shifts (o0.08 kcal mol À1 ) indicate that our pericyclic reactions are not complicated by surface crossings near the transition states.

Density functional theory geometric structures and relative energies
The key stationary points on the potential energy surfaces of the pericyclic reactions were optimized with the generalized gradient approximation BP86, 26 which has been shown to give accurate geometric structures 27 and is frequently employed for pericyclic reactions. 4a,b,7d The BP86 functional was implemented with a hierarchical series of Slater-type (STO) basis sets (DZ, DZP, TZ2P, QZ4P) of double (D)-, triple (T), and quadruple (Q)-z quality with one (P), two (2P), or four (4P) sets of polarization functions. 28 To gauge dispersion effects on the structures, the stationary points were also optimized using Grimme's D3 dispersion correction and Becke-Johnson damping, 29 as denoted by BP86-D3(BJ)/QZ4P. The accuracies of the fit scheme (Zlm fit) 30a and the integration grid (Becke grid) 30b were set to VERYGOOD. As in the benchmark computations, the character of each DFT stationary point was verified by performing a harmonic vibrational analysis, 31 including careful inspection of the normal mode of imaginary frequency for each transition state. Finally, single-point energies of the stationary points were computed using the TZ2P and QZ4P basis sets in combination with a panoply of density functionals: the localdensity approximation (LDA) functional VWN; 32 35,40 mPW1K, 57 S12H, 58 Table S3 of the ESI. † The Diels-Alder reaction between 1,3-cis-butadiene (B) and ethylene (E) proceeds via a concerted synchronous transition state (TS-DA) to the six-membered cycloadduct (P-DA). The transition state is C s -symmetric with two nascent CÁ Á ÁC bond distances of 2.234 Å between the terminal carbon atoms of B and E. The CQC double bonds in B and E stretch by 0.167 and 0.194 Å to 1.509 and 1.531 Å, respectively, and become C-C single bonds in the cycloadduct P-DA; at the transition state, 26% and 28% of this evolution has occurred. Simultaneously, the central C-C bond of B shortens from 1.474 Å to 1.341 Å and becomes a CQC double bond in P-DA, this contraction being 50% complete in TS-DA. Twisting and pyramidalization of methylene groups are key motifs of the DA reaction. The twisting is best characterized by the dihedral angle [t twist = t(C 2 C 3 C 4 -H 9 C 4 H 10 ) = t(C 1 C 2 C 3 -H 5 C 1 H 6 )] between the terminal CH 2 plane and that of the connected C-C-C backbone in B. In particular, t twist starts at 1.81 in B and ends at 82.11 in P-DA, while displaying 33% of this change in TS-DA. The pyramidalization is quantified by the angles formed by the C-C bond vectors out of the CH 2 planes [g pyr (B) = g(C 2 C 1 H 6 H 5 ) = g(C 3 C 4 H 10 H 9 ) and g pyr (E) = g(C 12 C 11 H 14 H 13 ) = g(C 11 C 12 H 15 H 16 )]. As the terminal carbon atoms hybridizes from sp 2 to sp 3 , g pyr (B) rises from almost 01 in reactant B to 56.41 in P-DA, achieving 37% of this transformation in TS-DA. In remarkable accord, g pyr (E) likewise changes from 01 in E to 55.01 in P-DA, exhibiting 38% of this increase in TS-DA. Overall, the key geometric parameters reflective of bond rearrangement consistently indicate that TS-DA is an ''early'' transition state, nicely consistent with the Hammond Postulate and the large exothermicity (DH 0 = À40.6 kcal mol À1 ) of the reaction.
The 1,3-dipolar cycloaddition between methyl azide (MA) and ethylene (E) proceeds, in contrast to the Diels-Alder case, via a concerted asynchronous transition state (TS-1,3-DC) to the five-membered cycloadduct (P-1,3-DC). The nascent N 1 -C 8 and N 3 -C 9 bonds in the transition state (2.152 and 2.104 Å) differ in length by 0.048 Å, and this marked asymmetry persists in the N-C bonds of the cycloadduct (1.493 and 1.465 Å). The two NÁ Á ÁC distances in TS-1,3-DC were previously computed by Blavins and Karadakov 9c with a variety of methods (HF, MP2, B3LYP, CCD, QCISD, MP4) and small Pople basis sets [6-31G (d) or 6-31G(d,p)]. Their N 1 -C 8 and N 3 -C 9 results span the ranges 2.102-2.159 and 2.028-2.178 Å, respectively, with the QCISD distances coming within À0.006 and +0.011 Å of our CCSD(T)/ cc-pVTZ values. The N 1 -N 2 , N 2 -N 3 , and CQC multiple bonds of MA and E stretch by 0.108, 0.160, and 0.186 Å, respectively, in forming the cycloadduct; the percentage of this elongation found in the transition state is 30%, 23%, and 25%, in order. Regarding pyramidalization of the ethylene moiety, the out-ofplane angles [g(C 9 C 8 H 11 H 12 ), g(C 8 C 9 H 12 H 13 )] increase from 01 in the reactant to 46.301 and 47.731, respectively, in the product, while comprising 35% and 38%, respectively, of this transformation in TS-1,3-DC. Therefore, the key bond distances and pyramidalization angles show uniform trends that characterize TS-1,3-DC as a prototypical ''early'' transition state in accord with the Hammond Postulate and the reaction exothermicity (DH 0 = À23.2 kcal mol À1 ). The azide bending angle is incongruous in this respect, because it contracts from 173.51 in MA to 112.61 in the cycloadduct and exhibits more than half of this deformation in TS-1,3-DC. The spectator methyl substituent of the system actually shows nonmonotonic internal rotation, as the torsion angle t(N 2 N 3 C 4 H 5 ) rocks from that in MA by À61 in TS-1,3-DC but +141 in P-1,3-DC.
The electrocyclic rearrangement of 1,3-cis-butadiene is an intramolecular concerted cyclization in which a new C-C bond is formed between the terminal sp 2 -hybridized carbon atoms of B. The C s -symmetric transition state (TS-ER) has a nascent CÁ Á ÁC bond distance of 2.144 Å, which is reduced to 1.572 Å in the strained, four-membered cycloadduct (P-ER, C 2v ). These distances are 0.010 and 0.007 Å longer, respectively, than the CCSD/6-311++G(df,p) results of Minkin and coworkers. 9b During the reaction, the C-C single bond of B contracts from 1.474 Å to become a CQC double bond of length 1.347 Å in P-ER, and 76% of this decrease has occurred in the transition state. Simultaneously, the CQC double bonds of B elongate by 0.180 Å to become C-C single bonds in P-ER, exhibiting 51% of this change in TS-ER. The C-C-C angle changes from 124.51(B) -104.31(TS-ER) -94.21(P-ER) during the reaction, exhibiting 67% completion in the transition state. As in the DA reaction, the ER process requires twisting and pyramidalization of the terminal methylene groups in B. The aforementioned t twist angle starts at 1.81 in B, ends at 901 in P-DA, and has achieved 57% of this change in TS-DA. Because the ER reaction is endothermic (DH 0 = +9.6 kcal mol À1 ), a ''late'' transition state is anticipated, consistent with the trends seen in the characteristic geometric parameters.
The SR transformation of (Z)-penta-1,3-diene (P) is an intramolecular identity reaction comprising a [1,5]-sigmatropic hydrogen shift between the terminal carbon atoms. The transition state (TS-SR) has C s symmetry, with a 1.417 Å distance for both the breaking and forming CÁ Á ÁH bonds. The reorganization of the p-system causes the backbone bond distances, i.e., C 1 -C 2 , C 2 -C 3 , C 3 -C 4 , C 4 -C 5 , to go from 1.343, 1.474, 1.347, 1.503 Å in the reactant (P), to 1.420, 1.398, 1.398, 1.420 Å in TS-SR, and finally to the re-ordered values 1.503, 1.347, 1.474, 1.343 Å in the identity product (P). The corresponding evolutions of these backbone distances in TS-SR (48%, 59%, 41%, 52%) are symmetrically displaced from the 50% progression simplistically expected for an identity reaction. The backbone bond angles of P, [y(C 1 -C 2 -C 3 ), y(C 2 -C 3 -C 4 ), y(C 3 -C 4 -C 5 )] = 125.91, 127.21, 127.01, contract by 3.21, 7.31, and 4.31 to form TS-SR before relaxing back to the initial values in reverse order. As terminal carbon C 1 in the reactant transitions from sp 2 to sp 3 hybridization, the out-of-plane angle g(C 2 C 1 H 6 H 7 ) undergoes the progression À1.31 (P) -33.91 (TS-SR) -53.61 (P), whereas the methylene twist angle, t twist = t(C 3 C 2 C 1 -H 6 C 1 H 7 ), displays the simultaneous evolution 1.41(P) -52.01(TS-SR) -78.91(P). In both the pyramidalization and twisting cases, about 65% of the total reorganization occurs before the transition state is reached. Notably, the expectation of the Hammond Postulate of a central transition state for a thermoneutral reaction is not generally upheld for the various SR geometric transformations.
The DGT identity reaction between ethane (Et) and ethylene (E) proceeds in a concerted and synchronous fashion with two hydrogen atoms migrating simultaneously from Et to E. The corresponding transition state (TS-DGT) of D 2h symmetry features a six-membered ring in which the breaking C-H bonds are elongated by 0.272 Å to 1.363 Å and the newly forming bonds have the same distance. Comparing with the work from Radom and coworkers 9a reveals that CCSD(T)/6-31G(d,p) theory gives the CÁ Á ÁH distance and CÁ Á ÁHÁ Á ÁC angle of TS-DGT within 0.003 Å and 0.51, respectively, of our benchmark values. The DGT switches the C-C and CQC bonds in the fragments, and the dual 1.424 Å distances in the transition state correspond to 55% of the contraction and 45% of the elongation that occurs in the overall reaction. TS-DGT exhibits four equivalent pyramidalizations, which all have a 29.31 angle of the C-C bonds out of adjacent methylene plane; these angles represent 56% of the total shift occurring in the E(sp 2 ) -Et(sp 3 ) rehybridization.
The TS-DGT carbon-carbon distances and pyramidalization angles agree with simple expectations for an identity reaction. However, this accord is not met by the C-H bonds peripheral to the hydrogen exchange; in particular, the corresponding distances and methylene angles in TS-DGT respectively constitute only 29% of the increase and 34% of the decrease involved in the E(sp 2 ) -Et(sp 3 ) transformation.

Benchmark pericyclic reaction energetics
The focal-point analyses to determine the benchmark reaction energies and barriers for our five paradigmatic pericyclic reactions are presented in Tables 1 and 2. The two-dimensional grids detail the systematic convergence of the energetic predictions toward both the complete basis set (vertical) and electron correlation (horizontal) limits, as approached via the aug 0 -cc-pVXZ atomic-orbital series (X = D, T, Q, 5) and the HF -MP2 -CCSD -CCSD(T) -CCSDT(Q) wave function hierarchy, respectively. Attesting to the excellent basis set convergence in all eight FPA tables, the explicitly computed aug 0 -cc-pV5Z results for RHF relative energies and MP2 correlation increments display mean absolute deviations of only 0.02 and 0.09 kcal mol À1 , respectively, with respect to the extrapolated CBS values. The entry-level RHF/CBS method greatly overestimates the barriers by 13.5-32.3 kcal mol À1 , while the corresponding reaction energy errors range from +2.4 to +9.0 kcal mol À1 . The MP2 corrections to the barriers are all negative and overly large in magnitude, to the extent that MP2/ CBS theory underestimates the barriers by 2.8-9.4 kcal mol À1 ; concomitantly, the reaction energy errors vary from À5.2 to +3.5 kcal mol À1 . Oscillations of computed barrier heights continue in the next step of the electron correlation hierarchy, as the CCSD/CBS correlation increments are all positive and this method overestimates the barriers by 3.2-6.4 kcal mol À1 ; however, the corresponding mean absolute error in the reaction energies is now 1.2 kcal mol À1 . Predictably, the (T) increments for all barriers are negative, reducing the overestimation of the barriers to only 0.08-0.26 kcal mol À1 at the CCSD(T)/CBS level of theory; all reaction energies given by this method are accurate to better than 0.8 kcal mol À1 . The CCSD(T) -CCSDT and CCSDT-CCSDT(Q) contributions are opposite in sign and largely cancel one another for the barriers, giving very small net d[CCSDT(Q)] corrections between À0.07 and À0. 16  The 2019 paper by Karton 8 reports W3lite-F12 computations of three reaction barriers (without ZPVE corrections) related to our pericyclic reaction systems: TS-DA arising from the Diels-Alder reaction between E and the trans-isomer of B, ring opening reaction (ROR) of P-ER, that is, the reverse reaction of the electrocyclic rearrangement (TS-ER with respect to P-ER), and TS-SR originating from the sigmatropic rearrangement of the E-isomer of P. Note, the trans-isomer of B and the E-isomer of P are not reached by the intrinsic reaction paths connected to TS-DA and TS-SR, but these isomers are lower in energy than their cis and Z counterparts shown in Fig. 1. We expanded our computations to include B(trans) and P(E) and re-worked the FPA tables to obtain the following DE final barrier predictions: .85 (TS-SR) kcal mol À1 , so that the largest difference between the W3lite-F12 and FPA barriers is a mere 0.09 kcal mol À1 . Notwithstanding this close agreement and the fact that both approaches include correlation effects up to CCSDT(Q)/cc-pVDZ, our FPA computations are more favorable in several respects for the three cases of concern. (1) Geometric structures were optimized at the CCSD(T)/ cc-pVTZ level of theory rather than at B3LYP-D3(BJ)/cc-pVTZ.
(2) Aug 0 -cc-pV5Z HF and MP2 energies were used to assist the FPA extrapolations to the CBS limit, as opposed to omitting MP2 in the hierarchy and employing the cc-pVQZ-F12 approach for the largest explicit computations. (3) More complete basis sets with diffuse heavy-atom functions [aug 0 -cc-pV{T,Q}Z] were utilized to extrapolate CCSD increments via an established X À3 form to the CBS limit, in contrast to evaluating such increments by means of cc-pV{T,Q}-F12 extrapolations with a phenomenological X À5.94 function. (4) CCSD(T) increments were also obtained from aug 0 -cc-pV{T,Q}Z extrapolations to the CBS limit rather than by scaling cc-pVTZ-F12 results by an ad hoc factor of 0.987. (5) CCSD(T)/cc-pVTZ harmonic vibrational frequencies and ZPVEs of all stationary points were computed as opposed to a complete omission of all vibrational effects. (6) The DBOC terms were evaluated and used to confirm the absence of surface crossing phenomena in the pericyclic reactions. (7) Reaction energies were computed rather than forward barriers alone, thus providing results of the same rigor for the reverse reactions.

Assessment of DFT geometric structures
In this section, geometric structures of the reactants, transition states, and products of each pericyclic reaction were fully optimized with five DFT methods: BP86/DZ, BP86/DZP, BP86/ TZ2P, BP86/QZ4P, and BP86-D3(BJ)/QZ4P. Cartesian coordinates of these DFT structures are provided in Tables S4-S8  (ESI †). To assess the quality of these DFT geometries, they were compared to the CCSD(T)/cc-pVTZ benchmark geometries by means of a Cartesian root-mean-square deviation (RMSD) analysis. The procedure rigorously solved for the 3 origintranslation and 3 axis-rotation variables that provide the minimum possible RMSD between the Cartesian coordinates of a given DFT geometry and those of the corresponding CCSD(T)/ cc-pVTZ benchmark geometry. The solution to this minimization  problem in the translational space is always achieved by making the centroids exactly coincident for the two structures under comparison. However, the RMSD minimization in the rotational space is not achieved in the general case by assuming unit masses and aligning the principal axes of inertia, although this approach often yields an excellent approximation. Analytic solutions of the RMSD rotational minimization problem have been reported by Kabsch, 68 and we have derived improved analytic formulas from the Euler-Rodrigues representation of the associated unitary matrix. Finally, we wrote a Mathematica program to execute numerical searches for the desired solutions (see ESI, † for Mathematica code CartRMSD). All three approaches gave identical results for every set of molecular structures, making our RMSD minimizations definitive.
The global-minimum Cartesian RMSD values for the five DFT methods are plotted in Fig. 3 for all stationary points of our pericyclic reactions. Note that for all methods the largest deviations occur for the transition state of the 1,3-dipolar cycloaddition reaction (TS-1,3-DC). Interestingly, the BP86/DZP geometries generally display the smallest deviations from the CCSD(T)/cc-pVTZ benchmarks with an average RSMD of 0.027 Å (Table S9, ESI †). Expanding the BP86 basis set to TZ2P and QZ4P generally increases the RMSD and thus does not improve accuracy despite the added computational expense. Moreover, the inclusion of Grimme's D3 dispersion in conjunction with Becke-Johnson damping does not always improve the BP86/QZ4P geometries, and overall BP86-D3(BJ)/QZ4P is still substantially inferior to BP86/DZP. In order to pinpoint where the DFT geometries exhibit disparities, we considered all geometric parameters identified in Fig. 1 and 2; and calculated the average percentage error of the distances, angles, and dihedral angles with respect to the CCSD(T)/cc-pVTZ benchmarks (Table S10, ESI †). Consistent with Fig. 3, the percentage errors of the distances, angles, and dihedral angles have both the smallest average and the smallest maximum for the BP86/DZP method. Thus, a pragmatic approach to obtain geometries closest to CCSD(T)/cc-pVTZ structures is to optimize stationary points at BP86/DZP. Independent of the DFT level of theory, the largest percentage errors in the internal coordinates are found for dihedral angles, and these errors are responsible for most of the deviations seen in the Cartesian RMSD analysis. In particular, dramatic deficiencies occur for the dihedral angles of the B, TS-1,3-DC, P-1,3-DC, and P stationary points. In contrast, the bond distances and angles for all DFT levels of theory differ on average by less than 1% from the CCSD(T)/cc-pVTZ reference values.

Performance of DFT methods for energetics
Next, we examine the relative energies of the stationary points computed with the density functionals (XC) specified in the Theoretical methods section in conjunction with a QZ4P Slatertype basis set atop the CCSD(T)/cc-pVTZ geometries, denoted as XC/QZ4P//CCSD(T)/cc-pVTZ. The DFT reaction barriers and reaction energies for the five pericyclic reactions are collected in Table 3. The performance of the various density functional approximations is assessed by a systematic comparison of the resulting reaction barriers and energies with our most accurate CCSDT(Q)/CBS benchmark values (Table 1). For all 60 functionals, we have computed the mean error (ME), mean absolute error (MAE), maximum unsigned error (MUE), and standard deviation (SD) of the reaction barriers, reaction energies, as well as the combination of both, relative to the FPA CCSDT(Q)/CBS targeted energy values for all pericyclic reactions ( Table 4). Note that the errors in reaction barriers and reaction energies for all stationary points computed at XC/QZ4P//CCSD(T)/cc-pVTZ with respect to CCSDT(Q)/CBS//CCSD(T)/cc-pVTZ can be found in Table S11 (ESI †). Table 4 and Table S11 (ESI †) show that the LDA functional VWN yields too low reaction barriers with a MAE of 17.3 kcal mol À1 and in the worst case of the double group transfer an MUE of 28.3 kcal mol À1 , which is likely due to the well-known overbinding of this functional. Next to that, VWN also leads to extremely exothermic reactions with a MAE and MUE of 12.4 kcal mol À1 and 17.8 kcal mol À1 , respectively. Additionally, the GGA functional PBEsol and hybrid functionals O3LYP and BHandH perform nearly as poorly as LDA, having a total MAE of 10.6, 18.9, and 8.3 kcal mol À1 , respectively.
The best overall agreement with our ab initio benchmark reaction barriers and reaction energies is obtained by representatives from the meta-hybrid and double-hybrid family, with total MAEs of 1.1 (M06-2X), 1.4 (B2K-PLYP), and 1.5 kcal mol À1 (mPW2K-PLYP and revDSD-PBEP86). The best performing density functional approximation, M06-2X, also has the lowest MUE and SD of only 2.4 and 0.6 kcal mol À1 , respectively. The three best performing double-hybrid functionals (B2K-PLYP, mPW2K-PLYP, and revDSD-PBEP86), on the other hand, have a slightly higher MUE and SD, with a MUE around 4 kcal mol À1 and a SD around 1 kcal mol À1 . Nevertheless, the error distribution of these functionals is not uniform as they under-and overestimate some of the reaction barriers and energies depending on the reaction, i.e., the ME is, due to error cancelation, unequal to the MAE. For instance, mPW2K-PLYP underestimates, amongst others, the reaction barrier of the Diels-Alder reaction by 3.9 kcal mol À1 , but, in contrast, overestimates, e.g., the reaction barrier of the electrocyclic rearrangement by 0.7 kcal mol À1 . We want to highlight that, despite a somewhat higher MAE of 5.8 kcal mol À1 , the still popular density functional approximation BP86 4a,b,7d is able to describe the trends in reaction barriers and reaction energies qualitatively accurately.
Karton and Goerigk proposed PWPB95-D3(BJ) as best performing density functional approximation, with a root mean square deviation in reaction barrier heights of 1.0 kcal mol À1 with respect to their WnF12 reference data for four classes of pericyclic reactions (26 reactions in total). 5c The difference between our results and those of Karton and Goerigk could be ascribed to several factors. (1) The density functional approximations in this work are compared to more accurate reference data computed on more accurate geometries (vide supra).
(2) The density functional approximations are, in contrast to the work of Karton and Goerigk, not only assessed based on their performance of accurately describing the reaction barrier, but also the reaction energies and thus the reverse reaction.
(3) Karton and Goerigk exclusively assess the performance of dispersion-corrected density functional approximations, which, as we show later, artificially lower the reaction barrier for pericyclic reaction that involve small reactants with minimal dispersive intermolecular interactions. When solely concentrating on the reaction barriers, as done by Karton and Goerigk, we find three other density functional approximations that outperform, even with respect to our more accurate reference  Table 3 Reaction barriers (D E ‡ ) and energies (DE rxn ) (kcal mol À1 ) for the pericyclic systems of Scheme 1, computed using various density functional approximations at the XC/QZ4P//CCSD(T)/cc-pVTZ level XC Diels-Alder reaction 1, 3  data, the PWPB95-D3(BJ) functional proposed by Karton and Goerigk (vide infra). This is, however, not the reason why PWPB95-D3(BJ) is not included in this DFT performance study. The reason why PWPB95-D3(BJ) is not included in this DFT performance study is that this XC functional is not available in the ADF software package. Evaluating the performance of the functionals in describing the reaction barrier and reaction energy separately reveals that the top three best functionals for accurately describing the reaction barriers are the double-hybrid functionals mPW2-PLYP, B2T-PLYP, and B2-PLYP with a MAE, for all reactions together, of only 0.7, 0.8, and 0.9 kcal mol À1 , respectively, and a MUE of no more than 2.0 kcal mol À1 . The best functionals for evaluating the reaction energies are, on the other hand, the meta-hybrid functional M06-2X, with a MAE of only 0.8 kcal mol À1 and a MUE of 1.2 kcal mol À1 , followed by the range-separated hybrid functionals oB97X-D that has a MAE of 1.1 kcal mol À1 and a MUE of around 2.2 kcal mol À1 and CAM-B3LYP that has a MAE of 1.4 kcal mol À1 and a MUE of around 2.3 kcal mol À1 . Although the double-hybrid functionals are able to accurately describe the reaction barriers and energies, a disadvantage of these functionals, up to this point, is that, in the context of an STO-based approach, they are computationally less efficient for geometry optimization or transition state search calculations. However, an effective workaround is to perform single-point energy calculations on stationary point geometries calculated using a more computationally more efficient functional (vide infra).
Next, we examine how the various DFT approximations perform for each individual one of the five pericyclic reactions (Table S11, ESI †). Interestingly, we find that, despite M06-2X, B2K-PLYP, mPW2K-PLYP, and revDSD-PBEP86 are able to accurately describe pericyclic reactions in general, they are not per definition the best functionals for describing the individual potential energy surfaces of the pericyclic reactions studied in this work. The most accurate energies for the Diels-Alder reaction are obtained using the M06, B1PW91, and OPBE functionals, with an error of less than 1.0 kcal mol À1 with respect to the CCSDT(Q)/CBS reference energies. The 1,3-dipolar cycloaddition, on the other hand, is described the best using the functionals S12H, PBE0, and revDSD-PBE, which have an error between 1.0-2.0 kcal mol À1 . The DFT approximations B2K-PLYP, mPW2K-PLYP, and revDSD-PBEP86 are, with their error of less than 1.1 kcal mol À1 , the best for describing the electrocyclic rearrangement. The reaction barrier of the sigmatropic rearrangement can be approximated within an error of only a few tenths of a kcal mol À1 with B2T-PLYP, mPW2-PLYP, and M06-2X. At last, mPW1K, oB97X-D, and KMLYP are able to approach the reference energies of the double group transfer with 0.0-0.1 kcal mol À1 .
Adding Grimme's empirical D3 dispersion correction, in combination with the Becke-Johnson damping, denoted as D3(BJ), results in a significant worsening of the performance of the GGA, hybrid, and meta-hybrid functionals with a maximum increase in total MAE of 6 kcal mol À1 . The D3(BJ) dispersion correction leads to an additional stabilization of the stationary points, which increases the deviation from the CCSDT(Q)/CBS reference energies when a functional is already underestimating the stationary point energy. For instance, BP86 underestimates the reaction barrier of the pericyclic reaction with errors ranging from À2.7 kcal mol À1 for the 1,3-dipolar cycloaddition reaction to À12.1 kcal mol À1 for the double group transfer and hence has a MAE, for the reaction barriers, of 6.5 kcal mol À1 . Adding the D3(BJ) dispersion correction, i.e., BP86-D3(BJ), aggravates the underestimation of the reaction barriers by this functional, increasing the errors to À7.1 kcal mol À1 for the 1,3-dipolar cycloaddition reaction and À15.5 kcal mol À1 for the double group transfer and, therefore, increases the MAE, for the reaction barriers, to 9.4 kcal mol À1 . This shows that including dispersion corrections to the evaluation of the reaction barriers of pericyclic reactions involving small reactants that exhibit minimal dispersive intermolecular interactions, such as steric attraction or hydrogen bonds, leads to a worsening (not improvement) of the accuracy of the computed energies. The D3(BJ) dispersion correction does, however, have a positive effect on the reaction energies. As these functionals generally underestimate the stability of the adducts, they can be partly compensated by including an additional stabilizing dispersion correction and hence reduces the MAE of the reaction energies. Notably, adding Grimme's D3 dispersion correction to the overall best performing functional M06-2X, i.e., M06-2X-D3, leads to a slight worsening of the accuracy of this functional by increasing the MAE to 1.2 kcal mol À1 .
Finally, we have additionally evaluated the performance of the 52 functionals and eight dispersion-corrected functionals using (i) a TZ2P basis set on CCSD(T)/cc-pVTZ geometries and (ii) a QZ4P and TZ2P basis set on DFT geometries optimized at BP86/DZ, BP86/DZP, BP86/TZ2P, BP86/QZ4P, BP86-D3(BJ)/ QZ4P (vide supra). We find that, for the CCSD(T)/cc-pVTZ geometries, there is almost no difference between energies computed with the TZ2P and QZ4P basis set, since the variation in MAE from TZ2P is QZ4P is only 0.2 kcal mol À1 . Hence, the energies are converged with a TZ2P basis set (Tables S13-S15, ESI †). Notably, the energies evaluated using the density functional approximations in this work do not always have the best agreement with the CCSDT(Q)/CBS benchmark values when employing the best, i.e., largest, basis set. This indicates an increased error cancelation for these functionals with a smaller basis set. Evaluating the performance of the density functional approximation using DFT geometries optimized at BP86/BS, where BS = DZ, DZP, TZ2P, QZ4P, and BP86-D3(BJ)/QZ4P give results with similar accuracy as the evaluation based on the ab initio reference geometries at CCSD(T)/cc-pVTZ (Tables S11-S45, ESI †). The differences in total MAE between the computations performed on the DFT and ab initio geometries are within a few tenths of a kcal mol À1 . Even the earlier established poor BP86/DZ geometries give a total MAE that differs only a few tenths of a kcal mol À1 with respect to the very accurate CCSD(T)/cc-pVTZ geometries. In analogy with the evaluation of the functionals on the ab initio geometries, we find that for the DFT optimized geometries, the best overall performance is again obtained by members of the meta-hybrid and doublehybrid family, namely, M06-2X, B2K-PLYP, mPW2K-PLYP, and revDSD-PBEP86, having a total MAE of no more than 1.5 kcal mol À1 (maximum total MAE for BP86/DZ geometries is 2.3 kcal mol À1 ) and an MUE of a few kcal mol À1 with respect to the CCSDT(Q)/CBS reference energy values. These results unequivocally show that one can obtain very accurate pericyclic reaction barriers and reaction energies by optimizing the geometries of the stationary points using an affordable DFT approach, such as BP86/DZP, followed by a refinement of the energy by performing single-point calculations using the more accurate, but also more expensive, hybrid functional M06-2X or a double-hybrid functionals B2K-PLYP, mPW2K-PLYP, and revDSD-PBEP86.

Conclusions
We have performed a comprehensive ab initio benchmark study yielding reference data of unprecedented accuracy for various pericyclic reactions, including the Diels-Alder reaction, 1,3dipolar cycloaddition, electrocyclic rearrangement, sigmatropic rearrangement, and double group transfer. Our reference data for reaction barriers and reaction energies emerge from a hierarchical series of ab initio quantum mechanical methods, up to CCSDT(Q), in combination with a series of correlationconsistent Gaussian-type basis sets, up to the complete basis set (CBS) limit, along which they converge within a few tenths of a kcal mol À1 .
We have analyzed the performance of 60 density functional approximations for describing the above-mentioned pericyclic reactions against our best ab initio benchmark data at CCSDT(Q)/ CBS and find that the meta-hybrid functional M06-2X is the best performing density functional approximation, with a mean absolute error (MAE) of only 1.1 kcal mol À1 and a maximum unsigned error (MUE) of 2.4 kcal mol À1 . The mean error (ME) and standard deviation (SD) of M06-2X are also small, only À0.3 and 0.6 kcal mol À1 , respectively. This meta-hybrid functional is closely followed by representatives from the double-hybrid family, namely, B2K-PLYP (MAE = 1.4 kcal mol À1 ), mPW2K-PLYP, and revDSD-PBEP86 (MAE = 1.5 kcal mol À1 ). The BP86 functional, which has frequently been employed to study pericyclic reactions, has a higher MAE of 5.8 kcal mol À1 . Inspection, however, shows that BP86 trends in barriers and reaction energies are qualitatively accurate. Furthermore, we found that, for practical application, the energies computed using most of the 60 DFT functionals reach their convergence with the TZ2P basis set. In most cases, the variation in MAE from TZ2P to QZ4P is 0.2 kcal mol À1 or less (e.g., 0.1 kcal mol À1 for M06-2X and BP86).
Finally, we find that an efficient and still accurate protocol for exploring pericyclic reactions consists of computing BP86/DZP geometries followed by single-point corrections, using the TZ2P basis set, with the more sophisticated density functional approximations M06-2X, B2K-PLYP, mPW2K-PLYP, or revDSD-PBEP86.

Conflicts of interest
There are no conflicts to declare.