Efficient and Accurate Description of Diels-Alder Reactions using Density Functional Theory
Modeling chemical reactions using Quantum Chemistry is a widely used predictive strategy capable to complement experiments in order to understand the intrinsic mechanisms guiding the chemicals towards the most favorable reaction products. However, to do so, it is mandatory to use reliable and computationally tractable theoretical methods. In this work, we focus on six Diels-Alder reactions of increasing complexity and perform an extensive benchmark of middle-to low-cost computational approaches to predict the characteristic reactions energy barriers. We found that Density Functional Theory, using the ωB97XD, LC-ωPBE, CAM-B3LYP, M11 and MN12SX functionals, with empirical dispersion corrections coupled to an affordable 6-31G basis set, provides quality results for this class of reactions, at a small computational effort.Such efficient and reliable simulation protocol opens perspectives for hybrid QM/MM molecular dynamics simulations of Diels-Alder reactions including explicit solvation.
Introduction<!>Reactions under consideration<!>Reference calculations<!>DFT calculations<!>Wave-function and Semi-empirical Hamiltonians<!>Basis Set Effect<!>Activation energies on benchmark reactions<!>6-31G basis set performances<!>Semiempirical methods<!>Application to Strained Cyclic Allene Dienophiles<!>Conclusions<!>Supporting Material
<p>The Diels-Alder reaction, first observed in 1928 by Otto Diels and Kurt Alder, has long been established as a key transformation in the toolbox of chemists. This two-bond-forming pericyclic combination of a 1,3-diene (4π component) with a dienophile (2π component) allows the formation of six-membered rings with well-defined regio-and stereo-chemistries and tolerates a large variety of functionalities. It is therefore not surprising that its applications in synthesis are numerous, going from the access to small, highly functionalized building blocks to the construction of complex natural products or polymers. 1,2 The possibility of describing these reactions via chemical simulations is thus of crucial importance. For example, it makes possible to understand (and predict) solvation effects or to model the role of pressure in tuning stereo-selectivity. [3][4][5][6] At this end, computational chemists are often faced to deal with relatively large systems, either because the reactants are complex and extended, or because the solvent must be explicitly considered. One example of this last case, was reported recently by us showing that to understand pressure effects one needs to use molecular dynamics simulations within a mixed quantum mechanics -molecular mechanics (QM/MM) approach. [6][7][8] Also, chemical dynamics simulations are often performed to elucidate the nature of mechanisms in organic reactions, for example to distinguish between synchronous or step-wise mechanisms, 9,10 to elucidate post transition state dynamics, [11][12][13][14] or to understand kinetic isotopic effects, 15 and this technique needs a large statistical sampling to have converged and robust results.</p><p>To allow the simulations of extended systems, two approaches are tempting: (i) using small basis sets in conjunction with electronic structure, both using Density Functional Theory (DFT) or Hartree-Fock (HF); (ii) using semi-empirical Hamiltonians (SEHs). In fact, the basis set size remarkably affects the computational cost of electronic structure calculations, while large basis are required, e.g., in ab initio methods. To reach a sufficient accuracy, some DFT functionals might require very large basis sets. 16 Given this constraint, MD simulations rapidly start to be unfeasible for real life applications, where a large phase space needs to be sampled.</p><p>DFT was widely used to computationally study Diels-Alder reactions, and in a seminal work by Houk and co-workers it was pointed out that the popular B3LYP functional is not able to correctly describe single bond formation and double bond breaking. 17 Since then, the developments of range-separated functionals and inclusion of dispersion made DFT more reliable and other functionals (like ωB97XD or M06-2X) are largely used in physical organic chemistry. Later, Grimme and co-workers pointed out how results of popular B3LYP functional on a set of organic reactions can be improved adding dispersion corrections. 18 However, the search for reliable and/or computationally cheap basis sets is an important task in the field. 19 SEHs were largely used in the past and replaced by modern theories when DFT calculations became available. However, in last years research in optimizing them was still active and they are still useful to address complex chemical reactions when a huge number of calculations are needed, like is the case when studying unimolecular dissociation in conjunction with mass spectrometry. 20,21 Of course these two options do not cover all the spectrum of possibilities, for example tight-binding DFT shows recently that it is able to deal with chemical reaction dynamics. 22,23 Another possibility is also to train and/or use a machine learning algorithm for chemical reactions. 24 While some attempts to use this technique have been performed to predict activation energies, 25 it is at a preliminary stage and will need a huge data-set. Very recently, some reactive force-fields were developed with a first application also to a simple Diels-Alder reaction, but, at the present stage, this method needs a specific parametrization for each reaction. 10 Surely, this approach will have interesting future applications in the field of organic reactions, but it needs, at the present stage, further development and more grounded common and standard utilities.</p><p>Nowadays organic chemistry can maximize the efficiency of well-known and studied approaches like DFT, HF or SEHs. At this end, we investigated how accurate is the description of Diels-Alder reactions using such relatively fast methods. We have first studied a set of benchmark calculations on well-established Diels-Alder reactions and then applied to specific (4+2) cycloadditions where the dienophiles are strained cyclic allene for which both experimental and very recent calculations are reported. [26][27][28][29] Prototypical Diels-Alder cycloadditions (see reactions B1-B3 in Scheme 1) were selected for this benchmark. These small systems were chosen in order to be able to perform highly-correlated reference calculations (here at CCSD(T) level of theory). We thus chose three reactions involving cyclic or acyclic dienes, aromatic or non-aromatic substrates, and dienophiles bearing different types of classical electron-withdrawing groups (cyano, nitro and ester/anhydride). This led us to select B1 involving cyclopentadiene and acrylonitrile, a classical reaction that we had recently studied computationally under high pressure. 6 In B2 cycloaddition, butadiene reacts with nitrobenzene as aromatic dienophile, a reaction linked to our interest in (4+2) Diels-Alder cycloaddition with electron-poor aromatic compounds. [30][31][32][33][34][35] Note that in this case the cyclo-adduct will not differ by endo or exo approach, while the TS will. Here we consider the TS approach corresponding to the conformation which is lower in energy (the aromatic ring and the diene on the same side). Finally, B3 reaction involves the classical Diels-Alder between furan and maleic anhydride to form norcantharidin, whose endo/exo diastereoselectivity under kinetic/thermodynamic control has been debated over years. [36][37][38] Recently, Houk and co-workers have studied computationally a class of Diels-Alder reactions in which the dienophile has a cyclic allene structure. 26 They are particularly interesting because it is possible to synthesize several stereochemically-reach products in this way. This is possible thanks to the high endo:exo ratio as reported by many experiments. [27][28][29]39,40 Quantum chemistry calculations provide a key tool to understand and predict such endo:exo ratio. Our theoretical results can be compared with such experimental data in addition to the aforementioned calculations which use DFT with an extended basis set which would be computationally not affordable in further and likely QM/MM molecular dynamics simula-Scheme 1: Reactions under consideration in benchmark calculations.</p><p>tions. We have thus selected three reactions among those studied by Houk and co-workers and reported in Scheme 2.</p><p>The first reaction (R1) is one of the simplest possible: furan reacts with cyclohexa-1,2diene yielding 80% of products with an endo:exo ration of 11:1. 27 This selectivity is notable because it arises from "a new type of secondary orbital interaction that results from the near perpendicular approach of the diene to the dienophile and the twisted nature of the strained allene". 26 We have then chosen two other reactions, R2 where the reactants were slightly modified, resulting in similar endo:exo ratio (9.2:1), 28 and R3 where experimentally only the endo adduct was observed. 29 The aforementioned calculations reported recently, done with the ωB97XD and a relatively large basis set (6-311+G(d,p)) in implicit solvent, found a relatively small endo:exo ratio (3.3:1) 26 which is even smaller than what obtained for R1 and R2 reactions. Note that calculations were done with R=Me (methyl group), while experiments were done with R=Bn (benzyl group). The use of a smaller basis set would allow to easily check if the use of the benzyl group as in experiment in conjunction with the best performing functionals can bridge this gap. Scheme 2: Diels-Alder reactions with a strained cyclic allene dienophile considered in this study. R1 is done in THF, while R2 and R3 in acetonitrile (MeCN).</p><!><p>The reactions studied are reported in Schemes 1 and 2, where the nomenclature used here and hereafter for the different reactions is also depicted. All calculations described in the following were done with Gaussian 16 41 and Mopac16 42 (for semi-empirical Hamiltonians) software packages.</p><!><p>Reference energies have been obtained applying the following recipe: pre-complex (PreC) and transition state (TS) structures are optimized using the best level available, given the size of the system; the CCSD(T) energy at the Complete Basis Set (CBS) limit is then obtained performing a two-points extrapolation procedure on the correlation energy only, 43 using those fixed structures. The largest affordable basis sets for the two points extrapolation have been used, typically the cc-pVDZ and cc-pVTZ Dunning's basis sets. The basis set superposition error (BSSE) counterpoise correction is applied in the CCSD(T) total energy calculations, from which a BSSE-corrected correlation energy is extracted. Just for Reaction B1 the calculations were affordable also using the aug-cc-pV(D,T)Z basis sets, and the resulting extrapolated CBS energy has been reported in the SI (see Table S2).</p><p>The final CCSD(T)/CBS reference energy is obtained adding such BSSE-corrected, extrapolated correlation energy to a reference state energy computed at the HF/aug-cc-pV5Z level, 44 so to reasonably account for any source of error coming from incomplete basis sets, including the BSSE also for the non-correlation part of the total electronic energy. This is shown to be a reasonable assumption for Reaction B1, where HF energies at aug-cc-pVQZ, aug-cc-pV5Z and aug-cc-pV6Z are compared, showing a quite low impact of basis set related errors already at the quadruple-ζ (see SI, Section S1) , and a quite small variation passing from quintuple-ζ to 6-ζ (see SI, Table S1). In Section S1 of the SI more details are reported on this topic and on the reference energies calculations in general.</p><!><p>We have used different functionals, listed in Table 1 where we classified them by type and where we specified if dispersion is also added and with which method. For dispersion, we used the D3BJ correction 45 when available for the given functional otherwise we used D2 46 or D3, 47 as specified in the same Table 1.</p><p>We used different basis sets, and namely: STO-3G, 3-21G, 6-31G, 6-31G(d,p) and 6-311++G(d,p). BSSE was calculated using the counterpoise method of Boys and Bernardi. 87 In the case of reactions involving strained cyclic allene dienophiles, we performed also calculations in implicit solvent using the SMD solvation model 88 as in the recent work by Houk and co-workers. 26 Table 1: Functional used in the benchmarking. We classified them in the Type column using the following abbreviations: local spind density approximation (LSDA), generalized gradient approximation (GGA), hybrid GGA (HGGA), meta-GGA (mGGA), hybrid-meta-GGA (HmGGA), range-separated hybrid (RS), non-separable gradient approximation (NGA), double hybrid (DH). When dispersion correction is added this is reported in the specific column, in parenthesis if the functional is used both with and without dispersion.</p><!><p>We have also considered Hartree-Fock (HF) and MP2 methods to compare with DFT and reference CCSD(T) calculations, varying also in those cases the basis set size, in the benchmarking calculations. Furthermore, we have tested some popular semi-empirical Hamiltonians (SEHs) which are computationally much faster than DFT and are thus tempting as method, in particular in molecular dynamics simulations. At this end we have considered the following methods: AM1, 89 AM1-D, RM1, 90 RM1-D, PM3, 91 PM6, 92 PM6-D, PM6-D3, PM6-D3H4 93 and PM7, 94 where D and D3 stand for the D 95 and D3 47 methods by Grimme, respectively.</p><!><p>We first consider in detail the role of the basis set superposition error (BSSE) in particular concerning energy barriers. This will be important to understand also the role of the intrinsic error of a finitie basis set. DFT calculations are here compared with CCSD(T)/CBS results which are taken as reference for both method and basis set.</p><p>First we consider the activation energy of reactions B1 and B2. The activation energy</p><p>is calculated considering as the starting point the non-covalent complex formed by the reactants (PreC). Note that the PreC is obtained once from an intrinsic reaction coordinate (IRC) calculation from the TS structure corresponding to the reactant state and the resulting structure is then directly optimized each time with the different methods.</p><p>For Reaction B1, which is the simplest system treated in this work, additional high-level electronic structure calculations are performed, to estimate the reliability of the reference chosen. We notice that, as reported in Table S2 of the Supporting Information (SI), the CCSD(T) calculations with the larger basis sets used for the CBS extrapolation are almost identical to those which use the slightly smaller basis sets. We then use this last one extrapolation for the other reactions with good confidence.</p><p>It is relevant to observe here that the MP2/CBS (see Section S1 of SI for details) barrier provided in the SI is roughly half of the CCSD(T)/CBS one. Since the latter method is widely considered as a golden standard, the reason of such a huge mismatch can be ascribed to an erroneous prediction of MP2, which can be rooted in i) the more approximated treatment of correlation energy, particularly evident ii) in its known deficiency of describing accurately non-covalent interactions.</p><p>We present now how BSSE behaves for some functionals considered: we have chosen some between those who perform the best as discussed in the following subsections. Also other DFT functionals perform as well as the ones reported hereafter, and the exhaustive list of computed energy barriers including all tested functionals is reported in the SI (see Tables S5 and S6).</p><p>In Figure 1  This is due to the, by definition, always negative BSSE value, which stabilizes more the TS, for which BSSE values are larger (see Figure 1, BOTTOM panel) than the PreCs, for which the BSSE is smaller (same panel in Figure 1). For comparison, we also reported the BSSE corrected barrier differences (empty symbols), which are, as expected, always higher than the corresponding non-BSSE corrected ones. For the larger 6-311++G(d,p) basis, the difference between non-BBSE corrected and BSSE corrected barriers is obviously much smaller, so that for PreCs, compared to the smallest Pople's basis set tested, the 3-21G. This fact can be interpreted as a smaller overlap between basis from different fragments when using STO-3G.</p><p>We also remark here that the BSSE correction applied to 3-21G calculations (TOP and MIDDLE panels in Figure 1, empty symbols) brings the results of the selected functionals close, or at least much closer than the non-BSSE corrected ones, to the reference. On the opposite the STO-3G results show in almost all cases a quite bad agreement with the reference. All in all, the reduced BSSE (see BOTTOM panel in Figure 1) improves the intrinsic behavior of DFT, so that BSSE corrected and non-BSSE corrected values are generally improved in the 6-31G and the 6-311++G(d,p) basis sets over the STO-3G qnd the 3-21G. In particular for the 6-31G basis the reduced BSSE (still larger than those of 6-311++G(d,p))</p><p>compensates for the intrinsic error of the DFT functionals, surprisingly giving a very good agreement between the reference (CBS) and the non-BSSE corrected barriers.</p><p>It is worth spending few words to quantify the time needed to perform these calculations. In Table 2 we report the slow-down time factor for an SCF energy evaluation in gasphase using the 6-311++G(d,p) basis, with and without BSSE counterpoise correction, with respect to the 6-31G one. The calculations are performed for both TS and PreC structures of Reaction B1, and we also include 6-31G calculations with BSSE correction. We chose to use two of the best working functionals, as discussed later, the MN12SX and the B2PLYP including empirical dispersion, but the conclusions we draw, reasonably, are not significantly affected by this choice. For both functionals, a slow-down factor between 3 and 15 is found for the 6-31G basis with BSSE and 6-311++G(d,p) basis without BSSE correction. For the 6-311++G(d,p) basis zith BSSE correction the slowing factor is much higher in average, being ∼48 in the worst case, for the PreC treated with the B2PLYP/D3BJ method.</p><p>If one is interested in simulating the dynamics of the system, even an apparently small factor, as a 2 or 3 slower calculation, can determine a drastic reduction in the number of dynamics steps that one can afford, reducing the portion of the phase space of the system that one can explore and consequently the reliability of calculated quantities.</p><!><p>In light of this preliminary analysis on the BSSE effect on the reaction barriers computed for two of the Diels-Alder reactions in Scheme 1, we extended the set of tested functionals to find good cost-effective methods to predict a reliable energy barrier. First we have considered the gas-phase barrier for reactions B1 and B2 comparing the DFT results with reference CCSD(T)/CBS calculations. Note that, since CCSD(T)/CBS calculations provide only electronic energies, here we compare this quantity without considering any contribution from molecular vibrations (typically to free energies): they will be added when studying the reactions with strained cyclic allene dienophiles in the last subsection. As shown in the previous subsection, STO-3G basis set provides very poor results and thus we did not include it in further comparisons.</p><p>The full set of results obtained from 3-21G, 6-31G, 6-31G(d,p) and 6-311++G(d,p) basis sets with different functionals for Reactions B1 and B2 are reported in Tables S3 and S4 in the SI where we also compare with CCSD(T)/CBS values. Note that for reaction B1 we report, for many functionals excluding the Minnesota class, values with and without dispersion correction. In general, dispersion decreases the barriers: when they are overestimated without dispersion they tend to become closer to the reference value, but the improvement is not generalized.</p><p>When dealing with reaction B2, the structures of the PreC are often totally wrong if dispersion is not included. In Figure 2 we show as an example the PreC obtained at B3LYP/6-31G level of theory with and without dispersion correction. We should remind that here we define a PreC as a structure which is formed when the reactants approach to each other and technically obtained by following the intrinsic reaction coordinate (IRC) backward to reactants from the TS down to the first minimum. In other terms, the PreC by definition should have a structure where the reactants are close together "ready" to form the TS (corresponding substrates in almost parallel plans in our case dealing with a Diels-Alder cycloaddition), which in the present case is what obtained when dispersion correction is added. On the other hand, when the dispersion is not included, the geometry is totally unrealistic, with substrates lying in perpendicular plans in the present case. The same behavior is observed with almost all the functionals and basis sets, but the Minnesota families of functionals, which are able to catch, in our application, the good behavior without adding empirical dispersion. The qualitatively good description of dispersion of the Minnesota functionals can be ascribed to their specific functional form, and thanks to the detailed and extended fitting procedure applied to evaluate the parameters which are involved in their definition. 96,97 With this consideration, and verifying that PreC and TS structures do not undergo significant distorsion, we did not add dispersion correction when using Minnesota functionals. with such relatively small basis set. Note that also other Minnesota functionals are able to provide good results, in particular for Reaction B2. Notably, when we just add polarization functions, the situation is not improved, and the same classes of functionals provide the best results. The 6-31G(d,p) basis set does not improve the quality of the results while it is computationally more expensive, even if slightly. Note that this is particularly important for molecular dynamics simulations when many energy and gradient calculations are needed and thus this can result in a huge difference from a computational point of view.</p><p>The same picture arises with a larger basis set, 6-311++G(d,p), for which again rangeseparated functionals (in particular LC-ωPBE, M11 and MN12SX) provide the barriers in better agreement with the reference values.</p><p>We now move to results for reaction B3, for which we have calculated the transition states for both endo and exo forms, and thus we can compare those two values and the difference between them. Values obtained with different basis sets (here we did not consider the smallest 3-21G basis which, as we have shown, provides too poor results) and functionals are reported in Before discussing the energetic of two analogous Diels-Alder reactions in solution, we should pause and consider the global behavior of all these functionals at a given basis set.</p><p>Since we have a relatively large number of functionals, we can see how they behave on average and if the best functionals have a considerable deviation from such average. We report in</p><!><p>As shown and discussed in the previous subsection, 6-31G basis set seems to be a good compromise. We now discuss in more details the performances of the different functionals with this basis set in order to provide a limited set of functionals which will be more likely to be used in real applications.</p><p>In Figure 3 we report the absolute error obtained by the different functionals for the barrier of reaction B1 with respect to CCSD(T)/CBS results. Notably many functionals provide an error less than 1 kcal/mol with respect to the reference, and notably the double hybrid B2PLYP/D3BJ and the range separated ωB97XD, M11, M11L, LC-ωPBE/D3BJ, MN12SX, MN12L while CAM-B3LYP/D3BJ has a slightly bigger error (2 kcal/mol).</p><p>As discussed previously, for reaction B2 many functionals (actually all but Minnesota  Also in this case, a good number of functionals show a small error even with the small 6-31G basis set, and notably the double-hybrid B2PLYP/D3BJ and the range separated ωB97XD, LC-ωPBE/D3BJ, M11 and MN12SX (as previously also other functionals show good results, like M062X and SOGGA11X). As for reaction B1, CAM-B3LYP/D3BJ provides a slightly worst energy barrier (1.9 kcal/mol).</p><p>Finally, in Figure 5 we report the activation energy barrier for endo TS of reaction B3 as well as the difference with respect to exo TS, namely a |∆E| a (endo-exo), in absolute value, for the series of tested functionals with 6-31G basis set. For comparison, in the same Figure , we also show the reference CCSD(T)/CBS results.</p><p>Surprisingly, many functionals are able to correctly catch the small |∆E| a (endo-exo) Figure 4: Absolute activation energy differences (|∆E| a ) between tested functionals (Density Functional Approximations, DFA on the y-axis) and the reference CCSD(T)/CBS, reported in kcal/mol, for 6-31G calculations for Reaction B2.</p><p>(∼0.5 kcal/mol in the reference). On the other hand, few functionals are able to provide the endo energy barrier close to the reference value. Notably, τ HCTH barrier is very close to the reference, but this method provides, wrongly, a lower exo barrier. Other well behaving functionals are SOGGA11X, M06L, CAM-B3LYP/D3BJ and ωB97XD. Still in acceptable agreement with the reference, but slightly worse, are results obtained by the Minnesota MN12SX and the N12SX functionals, which give very similar results as the B2PLYP/D3BJ.</p><p>Overall the results on the B3 reaction show less functionals in very good (close to 1.0-0.5 kcal/mol) agreement with the reference compared to the B1 and B2 cases, if we consider the absolute reaction barrier. However, concerning the difference between endo and exo barriers the general agreement is more satisfactory, as discussed in the previous section.</p><p>Figure 5: Absolute activation energies (E a in blue) and differences between endo end exo activation energies (|∆E| a in red) for the B3 reaction, are reported in kcal/mol and computed with the many functionals reported on the x-axis always using the 6-31G basis set; the reference CCSD(T)/CBS corresponding values are also reported for comparison.</p><!><p>In principle, Hartree-Fock (HF) and Semi-Empirical Hamiltonians (SEH) methods can provide an efficient way to evaluate reaction energetics, thanks to the much more favorable scaling with the size of the system if compared to MP2. We thus tested a number of SEHs and also HF with small basis sets. The full set of results are reported in Table S10 of SI, where we also report MP2 values with relatively small basis sets and reference calculations for the three benchmark reactions. In Figure 6 we show the results for SEHs and HF/6-31G calculations. As one can clearly notice from these results, the different methods show a poor agreement with respect to reference CCSD(T)/CBS calculations. For Reaction B2, we found that few SEHs, as the AM1-D, RM1, PM6-D3 and PM6-D3H4, reproduce quite well the reference, which is not the case for Reaction B1 nor B3, being probably due to fortuitous error corrections. Globally, a slight preference for RM1 or PM7 could be drawn, but the deviation from the reference is too large to be used with confidence. Thus, we do not consider in the following neither HF nor SEH methods. Note that, HF in general shows probably the worst results between all the methods tested in this work.</p><!><p>Given the results on benchmark reactions, we finally found that range separated functionals (CAM-B3LYP, LC-ωPBE, ωB97XD, M11 and MN12SX) show the best performances with 6-31G basis set, having computing time which are similar to standard functionals. Also the double-hybrid functional B2PLYP provides very good results, but it is much more computationally demanding, in particular when frequency calculations are needed and thus it will not be considered here. We have thus consider these range separated functionals plus the recently proposed SOGGA11X functional, which provides good performances with 6-31G basis set, to study some selected Diels-Alder reactions with a strained cyclic allene dienophile (see Scheme 2).</p><p>These reactions were recently studied by Houk and co-workers by using ωB97XD functional with extended 6-311++G(d) basis set in THF (reaction R1) and acetonitrile (MeCN)</p><p>(reactions R2 and R3). 26 Here we used the same implicit solvation model (SMD) with the aforementioned functionals with 6-31G basis set. Furthermore the final endo:exo ratio based on free energy TS barriers of the two forms is with experiments reported on these reactions. [27][28][29] We first report in the top panel of Figure 7 the difference of the free energy barriers between endo and exo approaches, as obtained by present calculations with 6-31G basis set and compared with what was reported by Houk and co-workers. In Table 4 we list all the associated values. We first notice that the basis set has a small effect (as noticed in benchmark calculations), in fact ωB97XD calculations done with the 6-31G basis set are very similar to what reported by Houk and co-workers with the same functional but larger basis set. Furthermore, also other functionals report values quite close to those ones, where CAM-B3LYP shows a slightly larger difference.</p><p>Calculations reported with the extended basis set are done with the methyl group for reaction R3 (R=Me in Scheme 2) while experiments were done with a benzyl group (R=Bn in Scheme 2). We have thus performed calculations with 6-31G basis set with both R=Me and Bn and the corresponding reactions are noted as R3-model and R3-real. Free energy barriers and endo-exo differences vary relatively little moving from Me to Bn group, and we will discuss more in details this aspect comparing our results directly with experimental data.</p><p>In fact, the experimental endo:exo ratios were reported for these reactions. We have calculated them simply from the transition state free energy barriers of endo and exo forms.</p><p>Results obtained from experiments and calculations are reported in the bottom panel of Figure 7. We should recall here that, from the experimental results, only the endo form is obtained for R3. This experimental finding is translated into an endo:exo ratio larger than 20.</p><p>Results show that for reactions R1 and R2 the different functionals with 6-31G basis set are in quite good agreement with experiments (and previous calculations). Notably, for reaction R1, where ωB97XD/6-311++G(d) calculations overestimate the endo ration, our results are better or equivalent (with a much lower computational cost). Moving to reaction R2, experimentally a slight decrease of endo:exo ratio is noticed experimentally, while almost all calculations show the opposite (as well as the previously reported ωB97XD/6-311++G(d) ones), with the exception of M11, which however overestimates this ratio for both reactions.</p><p>We should notice that the endo:exo ratio is very sensitive to small modification in free energy difference, since it has an exponential dependence on ∆∆G.</p><p>For reaction R3, in the previous work of Ramirez et al. the methyl group was considered and the endo:exo ratio was considerably lower than 20 (even lower than the other reactions).</p><p>We tested for reaction R3 both methyl (called here and hereafter R3-model) and benzyl group (called here and hereafter R-real), this last being the same as in experiments. Results</p><p>show that in many cases the endo:exo ratio is increased with respect to R3-model. This is particularly evident for LC-ωPBE/D3BJ and MN12SX results which are best in agreement with experiments with an endo:exo ratio larger than 20.</p><!><p>The aim of the present work is to find a reliable and affordable electronic structure-based startegy to efficiently describe, real-life, Diels-Alder reactions in the most accurate way possible. We have thus investigated the opportunity of applying reduced-cost electronic structure methods to study such reactions in order to predict reliable characteristic energy barriers. On a series of six reactions, we extensively analysed the performances of a number of recent and traditional DFT functionals (Table 1), with different basis sets, and SEHs, including in many cases empirical dispersion corrections, due to the known deficiencies of such approximated methods in treating long-range electron correlations. 96 All employed methods are commonly available in general electronic structure codes widely used in the Computational chemistry community.</p><p>In a first stage, a specific attention has been put on the role of BSSE when using DFT, testing the quality of reaction barriers predicted using different small-medium size basis sets on two relatively small reactions (reactions B1 and B2 in Scheme 1), for which BSSE effects have been explicitly analyzed by comparing the DFT results to reference CCSD(T)/CBS calculations. In most cases, a good compromise has been found using the 6-31G Pople's basis set due to a non-negligible error compensation due to the BSSE. In particular, rangeseparated functionals with dispersion corrections provide good quality results despite the modest basis set size and therefore greatly reduce the computational cost of the simulations.</p><p>On the other hand, as one could expect, HF and SEHs are shown to be quite ineffective, even with the new SEH parametrizations and dispersion corrections.  We have then applied the best DFT methods to the prediction of three different Diels-Alder reactions in solution using the small 6-31G basis set. We found that they are able to correctly describe the experimentally observed endo:exo ratios while free energy barriers appear in agreement with computations performed with an extended basis set. 26 Furthermore, the use of such affordable computational approach enables to consider the full complexity of the most extended considered system and therefore to obtain results in better agreement with experiments than those using a slightly reduced model (used in the Literature data).</p><p>This paves the way for studying either very large systems and/or to use QM/MM dynamics in solutions for Diels-Alder reactions, achieving an high level of accuracy and thus reliable, but cost-effective, predictions.</p><p>Concluding, we suggest that the best functionals to be used in future studies on Diels-Alder reactions (in association with the 6-31G basis) should be : ωB97XD, LC-ωPBE/D3BJ, CAM-B3LYP/D3BJ, M11 and MN12SX. It will be surely interesting to understand how (and if) they can be successfully applied also to model other organic reactions.</p><!><p>In the Supporting Information we report: (i) a pdf file with additional Tables mentioned in the manuscript; (ii) a tar file with xyz structures used in CCSD(T)/CBS calculations of B1, B2 and B3 reactions.</p>