Kinetic Analysis for the Multistep Profiles of Organic Reactions

Subscriber access provided by UNIV OF WISCONSIN OSHKOSH

Article

Kinetic Analysis for the Multistep Profiles of Organic Reactions: Significance of the Conformational Entropy on the Rate Constants of the Claisen Rearrangement Yosuke Sumiya, Yutaka Nagahata, Tamiki Komatsuzaki, Tetsuya Taketsugu, and Satoshi Maeda J. Phys. Chem. A, Just Accepted Manuscript • DOI: 10.1021/acs.jpca.5b09447 • Publication Date (Web): 16 Nov 2015 Downloaded from http://pubs.acs.org on November 21, 2015

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry A is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 33

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Kinetic Analysis for the Multistep Profiles of Organic Reactions: Significance of the Conformational Entropy on the Rate Constants of the Claisen Rearrangement

Yosuke Sumiya,1 Yutaka Nagahata,2 Tamiki Komatsuzaki,2,3 Tetsuya Taketsugu4 and Satoshi Maeda*4

1

Graduate School of Chemical Sciences and Engineering, Hokkaido University, Kita 13,

Nishi 8, Kita-ku, Sapporo 060-8628, Japan 2

Graduate School of Life Science, Hokkaido University, Kita 10, Nishi 8, Kita-ku,

Sapporo 060-0810, Japan 3

Molecule and Life Nonlinear Sciences Laboratory, Research Institute for Electronic

Science, Hokkaido University, Kita 20 Nishi 10, Kita-ku, Sapporo 001-0020, Japan 4

Department of Chemistry, Faculty of Science, Hokkaido University, Kita 10, Nishi 8,

Kita-ku, Sapporo 060-0810, Japan

ABSTRACT The significance of kinetic analysis as a tool for understanding the reactivity and selectivity of organic reactions has recently been recognized. However, conventional simulation approaches that solve rate equations numerically are not amenable to multistep reaction profiles consisting of fast and slow elementary steps. Herein, we present an efficient and robust approach for evaluating the overall rate constants of multistep reactions via the recursive contraction of the rate equations to give the overall

1

ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

rate constants for the products and byproducts. This new method was applied to the Claisen rearrangement of allyl vinyl ether, as well as a substituted allyl vinyl ether. Notably, the profiles of these reactions contained 23 and 84 local minima, and 66 and 278 transition states, respectively. The overall rate constant for the Claisen rearrangement of allyl vinyl ether was consistent with the experimental value. The selectivity of the Claisen rearrangement reaction has also been assessed using a substituted allyl vinyl ether. The results of this study showed that the conformational entropy in these flexible chain molecules had a substantial impact on the overall rate constants. This new method could therefore be used to estimate the overall rate constants of various other organic reactions involving flexible molecules.

1. INTRODUCTION In recent years, there have been a large number of theoretical studies directed towards the mechanisms of organic reactions.1-3 Many of these previous studies have discussed reactivity and selectivity in terms of the energy barriers associated with important steps in the reaction such as the rate- and selectively-determining elementary steps. Furthermore, the rate constants for these important steps have generally been calculated based on the transition state theory (TST).4,5 However, this approach is sometimes considered to be too simplistic, because the reactivity and selectivity may be affected by the whole reaction profile. Quantitative discussion concerning these systems would therefore require overall rate constants that take into account the entire reaction profile. In principle, an overall rate constant can be evaluated by calculating the rate constants for all of the elementary steps and solving the resulting rate equations in a 2


Page 2 of 33

Page 3 of 33

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


