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

1Guangdong 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 2Department of Reproductive Health and Infertility, Guangdong Women and Children Hospital, Guangzhou 510010, China 3The First Clinical School of Guangzhou Medical University, Guangzhou 510275, China 4College of Chemistry and Molecular Sciences, Wuhan University, Wuhan 430072, China


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][4][5][6][7][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.

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

Results and Discussion
3.1.Bond Order Analysis.Along the internal reaction coordinate (IRC), we simulated the change in bond order ( Mayer ) from the RC to the PC.Mayer defined  Mayer based on Mulliken population analysis and linear combination of atomic orbitals [18], combining both to derive the formula for  Mayer 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 1 and C 6 ; the value of  Mayer is close to 0. When the reaction proceeds to the PC, the value of  Mayer between C 1 and C 6 increases to approximately 1, whereas when C 5 =C 6 is changed to C 5 -C 6 , the corresponding value is reduced from 2 to 1.At the TS, the  Mayer values of C 5 -C 6 and C 1 -C 6 are approximately 0.4 and 1.4, respectively (not 1.0).The changes in the  Mayer value of C 5 =C 6 in the four reactions are relatively similar, and the changes in the bond orders between C 1 and C 6 are different.This is due to their different structures; that is to say, the AM and AME are affected by the electrondonating group, whereas the MA is affected by the electronwithdrawing group.Meanwhile, the trend in Figure 2(b) is relatively similar.

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 S N 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  Mayer 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  Mayer 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.

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 |∇()|/ 3 √() 4 , 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 longrange 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 1 -C 6 and C 4 -C 5 in BD + ET.As is well known, the RDG value is closer to 0, the () value is smaller, and the weak interaction is stronger.
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 1 -C 6 and C 4 -C 5 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 1 -C 6 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 4 -C 5 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 6 atoms here are substituted by electron-donating groups.For BD + MA, a continuous decrease in the RDG value between C 4 -C 5 was observed, while the change in the RDG between C 1 -C 6 is the same as that of BD + ET.It seems that the RDG value of C 1 -C 6 changes irregularly, while that of C 4 -C 5 changes regularly.It is noteworthy that C 5 is far away from the substituent based on our molecular design.We believed that the RDG value of C 5 , as a distal end, is more susceptible to the substituent than that of C 6 , as a near end.The electronwithdrawing group monotonically decreases the RDG value of C 4 -C 5 , while the electron-donating group produces the opposite effect.
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 1 -C 6 increases in the first stage and decreases to approximately the original value, while the RGD value of C 4 -C 5 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 6 , as a near end, is more susceptible to the substituent than that of C 5 , as a distal end.It is noteworthy that the RGD values of C 4 -C 5 in BD + AM and BD + AME are nearly the same.We found that the RGD value of C 4 -C 5 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 4 -C 5 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.

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 1 -C 6 are similar and negative, and the trends of C 4 -C 5 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  Mayer 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.

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 (VTST) are discussed [21].The temperature-dependent reaction rates,  TST (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  VTST () is where  represents the symmetry number and Δ(,  VTST * ) 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  TST and  VTST were calculated using different TSs.
IRC coordinate (；ＧＯ   Table 1 shows that the  TST values were similar in magnitude to the  VTST 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  TST and  VTST 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  TST value is always of the same order as the  VTST value but lower in the same reaction when determined by B97X-D, while the  TST value is not always lower than the  VTST 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 3 ⋅ (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.

Figure 3 :
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 1 -C 6 .(a) B3LYP and (b) B97X-D.

Table 1 :
The rate constants and tunneling factors.