Kinetics of the Strain-Promoted Oxidation-Controlled Cycloalkyne-1,2-quinone Cycloaddition: Experimental and Theoretical Studies.

Stimulated by its success in both bioconjugation and surface modification, we studied the strain-promoted oxidation-controlled cycloalkyne-1,2-quinone cycloaddition (SPOCQ) in three ways. First, the second-order rate constants and activation parameters (ΔH⧧) were determined of various cyclooctynes reacting with 4-tert-butyl-1,2-quinone in a SPOCQ reaction, yielding values for ΔH⧧ of 4.5, 7.3, and 12.1 kcal/mol, for bicyclo[6.1.0]non-4-yne (BCN), cyclooctyne (OCT), and dibenzoazacyclooctyne (DIBAC), respectively. Second, their reaction paths were investigated in detail by a range of quantum mechanical calculations. Single-configuration theoretical methods, like various DFT and a range of MP2-based methods, typically overestimate this barrier by 3-8 kcal/mol (after inclusion of zero-point energy, thermal, and solvation corrections), whereas MP2 itself underestimates the barrier significantly. Only dispersion-corrected DFT methods like B97D (yielding 4.9, 6.4, and 12.1 kcal/mol for these three reactions) and high-level CCSD(T) and multireference multiconfiguration AQCC ab initio approaches (both yielding 8.2 kcal/mol for BCN) give good approximations of experimental data. Finally, the multireference methods show that the radical character in the TS is rather small, thus rationalizing the use of single-reference methods like B97D and SCS-MP2 as intrinsically valid approaches.

Introduction

Since its coining in 2001,[1] click chemistry has played an important role in a wide variety of research areas including organic synthesis,[2] supramolecular chemistry,[3] drug development,[4] materials science,[5] and bioorthogonal labeling.[6] Among these reactions, the copper(I)-catalyzed <span class="Chemical">alkyne–azide cycloaddition (CuAAC) has been extensively studied owing to its bio-orthogonality and high reaction rates.[7] However, the use of copper(I) species is not ideally suited for labeling living systems without compromising cell function.[8,9] Consequently, in the past decade, increasing focus has been oriented toward developing novel click reactions that do not require the use of toxic metal catalysts.[10] Among the diversity of metal-free click reactions that have been identified, the strain-promoted cycloaddition of cyclooctynes (SPAAC) has been widely studied in life and material sciences.[11,12] Despite the wide applicability of the SPAAC reaction, the quest for new and faster metal-free click cycloaddition reactions is ongoing. Especially in the past decade, the inverse-electron-demand Diels–Alder reaction (IEDDA)[13,14] and the strain-promoted oxidation-controlled cycloalkyne–1,2-quinone cycloaddition (SPOCQ)[15,16] have emerged with extremely fast kinetics. Second-order rate constants up to 400–1200 M–1 s–1 have been measured for the IEDDA, while for the SPOCQ reaction with bicyclo[6.1.0]non-4-yne (BCN) k2 = 500 M–1 s–1, respectively. These values are high compared to those obtained for the SPAAC reaction of BCN and azides (0.05–1 M–1 s–1). Such rates have also been shown to be highly dependent on the nature of the alkyne; e.g., dibenzoazacyclooctyne (DIBAC) reacts in the SPOCQ reaction with a rate constant of only 0.19 M–1 s–1, i.e., 3 orders of magnitude slower than BCN. In this regard, the ability to control the polarization of strained cycloalkynes and their counterparts has also been shown to influence the rate and chemoselectivity in 1,3-dipolar cycloaddition reactions.[17−20]

The SPOCQ reaction between <span class="Chemical">cycloalkynes and 1,2-quinones is a Diels–Alder-type reaction that was shown to be highly useful for protein conjugation[15,21] and organic synthesis.[16] Very recently, we have reported the use and usefulness of the SPOCQ reaction on a solid surface,[22] showing for the first time a surface-bound metal-free “click reaction” true to the original definition of click chemistry, i.e., including the complete conversion of all quinone groups immobilized on the surface into the SPOCQ product. While full conversion is, of course, also highly desired for solution chemistry, any nonideality can still be remedied by, e.g., chromatography. In contrast, this safeguard is not present for organic conversions within firmly bound organic monolayers on the surface, where all products—wanted or not—stay at the place where they are formed (see for an analogy with polymers, ref (23)). Therefore, such 100% reaction efficiency is crucial for surface-bound reactions.[24−26] Since thus both the bio-organic syntheses and surface modification reactions have demonstrated SPOCQ to be a highly useful metal-free click reaction, we set out to explore the reaction in more detail.

