A classical potential-based framework for modeling mechanochemical reactivity via molecular distortion: demonstration for a Diels–Alder reaction

Abstract

Mechanochemical reactions are increasingly studied using molecular dynamics simulations to understand mechanically activated chemical transformations. However, accurately capturing reactivity under mechanochemical conditions using classical potentials remains a challenge because standard models inhibit force-induced distortion of reactant species. In this study, we used the REACTER protocol, a method for simulating reactive events via dynamic bond changes, with a classical potential modified to allow the molecular distortion observed in first-principles calculations of a 4 + 2 Diels–Alder cycloaddition reaction. The approach was used to simulate the reaction in non-mechanochemical conditions with a solvent and no external stress, as well as in mechanochemical conditions. Mechanochemical simulations were run at hydrostatic stresses of 0.1 MPa and 2.5 GPa, both with and without shear applied, to investigate how the stress state influences reactivity. Relative to the non-mechanochemical reference case, hydrostatic stress and shear stress increased reaction yield. This increase was due to molecular distortion, the primary mechanism by which mechanical force activates chemical reactions, that could only be modeled using the modified classical potential. However, some of the increase in reaction yield was attributable to secondary mechanochemical activation mechanisms. Specifically, hydrostatic stress decreased the distance between reactants and shear stress facilitated alignment of reactants in the direction of imposed shear. This work provides new insight into how the stress state affects mechanochemical reaction mechanisms and establishes a generally applicable framework for improving classical potential-based simulations for organic reactions.

---

Introduction

Mechanochemical reactions are chemical transformations accelerated by mechanical force.^1,2 These reactions often involve the direct application of mechanical energy to break, form, or rearrange chemical bonds, leading to new chemical products. The study of these processes^3–5 and their mechanisms^6–8 is called mechanochemistry.^6,9–11 Mechanochemistry is a powerful synthetic approach as it can access reaction pathways, mechanisms, selectivities, and rates distinct from those achieved through conventional thermal, photochemical, or electrochemical methods, often with less input energy and little to no use of wasteful solvent.^12–16

One means by which mechanical force can facilitate chemical reactions is distorting reactant molecules.^17–19 Mechanical force can stretch, compress, twist, or bend molecular bonds and bond angles, thereby changing the activation energy for a reaction to proceed. Experimentally, mechanochemical processes can be achieved through various methodologies, including ball milling,^2,20–22 ultrasonic cavitation,^23–25 molecular pulling techniques,^26–28 and local probe methods.^29,30 Recently, continuous-flow methods such as twin-screw extrusion³¹ and resonant acoustic mixing^32–34 have emerged as powerful tools for scaling up mechanosynthesis, enabling solvent-free processing and industrial translation.³⁵ In parallel, the development of in situ³⁶ and in operando techniques, including X-ray diffraction,^37,38 Raman spectroscopy,^39–41 and synchrotron-based^37,42 methods, have provided critical insights into reaction intermediates and mechanisms, allowing chemical transformations under mechanical activation to be monitored in real time.^43–45 However, the direct measurement of molecular distortion remains a challenge with such techniques.

Computational tools have emerged as valuable resources for advancing this field.^46,47 Tools such as Constrained Geometries simulate External Force⁴⁸ (CoGEF), Force-Modified-Potential Energy Surfaces⁴⁹ (FMPES), External Force Explicitly Included⁵⁰ (EFEI), and eXtended Hydrostatic Compression Force Field⁵¹ (X-HCFF) extend the capabilities of Density Functional Theory^52,53 (DFT) by simulating the mechanical force-induced distortion of molecules at the electronic structure level.^54,55 Recent investigations have extended DFT-based methods by introducing a mechanochemical reaction constant to rationalize reactivity trends⁵⁶ and by combining DFT with microkinetic modeling⁵⁷ to reproduce experimental observations of ball-milling reactivity. Such quantum mechanical methods allow precise modeling of bond elongation,^58,59 rupture,^60,61 and formation,^62,63 and mechanochemical pathways.^63–65 While quantum mechanical methods provide mechanistic understanding into force-induced molecular-level processes,⁴⁶ they are limited to relatively small system sizes and typically provide information only about static or short-timescale phenomena. This is a challenge for computational mechanochemistry because it means that bulk-scale phenomena involving large numbers of molecules, such as solvent interactions, collective thermal response, slow molecular rearrangement, cannot be modeled.

Bulk-scale effects can instead be captured within the tens of nanometer and nanosecond scales accessible to classical molecular dynamics (MD) simulations. However, chemical reactions can only be captured by MD simulations if a reactive potential, i.e., an empirical potential that allows the formation and dissociation of covalent bonds, is used. Reactive MD simulations can model bulk material deformation, distortion of individual molecules, and the resulting chemical reactions. This has been done previously using the ReaxFF potential for reactions including oligomerization reactions of allyl alcohol,⁶⁶ α-pinene,^67,68 and cyclohexene,⁶⁹ oxygen release from polymethacrylate-endoperoxide,⁷⁰ perfluoropolyether (PFPE) degradation,⁷¹ and tribochemical wear of silicon.⁷² While ReaxFF has provided valuable insight into force-induced reaction mechanisms, its applicability for mechanochemistry is still constrained by its lack of transferability, that is, its reduced accuracy when applied to chemical systems or conditions other than those used in its original parameterization. This limitation is particularly relevant in mechanochemical systems, since ReaxFF potentials are typically not parameterized on systems where chemical transformations are triggered by mechanical force.