numerical manner for all of the local states. However, it can be difficult to achieve a time evolution covering all of the timescales buried in an entire reaction profile when the elementary steps have a wide range of different rate constants. This situation occurs in the reaction profiles of most organic reactions. For example, the rate constants of local bond rotations are around 1012 s−1, whereas those of bond reorganization steps are around 10−4 s−1. Hence, approaches that simplify rate equations by adopting quasisteady state approximation (QSSA) and/or partial equilibrium assumptions have been employed in recent studies on the kinetics of organic reactions.6,7 Similar problems to this have also been studied in several other fields, including combustion chemistry, where the lumping approach8,9 and its extensions10-12 have been used to address these issues. The lumping approach allows for the reduction of a complicated reaction profile by eliminating unimportant reaction channels and states that rapidly decay to give more stable states. Approaches involving the grouping together of rapidly interconverting states have also been developed and used to evaluate global conformational rearrangements in peptides, as well as phase transition processes in clusters.13,14 The overall rate constants between the resulting superbasins can then be estimated using other methods such as the pre-equilibrium approximation (PEA).15,16 The kinetic Monte Carlo simulation is another powerful approach that provides efficient access to the overall rate constant of a chemical reacion.17,18 This approach effectively samples states that would otherwise require a long-time simulation to be reached by a Monte Carlo approach. In this paper, we have proposed a recursive method for calculating the overall rate constants. The aim of this study was to develop a simple and robust numerical algorithm, and investigate its application to the reaction profiles of organic reactions. This new 3



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

method only requires a reaction profile and a single parameter, namely kMAX, which represents the upper limits of the rate constants in the resulting rate equations. A reaction profile should be prepared prior to the analysis, whilst the kMAX value can be determined easily by adjusting the value based on results of the recursive contraction, as demonstrated below. In this study, reaction profiles were obtained by the singlecomponent artificial force induced reaction (SC-AFIR) method.19 We would like to emphasize that approaches for the construction of reaction profiles are not the main focus of this paper and that several reviews on path searching methods have been provided elsewhere in the literature.20-22

2. THEORY Precondition. A reaction profile consists of a series of reaction pathways with each reaction pathway connecting two different local states. In this paper, we have characterized the local states and reaction pathways as energy minima on the potential energy surface (PES) and intrinsic reaction coordinate (IRC) pathways.23 Let us consider a recursive contraction of the rate equations for a given reaction profile. When a reaction profile has a single product, its final (contracted) rate equation will have only two superstates, i.e., the reactant and the product. When a reaction profile has a distinct product, it will, in principle, remain as another superstate (depending on the timescale of the observation). During each loop of the recursive contraction, one state is contracted to the other neighboring states. In this way, the number of states in the profile can be reduced by one. In this study, we selected the elementary step with the maximum rate constant from the reaction profile, and the contraction was subsequently applied to the initial state 4


Page 4 of 33

Page 5 of 33

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


following this elementary step. The QSSA method is one of the most widely used approaches for the simplification of a reaction profile. The QSSA method was also used in this study for the first step of the contraction. Given that the QSSA simply deletes one state from the reaction profile, the solution to the rate equations at t = ∞ (i.e., the population of each state at thermal equilibrium) would not be conserved. For this reason, we introduced a simple correction scheme (see the Algorithm section) to ensure that the population at thermal equilibrium was conserved throughout the entire contraction procedure. In this study, the rate constant for an elementary step from state X to state Y via transition state TSX‒Y was estimated by the following equation of the TST. k X→Y = Γ

k BT −(∆GTS − ∆GX ) RT e h

1  hν ‡   Γ = 1 +  24  k BT 

(1)

2

(2)

In eq. (1), ∆GX and ∆GTS are the relative Gibbs free energy values for X and TSX‒Y, respectively, kB is the Boltzmann constant, T is the temperature, h is the Planck constant and R is the gas constant. In several applications, the transmission coefficient Γ in Eq. (1) has been set to unity. Several other representations of Γ have also been used to correct for the effects that are neglected in the TST, such as quantum tunneling and trajectory recrossing. The Wigner correction in Eq. (2) has been widely used as a onedimensional tunneling correction,4 where ν‡ is the magnitude of the imaginary frequency at the TS. When Eq. (1) is used for first-order reactions, the population of each state at t = ∞ is equal to the Boltzmann distribution. With this in mind, the Boltzmann distribution was conserved in our scheme (see the Algorithm section).

5



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 33

Algorithm. In the following algorithm, the set for the states obtained at the n-th loop is denoted by N ( n ) ; the i, j, k, l and m states correspond to the components of N ( n ) ( i, j , k , l , m ∈ N ( n ) ); the rate constant for the elementary steps from state i to state