The mechanism of the SPOCQ reaction has not been investigated until now, and in this regard, computational chemistry can help us to understand chemical systems at the molecular scale. In this vein, in the past decade, theoretical calculations have been employed to better understand <span class="Chemical">1,3-dipolar azide cycloadditions of cyclooctynes[27,28] and strained alkenes,[29] the IEDDA reaction of 1,2,4,5-tetrazines with strained alkenes and alkynes,[30−32] and other metal-free click reactions.[33,34] Most of the theoretical studies performed on this area are for reasons of efficiency based on density functional theory (DFT) methods, yielding reasonable agreement with experimental data. However, in some cases, it becomes necessary to use more precise methods that describe properly electron correlation effects, as DFT methods have been proven to be somewhat limited for reactions involving the transformation of electrons in π bonds to electrons in σ bonds.[35] (For a systematic benchmarking of the performance of ab initio and DFT methods for the prediction of activation barriers in pericyclic reactions see ref (36).) Regarding alternative methods to DFT, the use of second-order Møller–Plesset perturbation theory (MP2)[37] becomes a potential candidate, as MP2 is the least computationally intensive alternative method incorporating electron correlation, giving very good results for electrostatically dominated interactions and qualitatively accurate stabilization energies for dispersion-bound complexes. To its disadvantage, MP2 methods require more computational power and yield to an overestimation of the dispersion contribution to the correlation energy,[38,39] with as a result an underestimation of the transition-state energies.[40] In an attempt to compensate for these deficiencies, spin-component scaling (SCS-MP2)[41] and scaled opposite spin (SOS-MP2)[42] approaches have been developed in the past years. These two scaling factor methods have been found to predict activation barriers close to those determined by the high-accuracy G3B3 model chemistry for those cases in which single-determinant methods are known to work well. However, if the transition state starts to display significant open-shell character, perturbation methods are expected to break down, and nonperturbative methods like the coupled cluster methods (CCSD(T))[43,44] provide an interesting alternative to MP2. These have recently attracted strong, renewed attention due to an efficient program implementation in terms of the domain-based local pair natural orbital (DLPNO–CCSD(T)) approach.[45] In most demanding cases, multireference methods[46] like the multireference averaged quadratic coupled cluster (MR-AQCC) method[47] might be needed for a proper representation of especially the transition state (TS),[48] as methods like complete active space self-consistent field (CASSCF) have been shown to be insufficiently accurate.[36] The MR-AQCC method, while computationally expensive, is useful to analyze the radical character of the TS and, thus, the potential of single-reference methods to provide a proper description of the reactive system. We note that until now only the parent Diels–Alder reaction (butadiene + ethene)[36,48,49] and a relatively simple hetero-Diels–Alder reaction (butadiene + glyoxylates)[50] have been studied using multireference methods, but to our knowledge, no larger and or alkyne-based Diels–Alder reactions have been studied.

In view of the above, we report in this paper the study of the reactivity of different cyclooctyne derivatives in the <span class="Chemical">metal-free strain-promoted oxidation-controlled cyclooctyne–1,2-quinone cycloaddition (Scheme ). First, the scope of the reaction is investigated by studying the reactivity of a range of dienophiles, including cyclooctyne (OCT), BCN, and DIBAC. Second, kinetic data were obtained for these reactions at various temperatures, so as to deduce activation enthalpies. Finally, the reaction was studied theoretically by a combination of aforementioned DFT, MP2-based, multiconfiguration, and multireference methods, noting that until now no theoretical work has been published on this reaction. Here we aimed at two things: to clarify what level of theory is appropriate for the study of such metal-free cycloaddition–type “click” reactions, specifically with an eye on the question “how accurate are DFT methods?”, and to deduce detailed mechanistic information.

Scheme 1

Structures of the 1,2-Quinones (in Experiments: tert-Butyl-QUIN; in Theoretical Study: QUIN) and Different Dienophiles (Cyclooctyne Derivatives and Cyclic Alkenes) Used in This Study

Results and Discussion

Scope of the SPOCQ Reaction

Given the potential of being an extremely fast biorthogonal reaction, we explored the scope of the recently discovered SPOCQ reaction via the reactivity of the monosubstituted <span class="Chemical">4-tert-butyl-1,2-benzoquinone toward different cyclooctynes. For that, 4-tert-butyl-1,2-quinone was prepared by oxidation of 4-tert-butyl-1,2-catechol with sodium periodate and isolated as a pure solid.

We optimized our studies using the SPOCQ reaction of <span class="Chemical">cyclooctyne (OCT) and 4-tert-butyl-1,2-benzoquinone (Figure A). For that, the cycloaddition was investigated by adding 1.5 equiv of OCT to a solution of 4-tert-butyl-1,2-benzoquinone in 1,2-dichloroethane (DCE). This reaction leads to significant changes in the UV–vis spectrum due to conversion of the quinone (λmax (in 1,2-dichloroethane) = 380 nm) into a cycloaddition adduct that showed a broad band at 443 nm. This conversion yields the largest absorption differences at 384 nm (Figure B), where we therefore monitored the exponential decay of the absorbance of 4-tert-butyl-1,2-benzoquinone over time upon reaction with cyclooctyne (Figure C). On the basis of such measurements, the reaction in 1,2-dichloroethane at ambient temperature (25 °C) was determined to be complete in around 3 h. Similar studies for other cyclic octynes, such as, endo-BCN, exo-BCN, and DIBAC yielded under analogous conditions, reaction times of 0.1, 0.5, and 10 h, respectively.

