# Study of the Deformation/Interaction Model: How Interactions Increase the Reaction Barrier

Authors: Zhiling Liang; Houhe Liu; Nianjun Su; Dandan Song; Yun Zhang; Hong Huang; Jianqi Zheng; Cheng Zhong; Guodong Ye
Journal: Journal of Chemistry
 (2018)

---

## Abstract

The interactions (including weak interactions) between dienophiles and dienes play an important role in the Diels-Alder reaction. To elucidate the influence of these interactions on the reactivity, a popular DFT functional and a variational DFT functional corrected with dispersion terms are used to investigate different substituent groups incorporated on the dienophiles and dienes. The bond order is used to track the trajectory of the cycloaddition reaction. The deformation/interaction model is used to obtain the interaction energy from the reactant complex to the inflection point until reaching the saddle point. The interaction energy initially increases with a decrease in the interatomic distance, reaching a maximum value, but then decreases when the dienophiles and dienes come closer. Reduced density gradient and chemical energy component analysis are used to analyse the interaction. Traditional transition state theory and variational transition state theory are used to obtain the reaction rates. The influence of tunneling on the reaction rate is also discussed.


---

## Body

## 1. Introduction

Interactions play an important role in reactions. Notably, the long-distance weak interaction force, which is known as the weak interaction, has been examined extensively in the study of multicomponent systems, like dimmers, clusters, and so on [1, 2]. However, many chemical reactions involve the breakage and formation of chemical bonds as well as three-dimensional conformational arrangement. Traditional examination methods are unable to cope with such increasing complication. The main reason for the complexity is that two molecules usually react from the reactant complex (RC) to the saddle point (SP) until the product complex (PC) is formed. For example, conformational distortion includes steric strain, torsion, and repulsions in dynamic systems, and the weak interaction also includes hydrogen bonding, van der Waals forces, salt bonding, and hydrophobic interaction forces. This requires expensive treatments (like topological analysis, basis set superposition error correction, etc. [3–8]) to track the trajectory of the reaction and reveal the relationship between the interaction energy (Ei) and the conformational change during the breakage and formation of chemical bonds. Luckily, the deformation/interaction model [9] seems to avoid this complication, compared with the above method, without losing accuracy in the description of Ei.We can describe the change in chemical bonding during formation and dissociation using the Mayer bond order [10]. We can use the reduced density gradient (RDG) function isosurface to analyse the weak interaction region and determine the source of interaction [11, 12]. RDG is an extension of the atoms in molecules (AIM) theory, which gives clearer results. The relationship between energies and intermolecular forces was also elucidated by the energy parameters, according to chemical energy component analysis (CECA) [13], and other kinetic parameters. To better match our experimental results, we used the hybrid DFT functional B3LYP and the variational DFT functional ωB97X-D to obtain the results. The latter contains an empirical dispersion correction to the ωB97X functional, as this provides long-range van der Waals interactions without additional computational cost. It also emerges that optimization of the ωB97 functional with empirical dispersion corrections leads to essentially zero dispersion correction.The Diels-Alder reaction (DA) has a well-known secondary orbital interaction between the diene and the dienophile, which can influence the reaction [14]. The addition of dienophiles to dienes can proceed via a concerted mechanism, meaning that the two new CC bonds form “in concert” rather than sequentially. This popular reaction provides a good and standard example for us to examine the abovementioned interaction. We hope to find good correlation between the activation barrier and activation strain for many systems based on this exploration.


## 2. Calculation Methods

The four reactions with ethylene (ET), amylene (AM), methacrylate (MA), and allyl methyl ether (AME) as dienophiles with butadiene (BD) as the diene were studied. The structures are shown in Figure1. The reactants, transition state (TS) and products were optimized by the Gaussian 09 suite of programs [15] using the B3LYP and ωB97X-D functionals and the 6-311++G (d,p) basis set. Graphical presentation of the weak interaction was achieved by Multiwfn [16], and Ei was decomposed into a variety of energy matrices by APEX4. Reaction rates and tunneling factors were calculated by KiSThelP [17].Figure 1
Schemes of the four Diels-Alder reactions.