An alternative approach that combines the accuracy of DFT with the ability of reactive potentials like ReaxFF to model bulk-scale effects is the use of classical potentials with parameters modified based on DFT calculations of mechanochemical reactions. Some prior studies have implemented this approach with custom scripting within the large-scale atomic/molecular massively parallel simulator (LAMMPS)⁷³ to enable the classical potential bond formation or breaking.^74,75 In custom scripted methods, bond changes are typically hard-coded and tied to specific simulation steps or distance cutoffs, requiring manual updates to topology.

The REACTER protocol can overcome this limitation as it enables dynamic topology changes in LAMMPS simulations based on user-defined geometric and probabilistic constraints.^76–79 The protocol, which supports arbitrary reactions for any classical force field, directly converts pre-reaction topologies to post-reaction ones without extensive parameterization of reaction pathways. REACTER has been validated for a wide range of a non-mechanochemical bulk-scale reactions, including polymerization,⁸⁰ catalysis,⁸¹ and nanoscale morphology studies,⁸² demonstrating its versatility and efficiency in modeling chemical transformations. However, standard implementations of REACTER rely on classical potentials that do not allow for significant molecular distortion under applied force because the repulsive part of conventional non-bonded interaction models restrict atoms from coming into bonding proximity and standard dihedral potentials are too stiff to permit the angular distortions beyond the harmonic range. As a result, standard implementations of REACTER may not fully capture the role of molecular distortion in driving mechanochemical activation.

Here, we address this challenge by modifying the classical potential underlying REACTER to facilitate force-induced distortion of bonds and angles, allowing for a more explicit representation of distortion-induced reactivity within a classical potential-based framework. We use a Diels–Alder reaction as a testbed to demonstrate the approach and investigate how mechanochemical reactivity depends on the stress state. Specifically, we consider the mechanochemical 4 + 2 Diels–Alder cycloaddition reaction (Fig. 1) previously studied by Zhang et al.⁸³ They reported that the 4 + 2 Diels–Alder cycloaddition of cyclopentadiene (diene) with N-4-tosylmaleimide (dienophile) resulted in more than 99% yield when carried out under controlled mechanochemical milling conditions, significantly higher than the 14% yield observed for the solution-phase reaction conducted under non-mechanochemical conditions with hexane as the solvent.

Fig. 1 4 + 2 cycloaddition reaction between cyclopentadiene and N-4-tosylmaleimide.

To model this reaction, first, CoGEF calculations were performed on a single reactant pair to characterize molecular distortion during the reaction. Then, a classical MD potential was modified to enable molecular distortion by adjusting parameters for non-bonded interactions and dihedral angles containing atoms that participate in the reaction. Using the modified potential, we performed large-scale molecular dynamics simulations using the REACTER protocol in LAMMPS. We simulated the above-mentioned Diels–Alder cycloaddition reaction in non-mechanochemical conditions with and without solvent, and mechanochemical conditions at three different stress states. Reactivity was quantified by counting bond formation events during the simulations. The comparison between reaction yield in mechanochemical and non-mechanochemical solvent conditions enabled qualitative comparison to experimental results. Comparison of the results from mechanochemical simulations in various stress states enabled analysis of how hydrostatic stress and shear stress can influence reaction yield. To understand the mechanistic basis of reactivity trends under mechanochemical conditions, we analyzed the primary effect of force in distorting reactant species, as well as the secondary effects of force to increase density and shear to induce molecular alignment. The results demonstrate the effectiveness of the new approach for modeling mechanochemical reactions at the bulk scale and provide insight into the mechanisms by which mechanical force activates chemical reactions.

Computational details

To explore the molecular-level reaction energy profile for the Diels–Alder cycloaddition reaction, first, DFT calculations were performed using the Gaussian 16 program package.⁸⁴ The M06-2X^85,86 functional was employed with the 6-311+G(d,p) basis set. The CoGEF method was used to constrain the distance between the reacting carbon atoms of cyclopentadiene and N-4-tosylmaleimide. These calculations provided the energy profile as well as information about molecular distortion that occurred along the reaction pathway. The molecular distortion was then used to inform a second set of calculations performed with constraints on the C–C=C–H dihedral angle of the N-4-tosylmaleimide and C–C=C–C of the cyclopentadiene.

Bulk-scale MD reactions were then performed using the REACTER protocol^76–79 implemented through the fix bond/react command (which enables formation or breaking of predefined chemical bonds) in LAMMPS.⁷³ The molecular templates and the reaction map file required for the REACTER protocol were prepared using the LUNAR program.⁸⁷ Within the REACTER protocol, a bond was defined to form when the distance between the reacting carbon atoms of cyclopentadiene and N-4-tosylmaleimide was smaller than 2.6 Å (based on the DFT/classical fitting results; discussed later). Trial runs were performed with different bond definition distances and the number of reactions increased with increasing distance, but trends were unaffected (Fig. S1).

The MD simulations utilized the PCFF-IFF⁸⁸ potential, which includes all relevant molecular interactions for organic systems, including bond, angle, dihedral, and improper terms. The potential describes non-bonded interactions using Lennard-Jones and coulombic potentials with a cutoff distance of 12.0 Å. Long-range electrostatic interactions were computed using the particle–particle particle–mesh (PPPM) method with a relative accuracy of 10⁻⁴.

The system consisted of 200 cyclopentadiene molecules and 200 N-4-tosylmaleimide molecules. Periodic boundary conditions were applied in all three dimensions. To integrate Newton's equations of motion, the velocity-Verlet algorithm was used with a 1 fs timestep. The equilibration procedure was performed in two stages. First, a simulation was run in the NPT ensemble (constant number of atoms, pressure, and temperature) at 0.1 MPa (atmospheric pressure) and 300 K for 3 ns to allow the system to reach a stable density. This was followed by a 0.5 ns NVT (constant number of atoms, volume, and temperature) simulation run at 300 K to allow the system to reach a stable energy while maintaining the temperature. After equilibration, the REACTER protocol was activated so that the target chemical reaction could occur.

