On the method-dependence of transition state asynchronicity in Diels–Alder reactions

Abstract

This work discusses the dependence of transition state geometries on the choice of quantum chemical optimization method for the extensively studied Diels–Alder reaction. Rather significant differences are observed between post-Hartree–Fock methods and (hybrid) density functional theory, where the latter predicts larger asynchronicities. The results show that the low MP2 asynchronicity observed is likely artificial. Still, there are significant discrepancies between hybrid and pure density functionals. The role of the exchange functional seems to be most prominent in less activated reacting systems, while the importance of the correlation functional seems to increase as they become more activated by, e.g., an electron-donating group on the diene. To correct the dubious MP2 geometries, we employed the SCS-MP2 protocol for transition state optimization, which leads to significantly better results with respect to CCSD/6-31+G(d) level calculations. We conclude that in order for hybrid functionals to give descriptions consistent with the sample post-Hartree–Fock methods, a balanced combination of both Hartree–Fock exchange (with a couple of exceptions) and a well-behaved correlation functional is required. Given that the benchmark CCSD/6-31+G(d) geometries are sufficient representations, the best geometries were obtained using ωB97X(D), B2PLYP(D) and M06-2X.

---

1 Introduction

Becke's 3-parameter exchange–correlation functional (B3LYP)^1–3 has been the most popular choice for Density Functional Theory (DFT) calculations in theoretical organic chemistry during the last decade.⁴ Now, close to its 20th anniversary, problems for the functional to treat several systems have been highlighted.⁵ A pertinent example is cycloadditions.⁶ It was recently demonstrated that the large errors in reaction energy originate from an erroneous treatment of π → σ transitions, which increase with basis set size.⁷

The novel meta-hybrid-GGA functional M06-2X⁸ has recently been employed in a number of theoretical studies of Diels–Alder and related mechanisms.^9–15 In most cases, B3LYP has still been used for optimizations,^7,9,11–13 owing to its good reputation for reproducing geometries, while M06-2X has typically been employed for energy calculations.¹⁶ Some recent studies have however begun to use M06-2X for geometries as well.^10,14,15,17 We note that in contrast to the ongoing challenge from M06-2X and other functionals, a recent study concluded that B3LYP performed very well for predicting geometries.¹⁸

Despite the ‘reliability’ of the M06-2X//B3LYP approach, we have noticed^14,17 that B3LYP sometimes produces significantly different Diels–Alder transition state (TS) geometries compared to M06-2X and MP2.²⁰ The tendency is for B3LYP to yield very asynchronous TSs for activated reactants (see Fig. 1a for a definition of the asynchronicity Δd). A rather dramatic example is given in Fig. 2:¹⁷ the B3LYP description of this ionic TS is hardly reminiscent of the conventional picture of a Diels–Alder reaction. Indeed, a B3LYP-optimized Diels–Alder reaction between ionic reactants, or catalyzed by a Lewis acid, is generally seen as a stepwise process^21,22 (although older wave-function calculations give a different description^23,24). However, in the example from ref. 18, an IRC²⁵ calculation at the MP2/6-31+G(d) level indicated a concerted process, consistent with the less asynchronous TS in Fig. 2. Hence the accepted picture of polar Diels–Alder reactions being increasingly asynchronous with increasing activation and charge transfer^19,21,26 may be biased by the choice of computational method.

Fig. 1 (a) Atom enumeration and definition of asynchronicity used in this work. All considered transition states have endo-cis geometry, as this is normally found to be the one with lowest energy.¹⁹ (b) Reactants used in this study. The methyl group in 3 and the methoxy group in 4 activate the diene's 4π system by being weak electron donors, while the dimethylamino group in 5 is a relatively strong donor. A hydrogen bond donor (such as water) conversely activates the dienophile by withdrawing electrons from the conjugated system, acting as a weak Lewis acid.

Fig. 2 Transition state geometries of the Diels–Alder reaction between acrolein and the acid-coordinated enolate of cyclopentenone, optimized using three different methods (see ref. 18).

While such details may not always be the main interest in a computational study aiming at reproducing experimental energies, they become central in mechanistic studies, for example in the elucidation and design of molecular (or biomolecular²⁷) catalysts. Another aspect is the rising awareness that stepwise Diels–Alder is a rather common phenomenon,^13,14,28 and distinction between a concerted and two-step mechanism of course requires a reliable method.

With regard to the overwhelming dominance of B3LYP in computational organic chemistry and the emergence of new functionals such as M06-2X, we find it motivating to compare the geometric details of Diels–Alder TSs obtained using a range of popular and/or novel density functionals. For this purpose, we have used a number of common reactants (Fig. 1b) with electron-donating substituents and/or a water molecule gradually activating the Diels–Alder reaction. The results are discussed in terms of variations in asynchronicity with activation, and are compared to a few wave-function based methods. Finally, we comment on the accuracy in predicting activation energies in comparison to geometry.

