Unveiling the reaction mechanism and the regioselectivity of Diels–Alder cycloadditions between 2-methoxy-3-thiophenylbutadiene and methyl vinyl ketone using topological tools

Abstract

The [4 + 2] Diels–Alder cycloaddition reaction between 2-methoxy-3-thiophenylbutadiene 1 and methyl vinyl ketone 2 has been studied at the density functional theory level using a panoply of tools to unravel the regioselectivity and the reaction mechanisms. From the analysis of the CDFT reactivity indices, 2 behaves as an electrophile, while 1 behaves as a nucleophile. This cycloaddition is characterized by a non-polar character as evaluated by a low value of the electrophilicity difference between the reactants as well as a low value of the global electron density transfer at the transition state. Futher, the bonding evolution theory shows that this reaction takes place via a one-step asynchronous mechanism. QTAIM descriptors also highlight the asynchronicity, which confirms the absence of the formation of a new covalent bond. Generally, the ortho/endo pathway is both thermodynamically and kinetically favored over the other routes as revealed by the Distortion/Interaction-Activation Strain (DIAS) analysis, which shows that this preference is primarily driven by differences in interaction energies.

---

1. Introduction

Diels–Alder (DA) cycloadditions are fundamental [4 + 2] reactions that efficiently construct six-membered rings from conjugated dienes and unsaturated dienophiles.^1,2 They provide rapid access to highly substituted cyclohexenes with multiple stereogenic centers, making them essential for complex molecule synthesis.^3,4 The use of heteroatom-substituted dienes further expands their scope by introducing functional groups that can be retained or transformed in subsequent steps.^5–7

Regioselectivity, largely determined by substituents on the diene and dienophile, is a key factor in these reactions. In unsymmetrical dienes, one substituent often acts as a regio-director, controlling the major product. For example, the cycloaddition of 2-methoxy-3-thiophenylbutadiene 1 with methyl vinyl ketone 2 predominantly yields 1-thiophenyl-2-methoxy-4-acetylcyclohex-1-ene (3, ∼80%) due to the directing effect of the SPh group (see Scheme 1).⁸

Scheme 1 DA cycloaddition between 2-methoxy-3-thiophenylbutadiene 1 with methyl vinyl ketone 2.

Building on this, the present study uses computational methods to investigate the electronic factors controlling regioselectivity, providing mechanistic insight into the formation of the major regioisomer.

Quantum chemical methods have become indispensable for gaining insight into DA reactions, as they allow a detailed examination of activation barriers,^9,10 reaction pathways,^11,12 and stereochemical outcomes.^13,14 Computational studies have been particularly effective in identifying reactivity trends, distinguishing between concerted and stepwise mechanisms, and rationalizing endo/exo pref. 15 and 16 Within this context, Frontier Molecular Orbital (FMO)¹⁷ considerations have traditionally been used to interpret cycloaddition reactivity by correlating activation energies with the energetic matching of interacting orbitals. A reduction in the HOMO–LUMO energy separation, achieved for instance by introducing electron-donating substituents on the diene or electron-withdrawing groups on the dienophile, generally leads to an enhancement of the reaction rate.¹⁸

Conceptual Density Functional Theory (CDFT) was developed to provide a more quantitative understanding of chemical reactivity by introducing global descriptors such as the electronic chemical potential (µ), electronegativity (χ), chemical hardness (η), softness (S), electrophilicity (ω), and nucleophilicity (N).^19–25 In particular, the electrophilicity and nucleophilicity indices have been widely employed to assess the polar character of cycloadditions reactions.^26–29 At a more localized level, Parr and Yang proposed molecular Fukui functions,^30,31 which allocate global reactivity indices to individual atoms within a molecule. These local descriptors allow the identification of the most electrophilic and nucleophilic sites in a reagent, offering a detailed picture of atomic reactivity and providing insights into the regioselectivity observed in DA cycloadditions.^32,33

In parallel, approaches centered on electron density analysis have provided complementary insight into the nature of chemical bonding and its evolution along reaction coordinates. The Quantum Theory of Atoms in Molecules (QTAIM)³⁴ and the Electron Localization Function (ELF)³⁵ allow changes in electron density organization to be monitored during the course of a reaction.^36,37 On this basis, the Bonding Evolution Theory (BET)³⁸ was developed to describe bond formation and cleavage in terms of topological changes in ELF basins.^39,40 More recently, Molecular Electron Density Theory (MEDT)⁴¹ has emphasized that variations in electron density, rather than orbital interactions alone, are the driving force behind organic reactivity. Within this framework, the combined use of CDFT descriptors, electron density topology, and non-covalent interactions (NCI)⁴² analyses offers a consistent strategy for understanding both reactivity and selectivity in DA cycloadditions.^43,44

This study provides a comprehensive mechanistic and selectivity analysis of the DA cycloaddition between 2-methoxy-3-thiophenylbutadiene 1 and methyl vinyl ketone 2 using MEDT framework. Owing to the non-symmetry of the reactants, the reaction can proceed along four competing pathways, corresponding to two regioisomeric routes (ortho and meta) and two stereoisomeric routes (endo and exo). The four associated transition states (TSs) (TS-on, TS-ox, TS-mn, and TS-mx) leading to the corresponding cycloadducts (CA-on, CA-ox, CA-mn, and CA-mx) have been located and fully characterized (Scheme 2). By applying BET, we reveal how bonds form and break, how electron density is redistributed, and whether the reaction occurs via a concerted yet asynchronous pathway. Complementary analyses using NCI, QTAIM, and ELF further identify the key electronic and structural factors governing the regio- and stereoselectivity observed experimentally.⁹ Importantly, this work provides novel insight into the role of substituent-dependent interaction energies in directing both kinetic and thermodynamic preferences, offering a unified topological-energetic framework that connects asynchronicity, selectivity, and reactivity in DA cycloaddition between 2-methoxy-3-thiophenylbutadiene 1 with methyl vinyl ketone 2 a contribution not addressed in previous studies.

Scheme 2 Possible regio- and stereoisomeric pathways for the DA cycloaddition between 2-methoxy-3-thiophenylbutadiene 1 with methyl vinyl ketone 2.

2. Computational methods

All quantum chemical calculations, including geometry optimizations and harmonic vibrational frequency analyses for reactants, TSs, and products, as well as intrinsic reaction coordinate (IRC) computations,¹⁷ were carried out within the framework of DFT using the Gaussian 16 software suite.⁴⁵ The M06-2X hybrid meta-GGA functional developed by the Minnesota group⁴⁶ was employed in conjunction with the 6-311 + G(d,p) basis set. This computational protocol has proven to be well suited and reliable for describing similar reaction mechanisms. The D3 correction was also considered for the role of dispersion correction energies.⁴⁷ Solvent [diethylether, ε = 4.33] effects were modelled with the IEFPCM (integral equation formalism of the polarizable continuum model) method.⁴⁸ Calculations were performed for T = 298.15 K and, initially, for P = 1.0 atm. The stationary points corresponding to reactants, products, and TSs were identified and validated through frequency calculations, ensuring the absence of imaginary frequencies for minima and the presence of a single imaginary frequency for each TS. To confirm the proper connection between reactants and products, IRC calculations were performed for all TSs using the standard IRC procedure.