Five different conditions were modeled in the production simulations. First, non-mechanochemical simulations at 0.1 MPa without imposed shear strain were performed both with and without the addition of 400 hexane molecules (hexane was the solvent used in the previous experimental study⁸³). Then, mechanochemical simulations were run with the system subjected to hydrostatic stress and shear stress. Hydrostatic stress was applied by gradually compressing the system in all three spatial directions over a period of 4 ns until the target pressure was achieved, and shear stress was applied by deforming the xy plane of the simulation cell at a constant shear rate of 10⁸ s⁻¹. While this shear rate is higher than in some mechanochemical experiments, there are examples of recent mechanochemical studies approaching this rate.⁸⁹ High shear rates are typical in nonequilibrium molecular dynamics (NEMD) simulations where the short time scale necessitates fast deformation rates to observe measurable responses within computationally feasible simulation durations. However, such simulations can still provide qualitative mechanistic insight into bond distortion and reaction pathways.

Three different combinations of hydrostatic and shear stress were simulated, as shown in Fig. 2: 0.1 MPa hydrostatic stress without imposed shear strain, 2.5 GPa hydrostatic stress without imposed shear strain, and 2.5 GPa hydrostatic stress with shear strain. Each simulation was repeated twice starting from the initial step of packing molecules into the periodic simulation cell. These simulations can be considered to represent, for example, a nanoscale portion of the material that is compressed and sheared during ball milling, perhaps at a contact between two protruding asperities. While the exact stress distributions in ball milling are difficult to measure, applying controlled hydrostatic and shear stresses in simulations allows us to systematically probe their individual and combined effects on molecular distortion and reactivity.

Fig. 2 Simulation periodic boxes containing cyclopentadiene and N-4-tosylmaleimide at three conditions: (A) 0.1 MPa hydrostatic stress without imposed shear, (B) 2.5 GPa hydrostatic stress without imposed shear, and (C) 2.5 GPa hydrostatic stress with imposed shear.

Results

First principles calculations of molecular distortion

We performed DFT calculations to obtain the relaxed energy profile in Fig. 3A by systematically varying the distance between reacting carbon atoms (dashed lines in Fig. 3B). The reaction proceeds through an energy barrier of 13.64 kcal mol⁻¹ and ultimately leads to the product formation.

Fig. 3 (A) Energy profile diagram as a function of the reacting C–C distance from DFT for equilibrium conditions (maroon circles) and distorted conditions (purple triangles). The distorted conditions correspond to fixing the C–C=C–H dihedral angles of N-4-tosylmaleimide and C–C=C–C of cyclopentadiene, highlighted in orange in the inset, to 159° and 10°, respectively. The energy from the modified classical potential is shown as a solid purple line. (B) Optimized geometries of the (Ia) reactant state without molecular distortion, (Ib) transition state without molecular distortion, and (II) reactant state with molecules constrained to a distorted conformation. The dashed line indicates the reacting C–C distance. The pink lines highlight the bending of the molecular planes near the reactive sites.

To visualize the geometric changes associated with this reaction, the optimized geometries of the reactant state (C–C = 3.3 Å) and the transition state (C–C = 2.25 Å) are shown in Fig. 3B(Ia) and (Ib). The bond formation between the reacting carbon atoms is accompanied by distortion in the five-membered cyclopentadiene ring and changes in dihedral angles of the N-4-tosylmaleimide. In the distorted geometry (Fig. 3B(Ib)), the molecular plane of the diene becomes nonplanar and bends toward the dienophile (emphasized by the pink lines adjacent to the snapshots), while the hydrogen atoms of the reacting carbon atoms on the dienophile are displaced in the opposite direction of the diene. This is quantified by the C–C=C–H dihedral angle of the N-4-tosylmaleimide molecule and the C–C=C–C dihedral angle of cyclopentadiene molecule (highlighted in orange in the inset of Fig. 3A). Both angles increase with decreasing C–C distance; shown for the C–C=C–H dihedral angle of the N-4-tosylmaleimide in Fig. S3. At the transition state, the C–C=C–H dihedral angle of N-4-tosylmaleimide is 153° and the C–C=C–C dihedral angle of cyclopentadiene is 17°, distorted from their equilibrium values of 180° and 0°, respectively.

The primary effect of mechanical force to activate chemical reactions has been reported to be distortion of the reactant species.^17,90 To quantify this effect in our calculations, a second set of DFT calculations was performed to obtain the energy profile, but with the dihedral angles of the reactants constrained away from their equilibrium values. The dihedral angles were fixed at values just below the angles at the transition state. Particularly, the C–C=C–H dihedral angle of N-4-tosylmaleimide was set to 159° and the C–C=C–C dihedral angle of cyclopentadiene molecule was set to 10°, as illustrated in Fig. 3B(II). These values correspond to distortions of 21° and 10°, respectively, and were chosen to match the distorted conformations of the molecules just before the transition state. As shown in Fig. 3A, this distortion increases the reactant energy and therefore decreases the energy required for the reaction to proceed, which mimics the effects of external mechanical force on the reaction pathway.

Classical potential modified for mechanochemistry

