Solvent-Mediated Rate Deceleration of Diels–Alder Reactions for Enhanced Selectivity: Quantum Mechanical Insights

Abstract

Solvents can have a tremendous influence on the rate and selectivity of chemical reactions, but their effects are not always well accounted for. In the present work, density functional theory computations are used to investigate the influence of solvent on the Diels–Alder reactions of 9-methylanthracene with (5-oxo-2H-furan-2-yl) acetate and different anhydrides considering the overall reaction rates as well as selectivity between possible isomeric products. Crucially, we find that overall reaction rates are higher in non-polar toluene, whereas selectivity is enhanced in the polar solvent acetone. In the case of (5-oxo-2H-furan-2-yl) acetate, the difference in the reaction barriers is enhanced by 2.4 kJ/mol in acetone as compared to the gas phase, halving the yield of the side product. Similar results are found for the reaction of 9-methylanthracene with chloro-maleic anhydride and cyano-maleic anhydride, highlighting the generality of the trends observed. After presenting the energetics, a detailed discussion of the reactivity is given using electrostatic potentials, frontier orbitals, reactivity indices and Fukui functions. In summary, this study highlights the importance of solvent in influencing reaction rates and illustrates the possibility of studying its effects computationally.

1. Introduction

The Diels–Alder (DA) reaction is an indispensable chemical reaction in modern as well as traditional organic chemistry and it has found many applications in the production of several important chemical compounds [1,2,3]. For instance, the DA reaction is used in the synthesis of bioactive natural products [4,5,6], in polymer chemistry [7,8] and medicinal chemistry [9,10] and for the functionalization of fullerenes [11]. Furthermore, this reaction has a major role in biological systems [12], and the DA reactions involving furans as dienes are a very important synthetic tool in green chemistry [13,14]. Keeping in view the aforementioned and many other uses of the DA reaction in a wide range of sectors, it is necessary to have a comprehensive understanding of the mechanisms of this reaction. Since the discovery of the DA reaction, its mechanism has been investigated by a large number of researchers, but even today there is an enthusiastic debate about it. This reaction is a facile approach for the formation of new C-C bonds by the addition of an electron-rich diene to an electron-poor dienophile [15,16,17,18]. The processes involving bond formation and bond breaking in the DA reactions are thought to be concerted, but they are not always synchronous. In extreme examples, it may become a two-step chemical reaction with biradical or zwitterionic intermediates [19]. Even though the mechanism of the DA reaction has been much explored both theoretically and experimentally, some of its classes, such as the DA reactions of furoic acids, the aza-Diels–Alder reaction and the hetero-Diels–Alder (HDA) reaction, are not investigated enough [9,14,20].

Solvents are generally supposed to be a medium in which the chemical reaction between the dissolved molecules of solute takes place [21]. However, molecules of the solvents used in a chemical reaction also have a great influence on its mechanism and kinetics, in addition to bringing the reacting species into contact [22,23]. Changing the solvent for a given reaction can modify its activation barrier [24], thermodynamics, equilibrium composition [25], selectivity, rate constant [26] and its enantioselectivity [27]. As an extreme example, the reaction rate of a hydroxide anion with methyl formate is decreased by a factor of one billion when the reaction is carried out in water instead of the gas phase due to its stabilization and high solvation [28]. Solvents play a major role in pharmaceutical, electrochemical, photochemical and many other industrial processes. Hence, it is necessary to have an excellent solvent for maximizing the yield and selectivity of the desired product and modifying the reaction rate [29]. Experimental screening of the most effective and efficient solvent for a particular reaction is usually very expensive and restricted to a few solvents. Therefore, the theoretical approach is very useful for accelerating and improving the search for the best solvent [30,31]. Advancement in the field of computational chemistry has allowed the deep examination of the transition-state structures involved in the reactions occurring in solvents, and the computational studies on solution-phase reactions have revealed that solvents have a far more critical role in chemical reactions.