j is denoted by k i(→n ) j ; and the population of state i at t = ∞, i.e., the Boltzmann distribution of state i, is denoted by Pi (n ) . It follows that N (0 ) , ki(→0 ) j and Pi (0 ) correspond to the original sets of all states, rate constants and populations at t = ∞, respectively. The algorithm consists of the following twelve steps. 1.

Input a reaction profile and the parameter kMAX, and initialize n (n = 0).

2.

Compute k i(→0 ) j for all of the elementary steps on the reaction profile using eq. (1)

and obtain the rate equations. When there are two or more elementary steps directly connecting states i and j, k i(→0 ) j represents the sum of the rate constants for these elementary steps. When there is no elementary step directly connecting states i and j,

ki(→0 ) j is zero. The rate constants between the identical states, k i(→0 )i , is also zero. 3.

Compute the population at t = ∞, Pi (0 ) , for all of the states according to their

Boltzmann distribution. ( n)

4.

Identify the state pair i and j with the maximum rate constant ki→ j .

5.

If k i(→n )j < k MAX , then exit from the recursive contraction loop.

6.

(n ) Apply the QSSA to state i. In this step, all of the rate constants kk→l are updated

(n+1) to kk′→l using equations (3) and (4), which are given below. A derivation of these

equations based on the QSSA method is shown in Appendix A:

6


Page 7 of 33

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


k k′ (→n +l 1 ) = k k(n→) l + k k(n→) i k i(→n )lσ i( n )

σ i( n ) =

1

∑

m∈N ( n )

ki(→n )m

(3) (4)

7.