To enable force-induced distortions in both the diene and dienophile, we modified the classical potential to bias reactive moieties toward strained geometries prior to bond formation. These modifications were necessary because the repulsive part of conventional non-bonded interaction models does not allow atoms to come into bonding proximity under mechanical force, and standard dihedral potentials are too stiff to permit the angular distortions that promote reactivity. To address the proximity issue, we employed a hybrid/overlay potential. Specifically, for C–C pairs between 0.5 and 5.0 Å, a Morse-style interaction⁹¹ was added to the Lennard-Jones and coulombic models that capture long-range van der Waals and electrostatic effects in the classical potential. Beyond ∼5 Å, the Morse potential flattens, approaching its asymptotic limit where interatomic forces vanish and the atoms are effectively non-interacting. This choice is also supported by DFT calculations, which showed no significant interactions between the two molecules at separations greater than ∼3.3 Å. The Morse model is only attractive at intermediate distances, which lowers the energy barrier for close-range approach and allows molecules to more readily come into bonding proximity. Beyond 5.0 Å (up to the global cutoff of 12 Å), all atom pairs interacted through standard Lennard-Jones and coulombic forces.

To address the stiffness issue, the force constants in the classical potential were modified for dihedral angles containing atoms that participate in the reaction in both reacting molecules (shown in Fig. S2). Particularly, the amplitude and phase angle of the classical potential dihedrals were iteratively adjusted such that the energy calculated by the classical potential was consistent with the DFT (distorted) values until a C–C distance of approximately 2.6 Å. Fitting was only necessary up to an atomic separation of 2.6 Å because, below that separation, a covalent bond forms within the REACTER protocol and the reaction is complete. The modified classical potential reproduced the distorted DFT energy profile with an root mean square error (RMSE) of 1.482 kJ mol⁻¹, an average absolute error of 1.170 kJ mol⁻¹, and a maximum deviation of 2.716 kJ mol⁻¹. These values are on the order of the energy due thermal fluctuations at room temperature (k_BT ∼ 2.5 kJ mol⁻¹), suggesting the modified potential captured DFT-level distortion energetics with sufficient accuracy for use in classical MD simulations. The modified potential file is provided in the SI.

MD simulations with REACTER were initially run both with and without the modifications to the underlying classical potential. Distortion was quantified as the deviation from equilibrium of the C–C=C–H dihedral angle for all N-4-tosylmaleimide molecules in the last 2 ns of each MD simulation (although both the cyclopentadiene and N-4-tosylmaleimide underwent distortion prior to reaction, we used N-4-tosylmaleimide to quantify distortion, as it is more distorted than the cyclopentadine at the transition state). Recall from Fig. 3 that the deviation of this angle from its equilibrium value at the transition state, i.e., distortion was 180° − 153° = 27°. As shown in Fig. 4, with the unmodified classical potential, the dihedral angles of most N-4-tosylmaleimide molecules were less than 10° and no molecules reached the 27° distortion expected for the transition state based on DFT results. In contrast, most molecules in the simulations using the modified potential had a dihedral angle exceeding 27°. This indicates that the classical potential lacks the ability to reproduce the extent of distortion required for the Diels–Alder reaction to proceed, while the modified potential successfully allows the necessary dihedral distortions.

Fig. 4 Distribution of distortion in the C–C=C–H dihedral of N-4-tosylmaleimide molecules under hydrostatic stress of 2.5 GPa with imposed shear strain computed using classical MD simulations. Distortion is quantified as the deviation of the dihedral angle from equilibrium, based on the last 2 ns of the trajectory. Histograms compare results from simulations with the unmodified (maroon solid) and modified (purple hatched) classical potentials. The vertical dashed line at 27° represents the transition-state distortion predicted by DFT, serving as a reference for the amount of conformational distortion needed to activate the reaction. Each bar reflects the number of molecules exhibiting a given level of distortion, with 200 molecules analyzed per simulation.

Reactive MD simulations of mechanochemistry

Using the modified potential, we conducted simulations at five different conditions: in non-mechanochemical conditions at 0.1 MPa (atmospheric pressure) without imposed shear strain, either neat or with hexane solvent, and in mechanochemical conditions at 0.1 MPa without shear strain, 2.5 GPa without imposed shear strain, and 2.5 GPa with shear strain. In the cases with shear strain applied, the shear stress initially increased as the system responded to the imposed strain, then reached a steady state within approximately 3–5 ns (Fig. S4). The steady-state shear stress was 104.6 MPa (standard deviation of 62.2 MPa) for the 0.1 MPa hydrostatic stress case and 339.3 MPa (standard deviation of 86.7 MPa) for the 2.5 GPa hydrostatic stress case.

The number of reactions was tracked over time at each condition as shown in Fig. 5. At 0.1 MPa without imposed shear, the number of reactions was very low, both for neat and hexane solvent conditions (yields of 1.5% and 2.5% of the molecules reacting, respectively). Relative to these non-mechanochemical conditions, the application of shear alone increased the number of reactions slightly (5.5% yield), hydrostatic stress increased yield even more (10%), and both shear and hydrostatic stress increased the yield markedly (27.5%). Under solvent-based mechanochemical conditions (hexane), we observed the same qualitative ordering, with only small shifts in absolute yields (Fig. S5). Notably, no reactions occurred in simulations with the unmodified potential at any of these conditions because the potential inhibits the molecular distortion necessary for reacting atoms to be close enough to bond.

Fig. 5 Number of reactions observed under non-mechanochemical and mechanochemical conditions, simulated using the REACTER protocol with a modified classical force field. A total of 200 reactions is possible at each condition. Lines and symbols represent the average number of reactions across two independent simulations and shaded regions indicate the corresponding standard deviation. Dashed lines represent 0.1 MPa conditions, solid lines represent 2.5 GPa; open and filled circles indicate non-shear and shear simulations, respectively; the dashed black line corresponds to simulations with hexane solvent.

Force-induced molecular distortion is a primary effect of mechanochemical conditions.^17,90 In our simulations, this is modeled by modifying classical potentials to favor distorted geometries relevant to mechanochemical reactions. However, since these modifications are applied uniformly across all conditions, distortion alone cannot explain the observed differences in reactivity between shear and non-shear or low- and high-pressure environments (Fig. 5). This suggests that secondary mechanochemical mechanisms may be contributing. To isolate these mechanisms, we examined the effects of hydrostatic stress (pressure) and shear stress independently.

