Molecular Origin of Mechanical Sensitivity of the Reaction Rate in

Subscriber access provided by CORNELL UNIVERSITY LIBRARY

Article

Molecular Origin of Mechanical Sensitivity of the Reaction Rate in Anthracene Cyclophane Isomerization Reveals Structural Motifs for Rational Design of Mechanophores Nikolay V. Plotnikov, and Todd J. Martínez J. Phys. Chem. C, Just Accepted Manuscript • DOI: 10.1021/acs.jpcc.6b04924 • Publication Date (Web): 26 Jul 2016 Downloaded from http://pubs.acs.org on July 28, 2016

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry C is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 24

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Molecular Origin of Mechanical Sensitivity of the Reaction Rate in Anthracene Cyclophane Isomerization Reveals Structural Motifs for Rational Design of Mechanophores Nikolay V. Plotnikov†* and Todd J. Martinez* Department of Chemistry and the PULSE Institute, Stanford University, 333 Campus Drive, Stanford, CA 94305 and SLAC National Accelerator Laboratory, 2575 Sand Hill Road, Menlo Park, CA 94025 Abstract

The observed pressure sensitivity of the isomerization reaction rate of bis-anthracene cyclophane photoisomer has attracted significant attention in the rational design of mechanically sensitive materials. However, the molecular origin of this sensitivity remains unclear. We developed an ab initio molecular model to quantify the effect of pressure on the reaction rate and to elucidate its molecular origin. Pressure-induced deformations and changes along the reaction free energy surfaces are estimated from ab initio molecular dynamics trajectories. Our model predicts a barrier reduction of ~2 kcal/mol at 0.9 GPa (in agreement with experiment). The barrier reduction is linear in the low pressure regime (up to 2GPa), but has a nonlinear dependence at higher pressures. We find that pressure alters the reaction path and that the mechanical sensitivity of the reaction rate is caused by an uneven distribution of the free energy increase along the reaction surface. The uneven distribution primarily results from destabilization of the reactant, which has a lower mechanical rigidity along a particular deformation mode (flattening of anthracene rings as they are pushed towards each other in the cyclophane framework). Exploiting similar structural motifs to maximize rigidity differences along the reaction coordinate represents a promising rational design strategy.

Introduction The dependence of chemical reaction rates on mechanical stimuli such as pressure,1 tensile stress2-4 and force has been of great interest in chemistry.5-12 A full understanding of the coupling between chemical reactivity and mechanical stress would open the door to ab initio prediction of material mechanical properties and development of novel materials with enhanced sensing, self-healing, shock absorbing and energy conversion properties. The properties of materials are a function of their molecular structure. Thus, it is valuable to identify

structural motifs that grant materials the desired properties. Mechanically sensitive molecules, known as mechanophores, might have some common structural motifs underlying their mechanical sensitivity. Once identified, libraries of these structural motifs can be used in the rational design of new mechanophores. Notable examples of the identification of such structural motifs in mechanophores are two recent studies by the groups of Bardeen and Chronister.13,14 The anthracene cyclophane photodimer (Figure 1) represent an example of a successful molecular design.13 Unfortunately, the subsequent design

Plotnikov and Martínez – Page 1 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

attempt did not show high mechanical sensitivity.14 This indicates a failure of the employed design strategy, based on a perturbative picture of the mechanochemical coupling.1,3 Such a strategy seeks to minimize the reaction activation volume, while assuming it is independent of pressure. The activation volume estimated from crystallographic data for reactants and products or from electronic structure calculations in vacuum15 might fail to provide accurate predictions when the activation volume depends on the applied pressure. Obtaining reliable experimental estimates for pressure-induced barrier reduction can be quite challenging. For the bis-anthracene isomerization, two employed techniques provide significantly different estimates.13 While the Arrhenius plots (lnk vs 1/T) at normal pressure and at 0.9 GPa provided an estimate of 20.8 kcal/mol, the rate constant dependence on pressure (lnk vs p) gave an estimate of ~ 2 kcal/mol (see Table S1 for details). Another question from the experiment, which remains unanswered, is the nature of two molecular populations (one sensitive to the applied pressure and a second which is insensitive). To quantify the reaction rate pressuredependence and to elucidate its molecular origin, we compute Gibbs free energy surfaces from ab initio molecular dynamics over a range of pressures. Next, we extract the entropic, mechanical and internal energy contributions to the free energy along the reaction path. Finally, we establish correlations between pressure-induced energetic and structural changes to identify the structural motifs for the rational design of mechanophores. The paper is organized as follows. In the “Methods” section, we summarize the key components of our computational algorithm. We show how it is used to establish the effect of pressure on: 1) energetics along the reaction path; 2) topological changes of the reaction free energy surfaces (that is, the reaction mechanism); and, 3) structures of critical points along the reaction path. In the “Results” and “Discussion”, we report and compare our theoretical predictions with experimental data and establish the relationship between the observed energetic and structural changes. In the “Conclusions” section, we summarize the key findings of this study, and discuss application of molecular modeling in the rational design of mechanophores.