Figure 1

(A) SPOCQ reaction of 4-tert-butyl-1,2-benzoquinone with cyclooctyne. (B) UV–vis spectrum (in DCE) of 4-tert-butyl-1,2-benzoquinone (red), cyclooctyne (blue) and SPOCQ product (green). (C) Variation of UV–vis absorption of a mixture of 4-tert-butyl-1,2-quinone (100 μM) and cyclooctyne (200 μM) in DCE upon SPOCQ reaction at 25 °C.

(A) SPOCQ reaction of <span class="Chemical">4-tert-butyl-1,2-benzoquinone with cyclooctyne. (B) UV–vis spectrum (in DCE) of 4-tert-butyl-1,2-benzoquinone (red), cyclooctyne (blue) and SPOCQ product (green). (C) Variation of UV–vis absorption of a mixture of 4-tert-butyl-1,2-quinone (100 μM) and cyclooctyne (200 μM) in DCE upon SPOCQ reaction at 25 °C.

Kinetic Study

The reaction-rate constant for the SPOCQ reaction of <span class="Chemical">4-tert-butyl-1,2-benzoquinone and OCT was derived to be 13 ± 2 M–1 s–1 (Figure A) at 25 °C in 1,2-dichloroethane, as determined by UV–vis measurements (see the Supporting Information for details). This reaction rate is similar to that of a typical IEDDA reaction of tetrazines with cyclooctynes[51] and higher than that of most strain-promoted azido-alkyne cycloadditions. Interestingly, the reaction rate could be slightly increased by performing the cycloaddition in methanol (51 ± 4 M–1 s–1), showing a reaction rate enhancement due to the protic solvent.

Figure 2

Determination of reaction rate constant in 1,2-dichloroethane at 25 °C measuring the decay in absorbance of 4-tert-butyl-1,2-benzoquinone at 384 nm for (A) OCT and (B) endo-BCN.

Determination of reaction rate constant in 1,2-dichloroethane at 25 °C measuring the decay in absorbance of <span class="Chemical">4-tert-butyl-1,2-benzoquinone at 384 nm for (A) OCT and (B) endo-BCN.

Next, we investigated the reactivity of endo-BCN, which was found to react much faster than its homologous <span class="Gene">OCT, as indicated by an immediate color change from orange to yellow. We found the second-order rate constant for the reaction between endo-BCN and 4-tert-butyl-1,2-benzoquinone (in 1,2-dichloroethane) to be 219 ± 14 M–1 s–1 (Figure B). This rate is 2-fold slower than that reported in methanol/water (1:1)[15] but on the same order of magnitude as the reaction of BCN with tetrazines.[52] However, when the reaction was performed in methanol, the rate constant could be increased up to 838 ± 22 M–1 s–1. On the other hand, under the same reaction conditions, exo-BCN reacted more slowly (99 ± 3 and 298 ± 17 M–1 s–1, in 1,2-dichloroethane and methanol, respectively). Analogous reaction rates for the DIBAC have been reported to be 0.12 ± 0.02 and 0.51 ± 0.06 M–1 s–1 in 1,2-dichloroethane[15] and methanol, respectively, i.e., orders of magnitude slower. While this makes DIBAC less relevant from a synthetic point of view, it presents a good test case for further theoretical studies.

The reaction rates for these three cyclooctyne derivatives were also studied at different temperatures (from 15 to 65 °C, in <span class="Chemical">1,2-dichloroethane) to extract the thermodynamic parameter of activation (ΔH⧧) by Eyring plot analysis. From the slope of the least-squares fit of plots of ln (k/T) versus 1/T, we obtained the experimental values of ΔH⧧ = 7.3 ± 0.4 kcal/mol for OCT and 4.5 ± 0.3 kcal/mol for endo-BCN. Finally, for the bulky and slower DIBAC derivative, we obtained an ΔH⧧ value of 12.1 ± 0.5 kcal/mol. For some of these reactions, the plot is not perfectly linear (see, e.g., the Figure A inset), as small amounts of byproducts are observed. We therefore regard the experimental error to be about 1 kcal/mol, i.e., slightly larger than the standard deviations given above.

The experimentally observed difference in reactivity between 1,2-benzoquinone and the set of <span class="Chemical">cyclooctyne derivatives can be qualitatively rationalized based on the HOMO and LUMO energy levels. In this regard, a closer inspection of the LUMOdienophile–HOMOquinone energy gap (from M06-2X ground-state calculations) reveals an energy gap of 10.2 eV for the reaction between cyclooctyne and 4-tert-butyl-1,2-benzoquinone. However, when endo-BCN was used as dienophile, the LUMO energy was decreased by 0.27 eV, giving this an increased reactivity, as the LUMOdienophile–HOMOquinone gap is decreased.[53] Finally, for the DIBAC compound, the LUMOdienophile–HOMOquinone gap increased as a consequence of an increase in the LUMO energy for DIBAC. Accordingly, the smaller orbital energy gap between BCN and OCT makes the favorable orbital interaction in TS2 (BCN) stronger than that in TS1 (OCT) or TS3 (DIBAC) (see also Figure ).