NCI analyses were conducted based on the reduced density gradient and low-gradient isosurfaces,^42,49 employing the NCIPlot program.⁵⁰ In addition, QTAIM analysis was carried out using the Multiwfn package,⁵¹ relying on the mono-determinantal wavefunctions obtained at the M06-2X-D3/6-311 + G(d,p) level of theory.

For the BET investigation along each reaction pathway, topological analyses were performed by extracting the wavefunction at successive points along the IRC. The ELF was subsequently computed using the TopMod package,⁵² adopting a cubic grid with a spacing smaller than 0.2 bohr. The evolution of electron populations associated with atoms involved in bond formation along the reaction coordinate was visualized using the DrawProfile program,⁵³ whereas ELF isosurfaces were rendered with DrawMol.⁵⁴ Within the BET framework, two categories of ELF basins were identified: monosynaptic basins V(A), associated with localized lone electron pairs on atom A, and disynaptic basins V(A,B), corresponding to shared electron density between atoms A and B.

3. Results and discussion

3.1. Reactants properties

3.1.1. ELF and NPA analysis. The topological analysis of the ELF is a valuable method to study the distribution of electron density in a molecule, thereby facilitating the prediction of its chemical reactivity.^55,56 Thus, an ELF topological analysis of the electronic structures of 2-methoxy-3-thiophenylbutadiene 1 and methyl vinyl ketone 2 was first performed. ELF basin attractor positions, together with the most relevant valence basin populations, are shown in Fig. 1.

Fig. 1 ELF basin attractor positions of 2-methoxy-3-thiophenylbutadiene 1 and methyl vinyl ketone 2, in vacuo.

The ELF topology of 2-methoxy-3-thiophenylbutadiene 1 is defined by three disynaptic basins, V(C1, C2), V(C2, C3), and V(C3, C4), with electron populations of 3.64, 2.21, and 3.48e, respectively, revealing a markedly delocalized π system (Fig. 1). Notably, the significantly reduced population of the central V(C2, C3) basin indicates pronounced electronic depletion and polarization within the conjugated framework, which enhances the nucleophilic character of the terminal carbon atoms and facilitates their participation in the cycloaddition process.

In contrast, the ELF analysis of methyl vinyl ketone 2 exhibits a well-localized C5=C6 disynaptic basin with a population of 3.37e, characteristic of an electron-poor π bond (Fig. 1). This localization reflects the strong electron-withdrawing effect of the carbonyl group, reinforcing the electrophilic nature of the dienophile. The complementary ELF signatures of 1 and 2 therefore support a polar DA reaction, in which electron density flows from the delocalized π system of the diene toward the activated C=C bond of the dienophile, consistent with the observed reactivity trends.

Scheme 3 illustrates the Lewis-type representations derived from the ELF valence basin populations, together with the corresponding Natural Population Analysis (NPA) charges.⁵⁷ The NPA results reveal a marked charge separation within compound 1, with the C1, C3, and C4 carbon atoms bearing negative charges of −0.50, −0.20, and −0.31e, respectively, whereas the C2 carbon atom carries a positive charge of +0.30e. In comparison, methyl vinyl ketone 2 displays a more homogeneous charge distribution along the C5=C6 fragment, with both carbon atoms exhibiting identical negative charges of −0.30e.

Scheme 3 Proposed Lewis-like structures together with the natural atomic charges in average number of electrons e of 2-methoxy-3-thiophenylbutadiene 1 and methyl vinyl ketone 2. Negative and positive charges are shown in red, and blue colors, respectively, in vacuo.

3.1.2. CDFT analysis. CDFT provides a reliable approach to analyse reactivity in DA reactions.^58,59 The global reactivity indices, namely, the electronic chemical potential µ, chemical hardness η, global electrophilicity ω, and global nucleophilicity N of 2-methoxy-3-thiophenylbutadiene 1 and methyl vinyl ketone 2 are given in Table 1. The chemical potential µ of 2-methoxy-3-thiophenylbutadiene 1 (−3.72 eV) is higher than that of methyl vinyl ketone 2 (−4.79 eV), indicating that the GEDT is transferred from 1 to 2, consistent with a forward electron density flux (FEDF) mechanism.⁶⁰ The chemical hardness η of 2 (8.42 eV vs. 7.35 eV for 1) reflects greater resistance to electronic deformation and a more rigid π-system, in line with its electrophilic character. This is further supported by the global electrophilicity index ω, which is higher for 2 (1.36 eV) than for 1 (0.94 eV), whereas the nucleophilicity index N of 1 (2.09 eV) surpasses that of 2 (0.49 eV), highlighting the moderate nucleophilic character of the diene. Taken together, these descriptors classify 1 as a moderate nucleophile and weak electrophile, while 2 behaves as a moderate electrophile and a weak nucleophile, consistent with standard electrophilicity and nucleophilicity scales.⁶¹ Moreover, the studied DA reaction has a non-polar character, as substantiated by the low value of the difference of electrophilicity indices between the two reactants (Δω = 0.4 eV).⁶²

Table 1 DFT-based reactivity indices (in eV) of reactants, in vacuo

	µ	η	ω	N
1	−3.7	7.4	0.9	2.1
2	−4.8	8.4	1.4	0.5

---

In polar cycloaddition reactions that involve non-symmetrical reactants, regioselectivity is generally governed by the interaction between the most electron-rich site of the nucleophile and the most electron-deficient site of the electrophile. In this context, Domingo introduced the nucleophilic P⁺_k and electrophilic P⁺_k Parr functions,⁶³ which are obtained from electron-density variations associated with the GEDT from the nucleophile to the electrophile. These functions are currently regarded as one of the most reliable descriptors for evaluating local reactivity in polar and ionic reaction mechanisms.

Accordingly, the nucleophilic P⁻_k Parr functions of 1 and the electrophilic P⁺_k Parr functions of 2 were examined to identify their most reactive atomic centers (Fig. 2 and Table S3).

Fig. 2 3D representations of the spin densities (η = 0.02) of 1 and 2 together with the nucleophilic P⁻_k Parr functions for 1, and the electrophilic P⁺_k Parr functions for 2, in vacuo.

The analysis reveals that, within 1, the C4 carbon atom (see Scheme 2 for atom numbering) exhibits the highest nucleophilic character, with a P⁻_k value of 0.08. Conversely, in 2, the P⁺_k Parr function reaches its maximum at the C5 carbon atom, which displays a markedly high electrophilic value of 0.50. These findings indicate that the preferred interaction occurs between C4 of the nucleophile and C5 of the electrophile, promoting a regioselective C4–C5 bond formation. This interaction pathway leads to the formation of the ortho regioisomer, in agreement with the experimentally reported predominance of the corresponding ortho cycloadduct (∼80%).⁸

3.2. Thermodynamic analysis

