Many-Body Mechanochemistry: Intra-molecular Strain in Condensed Matter Chemistry

Mechanical forces acting on atoms or molecular groups can alter chemical kinetics and decomposition paths. So called mechanochemistry has been proposed to influence a variety of processes, from the formation of prebiotic compounds during planetary collisions to the shockinduced initiation of explosives. It has also been harnessed in various engineering applications such as mechanophores and ball milling in industrial applications. Experimental and computational tools designed to characterize the effect of mechanics on chemistry have focused exclusively on simple linear forces between pairs of atoms or molecular groups. However, the mechanical loading in condensed matter systems is significantly more complex and involves many-body deformations. Therefore, we propose a methodology to characterize the effect of many-body intra-molecular strains on decomposition kinetics and reaction pathways. We combine four-body external potentials with reactive molecular dynamics and show that many body strains that mimic those observed in condensed matter encourage bond rupture in a spiropyran mechanophore and accelerate thermal decomposition of condensed TATB, an energetic material. The approach is generalizable to a variety of systems and can be used in conjunction with ab initio molecular dynamics, and the two examples utilized here illustrates both the versatility of the method and the importance of many-body mechanochemistry. * Corresponding author: strachan@purdue.edu


Introduction
Mechanochemistry, the use of mechanical loads to trigger or influence chemistry, can enable otherwise forbidden reactions 1 , strengthen aerogels without limiting their elongation 2 , and prompt mechano-chromic/catalytic/fluorescent reactions in mechanophores 3,4 . Mechanochemistry is also believed to play a role when materials are subjected to extreme dynamical loads such as during planetary collisions 5,6 and shock initiation of explosives [7][8][9] . Over the last few decades, significant progress has occurred towards an understanding of covalent mechanochemistry in organic and molecular solids 10 . Experiments with exquisite resolution have revealed the strong effect of mechanical forces on chemistry. For example, atomic force microscopy (AFM) determined the force-extension response of individual covalent bonds and opened the door to a variety of experimental and computational studies on the effects of applied forces to chemical reactions 11 . Single molecule force clamp spectroscopy studies of reduction reactions of disulfide bonds showed an exponential dependence of the reaction rate on the applied force 12 . Numerous experimental and theoretical studies have explored and characterized mechanoluminescence, mechanocatalysis, and mechanochromic responses. 3,4,[13][14][15] Electronic structure calculations have shown that elongation forces on mechanophores can induce isomerization and ring-opening reactions, 3,4,16,17 the reduction in strength and distribution of stresses in knotted/entangled polymers, 18-21 explored the mechanical stresses in protein folding, 22 and assessed the rupture force and electronic properties for extension and scission of metal-organic junctions. [23][24][25] These types of studies use specialpurpose methods to characterize the effect of pair-wise external forces on chemistry including steered molecular dynamics (SMD) 26-30 , constrained geometries simulate external force (CoGEF) 31 and explicit force methods such as FMPES 17 (force-modified potential energy surface), EFEI 32 (external force explicitly included), and EGO 33 (enforced geometry optimization).
To date, computational and experimental studies designed to single out the effects of mechanochemistry have focused on elongation (two-body forces), while in most cases of extemporaneous mechanochemistry in condensed matter, such as shock-induced chemistry or the deformation of polymers, the mechanical loads are complex and many-bodied. 3,10,34,35 The importance of many-bodied mechanochemistry (MB-MC) effects has been revealed by several studies that tracked the reaction kinetics and pathways for organic systems under uniform shear loads 5,36,37 . These studies have demonstrated that dynamic mechanical loading can reduce activation energies and speed up chemistry. However, the results lack specificity and the individual deformations responsible for the observed, overall changes in kinetics and reactions paths cannot be directly determined.
In this paper, we propose a new computational approach to characterize MB-MC designed to capture the conditions experienced in condensed matter via the application four-body external potentials that affects torsional (proper dihedrals) and out-of-plane deformations (improper dihedrals). Our approach demonstrates the importance of many-body deformations on the activation energy of mechanophore isomerization (spiropyran) and the thermal decomposition of 1,3,5-trinitro-2,4,6-triaminobenzene (TATB). In the first example, we couple a four-body deformation around the spiro atom and apply the widely used CoGEF 31 method to load the mechanophore. We show that straining the central torsional angle by ~50 reduces the energy barrier associated with the mechanical activation of the mechanophore by a factor of two. The second example couples external many-body potentials that mimic the molecular strains observed under shock loading 35 using reactive molecular dynamics in the molecular crystal TATB. A Kissinger analysis of the MD simulations shows marked changes in chemical kinetics and decomposition paths that can explain previous inconsistencies between mechanical and thermal insults in this class of materials 8,9,38 . Importantly, this second example cannot be tackled at all with any prior approaches.