Theoretically, solvation effects can be taken into account via either implicit or explicit models. Explicit solvation models treat solvent molecules individually and thus can provide details of different interactions, including hydrogen bonding between solute and solvent molecules. However, modelling the solvent effects explicitly can significantly increase computation cost [32] while also posing challenges in terms of accurate sampling. On the other hand, in implicit solvation models, solvents are treated as a continuous medium which surrounds the solute molecules. These models offer the advantage of low computational cost but lack atomistic detail [33]. In this paper, an implicit solvation model, the integral equation formalism polarisable continuum model (IEFPCM), is employed, as it has been used extensively for studying Diels–Alder reactions [34,35,36,37,38,39], providing reliable results [40].

In the present work, the mechanism of the hetero-DA reactions between 9-methyl-anthracene (1) and four different dienophiles, i.e., maleic anhydride (2a), chloro-maleic anhydride (2b), cyano-maleic anhydride (2c) and (5-oxo-2H-furan-2-yl) acetate (2d) (see Scheme 1), is theoretically investigated in the gas phase, toluene and acetone. According to previous studies, some Diels–Alder reactions are accelerated in the presence of a solvent while others are decelerated [23]. The purpose of the current study is to evaluate whether the rate of the selected reactions is promoted or retarded by the use of these solvents and, specifically, how selectivity can be modulated. The activation barrier of the reactions is compared in the solvents of different polarities. Furthermore, analyses for molecular electrostatic potential (MESP), frontier molecular orbitals (FMO), global reactivity indices (GRI), natural bond orbitals (NBO) and Fukui functions are performed for reactions of 2a and 2d. We focus our study on 2a and 2d while selected results are also presented for 2b and 2c to test the generality of the results. This work will be helpful for gaining insight into the factors responsible for changing the properties of DA reactions due to variations in the environment of chemical species involved in the reaction. This understanding helps in deciding the chemical environment for the synthesis of a wide range of chemicals including biologically important molecules, agrochemicals, polymers and hybrid materials.

2. Computational Details

Density functional theory (DFT) calculations on all the chemical species were conducted using the Gaussian 09 program, and, for the visualisation of structures, Gauss-View 6.0.16 was employed [41]. The geometry of all structures was optimised at the M06-2X/6-31+G(d) level of theory [42,43], and frequency analysis was performed at the same level while single-point energy calculations were carried out with the 6-311+G(d,p) basis set. Optimisations and single-point computations were performed in the respective solvent using the integral equation formalism polarisable continuum model (IEFPCM) [44]. The M06-2X functional was previously shown to provide a reliable description of the energetics and geometries for Diels–Alder reactions, performing better than, e.g., the popular B3LYP functional [45,46].

Analysis of normal vibrational modes confirmed that all the reactants and products are minima with real frequencies only and all the transition state structures have single imaginary frequencies. Intrinsic reaction coordinate (IRC) [47] calculations were performed on the transition states in the gas phase to ensure they connected the reactants to the products (Figures S1–S3). Free energies were computed using a standard rigid rotor/harmonic oscillator model based on harmonic frequencies at the equilibrium geometries. Rates were computed using standard transition-state theory.

For molecules 3d and 3d’, a conformer search was performed using CREST version 2.12 (conformer–rotamer ensemble sampling tool) [48,49] and then the lowest-energy structures were selected for further optimisation. For the other molecules, no conformer analysis was performed due to their relative rigidity.

Energies of frontier molecular orbitals (FMO) were computed for all the molecules, and their energies were used to derive chemical hardness ($\eta$), chemical softness (S), the global electrophilicity index ($\omega$) and chemical potential ($\mu$) [50]. These parameters are computed using the following expressions [51,52,53]:

$$\eta ={\displaystyle \frac{E_{LUMO}-E_{HOMO}}{2}}$$

(1)

$$S={\displaystyle \frac{1}{\eta }}$$

(2)

$$\mu ={\displaystyle \frac{1}{2}}(E_{LUMO}+E_{HOMO})$$

(3)

$$\omega ={\displaystyle \frac{{\mu }^2}{2\eta }}$$

(4)