First, we analyzed the effect of hydrostatic stress by comparing simulations run at 0.1 MPa and 2.5 GPa without imposed shear. It is well established that elevated pressure can enhance chemical reactivity by bringing molecules into closer proximity, thereby increasing local density and the likelihood of reactive collisions.^92,93 To determine if this mechanisms is relevant in our system, we analyzed the distribution of cavity radii under different pressure conditions, comparing non-shear simulations at 0.1 MPa and 2.5 GPa. Cavity radius distributions were calculated by identifying the largest empty spherical region surrounding each atom without overlapping neighboring atoms, providing a measure of local free volume.⁹⁴ As shown in Fig. 6, there is a clear distinction between low- and high-pressure environments. At 0.1 MPa, the system exhibits relatively larger cavities, implying lower local density and greater intermolecular separation, which can reduce the probability of reactive collisions. In contrast, at 2.5 GPa, the cavities are smaller, reflecting a more compact molecular arrangement that enhances the likelihood of reactive encounters. This observation is further supported by the calculated densities: the system at 0.1 MPa has a density of 0.989 g cm⁻³, whereas at 2.5 GPa, the density increases to 1.223 g cm⁻³. The pressure-induced reduction in free volume thus contributes to increased reactivity by bringing reactive moieties closer.

Fig. 6 Histograms of the cavity radius in the system at 0.1 MPa and 2.5 GPa without shear, computed over the last 0.5 ns of the NVT stage of the equilibration process.

In addition to a secondary effect due to pressure-induced proximity, another secondary effect arises from shear strain, which in this case enhances reactivity by promoting favorable molecular orientations for reactants. This can be observed qualitatively from the snapshots in Fig. 7A that visually show increased alignment of molecules in the shear simulation vs. the simulation without shear. We quantified this by defining alignment as the absolute value of the dot product between the normal vector of a molecule's plane (planar five member ring) and the unit vector along the y-axis. A value of 1 indicates that the molecular plane is aligned along the shear direction (x-axis), while a value of 0 corresponds to a molecular plane oriented perpendicular to it. The ring plane is refit to the positions of the atoms in the five-membered ring of both molecules, so any ring bending is captured in the alignment metric. The average per-frame alignment results for cyclopentadiene and N-4-tosylmaleimide under 2.5 GPa non-shear and shear conditions are shown in Fig. 7B. In the absence of shear, there is no preferred orientation of the molecular planes. Under shear, however, molecules reorient such that their molecular planes become more aligned with the shear direction, promoting more favorable alignment of the reactants to each other. This shear-induced alignment between molecular planes facilitates productive intermolecular contacts which contributes to the enhanced reactivity observed under 2.5 GPa shear conditions compared to non-shear conditions. As shown in Fig. S6, however, this shear-induced alignment effect is negligible at 0.1 MPa, indicating that the orientational influence of shear becomes significant only under higher pressure conditions.

Fig. 7 (A) Qualitative evidence of molecular alignment due to shear is seen in side-view snapshots of the simulations at 30 ns with and without imposed shear at 2.5 GPa. Each molecule is colored based on its alignment, as defined in the color scale legend. (B) Average molecular alignment relative to the shear direction at 2.5 GPa without imposed shear (green) and with imposed shear (red). Each symbol represents the per-frame average alignment across all molecules of a given type: circles for cyclopentadiene and triangles for N-4-tosylmaleimide. Lines show the running average of these per-frame values over time: solid for cyclopentadiene and dashed for N-4-tosylmaleimide.

Recall that the previous experimental study of the reaction in our simulations reported 14% yield for the solution-phase reaction with hexane solvent and 99% yield under controlled mechanochemical milling.⁸³ While the simulations do not directly model the experimental conditions, and the term “yield” should not be compared quantitatively between simulations and experiments, there is qualitative agreement in the observed trends. Specifically, the simulations show a significant increase in the number of reactions with hydrostatic stress and shear and without solvent compared to at 0.1 MPa without imposed shear in hexane. This qualitative agreement provides evidence that the approach developed here using REACTER with a modified classical potential can effectively capture the enhanced reactivity and yield of mechanochemical conditions in previous experiments.

There is no previous experimental study comparing different mechanochemical conditions for the reaction modeled here. However, experimental measurements of film growth from zinc dialkyldithiophosphate (ZDDP) have shown that shear strongly drives mechanochemical film formation,^89,95 consistent with our findings. In one study, shear and compressive stresses were isolated by testing with different ratios of high-viscosity fluids, with shear rates calculated from data given in that paper to be as high as approximately 2 × 10⁷ s⁻¹.⁸⁹ The results showed that shear stress promoted film growth while compressive stress had the opposite effect, inhibiting film growth. This differs from our findings in which compression led to a small increase in yield for a cycloaddition reaction. The difference may be due to the fact that the film formation reaction starts with bond-breaking reactions that instigate adsorption of ZDDP precursor molecules on the sliding surface. Bond breaking is inhibited by compression but this effect, which is important for ZDDP film formation, is negligible for the cycloaddition reaction that occurs through bond formation only. Therefore, our observation that compression increases yield slightly while shear increases it dramatically is qualitatively consistent with the previous experimental study, further supporting the robustness of the simulation method reported here.

Conclusions

