Bouncing off walls - widths of exit channels from shallow minima can dominate selectivity control.

A selectivity model based on the widths of pathways to competing products, rather than barrier heights, is formulated for the butadiene + allyl cation reaction. This model was arrived at via analysis of stationary points, intrinsic reaction coordinates, potential energy surface shapes and direct dynamics trajectories, all determined using quantum chemical methods. This journal is © The Royal Society of Chemistry.

Introduction

Reactions between butadienes and allyl cations have been studied extensively from both mechanistic and synthetic perspectives. Reactions of these two types of components can lead to 4-, 5-, 6- or 7-membered rings via concerted or stepwise (formal) (2 + 2)/[π2 + π2],[1,2] (2 + 3)/[π2 + π2],[3,4] (4 + 2)/[π4 + π2][5-9] (the so-called “ionic Diels–Alder reaction”) or (4 + 3)/[π4 + π2][3,4,10-12] cycloaddition reactions, respectively (Scheme 1). Which products are observed can be influenced by which substituents are attached to each component and the environment in which the reaction is run. Here we reinvestigate the parent reaction of butadiene + allyl cation (in the gas phase and DMSO), adding results from density functional theory (DFT) calculations to Pascual-Teresa and Houk's previously reported HF and MP2 results,[13] and including an expanded view of the butadiene + allyl cation potential energy surface (PES). In addition, we have examined variational transition states (VTSs)[14] and subjected this reaction to ab initio molecular dynamics (AIMD) simulations, of both downhill and uphill varieties,[15] to assess whether or not non-statistical dynamic effects play a role in determining product selectivity.[16-19] Our results have led to the conclusion that the widths of pathways to products can sometimes be more important than barrier heights—an entropic effect that has often been overlooked but has potential utility in constructing useful models of selectivity and in designing new selective reactions.

Scheme 1

Possible modes of cycloaddition for butadiene + allyl cation.

Computational methods

The Gaussian 09 software suite[20] was used for all DFT calculations. The geometries and frequencies of minima and transition-state structures (TSSs) were optimized at the B3LYP-D3/6-311+G(d,p),[21-25] M06-2X/6-311+G(d,p),[26-28] M06-2X/6-311G(d), and ωB97X-D/6-311G(d)[29] levels of theory. TSSs were defined as such by frequency calculations that predicted exactly one imaginary frequency,[30] while minima showed no imaginary frequencies. Molecular structures were visualized using CYLview.[31] Intrinsic reaction coordinate (IRC)[32-34] calculations were performed to confirm that minimum energy pathways (MEPs) connected TSSs to the reactants and the products indicated. To explore the effect of solvation, calculations were also performed with M06-2X/6-311G(d) and the PCM solvation model for DMSO.[35] In addition, single-point energies were computed at the CCSD(T)-F12a/AVDZ[36] level using the MOLPRO 2015 package.[37] As shown in Fig. 1, the relative energies of stationary points are similar among different levels. VTSs and associated reaction paths were determined using Gaussrate[38]/Polyrate[39] with application of the Reorientation Of the Dividing Surface (RODS) algorithm at the M06-2X/6-311G(d) level of theory.[14,40]

Fig. 1

Relative free energies (single point electronic energies for CCSD(T) results) of stationary points on the allyl cation + butadiene PES (in kcal mol−1 relative to Rc) at 298.15 K. The values shown are from top to bottom: B3LYP-D3/6-311+G (d,p) (black), M06-2X/6-311+G(d,p) (red), M06-2X/6-311G(d) (blue), ωB97X-D/6-311G(d) (green), PCM(DMSO)-M06-2X/6-311G(d) (purple), CCSD(T)-F12a/AVDZ//M06-2X/6-311G(d) (yellow). Distances shown (Å) are from structures optimized with M06-2X/6-311G(d). Structure labels are defined in the text below.