Table 2 indicates that the investigated Diels–Alder reaction is thermodynamically favorable in the gas phase, as evidenced by negative enthalpy and Gibbs free energy changes for all cycloadducts. The calculated free energy variations span from −24.2 to −31.1 kcal mol⁻¹, confirming that cycloadduct formation is both exothermic and exergonic. Among the four possible products, CA-mn corresponds to the global thermodynamic minimum, exhibiting the largest free energy release (ΔG = −31.1 kcal mol⁻¹), while CA-on is the second most stable cycloadduct (ΔG = −26.1 kcal mol⁻¹).

Table 2 Thermochemical state functions (ΔE, ΔH, ΔG, in kcal mol⁻¹ and ΔS in cal mol⁻¹ K⁻¹) of the DA reaction between 1 and 2. Solvent (diethylether) effects were described using IEFPCM

	ΔE	ΔH	ΔS	ΔG
In vacuo
TS-on	8.8	9.9	−51.0	25.1
TS-ox	12.5	13.4	−46.8	27.3
TS-mn	9.6	10.9	−51.9	26.3
TS-mx	13.2	14.1	−46.9	28.1
CA-on	−44.5	−40.6	−48.4	−26.1
CA-ox	−42.2	−38.3	−47.1	−24.2
CA-mn	−48.8	−44.2	−46.6	−31.1
CA-mx	−42.2	−38.3	−46.4	−24.5
Diethylether (ε = 4.33)
TS-on	9.1	10.3	−50.9	25.4
TS-ox	12.4	13.4	−46.7	27.3
TS-mn	9.8	10.9	−50.9	26.2
TS-mx	12.9	13.7	−46.8	27.6
CA-on	−43.9	−39.6	−48.6	−25.2
CA-ox	−41.7	−37.4	−46.9	−23.5
CA-mn	−42.8	−38.2	−45.0	−25.3
CA-mx	−41.6	−37.2	−44.9	−24.3

---

The remaining products, CA-mx and CA-ox, are less stabilized but remain thermodynamically accessible (Fig. 3). The analysis of the activation Gibbs free energies highlights a clear differentiation between the regio- and stereochemical pathways. Along the ortho approach, the TS-on is kinetically favored over its exo counterpart TS-ox by 2.2 kcal mol⁻¹. A similar trend is observed for the meta pathway, where TS-mn exhibits a lower activation barrier than TS-mx. These results consistently point to a pronounced preference for endo selectivity, regardless of the regiochemical mode of approach. A comparison of all TSs further reveals that TS-on possesses the lowest activation free energy among the accessible pathways. Consequently, despite the theoretical accessibility of multiple reaction channels, the reaction is kinetically dominated by the TS-on pathway, which accounts for more than 85% of the TS Maxwell–Boltzmann populations.⁶⁴ This strong kinetic bias dictates the initial regioselectivity of the reaction.

Fig. 3 Relative Gibbs free energy diagram (298.15 K) of the DA reaction between 1 and 2 in vacuo.

In contrast, the relative stability of the cycloadducts follows a different ordering, with CA-mn being overwhelmingly favored on thermodynamic grounds. This apparent divergence between kinetic and thermodynamic preferences indicates that the reaction operates under kinetic control at the TS level, followed by thermodynamic stabilization of the products. Such a scenario consistently explains the experimentally observed product distribution, where CA-on is formed as the major product (80%), while CA-mn appears as a significant minor component (20%).⁹ To further assess the interplay between kinetics and thermodynamics, a linear correlation between the calculated activation energies (ΔE) and reaction energies (ΔE°) was examined within the framework of the Bell–Evans–Polanyi (BEP) principle.^65,66 The resulting relationship, ΔE^≠ = 0.5 Δ_rE + 33.24 (R² = 0.52) (Fig. S1), displays a slope close to the classical BEP value of ∼0.5,⁶⁷ indicating a moderate thermodynamic contribution to the activation barriers. However, the relatively modest correlation coefficient suggests that barrier heights are not governed exclusively by reaction exergonicity but are significantly influenced by additional electronic and structural factors specific to each competing pathway. This interpretation is further supported by the DIAS analysis (see Section 3.6), which demonstrates that differences in activation barriers primarily arise from pathway-dependent interaction energies rather than from distortion effects alone. Therefore, while a general BEP-type tendency can be identified, the noticeable deviations arise from the distinct electronic character and structural organization of the competing TSs. This ultimately results in kinetic selectivity that does not directly mirror the thermodynamic stability of the final products.

The inclusion of diethyl ether slightly reduces the exergonicity of cycloadduct formation (e.g., ΔG = −25.3 kcal mol⁻¹ for CA-mn and −25.2 kcal mol⁻¹ for CA-on) without significantly affecting the activation barriers (ΔG ≈ 25.4–27.7 kcal mol⁻¹), pointing to a modest solvent effect that preserves the regio- and stereoselectivity trends. Overall, these results demonstrate that the outcome of the reaction is predominantly controlled by kinetic factors, with product distributions determined mainly by the relative heights of the activation barriers rather than by the thermodynamic stability of the cycloadducts.

Methyl vinyl ketone exhibits two distinct geometric conformers, namely s-cis and s-trans. DFT calculations show that the s-trans conformer is slightly more stable in terms of electronic energy and enthalpy than the s-cis one (Table S1). However, when Gibbs free energies are considered, the s-cis conformer becomes marginally favored (ΔG = −0.17 kcal mol⁻¹), indicating a favorable entropic contribution. As a result, the Maxwell–Boltzmann⁶⁷ populations at 298 K predict a predominance of the s-cis conformer in solution (P ≈ 57%), while the s-trans conformer accounts for about 43% of the population. Owing to this small free-energy difference, both conformers are significantly populated, and their relative reactivities are therefore worth considering. Unless otherwise stated, the following Sections and Subsections implicitly refer to the s-trans conformation.

The TSs geometries are depicted in Fig. 4 together with the length of new forming C1–C6 and C4–C5 bonds for the ortho pathway and C1–C5 and C4–C6 for the meta pathway (see Scheme 2 for the nomenclature). To further elucidate the asynchronous character of the studied DA reaction, the interatomic distances (r), bond development indices (l), and the asymmetry parameter (Δl) were examined for both the ortho and meta pathways (Table 3).⁶⁸ In the ortho approach, TS-on displays l values of 0.671 and 0.455 for the for the C4–C5 and C1–C6 bonds, respectively, resulting in a pronounced asymmetry (Δl = 0.216). Similarly, TS-ox shows slightly lower but still significant asynchronicity, with l values of 0.662 (C4–C5) and 0.449 (C1–C6), and Δl = 0.168. In both TSs, the formation of the C4–C5 bond is clearly more advanced than that of the C1–C6 bond, confirming a strongly asynchronous bond-forming process along the ortho pathway. In contrast, For TS-mn, the forming bond distances are 2.11 Å (C1–C5) and 2.35 Å (C4–C6), corresponding to l values of 0.639 and 0.494, with Δl = 0.145. For TS-mx, the respective distances are 2.13 Å (C1–C5) and 2.34 Å (C4–C6), with l values of 0.626 and 0.471, giving Δl = 0.155. In both TSs, the C1–C5 bond formation is more advanced than that of the C4–C6 bond, although the smaller Δl values indicate a significantly lower degree of asynchronicity compared to the ortho pathway. They indicate that the studied DA cycloaddition take place via a one-step asynchronous mechanism, with the ortho pathway being more asynchronous than the meta one.