This study investigated the mechanochemical reactivity of a 4 + 2 Diels–Alder type cycloaddition reaction with MD simulations. The simulations used the REACTER protocol with a classical potential modified to facilitate force-induced molecular distortion. Simulations using modified potentials run under mechanochemical conditions resulted in higher reaction yield than those in non-mechanochemical conditions, in qualitative agreement with previous experiments. Simulations run at different pressures with and without imposed shear enabled differentiation of the effects of hydrostatic stress and shear stress. The increase in reaction yield under hydrostatic stress conditions was attributed to reduced cavity size and increased molecular density, which improved the likelihood of encounters between reactants. Shear further amplified this effect by inducing directional molecular alignment, facilitating favorable spatial arrangements of reactive moieties. These secondary mechanochemical effects can only be captured in MD simulations that enable modeling of many-molecule systems. Further, as demonstrated here, the MD simulations must use modified classical potentials to allow the primary mechanochemical effect of molecular distortion.

Looking ahead, these insights underscore the importance of integrating mechanical stress factors, such as pressure-induced compaction and shear-induced alignment, into computational models for reactive systems. Further, the approach demonstrated here can be applied to any mechanochemical reaction by modifying a classical potential based on distortions obtained from DFT calculations. This approach not only improves the predictive capability of molecular simulations but also provides a rational framework for designing stress-responsive materials and optimizing conditions for solvent-free or low-solvent mechanochemical synthesis.

Conflicts of interest

There are no conflicts to declare.

Data availability

The force field parameter files used in this study are provided as supplementary information (SI). All additional simulation-related files are available from the corresponding author upon reasonable request (E-mail: amartini@ucmerced.edu; ). Supplementary information: additional figures and analyses detailing the influence of reaction cutoff on reactivity, dihedral parameters for the modified potential, dihedral distortion as a function of C–C distance, shear stress evolution with time, solvent effects on reactivity, and molecular alignment under low-pressure conditions. See DOI: https://doi.org/10.1039/d5mr00099h.

Acknowledgements

This work was supported by the NSF Center for the Mechanical Control of Chemistry (CMCC), CHE-2303044. The CMCC is part of the Centers for Chemical Innovation Program. The authors thank Prof. Adam B. Braunschweig from CUNY Advanced Science Research Center for his helpful insights.

References