2 Computational methods

Our survey mainly concerns methods available in the 2009 release of the Gaussian package.²⁹ Apart from B3LYP, we have considered a number of related methods based on Becke's GGA exchange functionals: BLYP,^2,3 B3P86,³⁰ B3PW91³¹ B97D,³²ωB97X, ωB97XD,³³ B2PLYP and B2PLYPD,³⁴ (functionals with a ‘D’ index contain an empirical dispersion correction³⁵). Within the meta-GGA class we evaluated M06-2X, M06L³⁶ and HCTH/407.³⁷ In addition, we employed two analytic functionals of Perdew and coworkers, mPW1PW9^31,38 and PBE1PBE (also known as PBE0).³⁹ Note that BLYP, B97D, M06L, and HCTH/407 are ‘pure’ exchange–correlation functionals, viz., they do not contain any exact exchange. The choice of functionals allows us to some extent to trace the origin of certain geometric features. For example, by comparing results from the B3LYP, mPW1PW91 and B3PW91 functionals, one can determine whether a potential anomaly stems from the exchange or the correlation part. Comparing pure and hybrid functionals may furthermore illustrate the importance of including exact exchange.

The functionals were evaluated based on the optimized TS deviations from geometries calculated using coupled cluster with single and double substitutions (CCSD).⁴⁰ TSs for all systems except 1·W + 3, 1·W + 4 and 1·W + 5 were optimized using CCSD. For reference, geometries were also optimized at the MP2 and HF levels. The calculations were performed using Gaussian09. The 6-31+G(d,p) basis set was used with all methods except CCSD, for which we employed the smaller 6-31G+(d) basis set. Frequency calculations were performed for each optimized TS to verify that it contained exactly one imaginary frequency corresponding to the reaction coordinate. Although a small basis set may increase the basis set superposition error's (BSSE) effect on the geometries for high-level methods, it should be noted that diffuse functions, such as present in 6-31+G(d), are known to significantly reduce the BSSE.⁴¹

Since SCS-MP2⁴² has been shown to perform well for Diels–Alder reaction and barrier energies,^6,7 we employed the Turbomole suite⁴³ to perform complementary calculations of all systems (vide infra). These optimizations were performed using the RICC2 module and the 6-31+G(d,p) basis set to obtain both MP2 and SCS-MP2 geometries. As an auxiliary basis set we employed the Alrichs type QZVPP.⁴⁴ In contrast to the Gaussian calculations, the Turbomole computations did not invoke the frozen core approximation and used spherical d-functions. Despite the slightly different approach, however, the RI-MP2 geometries are virtually identical to those produced using the Gaussian protocol, as seen below (Tables 1 and 2). The Turbomole-derived SCS-MP2 geometries can thus be directly compared to the results obtained with Gaussian.

Table 1 Incipient bond distances for the uncatalyzed systems^a

Method	1+2		1+3		1+4		1+5
	d β	d α	d β	d α	d β	d α	d β	d α
a All distances are given in Ångströms. b Optimized using the Turbomole suite and the 6-31+G(d,p) basis set (vide infra).
CCSD	2.026	2.522	2.003	2.597	1.945	2.729	2.032	3.000
HF	2.032	2.413	1.992	2.495	1.881	3.209	1.979	3.204
BLYP	2.013	2.800	2.021	2.859	1.990	3.004	2.117	3.159
B3LYP	2.028	2.670	2.014	2.764	1.979	2.885	2.093	3.083
B3P86	2.084	2.686	2.074	2.754	2.054	2.832	2.183	3.020
B3PW91	2.080	2.682	2.069	2.737	2.042	2.845	2.162	3.039
B97D	2.034	2.708	2.040	2.790	2.020	2.890	2.168	3.003
ωB97X	2.029	2.574	2.005	2.673	1.967	2.884	2.081	3.025
ωB97XD	2.044	2.616	2.021	2.724	1.990	2.850	2.108	3.008
PBE1PBE	2.093	2.673	2.080	2.730	2.056	2.828	2.181	3.017
mPW1PW91	2.086	2.661	2.070	2.718	2.046	2.818	2.165	3.015
HCTH/407	1.996	2.924	1.994	3.011	1.982	3.205	2.094	3.282
M06-2X	2.072	2.523	2.031	2.632	1.967	2.814	2.098	2.992
M06-L	2.036	2.821	2.047	2.904	2.054	2.999	2.162	3.020
B2PLYP	2.042	2.654	2.036	2.723	2.000	2.829	2.113	3.013
B2PLYPD	2.045	2.639	2.039	2.720	2.010	2.830	1.967	2.814
MP2	2.121	2.617	2.124	2.632	2.076	2.655	2.160	2.832
MP2^b	2.122	2.614	2.125	2.636	2.077	2.653	2.163	2.838
SCS-MP2^b	2.048	2.571	2.048	2.601	2.002	2.633	2.055	2.869