Fig. 4 Geometries of the TSs of the DA reactions between 1 and 2 in vacuo.

Table 3 M06-2X-D3/6-311 + G(d,p) key parameters of critical structures for the studied DA reaction, in vacuo

	rC4–C5 (Å)	lC4–C5	rC1–C6 (Å)	lC1–C6	Δl
TS-on	2.06	0.671	2.41	0.455	0.216
CA-on	1.55		1.56
TS-ox	2.06	0.662	2.42	0.449	0.168
CA-ox	1.54		1.56

---

	rC1–C5 (Å)	lC1–C5	rC4–C6 (Å)	lC1–C6	Δl
TS-mn	2.11	0.639	2.35	0.494	0.145
CA-mn	1.55		1.56
TS-mx	2.13	0.626	2.34	0.471	0.155
CA-mx	1.55		1.53

---

The amount of GEDT at the TSs was estimated as the sum of the natural atomic charges of the atoms belonging to the 2 moieties. GEDT values above 0.20e are characteristic of polar reactions, while values below 0.05e indicate non-polar reactions.⁶⁹ The computed GEDT values at the TSs range between 0.02 and 0.08e (Fig. 4), indicating a non-polar character of the studied DA reaction, which is in good agreement with the CDFT analysis.

3.3. Topological analysis of the ELF at the TSs

The Fig. 5 illustrates the ELF localization domains for the four competitive TSs, TS-on, TS-ox, TS-mn, and TS-mx, involved in the studied DA reaction.

Fig. 5 ELF localization domains and the positions of TS attractor basins (η = 0.80). Monosynaptic basins are shown in red, protonated basins in blue, disynaptic basins in green, and attractor positions in purple, in vacuo.

The ELFs of all TSs display the presence of a three disynaptic basins, V(C1, C2), V(C2, C3), and V(C3, C4), associated with the diene framework, and a single disynaptic basin, V(C5, C6), at the dienophile framework. V(C1, C2), V(C3, C4), and V(C5, C6) exhibit integrated populations of 3.30–3.34e, 3.21–3.25e, and 2.94–3.14e, corresponding to the C1–C2, C3–C4, and C5–C6 double bonds, respectively, whereas V(C2, C3) has an integrated population of 2.65–2.67e, consistent with a single C2–C3 bond. Notably, the ELF of TS-on and TS-ox displays a one monosynaptic basin, V(C4), with populations of 0.19 and 0.20e, respectively. The absence of the disynaptic basins, V(C1, C6) and V(C4, C5), indicates that the formation of these new C–C bonds has not yet started at these TSs.

3.4. QTAIM analysis of the TSs

The QTAIM analysis provides insight into the nature of the interactions at the TSs through the electron density ρ(r) and its Laplacian ∇²ρ(r) at the bond critical points (BCPs) (Fig. 6 and Table 4). For TS-on, TS-ox, TS-mn, and TS-mx, the ρ(r) values at all forming C–C BCPs remain below 0.1 a.u., and the corresponding ∇²ρ(r) values are positive, which are characteristic features of closed-shell interactions rather than fully developed covalent bonds. These results indicate that the new C–C bonds are not yet formed and that bond formation is still at an early stage at these TS structures, in agreement with the analysis of their geometric structures and ELF distributions.

Fig. 6 AIM representation of TSs. The violet ellipses highlight the bond critical points of new forming bonds, in vacuo.

Table 4 M06-2X-D3/6-311 + G(d,p) QTAIM parameters (in a.u.) of the bond critical points presented in the four competitive TSs associated with the studied DA reaction, in vacuo

TSs	BCPs
	C1–C6		C4–C5		C1–C5		C4–C6
	ρ	∇^2ρ(r)	ρ	∇^2ρ(r)	ρ	∇^2ρ(r)	ρ	∇^2ρ(r)
TS-on	0.0665	0.0371	0.0435	0.0499	—	—	—	—
TS-ox	0.0725	0.0321	0.0389	0.0483	—	—	—	—
TS-mn	—	—	—	—	0.0647	0.0385	0.0445	0.0486
TS-mx	—	—	—	—	0.0723	0.0321	0.0392	0.0500

---

3.5. NCI analysis of the TSs

The NCI analysis of the TSs reveals distinct patterns of weak interactions (Fig. 7). TS-on is characterized by predominantly green surfaces, indicating weak van der Waals interactions, with small blue patches corresponding to minor attractive contacts (e.g., C–H⋯O) and scattered red regions reflecting localized steric repulsion, suggesting moderate stabilization but less favorable than TS-mn and TS-mx. TS-ox shows a similar pattern with slightly enhanced green regions, indicating additional van der Waals stabilization. In contrast, TS-mn and TS-mx display more extensive blue regions, indicating stronger attractive interactions (e.g., C–H⋯O and π⋯π contacts), along with reduced red regions and moderate green areas, highlighting lower steric hindrance and overall greater stabilization, with TS-mx being the most stabilized.

Fig. 7 NCI isosurfaces (0.4 isovalue) associated with the density overlap from −0.02 < sign(λ₂) ρ(r) < 0.02 a.u. at the TSs involved in the studied DA reaction, in vacuo.

3.6. Distorsion/interaction-activation strain (DIAS) analysis

Considerable progress over the past two decades has provided a clear mechanistic understanding of asynchronicity in DA reactions. Within the activation strain model, this feature arises from the interplay between strain and interaction energies, with reduced Pauli repulsion identified as the dominant factor.^70–72 To quantitatively assess the energetic factors governing reactivity and selectivity, the reaction profiles were analysed using the Distortion/Interaction-Activation Strain (DIAS) framework.^73,74 This approach decomposes the total electronic energy (ΔE = ΔE_Dis + ΔE_Int) into distortion (ΔE_Dis) and interaction (ΔE_Int) components, thereby providing direct insight into the physical origin of the activation barriers. The calculations were performed with the autoDIAS program.⁷⁵

In this context, DIAS analyses were carried out for the ortho (Fig. 8) and meta (Fig. 9) approaches, considering both endo and exo pathways. For the ortho pathway (Fig. 8), comparison of the ΔE_Int and ΔE_Dis curves along the forming C1⋯C4 bond reveals that the TS-on pathway exhibits a significantly more stabilizing interaction energy (ΔE_Int) than TS-ox, while the distortion energy (ΔE_Dis) profiles remain very similar for both TSs. When both contributions are considered, the activation energy is lower for TS-on than for TS-ox, and the corresponding TS occurs earlier along the reaction coordinate, in agreement with Hammond's postulate.

Fig. 8 Comparative DIAS diagrams for the endo (solid lines) and exo (dashed lines) approaches of the ortho pathway of the DA reaction between 2-methoxy-3-thiophenylbutadiene 1 and methyl vinyl ketone 2 along the reaction coordinate projected onto the forming C4⋯C5 bond distance as determined with the M06-2X-D3/6-311 + G(d,p) level. The TSs are indicated by a dot.