Methods Potential energy surface description with density functional theory The reaction potential energy surface was computed by constrained energy minimization with QChem16 along two interatomic distances, which also correspond to the breaking C-C chemical bonds in the reaction (denoted by R1 and R2 in Figure 1). The classical vibrational entropy at 300K and zero point energy in vacuum were estimated from harmonic analysis for the reactant (RS), product (PS), and transition state (TS) geometries (details in SI). The TS geometry was confirmed to have a single negative frequency. The minimum energy reaction path was computed using the Nudged Elastic Band (NEB) method17 as implemented in TeraChem.18 We used spin-unrestricted density functional theory with exact exchange to describe the reaction potential energy surface. Unless otherwise stated, calculations presented in the main text were carried out using the spin-unrestricted B3LYP19 functional with the 6-31G basis set and an empirical dispersion correction.20 We compared B3LYP results for the activation energy, heat of reaction, and structural parameters with various other functionals and experimentally available data (Table S2). We note that in this work we focus on differences in the activation energy. The differences between the activation energies at different pressures are expected to be more accurate than absolute activation energies due to error cancellation. Computing the free energy surfaces The effect of stress on the reaction rate constant is estimated through the activation free energy, which is related to the rate constant by transition state theory:21 (1)

k=

κ exp ( − β ∆g ‡ ) βh

where κ is the transmission coefficient; β=1/kBT; h is Planck’s constant; kB is Boltzmann’s constant; T is the absolute temperature; and ∆g‡ is the difference between the Gibbs free energy function along the reaction coordinate, ξ, at the TS and at the RS. While the transmission coefficient depends on the choice of the reaction coordinate, we assume that it is not affected by the applied mechanical stress. Under this assumption, the ratio of two rate constants for the same reaction at different mechanical conditions is determined by the difference in computed free energy barriers from Potentials of Mean Force (PMF).


Page 2 of 24

Page 3 of 24

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The free energy function (or PMF) is computed from the probability distribution function, ρ(ξ ξ), obtained by sampling with a thermodynamic potential H=U+PV:

g p ( ξ ) = −β −1 ln ρ ( ξ )

(2) where

...

H

H

designates the configurational average

obtained with the potential energy U at pressure p. To enhance the configurational sampling of highenergy regions of the reaction free energy surface, we introduce a set of potentials biased to particular sections of the reaction path: (3)

Tm = H + wm ( ξ m0 )

Here

wm ( ξ 0m ) is a harmonic potential, which

restrains components of the reaction coordinate vector to the R1 and R2 values, determined from the NEB path in vacuum. That is, for the m-th bead along the NEB path, it is given by

( )

(

(4) wm ξ 0m = K R1 − R1,0 m

)

2

+ K ( R2 − R2,0 m )

2

To remove the bias from the probability distribution we used a combination22 of the free energy perturbation23,24 and umbrella sampling25 methods:

(5)

Where

ρ (ξ )

H

=

δ ( ξ′ − ξ ) exp [ β wm ] T exp [ β wm ]

m