## 3. Results and Discussion

### 3.1. Bond Order Analysis

Along the internal reaction coordinate (IRC), we simulated the change in bond order (nMayer) from the RC to the PC. Mayer defined nMayer based on Mulliken population analysis and linear combination of atomic orbitals [18], combining both to derive the formula for nMayer as (1)nMayer=∑μ∈A∑ν∈BPSμνPSνμ,where the off-diagonal matrix elements of PS are known as the Mulliken overlap populations.Figure2(a) shows that these trends are consistent with the empirical view. At the beginning, there is no bond formed between C1 and C6; the value of nMayer is close to 0. When the reaction proceeds to the PC, the value of nMayer between C1 and C6 increases to approximately 1, whereas when C5=C6 is changed to C5-C6, the corresponding value is reduced from 2 to 1. At the TS, the nMayer values of C5-C6 and C1-C6 are approximately 0.4 and 1.4, respectively (not 1.0). The changes in the nMayer value of C5=C6 in the four reactions are relatively similar, and the changes in the bond orders between C1 and C6 are different. This is due to their different structures; that is to say, the AM and AME are affected by the electron-donating group, whereas the MA is affected by the electron-withdrawing group. Meanwhile, the trend in Figure 2(b) is relatively similar.Figure 2
Analysis of the bond order. Bond order curves of BD + ET (black square), BD + AM (cyan star), BD + MA (red diamond), and BD + AME (blue triangle) from the RC to the PC. C1-C6 (- - -), C5=C6 (—). (a) B3LYP and (b) ωB97X-D.


### 3.2. Energy Analysis

We used the deformation/interaction model to analyse dynamic changes in the four DA reactions. The formula for the activation energy is(2)Ea=Ei+Ed,where Ea is the activation energy, Ed is the energy required to distort the olefin reactants into their TS geometries, and Ei arises from a combination of closed-shell repulsion, charge transfer involving occupied and vacant orbital interactions, electrostatic interactions, and polarization effects [19].As shown in Figures3(a) and 3(b), although the energies obtained by ωB97X-D are smaller than those obtained by B3LYP, the two functionals showed similar trends in the calculated energy, and both Ea and Ed increased gradually, while Ei increased initially, reaching a maximum value, but then decreased at a smaller distance. There is an inflection point (IP) in the curve of Ei. Therefore, the process can be divided into two stages: the first stage is from RC to IP, and Ei plays an important role in the activation barrier at this stage, while the second stage is from IP to SP, and Ed exceeds the other chemical forces and governs the mechanism of the reaction. During the course of the reaction, Ea is much higher than Ed, therefore resulting in a positive Ei. This is different from SN2 reactions, in which Ei is negative. It has been proven that the activation barrier decreases relative to Ed in all cases [20]. In particular, the values of Ea, Ei, Ed, and nMayer of BD + AM and BD + AME are similar and higher than those of BD + ET. This result confirms that the dienophile containing an electron-donating group is not conducive to the DA reaction, and while the Ea, Ei, and nMayer values of BD + MA are the smallest, they seem to be affected by the presence of an electron-withdrawing group, which is relatively favourable for the reaction. The symmetry or asymmetry of the structure often affects the Ed value of atoms, and the Ed values of BD in the four reactions are not the same. We found that the Ed value of diene is always greater than that of the dienophiles.Figure 3
Energy analysis by the deformation/interaction model. Changes in the activation energy (…), deformation energy (—), and interaction energy of BD + ET (black square), BD + AM (cyan star), BD + MA (red diamond), and BD + AME (blue triangle) were analysed by the deformation/interaction model and projected onto the distance between C1-C6. (a) B3LYP and (b) ωB97X-D.


### 3.3. Analysis of the Weak Interaction

