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.
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 () 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 .
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 Figure 1. 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 was decomposed into a variety of energy matrices by APEX4. Reaction rates and tunneling factors were calculated by KiSThelP [17].
3. Results and Discussion
3.1. Bond Order Analysis
Along the internal reaction coordinate (IRC), we simulated the change in bond order () from the RC to the PC. Mayer defined based on Mulliken population analysis and linear combination of atomic orbitals [18], combining both to derive the formula for as where the off-diagonal matrix elements of PS are known as the Mulliken overlap populations.
Figure 2(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 is close to 0. When the reaction proceeds to the PC, the value of 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 values of C5-C6 and C1-C6 are approximately 0.4 and 1.4, respectively (not 1.0). The changes in the 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.
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 iswhere is the activation energy, is the energy required to distort the olefin reactants into their TS geometries, and 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 Figures 3(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 and increased gradually, while increased initially, reaching a maximum value, but then decreased at a smaller distance. There is an inflection point (IP) in the curve of . Therefore, the process can be divided into two stages: the first stage is from RC to IP, and plays an important role in the activation barrier at this stage, while the second stage is from IP to SP, and exceeds the other chemical forces and governs the mechanism of the reaction. During the course of the reaction, is much higher than , therefore resulting in a positive . This is different from 2 reactions, in which is negative. It has been proven that the activation barrier decreases relative to in all cases [20]. In particular, the values of , , , and 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 , , and 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 value of atoms, and the values of BD in the four reactions are not the same. We found that the value of diene is always greater than that of the dienophiles.
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:This part only focuses on , where is the gradient operator and is the modulus of the electron density gradient. By combining the RDG and , 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, and are smaller, while 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 value is smaller, and the weak interaction is stronger.
(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 Figure 5. When is plotted along the -axis, it gradually increases during the reaction. When the RDG values are plotted along the -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.
(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 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 Figure 6. 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 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.
3.5. Reaction Rate Results
In this part, the calculation results of the reaction rate constants () based on transition state theory (TST) and variational transition state theory (TST) are discussed [21]. The temperature-dependent reaction rates, (T), for bimolecular reactions were calculated by applying TST: where , , and represent the total partition function of the TS and the reactants and , represents Boltzmann’s constant, represents the gas constant, represents the temperature, represents Planck’s constant, and is the activation energy. The constant expression of iswhere represents the symmetry number and coincides with the maximum free energy, .
Since the tunneling factor is capable of affecting , it is necessary to calculate it when 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. can be obtained from the Wigner formula, which is shown below:where is the imaginary frequency of the TS. Two virtual frequencies were found here because and were calculated using different TSs.
Table 1 shows that the values were similar in magnitude to the values. This means that the free and electron potential energy surfaces were nearly the same. However, the obtained by ωB97X-D was greater than that obtained by B3LYP, which matches the results in Figure 3. The and 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 . However, the value is always of the same order as the value but lower in the same reaction when determined by ωB97X-D, while the value is not always lower than the 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 of all reactions were found to range from 1.2 to 1.4. This increased to almost the same level, which implies that does not significantly affect . 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 obtained by ωB97X-D was close to the experimental result.
4. Conclusions
The incorporation of an electron-withdrawing group affects the change in . 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, was related to , but increased in equal proportion to .
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Authors’ Contributions
Zhiling Liang, Houhe Liu, and Nianjun Su contributed to this work equally and should be regarded as co-first authors.
Acknowledgments
The authors’ works are supported by the National Natural Science Foundation of China (Grant no. 21274032) and the Natural Science Foundation of Guangdong Province (Grant no. 2014A030313500). It was also funded by the project of the high level university of Guangdong province. The corresponding author, Guodong Ye, thanks the above funding sources.
Supplementary Materials
The coordinates of the RCs, PCs, and their TSs are given in the Supplementary Material. (Supplementary Materials)