Fig. 9 Comparative DIAS diagrams for the endo (solid lines) and exo (dashed lines) approaches of the meta pathway of the DA reaction between 2-methoxy-3-thiophenylbutadiene 1 and methyl vinyl ketone 2 along the reaction coordinate projected onto the forming C1⋯C5 bond distance as determined with the M06-2X-D3/6-311 + G(d,p) level. The TSs are indicated by a dot.

For the meta pathway (Fig. 9), along the C1⋯C5 forming bond, ΔE_Int is even more stabilizing for the TS-mn pathway compared to TS-mx. The ΔE_Dis curves remain very close throughout the reaction coordinate, particularly beyond a C1⋯C5 distance of 2.2 Å, where only a slight elevation is observed for TS-mx. Consequently, the activation barrier is lower for TS-mn than for TS-mx, and its TS is also reached earlier, again consistent with Hammond's postulate.

Notably, all four competitive pathways display comparable distortion energies (ΔE_Dis), which mainly originate from deformation of the 2-methoxy-3-thiophenylbutadiene 1 fragment. Therefore, regioselectivity and endo/exo selectivity are primarily governed by differences in interaction energies rather than by distortion effects. The endo selectivity arises from a slightly more stabilizing interaction energy, despite a marginal distortion penalty relative to the exo pathway.

For the meta-approach, the higher ΔE_Dis value observed for the endo pathway originates from the closer spatial proximity between the carbon atom of the come group in 2 and the SPh group in 1, with d(C⋯SPh) distances of 3.81 Å (endo) and 5.01 Å (exo). Similarly, for the ortho approach, the slightly higher ΔE_Dis value for the endo pathway is attributed to closer steric contact between these groups, with d(C⋯SPh) distances of 4.72 Å (endo) and 6.55 Å (exo).

Overall, the DIAS results demonstrate that deviations from a simple BEP relationship (Section 3.2) originate from pathway-dependent interaction energies rather than from distortion effects alone. The SPh substituent plays a key directing electronic role by enhancing stabilizing interactions at the TS level, thereby governing regioselectivity, kinetic preference, and the degree of asynchronicity within a consistent energetic-topological framework. These findings indicate that the regioselectivity of the reaction is primarily dictated by the directing electronic effect of the SPh group, which enhances stabilizing interaction energies at the TS level.

3.7. BET analysis for the 1 + 2 DA reaction along the ortho- and meta- path

3.7.1. BET analysis along the TS-on of the studied DA reaction. The Fig. 10 presents the changes in electron populations of the key basins involved in the formation of CA-on, while the corresponding ELF isosurfaces are shown in Fig. 11.

Fig. 10 Electron populations (e) evolution for selected basins along the IRC for the TS-on path of the 1 + 2 DA reaction in vacuo. The relative potential energy along the IRC corresponds to the black dotted curve.

Fig. 11 ELF (η = 0.72) snapshots of selected basins related to the bond forming process along the IRC for the TS-on path of the 1 + 2 DA reaction, in vacuo.

The joint examination of these figures reveals that the reaction coordinate can be partitioned into seven successive structural stability domains (SSDs), each reflecting a specific topological reorganization associated with the stepwise creation of the two σ bonds C4–C5 and C1–C6 (atom numbering given in Scheme 2). In SSD-I, the ELF topology corresponds to that of the isolated reactants. For 1, three disynaptic basins, V(C1, C2), V(C2, C3), and V(C3, C4) are observed, with electron populations of 3.62, 2.28, and 3.48e, respectively. 2 exhibits a single disynaptic basin related to the C5=C6 double bond, V(C5, C6), integrating 3.35e. As SSD-I progresses, slight reductions in electron density occur in the V(C1, C2), V(C3, C4), and V(C5, C6) basins, with losses of 0.29, 0.22, and 0.20e, respectively, indicating a gradual attenuation of π bonding in favor of σ character. Simultaneously, the V(C2, C3) basin becomes more populated, reaching 2.59e, which signals the incipient development of double-bond character in the C2–C3 linkage.

The passage from SSD-I to SSD-II and SSD-III is marked by the appearance of two monosynaptic basins, V(C5) and V(C4), localized on the C5 and C4 atoms and containing 0.20 and 0.37e, respectively. These basins correspond to pseudoradical centers generated mainly through electron density withdrawal from the V(C3, C4) and V(C5, C6) bonding regions, which experience comparable depopulation. In the subsequent SSD-IV domain, a further pseudoradical center emerges on atom C6 of 2, characterized by a population of 0.16e. This feature results from the continued depletion of the C5–C6 bonding basin, amounting to an overall loss of approximately 0.72e. At the beginning of SSD-V, the pseudoradical centers located on C4 and C5 combine to form a new disynaptic basin, V(C4, C5). This constitutes the first significant topological event, identified as a cusp-type catastrophe, and corresponds to the formation of the C4–C5 σ bond. At this point, the newly created basin holds 1.03e, while the electron populations of the V(C1, C2), V(C3, C4), and V(C5, C6) basins continue to decrease. Further electronic rearrangements take place in the SSD-VI region, where an additional monosynaptic basin, V(C1), appears. This basin identifies C1 as a new pseudoradical center and represents a necessary precursor for the establishment of the second σ bond. Ultimately, at the onset of SSD-VII, the monosynaptic basins V(C1) and V(C6) merge, giving rise to the disynaptic basin V(C1, C6), which signals the formation of the C1–C6 bond. During SSD-VII, the population of this basin increases steadily from 1.13 to 1.76e, confirming the progressive stabilization of the newly formed σ interaction.

3.7.2. BET analysis along the TS-mn of the studied DA reaction. The Fig. 12 illustrates the variation of electron populations within the principal basins involved in the formation of CA-mn, whereas the corresponding ELF isosurfaces are shown in Fig. 13.

Fig. 12 Electron populations (e) evolution for selected basins along the IRC for the TS-mn path of the 1 + 2 DA reaction in vacuo. The relative potential energy along the IRC corresponds to the black dotted curve.

Fig. 13 ELF (η = 0.72) snapshots of selected basins related to the bond forming process along the IRC for the TS-mn path of the 1 + 2 DA reaction, in vacuo.

The combined analysis of these figures indicates that the reaction path can be segmented into five successive SSDs, which account for the stepwise topological changes associated with the formation of the two new σ bonds, C4–C6 and C1–C5 (atom labelling according to Scheme 2). SSD-I corresponds to the initial electronic structures of the separated reactants. The passage to SSD-II is marked by the emergence of three monosynaptic basins localized on C1, C5, and C6 atoms, with initial electron populations of 0.39, 0.23, and 0.16e, respectively. These populations arise from the depletion of the V(C1, C2) and V(C5, C6) disynaptic basins, which lose 0.54 and 0.20e, respectively.

In SSD-III, the pseudoradical centers on atoms C1 and C5 merge, leading to the appearance of the V(C1, C5) disynaptic basin. This event represents the first major topological transformation, classified as a cusp catastrophe, and corresponds to the formation of the first C1–C5 σ bond. At the onset of SSD-III, this basin holds an electron population of 0.99e, while the electron densities of the V(C1, C2), V(C3, C4), and V(C5, C6) basins continue to decrease. At the beginning of SSD-IV, further depletion of the V(C3, C4) basin induces the formation of an additional pseudoradical center on atom C4, with a population of 0.27e.