We used a new method to visualize the weak interaction, which depends on the RGD via the arithmetic expression as shown below:(3)RDG=I023π23∗∇ρrρr43.This part only focuses on ∇ρr/ρr43, where ∇ is the gradient operator and ∇ρr is the modulus of the electron density gradient. By combining the RDG and ρ(r), we can determine which regions of the molecule are clearly associated with weak interactions. The RDG was used to analyse the interactions, which can be divided into medium- and long-range interactions based on the distance or into hydrogen bonding and electrostatic and van der Waals interactions according to the type of bond. A visualization of the weak interaction is shown in Figure 4, using BD + ET as an example. The points on the right side within the abscissa represent the nuclear area. The peak that extends to the upper left corner represents the molecular edge area. When it is farther away from the molecule, ∇ρr and ρ(r) are smaller, while ρ(r) decreases faster, and the RDG value is stronger. The vertical bar, called a “spike” [11], on the leftmost side close to 0 indicates that there is a weak interaction (van der Waals interaction) between C1-C6 and C4-C5 in BD + ET. As is well known, the RDG value is closer to 0, the ρ(r) value is smaller, and the weak interaction is stronger.Figure 4
The RDG region of BD + ET. RC (black square) to IP (red circle) to SP (blue triangle). C1-C6 (a) and C4-C5 (b).
(a)
(b)We also investigated the RDG values of different reactions from the RC to the IP until reaching the SP. The result is shown in Figure5. When ρ(r) is plotted along the X-axis, it gradually increases during the reaction. When the RDG values are plotted along the Y-axis, we found that the RDG values of C1-C6 and C4-C5 in BD + ET are similar. This is because the structure of BD + ET is symmetric; specifically, the RDG value begins to grow slightly in the first stage and decreases to almost the initial value in the second stage, as shown in Figure 5(a). However, the RDG value of C1-C6 in BD + AM and BD + AME decreases in the first stage and increases back to the original value in the second stage. At the same time, the RDG value of C4-C5 in BD + AM increases rapidly, while that of BD + AME increases slowly as the two atoms come closer. It should be noted that the two C6 atoms here are substituted by electron-donating groups. For BD + MA, a continuous decrease in the RDG value between C4-C5 was observed, while the change in the RDG between C1-C6 is the same as that of BD + ET. It seems that the RDG value of C1-C6 changes irregularly, while that of C4-C5 changes regularly. It is noteworthy that C5 is far away from the substituent based on our molecular design. We believed that the RDG value of C5, as a distal end, is more susceptible to the substituent than that of C6, as a near end. The electron-withdrawing group monotonically decreases the RDG value of C4-C5, while the electron-donating group produces the opposite effect.Figure 5
The change in the spike value among the IRC from the RC to the IP until reaching the SP for the four reactions. (a) BD + ET, (b) BD + AM, (c) BD + MA, (d) BD + AME, C1-C6 (white circle), C4-C5 (red circle). B3LYP (—), ωB97X-D (- - -).
(a)
(b)
(c)
(d)In the calculation using theωB97X-D functional, the RGD value of BD + ET increases in the first stage and decreases in the second stage, which is similar to the trend obtained with the B3LYP functional. The RGD value of C1-C6 increases in the first stage and decreases to approximately the original value, while the RGD value of C4-C5 increases monotonously during the two stages. Compared with the results calculated by B3LYP, the calculation using the ωB97X-D functional accounts for the dispersion force to make the results more accurate. These results also prove that the RDG value of C6, as a near end, is more susceptible to the substituent than that of C5, as a distal end. It is noteworthy that the RGD values of C4-C5 in BD + AM and BD + AME are nearly the same. We found that the RGD value of C4-C5 in BD + MA also increases. This is in contrast to the results found using B3LYP and is different from the results of BD + AM and BD + AME. Its amplitude is obviously smaller than that of BD + AM and BD + AME. This validated our hypothesis that the influence of the electron-donating group on increasing the RDG value of C4-C5 is greater than that of the electron-withdrawing group. In brief, if the dienophile is substituted by an electron-withdrawing group, the spike value of a distal end determined by ωB97X-D increases slowly, while value determined by B3LYP shows a reverse trend from the RC to the IP until reaching the SP. This is in contrast to the spike value when the dienophile is substituted by an electron-donating group, which increases quickly as the reaction progresses.