the 1D representation using the R1+R2 value (which is a projection on the concerted path, corresponding to the simultaneous rupture of the two bonds). The path defined this way implies the approximation that the reaction path is not shifted by stress (in terms of R1 and R2 values) along the dividing surface or at least adequately represents the minimum free energy path within an acceptable error. To quantify this assumption, we computed 2D free energy surfaces in a range of pressures. Finally, we analyzed the curvature of the transition state region in the direction defined by R1-R2, which is perpendicular to the reaction path. Projection of the free energy surface on this coordinate illustrates the effect of pressure on the difference between concerted and sequential reaction paths in this isomerization. The high computational cost of the required configurational sampling is circumvented by solving the most computationally expensive parts of the SCF procedure on general-purpose NVidia graphical processing units using the development version of TeraChem,18 which is interfaced with our implementation of enhanced sampling MD in an ideal gas pressure bath (see below).

Computing the internal energy, entropic and mechanical contributions to the free energy profile To construct the internal energy profile (for a given temperature and pressure) we compute the ensemble average of the SCF energy:

Tm

δ ( ξ′ − ξ ) is the delta function of a vector

of the reaction coordinate. The denominator of Eq. (5) corresponds to the free energy penalty of introducing the bias. This penalty was estimated with several free energy estimators: the mid-point exponential average,26 the average of the forward and backward exponential averages, the linearresponse approximation27 and the acceptance ratio.28 Estimates of Eq. (5) for different biasing points were combined with the weighted average as described previously,29 and with the weighted histogram analysis method30,31 (WHAM), which provided another independent estimate. The deviation between different free energy estimators was used to evaluate the error in the free energy calculations. To define the reaction free energy path, we used R1 and R2 values from the gas-phase NEB path. In molecular dynamics these values were constrained with harmonic potentials, given in Eq. (4). We found it convenient to define the reaction path in

(6)

U ξ =ξ 0 = E AI m

H

=

E AI exp  β wm  exp  β wm 

Tm

Tm

where EAI is the ab initio total energy in a given configuration (the self-consistent field energy obtained from the electronic structure calculation). With eq. (6) we estimate the internal energy profile at a given pressure and the internal energy increase for every point along the reaction coordinate relative to energies of equilibrium geometries in vacuum. Thus, we position the internal energy profiles at different pressures relative to each other. The difference between the free energy profile and the internal energy profile provides an estimate for the sum of entropic and mechanical contributions to the free energy profiles: (7)

∆G − ∆U = P∆V − T ∆S

Note that these differences are computed relative to the relevant values for the reactant state. To further separate the entropic and mechanical terms we evaluate the entropic contribution for



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

each MD trajectory using the quasi-harmonic approximation.32-34 Thus, we compute a profile along the reaction coordinate for the entropic contribution to the free energy (-T∆S). This let us finally separate out the mechanical contribution to the free energy profile (P∆V), and compute the corresponding profile for the volume change at a given pressure. In other words, we estimate the difference in the volume reduction relative to the volume reduction of the corresponding reactant state. Finally, we identify which molecular dimensions are affected by pressure along the reaction path. For this we computed the average geometry for RS, TS and PS from propagated MD trajectories.

Ab initio molecular dynamics in the ideal gas barostat The effect of pressure in ab initio MD was simulated using our adaptation of the statistical ideal gas pressure bath35,36 which represents the effect of mechanical stress (hydrostatic pressure). The simulation system, treated quantum mechanically (QM region), is immersed in the reservoir of particles (the detailed description of our implementation is available in SI). The ideal gas particles obey ideal gas statistics (i.e. particle fluctuation in a volume element, the equation of state and the velocity distribution) and do not interact with each other but do interact with the QM region via the repulsive term of the LennardJones potential, thus exerting pressure on the simulated system. In this work, we explicitly capture the dependence of all degrees of freedom of a molecule on the applied mechanical stress by performing molecular dynamics simulations in the pressurized simulation box. Thus, there is no need to postulate a priori how the mechanical stress is applied (i.e. by exerting force on selected atoms). While in some cases, such as bond rupture in polymers, the atoms to which a mechanical stress (e.g. force in the atomic force microscopy) is applied can be relatively easy identified,7,37 in the case of applied pressure, it is not obvious which atoms experience mechanical stress and to what degree. Therefore, formulation in terms of force-modified potential energy surfaces38 becomes non-trivial. An extension of the forcemodified potential energy surface has been recently proposed39 which creates “pseudo-hydrostatic pressure” by applying force directed to the molecular centroid on all atoms. Here we use ab initio MD, (Born-Oppenheimer MD) in which the electronic density is relaxed for