---

Table 2 Incipient bond distances for the water-catalyzed systems^a

Method	1·W+2		1·W+3		1·W+4		1·W+5
	d β	d α	d β	d α	d β	d α	d β	d α
a All distances are given in Ångströms. b Optimized using the Turbomole suite and the 6-31+G(d,p) basis set (vide infra).
CCSD	1.997	2.602
HF	1.983	2.501	1.963	2.597	1.907	3.264	2.029	3.258
BLYP	2.004	2.865	2.045	2.920	2.027	3.124	2.215	3.227
B3LYP	2.006	2.758	2.036	2.815	2.007	3.024	2.186	3.168
B3P86	2.076	2.762	2.113	2.796	2.088	2.959	2.312	3.107
B3PW91	2.068	2.761	2.098	2.794	2.070	2.967	2.261	3.123
B97D	2.028	2.771	2.080	2.810	2.060	2.980	2.316	3.055
ωB97X	2.004	2.666	2.028	2.716	2.001	3.019	2.159	3.091
ωB97XD	2.021	2.716	2.052	2.765	2.026	2.991	2.194	3.076
PBE1PBE	2.082	2.753	2.115	2.775	2.087	2.956	2.300	3.106
mPW1PW91	2.072	2.743	2.102	2.768	2.072	2.948	2.270	3.104
HCTH/407	1.994	2.989	2.024	3.029	2.018	3.290	2.149	3.349
M06-2X	2.029	2.618	2.026	2.698	1.994	2.964	2.192	3.049
M06-L	2.039	2.846	2.069	2.895	2.082	3.034	2.211	3.049
B2PLYP	2.028	2.727	2.066	2.751	2.024	2.958	2.202	3.078
B2PLYPD	2.034	2.710	2.081	2.740	2.040	2.930	2.190	3.015
MP2	2.118	2.663	2.133	2.681	2.079	2.747	2.185	2.893
MP2^b	2.119	2.662	2.136	2.682	2.078	2.689	2.189	2.894
SCS-MP2^b	2.035	2.619	2.057	2.640	1.997	2.682	2.092	2.950

---

3 Results and discussion

Incipient bond distances (d_β and d_α) are reported for the uncatalyzed systems in Table 1. Corresponding values for the water-mediated systems are given in Table 2. The resulting asynchronicities are presented in Fig. 3 as colored bars.

Fig. 3 Transition state asynchronicity per system, represented by colored bars. No CCSD geometries were optimized for the 1·W+3, 1·W+4 and 1·W+5 systems. Note that the SCS-MP2 geometries have been obtained using Turbomole.

As seen from the tables and figures, all methods predict the same overall trend of increasing asynchronicity with activation. Steric factors in 4 and 5 make it difficult to determine the purely electronic effects, and the 1·W + 5 TS is in fact less asynchronous than 1·W + 4. The general trend should however be obvious. MP2 asynchronicities are rather small compared to all DFT methods and, notably, increase only marginally with activation. In contrast, the CCSD description yields dramatic changes in going from 2 to 5, which makes it difficult to determine which DFT method gives the ‘best’ corresponding geometries. We note by comparing with geometries at the 6-31G(d) level (not shown) that the diffuse functions in 6-31+G(d) have a significant effect on the CCSD asynchronicity, which indicates flat potential energy surfaces with respect to the incipient bonds. CCSD/6-31G(d) asynchronicities were found to be similar to those obtained with MP2/6-31G(d), in stark contrast to the CCSD/6-31+G(d) ones shown in Fig. 3. MP2 appears much less sensitive to the diffuse basis set; the 6-31+G(d) geometries are only slightly more asynchronous than the 6-31G(d) ones. We can therefore not conclude if the CCSD geometries are entirely converged. We will henceforth assume that CCSD/6-31+G(d) yields the best available geometries and use it as a benchmark in the following. However, MP2 calculations with larger basis sets (vide infra) indicate a slight decrease in asynchronicity due to the additional polarization functions.