For calculating the nucleophilicity index (N) of the molecules, the energy for the highest occupied molecular orbital of tetracyanoethylene (TCNE), $E_{HOMO}^{\mathrm{TCNE}}$, is subtracted from the $E_{HOMO}^{\mathrm{nuc}}$ of the selected molecule [54]. TCNE is taken as a reference because of its very low-lying HOMO energy [55].

$$N=E_{HOMO}^{\mathrm{nuc}}-E_{HOMO}^{\mathrm{TCNE}}$$

(5)

Condensed-to-atom Fukui functions for molecules were calculated following refs. [56,57] using

$$f_w^{+}=q_w\left(N\right)-q_w(N+1)$$

(6)

for the attack of the nucleophile and

$$f_w^{-}=q_w(N-1)-q_w\left(N\right)$$

(7)

for the attack of electrophile. In these equations, $q_w$ represents the charge on the site w while N, $(N+1)$ and $(N-1)$ represent neutral, anionic and cationic systems, respectively [58].

3. Results and Discussion

3.1. Reaction Thermodynamics

The thermodynamic properties of the chemical species involved in the reaction can provide assistance in predicting whether a chemical reaction will occur spontaneously under a given set of conditions. Moreover, thermodynamic parameters help in deciding the reaction mechanism, determining the relative stability of molecules and measuring the effect of solvents, catalysts, etc., on chemical reactions [59]. Figure 1 displays the activation barriers for the reaction of 9-methylanthracene (1) with the four dienophiles. We start the discussion with 2a and 2d, which are the focus of this study. The DA reaction between 1 and 2a has a reasonably low gas-phase single-point free-energy barrier of 80.2 kJ/mol. The reaction is also strongly exergonic (–78.8 kJ/mol), highlighting that there is a strong thermodynamic driving force. A significant increase in the reaction barrier is observed when considering 2d as the dienophile. This dienophile can approach the diene in different orientations, resulting in the formation of isomeric products 3d and 3d’. In product 3d, the methyl group on the diene and the acetoxy group on the dienophile are in the close vicinity of one another, while, in 3d’, these two groups are far from one another. The orientation of the dienophile in the case of 3d’ is kinetically more favourable, as indicated by the lower activation free energy for its synthesis compared to that of 3d. This is consistent with the experimental results when the reaction was carried out in toluene and the yield of 3d’ was considerably higher than that of 3d [60]. The difference in the value of the activation barrier for these two isomers in the gas phase is about 3.1 kJ/mol. This increase in the reaction barrier can be attributed to steric hindrance between two bulky groups, i.e., methyl and acetoxy [61]. For these reactions, the enthalpy and free energy of products are lower than those of the reactants; hence, these are exothermic and exergonic at room temperature in the gas phase [62].

For theoretical evaluation of the solvent effects on the yield of 3d and 3d’, we have determined it in five different solvents with different dielectric constants, i.e., DMSO (dimethylsulfoxide), acetone, DCM (dichloromethane), diethyl ether, toluene and n-hexane. The results are presented in Figure 2. The qualitative effect of all solvents is the same, i.e., increased free energy of activation when compared to the gas phase, as presented above. This overall effect may be because reactants have a higher accessible surface, and hence their stabilization by solvents is greater than that of the products. Crucially, solvation not only raises the barrier heights, but it also raises the difference in barrier heights. In the gas phase, the $\Delta \Delta G^{\ne }$ for the two reaction pathways is only 3.1 kJ/mol, which corresponds to a product distribution of 78:22 according to standard Boltzmann statistics. The selectivity is slightly improved in toluene ($\Delta \Delta G^{\ne }=3.7$ kJ/mol) to 82:18, while acetone gives ($\Delta \Delta G^{\ne }=5.5$ kJ/mol) a significant further improvement to 90:10 (Figure 3). In DMSO solvent, which has a higher dielectric constant than acetone, this selectivity is enhanced to 91.5:8.5. Figure 2 highlights a roughly linear relationship between the inverse of the dielectric constant and the yield of the side product 3d. Conversely, more erratic behaviour and no clear trends are observed for the overall barrier heights. In summary, the results highlight that the yield of the side product 3d is reduced by more than half when changing the solvent from n-hexane to acetone or DMSO.

