Journal of Chemistry
Volume 2018 (2018), Article ID 3106297, 8 pages
https://doi.org/10.1155/2018/3106297

Research Article

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

Zhiling Liang,¹ Houhe Liu,¹ Nianjun Su,² Dandan Song,¹ Yun Zhang,¹ Hong Huang,¹ Jianqi Zheng,³ Cheng Zhong,⁴ and Guodong Ye¹

¹Guangdong Provincial Key Laboratory of Molecular Target & Clinical Pharmacology, School of Pharmaceutical Sciences, Guangzhou Medical University and The Fifth Affiliated Hospital of Guangzhou Medical University, Guangzhou 510700, China
²Department of Reproductive Health and Infertility, Guangdong Women and Children Hospital, Guangzhou 510010, China
³The First Clinical School of Guangzhou Medical University, Guangzhou 510275, China
⁴College of Chemistry and Molecular Sciences, Wuhan University, Wuhan 430072, China

Correspondence should be addressed to Guodong Ye; gzhygd@gzhmu.edu.cn

Received 9 June 2017; Revised 28 September 2017; Accepted 31 October 2017; Published 14 January 2018

Academic Editor: Arturo Espinosa Ferao

Copyright © 2018 Zhiling Liang et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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].

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 () 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 C₁ and C₆; the value of is close to 0. When the reaction proceeds to the PC, the value of between C₁ and C₆ increases to approximately 1, whereas when C₅=C₆ is changed to C₅-C₆, the corresponding value is reduced from 2 to 1. At the TS, the values of C₅-C₆ and C₁-C₆ are approximately 0.4 and 1.4, respectively (not 1.0). The changes in the value of C₅=C₆ in the four reactions are relatively similar, and the changes in the bond orders between C₁ and C₆ 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. C₁-C₆ (- - -), C₅=C₆ (—). (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

where 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.

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 C₁-C₆. (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:

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 C₁-C₆ and C₄-C₅ 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)

Figure 4: The RDG region of BD + ET. RC (black square) to IP (red circle) to SP (blue triangle). C₁-C₆ (a) and C₄-C₅ (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 C₁-C₆ and C₄-C₅ 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 C₁-C₆ 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 C₄-C₅ 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 C₆ atoms here are substituted by electron-donating groups. For BD + MA, a continuous decrease in the RDG value between C₄-C₅ was observed, while the change in the RDG between C₁-C₆ is the same as that of BD + ET. It seems that the RDG value of C₁-C₆ changes irregularly, while that of C₄-C₅ changes regularly. It is noteworthy that C₅ is far away from the substituent based on our molecular design. We believed that the RDG value of C₅, as a distal end, is more susceptible to the substituent than that of C₆, as a near end. The electron-withdrawing group monotonically decreases the RDG value of C₄-C₅, while the electron-donating group produces the opposite effect.

(a)

(b)

(c)

(d)

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, C₁-C₆ (white circle), C₄-C₅ (red circle). B3LYP (—), ωB97X-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 C₁-C₆ increases in the first stage and decreases to approximately the original value, while the RGD value of C₄-C₅ 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 C₆, as a near end, is more susceptible to the substituent than that of C₅, as a distal end. It is noteworthy that the RGD values of C₄-C₅ in BD + AM and BD + AME are nearly the same. We found that the RGD value of C₄-C₅ 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 C₄-C₅ 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 C₁-C₆ are similar and negative, and the trends of C₄-C₅ 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.

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). C₁-C₆ (…), C₄-C₅ (—). (1) B3LYP; (2) ωB97X-D.

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 is

where 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 cm³ · (molec⁻¹·s⁻¹) [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.

Table 1: The rate constants and tunneling factors.

Reactant
Functional

(cm³/molec/s)
kJ/mol
(cm³/molec/s)

BD + ET
B3LYP
5.52 × 10⁻³⁴
151.66
5.40 × 10⁻³⁴
−544.36
1.29

B97X-D
4.46 × 10⁻³²
138.01
6.94 × 10⁻³²
−613.41
1.37

BD + AM
B3LYP
4.37 × 10⁻³⁵
167.16
4.69 × 10⁻³⁵
−540.58
1.28

B97X-D
3.60 × 10⁻³³
146.69
4.23 × 10⁻³³
−612.82
1.36

BD + MA
B3LYP
2.70 × 10⁻³¹
137.97
2.51 × 10⁻³¹
−475.59
1.22

B97X-D
2.62 × 10⁻²⁷
109.43
1.28 × 10⁻²⁶
−543.45
1.29

BD + AME
B3LYP
1.76 × 10⁻³⁵
162.37
1.44 × 10⁻³⁵
−534.06
1.28

B97X-D
1.81 × 10⁻³²
143.03
2.07 × 10⁻³²
−594.03
1.34

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)

References

E. F. Kiefer and W. Gericke, “Nuclear magnetic resonance investigation of secondary valence forces. II. Intramolecular mercury-halogen coordination in 3-halopropylmercury compounds,” Journal of the American Chemical Society, vol. 90, no. 19, pp. 5131–5136, 1968.

X. An, H. Zhuo, W. Tian, Q. Li, W. Li, and J. Cheng, “Novel non-covalent interactions involved with the Al13M cluster (M = Li, Na, K, Cu, Ag, Au),” Molecular Physics, vol. 112, no. 22, pp. 2954–2962, 2014.

R. Chaudret, B. De Courcy, J. Contreras-García et al., “Unraveling non-covalent interactions within flexible biomolecules: From electron density topology to gas phase spectroscopy,” Physical Chemistry Chemical Physics, vol. 16, no. 21, pp. 9876–9891, 2014.

D. J. R. Duarte, G. L. Sosa, and N. M. Peruchena, “Nature of halogen bonding. A study based on the topological analysis of the Laplacian of the electron charge density and an energy decomposition analysis,” Journal of Molecular Modeling, vol. 19, no. 5, pp. 2035–2041, 2013.

D. Mani and E. Arunan, “The X-C⋯Y (X = O/F, y = O/S/F/Cl/Br/N/P) 'carbon bond' and hydrophobic interactions,” Physical Chemistry Chemical Physics, vol. 15, no. 34, pp. 14377–14383, 2013.

D. Hugas, L. Guillaumes, M. Duran, and S. Simon, “Delocalization indices for non-covalent interaction: Hydrogen and DiHydrogen bond,” Computational and Theoretical Chemistry, vol. 998, pp. 113–119, 2012.

F. Fuster and S. J. Grabowski, “Intramolecular hydrogen bonds: The QTAIM and ELF characteristics,” The Journal of Physical Chemistry A, vol. 115, no. 35, pp. 10078–10086, 2011.

S. Grimme, C. Mück-Lichtenfeld, and J. Antony, “Noncovalent interactions between graphene sheets and in multishell (hyper)fullerenes,” The Journal of Physical Chemistry C, vol. 111, no. 30, pp. 11199–11207, 2007.

S. Osuna and K. N. Houk, “Cycloaddition reactions of butadiene and 1,3-dipoles to curved arenes, fullerenes, and nanotubes: Theoretical evaluation of the role of distortion energies on activation barriers,” Chemistry - A European Journal, vol. 15, no. 47, pp. 13219–13231, 2009.

A. J. Bridgeman, G. Cavigliasso, L. R. Ireland, and J. Rothery, “The Mayer bond order as a tool in inorganic chemistry,” Journal of the Chemical Society, Dalton Transactions, no. 14, pp. 2095–2108, 2001.

E. R. Johnson, S. Keinan, P. Mori-Sánchez, J. Contreras-García, A. J. Cohen, and W. Yang, “Revealing noncovalent interactions,” Journal of the American Chemical Society, vol. 132, no. 18, pp. 6498–6506, 2010.

N. Gillet, R. Chaudret, J. Contreras-García, W. Yang, B. Silvi, and J.-P. Piquemal, “Coupling quantum interpretative techniques: Another look at chemical mechanisms in organic reactions,” Journal of Chemical Theory and Computation, vol. 8, no. 11, pp. 3993–3997, 2012.

A. Hamza and I. Mayer, “Physical analysis of the diatomic "chemical" energy components,” Theoretical Chemistry Accounts, vol. 109, no. 2, pp. 91–98, 2003.

E. J. Corey and T. Loh, “First application of attractive intramolecular interactions to the design of chiral catalysts for highly enantioselective Diels-Alder reactions,” Journal of the American Chemical Society, vol. 32, no. 23, pp. 6289–6292, 1991.

M. J. Frisch, G. W. Trucks, and H. B. Schlegel, 2010, Gaussian 09. Gaussian Inc., Wallingford.

T. Lu and F. Chen, “Multiwfn: a multifunctional wavefunction analyzer,” Journal of Computational Chemistry, vol. 33, no. 5, pp. 580–592, 2012.

S. Canneaux, F. Bohr, and E. Henon, “KiSThelP: A program to predict thermodynamic properties and rate constants from quantum chemistry results,” Journal of Computational Chemistry, vol. 35, no. 1, pp. 82–93, 2014.

I. Mayer, “On bond orders and valences in the Ab initio quantum chemical theory,” International Journal of Quantum Chemistry, vol. 29, no. 1, pp. 73–84, 1986.

C. A. Sader and K. N. Houk, “Distortion/interaction analysis of the reactivities and selectivities of halo- and methoxy-substituted carbenes with alkenes,” Arkivoc Archive for Organic Chemistry, vol. 77, no. 77, pp. 843–850, 2014.

D. H. Ess and K. N. Houk, “Theory of 1,3-dipolar cycloadditions: Distortion/interaction and frontier molecular orbital models,” Journal of the American Chemical Society, vol. 130, no. 31, pp. 10187–10198, 2008.

A. Fernandez-Ramos, B. A. Ellingson, B. C. Garrett, and D. G. Truhlar, “Variational transition state theory with multidimensional tunneling,” Reviews in Computational Chemistry, vol. 23, pp. 125–232, 2007.

L. L. Thomas, J. Tirado-Rives, and W. L. Jorgensen, “Quantum mechanical/molecular mechanical modeling finds diels alder reactions are accelerated less on the surface of water than in water,” Journal of the American Chemical Society, vol. 132, no. 9, pp. 3097–3104, 2010.

J. Sauer and R. Sustmann, “Mechanistic Aspects of Diels‐Alder Reactions: A Critical Survey,” Angewandte Chemie International Edition, vol. 19, no. 10, pp. 779–807, 1980.