The most spurious results are apparently found with HF and HCTH/407, which are methods with either no correlation or exact exchange included, respectively. Two other pure DFT functionals, namely BLYP and M06-L, similarly yield spuriously large asynchronicities for the less activated systems. It thus seems like a balanced treatment of both exchange and correlation is vital to an accurate and consistent treatment of Diels–Alder TS geometries. Notably, the pure functionals B97D and mPW1PW91 render descriptions more in line with the hybrid functionals.

Among the family of Becke hybrid functionals, B3LYP consistently gives the largest asynchronicity. Somewhat surprisingly, the discrepancies with respect to CCSD are largest in the least activated systems; the B3LYP and CCSD asynchronicities are in fact very close for 1 + 5. Nevertheless, a comparison of the related hybrid functionals allows us to discuss the influence of exchange and correlation functionals, respectively. First, we note that B3P86 and B3PW91 give similar asynchronicities as the PBE1PBE and mPWPW91 functionals. Hence, the main problem seems to lie in the LYP correlation and not in the choice of exchange functional. Second, B2PLYP and B2PLYPD produce results similar to and consistent with the mentioned functionals, and one can surmise that inclusion of perturbative correlation energy offsets the problems caused by LYP.

M06-2X reproduces the CCSD geometries well, in stark contrast to the meta-GGA functionals M06-L and HCTH/407. The functionals ωB97X and ωB97XD by Head-Gordon and coworkers are also designed to include long-range corrections, and yield consistent results for all systems. The effect of an empirical dispersion term can be further checked by comparing the ωB97X/ωB97XD and B2PLYP/B2PLYPD geometries. Although the difference between ωB97X and ωB97XD is marginal, it is interesting to note that the asynchronicity grows faster with ωB97X. The dispersion contribution is even smaller between B2PLYP and B2PLYPD, which indicates that an adequate correlation representation diminishes the importance of the dispersion term.

Since all methods with a large extent of HF exchange (including HF!) come fairly close to the CCSD reference for systems including 2 and 3, but much less so for 4 and 5, it seems like electron correlation effects become increasingly important for TSs with a large amount of charge-transfer. For example, whereas M06-2X (54% HF exchange) yields smaller asynchronicities than PBE1PBE (25% HF exchange) for the two former systems, PBE1PBE gives the smaller value for 4 and 5. The trend is repeated in the catalyzed systems, but here the two functionals' relationship is reversed earlier due to increased activation.

At this point, it is clear that it is difficult to determine what is the best ‘recipe’ for Diels–Alder TS optimization. Specifically, the method closest to the CCSD reference varies across the series of systems. It could be that the Δd definition of asynchronicity is suboptimal. One could instead use the ‘relative asynchronicity’, which we define as Δd/(d_α + d_β), because it would take into account the fact that methods developed with medium- and long-range interactions in mind render TSs somewhat ‘earlier’ on the reaction coordinate. An early TS with a Δd of 1 Å would thus be less asynchronous than one with the same Δd but shorter distances. However, we find that there is no change in the apparent differences between the methods provided by this metric (see Fig. S1, ESI†). It can furthermore be deduced from Tables 1 and 2 that the main reason for the differences in asynchronicity is d_α. d_β does not differ much between any DFT method and CCSD, whereas d_α varies considerably across the series of reactants and methods. An illustration of this phenomenon is provided in Fig. 4, which shows all TSs of three systems superimposed on the CCSD geometry.

Fig. 4 Transition states of 1+3, 1+5 and 1·W+4, optimized using all methods discussed in this work, and superimposed on the CCSD/6-31G+(d) geometry. Note the larger discrepancies in the C₁ and C_α atom positions.

As an additional measure of geometric similarity, we determined the carbon-only root-mean-square deviations (RMDSs) between all methods for the systems that were optimized using CCSD. Two-dimensional matrix representations of the RMSD calculation are given in Fig. S2 (ESI†), from which it appears that the methods most consistently similar to the overall CCSD geometries are ωB97X, ωB97XD, B2PLYP, B2PLYPD and M06-2X. In addition, no method reproduces the CCSD incipient carbon–carbon distances better than M06-2X (cf.Tables 1 and 2). With a mean absolute deviation of 0.034 Å for the ten distances present at the CCSD level, it precedes ωB97X at 0.046 Å. In addition to the central Δd metric, these two measures corroborate the initial assumption that M06-2X is well suited for Diels–Alder optimizations – along with the less widely spread Head-Gordon functionals.

Our conclusions so far are consistent with those obtained in a recent study by Fleurat-Lessard and coworkers, who studied a nucleophilic aromatic substitution and a nucleophilic substitution to a carbonyl using a range of density functionals.⁴⁵ The two stepwise reactions form charge-separated intermediates with a π → σ transition, which is reminiscent of what happens in a Diels–Alder TS.⁷ The authors found that a large amount of HF exchange was required to reproduce the stepwise mechanism, and that the LYP correlation functional was a poor choice for these reactions.