### 3.4. Analysis of the Chemical Energy Component

To explain the above unusual RDG of the BD + MA system, we analysed the chemical energy component. Mayer proposed an energy decomposition method called CECA, in which the sum of the one- and two-center energy components should reproduce the total molecular HF energy. The interaction terms are combined and interpreted as electrostatic, exchange, overlapping, and kinetic energy contribution according to CECA [13]. Note that the two-centered energy here is not equivalent to Ei because it contains other energies (e.g., covalent components). However, we believe it can also elucidate the nature of the molecular interactions.We compare the results of the two functional calculations in Figure6. In Figures 6(a1), 6(b1), 6(c1), and 6(d1), the trends of the energies as the reaction proceeded are almost similar, except for the electrostatic contribution. The trends of C1-C6 are similar and negative, and the trends of C4-C5 are similar and positive in the four systems. However, we still found that BD + MA showed the largest difference (Figure 6(a1)). The electron-withdrawing/donating group influences the inductive effect. Thus, the electrostatic contribution was proven to be related to the nMayer and RDG. The contribution from the exchange and overlap terms continued to decline. We found that the exchange and overlap contributed a large portion of the two-center energy and had a dominating effect. In Figures 6(a2), 6(b2), 6(c2), and 6(d2), through calculation of the dispersion correction, we found that the trends of BD + MA were not consistent in Figures 6(b2), 6(c2), and 6(d2), but this did not affect the system as a whole.Figure 6
Analysis of the electrostatic (a), exchange (b), overlap (c), and kinetic energy (d). BD + ET (black square), BD + AM (cyan star), BD + MA (red diamond), BD + AME (blue triangle). C1-C6 (…), C4-C5 (—). (1) B3LYP; (2) ωB97X-D.


### 3.5. Reaction Rate Results