Figure 3

M06-2X/6-311+G(d,p)-optimized transition-state structures for the SPOCQ reaction of QUIN with OCT (TS1), BCN (TS2), and DIBAC (TS3).

M06-2X/6-311+G(d,p)-optimized transition-state structures for the SPOCQ reaction of <span class="Chemical">QUIN with OCT (TS1), BCN (TS2), and DIBAC (TS3).

Theoretical Investigations

In order to get new insights in the mechanism of this novel click reaction, we next studied theoretically the reactivity difference for the SPOCQ between <span class="Chemical">BCN and OCT with 1,2-quinone by quantum chemical means. To simplify, the monosubstituted 4-tert-butyl-1,2-benzoquinone, as used experimentally, was replaced by 1,2-benzoquinone. Comparative M06-2X/6-311+G(d,p) and B3LYP/6-311+G(d,p) calculations did not show a significant influence of the 4-tert-butyl group in the thermodynamics of the cycloaddition.

Density Functional Theory calculations

First, a variety of popular DFT[54] was tested because of the combination of their potentially high accuracy and fast computational performance. For this preliminary study, methods based on Becke’s GGA exchange functionals (BLYP,[55] B97D,[56] and B3LYP[57]) and meta-GGA class functionals (M06-2X,[58] M06L,[59] M11,[60] <span class="Mutation">M11L,[61] and MN15[62]) with the 6-311+G(d,p) basis set were used and compared. The three SPOCQ cycloadditions under study were found to be highly exothermic. As an example, the computed reaction free energies (ΔGrxn) at the M06-2X level were found to be −38.4, −38.6, and −42.1 kcal/mol for OCT, BCN, and DIBAC, respectively. [Most DFT and MP2-based calculations (see below) typically yield highly similar figures.] Typically, the computed ΔGrxn for the DIBAC reaction are 2–5 kcal/mol more exothermic than the reactions of BCN and OCT (see Table S7). In the reaction of 1,2-quinone with the cyclooctynes, the activation free energies increase from BCN to OCT, in agreement with the higher reaction rate measured for BCN (see Table S6).

In nearly all cases, the computed activation barriers were significantly higher than the experimental values, with the exception of B97D (and some ab initio methods; see the next section) which gives good agreement with experiment. Among the Becke’s GCA exchange functionals, B97D gave a better agreement (lower activation barriers) than the pure BLYP and the hybrid B3LYP functional (see Table S12). Inclusion of the dispersion energy (e.g., via D3) has been shown to be crucial in the case of intermolecular complexes or transition structures,[63] and in addition, in our study this addition is consistently improving the agreement with experiment: the computed activation energy is typically reduced by ∼1 kcal/mol in going from the BLYP or B3LYP functional to the dispersion-corrected BLYP+D3 and B3LYP+D3 calculations (Tables S1 and S2). For the meta-GGA class functionals, it was observed that M06-2X and M06L gave slightly lower activation barriers for <span class="Chemical">BCN than M11 and M11L; however, for OCT, this difference was almost negligible (0.3 kcal/mol). For the hybrid meta-NGA MN15 functional, activation energies for BCN and OCT were similar to those computed by the M06-2X functional, although for DIBAC this value was similar that obtained by M11 functional. Overall, the calculated activation energies with B3LYP and the entire set of Minnesota density functionals under study were roughly in the same range and overestimating the activation enthalpy beyond the experimental error (typical overestimation ∼4 kcal/mol), and only B97D and dispersion-corrected BLYP+D3, yielded lower activation energies for the three cyclooctynes, with B97D being the best across the board (see Figure ).

Figure 4

Comparison of DFT-calculated activation enthalpy for the SPOCQ reaction of BCN, OCT, and DIBAC with QUIN in 1,2-dichloroethane. [All calculations were performed using the 6-311+G(d,p) basis set, and values are given in kcal/mol.]

Comparison of DFT-calculated activation enthalpy for the SPOCQ reaction of <span class="Chemical">BCN, OCT, and DIBAC with QUIN in 1,2-dichloroethane. [All calculations were performed using the 6-311+G(d,p) basis set, and values are given in kcal/mol.]