The systematically lower estimation of TS asynchronicity by MP2 is consistent with Fig. 2. It has further been argued that MP2 is inferior to M06-2X in predicting π–π stacking interactions,⁴⁶ and it is well-known that MP2 tends to underestimate Diels–Alder activation energies (cf.Table 3).⁶ On the other hand, the spin-component-scaled MP2 method (SCS-MP2)⁴² has been shown to yield realistic energies of pericyclic reactions,^6,7 and has been used by us^12,14 and others⁴⁷ for accurate energy calculations. It thus appears motivating to employ SCS-MP2 to calculate more realistic TS geometries. However, it does not seem like the method has been used extensively for geometry optimization apart from a few tests of non-covalent interactions.⁴⁸

Table 3 Calculated activation energies using the various methods at the level of optimization^a

Method	1+2	1+3	1+4	1+5	1·W+2	1·W+3	1·W+4	1·W+5
a All energies are given in kcal mol^−1. b SCS-MP2 energies recalculated from the MP2/6-31+G(d,p) geometries. c Computed as a single point CCSD(T)/6-31+G(d) energy on the CCSD/6-31+G(d) geometries, with a MP2/6-311++G(2df,2pd) basis set extrapolation. d SCS-MP2 energies obtained from the MP2/6-311++G(2df,2pd) single point calculations on the CCSD geometries.
HF	43.7	43.3	43.4	33.9	41.5	40.8	39.4	30.3
BLYP	18.9	17.9	17.5	10.2	16.5	14.9	15.3	7.4
B3LYP	20.4	19.5	19.8	12.0	17.8	16.5	17.4	9.2
B3P86	14.6	14.0	14.6	7.3	12.0	11.1	12.3	5.2
B3PW91	17.2	16.6	17.0	10.0	15.0	14.1	14.9	7.8
B97D	9.5	8.3	7.9	0.8	6.2	5.1	5.3	−2.0
ωB97X	18.8	17.9	19.0	10.7	15.5	15.0	15.9	7.8
ωB97XD	15.2	14.1	14.8	6.9	11.5	10.7	11.6	3.7
PBE1PBE	14.2	13.6	14.3	7.2	11.6	10.9	12.0	5.0
mPW1PW91	15.9	15.3	15.9	8.9	13.5	12.8	13.7	6.7
HCTH	20.4	18.9	18.0	11.4	17.9	16.3	15.7	8.3
M06-2X	13.8	13.2	14.8	6.6	10.2	10.6	11.9	3.6
M06-L	11.0	10.0	9.8	1.8	7.5	7.3	7.2	−1.3
B2PLYP	17.5	16.4	17.2	9.8	14.7	13.8	14.6	7.3
B2PLYPD	10.8	11.4	12.2	4.7	7.4	8.5	9.3	1.7
SCS-MP2^b	18.4	17.0	19.0	12.7	15.7	15.0	16.5	9.9
MP2	10.4	8.9	10.8	5.3	7.5	6.8	8.3	2.5
CCSD(T)^c	15.1	14.2	16.1	9.0	12.3
SCS-MP2^d	16.2	15.5	18.2	11.5	13.8

---

Hence, we optimized TS geometries using MP2 and SCS-MP2 with the RICC2 module in Turbomole.⁴³ The incipient bond distances resulting from these calculations are provided in Tables 1 and 2, and a graphical comparison of the Turbomole TS asynchronicities is given in Fig. 5. The Turbomole MP2 geometries are, within an acceptable error, identical to those produced using Gaussian09. Encouragingly, the SCS-MP2 optimization increases the asynchronicity compared to MP2. The Δd values are however still significantly lower than for most methods in Fig. 3. From Tables 1 and 2 it may be concluded that the main difference between MP2 and most other geometries is that the former predicts comparatively large d_β, which in turn leads to more synchronous TSs. SCS-MP2 mitigates this phenomenon, but also lowers d_α with respect to the DFT functionals. It is clear, however, that the spin-component scaling does provide TS geometries more similar to those derived from CCSD and ‘well-behaved’ DFT methods. This indicates that the balance in correlation treatment between parallel and antiparallel electron pairs is a key factor for Diels–Alder TS geometries.

Fig. 5 Transition state asynchronicity per system using MP2 (red) and SCS-MP2 (blue) optimization in Turbomole and the 6-31+G(d,p) basis set. The thin black and grey bars represent geometries produced using the TZVPP and TZVPP+basis sets, respectively (see the text).