Methods
Computational mechanochemistry studies use various techniques to deform systems of interest and explore the evolution of the system as an external force or displacement is increased. For example, SMD anchors atoms to an evolvable position, with restoring forces on the anchor site, using a simple potential, typically harmonic. CoGEF performs geometry optimizations while constraining the desired atoms. Evolving these constraints (e.g. increasing linear distance between two atoms) can simulate molecular deformation to incite mechanochemistry. Such approaches cannot readily capture complex many body molecular distortions observed in condensed matter under mechanical loads 9,35,39 . To address this gap, we develop Many-Bodied Steered Molecular Dynamics (MB-SMD) in which we use four-body external biasing potential to deform molecular groups of interest along torsional and out-of-plane degrees of freedom. Four-body terms are motivated by our observation of the large molecular distortions observed under shock loading 35 as well as their general applicability to explore the conformational space of a molecule. This external biasing potential is implemented in the form of four-body harmonic inter-atomic potential, = ( − ) 2 (1) where is either a dihedral angle measuring torsions around the central jk bond or an improper dihedral angle measuring out of plane distortions, see Fig. 1. MB-SMD implementation details, as well as statistical characterizations of deformed molecules, are available in Sections S1 and S2 of the supplementary information. We apply torsional and out-of-plane deformations in independent simulations for a variety of values for the spring constant ( ) and target angle ( ) to parametrically study the effects of each deformation on reaction kinetics and first step pathways. These potentials will balance with the internal energy expression of the system of interest and will result in equilibrium configurations than can be controlled via the external potential's spring constant and target angle, as shown in S3.
All simulations use all-atom representations and molecular dynamics (MD) simulations were performed using LAMMPS 40 with a 0.025 ps timestep under NVE dynamics. The ReaxFF 41 reactive force field was used and partial atomic charges were calculated using the charge equilibration scheme 42 at every step in the simulations with a tolerance of 1x10 -6 .

Many-body distortions of spiropyran
For mechanophore simulations, a gas phase spiropyran molecule, Figure 1a, is employed to act as a truncated mechanophore/polymer system. Atomic forces were calculated with the ReaxFF-LG 43 parametrization. We employ the CoGEF 31 algorithm, with anchored H atoms (cyan in Figure  1a), using a constrained structural relaxation via energy minimization. Geometry optimization was conducted via the conjugate gradient method 44 with energy and force tolerances of 1 x 10 -4 (unitless) and 1 x 10 -6 kcal/molA, respectively. To better explore the underlying energy landscape and being trapped in a local minimum, we combined CoGEF with MD, performing 10 fs of NVE dynamics between geometry optimizations.
To explore the effect of many-body deformations on the activation of the mechanophore, we include a four-body potential to control the dihedral angles around the spiro atom and one the C atoms in the six-member ring, yellow highlight in Fig. 1a. The equilibrium dihedral angle predicted by ReaxFF is 72. A spring constant of 50 kcal/mol is used and the target angle is varied from 70 to 25 across independent CoGEF simulations. We determine the energy barrier for isomerization as the difference in total energy (not including the biasing potential) between the initial configuration and that under which the opening of the C-O bond is observed. Thermalizations (20 ps NVT) were conducted with the external potential prior to CoGEF runs to "pre-deform" the mechanophore to the desired angle of the external potential and keep the initial state at ~300 K.

Many body distortions of TATB
For TATB, see Figure 1b, we used a simulation cell with 250 molecules (6000 atoms) in the triclinic crystalline lattice setting of Cady and Larson 45 . Atomic interactions were calculated using the ReaxFF-2018 46 parametrization that has previously been used to study the TATB system for shock loading and thermal decomposition 47,48 . Thermal decomposition under the presence of the external potential is simulated by linearly heating the system from 1000 K to 3000 K at various rates (from 4 K/ps to 40 K/ps) 49 .
We applied two external potentials, in two independent sets of simulations: i) proper dihedrals of C-C-N-O atoms that controls the torsional rotation around the CN bond of the nitro group (atoms are highlighted in yellow on Fig. 1 (b)), and ii) improper dihedrals of C-C-C-N (nitro group only, atoms are highlighted in yellow with a red outlines) that bends the nitro groups away from the plane of the ring, an 'out-of-plane' deformation. Both TATB deformations (which are studied one at a time) are applied simultaneously to all three nitro groups in all 250 molecules in a crystalline lattice with fully periodic boundaries. For both cases, the system is equilibrated for 20 ps at 300 K under the presence of the external potential to "pre-deform" the system prior heating and decomposition. We utilize multiple spring constants for each set angle of each respective deformation to assess to influence of both external potential parameters. Spring constant effects, which are minimal, are displayed in Section S3 of the SI.