n+1) Update rate constants involving state i to zero ( ki′→(nm+1) = km′(→ i = 0 ).

8.

Set the population of state i at t = ∞ to zero, Pi (n +1 ) = 0 ; in the step 7, the reaction

flux entering state i was set to zero, and Pi (n +1) should therefore equal to zero. 9.

Update Pk(n +1) . In this study, Pi ( n ) , which was updated to zero in step 8, is

distributed to the neighboring states. In contrast, states that rapidly decay to give more stable states were simply grouped or deleted in the previously published approaches.12-14 In this study, by referring the reaction flux per unit time from state i to state k at t = ∞, i.e., k i(→n )k Pi ( n ) , and the normalization constant σ i( n ) in equation (4), part of the Pi ( n ) value, i.e., k i(→n )k σ i( n ) Pi (n ) , is added to Pk(n ) . Pk( n +1) can therefore be obtained by the following equation: Pk( n +1) = Pk( n ) + k i(→n )k σ i( n ) Pi (n )

10.

(5)

(n+1) (n+1) Update kk′→l to kk →l using equation (6), which is provided below. The

(n+1) solutions to the rate equations using kk′→l at t = ∞ do not generally reproduce Pk( n +1) . On

(n+1) the other hand, the solution to the rate equations at t = ∞ using kk →l reproduces Pk(n +1) ,

as demonstrated in Appendix B:

kk(n→+l1) =

1 ( n) (n ) i k →i

1+ σ k

kk′(→n+l 1)

(6)

11.

Define N (n +1) as the set for all of the states except the contracted state i.

12.

Increase n by 1 (n → n + 1), and return to step 4.

7



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 1. A model reaction profile.

An illustration of the procedure described above is provided in Figure 1 for a model reaction profile. In Figure 1, the elementary step from B to C has the maximum rate constant. The QSSA in step 6 would therefore be applied to state B. For this example, k A(0→) C , which would be zero, would be updated to k A′ (1→) C using the following equations. k A′ (1→) C = k A(0→) C + k A(0→) B k B(0→) Cσ B(0 )

σ B(0 ) =

1 (0 )

(0 )

k B → A + k B → C + k B(0→) D

(7) (8)

All of the rate constants concerning the states directly connected to state B would be 1) ′ (1 ) A , k A′ (1→) D , k D→ ′ (1) A , k C→ ′ (1 ) D and k D′ (1→) C ). The updated in a similar manner (i.e., k A′ (→ , k C→ C

other rate constants would simply be copied because the second term in eq. (3) is zero. 1) ′ (1 ) A , k B→ ′ (1 ) C , k C′ (→ ′ (1 ) In the step 7, the rate constants related to state B (i.e., k A′(1→) B , k B→ B , k B→ D 1) and k D′ (→ ) would be set to zero. In the step 8, PB(1) would also be set to zero. PB(0 ) would B

then be distributed to the neighboring states (i.e., states A, C and D). For example, PA(1) would be obtained as follows: PA(1) = PA(0 ) + k B(0→) Aσ B(0 ) PB(0 )

(9)

8


Page 8 of 33

Page 9 of 33

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


This equation therefore indicates that a part of PB(0 ) has been added to PA(0 ) with the (0 ) multiplied by the normalization constant σ B(0) . PC(1) and weight of the rate constant k B→ A

PD(1) can be updated in a similar manner. The populations of the E, F and G states at t = ∞ would not change because the second term in Eq. (5) is zero. In step 10, the rate constants could be further updated to allow for the solutions to the rate equations at t = ∞ to reproduce the populations at t = ∞ in step 9. For example, k A(1→) E would be obtained by the following equation. ) k A(1→ E =

1 (0 ) ( 0 )

1 + σ B kA→B

k A′(1→) E

(10)

A brief explanation of this correction has been provided below (see Appendix B for a detailed explanation). In step 9, the population of state A at t = ∞ (Boltzmann distribution) would increases. This increase would correspond to a decrease in the relative free energy of state A. As a result, the difference in the free energies of state A and the TS would increase and the corresponding rate constants would decrease. In eq. (10), the factor multiplied by k A′(1→) E represents the decrease in the rate constant. In the (1 ) (1 ) step 10, the other rate constants (i.e., k A(1→) C , k A(1→) D , k C→ , (1) , k C→ , k D(1→) A , k D(1→) C and A k C→ D F (1 ) ) would also be updated. Finally, a new state set N k D→ G

(1)

would be obtained by

removing the contracted state B.

3. RESULTS AND DISCUSSION The Claisen rearrangement of allyl vinyl ether CH2-CHO-CH2-CH-CH2 is one of the simplest of all of the known organic rearrangement reactions, and the mechanism of this reaction has been fully elucidated.24 Theoretical calculations have shown that the C-C 9



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

bond formation and C-O bond cleavage steps in this rearrangement occur in a concerted manner with a single TS.25 The results of a gas-phase kinetic study showed that the rate constant of the Claisen rearrangement of allyl vinyl ether follows the Arrhenius equation.26 In this study, we used a reaction profile consisting of 10 and 13 conformers for the wells of the reactant and the product, respectively. From this point onward, the terms “reactant” and “product” shall be used to mean sets of 10 and 13 conformers in the corresponding wells, respectively. The reaction profile was obtained using density functional theory at the B3LYP/6-31G level. All of the local minima (MINs) and transition states (TSs) were subsequently reoptimized at the M062X/6-311+G(2d,p) level.27 Finally, we obtained 23 MINs and 66 TSs at the M062X level. The IRC paths were computed at the M062X level starting from all 66 of the different TSs. Among the 66 IRC paths generated in this way, two corresponded to pathways directly connecting the reactant and the product. These two IRC paths also corresponded to pathways that proceed through chair- and boat-type TSs.25 The other 64 IRC paths corresponded to pathways associated with local conformational changes. To achieve high levels of chemical accuracy, single point CCSD(T)-F12a/jul-cc-pVTZ calculations28,29 were performed for all of the MINs and TSs obtained in the current study. Finally, the Gibbs free energy values were estimated at the experimental temperature (469.1 K) using the electronic energy values at the CCSD(T) level, a rigid rotor approximation based on structures at the M062X level, harmonic vibrational frequencies at the M062X level and the ideal gas approximation. The rate constants for all of the elementary steps were estimated using Eq. (1)

10


Page 10 of 33

Page 11 of 33

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 2. Elementary steps via the chair-type TS (a) and the boat-type TS (b), and the six most stable conformers in the well of the reactant (c) for the Claisen rearrangement of allyl vinyl ether. The relative Gibbs free energy values at 469.1 K are shown in parentheses in kJ mol−1.

The two elementary steps that directly connected one conformer of the reactant with one conformer of the product are shown in Figure 2 (a) and (b). The six most stable conformers in the well of the reactant are also shown in Figure 2 (c). The relative Gibbs free energy values of the different conformers at 469.1 K are shown in parentheses in kJ mol−1. The rate constants for the two elementary steps in Figure 2 (a) and (b) were determined to be 1.454×10−2 and 3.209×10−4 s−1, respectively. When these two direct paths were the only pathways considered for the Claisen rearrangement, the total rate constant was estimated as the sum of the rate constants of these two steps. In this study, we have referred to this approach as the rate-determining elementary step (RDES)

11



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

model, because this model only takes the rate-determining elementary steps into consideration. The rate constant of the RDES model was determined to be 1.486×10−2 s−1, which was very different from that of the experimental value, 2.875×10−3 s−1.26 It was envisaged that the rate constant could be improved by considering the most stable conformer in the well of the reactant. With this in mind, we considered the use of the lowest conformer to transition state (LCTS) model. For the LCTS model, we assumed that the most stable conformer in the well of the reactant shown in Figure 2 (c) was the initial state of the two elementary steps in Figure 2 (a) and (b). In other words, the rate constant was calculated using Eq. (1) based on the Gibbs free energy gaps between the TSs in Figure 2 (a) and (b) and the lowest energy conformer in Figure 2 (c). This model is typically used for the analysis of organic reactions when the most stable (or likely most stable) conformer in the well of the reactant is available. The rate constant determined using the LCTS model, which is the sum of the two LCTS rate constants for the two elementary steps in Figure 2 (a) and (b), was 8.881×10−3 s−1. The value obtained by the LCTS model was therefore closer to the experimental value than that of the RDES model.

Table 1. Rate constants k and the corresponding activation free energies ∆G‡. k / s−1

∆G‡ / kJ mol−1

RDES model a

1.486×10−2

133.1 (133.5) c

LCTS model b

8.881×10−3

135.1 (135.5) c

This work

1.794×10−3

141.3 (141.8) c

Experiment

2.875×10−3

139.5

a

The rate-determining elementary step model (see text). 12


Page 12 of 33

Page 13 of 33

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


b

The lowest conformer to transition state model (see text).

c

Values in parentheses are obtained by Eq. (1) with Γ = 1.

A detailed discussion of the rate constant obtained by the current method is provided below. This rate constant accounts for all of the conformational rearrangement pathways in the well of the reactant. The rate constant obtained using the current method was 1.794×10−3 s−1. Notably, this value was in very good agreement with the experimental value. In this study, the kMAX value was set as 6.931×10−1 s−1; in an irreversible reaction with a k value of 6.931×10−1 s−1, the population of the reactant would decrease to 50% of its original value in one sec with the total number of states being two. This kMAX value was reasonable because only two superstates remained after the recursive contraction and these two corresponded to those for the reactant and product. It is noteworthy that the overall rate constant of the current method did not change when the kMAX value was varied from 1.794×10−3 to 6.061×1010 s−1. The consistency of the rate constant was attributed to the fact that the number of final superstates did not change for any of these variations in the kMAX value for this reaction profile. This result therefore indicated that the reaction profile was not sensitive to the kMAX value. The kMAX value can therefore be decided easily based on the final superstates obtained with a few different kMAX value. The rate constants k obtained using the different models are shown in Table 1, together with the corresponding activation free energies ∆G‡, which were estimated using Eq. (11).

 hk   ∆G ‡ = − RT ln  Γk BT 

(11)

13



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The activation free energies based on the rate constants obtained using Eq. (1) with an Γ value of Eq. (2), and those calculated with an Γ value of 1, were very close to each other. This result indicated that the one-dimensional tunneling correction in Eq. (2) was having very little impact on the rate constant. The rate constants were then systematically improved by considering the most stable conformer (LCTS model vs. RDES model), as well as all of the conformational rearrangement paths (i.e., the results of this work vs. those of the LCTS model). With regard to the activation free energy, the current method reproduced the experimental value with an error of only 1.8 kJ mol−1. The difference between the activation free energies of the current method and the RDES model was 8.2 kJ mol−1, which was attributed to the conformational entropy in the well of the reactant.

Figure 3. Elementary steps via the chair-type TS (a) and the boat-type TS (b), and the six most stable conformers in the well of the reactant (c) for the Claisen rearrangement 14


Page 14 of 33

Page 15 of 33

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


of a substituted allyl vinyl ether. The relative Gibbs free energy values at 469.1 K are shown in parentheses in kJ mol−1

We subsequently applied the current method to the Claisen rearrangement of the substituted allyl vinyl ether Ph-CH-CHO-CH2-CH-CH-Me, where a terminal H atom on the allyl and vinyl groups was substituted with a Me and Ph group, respectively. The reaction profile for this rearrangement was obtained at the M062X/6-311+G(2d,p) level.27 CCSD(T) single point calculations were not performed in this case because an experimental rate constant was not available. Furthermore, the purpose of this part of the study was to demonstrate the applicability of the current method to a reaction profile that also results in the formation of a byproduct. This reaction gives two stereoisomers, including the syn- and anti-products. Using the SC-AFIR, we obtained 17, 35 and 32 conformers for the reactant, syn-product and anti-product, respectively, as well as 278 TSs. Among the 278 TSs, two corresponded to the Claisen rearrangement, whilst the remaining 276 corresponded to conformational rearrangements. This result was confirmed by computing all of the IRC paths starting from the 278 different TSs. Figure 3 (a) and (b) shows the two elementary steps responsible for the formation of the synand anti-products via the chair- and boat-type TSs, respectively. The six most stable conformers in the well of the reactant are also shown in Figure 3 (c). The relative free energy values at 469.1 K are shown in kJ mol−1 in parentheses in Figure 3. These results indicated that the preferred path to the syn-product proceeded via a chair-type TS, and were therefore in agreement with our general understanding of the Claisen rearrangement.

15



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 33

Table 2. Rate constants k and the corresponding activation free energies ∆G‡. Method

k c / s−1

∆G‡ c / kJ mol−1

k d / s−1

∆G‡ d / kJ mol−1

RDES model a

2.231×10−2

131.5 (131.9) e

4.848×10−4

146.4 (146.8) e

LCTS model b

8.327×10−3

135.3 (135.8) e

1.526×10−4

150.9 (151.4) e

This work

1.434×10−3

142.2 (142.6) e

2.601×10−5

157.8 (158.2) e

a

The rate-determining elementary step model (see text).

b

The lowest conformer to transition state model (see text).

c

Values for the syn-product.

d

Values for the anti-product.

e

Values in parentheses are obtained by Eq. (1) with Γ = 1.

The rate constants from the well of the reactant to the two products were estimated by the three different models described above, and the resulting values are shown in Table 2 for comparison. The corresponding activation free energies were estimated using Eq. (11), and the resulting values are also shown in Table 2. The kMAX parameter for these calculations was once again set to 6.931×10−1 s−1. This kMAX value was reasonable also in this system because three superstates remained and these three corresponded to those for the reactant, product (syn-product) and byproduct (antiproduct). The general trends observed in Table 2 were similar to those observed in Table 1. The rate constants obtained by the RDES model were the largest, followed by those of the LCTS model and the present method, which were the smallest. The results with Γ = 1 indicated that the one-dimensional tunneling correction was making a negligible contribution to the calculations. It is noteworthy, that the activation free 16


Page 17 of 33

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


energy determined using the RDES model was significantly lower than the value obtained by the current method, with a difference of 10.7 kJ mol−1. This difference in the two values was attributed to differences in the conformational entropy in the well of the reactant. The impact of this difference on the overall rate constant was accounted for in the current method by considering all of the conformers in the well of the reactant, as well as their interconversion pathways. The importance of sampling the different molecular conformations of the molecules involved in an organic reactions has recently been highlighted in the literature.5,30 The importance of different molecular conformation has also been recognized in enzymatic reactions.31,32 The results of the current study have demonstrated that conformational effects can also have a significant impact on the activation free energy of simple organic molecules. In contrast, the impact of the conformational entropy on the selectivity of this reaction was small. The differences in the activation free energies between the two different products were 14.9 and 15.6 kJ mol−1 in the RDES model and current method, respectively. To let the discussion be more accurate, we estimated the ∆G‡ values after the CCSD(T)-F12a single-point calculations by applying a very simple correction to the profile obtained by the M062X calculations. For this purpose, the relative energy of the two key TSs shown in Fig. 3 (a) and (b) were corrected by the difference between the relative energy obtained for the analogous TSs shown in Fig. 2 (a) and (b) with the M062X and CCSD(T)-F12a levels. This treatment relies on an assumption that level dependences of computational errors in these bond reorganization TSs are similar between the two different systems. The ∆G‡ values were then estimated by the current method or the RDES model using the reaction profile of the M062X calculations with this correction. The correction reduced the ∆G‡ values for reactions giving the syn- and 17



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

anti-product by 6.3 and 8.9 kJ mol−1, respectively, regardless of models (the current method or the RDES model), and our best estimates of the ∆G‡ values for these two channels were therefore 135.9 and 148.9 kJ mol−1, respectively. In each contraction loop, the population at t = ∞ of the contracted state was distributed to the adjacent superstates. It was therefore possible to identify the contribution of each original state to each superstate. For the Claisen rearrangement of allyl vinyl ether, the contributions of all of the original states to the two superstates are shown in Appendix C. The Boltzmann distribution for all of the original states was perfectly reproduced using the equilibrium constant obtained from the overall rate constant and the contributions of the original states to the two superstates (see Appendix C). This condition was explicitly imposed in steps 9 and 10 of the algorithm used in the current study. The current method therefore reproduced the thermodynamic equilibrium despite the fact that is used the QSSA. The correction in Eq. (6), which was introduced to impose this condition, reduced the errors caused by the QSSA. A numerical illustration of this improvement is provided in Appendix D for a simple profile. Further investigation towards developing a deeper understanding of the importance of this correction is currently underway in our group. The accuracy of any method of quantum chemical calculation should be considered as an important performance parameter because poor accuracy can lead to significant errors in the calculations. The approximations used in the Gibbs free energy evaluations, including the harmonic, rigid rotor and ideal gas approximations, can also be significant error sources. In this study, we assumed that the IRC was the reaction path and that the two endpoints of the IRC were the initial and final states of each elementary step. An assumption of this type can lead to large errors in molecular systems that strongly 18


Page 18 of 33

Page 19 of 33

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


deviate from the IRC in terms of their dynamics or systems with dynamic bifurcation pathways.22 The use of Eqs. (1) and (2) to determine the rate constant of a specific transformation, the Born-Oppenheimer approximation, as well as the assumption that the system follows the Markov property, can also result in errors. Whilst we appreciate that there are many different potential sources of error in these calculation, we feel that a detailed study of these different error sources and their impact on the calculations conducted in this study is beyond the scope of this paper. The current method gave a few superstates that were connected by a slow process. This result was similar to those of several previous methods where the original states were grouped together to give only a few superbasins.13,14 The results of the current study were also similar in some ways to those of the reduction methods, where the rapidly converting channels and states are eliminated.10-12 The current method can therefore be considered as an extension of these methods. However, the definition of a superstate used in the current method differed from the definitions used in previous methods, in that the superstates were expressed as the weighted sums of the original states. In other words, all of the original states in the current study made non-zero contributions to all of the superstates. This representation of the superstate would therefore provide a reasonable chemical picture for intermediate states that provide an equal contributing to two superstates. One advantage of the current method is that the overall rate constants can be obtained by simple inputs, i.e., a rough reaction timescale (e.g., 1/kMAX) and a reaction profile. As shown in the second example, the current method allowed for several superstates to remain and provided a platform for discussing the selectivity. That is, when the value of timescale 1/kMAX scans one can capture the hierarchical organization 19



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

among the different states with respect to time. The hierarchical organization of complex reactions with respect to time in reaction networks will be published elsewhere in a separate paper with a new definition of transition state throughout a reaction profile/network.33 The computational cost of the present algorithm was found to be proportional to

Kinetic Analysis for the Multistep Profiles of Organic Reactions

Recommend Documents