Post-Hartree–Fock methods, such as CCSD and MP2, are sensitive to the choice of basis set. However due to limitation in computational resources, we were not able to optimize geometries with a larger basis set than 6-31+G(d) for the CCSD method. Instead we decided to explore the basis set dependence on the MP2 and SCS-MP2 geometries. We optimized all geometries using the Alrichs-type TZVPP basis set.⁴⁹ This is larger than 6-31+G(d,p) and features a (2df,2pd) polarization but lacks diffuse functions. We therefore also investigated the effect on the geometries by adding the diffuse functions from 6-31+G(d,p) to TZVPP; we refer to this basis set as TZVPP+. It can first be noted that there are very small differences in geometries between TZVPP and TZVPP+ (Fig. 5, black and grey bars, respectively), indicating that the geometries are nearly converged with respect to the basis set. Furthermore, we note that the asynchronicities (Δd) are consistently reduced by 0.08–0.12 with the larger basis sets for both MP2 and SCS-MP2; increasing the basis set does not lead to a relatively larger decrease of the asynchronicity of the more activated systems, rather the opposite. It can therefore be anticipated that the relative trend of the CCSD/6-31+G(d) geometries with strongly increasing asynchronicities along the series would hold also at the basis set limit, although their magnitudes would supposedly be somewhat smaller. This strengthens the conclusion that the DFT methods in general overestimate the asynchronicity of the less activated systems while providing more reliable results for those that are most strongly activated.

We finally note that several DFT methods predict activation energies in good comparison with the CCSD(T) level (Table 3) without necessarily reproducing the geometries. The B3P86, ωB97XD, PBE1PBE, mPW1PW91 and M06-2X energies are mostly within 1.5 kcal mol⁻¹ from the CCSD(T) values, while we observe the well-known overestimation by B3LYP. The empirical dispersion term seems to result in underestimated activation energies. Moreover, recalculating the usual underestimations by MP2 using the SCS-MP2 protocol gave slightly overestimated energies. It is normally affordable to improve the energy calculations using larger basis sets; Table 3 merely shows how accurate an estimate one can expect to obtain at the level of optimization. We note that the MP2/6-311++(2df,2pd) level energies, used for the basis set extrapolation of CCSD(T), effectively reduce the MP2 barriers by little over 3 kcal mol⁻¹. Applied with the SCS-MP2 protocol, the energies instead become overestimated by ca. 1.5 kcal mol⁻¹ with respect to CCSD(T). Interestingly, there appears to be no strong dependence on the geometry for the energy prediction of the DFT methods (comparing Table 3 with Fig. 3). This point manifests our previous note on the flatness of the potential energy surfaces, which seems to be a general feature of Diels–Alder TSs.²⁶

4 Conclusions

This study highlights the challenges one faces when optimizing Diels–Alder transition states using some quantum chemical or DFT approaches. First, the DFT geometries are sensitive to both the choice of exchange and correlation functionals. Second, it is difficult to establish a reliable benchmark since the post-Hartree–Fock results are sensitive to the basis set composition. On the basis of the results from the MP2 and SCS-MP2 optimizations with large basis sets, we believe that the CCSD/6-31+G(d) asynchronicities are somewhat exaggerated. Despite the benchmark uncertainty, it is clear that the Diels–Alder reaction is an exception from Simón and Goodman's conclusion¹⁸ that B3LYP reproduces geometries well, particularly for moderately activated systems. B3LYP is however not by far the most deviant method compared to CCSD.

Bearing in mind that the Diels–Alder reaction mechanism is thought to transit continuously^13,21 from concerted synchronous to a polar asynchronous, and in some cases a stepwise mechanism with increasing activation,^14,22 the results presented here suggest that predictions about such mechanistic distinctions should be made carefully. Perhaps the safest recommendation that can be made from this study is that one should use more than one method for geometry optimizations to obtain a consensus description of the reaction. It may well be that some of the generally accepted aspects regarding the Diels–Alder mechanisms, in particular those derived from the large number of studies employing B3LYP, need reconsideration.

While the M06-2X functional certainly seems to be advantageous over B3LYP for geometry optimizations of Diels–Alder TSs, the difference between the two methods decreases down the series of increasingly activated reactants. It appears that the choice of exchange functional is most critical for less activated systems, while the choice of correlation becomes crucial in the more activated examples. Specifically, the relative difference between B3LYP and M06-2X decreases with more activation (note also that for TSs containing 4 and 5, the asynchronicities of the deviant BLYP and M06-L methods are virtually in line with B3LYP and M06-2X, respectively). However, for the more polar TSs in our previous studies, B3LYP and M06-2X again yield quite different geometries (cf.Fig. 2).^14,17 Functionals containing P86 and PW91 correlation are well-behaved across the series of reactants, and also seem to provide good estimates of barrier heights. B2PLYP produces decent geometries and energies, but at a higher computational cost.