Role of internal rotation on Spiropyran isomerization
Extension of the spiropyran molecule, Figure 1(a), can induce the scission of the C-O bond at the spiro atom through an isomerization reaction that results in a chromatic response, as shown in Ref. 4. As was done in prior studies of mechanophores 4,17 , we performed CoGEF simulations to characterize this reaction but, in addition to the constraints applied to the anchor H atoms, we apply MB-SMD around a central bond marked with yellow in Fig. 1a. The associated dihedral has an equilibrium value of 72 and ten independent CoGEF simulations were performed for target torsional angles ranging from 70 to 25. Figure 2(a) shows the activation barrier (increase in ReaxFF internal energy from the initial configuration to isomerization) as a function of target torsion angle. The simulations reveal a significant decrease in the energy barrier needed to activate the isomerization reaction caused by the many-body molecular strain. A deformation of 30 in the selected dihedral angle reduces the activation energy by 1 eV, or one third. Renderings of the reaction pathway are available in S4 in the SI. The central rotation studied here does not greatly affect the reaction path or product state of the reaction, which an understanding of is critical to the engineering of mechanophores to perform a variety of capabilities.
Importantly, the absolute energy required for the isomerization reaction to occur, referenced by the system with no MB-SMD, decreases with the application of the many-body potential, see Fig. 2b. This indicates that the many-bodied deformation not only includes additional strain energy into the system, but also lowers the energy threshold needed. The decrease in activation barrier can be separated into the decrease in absolute energy threshold and the initial energy state increase due to deformation. Roughly two thirds of the activation barrier decrease from lowering the absolute energy needed, and one third of the decrease is from the torsional strain energy added to the system. Thus, many-body strains can not only reduce the energy required to activate mechanochemistry, but also drive it through a new lower energy path. Lastly, torsional deformations resulting in dihedral angles approaching 20 (change in >50) result in a small percentage of systems undergoing isomerization without applying a linear strain, at 300 K.