In order to test the general validity of the enhanced reaction barriers, we investigated the Diels–Alder reactions of the same diene with two additional asymmetric dienophiles, chloro-maleic anhydride (2b) and cyano-maleic anhydride (2c in acetone and toluene). Similarly to 2d, these two anhydrides can also approach the diene in two different directions, resulting in the formation of isomeric products. Figure 3 shows that acetone does increase their selectivity.

In summary, we mention that, in many applications of the DA reaction, maximizing the yield of one isomeric product over the other is an essential goal. Changing the solvent can be a convenient and simple way to move toward this objective.

3.2. Geometries

The gas-phase optimised geometries of all the minima structures and transition states are given in Figure S4. Among two dienophiles (2a and 2d), 2a has a slightly higher bond distance of the C3=C4 double bond (see Scheme 1 for the numbering scheme). In the transition states, there is a little increase of about 0.06 Å in this bond length, and, finally, in the products, it further increases to 1.54 Å from an initial length of almost 1.3 Å. This increase in the bond length is consistent with the fact that a double bond is converted into a single bond. In TSa, as the dienophile approaches C9 and C10 (see Scheme 1) of the diene, C9 maintains a bit more distance from the dienophile as compared to C10 (see Figure S4). It indicates that the methyl group on the diene does not offer favourable interaction for this cycloaddition reaction. Although, the methyl group is electron-donating and its presence on the diene increases the reaction rate, in the present case, the steric hindrance seems to be a dominating factor [63]. Analysis of the bonds which are being formed in TSd and TSd’ reveals that the C3, which is directly attached to C2 has a greater ability to approach the diene than C4. This is because the acetoxy group has an electron-withdrawing nature and the presence of such a group on dienophiles makes the interactions of the diene and dienophile more favourable by increasing the global electron density transfer from the nucleophile to the electrophile [63]. The extent of synchronicity for the reactions can be determined via the difference in the bond lengths of the two new sigma bonds being formed in the transition state [64].

The geometries of the transition states shown in Figure S4 indicate that the Diels–Alder reactions under investigation are not completely synchronous. However, the degree of asynchronicity in TSa, TSd and TSd’ is not high enough to propose a stepwise mechanism (thus, these reactions can be expected to proceed via the commonly followed concerted mechanism). The asynchronicity of TSd is less than TSa by a value of 0.07 Å, which is contrary to the usual observation that asymmetric reactants increase the asynchronicity of the transition state [65] as 2d is clearly more asymmetric than 2a. The asynchronicity of TSb’ and TSc’ is higher than that of the TSb and TSc, respectively. The lower barriers for the former transition states may be attributed to their higher asynchronicity [66].

In the solvents (toluene and acetone), the extent of asynchronicity is increased for all the transition sates, except for the case of TSb. For TSa, it is increased from 0.12 Å (in the gas phase) to 0.13 Å and 0.14 Å in toluene and acetone. Likewise, for other transition states, the increase in asynchronicity in acetone is higher than in toluene. TSb is an exception, where it is decreased from 0.048 Å in the gas phase to 0.041 Å and 0.035 Å in toluene and acetone, respectively. Although, for TSb, it is decreased in both solvents, the overall change in acetone remains greater than in toluene.

3.3. Molecular Electrostatic Potential Analysis