every nuclear configuration (at each MD time step) with the usual SCF procedure.18

Results Effect of pressure on the reaction rate Our calculated changes of the activation free energy (Table 1) agree well with the experimental data13 for lnk vs. p (Table S1). The activation volume, ∆V‡, is estimated to be -16 Å3 in the range of 0-1.5 GPa (see inset in Figure 2). The sensitivity of the rate constant to pressure in the original experiment was limited by compressibility of the polycycloolefin resin (ZEONEX Z480, with glass transition temperature Tg = 138−140 °C). Our calculations performed in the ideal gas barostat predict a continuing strong dependence of the free energy barrier up to ~3 GPa (Figure 2). This dependence starts decreasing significantly above 3 GPa. Thus, over a wide range of pressures (0-7.3 GPa) our calculations predict a nonlinear pressure dependence of the reaction rate: a strong (nearly linear) dependence from 0-3GPa and a weaker but still notable dependence from 3-7.3 GPa. The barrier is reduced by 11.6 kcal/mol at the maximum simulated pressure of 7.3 GPa. The activation volume of -16 Å3 corresponds to reduction of the activation energy by ~ 2 kcal/mol at 0.9 GPa. This estimate is significantly smaller than 20.8 kcal/mol, which was an experimental estimate obtained from the temperature dependence of rates13 at ambient pressure and at 0.9 GPa. A narrow temperature range 295~320 K used to construct the Arrhenius plot at 0.9 GPa might have resulted in a significant extrapolation error. Our calculations also predict a strong dependence of the reaction exothermicity on pressure (Table 1). The corresponding change in the reaction free energy is always larger than in the activation energy (Figure 3). Unfortunately, no experimental calorimetry data are available for comparison. The negative activation volume13,15 was proposed to be the driving factor in lowering the activation barrier (by p∆V‡). However, the difference in volumes of the transition state and the reactant state in vacuum was estimated15 to be -60/-70 Å3. The energy equivalent for this activation volume is ~9 kcal/mol at 0.9 GPa and ~24 kcal/mol at 2.6 GPa. Therefore, the activation volume corresponds to the difference between pressure-induced change in volumes of the reactants and of the transition state. This difference is due to different mechanical compliance of two states.40 In other words, the activation volume is the excessive volume reduction of the RS compared to TS, that is:


Page 4 of 24

Page 5 of 24

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


∆V ‡ = (V ‡ ( p ) − V ‡ ( p = 0 ) )

(7)

− (V RS ( p ) − V RS ( p = 0 ) )

Here the volume reduction is defined relative to vacuum. All states along the reaction path are compressed, but at low pressures the energies of each individual state are linear in the applied stress. While there are various definitions of molecular volume,41 we compute volume as the extensive conjugate variable of pressure. That is, changes in energy have a thermodynamic (energy) equivalent to the change in volume. Notably, even an applied pressure of 7.3 GPa is insufficient to cause a spontaneous reaction (the activation energy is 5.5 kcal/mol). The rate limiting step remained the rupture of the first C-C bond in the cyclobutane ring. Thus, our findings do not support the previously proposed two-step mechanism13 where flattening of the anthracene rings is followed by breaking C-C bonds. Actually, flattening of the anthracene rings from the equilibrium value in vacuum to that at 7.3 GPa would require significantly more energy than the original barrier: ~60 kcal/mol of energy (Figure S8) on top of the reduced 5.5 kcal/mol barrier. Absolute values of the reaction barriers with B3LYP functional are systematically underestimated (which is a well-known self-interaction problem in DFT-functionals42) and the reaction heat is overestimated. However, agreement between experimental data for the rate constants and calculated pressure-induced barrier reduction indicates error cancellation in estimating the relative values (see SI). Since we are interested in the pressure-induced difference in the free energy values, the error cancellation significantly improves the accuracy of our theoretically predicted pressureinduced free energy differences. This statement however, does not universally apply to a wide range of mechanical perturbation, e.g. when the unperturbed barrier is smaller than the experimentally observed reduction. Furthermore, as the reaction barrier gets smaller and the TS curvature significantly decreases, the assumption of a constant transmission factor and the neglected probability of tunneling43 might result in a higher error in the rate constant estimate.