Both quasi-classical and classical AIMD simulations were performed using M06-2X/6-311G(d). For the former, zero-point energies are included in initial sampling while for the latter, they are not. Quasi-classical AIMD simulations also were performed using ωB97X-D/6-311G(d). All trajectories were propagated using the Progdyn[41] script. To explore the effect of continuum solvent on dynamics, the PCM(DMSO)-M06-2X/6-311G(d) level of theory was employed. Uphill (initiated in the Rc region) and downhill (initiated in transition state regions) dynamics simulations were both carried out. For downhill dynamics, starting geometries were generated from a Boltzmann distribution of vibrations at 298 K and trajectories were propagated in both product and reactant directions with the Verlet algorithm, using Gaussian 09 to calculate forces at each time step (the time step was set as 1 fs).[42] For uphill dynamics simulations, starting geometries were generated from a Boltzmann sampling of vibrations at 298 K then propagated towards the product direction until one or the other of the products is formed. Starting points for uphill dynamics are overlaid in Fig. 2. Note that our analysis here is for a “relaxed” Rc; i.e., momentum for its formation is not included.[18,19] The following geometric criteria (atom numbers shown in Fig. 1) were employed to halt trajectories: when the C1–C4 and C5–C7 distances both dropped below 1.70 Å, a trajectory was labeled as the (4 + 3) adduct (7-membered ring product, Prod-1); when the C1–C4 and C2–C7 distances both dropped below 1.60 Å, a trajectory was labeled as the (2 + 3) adduct (5-membered ring product, Prod-2).

Fig. 2

Overlays of all starting geometries consistent with a Boltzmann distribution at 298 K in the region of Rc: (a) M06-2X/6-311G(d), (b) ωB97X-D/6-311G(d), (c) PCM(DMSO)-M06-2X/6-311G(d).

Results and discussion

Reaction coordinates

Energy profiles from different levels of theory for the allyl cation + butadiene reaction are shown in Fig. 1. No TSSs for concerted cycloaddition or formation of Rc, the product of addition of one end of the butadiene π-system to one end of the allyl π-system, were found; it appears that adduct formation is barrierless, at least in the absence of explicit solvent (appropriate modeling of which is beyond our current capabilities). This conclusion is consistent with previous studies on related reactions.[11] Cation Rc can then be converted to 7- or 5-membered ring-containing products via concerted processes involving TS-a and TS-b, respectively. Note that the CC group involved in bond formation (that which originated in the allyl fragment) rotates in opposite directions to form the 7- and 5-membered ring-containing products (i.e., for an overall pathway connecting Prod-1 to Prod-2, the migrating CH2 group interacts antarafacially with the allyl group over which it migrates). While tempting to consider this observation a result of “vestiges of orbital symmetry,”[43] that cannot be done with confidence without more detailed examination of the PES preceding Rc. At all DFT levels of theory (but not CCSD(T), for which the data correspond to single-point electronic energies), TS-b is predicted to be higher in energy than TS-a and Prod-2 was predicted to be higher in energy than Prod-1. Namely, in terms of either kinetic or thermodynamic control, formation of the 7-membered ring product, Prod-1, is favored. However, Rc is a shallow minimum, which can increase the lifetime and play a significant role in the dynamics, particularly at low temperatures.

Variational transition states

To address issues with standard computations of entropies, reaction pathways derived using variational transition state theory (VTST) calculations (using M06-2X/6-311G(d)) were computed; these are shown in Fig. 3. For both pathways, VTSs differ from the PES TSSs (compare the geometries of the VTSs in Fig. 4 with those of the TSSs in Fig. 1). The C5–C7 distance in VTS-a is 0.02 Å longer than that in TS-a, while the C2–C7 distance in VTS-a is 0.01 Å shorter than that in TS-a. The C5–C7 and C2–C7 distances in VTS-b are both shorter than those in TS-b, by 0.05 Å and 0.09 Å, respectively. Both VTSs are higher in free energy than the corresponding TSSs. As shown in Fig. 3, this effect is not simply an entropy effect, but also an effect of zero point energy. Note, however, that the difference in energy between the peaks of the free energy and potential energy curves is small—tenths of a kcal mol−1—consistent with flat surfaces, in terms of both potential energy and free energy, near the transition states. No evidence for “entropic intermediates” is found here.[44,45]

Fig. 3

Reaction coordinate plots (M06-2X/6-311G(d)) for pathways a (left) and b (right). The locations of VTS-a and VTS-b are highlighted.

Fig. 4

Geometries (M06-2X/6-311G(d)) of variational transition states. Selected distances are shown in Å.

Direct dynamics trajectories

Given the flatness of the PES near Rc and the transition states for cyclization, AIMD trajectory simulations were employed to evaluate (1) whether either transition state can lead to both products and (2) the factors influencing product selectivity.

Downhill trajectories initiated from TS-a did in fact lead to both products, but just barely (Table 1), while trajectories initiated from TS-b did not. In both cases, however, a large amount of recrossing (here including trajectories connecting the reactant to the reactant, one product to the same product, or one product to the other product) was observed, consistent with a flat region of the PES.