In this part, the calculation results of the reaction rate constants (k) based on transition state theory (TST) and variational transition state theory (vTST) are discussed [21]. The temperature-dependent reaction rates, kTST (T), for bimolecular reactions were calculated by applying TST: (4)kTSTT=λkbTh·Q≠QAQB·exp⁡-ERT,where Q≠, QA, and QB represent the total partition function of the TS and the reactants A and B, kb represents Boltzmann’s constant, R represents the gas constant, T represents the temperature, h represents Planck’s constant, and E is the activation energy. The constant expression of kvTST(T) is(5)kvTSTT=σhβ·exp⁡-β·ΔGT,s∗vTST,where σ represents the symmetry number and ΔGT,s∗vTST coincides with the maximum free energy, β=kb/T.Since the tunneling factorκ(T) is capable of affecting k, it is necessary to calculate it when k is considered in a precise way. In quantum chemistry, there is still a certain probability to achieve the reaction even when the energy is lower than the barrier. This may be due to penetration of a part of the particle through the middle of the barrier (also known as barrier penetration), which is represented by a correction term in the perturbative expansion of the parabolic barrier. κ(T) can be obtained from the Wigner formula, which is shown below:(6)κT=1+124hω≠kbT2,where ω≠ is the imaginary frequency of the TS. Two virtual frequencies were found here because kTST and kvTST were calculated using different TSs.Table1 shows that the kTST values were similar in magnitude to the kvTST values. This means that the free and electron potential energy surfaces were nearly the same. However, the k obtained by ωB97X-D was greater than that obtained by B3LYP, which matches the results in Figure 3. The kTST and kvTST values of BD + MA were found to be the highest among the four systems. This result agreed with the result that BD + MA had the smallest Ea. However, the kTST value is always of the same order as the kvTST value but lower in the same reaction when determined by ωB97X-D, while the kTST value is not always lower than the kvTST value when determined by B3LYP. It seems that the results obtained by ωB97X-D become more regular through the use of the DFT-D dispersion correction. The κ(T) of all reactions were found to range from 1.2 to 1.4. This increased k to almost the same level, which implies that κ(T) does not significantly affect k. The reaction rate of cyclopentadiene with methyl vinyl ketone was 1.46 × 10–23 cm3 · (molec−1·s−1) [22]. Experiments gave an Arrhenius activation energy for BD + ET of 115.5 kJ/mol at temperatures between 760 K and 921 K [23]. The Ea obtained by ωB97X-D was close to the experimental result.Table 1
The rate constants and tunneling factors.
Reactant
Functional
k
T
S
T
E
a
k
v
T
S
T
ω
≠
κ
(
T
)
(cm3/molec/s)
kJ/mol
(cm3/molec/s)
BD + ET
B3LYP
5.52 × 10−34
151.66
5.40 × 10−34
−544.36
1.29
ωB97X-D
4.46 × 10−32
138.01
6.94 × 10−32
−613.41
1.37
BD + AM
B3LYP
4.37 × 10−35
167.16
4.69 × 10−35
−540.58
1.28
ωB97X-D
3.60 × 10−33
146.69
4.23 × 10−33
−612.82
1.36
BD + MA
B3LYP
2.70 × 10−31
137.97
2.51 × 10−31
−475.59
1.22
ωB97X-D
2.62 × 10−27
109.43
1.28 × 10−26
−543.45
1.29
BD + AME
B3LYP
1.76 × 10−35
162.37
1.44 × 10−35
−534.06
1.28
ωB97X-D
1.81 × 10−32
143.03
2.07 × 10−32
−594.03
1.34


## 3.1. Bond Order Analysis

Along the internal reaction coordinate (IRC), we simulated the change in bond order (nMayer) from the RC to the PC. Mayer defined nMayer based on Mulliken population analysis and linear combination of atomic orbitals [18], combining both to derive the formula for nMayer as (1)nMayer=∑μ∈A∑ν∈BPSμνPSνμ,where the off-diagonal matrix elements of PS are known as the Mulliken overlap populations.Figure2(a) shows that these trends are consistent with the empirical view. At the beginning, there is no bond formed between C1 and C6; the value of nMayer is close to 0. When the reaction proceeds to the PC, the value of nMayer between C1 and C6 increases to approximately 1, whereas when C5=C6 is changed to C5-C6, the corresponding value is reduced from 2 to 1. At the TS, the nMayer values of C5-C6 and C1-C6 are approximately 0.4 and 1.4, respectively (not 1.0). The changes in the nMayer value of C5=C6 in the four reactions are relatively similar, and the changes in the bond orders between C1 and C6 are different. This is due to their different structures; that is to say, the AM and AME are affected by the electron-donating group, whereas the MA is affected by the electron-withdrawing group. Meanwhile, the trend in Figure 2(b) is relatively similar.Figure 2
Analysis of the bond order. Bond order curves of BD + ET (black square), BD + AM (cyan star), BD + MA (red diamond), and BD + AME (blue triangle) from the RC to the PC. C1-C6 (- - -), C5=C6 (—). (a) B3LYP and (b) ωB97X-D.


## 3.2. Energy Analysis

