Sanha
Lee
and
Jonathan M.
Goodman
*
Yusuf Hamied Department of Chemistry, University of Cambridge, Lensfield Road, Cambridge, CB2 1EW, UK. E-mail: jmg11@cam.ac.uk
First published on 6th April 2021
In recent years, a growing number of organic reactions in the literature have shown selectivity controlled by reaction dynamics rather than by transition state theory. Such reactions are difficult to analyse because the transition state theory approach often does not capture the subtlety of the energy landscapes the compounds traverse and, therefore, cannot accurately predict the selectivity. We present an algorithm that can predict the major product and selectivity for a wide range of potential energy surfaces where the product distribution is influenced by reaction dynamics. The method requires as input calculation of the transition states, the intermediate (if present) and the product geometries. The algorithm is quick and simple to run and, except for two reactions with long alkyl chains, calculates selectivity more accurately than transition state theory alone.
For example, Singleton et al. discovered the experimental product ratios in ozonolysis of alkyl vinyl ethers do not fit the predictions based on the transition state theory, Scheme 1.6 Mechanistically, the vinyl ether 4e and O3 react together to form the intermediate 4-INT. The intermediate can either fragment to products 4a and 4bvia TS2A or 4c and 4dvia TS2B. Singleton et al. discovered the comparison of barrier heights of TS2A and TS2B would predict the A/B ratio to be more than three orders of magnitude higher than the experimentally observed product distributions. The unexpected selectivity arises because the intermediate formed from TS1 has enough excess energy to favour the cleavage leading to products B over A. Therefore, TST does not give an accurate product ratio prediction.
The importance of reaction dynamics has been highlighted in many other organic reactions. In bifurcating reactions, a single TS is shared by two or more reaction pathways leading to different products. The product distribution is then governed by the shape of the potential energy surface (PES) and the resulting dynamic effects.7,8 Reactions with very shallow intermediates on the PES can also show selectivity dependent on nonstatistical dynamics. Examples of well-known organic reactions showing such behaviour include hydroboration,9,10 C2–C6 cyclisation,11 Garratt-Braverman cyclisations,12 and ketene cycloadditions.13 Shallow intermediates on a single reaction profile can also lead to multiple reaction paths due to the ‘concerted components’ and ‘stepwise components’. Acetone radical cation fragmentation,14 1,2,6-heptatriene rearrangement,15 cope rearrangement of 1,5-hexadiene16 and triplet di-π-methane rearrangement17 are examples showing this behaviour. Reactions involving a flat region on the PES called a ‘caldera’ also show dynamically controlled selectivities.18,19 Finally, Merrer et al. found the reaction between dichlorocarbene and 1,2-disubstituted cyclopropenes preferred the thermodynamically unfavoured reaction pathway due to the conservation of momentum.20
We have previously developed an algorithm ‘ValleyRidge.py’ which can quickly predict selectivity for organic reaction proceeding on bifurcating potential energy surfaces.21 The algorithm achieves this by importing the transition states and the product geometries and reducing the dimensionality of the PES by identifying the key bond difference between the products. The major product is determined by analysing the direction the imaginary eigenvector of the first transition state (TS1) points on the ridge. Finally, the product ratios are calculated using the harmonic approximation on the first transition state real normal modes.
ValleyRidge.py was developed for bifurcating organic reactions with a valley-ridge inflection point (VRI) on the PES, Fig. 1a. However, the study by Singleton et al.13 motivated us to investigate how this algorithm could be extended beyond these applications to PES without a VRI. When the cycloaddition of dichloroketene is modelled using the MPW1K functional, the reaction path does bifurcate as shown in Fig. 1a. When the same reaction is modelled using the B3LYP functional, a shallow intermediate is now present between the two products as illustrated. However, the reaction trajectories can easily bypass this shallow intermediate to form the products directly from TS1. Our new algorithm is able to predict the selectivity in such circumstances and to a range of other PES types with the selectivity controlled by nonstatistical dynamics.
Fig. 1 (a) The qualitative PES for cycloaddition of dichloroketene modelled by Singleton et al.13 for the MPW1K surface and the B3LYP surface. P1 and P2 are the two products, INT is the intermediate and ā is the imaginary eigenvector of TS1. (b) Fourty. |
We present this new algorithm, Valley Ridge Augmented Implementation Selectivity (VRAI-selectivity), which can model reactions with a shallow intermediate and diverging reaction pathways, as well as TST-controlled reactions. We summarise the different PES VRAI-selectivity can model in Fig. 1b.
Importantly, the algorithm automatically picks out the bond differences and, thus, the dimensionality reduction does not require any user input.21 When different users model the same reaction, this ensures that the same set of bond differences are always used.
The simplified 2D PES diagram for Fig. 1a, B3LYP PES with a shallow intermediate, is shown in Fig. 3. If we let vector ḡ be the separation of the intermediate from TS1, the angle ϕ becomes the angle the imaginary eigenvector of TS1 (ā) makes with vector ḡ. We can define vectors and to be the displacement vectors from TS2 to product 1 and product 2, respectively. By defining the product vectors, θ1 becomes the angle vector −ḡ makes with vector (product vector on the side of ḡ that vector ā points directly) and θ2 becomes the angle −ḡ makes with product vector on the opposite side of vector . VRAI-selectivity rejects the bond difference pair if ϕ is less than 0.1 degrees and selects the next highest pair in the rank. Qualitatively, the algorithm examines which side of the intermediate the imaginary eigenvector of TS1 points to decide the major product. The harmonic oscillator approximation is then used on the real eigenvectors of TS1 to estimate the width of the trajectory stream. The selectivity of the reaction is then estimated by considering the point of the closest approach of the intermediate and how much of the trajectory stream favour each product. We have found a strong linear correlation between the prediction error and [ϕ + (π − θ2)] (see ESI, section 4†). We use this relationship to correct the predicted selectivity. We will discuss all the reasons for these changes from the earlier algorithm in the following sections.
Fig. 4 Energy profile from Singleton et al. 200613 study at B3LYP/6-31G(d) level. |
The free energy profile for the diphenylketene reactant shows different characteristics to the dichloroketene case (Fig. 4b). The free energy drop from TS1 is smaller at 1.5 kJ mol−1. Furthermore, the barrier height for [4 + 2] and [2 + 2] additions are significantly higher than the 1.5 kJ mol−1 free energy gained from intermediate formation at 3.8 and 9.4 kJ mol−1 respectively. In such circumstances, we would expect the intermediate will be sufficiently relaxed before proceeding over TS2 and the dynamic effects will be negligible. VRAI-selectivity therefore, uses the free energy difference between the TS2A and TS2B to calculate the selectivity using the TST. The [4 + 2] addition product is predicted to be the major product with 90.4% selectivity. The MD prediction at B3LYP/6-31G(d) level is 85.7% and therefore VRAI-selectivity is able to replicate the MD ratio.
The results show two important conclusions. For nonstatistical predictions from VRAI-selectivity to be valid
1. The energy drop from TS1 should be sufficiently large that the intermediates gain significant energy.
2. The barriers for the subsequent reactions are small relative to the energy drop from TS1 and, therefore, the trajectories pass over the intermediates easily and the dynamic effects dominate.
For all the reactions we have analysed, if TS2A and TS2B are both lower in free energy than TS1 by more than 6.6 kJ mol−1, the dynamics analysis gives the correct selectivity predictions (see ESI, section 6†). We therefore approximate the behaviour of the molecules to be Arrhenius and the algorithm decides whether to use TST or the dynamics based on the barrier height difference between TS1 and the two TS2s. If either of the TS2 barrier heights is greater than TS1 free energy minus 5 kJ mol−1, the algorithm will use TST analysis to predict the selectivity. Otherwise, VRAI-selectivity will output the nonstatistical selectivity analysis as seen from dichloroketene example.
Considerable dynamic effects are observed in the hydroboration reaction studied by Singleton et al.9 (Scheme 2, reaction (3)). Truhlar et al. have reported a qualitative theory to account for the selectivity.22 The intermediate in this reaction is the reactive complex formed by molecules 3c and 3d coming together. Singleton et al. found the formation of the intermediate is barrierless on the enthalpic surface but does have a barrier on the free energy surface. Locating such point on the free energy surface is difficult. We therefore decided the TS1 geometry should simply be the free energy barrier geometry quoted by Singleton et al. The frequency calculation on this structure shows two imaginary frequencies. Performing the bond length constrained optimisation removes one of the imaginary frequencies as the alkene and the BH3 became staggered from eclipsed, Fig. 5. However, we believe the eclipsed structure is better because the IRC calculation required fewer steps and showed a smoother profile (see ESI, section 3†).
Fig. 5 Eclipsed and staggered TS1 geometries from Singleton et al. 2009.9 |
The computed energy profile showed a significant free energy drop of 19.7 kJ mol−1 from the eclipsed TS1 geometry to the intermediate complex. The subsequent anti-Markovnikov and Markovnikov activation barriers are 2.9 and 13.0 kJ mol−1 respectively. The TST selectivity for anti-Markovnikov product is 98.3%, higher than the MD simulation predicted 90.0%. Expressing the product selectivity as ratio highlights the difference better as TST ratio is 59:1 whereas the MD ratio is 9:1. The VRAI-selectivity prediction for anti-Markovnikov product is 91.4% with the intermediate activation, showing a better agreement with the MD simulation results.
As previously outlined, the MD study of ozonolysis reaction by Singleton et al.6 (reactions (4)–(7), Scheme 2) found that the experimental product ratios do not fit the expectations based on the TST. We investigated the same reaction and the calculated the energy profile for reaction (4) is shown in Fig. 6. The free energy drop from TS1 to INT is very large at 223.0 kJ mol−1. The barriers for fragmentation steps are less than a third of the energy gained by the intermediate and, therefore, dynamic effects are likely to prevail. The energy difference in subsequent fragmentation transition states is fairly large at 21.6 kJ mol−1 and TST would predict the major product ratio to be over 6000:1. This prediction is much higher than the experimental selectivity of 96.3%. The VRAI-selectivity predicts the major product selectivity to be 97.1%, in better agreement with the experimental result.
Fig. 6 Energy profile from Singleton et al. 20116 (reaction (4)). |
We then explored the effect of the chain length on the selectivity prediction. The modification of the R group to butyl, octyl and dimethyl octyl groups does not change the free energy drop from TS1 significantly since the values were 221.7, 221.7 and 221.8 kJ mol−1 respectively. The difference in TS2A and TS2B free energy gap does not change significantly either and the values were 19.0 kJ mol−1, 19.0 kJ mol−1 and 18.9 kJ mol−1 for butyl, octyl and dimethyloctyl groups respectively. The TST selectivity prediction still remains very high above 2000:1.
Modifying the R group from methyl to butyl does not change the experimental selectivity much. The VRAI-selectivity predicts the major product selectivity to be 95.0%, again showing a good agreement with the experimental 97.8%. Modifying the R group to octyl and dimethyl octyl groups increases the experimental selectivity to 98.2% and 98.3%, respectively. The high frequency vibrations in such reactions with long alkyl chains have very small real eigenvector components. We therefore modelled the octyl and the dimethyl octyl group reaction with VRAI-selectivity but ignored the real eigenvector components that have zero components up to 5 decimal places. The VRAI-selectivity results are less successful for these reactions since the product selectivity are predicted to be 88.7% and 89.0% for octyl and dimethyl octyl groups respectively. The long chain molecules have more normal modes, and so TS1 is broader. Following the VRAI-selectivity analysis, this leads to a lower selectivity. The experimental results have higher selectivity, suggesting a switch from the non-statistical process back to a TST controlled process because the relaxation of the intermediate is faster. Capturing the effect of the long alkane chain on the selectivity is currently beyond the capability of the algorithm but is something we are looking to investigate in our future work.
Schmittel et al.12 found significant nonstatistical dynamic effects in [1,5]-H shift of 8c (reaction (8)). We investigated the same reaction and the calculated energy profile is shown in Fig. 7. The intermediate, TS1 and the H-transfer TSs were modelled as triplet states using the UB3LYP open shell method. The intermediate would have plenty of excess energy from the 246.9 kJ mol−1 free energy drop from TS1. The Z-alkene forming H-transfer TS (TS2A) has slightly lower free energy barrier than E-alkene forming H-transfer TS. Therefore, TST would incorrectly predict the major product to be the Z-alkene. Furthermore, the free energy difference between the H-transfer TSs is very small at 0.5 kJ mol−1. The predicted TST major product percentage would be 45.0%, long way off from the experimental 91.1%. The VRAI-selectivity on the other hand correctly predicts the E-alkene major product with 92.2% selectivity. Therefore, the VRAI-selectivity result has better agreement with the experimental outcome.
We modelled the reaction using the VRAI-selectivity code without the intermediate activation. We used the TS connecting the bicyclobutane and butadiene product as the TS2 geometry since no stationary point is present at the branching point. The algorithm predicted the bicyclobutane to be the major product with 86:14 ratio, showing a good agreement with the experimental results.
We modelled this reaction using a combination of TST and VRAI-selectivity method without intermediate activation. The proportion of the reaction trajectory passing over TS1 and TS1d is estimated using the transition state theory. The increase in free energy going from the starting material to TS1 and TS1d is 84.0 and 84.6 kJ mol−1 respectively for reaction (10). Therefore, the use of TST is justifiable. The proportion of the reaction trajectories passing through TS1 reaching directly the products is then estimated using VRAI-selectivity. For reaction (2), 89.5% of the molecules are estimated to go directly from TS1 to 10a and remaining 10.5% will reach the intermediate. Applying VRAI-selectivity to TS1d estimated 84.3% of the trajectories to reach 10b directly whereas the remaining 15.7% will fall into the intermediate well. Finally, we used TST to estimate the proportion of the molecules in the intermediate well that will reach 10a by TS2 and 10b by TS2d. The intermediate should be sufficiently relaxed because the barrier heights of TS2A and TS2B are very similar to the free energy drop from TS1 to the intermediate.
The starting material can reach product 10a by three possible paths: [S → TS1 → P], [S → TS1 → INT → TS2 → P] or [S → TS1d → INT → TS2 → P]. We calculated the percentage of the trajectories reaching 10a by multiplying the probabilities at each stage and summing over the three possible routes. The method then predicts the 10a:10b ratio to be 53:47 which is in good agreement with the experimentally observed ratio of 60:40. The experimental selectivity increases for reactions (11) and (12) as 11a:11b ratio is 63:37 and 12a:12b ratio is 71:29. Our analysis predicted the 11a:11b ratio to be 64:36 and 12a:12b ratio to be 74:26. Therefore, VRAI-selectivity is able to predict the selectivity as well as the experimental trend as the R groups are changed.
In Fig. 1b, we previously summarised the four types of potential energy surfaces for which VRAI-selectivity analysis can correctly predict the selectivity. The first case is when the free energy drop from TS1 is small and both TS2A and TS2B have high barriers from the intermediate. The VRAI-selectivity selects TST to predict the selectivity as explained for the Singleton et al. diphenylketene example. VRAI-selectivity is also able to analyse the selectivity for PES of type (a) if TS1 is not present and the starting material connects the two products via TS2A and TS2B. The algorithm would simply use TST, only requiring TS2A and TS2B data as the input. The second case is when the free energy drop from TS1 is large and the intermediate has sufficient excess energy to go over the TS2 barriers. The selectivity is then controlled by the dynamic effects and TST often cannot predict the selectivity accurately. The VRAI-selectivity approach can accurately predict the MD and the experimental selectivity, Table 1. The third type of potential energy surface is the Merrer et al. example where the minimum energy path does not connect TS1 to the major product. TST approach is not applicable for such PES since no intermediate is present. However, VRAI-selectivity analysis can predict the correct selectivity. The final type of PES is bifurcating reaction profile. TST is again not applicable for such PES. However, the VRAI-selectivity analysis predicts the correct major product and the selectivity using the TS2 geometry instead of the intermediate geometry. VRAI-selectivity can, therefore, analyse selectivity on PES for which TST is applicable as well as other types of PES where TST will not give a correct outcome.
Study | Reaction | Exp/Traj % | TST % | VRAI % |
---|---|---|---|---|
Singleton 200613 | Cl | 83.3 | 52.4 | 82.9 |
Singleton 200613 | Ph | 85.7 | 90.4 | 90.4 |
Singleton 20099 | B | 90.0 | 98.3 | 91.4 |
Singleton 20116 | Me | 96.4 | 100.0 | 97.1 |
Singleton 20116 | But | 97.8 | 100.0 | 95.0 |
Singleton 20116 | Oct | 98.2 | 100.0 | 88.8 |
Singleton 20116 | DimeOct | 98.3 | 100.0 | 89.1 |
Schmittel 2014A12 | Triplet | 91.1 | 100.0 | 92.2 |
Merrer 200520 | Cl | 80.0 | 0.0 | 85.5 |
Schmittel 201424 | TSM/N(CH3)2 | 59.7 | 42.5 | 52.7 |
Schmittel 201424 | TMS/OCH3 | 63.2 | 58.7 | 63.3 |
Schmittel 201424 | t Bu/N(CH3)2 | 71.4 | 44.7 | 72.2 |
The latest version of VRAI-selectivity script is available for download on GitHub (https://github.com/sanha0213/VRAI-selectivity).
Footnote |
† Electronic supplementary information (ESI) available: VRAI-selectivity: better calculation of selectivity than transition state theory (PDF). VRAI-selectivity is available for download from Github: https://github.com/sanha0213/VRAI-selectivity. See DOI: 10.1039/d1ob00234a |
This journal is © The Royal Society of Chemistry 2021 |