Results of downhill dynamics simulations initiated from TS-a and TS-b (M06-2X/6-311G(d))

	Total trajectories	Trajectories that form Prod-1	Trajectories that form Prod-2	Trajectories that recross
TS-a	702	403	3	296
TS-b	699	0	451	248

Downhill trajectories were also initiated from VTS-a and VTS-b. Only Prod-1 was formed from VTS-a and only Prod-2 was formed from VTS-b (Table 2). While recrossing trajectories were much reduced for VTS-b, compared to TS-b, a high percentage of recrossing was still observed for VTS-a.

Results of downhill dynamics simulations initiated from VTSs (M06-2X/6-311G(d))

	Total trajectories	Trajectories that form product	Trajectories that recross
VTS-a	139	64	75
VTS-b	103	93	10

Results for uphill AIMD trajectories initiated at Rc are shown in Table 3. At all levels of theory, formation of the 7-membered ring-containing product, Prod-1, is much preferred. However, the predicted Prod-1/Prod-2 ratio is lower than expected based on differences in TSS free energies (Fig. 1). For example, at ωB97-XD/6-311G(d), the ratio from dynamics simulation is 90 : 10, while that from free energy differences for TSSs is 97 : 3. As described above, VTS free energies differ only slightly from TSS free energies. With PCM(DMSO)-M06-2X/6-311G(d), free energy differences predict essentially no Prod-2, but our dynamics simulations predict that 16.5% of trajectories lead to this product. What is the origin of this discrepancy? Is it technical or chemical in nature?

Results of uphill dynamics simulations initiated from Rc

Method	Total trajectories	Trajectories that form Prod-1	Trajectories that form Prod-2
M06-2X/6-311G(d)	206	171 (83%)	35 (17%)
ωB97X-D/6-311G(d)	154	138 (90%)	16 (10%)
PCM(DMSO)-M06-2X/6-311G(d)	91	76 (84%)	15 (16%)

Since our trajectories often involved thousands of fs to form products (see ESI†), we were concerned about ZPE leakage (generally a worry past ∼0.5 ps).[46-49] To address this issue, we also carried out classical AIMD simulations. Results (using M06-2X/6-311G(d)) for uphill AIMD trajectories initiated at Rc and for downhill AIMD trajectories initiated at VTS-a are shown in Tables 4 and 5, respectively. For uphill simulations, the predicted Prod-1/Prod-2 ratio does not change significantly between classical and quasi-classical simulations (compare Tables 3 and 4). For downhill simulations, a significantly smaller proportion of trajectories recrossed for classical vs. quasi-classical dynamics (compare Tables 2 and 5), but a large proportion of recrossing was still observed.

Classical dynamics results of uphill dynamics simulations initiated from Rc (M06-2X/6-311G(d))

	Total trajectories	Trajectories that form Prod-1	Trajectories that form Prod-2
Rc	90	75 (83%)	15 (17%)

Classical dynamics results of downhill dynamics simulations initiated from VTS-a (M06-2X/6-311G(d))

	Total trajectories	Trajectories that form product	Trajectories that recross
VTS-a	130	91	39

Potential energy surfaces

To explore how shape of the PES influences the selectivity of the reaction, PESs for gas phase (M06-2X/6-311G(d)) and DMSO (PCM-M06-2X/6-311G(d)) reactions were built (Fig. 5) using relaxed scans along C5–C7 and C2–C7 distances (see Fig. 1 for atom labels) as x- and y-axes and energy displayed as height and color (top surfaces) or color alone (bottom surfaces, projections of top surfaces into the x–y plane). While these are reduced-dimension PESs, the geometric parameters used correspond to those of the forming bonds.[50,51] Small C5–C7 distances correspond to the Prod-1 region, while small C2–C7 distances correspond to the Prod-2 region. The locations of the two TSSs are marked, as are the IRC paths. One observation is immediately apparent from these surfaces: the pathway toward Prod-1 is wider than the pathway toward Prod-2.

Fig. 5

Potential energy surfaces: (a) M06-2X/6-311G(d), (b) PCM(DMSO)-M06-2X/6-311G(d). IRCs for TS-a are shown as red lines and IRCs for TS-b are shown as green lines. TS-a and TS-b are indicated as stars and circles, respectively.