We used the deformation/interaction model to analyse dynamic changes in the four DA reactions. The formula for the activation energy is(2)Ea=Ei+Ed,where Ea is the activation energy, Ed is the energy required to distort the olefin reactants into their TS geometries, and Ei arises from a combination of closed-shell repulsion, charge transfer involving occupied and vacant orbital interactions, electrostatic interactions, and polarization effects [19].As shown in Figures3(a) and 3(b), although the energies obtained by ωB97X-D are smaller than those obtained by B3LYP, the two functionals showed similar trends in the calculated energy, and both Ea and Ed increased gradually, while Ei increased initially, reaching a maximum value, but then decreased at a smaller distance. There is an inflection point (IP) in the curve of Ei. Therefore, the process can be divided into two stages: the first stage is from RC to IP, and Ei plays an important role in the activation barrier at this stage, while the second stage is from IP to SP, and Ed exceeds the other chemical forces and governs the mechanism of the reaction. During the course of the reaction, Ea is much higher than Ed, therefore resulting in a positive Ei. This is different from SN2 reactions, in which Ei is negative. It has been proven that the activation barrier decreases relative to Ed in all cases [20]. In particular, the values of Ea, Ei, Ed, and nMayer of BD + AM and BD + AME are similar and higher than those of BD + ET. This result confirms that the dienophile containing an electron-donating group is not conducive to the DA reaction, and while the Ea, Ei, and nMayer values of BD + MA are the smallest, they seem to be affected by the presence of an electron-withdrawing group, which is relatively favourable for the reaction. The symmetry or asymmetry of the structure often affects the Ed value of atoms, and the Ed values of BD in the four reactions are not the same. We found that the Ed value of diene is always greater than that of the dienophiles.Figure 3
Energy analysis by the deformation/interaction model. Changes in the activation energy (…), deformation energy (—), and interaction energy of BD + ET (black square), BD + AM (cyan star), BD + MA (red diamond), and BD + AME (blue triangle) were analysed by the deformation/interaction model and projected onto the distance between C1-C6. (a) B3LYP and (b) ωB97X-D.


## 3.3. Analysis of the Weak Interaction

We used a new method to visualize the weak interaction, which depends on the RGD via the arithmetic expression as shown below:(3)RDG=I023π23∗∇ρrρr43.This part only focuses on ∇ρr/ρr43, where ∇ is the gradient operator and ∇ρr is the modulus of the electron density gradient. By combining the RDG and ρ(r), we can determine which regions of the molecule are clearly associated with weak interactions. The RDG was used to analyse the interactions, which can be divided into medium- and long-range interactions based on the distance or into hydrogen bonding and electrostatic and van der Waals interactions according to the type of bond. A visualization of the weak interaction is shown in Figure 4, using BD + ET as an example. The points on the right side within the abscissa represent the nuclear area. The peak that extends to the upper left corner represents the molecular edge area. When it is farther away from the molecule, ∇ρr and ρ(r) are smaller, while ρ(r) decreases faster, and the RDG value is stronger. The vertical bar, called a “spike” [11], on the leftmost side close to 0 indicates that there is a weak interaction (van der Waals interaction) between C1-C6 and C4-C5 in BD + ET. As is well known, the RDG value is closer to 0, the ρ(r) value is smaller, and the weak interaction is stronger.Figure 4
The RDG region of BD + ET. RC (black square) to IP (red circle) to SP (blue triangle). C1-C6 (a) and C4-C5 (b).
(a)
(b)We also investigated the RDG values of different reactions from the RC to the IP until reaching the SP. The result is shown in Figure5. When ρ(r) is plotted along the X-axis, it gradually increases during the reaction. When the RDG values are plotted along the Y-axis, we found that the RDG values of C1-C6 and C4-C5 in BD + ET are similar. This is because the structure of BD + ET is symmetric; specifically, the RDG value begins to grow slightly in the first stage and decreases to almost the initial value in the second stage, as shown in Figure 5(a). However, the RDG value of C1-C6 in BD + AM and BD + AME decreases in the first stage and increases back to the original value in the second stage. At the same time, the RDG value of C4-C5 in BD + AM increases rapidly, while that of BD + AME increases slowly as the two atoms come closer. It should be noted that the two C6 atoms here are substituted by electron-donating groups. For BD + MA, a continuous decrease in the RDG value between C4-C5 was observed, while the change in the RDG between C1-C6 is the same as that of BD + ET. It seems that the RDG value of C1-C6 changes irregularly, while that of C4-C5 changes regularly. It is noteworthy that C5 is far away from the substituent based on our molecular design. We believed that the RDG value of C5, as a distal end, is more susceptible to the substituent than that of C6, as a near end. The electron-withdrawing group monotonically decreases the RDG value of C4-C5, while the electron-donating group produces the opposite effect.Figure 5
The change in the spike value among the IRC from the RC to the IP until reaching the SP for the four reactions. (a) BD + ET, (b) BD + AM, (c) BD + MA, (d) BD + AME, C1-C6 (white circle), C4-C5 (red circle). B3LYP (—), ωB97X-D (- - -).
(a)
(b)
(c)
(d)In the calculation using theωB97X-D functional, the RGD value of BD + ET increases in the first stage and decreases in the second stage, which is similar to the trend obtained with the B3LYP functional. The RGD value of C1-C6 increases in the first stage and decreases to approximately the original value, while the RGD value of C4-C5 increases monotonously during the two stages. Compared with the results calculated by B3LYP, the calculation using the ωB97X-D functional accounts for the dispersion force to make the results more accurate. These results also prove that the RDG value of C6, as a near end, is more susceptible to the substituent than that of C5, as a distal end. It is noteworthy that the RGD values of C4-C5 in BD + AM and BD + AME are nearly the same. We found that the RGD value of C4-C5 in BD + MA also increases. This is in contrast to the results found using B3LYP and is different from the results of BD + AM and BD + AME. Its amplitude is obviously smaller than that of BD + AM and BD + AME. This validated our hypothesis that the influence of the electron-donating group on increasing the RDG value of C4-C5 is greater than that of the electron-withdrawing group. In brief, if the dienophile is substituted by an electron-withdrawing group, the spike value of a distal end determined by ωB97X-D increases slowly, while value determined by B3LYP shows a reverse trend from the RC to the IP until reaching the SP. This is in contrast to the spike value when the dienophile is substituted by an electron-donating group, which increases quickly as the reaction progresses.