The Diels–Alder reaction proceeds for <span class="Chemical">BCN with near-synchronous C–C bond formation for both new C–C bonds (both 2.32 Å, B97D data, optimized in DCE). Somewhat more dischronicity is found for OCT (2.29 and 2.38 Å), whereas for DIBAC, the TS is not symmetric at all, with bond lengths for the C–C bonds of 2.04 and 2.78 Å, respectively. While the values are numerically slightly different, the same trend was observed for M06-2X (numerical values: both 2.32 Å for BCN, 2.32 and 2.34 Å for OCT, and 2.22 and 2.35 Å for DIBAC), suggesting that this indeed represents the reactions rather than the methods. The method dependence of the dischronicity of the Diels–Alder reaction has also been investigated systematically for smaller unsymmetrical reactions.[64] It was found that, among other methods, the M06-2X approach gave the best results, which also confirms our findings.

To probe deeper into this trend in reactivity in the SPOCQ reaction, we calculated the interactions and distortion energies for the transition states.[65,66] To this end, each transition structure was separated into two fragments (the distorted <span class="Chemical">quinone and cyclooctyne compound), and B97D single-point energy calculations were performed on each of these fragments in their TS geometry to obtain the distortion energy of the quinone and cyclooctyne, calculated as the energy difference between the distorted and optimized ground-state structure of each component. In the case of BCN, the calculated distortion energy was 4.0, for OCT 5.1, and for DIBAC 7.4 kcal/mol (cf. calculated enthalpy barriers: 4.5, 6.4, and 9.3 kcal/mol, respectively). The interaction energy between the two distorted components in their TS geometry is the difference between the activation energy and the total distortion energy (sum of distortion energy of the quinone and cyclooctyne).[67] These data imply that the contributions from both geometrical distortion and TS interaction increase in going from BCN via OCT to DIBAC: in these latter two, slower reactions both terms are causing the slowdown. This can, for example, be understood from the increasing dischronicity of the two C–C bond formations but also from the C1–C2–C3–C4 dihedral angle in the TSs: for BCN this value is near zero, i.e., yielding the more favorable orbital overlap in the cycloaddition transition state, but it rises for OCT to become 8.5° for DIBAC (see Table S8 for more extensive tabulations). We thus also observe, in line with recent findings, that more dischronicity leads to both reduced distortion energies, but thus also to interaction energies.[68]

As described above, we observed in the experiments a significant solvent effect as shown by the increase in the reaction rate when moving from 1,2-dichloroethane to <span class="Chemical">methanol. Therefore, to gain more insights about this solvent effect, we also calculated the activation free energies in 1,2-dichloroethane (relative dielectric constant ε = 10.1) and methanol (ε = 32.6) using the Conductor-like Polarizable Continuum Model (CPCM) for the variety of dienophiles used in this study. As shown in Table S2, the (ΔH⧧calc)gas were typically 1–2 kcal/mol higher than (ΔH⧧calc)DCE or (ΔH⧧calc)MeOH. However, contrary to the significant solvent effect observed experimentally, no substantial differences were observed when comparing the calculated activation free energies in 1,2-dichloroethane and methanol. This small influence of the continuum solvent model on the activation enthalpies points to the limitations of this solvation method specifically for protic solvents. The calculations would require more sophisticated, combined models in which solvent molecules are included explicitly in the quantum chemical calculation,[69,70] beyond the scope of the current study.

Finally, in order to validate the practical use of DFT methods in predicting the rate of other SPOCQ reaction partners, we studied the correspondence between the B97D- and M06-2X-computed activation free enthalpies (ΔH⧧calc) and the experimental values (ΔH⧧exp). As shown in Figure , both B97D and M06-2X predict the correct reactivity pattern for the different dienophiles used in this study. However, whereas M06-2X can only be used for the more qualitative trends (and this is actually true for nearly all DFT methods under study), B97D actually yields a good agreement with experiment, at least over the range studied.

Figure 5

Plot of computed activation free energies (ΔH⧧comp) vs experimental observed values (ΔH⧧exp). Dotted line indicates x = y.

Ab Initio Calculations

DFT calculations need to be carefully used as the level of accuracy can differ significantly among different chemical systems and sometimes in a not straightforward manner (we do, for example, not observe a clear improvement in going from M06/M06-2x via M11 to MN15, whereas that was evidently the case for the large test sets used in their development).[62] Therefore, ab initio methods were also investigated, and we specifically used both single-reference methods (<span class="Gene">MP2 and variants, and DLPNO-CCSD(T); see below) and multireference methods (MR-AQCC) to evaluate the possibility of using single-reference methods (including DFT). In all cases, the −CH2OH moiety was not included in the calculations, based on the small effect thereof observed in both experiment and the DFT calculations, in order to reduce the computational requirements.