Next, we were interested in rationalising the results via molecular electrostatic potential (MESP) analysis. MESP analysis provides an efficient way of analysing distinct regions in molecules. It helps in understanding the effect of structural modifications of molecules on their reactivities [67]. It displays the distribution of electron density and indicates the sites suitable for the attack of nucleophiles and electrophiles [58]. MESP analysis for the reactants and products has been carried out and Figure 4 displays these maps. The blue and red colours in these plots indicate the regions of highest and lowest electrostatic potential, respectively. Higher electrostatic potential corresponds to a lower density of electrons, while a more negative value of electrostatic potential indicates an abundance of electrons. Intermediate values of electrostatic potential are represented by light-blue, yellow and green colours, depicting marginally electron-deficient, slightly electron-rich and neutral sites of molecules, respectively [68]. In DA reactions, the electron density is transferred from the diene to the dienophile; therefore, reactant 1 (diene) should behave like a nucleophile while the reactants 2a and 2d (dienophiles) must act as an electrophile. An examination of Figure 4 clearly reveals that carbon atoms of the diene have a higher electron density than the carbon atoms of dienophiles. If the relative electrostatic potential of the two dienophiles (2a and 2d) is compared, then it can be noticed that 2a has a comparatively greater value around its CC double bond. It means that 2a is a better electrophile as compared to 2d and this is consistent with the fact that the activation free energy for the reaction of 1 and 2a is smaller than that of the reactions between 1 and 2d. The MESP of products indicates that the blue regions of dienophiles and red regions of the diene have turned green. It verifies that electron density has been transferred from the diene to the dienophiles and the molecules have become less polar overall.

3.4. Frontier Molecular Orbital (FMO) Analysis

To gain more insight into the electronic properties of the reacting molecules, the energies of their HOMO (highest occupied molecular orbital) and LUMO (lowest occupied molecular orbitals) were computed in the gas phase, toluene and acetone [69]. According to molecular orbital theory, the energy of the HOMO indicates the electron-releasing ability of molecules, while that of the LUMO illustrates the potential of molecules to attract electrons [70]. As shown in Figure 5, the molecule of the diene (1) has a higher value of HOMO energy than the dienophiles, leading to its higher electron-donating ability.

In toluene and acetone, the HOMO of the diene is slightly reduced, while the LUMO of both dienophiles is somewhat increased in terms of energy compared to the gas phase. Despite this individual change in energies, their relative energies remain in the same order in both solvents as in the gas phase (a figure for FMO interactions in toluene and acetone is provided in Supplementary Information, i.e., Figure S5). The LUMO energy is lower for 2a than for 2d facilitating the reaction, in agreement with the trend in activation energies. The difference in the energies of these frontier molecular orbitals (FMOs) is an indicator of the kinetic stability of molecules. A higher energy gap corresponds to greater kinetic stability, and vice versa for a lower gap [71]. The kinetic stability of TSd is higher than that of the TSd’ (as indicated in Table S2), and hence later has greater ease of transition [72]. Depending upon the relative FMO energies of the diene and dienophile, the Diels–Alder reaction may proceed via normal electron demand (NED) or inverse electron demand (IED) [73]. In Figure 4, the interactions of the frontier orbitals of the reactants are represented in the gas phase. It shows that the studied reactions follow the normal electron demand pathway as the difference between the HOMO (diene) and the LUMO (dienophile) is smaller than the energy gap between the HOMO (dienophile) and the LUMO (diene) [74]. In toluene and acetone, these reactions still have normal electron demand; however, the value of the IED and NED gaps is changed slightly. For 1+2a, the NED is increased from 4.2 eV in the gas phase to 4.3 eV and 4.5 eV in toluene and acetone, respectively, while, for 1+2d, it is same in the gas phase and toluene and increased by 0.1 eV in acetone. The IED for them is reduced from 8.8 eV in the gas phase to 8.7 eV in toluene and 8.6 eV in acetone.

3.5. Global Reactivity Indices (GRI)