Representative gas phase trajectories are plotted in Fig. 6. As expected, long-time trajectories sample more of the Rc region on the PES before forming products. Importantly, some trajectories heading toward product “turn back” before crossing the dividing surface that separates reactant from products, i.e., they miss the exit channel and instead “bounce off the wall” flanking it. This is distinct from recrossing, where the dividing surface is passed before a trajectory turns back (at least for the dimensions we have plotted here).[16-19] This is also distinct from bouncing off the side walls of a valley en route to product(s), i.e., where deviations from an IRC still lead to product(s), not back to reactant(s).[16] The proportion of Prod-1 forming trajectories that arose from “turning back” (62/171 or 36%) is larger than the proportion of Prod-2 forming trajectories that arose from “turning back” (7/35 or 20%), i.e., more Prod-1 is formed from trajectories that initially miss the exit channel toward which they started, presumably because the exit channel toward Prod-1 is wider. The product ratio predicted on the basis of “turn back” trajectories alone is 90 : 10 (62 vs. 7 trajectories). The ratio predicted on the basis of direct trajectories (here there are 137, making for a total of 206 as shown in Table 3) is 80 : 20, and the overall ratio accounting for all product-forming trajectories is 83 : 17 (Table 3). Thus, the ratio predicted from direct uphill trajectories does not reflect the predicted free energies of transition states, be they TSSs or VTSs, and the “turning back” phenomenon increases selectivity for the product formed via the wider pathway. A similar effect was recently discussed for a different carbocation reaction, but from the perspective of downhill trajectories and post-transition state bifurcations.[52]

Fig. 6

(1) Projections of trajectories onto the M06-2X/6-311G(d) PES, based on trajectory length (time). Left to right: 0–500 fs, 500–1000 fs, 1000–1500 fs, 1500–2000 fs. The first row corresponds to Prod-1-forming trajectories. The second row corresponds to Prod-2-forming trajectories. (2) Projections of representative “turn back” trajectories onto the M06-2X/6-311G(d) PES. The third row shows representative Prod-1-forming “turn back” trajectories. The bottom row shows representative Prod-2-forming “turn back” trajectories.

A close-up view of the reactant region of the PES is shown in Fig. 7, which also shows the locations of the VTSs. TS-a and VTS-a occupy similar positions in this plot, consistent with similar amounts of recrossing in downhill AIMD. VTS-b is noticeably later (closer to product in terms of the geometric parameters plotted here) than is TS-b, consistent with reduced recrossing. Being later, it is further into the least wide region of the path to products, consistent with it being even more difficult for uphill trajectories to reach Prod-2 than if they only had to reach TS-b. This region of decreased width indicates that varying the C5–C7 distance happens to be particularly difficult near the transition state for formation of the C2–C7 bond. This phenomenon was unexpected, but is perhaps related to the fact that the forming σ-bond is much shorter in TS-b and VTS-b than in TS-a and VTS-a and conformational restrictions due to hyperconjugation with the forming σ-bond are expected, therefore, to be more severe.

Fig. 7

Close-up view of the reactant region of the surface from Fig. 5a.

This pathway width effect is clearly connected to entropy, i.e., a wider path is related to more flexibility and more possible ways to cross through a dividing surface.[40,45,53-55] This effect appears not to be captured fully in VTS free energies, however, at least for the system examined here, which involves a flat energy surface, although we cannot definitively rule out the possibility that the apparent discrepancy is due to the particulars of running each type of calculation. Nonetheless, interpreting results of this type in terms of pathway widths allows one to better tie specific structural features to entropy changes along reaction coordinates.

Conclusions

We have explored the mechanism of the reaction between butadiene and allyl cation using quantum chemical computations on PESs and direct AIMD simulations (quasi-classical and classical). Consistent with results from previous calculations at other levels of theory,[13] concerted cycloaddition pathways are not found. Instead, stepwise pathways leading to products of formal (4 + 3) and (2 + 3) cycloaddition are present, but these involve a shallow intermediate with low barrier exit channels. As for other reactions with shallow intermediates—twixtyls, calderas, mesas, para-intermediates[19,56-58]—non-traditional effects appear to play a key role in product formation. We propose that the product selectivity in the system examined here—and likely others[52] —is determined by dynamic effects related to the width of pathways to products. In particular, there is an entropy effect—transition states with wider energy wells orthogonal to the reaction coordinate have more possible paths through them—which is modulated, enhanced in this case, by an unequal probability of trajectories turning back before reaching transition states. We are currently expanding on this work to see if models based on pathway width can rationalize product selectivity for cases with existing experimental evidence for entropy-controlled selectivity where the origins of that entropy control are not known.

Conflicts of interest

The authors declare no competing financial interest.