One of the most common methodological approximations is the second-order Møller–Plesset perturbation theory (MP2)—the main advantage of which is relative computational efficiency—and the <span class="Chemical">SCS-MP2 and SOS-MP2 variants thereof. To investigate their usefulness, optimizations of the SPOCQ reactions of BCN and OCT, respectively, with QUIN were performed using these three perturbation methods with the 6-311+G(d,p) basis set (DIBAC was omitted for reasons of computational resources; all methods contain a ZPE and solvent correction at the standard MP2 level). Again, the SPOCQ cycloadditions under study were found to be highly exothermic. The results for the computed activation energies are summarized in Figure . From this figure one finds a significant underestimation of the activation barrier at MP2 level as compared to the DFT results (Figure ). In the case of BCN, the calculated barrier even yields a value approaching zero. This we attributed to overstabilization of the molecular interaction in the transition states as it was also found in case of stabilization energies of charge-transfer complex between aromatic systems and tetracyanoethylene.[71] This suggests that MP2 itself is indeed not a useful method for these Diels–Alder type of reactions due to the overstabilization of the transition-state complex.

Figure 6

Comparison of activation enthalpies for the SPOCQ reaction of BCN and OCT with 1,2-quinone as calculated by full geometry optimizations at the different MP2 approaches and single-point DLPNO–CCSD(T)/def2-TZVP//SOS-MP2/TZVP+sp and MR-AQCC/6-311G(2d)//SOS-MP2/TZVP+sp calculations. All values are given in kcal/mol.

Comparison of activation enthalpies for the SPOCQ reaction of <span class="Chemical">BCN and OCT with 1,2-quinone as calculated by full geometry optimizations at the different MP2 approaches and single-point DLPNO–CCSD(T)/def2-TZVP//SOS-MP2/TZVP+sp and MR-AQCC/6-311G(2d)//SOS-MP2/TZVP+sp calculations. All values are given in kcal/mol.

The SOS and <span class="Chemical">SCS modifications of the MP2 method compensate for these artifacts of the original MP2 method, with SCS-MP2 being slightly better than SOS-MP2. The resulting activation energies are, however, too high in comparison with experiment. Increasing the basis set from 6 to 311G(d,p) to TZVP+sp for the cases of BCN and OCT reduces the activation barrier by 1.1 and 1.7 kcal/mol, respectively, but they are still a bit too high. Next, SOS-MP2/TZVP+sp and SCS-MP2/TZVP+sp full optimizations were carried out and compared to SOS-MP2 and SCS-MP2 single-point calculations on the fully optimized structures from regular MP2/TZVP+sp calculations. As no solvent corrections were available at these levels, for reasons of consistency, the M06-2X values were taken (also for the DLPNO-CCSD(T) and MR-AQCC methods discussed in the following section). In particular, SCS-MP2 calculations perform well, with calculated barriers for BCN and OCT of 5.5 and 9.3 kcal/mol (exp: 4.5 and 7.3 kcal/mol). Only small differences (<1 kcal/mol; optimization reduces in all cases the barrier slightly) were found in the activation energies for OCT when comparing SCS-MP2//MP2 and SOS-MP2//MP2 single-point calculations with full geometry optimization with the SCS-MP2 and SOS-MP2 methods (0.4 kcal/mol in both cases). In the case of BCN, the difference energy between these two approaches was found to be just slightly higher, 0.7 and 0.8 kcal/mol, for the SCS-MP2 and SOS-MP2, respectively. Therefore, we conclude that full optimization at the SCS-MP2 and SOS-MP2 levels does improve the agreement with experiment, but the effect is rather small when optimization at the SCS-MP2 or SOS-MP2 levels of theory is available (as, e.g., in Orca[72] or Turbomole[73]); this is then to be preferred. However, if that is not the case (e.g., as in Gaussian 16[74]) then single-point SCS-MP2 or SOS-MP2 calculations on MP2-optimized geometries would be useful alternatives. We also note that single-point calculations with the larger aug-cc-pVTZ and aug-cc-pVQZ basis sets indicated that basis set effects beyond the TZVP+sp basis do not seem to be significant. Finally, also at, e.g., the SCS-MP2 level some slight dischronicity of the TS is observed (e.g., from BCN the C–C bond lengths of the bonds to be formed are 2.23 and 2.32 Å, which is analogous to what was observed with e.g. DFT methods; vide supra), thereby confirming the method-independence of the general conclusions.

Coupled-Cluster and Multireference Methods

The second-order perturbation character of the MP2 methods is certainly not sufficient to draw final conclusions about energy barrier heights. Therefore, two substantially more evolved methods have been used as well: (i) CCSD(T) and (ii) MR-AQCC. CCSD(T) represents an excellent standard in cases where the wave function is dominated by one (closed-shell) determinant. The MR-AQCC method has been chosen to test the biradicaloid character of the structures involved in the <span class="Chemical">SPOCQ reaction, in particular, at the transition state. This high-end method was chosen for a number of reasons: First, it allows the inclusion of quasi-degenerate configurations in the reference wave function and provides approximate size-extensive results as it takes account of single and double excitations explicitly when dealing with the dynamic electron correlation.[47] Second, whereas this method was previously only feasible for smaller systems, recently the scope of this method has been significantly extended[75] via the introduction of localized orbitals in combination with the weak pairs (WP) approximation[76] with additional freezing selected orbitals (e.g., the C–H bonds in this case) as described in the Experimental Section and the Supporting Information. This strongly reduces the computational requirements to manageable proportions on reasonably large computer clusters for standard organic chemistry reactions.