1. G.-W. Wang, Chem. Soc. Rev., 2013, 42, 7668–7700.
2. K. Kubota and H. Ito, Trends Chem., 2020, 2, 1066–1081.
3. K. J. Ardila-Fierro and J. G. Hernández, Angew. Chem., Int. Ed., 2024, 63, e202317638.
4. D. Tan and T. Friščić, Eur. J. Org Chem., 2018, 2018, 18–33.
5. M. Carta, E. Colacino, F. Delogu and A. Porcheddu, Phys. Chem. Chem. Phys., 2020, 22, 14489–14502.
6. F. Delogu and L. Takacs, J. Mater. Sci., 2018, 53, 13331–13342.
7. M. H. Barbee, T. Kouznetsova, S. L. Barrett, G. R. Gossweiler, Y. Lin, S. K. Rastogi, W. J. Brittain and S. L. Craig, J. Am. Chem. Soc., 2018, 140, 12746–12750.
8. P. A. Julien and T. Friscic, Cryst. Growth Des., 2022, 22, 5726–5754.
9. S. L. James, C. J. Adams, C. Bolm, D. Braga, P. Collier, T. Friščić, F. Grepioni, K. D. Harris, G. Hyett and W. Jones, et al. , Chem. Soc. Rev., 2012, 41, 413–447.
10. L. Takacs, Chem. Soc. Rev., 2013, 42, 7649–7659.
11. F. Cuccu, L. De Luca, F. Delogu, E. Colacino, N. Solin, R. Mocci and A. Porcheddu, ChemSusChem, 2022, 15, e202200362.
12. M. Wollenhaupt, M. Krupička and D. Marx, ChemPhysChem, 2015, 16, 1593–1597.
13. G. S. Kochhar, A. Bailey and N. J. Mosey, Angew. Chem., 2010, 122, 7614–7617.
14. J. L. Howard, M. C. Brand and D. L. Browne, Angew. Chem., 2018, 130, 16336–16340.
15. J. Batteas, K. G. Blank, E. Colacino, F. Emmerling, T. Friščić, J. Mack, J. Moore, M. E. Rivas and W. Tysoe, RSC Mechanochem., 2025, 2, 10–19.
16. S. Kumar, D. Button-Jennings, T. P. Hanusa and A. Martini, RSC Mechanochem., 2025, 2, 529–537.
17. Y. S. Zholdassov, L. Yuan, S. R. Garcia, R. W. Kwok, A. Boscoboinik, D. J. Valles, M. Marianski, A. Martini, R. W. Carpick and A. B. Braunschweig, Science, 2023, 380, 1053–1058.
18. Y. S. Zholdassov, R. W. Kwok, M. A. Shlain, M. Patel, M. Marianski and A. B. Braunschweig, RSC Mechanochem., 2024, 1, 11–32.
19. A. Martini and S. H. Kim, Tribol. Lett., 2021, 69, 1–14.
20. T. Friščić, C. Mottillo and H. M. Titi, Angew. Chem., 2020, 132, 1030–1041.
21. V. Martinez, T. Stolar, B. Karadeniz, I. Brekalo and K. Užarević, Nat. Rev. Chem., 2023, 7, 51–65.
22. I. R. Speight, K. J. Ardila-Fierro, J. G. Hernández, F. Emmerling, A. A. Michalchuk, F. García, E. Colacino and J. Mack, Nat. Rev. Methods Primers, 2025, 5, 29.
23. G. Cravotto, E. C. Gaudino and P. Cintas, Chem. Soc. Rev., 2013, 42, 7521–7534.
24. D. Yildiz, R. Göstl and A. Herrmann, Chem. Sci., 2022, 13, 13708–13719.
25. P. A. May and J. S. Moore, Chem. Soc. Rev., 2013, 42, 7497–7506.
26. J. Wang, T. B. Kouznetsova, Z. Niu, M. T. Ong, H. M. Klukovich, A. L. Rheingold, T. J. Martinez and S. L. Craig, Nat. Chem., 2015, 7, 323–327.
27. H. M. Klukovich, T. B. Kouznetsova, Z. S. Kean, J. M. Lenhardt and S. L. Craig, Nat. Chem., 2013, 5, 110–114.
28. R. T. O'Neill and R. Boulatov, Nat. Rev. Chem., 2021, 5, 148–167.
29. N. N. Gosvami, J. Bares, F. Mangolini, A. Konicek, D. Yablon and R. Carpick, Science, 2015, 348, 102–106.
30. J. R. Felts, A. J. Oyer, S. C. Hernández, K. E. Whitener Jr, J. T. Robinson, S. G. Walton and P. E. Sheehan, Nat. Commun., 2015, 6, 6467.
31. R. R. Bolt, J. A. Leitch, A. C. Jones, W. I. Nicholson and D. L. Browne, Chem. Soc. Rev., 2022, 51, 4243–4260.
32. A. A. Michalchuk, K. S. Hope, S. R. Kennedy, M. V. Blanco, E. V. Boldyreva and C. R. Pulham, Chem. Commun., 2018, 54, 4033–4036.
33. D. Kong, L. Yi, A. Nanni and M. Rueping, Nat. Commun., 2025, 16, 3983.
34. A. Nari, J. S. Ovens and D. L. Bryce, RSC Mechanochem., 2024, 1, 50–62.
35. E. Colacino, V. Isoni, D. Crawford and F. García, Trends Chem., 2021, 3, 335–339.
36. T. Rathmann, H. Petersen, S. Reichle, W. Schmidt, A. P. Amrute, M. Etter and C. Weidenthaler, Rev. Sci. Instrum., 2021, 92, 114102.
37. H. Petersen, S. Reichle, S. Leiting, P. Losch, W. Kersten, T. Rathmann, J. Tseng, M. Etter, W. Schmidt and C. Weidenthaler, Chem.–Eur. J., 2021, 27, 12558–12565.
38. N. Y. Gugin, K. V. Yusenko, A. King, K. Meyer, D. Al-Sabbagh, J. A. Villajos and F. Emmerling, Chem, 2024, 10, 3459–3473.
39. S. K. Hinojosa, T. Dogan, S. Fabig, S. Grätz and L. Borchardt, Chem.–Eur. J., 2025, e01336.
40. T. H. Borchers, F. Topić, M. Arhangelskis, M. Ferguson, C. B. Lennox, P. A. Julien and T. Friščić, Chem, 2025, 11.
41. L. Vugrin, C. Chatzigiannis, E. Colacino and I. Halasz, RSC Mechanochem., 2025, 2, 482–487.
42. I. Halasz, S. A. Kimber, P. J. Beldon, A. M. Belenguer, F. Adams, V. Honkimäki, R. C. Nightingale, R. E. Dinnebier and T. Friščić, Nat. Protoc., 2013, 8, 1718–1729.
43. I. Halasz, T. Friščić, S. A. Kimber, K. Užarević, A. Puškarić, C. Mottillo, P. Julien, V. Štrukil, V. Honkimäki and R. E. Dinnebier, Faraday Discuss., 2014, 170, 203–221.
44. S. Lukin, L. S. Germann, T. Friscic and I. Halasz, Acc. Chem. Res., 2022, 55, 1262–1277.
45. A. A. Michalchuk and F. Emmerling, Angew. Chem., Int. Ed., 2022, 61, e202117270.
46. T. Stauch and A. Dreuw, Chem. Rev., 2016, 116, 14137–14180.
47. J. Ribas-Arino and D. Marx, Chem. Rev., 2012, 112, 5412–5487.
48. M. K. Beyer, J. Chem. Phys., 2000, 112, 7307.
49. M. T. Ong, J. Leiding, H. Tao, A. M. Virshup and T. J. Martínez, J. Am. Chem. Soc., 2009, 131, 6377–6379.
50. J. Ribas-Arino, M. Shiga and D. Marx, Angew. Chem., Int. Ed., 2009, 48, 4190–4193.
51. T. Stauch, J. Chem. Phys., 2020, 153, 134503.
52. W. Kohn and L. J. Sham, Phys. Rev., 1965, 140, 1133–1138.
53. P. Hohenberg and W. Kohn, Phys. Rev., 1964, 136, 864–871.
54. S. Kumar, PhD thesis, Universität Bremen, 2023.
55. S. Kumar, G. Adabasi, M. Z. Baykara and A. Martini, J. Phys. Chem. C, 2024, 129, 369–378.
56. W. Sakai, L. Gonnet, N. Haruta, T. Sato and M. Baron, Phys. Chem. Chem. Phys., 2024, 26, 873–878.
57. B. S. Pladevall, A. de Aguirre and F. Maseras, ChemSusChem, 2021, 14, 2763–2768.
58. C. R. Wick, E. Topraksal, D. M. Smith and A.-S. Smith, Forces Mech., 2022, 9, 100143.
59. Z. Chen, X. Zhu, J. Yang, J. A. Mercer, N. Z. Burns, T. J. Martinez and Y. Xia, Nat. Chem., 2020, 12, 302–309.
60. I. M. Klein, C. C. Husic, D. P. Kovács, N. J. Choquette and M. J. Robb, J. Am. Chem. Soc., 2020, 142, 16364–16381.
61. S. Kumar and T. Stauch, RSC Adv., 2021, 11, 7391–7396.
62. S. Kumar, F. Zeller and T. Stauch, J. Phys. Chem. Lett., 2021, 12, 9470–9474.
63. S. Kumar, R. Weiß, F. Zeller and T. Neudecker, ACS Omega, 2022, 7, 45208–45214.
64. S. Aydonat, D. Campagna, S. Kumar, S. Storch, T. Neudecker and R. Göstl, J. Am. Chem. Soc., 2024, 146, 32117–32123.
65. L. J. Mier, G. Adam, S. Kumar and T. Stauch, ChemPhysChem, 2020, 21, 2402–2406.
66. J. Yeon, X. He, A. Martini and S. H. Kim, ACS Appl. Mater. Interfaces, 2017, 9, 3142–3148.
67. A. Khajeh, X. He, J. Yeon, S. H. Kim and A. Martini, Langmuir, 2018, 34, 5971–5977.
68. F. H. Bhuiyan, S. H. Kim and A. Martini, Appl. Surf. Sci., 2022, 591, 153209.
69. F. H. Bhuiyan, Y.-S. Li, S. H. Kim and A. Martini, Sci. Rep., 2024, 14, 2992.
70. J. Cobeña-Reyes, F. H. Bhuiyan and A. Martini, RSC Mechanochem., 2024, 1, 361–366.
71. H. Shekhar, S. Uchiyama, Y. Song, H. Zhang, K. Fukuzawa, S. Itoh and N. Azuma, RSC Mechanochem., 2025, 2, 745–755.
72. J. Wen, T. Ma, W. Zhang, G. Psofogiannakis, A. C. van Duin, L. Chen, L. Qian, Y. Hu and X. Lu, Appl. Surf. Sci., 2016, 390, 216–223.
73. A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. In’t Veld, A. Kohlmeyer, S. G. Moore and T. D. Nguyen, et al. , Comput. Phys. Commun., 2022, 271, 108171.
74. B. Koo, A. Chattopadhyay and L. Dai, Comput. Mater. Sci., 2016, 120, 135–141.
75. S. Kumar, B. Demir, A. Dellwisch, L. Colombi Ciacchi and T. Neudecker, Macromolecules, 2023, 56, 8438–8447.
76. J. R. Gissinger, B. D. Jensen and K. E. Wise, Polymer, 2017, 128, 211–217.
77. J. R. Gissinger, B. D. Jensen and K. E. Wise, Macromolecules, 2020, 53, 9953–9961.
78. J. R. Gissinger, B. D. Jensen and K. E. Wise, Comput. Phys. Commun., 2024, 304, 109287.
79. J. R. Gissinger, I. Nikiforov, Y. Afshar, B. Waters, M.-k. Choi, D. S. Karls, A. Stukowski, W. Im, H. Heinz and A. Kohlmeyer, et al. , J. Phys. Chem. B, 2024, 128, 3282–3297.
80. S. Zheng, J. Gissinger, B. S. Hsiao and T. Wei, ACS Appl. Mater. Interfaces, 2024, 16, 65677–65686.
81. M. S. Yeganeh, A. Jusufi, S. P. Deighton, M. S. Ide, M. Siskin, A. Jaishankar, C. Maldarelli, P. Bertolini, B. Natarajan and J. L. Vreeland, et al. , Sci. Adv., 2022, 8, eabm0144.
82. L. Alzate-Vargas, S. M. Blau, E. W. C. Spotte-Smith, S. Allu, K. A. Persson and J.-L. Fattebert, J. Phys. Chem. C, 2021, 125, 18588–18596.
83. Z. Zhang, Z.-W. Peng, M.-F. Hao and J.-G. Gao, Synlett, 2010, 2010, 2895–2898.
84. M. J. 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, H. P. Hratchian, J. V. Ortiz, A. F. Izmaylov, J. L. Sonnenberg, D. Williams-Young, F. Ding, F. Lipparini, F. Egidi, J. Goings, B. Peng, A. Petrone, T. Henderson, D. Ranasinghe, V. G. Zakrzewski, J. Gao, N. Rega, G. Zheng, W. Liang, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, K. Throssell, J. A. Montgomery Jr, J. E. Peralta, F. Ogliaro, M. J. Bearpark, J. J. Heyd, E. N. Brothers, K. N. Kudin, V. N. Staroverov, T. A. Keith, R. Kobayashi, J. Normand, K. Raghavachari, A. P. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, J. M. Millam, M. Klene, C. Adamo, R. Cammi, J. W. Ochterski, R. L. Martin, K. Morokuma, O. Farkas, J. B. Foresman and D. J. Fox, Gaussian Inc., Wallingford CT, 2016.
85. Y. Zhao and D. G. Truhlar, Theor. Chem. Acc., 2008, 120, 215–241.
86. M. Linder and T. Brinck, Phys. Chem. Chem. Phys., 2013, 15, 5108–5114.
87. J. Kemppainen, J. R. Gissinger, S. Gowtham and G. M. Odegard, J. Chem. Inf. Model., 2024, 64, 5108–5126.
88. H. Heinz, T.-J. Lin, R. Kishore Mishra and F. S. Emami, Langmuir, 2013, 29, 1754–1765.
89. L. Fang, S. Korres, W. A. Lamberti, M. N. Webster and R. W. Carpick, Faraday Discuss., 2023, 241, 394–412.
90. G. S. Kochhar and N. J. Mosey, Sci. Rep., 2016, 6, 23059.
91. P. M. Morse, Phys. Rev., 1929, 34, 57.
92. F.-G. Klärner and F. Wurche, J. Prakt. Chem., 2000, 342, 609–636.
93. C. Bolm and J. G. Hernández, Angew. Chem., Int. Ed., 2019, 58, 3285–3299.
94. A. Stukowski, Model. Simul. Mater. Sci. Eng., 2010, 18, 015012.
95. J. Zhang and H. Spikes, Tribol. Lett., 2016, 63, 24.

---