Pressure-induced changes in topology of the reaction surface Pressure significantly alters the topology of the reaction free energy surface: it affects both the relative position of the final states and curvature of the TS region (Figure 4). The concerted and

sequential reaction paths are shown with blue and black dots, connected with a solid line, in the top right panel in Figure 4. Notably, the difference between the reaction barriers of the concerted and of the sequential reaction paths decreases with increasing pressure (see the top left panel in Figure 4 and Figure S7). An applied pressure of 7.3 GPa reduces the difference in the activation energy between two mechanisms by ~6 kcal/mol compared to vacuum (Table S2). Thus, high pressure alters the reaction path and, at least in principle, can even alter the reaction mechanism. In this respect, the effect of stress on the reaction surface is similar to the effect of solvent1 (e.g. on the nucleophilic substitution mechanisms).

Pressure-induced structural changes Changes in the reaction free energy profile correlate with pronounced structural changes, which the mechanophore undergoes along the reaction path at different pressures (Figure 5). Average structures for reactant, product and transition states (see Figure 5) were computed from superposition of corresponding MD trajectories.34,44 We estimated the change in volume along the reaction path relative to the volume reduction of the RS (eq. (7)). This profile is the mechanical equivalent of the corresponding contribution to the free energy change, and it shows the relative mechanical stability of the mechanophores at different points of the reaction path. The RS is the least stable with respect to the applied pressure. We note that pressure significantly deforms RS, TS and PS compared to the vacuum equilibrium geometries (with the reactant state being most affected, as can be seen in Figure 5). Internal energy and entropic contribution The activation energy reduction in the pressures range 0-2.0 GPa is not accompanied by notable changes in the internal energy contribution (Figure S8). However, at higher pressures the internal energy of the RS increases faster than for TS and PS, thus contributing to lowering of the reaction barrier and higher exothermicity (Table 1). Destabilization of the RS reaches a maximum of ~ 5 kcal/mol at 2.7 GPa (Figure S8). No significant decreases in the activation entropy were detected at elevated pressure, even at 7.3 GPa (Figure 6 and Figure S14). The entropic contribution to the free energy barrier at 300K in vacuum was estimated at +0.3 kcal/mol (-0.4 kcal/mol to the reaction free energy). Thus, the entropic factor is not responsible for the observed reduction of the activation energy with pressure (within the error of 1 kcal/mol - see also Figure S14 - at 300 K, based on



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

the harmonic approximation in vacuum and quasiharmonic estimates at elevated pressures). The decrease in the activation energy at pressures below 3.0 GPa is caused by RS destabilization. The difference in the free energy increase between the RS and the PS also results in increased reaction exothermicity (see Table 1). Free energies of the TS and of the PS also increase, but at slower rates, since the RS is more compressible. The change in the activation free energy is due to the mechanical term in this pressure range. The internal energy and the entropic contributions are insignificant. The barrier reduction follows the change in reaction free energy much closer at pressures exceeding 3.0 GPa. Thus, the TS becomes increasingly more stable (more rigid) with respect to stress than the final states. This might indicate relative stabilization of the polyradical intermediate.

Discussion Strain in cyclobutane and aromaticity of polyradical intermediate Structural deformations increase strain in the cyclobutane rings of PI, which might result in an additional RS destabilization, and, counterintuitively, in a possible accelerated dissociation of the very same bonds. In PI, these C4 rings are less strained than in cyclobutane, and C-C bonds are elongated to 1.65 Å (compared to 1.55 Å in cyclobutane45). The strain energy in a cyclobutane ring45 is ~26 kcal/mol, while the difference between the activation energies of dissociation for PI and dianthracene is only about 13 kcal/mol,13 which translates to only ~ 6.5 kcal/mol of strain per ring (which is close to the strain energy in cyclopentane). Compression of R1 and R2 (or increasing planarity of rings) might cause additional strain in the RS and lower the activation barrier. While the TS contains one cyclobutane ring (Figure 1), the concerted TS and bisanthracene do not contain any strained cycles. The increase in strain, however, should correlate with the decrease in the internal energy of activation, which is observed in the intermediate pressure range (at pressures above 2.0 GPa). Flattening the anthracene rings in the TS might offset the increased strain since aromaticity is facilitated in a planar structure. Furthermore, increased planarity allows for a higher delocalization of the unpaired electrons over the anthracene rings, thus stabilizing the polyradical character. The resonance energy of benzene is about 36 kcal/mol while the resonance energy of anthracene is 83 kcal/mol. Thus, the total resonance