The energies of frontier molecular orbitals can be used to determine global reactivity descriptors, including chemical hardness and softness, chemical potential and the nucleophilicity index and electrophilicity index. Equations (1)–(5), provided in the computational details section, are used for the calculations of GRI. Global hardness estimates the extent to which charge transfer is prohibited in a molecule; conversely, global softness measures the electron-accepting ability of a molecule. The electrophilicity index calculated from the chemical hardness and chemical potential helps in determining the strength of the electrophile; the molecule may be a strong electrophile ($\omega$ > 1.5 eV), marginal electrophile ($\omega$ < 0.8 eV) or moderate electrophile (1.5 eV > $\omega$ > 0.8 eV). A smaller value of $\omega$ indicates a poor electrophile but a better nucleophile and vice versa [75,76,77,78]. Figure 6 indicates that 1 is a poorer electrophile as compared to the dienophiles, and, among the two dienophiles, 2a clearly outperforms 2d. This figure also shows that the solvents decreased the electrophilic character of the dienophiles but raised it for the diene, and the effect of acetone is more pronounced than that of toluene. From Table S2, it can be inferred that, among the reacting molecules, the dienophiles are better electrophiles than the diene. It also indicates that the softness of the dienophiles is decreased in solvents, and, as the role of the dienophile is to accept the electrons in the DA reaction, a decrease in its softness causes the reaction to proceed at a lower rate. A comparison of chemical potential indicates that the charge will be transferred from electron-rich species 1 to electron-deficient species (2a and 2d) [63]. The nucleophilicity index of 1 shows that it is a strong nucleophile (N > 3 eV), while 2a and 2d (N < 2 eV) are very weak nucleophiles [79]. As the diene in the DA reaction behaves as a nucleophile and the dienophile as an electrophile, the higher N of 1 and larger $\omega$ of 2a and 2d are favourable factors for these reactions.

3.6. Natural Bond Orbital (NBO) Analysis and Fukui Functions

For determining the atomic charges on the relevant atoms of reactants and products, NBO analysis is performed [80]. In 2a, the NBO charges on C3 and C4 are −0.279e, while, in 2d, the charge on C3 is −0.197e and on C4 is −0.320e. The negative charge on C3 is decreased due to the presence of an electron-withdrawing acetoxy group on its neighbouring carbon C2. In the products, the negative charges on C3 and C4 are increased compared to the dienophiles because they received the negative charge from the diene [56]. Fukui functions are local reactivity indices, and they define the selectivity of a particular site or atom in a given molecule. These functions can be described as the first differentials of electron density with respect to the number of electrons in a chemical system, keeping the external potential constant. These are determined from single-point energy calculations on optimised structures of reactants using Equations (6) and (7). $f_w^{-}$ are computed for the diene and $f_w^{+}$ for the dienophiles. Figure 7 indicates that the C9 and C10 atoms of the diene have the highest values of $f_w^{-}$, leading to their greater susceptibility to the attack of the electrophile [81]. On the other hand, $f_w^{+}$ is highest for atoms C3 and C4 of the dienophiles, making them the most reactive sites for the attack of nucleophiles. From these values, it can be predicted that new bonds will form between the aforementioned atoms of the reactants, in line with the observed reaction mechanism.

4. Conclusions

In the present work, the cycloaddition reactions of 9-methylanthracene with maleic anhydride, chloro-maleic anhydride, cyano-maleic anhydride and (5-oxo-2H-furan-2-yl) acetate were investigated. These reactions were explored using density functional theory in the gas phase, toluene and acetone. It was found that these reactions are faster in the gas phase than in toluene and in acetone because these solvents increase the gap between the LUMO of the dienophiles and the HOMO of the diene. The barrier for the reaction of 2a is lower as compared to that for the reaction of 2d, which is an indication of the superior dienophilic properties of maleic anhydride over the (5-oxo-2H-furan-2-yl) acetate. The reaction of 2d with the diene formed two products (3d and 3d’), and the formation of 3d’ was favoured over 3d in the gas phase as well as in the solvents, which is in agreement with the experimental observation. Crucially, the retardation in the rate of reaction by acetone is lower for 3d’ than for 3d, revealing the possibility of increasing the selectivity of 3d’ by using polar solvents. The same effect is also observed for the reaction of 1 with 2b and 2c. Values of their global softness, chemical potential and electrophilicity index were used to illustrate the underlying electronic structure. Furthermore, the values of Fukui functions showed that the atoms of the diene and dienophiles involved in bond formation are the most suitable sites for the attack of electrophiles and nucleophiles, respectively, in line with the observed reaction mechanism. This work would be helpful for designing optimum conditions for the Diels–Alder reaction in the future taking into account the crucial influence of the solvent on selectivity.