The final topological rearrangement occurs in SSD-V, characterized by the emergence of the V(C4, C6) disynaptic basin with an electron population of 1.01e, signalling the formation of the second σ bond. This population is nearly equal to the combined electron densities of the former monosynaptic basins V(C4) and V(C6) at the conclusion of SSD-IV. Additionally, throughout the progression from SSD-II to SSD-V, a continuous increase in the population of the V(C2, C3) basin is observed, reaching a total of 3.44e, which is consistent with the conversion of a single into a double bond.

4. Conclusion

In summary, the DA reaction cycloaddition between 2-methoxy-3-thiophenylbutadiene 1 with methyl vinyl ketone 2 has been investigated at the DFT level, using a broad range of descriptors and analysis schemes, leading to consistent and comprehensive results. Adopting the tools of CDFT, 1 is classified as a moderate nucleophile and weak electrophile, while 2 behaves as a moderate electrophile and a weak nucleophile. Moreover, the studied DA reaction has a non-polar character, as substantiated by the low value of the difference of electrophilicity indices between the two reactants (Δω = 0.4 eV) as well as by a small GEDT at the TS (0.02 < GEDT < 0.08e). The calculation of Parr functions reveals that the ortho regioisomeric pathway is more favorable than the meta one. Energy calculations in both the gas phase and diethyl ether indicate that the reaction proceeds under kinetic control, affording exclusive ortho regioselectivity and pronounced endo stereoselectivity. The solvent increases the activation barriers and attenuates the overall thermodynamic driving force by stabilizing the reactants more effectively than the TSs and products.

To rationalize these trends, the DIAS analysis was performed. The energetic decomposition reveals that the observed selectivity is governed primarily by differences in interaction energies among competing pathways, whereas distortion energies remain essentially comparable. The favored TS-on channel is characterized by more stabilizing TS interactions, which account for their lower activation barriers. Conversely, strain contributions do not provide significant discrimination between the pathways.

This interaction-controlled mechanism explains the departure from a simple Bell–Evans–Polanyi relationship, as the activation energies are modulated by pathway-dependent interaction terms rather than by structural deformation effects. In this framework, the SPh substituent plays a decisive electronic role by reshaping the interaction-energy profile at the TS level, thereby determining regioselectivity, stereoselectivity, and kinetic preference. The agreement between the ELF, QTAIM, and DIAS analyses establishes a coherent mechanistic framework, demonstrating that the reaction proceeds through a one-step asynchronous pathway.

Author contributions

Mohamed Chellegui, Sabir A. Mohammed Salih, and Lakhdar: writing, investigation, validation, methodology and Haydar A. Mohammad-Salimh, and Ali Ben Ahmed: editing, reviewing, and supervision. All authors have read and approved the final version of the manuscript.

Conflicts of interest

There are no conflicts to declare.

Data availability

The data supporting this article have been included as part of the supplementary information (SI). Supplementary information is available. See DOI: https://doi.org/10.1039/d6ra00118a.

Acknowledgements

M. C. is grateful to UNamur for financially contributing to his research stay.

References