Role of many-body strains in the thermal decomposition of TATB
When materials experience dynamic mechanical loading, such as during a high velocity impact, the resulting shockwave can generate a plethora of ultrafast thermal, mechanical, and chemical processes, 50 including phase transitions [51][52][53] , plasticity and failure [54][55][56] , and chemical reactions 6,57,58 . These processes are accelerated by energy localization into so-called hotspots. There is increasing evidence, from computational studies, that hotspots formed under shock loading react at rates significantly faster than equivalent ones taken to the same temperature and pressure states, but under equilibrium conditions 8,9 . We have recently shown that a possible origin of this increased reactivity is the presence of highly strained and deformed molecules in the hotspot. 35,39 This energy localization, in the form of intra-molecular potential energy, is persistent well into reaction timescales. The increase in strain energy in TATB is caused primarily by torsional and out-of-plane deformations of the nitro and amino groups, shown in Section S5 of the SI, and incurs a local amorphization of the material that affects kinetics 59 and heat transport 60 .
To assess the influence of these extemporaneous deformations on reactivity, we applied MB-SMD as four-body biasing potentials to all nitro groups in TATB molecules packed into its crystal structure with two sets of simulations: torsional deformations which rotate the nitro groups with respect to the carbon ring, and out-of-plane deformations which bends nitro groups upwards/inwards. To study the effect of these molecular strains, we performed a series of reactive MD simulations where we increase the temperature at various rates, varying stiffnesses and set angles of the external potential, until significant reactions are observed. Quite interestingly, the simulation results follow Kissinger 49 statistics quite closely, even for highly deformed systems. The characteristic reaction time, td, defined as the point of maximum internal energy, follows the expected rate dependence equation, see S6 of the SI. Figure 3a shows the resulting activation energy associated with decomposition as a function of the average value of the angle (dihedral or out of plane) obtained via the different set angles and spring constants. The deformation angle used in Figure 3a is the measured deformation level at 300 K, prior to heating-induced chemical reactions. The green star represents the value for a system with no applied external potential. Both external strains result is a reduction of the activation energy for decomposition, which explains the increased reactivity in energetic materials under a dynamical loading that deforms molecules. We find a significantly larger slope for the nitro outof-plane bending deformation, as compared to the torsional deformation. It is interesting to note that the application of MB-SMD using angle values close to the equilibrium values of zero for high spring constants results in an increase in the activation energy, as the molecules become more rigid. To further assess the MB-MC potency of each deformation, Figure 3b shows the decrease in activation energy, referenced from the undeformed case, as a function of the increase in internal energy caused by MB-SMD at 300 K (average energy of the deformed system without counting the energy pentalty of the external biasing potential). Clearly, the out-of-plane deformation of the nitro groups has a stronger effect on facilitating the decomposition of TATB. In S7 of the SI, we fit and compare these results to a simple mechanochemical kinetics model 61 . Lastly, we turn our attention whether MB-MC affects the decomposition path, specifically, the initial step. Several initiation mechanisms have been proposed and studied for TATB. Intramolecular hydrogen transfer from the nitro to amino groups is the primary reaction pathway, it occurs in nearly 100% of initial reactions under thermal 62 insults and the majority of initial reactions for shock insults 48 . Inter-molecular hydrogen transfer is the main alternative reaction pathway for undeformed TATB 62 and the study of chemistry on shear bands in TATB showed an increase of nitro scission reactions 9 .
Our simulation reveals that both molecular strains studied have a strong effect on decomposition paths, see Figure 4. While intra-molecular H transfers dominate initiation for no or small strains, torsional deformations around the C-N (nitro group) bond of over 40 opens a second decomposition path: bi-molecular hydrogen transfer. This bi-molecular path becomes dominant for angles above 70. The inset of Fig. 4a exemplifies this inter-molecular reaction path between adjacent amino and nitro groups in neighboring TATB molecules in the crystal. We find it interesting that for angles below the 40 deformation, amino groups increase their level of deformation as nitro groups deform. Above 40 the amino groups begin to significantly un-deform with increasing nitro deformation, i.e. they find that breaking their hydrogen bonds with the nitro groups and relaxing is energetically preferable to deforming alongside the nitro groups, as shown in S9 of the SI, which corresponds to the onset of the alternate reaction pathway.
For the out-of-plane bending, shown in Figure 4b, there is also a change in mechanism for strong enough many-body potentials. Surprisingly, this deformation opens up a different alternate reaction path: nitro group scission, which is higher in activation energy than both the intra-and inter-molecular H-transfer reactions 62 , and does not activate the inter-molecular H transfer at all. The availability of the new path coincides with a conformation change in the TATB. As described in detail in Section S8 of the SI, the inner C ring buckles and the nitro and amino groups may more freely rotate. This structural change causes the change in slope of activation energy vs. deformation shown in Figure 3 and Section S7 of the SI. The inset in Figure 4b shows instances of an intramolecular h-transfer (bottom arrow) and a nitro scission (top arrow). Hence, we find that not only does deforming a different DoF in the system result in different levels of activation barrier reduction, but it also opens different pathways for decomposition reactions to take place.

Summary
We developed a new method, Many-Bodied Steered Molecular Dynamics, to study the effect of four-body deformations on the chemistry in condensed, molecular systems. This approach captures the complex mechanochemical states experienced in condensed matter chemistry, including shock induced chemical initiation in high explosives. Our results indicate that that manybody strains play an important, sometimes dominant, role in covalent mechanochemistry.
In the isomerization reaction of the mechanophore spiropyran, torsional deformations around the central spiro atom led to a decrease in the activation barrier by 45% for rotations of 45. For crystalline TATB, MB-SMD strains designed to mimic the deformations found during shock loading also led to a significant decrease in activation energy and changes in the decomposition path. These results explain the increasing evidence of accelerated chemistry in highly deformed areas in shocked explosives and MB-SMD will likely be vital method for developing a predictive understanding of their initiation.
These two examples illustrate the versatility of the MB-SMD method and indicates the pervasiveness of many-body mechanochemistry in condensed matter. Our results show that manybody strains can reduce the energy required to activate mechanochemistry, thermal-induced chemistry, and drive chemical reactions.