energy of RS (due to four benzene rings) is ~144 kcal/mol, while the total resonance energy of PS (2 anthracenes) is ~166 kcal/mol. Therefore, the resonance energy of the TS is between 144 and 166 kcal/mol by these estimates, and the energy gain (due to the TS stabilization) should be smaller than the difference between the final states of ~22 kcal/mol. The B3LYP energy of the diradical state is too high, and we only observe a plateau in the TS region in vacuum, while the semi-empirical PM3 method predicts a stable diradical intermediate (see Figure S13). Given this and the magnitude of the internal energy change by high pressure (~3 eV) one might ask about the effect of pressure on the energy gap between electronic states. That is, whether pressure increases or reduces the energy gap between the ground and the excited state to the point where a single-determinant wave function or even the Born-Oppenheimer approximation fail. This would depend on the equilibrium geometries of the corresponding states, and their sensitivity to the applied stress. We qualitatively considered the effect of placing an electron from the highest occupied (HOMO) to the lowest unoccupied (LUMO) molecular orbital (Figure S12, bottom panel). While both of the final states have an anti-bonding HOMO and bonding LUMO with respect to the π-π interactions, the HOMO of the TS is of bonding character with respect to the π-π bonding interactions between the aromatic rings while its LUMO is anti-bonding. This means that energy of HOMO will increase upon compression in the final states and decrease in the TS. This kind of interaction was observed in the excimer of benzene dimers and studied theoretically, resulting in a tighter binding in the excited state (8.1 kcal/mol) relative to the ground state (2.6 kcal/mol), see Ref. 46 and references therein. While the pressure-induced strain increases the RS energy, the flattening of the anthracene rings increases the planarity and aromaticity and, possibly, π-π bonding interactions in the TS (which has a “third character” due to a hypothetical polyradical intermediate). Any stabilizing contribution offsets the destabilizing effects caused by other stress-induced deformations, leading to a higher relative stability with respect to other states. Thus, this contribution might result in lowering the TS energy when Z distances are decreased (in Figure 1) due to a higher orbital overlap. This also suggests that the energy gap between the ground and the excited state will decrease for the final states, and increase for the TS as force along the Z-coordinate makes the TS more aromatic. The concerted TS gains more from increased planarity (compared to


Page 6 of 24

Page 7 of 24

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


the sequential one) since flattening of the anthracene rings increases strain in the C4 ring of the latter. This could explain the decreasing free energy difference between the sequential and the concerted paths at high pressures (Figure 4).

Mechanical contribution and work dissipation due to uneven compressibility The uneven increase of the Gibbs free energy along the reaction path at pressures below 2 GPa is primarily due to the mechanical term. Higher pressures also cause uneven increases of the internal energy. Pressures below 2 GPa cause significantly larger deformations of the RS (Figure 5), which result in higher compression of the RS (Figure 7). The free energy dependence is linear under low pressures. Higher pressures notably compress the PS as well. The free energy dependence in the high pressure regime becomes non-linear. High pressures induce substantial changes in the electronic density and structure, which result in an uneven increase of the internal energy along the reaction path. Thus, we can identify two mechanisms for barrier reduction. Initially, the mechanism is a predominant destabilization of the RS. The TS is destabilized less than the RS since it mechanically resembles the more stable PS (the “late” TS). In this regime, the change in the activation energy is smaller than the change in the reaction free energy. Higher pressures deform the PS and the TS as well, primarily by flattening the anthracene rings, what increases the resonance energy of the TS. In the high-pressure regime, the change in reaction free energy closer follows the activation energy change. Deformations at low pressures are rather elastic, and the free energy change is due to the difference in compressibility of the molecule along the reaction path. At high pressures deformations are inelastic; they trigger substantial changes in the electronic structure, and in the internal energy profile. Thermodynamic origin of mechanical sensitivity of chemical reaction rates The work, which is done on a system by external forces, raises its free energy. The reaction rate constant is mechanically sensitive if the free energy increase is unevenly distributed along the reaction path (or, more generally, across the free energy surface). If the free energy increases evenly along the reaction path, the reaction rate is not sensitive to the external force. These two scenarios are shown in Figure 8. Differences in relative stability of critical points on the reaction path (RS, TS, and PS) might result in two possible scenarios for the reaction barrier reduction:

a) RS is destabilized more (the free energy increase is higher) than PS or TS b) both final states are destabilized more than the TS. This destabilizing effect of mechanical stress is conceptually similar to the stabilizing effect of solvent (e.g. on SN1-SN2 reaction mechanisms). Thus, empirical arguments of physical organic chemistry47,48 are often used to rationalize correlation between the energetic and topological changes of the reaction surface. Case (a) requires that the TS is closer to the PS, which, in turn, is more stable than the RS. Case (b) requires that the TS has some other electronic character, due to an intermediate state, which is more stable than PS and RS with respect to a mechanical perturbation. This qualitative picture can be quantified by both experimentally measurable and computationally accessible thermochemical parameters. Namely: a) stability of PS relative to RS is given by the change in reaction free energy, ∆∆G b) stability of TS relative to RS is given by the barrier reduction, ∆∆g‡ Therefore, by comparing these two parameters we can determine the thermodynamic stability of critical points with respect to the applied stress and relative to each other. Our calculations show that the barrier reduction is always somewhat smaller than the reaction free energy change: |∆∆g‡| < |∆∆G|, however, the difference between the two is larger at small pressures and decreases as the pressure increases. The same line of arguments can be applied to every individual free energy component to further clarify the origin of mechanochemical sensitivity. The Gibbs free energy decomposition into internal energy, mechanical and entropic contributions demonstrates that the origin of the mechanical sensitivity of this reaction rate is almost entirely due to the mechanical term under small pressures (up to 2 GPa). Higher pressures also cause an increasingly uneven distribution of the internal energy. Thus, we can identify two regimes corresponding to low and high pressure and linear and non-linear free energy dependence on pressure. In the linear regime, electronic structures of molecules remain close to the unperturbed ones, while high pressures cause significant changes in the electronic structure of mechanophores. The stress-induced change in the reaction free energy, therefore, can be used as an upper limit for the barrier reduction, at least under small to intermediate pressures. The reaction free energy



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

change is about -17 kcal/mol at 7.3 GPa and is only about -3 kcal/mol at 0.9 GPa. Thus, we conclude that the estimate of 2 kcal/mol for the barrier reduction at 0.9 GPa made from one of the experiments is reasonable, while the estimate of 20.8 kcal/mol made from another experiment is an overestimate. Perhaps, the higher estimate could be reasonable for the maximum reaction sensitivity under high pressure. The electronic structure method which we used underestimates the activation energy by about 5-7 kcal/mol and overestimates the reaction free energy change by 1315 kcal/mol. However, pressure-induced differences in the activation energies are in excellent agreement with experimentally derived rate constants under small pressures due to error-cancellation. In the high pressure range, the reliability of the quantitative estimate decreases due to the method’s intrinsic error. However, quantitative reliability of these free energy estimates can be further perturbatively improved29 with high level theoretical methods (e.g. using the GPU-accelerated version of multi-reference49 approach).

energy of atoms, which is dissipated as heat. Thus, the mechanical sensitivity is dependent upon the difference in the rigidity (or the force compliance40) of the system at critical points along the reaction path. Here it should be noted that neither the critical points, nor the reaction path under stress are guaranteed to remain the same as the unperturbed counterparts. As long as Hooke’s law holds (which is typically true for relatively small deformations) the difference in the absorbed mechanical energy by a system in two different states will be linearly proportional to the applied force, F, with the slope ∆xRS-∆xPS (note that if kmin

Molecular Origin of Mechanical Sensitivity of the Reaction Rate in

Recommend Documents