Because of the still considerable computational effort (e.g., precluding the calculation of the vibrational frequencies) for both types of computational methods, only the practically most interesting case, the reaction of BCN with <span class="Chemical">1,2-benzoquinone, was considered. The optimized geometries were taken from the SOS-MP2 geometry optimization calculations. In the case of DLPNO–CCSD(T), the def2-TZVP basis set[77] was selected, which yielded a gas-phase enthalpy barrier of 10.5 kcal/mol. The importance of the triples correction becomes obvious from the 15.2 kcal/mol gas-phase enthalpy barrier found at the lower DLPNO–CCSD level. Next, the DLPNO–CCSD(T) gas-phase value was corrected for ZPE and thermal contributions to the activation barrier ΔH⧧ (−0.55 kcal/mol) and solvent effects (−1.8 kcal/mol), both taken from the M06-2X calculations. This thus leads for DLPNO–CCSD(T) to an estimated enthalpy of activation ΔH⧧ in 1,2-dichloroethane of 8.2 kcal/mol.

In the MR-AQCC calculations, the range of computed gas-phase barriers lies between 7.6 kcal/mol (obtained with the 6-31G* basis) and 10.7 kcal/mol (computed with the larger 6-311G(2d)-red basis; the label “red” in this basis indicates that for the hydrogen atoms and the <span class="Chemical">carbon atoms not directly involved in the reaction only the 6-31G* basis was used; Table S13). Our best estimate amounts to 10.5 kcal/mol, which is exactly the same as the DLPNO–CCSD(T) value. To this value, the same ZPE, thermal and solvent corrections were applied as used for the DLPNO–CCSD(T) calculation to obtain the estimated activation enthalpy of 8.2 kcal/mol.

While, without any doubt, the CCSD(T) and MR-AQCC calculations are of the highest levels possible, there is still a 3–4 kcal/mol difference with experiment. We can attribute this to a few factors. First, the geometries at these levels are optimized at the SOS-<span class="Gene">MP2 level, not at the CCSD or MR-AQCC level, so these geometries do not really correspond to the stationary points on the respective potential energies. This likely yields a difference on the order of 1–2 kcal/mol, as was, e.g., observed on the SOS-MP2 and SCS-MP2 single-point and full optimization computations discussed above. In addition, the solvent effect via PCM methods may not be perfect, and we removed the −CH2OH moiety in BCN, the combined effects of which likely also contribute to an uncertainty on the order of a kcal/mol. Finally, the experimental uncertainty is ∼1 kcal/mol. Without aiming to bring these data into accordance with each other, it is clear that the single-point CCSD(T) and MR-AQCC function are as good as one might expect, whereas we regard a better than 2 kcal/mol agreement of experiment and theory (as, e.g., seen for B97D) as also somewhat fortuitous.

There is one more important conclusion to be drawn from these calculations. Because of computational efficiency, DFT has become the method of first choice to study Diels–Alder and related pericyclic reactions. This approach hinges on the assumption that a single-reference approach is correct. That observation that the highly correlated, but single-reference method CCSD(T) gives precisely the same barrier as the equally highly correlated, but multireference-method MR-AQCC indicates that the single-reference methods are good enough for the study of these reactions as long as the TS is more or less symmetric (near-equal lengths of the two C–C bonds that are being formed). Moreover, the largest deviations in the occupation numbers of natural orbitals from the limiting closed-shell case of two and zero, respectively, amount to about 0.12 e in case of the TS structure and show very similar amounts of around 0.1–0.13 e for the reactants. Chemically, this means that there is little biradicaloid character to be found in the TS for both <span class="Chemical">BCN and OCT; computationally, this means that also other single-reference approaches, like DFT or SCS-MP2, are viable. For the larger DIBAC system, with its asymmetric TS (clearly unequal lengths of the two C–C bonds that are being formed), this biradical character might be higher, and this system would thus be an interesting target for future MR-AQCC calculations.

Conclusions