## 3.4. Analysis of the Chemical Energy Component

To explain the above unusual RDG of the BD + MA system, we analysed the chemical energy component. Mayer proposed an energy decomposition method called CECA, in which the sum of the one- and two-center energy components should reproduce the total molecular HF energy. The interaction terms are combined and interpreted as electrostatic, exchange, overlapping, and kinetic energy contribution according to CECA [13]. Note that the two-centered energy here is not equivalent to Ei because it contains other energies (e.g., covalent components). However, we believe it can also elucidate the nature of the molecular interactions.We compare the results of the two functional calculations in Figure6. In Figures 6(a1), 6(b1), 6(c1), and 6(d1), the trends of the energies as the reaction proceeded are almost similar, except for the electrostatic contribution. The trends of C1-C6 are similar and negative, and the trends of C4-C5 are similar and positive in the four systems. However, we still found that BD + MA showed the largest difference (Figure 6(a1)). The electron-withdrawing/donating group influences the inductive effect. Thus, the electrostatic contribution was proven to be related to the nMayer and RDG. The contribution from the exchange and overlap terms continued to decline. We found that the exchange and overlap contributed a large portion of the two-center energy and had a dominating effect. In Figures 6(a2), 6(b2), 6(c2), and 6(d2), through calculation of the dispersion correction, we found that the trends of BD + MA were not consistent in Figures 6(b2), 6(c2), and 6(d2), but this did not affect the system as a whole.Figure 6
Analysis of the electrostatic (a), exchange (b), overlap (c), and kinetic energy (d). BD + ET (black square), BD + AM (cyan star), BD + MA (red diamond), BD + AME (blue triangle). C1-C6 (…), C4-C5 (—). (1) B3LYP; (2) ωB97X-D.


## 3.5. Reaction Rate Results