1. O. Diels and K. Alder, Synthesen in der hydroaromatischen Reihe, Justus Liebigs Ann. Chem., 1928, 460, 98–122.
2. O. Diels and K. Alder, Synthesen in der hydro-aromatischen Reihe, II. Mitteilung: Über Cantharidin, Ber. Dtsch. Chem. Ges., 1929, 62, 554–562.
3. T. Lazar, Modern Methods of Organic Synthesis, Fourth Edition, Synthesis, 2006, 4269.
4. Cycloaddition Reactions in Organic Synthesis, ed. F. S. J Ĺ Baldwin, S. FR and Ń. D. Magnus, Maxwell Macmillan Pergamon Publishing Corporation, 2002.
5. P. J. Proteau and P. B. Hopkins, Sulfur-directed Diels-Alder reactions. Synthesis of 1,5-disubstituted cyclohexene derivatives, J. Org. Chem., 1985, 50, 141–143.
6. O. de Lucchi and L. Pasquato, The role of sulfur functionalities in activating and directing olefins in cycloaddition reactions, Tetrahedron, 1988, 44, 6755–6794.
7. T. H. Chan and C. V. C. Prasad, Chemistry of 1-(trimethylsiloxy)-1-methoxy-3-[alkyl(aryl)thio]-1,3-butadienes and the synthesis of aryl sulfides via a cycloaromatization reaction, J. Org. Chem., 1986, 51, 3012–3016.
8. B. M. Trost, W. C. Vladuchick and A. J. Bridges, Sulfur as a regiochemical control element. Synthesis of 2-alkoxy(acyloxy)-3-alkyl(aryl)thiobuta-1,3-dienes, J. Am. Chem. Soc., 1980, 102, 3548–3554.
9. J. S. Fell, B. N. Martin and K. N. Houk, Origins of the Unfavorable Activation and Reaction Energies of 1-Azadiene Heterocycles Compared to 2-Azadiene Heterocycles in Diels–Alder Reactions, J. Org. Chem., 2017, 82, 1912–1919.
10. M. Chellegui, B. Champagne and M. Trabelsi, Lewis acid-catalyzed Diels–Alder cycloaddition of 2,5-dimethylfuran and ethylene: a density functional theory investigation, Theor. Chem. Acc., 2022, 141, 21.
11. E. Cherni, A. I. Adjieufack, B. Champagne, M. Abderrabba, S. Ayadi and V. Liégeois, Density Functional Theory Investigation of the Binding of ThioTEPA to Purine Bases: Thermodynamics and Bond Evolution Theory Analysis, J. Phys. Chem. A, 2020, 124, 4068–4080.
12. F. Teixeira, J. E. Rodríguez-Borges, A. Melo and M. N. D. S. Cordeiro, Stereoselectivity of the aza-Diels–Alder reaction between cyclopentadiene and protonated phenylethylimine derived from glyoxylates. A density functional theory study, Chem. Phys. Lett., 2009, 477, 60–64.
13. T. A. Hamlin, F. M. Bickelhaupt and I. Fernández, The Pauli Repulsion-Lowering Concept in Catalysis, Acc. Chem. Res., 2021, 54, 1972–1981.
14. S. Portela and I. Fernández, Origin of Catalysis and Selectivity in Lewis Acid-Promoted Diels–Alder Reactions Involving Vinylazaarenes as Dienophiles, J. Org. Chem., 2022, 87, 9307–9315.
15. M. Chellegui, M. Trabelsi, B. Champagne and V. Liégeois, DFT Investigation of the Stereoselectivity of the Lewis-Acid-Catalyzed Diels–Alder Reaction between 2,5-Dimethylfuran and Acrolein, ACS Omega, 2025, 10, 833–847.
16. M. Chellegui, A. I. Adjieufack, M. Trabelsi, V. Liégeois and B. Champagne, Unveiling the Reaction Mechanism of Diels-Alder Cycloadditions between 2,5-Dimethylfuran and Ethylene Derivatives Using Topological Tools, ChemPhysChem, 2025, 26(6), e202400896.
17. K. Fukui, The path of chemical reactions – the IRC approach, Acc. Chem. Res., 1981, 14, 363–368.
18. K. Fukui, Formulation of the reaction coordinate, J. Phys. Chem., 1970, 74, 4161–4163.
19. R. G. Parr and W. Yang, Density-Functional Theory of the Electronic Structure of Molecules, Annu. Rev. Phys. Chem., 1995, 46, 701–728.
20. W. Kohn, A. D. Becke and R. G. Parr, Density Functional Theory of Electronic Structure, J. Phys. Chem., 1996, 100, 12974–12980.
21. P. W. Ayers and R. G. Parr, Variational Principles for Describing Chemical Reactions: The Fukui Function and Chemical Hardness Revisited, J. Am. Chem. Soc., 2000, 122, 2010–2018.
22. F. De Proft and P. Geerlings, Conceptual and Computational DFT in the Study of Aromaticity, Chem. Rev., 2001, 101, 1451–1464.
23. P. Geerlings, F. De Proft and W. Langenaeker, Conceptual Density Functional Theory, Chem. Rev., 2003, 103, 1793–1874.
24. R. G. Parr, R. A. Donnelly, M. Levy and W. E. Palke, Electronegativity: The density functional viewpoint, J. Chem. Phys., 1978, 68, 3801–3807.
25. W. Yang and R. G. Parr, Hardness, softness, and the fukui function in the electronic theory of metals and catalysis, Proc. Natl. Acad. Sci., 1985, 82, 6723–6726.
26. M. Chellegui, S. Benmetir, R. N. Salih, H. A. Mohammad-Salim, J. V. de Julián-Ortiz and A. Ben Ahmed, A molecular electron density theory investigation of the mechanism of intramolecular [3+2] cycloaddition (32CA) with the participation of nitrile N -oxide and ethene molecular segments, New J. Chem., 2025, 49, 7302–7313.
27. M. Chellegui, I. Salhi, A. Ben Ahmed, S. Benmetir, R. N. Salih, H. A. Mohammad-Salim and J. V. de Julián-Ortiz, Mechanistic study of the [2 + 2] cycloaddition of ethylene with ketene derivatives via MEDT, Struct. Chem., 2026, 37, 111–124.
28. S. Benmetir, M. Chellegui, L. Benhamed, R. Mowafak Al-Mokhtar, R. N. Salih, M. A. Algso, J. V. de Julián-Ortiz and H. A. Mohammad-Salim, Mechanistic insights into the regio- and stereoselectivity of [3+2] cycloaddition reactions between N -methyl-phenylnitrone and trans -1-chloro-2-nitroethylene within the framework of molecular electron density theory, New J. Chem., 2025, 49, 11191–11202.
29. N. N. Cyrille, A. A. Idrice, M. B. Maraf, F. A. Charles, M. N. Ibrahim, M. Ríos-Gutierrez and L. R. Domingo, Understanding the Mechanism of Nitrobenzene Nitration with Nitronium Ion: A Molecular Electron Density Theory Study, ChemistrySelect, 2019, 4, 13313–13319.
30. R. G. Parr and W. Yang, Density functional approach to the frontier-electron theory of chemical reactivity, J. Am. Chem. Soc., 1984, 106, 4049–4050.
31. W. Yang and W. J. Mortier, The use of global and local molecular parameters for the analysis of the gas-phase basicity of amines, J. Am. Chem. Soc., 1986, 108, 5708–5711.
32. L. R. Domingo, M. J. Aurell, P. Pérez and R. Contreras, Quantitative Characterization of the Local Electrophilicity of Organic Molecules. Understanding the Regioselectivity on Diels-Alder Reactions, J. Phys. Chem. A, 2002, 106, 6871–6875.
33. P. Pérez, L. R. Domingo, M. Duque-Noreña and E. Chamorro, A condensed-to-atom nucleophilicity index. An application to the director effects on the electrophilic aromatic substitutions, J. Mol. Struct.:THEOCHEM, 2009, 895, 86–91.
34. R. F. W. Bader, Molecular fragments or chemical bonds, Acc. Chem. Res., 1975, 8, 34–40.
35. A. D. Becke and K. E. Edgecombe, A simple measure of electron localization in atomic and molecular systems, J. Chem. Phys., 1990, 92, 5397–5403.
36. J. Andres, S. Berski, L. R. Domingo, V. Polo and B. Silvi, Describing the Molecular Mechanism of Organic Reactions by Using Topological Analysis of Electronic Localization Function, Curr. Org. Chem., 2011, 15, 3566–3575.
37. L. R. Domingo, M. Ríos-Gutiérrez, E. Chamorro and P. Pérez, Aromaticity in Pericyclic Transition State Structures? A Critical Rationalisation Based on the Topological Analysis of Electron Density, ChemistrySelect, 2016, 1, 6026–6039.
38. X. Krokidis, S. Noury and B. Silvi, Characterization of Elementary Chemical Processes by Catastrophe Theory, J. Phys. Chem. A, 1997, 101, 7277–7282.
39. A. I. Adjieufack, I. M. Ndassa, I. Patouossa, J. K. Mbadcam, V. S. Safont, M. Oliva and J. Andrés, On the outside looking in: rethinking the molecular mechanism of 1,3-dipolar cycloadditions from the perspective of bonding evolution theory. The reaction between cyclic nitrones and ethyl acrylate, Phys. Chem. Chem. Phys., 2017, 19, 18288–18302.
40. A. I. Adjieufack, M. Mbah Bake, J. Ketcha Mbadcam, I. Mbouombouo Ndassa, J. Andrés, M. Oliva and V. S. Safont, How effectively bonding evolution theory retrieves and visualizes curly arrows: The cycloaddition reaction of cyclic nitrones, Int. J. Quantum Chem., 2019, 119(19), e25985.
41. L. Domingo, Molecular Electron Density Theory: A Modern View of Reactivity in Organic Chemistry, Molecules, 2016, 21, 1319.
42. E. R. Johnson, S. Keinan, P. Mori-Sánchez, J. Contreras-García, A. J. Cohen and W. Yang, Revealing Noncovalent Interactions, J. Am. Chem. Soc., 2010, 132, 6498–6506.
43. L. Domingo, M. Ríos-Gutiérrez and P. Pérez, Applications of the Conceptual Density Functional Theory Indices to Organic Chemistry Reactivity, Molecules, 2016, 21, 748.
44. L. R. Domingo, M. Ríos-Gutiérrez, B. Silvi and P. Pérez, The Mysticism of Pericyclic Reactions: A Contemporary Rationalisation of Organic Reactivity Based on Electron Density Analysis, European J. Org. Chem., 2018, 2018, 1107–1120.
45. J. B. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, G. A. Petersson, H. Nakatsuji, X. Li, M. Caricato, A. V. Marenich, J. Bloino, B. G. Janesko, R. Gomperts, B. Mennucci and H. P. Hratchian, Gaussian 16, Revision B.01, 2016, https://gaussian.com/gaussian16/.
46. Y. Zhao and D. G. Truhlar, The M06 suite of density functionals for main group thermochemistry, thermochemical kinetics, noncovalent interactions, excited states, and transition elements: two new functionals and systematic testing of four M06-class functionals and 12 other function, Theor. Chem. Acc., 2008, 120, 215–241.
47. S. Grimme, Density functional theory with London dispersion corrections, WIREs Comput. Mol. Sci., 2011, 1, 211–228.
48. J. Tomasi, B. Mennucci and R. Cammi, Quantum Mechanical Continuum Solvation Models, Chem. Rev., 2005, 105, 2999–3094.
49. J. R. Lane, J. Contreras-García, J.-P. Piquemal, B. J. Miller and H. G. Kjaergaard, Are Bond Critical Points Really Critical for Hydrogen Bonding?, J. Chem. Theory Comput., 2013, 9, 3263–3266.
50. J. Contreras-García, E. R. Johnson, S. Keinan, R. Chaudret, J.-P. Piquemal, D. N. Beratan and W. Yang, NCIPLOT: A Program for Plotting Noncovalent Interaction Regions, J. Chem. Theory Comput., 2011, 7, 625–632.
51. T. Lu and F. Chen, Multiwfn: A multifunctional wavefunction analyzer, J. Comput. Chem., 2012, 33, 580–592.
52. S. Noury, X. Krokidis, F. Fuster and B. Silvi, Computational tools for the electron localization function topological analysis, Comput. Chem., 1999, 23, 597–604.
53. V. Liégeois, DrawProfile, UNamur, 2024 https://www.unamur.be/en/sciences/chemistry/research/ucpts/lct/drawsuite/drawprofile.
54. V. Liégeois, DrawMol, UNamur, 2024, https://www.unamur.be/en/sciences/chemistry/research/ucpts/lct/drawsuite/drawmol.
55. M. Chellegui, R. M. Al-Mokhtar, R. N. Salih, L. Benhamed, S. Benmetir, J. Vicente de Julián-Ortiz, H. A. Mohammad-Salim and A. Ben Ahmed, Mechanistic insights into the (3 + 2) cycloaddition of azomethine ylide with dimethyl acetylenedicarboxylate via bond evolution theory, RSC Adv., 2025, 15, 29666–29679.
56. M. Chellegui, L. Benhamed, R. N. Salih, I. Salhi, S. Benmetir, A. Ben Ahmed, H. A. Mohammad-Salim and J. V. de Julián-Ortize, MEDT insights into the mechanism and selectivity of the (3 + 2) cycloaddition of (Z)- N -methyl- C -(2-furyl)-nitrone with but-2-ynedioic acid and the bioactivity of the reaction products, RSC Adv., 2025, 15, 32271–32283.
57. A. E. Reed, L. A. Curtiss and F. Weinhold, Intermolecular interactions from a natural bond orbital, donor-acceptor viewpoint, Chem. Rev., 1988, 88, 899–926.
58. M. Chellegui, S. Koudjina, I. Salhi, S. Benmetir, R. N. Salih, H. A. Mohammad-Salim, G. Y. S. Atohoun and J. V. de Julián-Ortiz, Mechanism and stereoselectivity of a [3 + 2] cycloaddition involving a glucosyl nitrone: a MEDT study, Org. Biomol. Chem., 2025, 23, 5016–5032.
59. R. N. Salih, H. Mohammad-Salim and M. Algso, Mechanistic study of N-t-butyl nitrone and cyanoacetylene [3 + 2] cycloaddition: a combined DFT, docking, and ADMET approach, Monatshefte für Chemie - Chem. Mon., 2025, 156, 443–456.
60. L. R. Domingo, M. Ríos-Gutiérrez and P. Pérez, A molecular electron density theory study of the participation of tetrazines in aza-Diels–Alder reactions, RSC Adv., 2020, 10, 15394–15405.
61. L. R. Domingo, M. Ríos-Gutiérrez and P. Pérez, Electrophilicity w and Nucleophilicity N Scales for Cationic and Anionic Species, Sci. Radices, 2025, 4, 1–17.
62. L. R. Domingo, M. J. Aurell, P. Pérez and R. Contreras, Quantitative characterization of the global electrophilicity power of common diene/dienophile pairs in Diels–Alder reactions, Tetrahedron, 2002, 58, 4417–4423.
63. L. R. Domingo, P. Pérez and J. A. Sáez, Understanding the local reactivity in polar organic reactions through electrophilic and nucleophilic Parr functions, RSC Adv., 2013, 3, 1486–1494.
64. L. Boltzmann, Ableitung des Stefan’schen Gesetzes, betreffend die Abhängigkeit der Wärmestrahlung von der Temperatur aus der electromagnetischen Lichttheorie, 1884, vol. 258, pp. 291–294.
65. R. P. Bell, The theory of reactions involving proton transfers, Proc. R. Soc. London. Ser. A – Math. Phys. Sci., 1936, 154, 414–429.
66. M. G. Evans and M. Polanyi, Inertia and driving force of chemical reactions, Trans. Faraday Soc., 1938, 34, 11.
67. M. G. Evans and M. Polanyi, Further considerations on the thermodynamics of chemical equilibria and reaction rates, Trans. Faraday Soc., 1936, 32, 1333.
68. R. Jasiński, A stepwise, zwitterionic mechanism for the 1,3-dipolar cycloaddition between (Z)-C-4-methoxyphenyl-N-phenylnitrone and gem-chloronitroethene catalysed by 1-butyl-3-methylimidazolium ionic liquid cations, Tetrahedron Lett., 2015, 56, 532–535.
69. L. R. Domingo and M. Ríos-Gutiérrez, A Useful Classification of Organic Reactions Based on the Flux of the Electron Density, Sci. Radices, 2023, 2, 1–24.
70. P. Vermeeren, T. A. Hamlin, I. Fernández and F. M. Bickelhaupt, Origin of rate enhancement and asynchronicity in iminium catalyzed Diels–Alder reactions, Chem. Sci., 2020, 11, 8105–8112.
71. P. Vermeeren, S. C. C. van der Lubbe, C. Fonseca Guerra, F. M. Bickelhaupt and T. A. Hamlin, Understanding chemical reactivity using the activation strain model, Nat. Protoc., 2020, 15, 649–667.
72. P. Vermeeren, T. A. Hamlin and F. M. Bickelhaupt, Origin of asynchronicity in Diels–Alder reactions, Phys. Chem. Chem. Phys., 2021, 23, 20095–20106.
73. W.-J. van Zeist and F. M. Bickelhaupt, The activation strain model of chemical reactivity, Org. Biomol. Chem., 2010, 8, 3118.
74. F. M. Bickelhaupt and K. N. Houk, Analyzing Reaction Rates with the Distortion/Interaction-Activation Strain Model, Angew. Chemie Int. Ed., 2017, 56, 10070–10086.
75. D. Svatunek and K. N. Houk, autoDIAS: a python tool for an automated distortion/interaction activation strain analysis, J. Comput. Chem., 2019, 40, 2509–2515.

---