From temperature-dependent kinetics and high-end quantum chemical calculations we determined the activation enthalpies ΔH⧧ for the reaction of several cyclooctyne derivatives in the strain-promoted oxidation-controlled cyclo<span class="Chemical">alkyne–quinone reaction with 1,2-benzoquinones in 1,2-dichloroethane and methanol. Experimentally, we found that ΔH⧧ = 7.3 ± 0.4 kcal/mol for OCT and 4.5 ± 0.3 kcal/mol for endo-BCN, and 12.1 ± 0.5 kcal/mol for the more bulky DIBAC derivative. Quantum chemically, this was studied by a wide range of DFT and ab initio methods as well. Here, we found that DFT methods with explicit inclusion of dispersion (especially B97D) mimic this well, whereas DFT methods that do not treat dispersion explicitly (as, e.g., the Minnesota family of functionals) yield barriers that are typically 4 kcal/mol or more too high. In contrast, MP2 calculations yield barriers that are much too small, an effect that is corrected quite nicely by especially SCS-MP2 (and somewhat less by SOS-MP2). Using CCSD(T) and the multireference method AQCC we were able to show that the transition state of this Diels–Alder-type reaction does not display significant biradicaloid character for BCN (and likely also OCT) and as such can be described well by DFT methods. This holds promise for especially DFT and SCS-MP2 approaches, as those start off from the single-reference, closed-shell approximation that is hereby confirmed to be allowed. Finally, the appreciable reactivity differences between the various cycloalkynes can be seen as the effect of three factors: for the fastest cyclooctyne under study (BCN), the TS was most synchronized in terms of formation of the two C–C bonds that are formed, the TS geometry required the smallest distortion energy to be formed, and the interaction between the components was most optimal as observed from both the interaction energies and, e.g., the orbital overlap in the TS. For slower SPOCQ reactions, all these factors were less favorable. These data therefore provide unprecedented mechanistic insight in this highly useful reaction that has displayed unique characteristics in terms of both rate and yield in both bioconjugation and surface reactions. Furthermore, they provide a route for the study of other, related reactions that convert π bonds into σ bonds.

Experimental Section

Kinetic Measurements

UV/vis measurements were performed on a Varian Cary-50 spectrophotometer using a 1 cm quartz cuvette. A stock solution of OCT was prepared in <span class="Chemical">DCE or MeOH to a final concentration of 400 μM. 4-tert-Butyl-1,2-benzoquinone was also dissolved in DCE or MeOH to a final concentration of 200 μM. For a total volume of 2 mL, 1 mL of the OCT solution (final concentration 200 μM) was added to the quartz cuvette containing 1 mL of the 4-tert-butyl-1,2-benzoquinone solution (final concentration 100 μM). The reaction was followed at a wavelength of 384 nm with a scanning interval of 0.01 min and was performed in triplicate.

Determination of Activation Parameters

The activation parameter enthalpy of activation (ΔH⧧) was obtained using the Eyring-Polanyi equation by a least-squares fit of plots of ln (k/T) versus 1/T. Values for ΔH⧧ and ΔS⧧ were calculated from the slope and intercept of the aforementioned plotwhere k = measured rate constant (in s–1); T = absolute temperature (in K); K = transmission coefficient (which was set to be 1); kB = Boltzmann’s constant; h = Planck’s constant; ΔS⧧ = entropy of activation; R = gas constant, and ΔH⧧ = enthalpy of activation.

Computational Details

The calculations were performed with the Gaussian 16,[74] Turbomole,[73] Orca,[72] and COLUMBUS[78,79] suite of programs. The stationary points were characterized by means of harmonic vibrational frequency analysis at DFT and the different MP2 levels. For all of the transition structures, the normal mode related to the imaginary frequency corresponds to the nuclear motion along the reaction coordinates under study. Additionally, in some cases, we carried out intrinsic reaction coordinate calculations (IRC) to verify that the transition structures connect with reactants and products. DFT calculations were performed using BLYP, B97D, B3LYP, M06-2X, <span class="Mutation">M06L, M11, M11L, and MN15 functionals, all with the 6-311+G(d,p) basis set. The MP2, SOS-MP2, and SCS-MP2 calculations have been performed using the TZVP+sp basis, in which diffuse functions (s,p) are added to the triple-z valence polarized (TZVP).[80] The DLPNO-CCSD and DLPNO-CCSD(T) calculations have been performed with the larger def2-TZVP basis set.[81]

A complete active space (CAS) of eight electrons and eight orbitals (CAS(8,8)) was used in the CAS self-consistent field (CASSCF) calculation and as reference space in the MR-AQCC calculation. The resulting shape of these eight active orbitals is displayed in Figure S1 for the transition state. The plots show that these orbitals describe the six orbitals considered in the standard 2s+4s cycloaddition of ethylene + <span class="Chemical">butadiene plus the bonding and antibonding orbitals of the orthogonal acetylenic π bond in BCN. Full MR-AQCC(8,8) calculations correlating all valence orbitals are very time-consuming. Therefore, the calculations were performed in terms of localized orbitals which were obtained for the reference doubly occupied orbitals. The localization allowed freezing of types of orbitals not participating in the cycloaddition reaction. Furthermore, the weak-pair approximation[75,76,82] was used as implemented into the COLUMBUS program system, which allowed the omission of interaction of spatially distant electron pairs. More details of this approach are described in the Computational Details in the Supporting Information. In most of the calculations, the 6-31G* basis was used. However, for the best estimate the 6-311G(2d) basis set was placed on the atoms located in the region of the reaction. All CASSCF, MR-CISD, MR-CISD+Q, and MR-AQCC calculations, including those with the weak pairs approximation, were performed with the COLUMBUS program system.