In conclusion, M06-2X appears to be a better choice for Diels–Alder TS optimization than B3LYP, especially when weighing absolute distances and RMSD values. However, this study shows that there are several other functionals which seem to be equally apt to perform this task.

Acknowledgements

This work has been supported by the Cambridge Crystallographic Data centre (CCDC).

References

1. P. J. Stephens, F. J. Devlin, C. F. Chabalowski and M. J. Frisch, J. Phys. Chem., 1994, 98, 11623–11627.
2. (a) A. D. Becke, Phys. Rev. A: At., Mol., Opt. Phys., 1988, 38, 3098–3100; (b) A. D. Becke, J. Chem. Phys., 1993, 98, 1372–1377; (c) A. D. Becke, J. Chem. Phys., 1993, 98, 5648–5652.
3. C. Lee, W. Yang and R. G. Parr, Phys. Rev. B: Condens. Matter Mater. Phys., 1988, 37, 785–789.
4. P. H.-Y. Cheong, C. Y. Legault, J. M. Um, N. Celebi-Ölcüm and K. N. Houk, Chem. Rev., 2011, 111, 5042–5137.
5. P. Schreiner, Angew. Chem., Int. Ed., 2007, 46, 4217–4219.
6. T. P. M. Goumans, A. W. Ehlers, K. Lammertsma, E.-U. Würthwein and S. Grimme, Chem.–Eur. J., 2004, 10, 6468–6475.
7. S. Pieniazek, F. Clemente and K. Houk, Angew. Chem., Int. Ed., 2008, 47, 7746–7749.
8. Y. Zhao and D. G. Truhlar, Theor. Chem. Acc., 2008, 120, 215–241.
9. B. H. Kirk and D. H. Ess, Tetrahedron Lett., 2011, 52, 1245–1249.
10. R. S. Paton, S. Kim, A. G. Ross, S. J. Danishefsky and K. N. Houk, Angew. Chem., Int. Ed., 2011, 50, 10366–10368.
11. S.-Y. Tang, J. Shi and Q.-X. Guo, Org. Biomol. Chem., 2012, 10, 2673–2682.
12. M. Linder, A. J. Johansson and T. Brinck, Org. Lett., 2012, 14, 118–121.
13. H. V. Pham, D. B. C. Martin, C. D. Vanderwal and K. N. Houk, Chem. Sci., 2012, 3, 1650–1655.
14. M. Linder and T. Brinck, J. Org. Chem., 2012, 77, 6563–6573.
15. X.-F. Gao, C.-X. Cui and Y.-J. Liu, J. Phys. Org. Chem., 2012, 850–855.
16. Y. Zhao and D. G. Truhlar, Chem. Phys. Lett., 2011, 502, 1–13.
17. M. Linder, A. J. Johansson, B. Manta, P. Olsson and T. Brinck, Chem. Commun., 2012, 48, 5665–5667.
18. L. Simón and J. M. Goodman, Org. Biomol. Chem., 2011, 9, 689–700.
19. S. Kong and J. Evanseck, J. Am. Chem. Soc., 2000, 122, 10418–10427.
20. C. Møller and M. S. Plesset, Phys. Rev., 1934, 46, 618–622.
21. L. R. Domingo and J. A. Saéz, Org. Biomol. Chem., 2009, 7, 3576–3583.
22. L. R. Domingo, J. Andrés and C. N. Alves, Eur. J. Org. Chem., 2002, 2557–2564.
23. D. M. Birney and K. N. Houk, J. Am. Chem. Soc., 1990, 112, 4127–4133.
24. J. Garcia, J. Mayoral and L. Salvatella, J. Am. Chem. Soc., 1996, 118, 11680–11681.
25. (a) C. Gonzalez and H. B. Schlegel, J. Chem. Phys., 1989, 90, 2154–2161; (b) C. Gonzalez and H. B. Schlegel, J. Phys. Chem., 1990, 94, 5523–5527.
26. M. Linder and T. Brinck, Org. Biomol. Chem., 2009, 7, 1304–1311.
27. J. B. Siegel, A. Zanghellini, H. M. Lovick, G. Kiss, A. R. Lambert, J. L. St Clair, J. L. Gallaher, D. Hilvert, M. H. Gelb, B. L. Stoddard, K. N. Houk, F. E. Michael and D. Baker, Science, 2010, 329, 309–313.
28. (a) S. Lakhdar, F. Terrier, D. Vichard, G. Berionni, N. El Guesmi, R. Goumont and T. Boubaker, Chem.–Eur. J., 2010, 16, 5681–5690; (b) D. V. Steglenko, M. E. Kletsky, S. V. Kurbatov, A. V. Tatarov, V. I. Minkin, R. Goumont and F. Terrier, Chem.–Eur. J., 2011, 17, 7592–7604.
29. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, B. Mennucci, G. A. Petersson, H. Nakatsuji, M. Caricato, X. Li, H. P. Hratchian, A. F. Izmaylov, J. Bloino, G. Zheng, J. L. Sonnenberg, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, J. A. MontgomeryJr, J. E. Peralta, F. Ogliaro, M. Bearpark, J. J. Heyd, E. Brothers, K. N. Kudin, V. N. Staroverov, R. Kobayashi, J. Normand, K. Raghavachari, A. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, N. Rega, J. M. Millam, M. Klene, J. E. Knox, J. B. Cross, V. Bakken, C. Adamo, J. Jaramillo, R. Gomperts, R. E. Stratmann, O. Yazyev, A. J. Austin, R. Cammi, C. Pomelli, J. W. Ochterski, R. L. Martin, K. Morokuma, V. G. Zakrzewski, G. A. Voth, P. Salvador, J. J. Dannenberg, S. Dapprich, A. D. Daniels, O. Farkas, J. B. Foresman, J. V. Ortiz, J. Cioslowski and D. J. Fox, Gaussian 09, Revision A.02, 2009.
30. J. P. Perdew, Phys. Rev. B: Condens. Matter Mater. Phys., 1986, 33, 8822–8824.
31. J. P. Perdew, K. Burke and Y. Wang, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 16533–16539.
32. A. Becke, J. Chem. Phys., 1997, 107, 8554–8560.
33. J.-D. Chai and M. Head-Gordon, Phys. Chem. Chem. Phys., 2008, 10, 6615–6620.
34. S. Grimme, J. Chem. Phys., 2006, 124, 034108.
35. (a) S. Grimme, J. Comput. Chem., 2004, 25, 1463–1473; (b) S. Grimme, J. Comput. Chem., 2006, 27, 1787–1799.
36. Y. Zhao and D. G. Truhlar, J. Chem. Phys., 2006, 125, 194101.
37. A. D. Boese and N. C. Handy, J. Chem. Phys., 2001, 114, 5497–5503.
38. C. Adamo and V. Barone, J. Chem. Phys., 1998, 108, 664–675.
39. (a) J. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868; (b) C. Adamo and V. Barone, J. Chem. Phys., 1999, 110, 6158–6170.
40. (a) G. E. Scuseria, C. L. Janssen and H. F. Schaefer III, J. Chem. Phys., 1988, 89, 7382–7387; (b) G. E. Scuseria and H. F. Schaefer III, J. Chem. Phys., 1989, 90, 3700–3703.
41. P. Salvador, M. Duran and X. Fradera, J. Chem. Phys., 2002, 116, 6443–6457.
42. S. Grimme, J. Chem. Phys., 2003, 118, 9095–9102.
43. (a) R. Ahlrichs, M. Bär, M. Häser, H. Horn and C. Kölmel, Chem. Phys. Lett., 1989, 162, 165–169; (b) TURBOMOLE, v. 6.0, 2009, http://www.turbomole.com.
44. (a) C. Hättig, Phys. Chem. Chem. Phys., 2005, 7, 59–66; (b) F. Weigend, A. Köhn and C. Hättig, J. Chem. Phys., 2002, 116, 3175–3183.
45. N. Cheron, D. Jacquemin and P. Fleurat-Lessard, Phys. Chem. Chem. Phys., 2012, 14, 7170–7175.
46. J. Gu, J. Wang, J. Leszczynski, Y. Xie and H. F. Schaefer III, Chem. Phys. Lett., 2008, 459, 164–166.
47. I. Held, A. Villinger and H. Zipse, Synthesis, 2005, 1425–1426.
48. (a) R. A. Distasio Jr. and M. Head-Gordon, Mol. Phys., 2007, 105, 1073–1083; (b) T. Takatani and C. David Sherrill, Phys. Chem. Chem. Phys., 2007, 9, 6106–6114.
49. (a) A. Schäfer, H. Horn and R. Ahlrichs, J. Chem. Phys., 1992, 97, 2571–2577; (b) A. Schäfer, C. Huber and R. Ahlrichs, J. Chem. Phys., 1994, 100, 5829–5835.

---

Footnote

† Electronic supplementary information (ESI) available: Additional figures, including 2D-RMSD matrices of transition state geometries. See DOI: 10.1039/c3cp44319a

---