In this part, the calculation results of the reaction rate constants (k) based on transition state theory (TST) and variational transition state theory (vTST) are discussed [21]. The temperature-dependent reaction rates, kTST (T), for bimolecular reactions were calculated by applying TST: (4)kTSTT=λkbTh·Q≠QAQB·exp⁡-ERT,where Q≠, QA, and QB represent the total partition function of the TS and the reactants A and B, kb represents Boltzmann’s constant, R represents the gas constant, T represents the temperature, h represents Planck’s constant, and E is the activation energy. The constant expression of kvTST(T) is(5)kvTSTT=σhβ·exp⁡-β·ΔGT,s∗vTST,where σ represents the symmetry number and ΔGT,s∗vTST coincides with the maximum free energy, β=kb/T.Since the tunneling factorκ(T) is capable of affecting k, it is necessary to calculate it when k is considered in a precise way. In quantum chemistry, there is still a certain probability to achieve the reaction even when the energy is lower than the barrier. This may be due to penetration of a part of the particle through the middle of the barrier (also known as barrier penetration), which is represented by a correction term in the perturbative expansion of the parabolic barrier. κ(T) can be obtained from the Wigner formula, which is shown below:(6)κT=1+124hω≠kbT2,where ω≠ is the imaginary frequency of the TS. Two virtual frequencies were found here because kTST and kvTST were calculated using different TSs.Table1 shows that the kTST values were similar in magnitude to the kvTST values. This means that the free and electron potential energy surfaces were nearly the same. However, the k obtained by ωB97X-D was greater than that obtained by B3LYP, which matches the results in Figure 3. The kTST and kvTST values of BD + MA were found to be the highest among the four systems. This result agreed with the result that BD + MA had the smallest Ea. However, the kTST value is always of the same order as the kvTST value but lower in the same reaction when determined by ωB97X-D, while the kTST value is not always lower than the kvTST value when determined by B3LYP. It seems that the results obtained by ωB97X-D become more regular through the use of the DFT-D dispersion correction. The κ(T) of all reactions were found to range from 1.2 to 1.4. This increased k to almost the same level, which implies that κ(T) does not significantly affect k. The reaction rate of cyclopentadiene with methyl vinyl ketone was 1.46 × 10–23 cm3 · (molec−1·s−1) [22]. Experiments gave an Arrhenius activation energy for BD + ET of 115.5 kJ/mol at temperatures between 760 K and 921 K [23]. The Ea obtained by ωB97X-D was close to the experimental result.Table 1
The rate constants and tunneling factors.
Reactant
Functional
k
T
S
T
E
a
k
v
T
S
T
ω
≠
κ
(
T
)
(cm3/molec/s)
kJ/mol
(cm3/molec/s)
BD + ET
B3LYP
5.52 × 10−34
151.66
5.40 × 10−34
−544.36
1.29
ωB97X-D
4.46 × 10−32
138.01
6.94 × 10−32
−613.41
1.37
BD + AM
B3LYP
4.37 × 10−35
167.16
4.69 × 10−35
−540.58
1.28
ωB97X-D
3.60 × 10−33
146.69
4.23 × 10−33
−612.82
1.36
BD + MA
B3LYP
2.70 × 10−31
137.97
2.51 × 10−31
−475.59
1.22
ωB97X-D
2.62 × 10−27
109.43
1.28 × 10−26
−543.45
1.29
BD + AME
B3LYP
1.76 × 10−35
162.37
1.44 × 10−35
−534.06
1.28
ωB97X-D
1.81 × 10−32
143.03
2.07 × 10−32
−594.03
1.34


## 4. Conclusions

The incorporation of an electron-withdrawing group affects the change innMayer. Ei is the dominating energy when the cycloaddition reaction proceeds from the RC to the IP. The RDG analysis found that a dienophile substituted by an electron-withdrawing group decreased the spike value of the distal end. CECA also revealed that exchange and overlap terms provided the largest contributions to the two-center energy. Finally, k was related to Ea, but κ(T) increased in equal proportion to k.



---

*Source: 3106297-2018-01-14.xml*