The global kinetic–thermodynamic relationship derived from first principles

Abstract

What governs the relationship between the reaction rate and thermodynamic driving force? Despite decades of rate theory, no general physically grounded equation exists to relate rate and driving force across all regimes. Classical models, such as the Marcus equation and Leffler equations, either rely on under-realistic assumptions or only capture the local behaviour, failing outside narrow regimes. We derive a general, non-linear equation from microscopic reversibility, arriving at three physically meaningful parameters: a minimum preorganisational barrier (E_min), a reaction symmetry offset (E_eq), and a kinetic curvature factor (θ). This model captures global behaviour, recovers known limits, and explains why classical models like Leffler equation exhibit the rate–driving force responsiveness (the Brønsted slope) as they do, by revealing their physical origin, not by fitting them. The model enables a causal reinterpretation of experimentally observed curved rate–driving force plots, such as in hydrogen atom transfer to Fe(IV)=O. Importantly, this model does not require replacing existing models, it explains their physical foundations, enabling chemists to continue using them while understanding when and why they apply, and where they break down. Beyond case studies including hydride shifts, rearrangements, and cyclisations, the framework's strength lies in its deductive foundation enabling the physically grounded design of reactions with desired kinetics across diverse chemical systems. By revealing the global structure of the rate–driving force relationship, this framework enables chemists to recognise, predict, and design reactivity that would otherwise appear anomalous or inaccessible, better clarifying the unknowns. Examples include highly exergonic regimes near E_min, where further increases in exergonicity offer little rate improvement, and control shifts to structural factors; even when the rate–driving force plot appears linear, the model uncovers hidden curvature and deeper physical meaning.

---

1 Introduction

What determines how reaction rates respond to a thermodynamic driving force? Why do some reactions speed up dramatically with slight changes in driving force, while some others appear much less responsive? The interplay between kinetics and thermodynamic driving forces has been a cornerstone of physical organic chemistry. The most elegant aspect of Marcus theory lies in dissecting the reaction barrier into thermodynamic-independent and thermodynamic-dependent components.¹ This dissection immediately raises a critical question: what are the “thermodynamic-independent” parameters defining the inherent nature of a family of reactions?² This question forms the foundation of different attempts to describing the interplay, leading to different models. However, many previous models rely either on excessively strong assumptions,^1–3 or, in attempts to be more comprehensive, on extensive parametrisation of reaction-specific concepts,⁴ which greatly restricts the scope and the ability to derive meaningful physical insights from these frameworks due to their more statistical nature. These include models that incorporate bond breaking and bond formation to estimate the energy of the transition state, which often require treatments very specific to one reaction class.⁴ Existing models interpret these cases as breakdowns or edge effects, but perhaps they are evidence of a deeper, universal relationship. Addressing these limitations requires a model that captures the global relationship between activation and reaction energies across the entire thermodynamic range that relies on minimal assumptions. Existing models are often thermodynamically narrow, too empirical, or too dependent on system-specific parametrisation. What has been missing is a unified, physically interpretable framework that recovers known behaviour as limits, and explains why existing models succeed or fail.²

For example, the most experimentally relevant is the Leffler equation (Fig. 1a), which expresses the energy of activation ΔE^‡ by a combination of the reaction energy ΔE and the intrinsic barrier ΔE^‡₀, which corresponds to ΔE^‡ for ΔE = 0, in the linear approximation, ΔE^‡ = ΔE^‡₀ + αΔE.² The Brønsted slope α indicates the sensitivity of the reaction rates to changes in driving force. ΔE^‡₀ and α are intrinsic properties of each family of reactions, presumed to be independent of the thermodynamic driving force. For constant ΔE^‡₀ and α, the Leffler equation is reduced to the Bell–Evans–Polanyi (BEP) principle, describing the kinetic–thermodynamic relationship within the same family of reactions.⁵ The value of α is often used to probe how early or late in a reaction the transition state develops, with values lower than 0.5 corresponding to an earlier transition state. Although the linear kinetic–thermodynamic relationship is commonly observed experimentally, it is the approximation over a limited range of thermodynamic driving force. As the driving force becomes highly exergonic or endergonic, the linearity breaks down.⁶ In the extreme limits, the behaviour of α reveals distinct kinetic regimes. For highly exergonic reactions, α approaches 0; conversely, for very endergonic reactions, α approaches 1.⁷ BEP plots often show deviations from linearity.⁸ These outliers prompt a critical question: do they indicate curvature due to driving forces exceeding the linear range, or do they reflect reactions from distinct families with different intrinsic properties? Analysing energy relationships in their exact non-linear form is essential to addressing this, as true chemical information can only be found when moving beyond parameters for a local interpolation. A universal, deductive relationship between rate and reaction energy remains conspicuously absent. This absence has not gone unnoticed, but it is often obscured by models that appear to work locally while failing to generalise, with all current models either empirical, local, or reliant on idealised reaction assumptions. Marcus derived a curved equation from the analysis of crossing the reactant state parabola and the product state parabola (Fig. 1b), and a follow-up work offers a thermodynamically rigorous description.⁹ This model predicts a Brønsted slope α near 0.5 for reactions close to zero driving force. A serious assumption in the Marcus formation of the curve is that the parabolas of the reactant and the product states have identical shape. These shapes are influenced by factors such as bonding, conjugation, steric properties, solvation, and non-covalent interactions.¹⁰ This oversimplification obscures the true origins of deviations from linear BEP relationships. Koeppl and Kresge demonstrated that allowing bond strengths and reaction distances to vary realistically produces a sigmoid Brønsted slope, better reflecting real systems,¹¹ which would deviate from Marcus predictions.¹² Similarly, Richard and Jencks found significant deviations in Brønsted slopes for general base catalysis, attributing them to differential charge development in the transition state and reactant-product asymmetry.¹³ These studies introduce reaction-specific corrections. Each of these formulations and models rely on distinct assumptions and, while having advanced our understanding of kinetic–thermodynamic relationships, still often address aspects in isolation.

Fig. 1 Models in describing the interplay between thermodynamic-independent and thermodynamic-dependent components of chemical reactivity: (a) Leffler equation as a local linear fitting with an empirically derived slope α treated as a constant, (b) Marcus model predicting a curved relationship, but with inaccurate curvature and limiting behaviour, (c) our derived exact model yielding three key parameters with physical meaning: a minimum preorganisational barrier (E_min), a reaction symmetry offset (E_eq), and an energetic coupling factor (θ). Our causal model explains why the Leffler equation yields the Brønsted slope value as it does and clarifies the physical basis of the observed curvature.

Herein, we present a model that originates from a distinct foundation and goes significantly beyond the scope of existing frameworks. In this study (Fig. 1c), we derive a non-linear equation grounded in the principle of microscopic reversibility, with minimal assumptions, which captures (or “arrives at”, given its deductive nature) all these essential physical features: reaction asymmetry, reactant preorganisation, product reorganisation, and the sensitivity of kinetic–thermodynamic response at each thermodynamic driving force. Together, these portray the full reactivity landscape. In comparison to G^‡₀ and α that represent the local thermodynamic-independent properties of a family of reactions, we introduce a new set of thermodynamic-independent parameters overarching the global thermodynamic range. Each parameter has a distinct physical interpretation, grounded in common chemical meaning (in contrast to statistical in nature). We demonstrate the applicability of this equation through diverse chemical reactions, illustrating its capacity to describe non-linearity, identify valid linear approximations, and elucidate the role of each thermodynamic-independent parameter. Experimental observations show that some families of reactions exhibit high sensitivity to rate–driving force responsiveness, while others display relatively insensitive correlations.¹⁴ Some sets of reactions exhibit crossing behaviour: one is more reactive in the endergonic region, while the other is more favoured in the exergonic region of thermodynamic driving forces. These distinctions raise practical questions: What should we do to accelerate a given reaction? To what extent will making a reaction more thermodynamically downhill improve its rate? Our derived equation also provides a lens through which to examine these questions. This framework reveals the physical origins of observed Brønsted slopes, showing why they have the values they do across different regimes, and recovers them as limiting behaviour under a broader kinetic–thermodynamic relationship.

2 Results and discussion

Finding a global model for the kinetic–thermodynamic relationship generalisable for any chemical reaction must arise from the most fundamental characteristics. Any approach from linearity or parametrisations specific to one reaction class, even when not ad hoc, hamper both the physicality and generality of the resulting equations. For these reasons, we investigated what mathematical features should be present in the global energy responses, followed by a deductive approach in agreement with the former and stemming from an interpretation of microscopic reversibility.

2.1. Derivation of the constraints and properties of the non-linear relationship

In order to reach a model dependent on physically meaningful parameters, the question must be approached in the most general way possible by considering the properties common to any kind of chemical reaction. Through this, we should deductively arrive at a set of constraints from characteristics general to all energy responses, rather than assuming that a particular behaviour should apply to the non-linear relationship. For a given family of reactions, it would be possible to determine ΔE^‡ as a function of ΔE_r, where ΔE_r is the reaction energy and ΔE^‡ is the activation energy. In any chemical reaction, the reaction energy can be expressed as the difference between the forward and backward energy barriers (Fig. 2).

Fig. 2 The activation energies and reaction energy of a reverse reaction are fundamentally connected to the forward reaction.

This non-linear function would have to satisfy a number of physically meaningful requirements that define its properties, particularly those concerning the symmetry of the function and its derivatives. As there is a fundamental relationship between the energies of a given reaction and its reverse, the function linking the reaction energy to the forward energy barrier, ΔE^‡(ΔE_r) and the one linking the (negative) reaction energy to the reverse energy barrier, ΔE^‡₋₁(−ΔE_r), have to be connected somehow. This feature of all reactions constrains the nature of the function. By considering how the energy barriers and reaction energies invert between a general forward and reverse reaction, the following constraints on the nature of the non-linear relationship are obtained:

ΔE^‡₋₁(−ΔE_r) = ΔE^‡(ΔE_r) − ΔE_r (1)

−[ΔE^‡₋₁]′(−ΔE_r) = −1 + [ΔE^‡]′(ΔE_r) (2)

[ΔE^‡₋₁]′′(−ΔE_r) = [ΔE^‡]′′(ΔE_r) (3)

Eqn (1) relates the activation energy of a reverse process to the activation energy of the forward process, which are always different by the reaction energy, ΔE_r. Therefore, given a constant horizontal displacement in ΔE^‡₋₁, a vertical displacement of the same magnitude would appear in ΔE^‡. By differentiating eqn (1) once and twice, we obtain eqn (2) and (3), respectively. The first derivative shows that the gradients of the forward and reverse reactions are complementary, adding up to one. The second derivative, the curvature, of the backward reaction is the mirror image about the y axis of that of the forward reaction (as −ΔE_r,−1 = ΔE_r) (Fig. 3).

Fig. 3 Graphical representation of (a) the non-linear kinetic–thermodynamic response, eqn (1); (b) its gradient, eqn (2); and (c) its curvature, eqn (3); for a set of forward reactions S and the set of reverse reactions S₋₁.

The additional conceptually derivable constraints on a function that models the kinetic–thermodynamic relationship involve the limits at extreme driving forces. Due to the TS interpolating between the minima, the regimes at extreme driving forces must be asymptotic to the energies of the minimum closest in energy. As the reaction energy approaches negative infinity, the factors stabilising the product cannot be present in the TS and so the energy of the TS must only depend on that of the reactant. In this limit, the gradient must approach zero, with the barrier approaching , where E_preorg represents the ΔE_r-independent minimum barrier to be overcome regardless of the driving force. The reverse is also true for extremely endergonic reactions. In the limit, the barrier must be equal to the reaction energy plus a reactant-independent offset and, hence, independent of ΔE_r. This is represented as ΔE_r + E_reorg, where E_reorg is the ΔE_r-independent reorganisation energy of the product. As these limits are linear, the gradient of the non-linear relationship must transition between its asymptotic limits at 0 and 1. For an idealised reaction following Hammond's postulate, and so in which the geometry of the TS is fully and exclusively defined by its connected minima, only the factors influencing the TS are present to various degrees in the connected minima, and so the response of ΔE^‡ must vary with the energetic proximity to the minima. As the position of the energy barrier stems from an interpolation of stabilising factors that change monotonically in elementary processes, the response of ΔE^‡ thereto must be of the same sign as the response of ΔE_r. In cases where this does not hold, the TS cannot be seen as an interpolation of factors. This situation only occurs if the potential energy surface shows extensive asynchronicity, caldera-like regions or bifurcation-like behaviour.¹⁵ These effects of asynchronicity can lead to incorporation of extraneous stabilisation factors, as one of the minima would be less geometrically related to the TS, and presumably significant changes in reaction energy unrelated to the factors affecting the barrier. Therefore, for any single synchronous elementary process, must always be positive and monotonically increasing, with the second derivative positive for all values of ΔE_r.

We now consider an identity reaction rⁱ, with its reactants and products being identical. There are many examples of such identity reactions, including degenerate Cope rearrangements, degenerate deprotonations with a conjugate base, and type I and type II dyotropic reactions. This reaction belongs to a family of analogous reactions, A, such that the corresponding family of backward reactions A₋₁ is equivalent. Such reactions may be labelled isodromic, as the reverse reactions show the same features as the forward reactions, stemming from how the reactants form products of similar nature.¹⁶ The energy response of these systems perfectly aligns with the notion of early or late TS, with exergonic reactions showing earlier transition states and thermoneutral reactions having the TS at the middle of the process. Being so closely related to the identity reaction, these isodromic reactions would energetically behave like the identity reaction upon the appearance of a driving force due to changes in reaction energy, the only difference between forward and backward reactions being the reaction energy. As these are indistinguishable, the kinetic–thermodynamic relationship of the forward and backward isodromic reactions, where an identity reaction belongs by definition to both sets, must be modelled by the same function as by definition they show the same energy-related features. The identity reaction is, by definition, thermoneutral, and its reverse is itself, and so ΔE^‡(ΔE_r) = ΔE^‡₋₁(ΔE_r) for any idealised isodromic family containing rⁱ. Following eqn (2), the gradient [ΔE^‡]′(0) = 0.5 and so the function is symmetric about (0, 0.5).¹⁷

Conversely, for a family of anisodromic reactions, B, such that no reaction of B₋₁ belongs to B, the curves ΔE^‡ and ΔE^‡₋₁ would intersect at ΔE_r = 0 by definition. However, their gradients at ΔE_r = 0 could be different from 0.5. This is clearly identifiable as thermoneutral reactions would still be found to be early or late in spite of having no apparent thermodynamic bias, resulting in a mismatch with any symmetric model.¹⁸ The TS of thermoneutral anisodromic reactions could thus in principle be displaced from the centre of the reaction coordinate and so early or late TSs switch at a reaction energy different from 0 (Fig. 4). Applying eqn (2), which provides the symmetry constraints of the gradient, shows that there is an energy E_eq for which [ΔE^‡]′(ΔE_r + E_eq) = [ΔE^‡₋₁]′(ΔE_r − E_eq) = 0.5 (Fig. 5). Hence, [ΔE^‡]′ and [ΔE^‡₋₁]′ are symmetric about and (−E_eq, 0.5), respectively. E_eq would represent the point at which the influence of reactants and products on the energy barrier is equal. This is already unaccounted for in the Marcus quadratic equation, for which the gradient at 0 is invariant since, by construction, the Marcus equation assumes perfectly isodromic behaviour. This can cause large deviations for anisodromic reactions.¹⁷

Fig. 4 Reaction coordinate for thermoneutral (plain) and non-thermoneutral (dashed) reactions of differing isodromicity. A vertical line shows the perfect symmetry for the thermoneutral isodromic reaction, as it would be an identity reaction.

Fig. 5 Relationship between the gradient of the kinetic–thermodynamic response of a set of anisodromic reactions B and the set of their reverse counterparts B₋₁. The points where their gradient is 0.5 are separated from the y axis by E_eq.

Consequently, as both [ΔE^‡]′(ΔE_r + E_eq) and [ΔE^‡₋₁]′(ΔE_r − E_eq) are symmetric about (0, 0.5), and their curvature is identical (eqn (3)), they must be equal. We can of course integrate these with respect to the reaction energy to recover the kinetic–thermodynamic relationship. This means that the kinetic–thermodynamic relationship functions themselves, ΔE^‡(ΔE_r + E_eq) and ΔE^‡₋₁(ΔE_r − E_eq), can only differ by a single constant (the integration constant) for any ΔE_r. The difference between their limits at negative infinity, E_preorg − E_preorg,−1, and at positive infinity, (E_reorg + E_eq) = (E_reorg,−1 − E_eq), must be equal. As the reorganisation of the product from the TS is, by definition, equivalent to the preorganisation of the reverse reaction, E_reorg = E_preorg,−1, E_eq is obtained from the limits of the function:

E_eq = E_preorg − E_reorg (4)

The value of the barrier at E_eq can only depend on the curvature, [ΔE^‡]′, which itself controls the transition between the limiting linear behaviour of the energy barrier at extreme driving forces. Owing to the symmetry constraints (eqn (3)), the curvature reaches a maximum at E_eq, while the gradient of the function was shown to be limited to the range between 0 and 1. For a self-similar curvature profile, invariant under rescaling as would be expected for a single-mechanism interpolation, the value of ΔE^‡(E_eq) is therefore inversely proportional to the curvature at that point.¹⁹ The curvature is responsible for describing the coupling of the two minima to the TS and is the defining feature of the non-linear equation, with the barrier at E_eq being a direct consequence thereof.

2.2. Derivation of equation

We return to ideal isodromic processes, where the stabilising factors affecting the reactants and products would be expected to impact the energy of the transition state proportionately to their influence on the geometry of the transition state. From the previous explored constraints, this influence must be dependent on the energy of the minima, consequently coupling the energy of the TS to those of the intermediates. Ultimately, an exponential dependence on energy can be expected, as it is found in analogous questions such as transition probabilities, rate constants and other aspects in non-equilibrium transition-state theory.²⁰ This would be derived from the energy difference, not unlike a Boltzmann distribution, albeit with a scaling parameter different to inverse temperature. For an isodromic reaction set, in which the reactants and products are of the same nature, any changes in energy on either of the minima would result in an identical response on the transition state (TS). Drawing an analogy to microscopic reversibility, as the forward and backward factor transmission influences constitute the total origin of all factors acting on the TS from either minimum, they must add to 1. This can be expressed as a normalised sum of exponentials, where θ is the reaction-specific scaling factor related to the responsiveness of the barrier (eqn (5)):

The value of θ captures the partition of influences and, hence, must be a positive real number. Larger values would correspond to a stronger coupling, where the TS is highly influenced by the minima, whereas values closer to zero would mean that the transition state responds less rapidly to thermodynamic changes. This can be rationalised as the stabilisation of the TS arising from the competing factors modifying the reactants and the products. The ease of transmission of stabilising factors to the TS would thus be modulated by θ: higher values would lead to a TS more responsive to the factors that stabilise the minima. Conceptually, this would represent the information about the reaction hitherto usually ascribed to the Brønsted slope. Re-expressing ΔE^‡₋₁ in eqn (5) as (ΔE^‡ − ΔE_r) and rearranging²⁰ shows the energetic dependence of the barrier with the minima it connects as an explicit function, for the idealised set of reactions:

From eqn (5), for the thermoneutral isodromic identity reaction (with E_eq = 0), the forward and backward barriers are the same, . This inverse proportionality to θ matches the constraint we found regarding the connection between curvature and the barrier at E_eq, with the curvature at that point being . As explored earlier, the asymptotic minimum energy of preorganisation could be non-zero. In addition, for a set of anisodromic reactions B, the expression would be displaced by E_eq, and so by the point with a gradient of 0.5 for reactions with unequal responses, giving rise to the final expression for the function:

The limits of this function at negative and positive infinity are E_min and ΔE_r + (E_min − E_eq) respectively, which match the derived constraints for the asymptotic preorganisation and reorganisation energies, as well as the definition of E_eq. Note that the constraints were derived from the general behaviour of forward and backward reactions (vide supra), whereas the equation was derived independently from the balance of influences in eqn (7). Hence, both the correct asymptotic behaviour and the symmetry of the gradient expected from the constraints are reached separately and serve as additional validation of the energy relationship. The effects of each parameter energy response are shown in Fig. 6.

Fig. 6 Graphical representation of the effects of each parameter on the behaviour of the equation. Note that the gradient at E_eq is always equal to 0.5.

While the kinetic–thermodynamic plots of some reactions exhibit observed slopes outside of the range 0–1 (such as the nitroalkane anomaly or the Marcus inverted region),²¹ these observed local gradients do not contradict our model but instead arise from parallel system-dependent effects. In most known examples, the disruption of the monotonicity of the kinetic–thermodynamic response occurs because of the weakening of Hammond's postulate, with the outcome of such a treatment akin to a Hammett-type correlation instead of a causal relationship, which sometimes results in the inverted region. Experimental evidence supports this, as in most observed inverted regions, the inverted Brønsted slopes vary significantly in magnitude,²² deviating from the expected same-magnitude-inverse-sign predicted by the quadratic equation in Fig. 1b. In fact, and in any case as explored in Section I, a non-monotonic energy response cannot be compatible with Hammond's postulate. In ET and many other reactions, asynchronicity is common,²³ which limits the validity of Hammond's postulate for the particular reaction, while in intersystem crossings there is no required smooth interpolation of the potential energy surface. With asynchronicity, as more than one process occurs, such trends are context-dependent and have also been related to the nitroalkane anomaly,²⁴ reinforcing the Hammett-like behaviour observed in these cases. Only in this scenario, perhaps closer in nature to a collection of independent steps than to a concerted process, would the TS correlate weakly enough for monotonicity not to be preserved: this can lead to cases where local observed Brønsted slopes were negative or higher than 1, a few examples of which have been presented as outliers in the literature.²⁵ In this light, our equation must always hold unless the reaction is highly asynchronous or the reaction class is ill-defined due to external factors. It offers a more fundamental framework for distinguishing when additional factors shape reaction behaviour and prompts the need to capture the evolving phases of a reaction, which in turn also enables conceptually more accurate and predictive kinetic analysis. Anomalous kinetic behaviour can only then be systematically interpreted, instead of heuristically fitted.

2.3. Demonstrations on chemical reactions

A particularly clear example of experimentally found non-linearity in kinetic–thermodynamic relationships is that by Groves and coworkers involving the hydrogen atom transfer from a number of organic molecules to a key iron(IV) oxo complex in a peroxygenase (Fig. 7).²⁶ While treating all C–H bonds amenable to activation in each substrate as equivalent and observing the average effect, as well as including substrates of a very divergent nature, with the consequent spread in the values, the results displayed visible non-linearity. The pseudo-barriers derived from the rate constants shifted from an essentially asymptotic unresponsive region to a more responsive section. Our derived equation provides a key causal link that accounts for these prior experimental findings. Practically, this implies that further increasing the thermodynamic driving force would no longer accelerate the reaction in this unresponsive region, reflecting an insight that the Leffler's equation cannot tell.

Fig. 7 Experimental observation of a non-linear kinetic–thermodynamic response in hydrogen atom transfer to an iron(IV) oxo complex.²³ The horizontal axis is the C–H bond dissociation energy (BDE) of substrates, the difference of which gives the variation in reaction enthalpy. The vertical axis is the pseudo-barrier obtained by the rate constant, where C = RT ln(k_BT/h), with R as the gas constant, k_B as the Boltzmann constant, h as Planck's constant, and T = 277 K in this example.

In order to better explore the non-linear responses to demonstrate the emergence of the derived behaviour, we studied a series of model systems with tunable reaction energies while retaining a closer structural similarity among substrates. In this way, the kinetic–thermodynamic relationship of the family of reactions would clearly manifest the three parameters we defined. To see non-linearity clearly, the reactions should display easily modulated driving forces over a wide range of energies while retaining the same orbital makeup. For these reasons, we first explored a 1,2-hydride shift on a vinyl cation, a system adapted from a reported reaction.²⁷ In order to determine the energies associated with this process, we performed density functional theory (DFT) calculations on the system. The reported energies were calculated, in the gas phase, in ORCA 6.0 (ref. 28) using the dispersion-corrected²⁹ ωB97X-D3 functional³⁰ and the def2-TZVPPD basis set³¹ for all atoms, with the RIJCOSX³² approximation.³³ A family of 1,2-diarylethenyl cations was calculated with substitution on one of the aryl rings, together with the corresponding 1,2-hydride shift transition states and their products (Scheme 1).

Scheme 1 Substrates studied in the 1,2-hydride shift of 1,2-diarylethenyl cations. Non-linear relationship between the activation and reaction potential energies (circles) and Gibbs energies (triangles) of 1,2-hydride shifts on diarylethenyl cations listed in Scheme 1. R² = 0.9625 (ΔE), 0.945 (ΔG).

The geometry of the transition state showed an essentially linear geometry of both central carbon atoms, best represented as formal η² coordination of a diarylacetylene to the migrating hydrogen atom. In accordance with Hammond's postulate, the more favourable reactions (substrates A17 and A18) showed an earlier transition state,³⁰ with the C–H distances reflecting the position of the transition state across all substrates. The corresponding plot showed a well-defined curve in a driving force range of around 25 kcal mol⁻¹ (Scheme 1b, circles).³⁴

The relationship evidences a large, thermodynamically independent contribution to the energy barrier of 12.5 kcal mol⁻¹. This is presumably linked to the change in hybridisation and consequent weakening of the C–H bond. In addition, the curvature of the plot, defined by θ, is markedly visible. This would stem from the geometric proximity between the TS and the minima and an ensuing stronger impact of the same factors on the three stationary points. Such chemical proximity is captured by the coupling of the energy of the transition state to that of the connecting minima (θ = 0.0945 mol kcal⁻¹), which is stronger than in the other studied reactions in this work. Finally, the E_eq of 7.7 kcal mol⁻¹ hints at a degree of asymmetry in influence, showing how a hydrogen migration from a variable aryl-substituted sp² carbon to a phenyl-substituted one responds differently because of the modified fragment. This results in the reorganisation, or the E_min of the reverse reaction (E_min − E_eq), being only 4.8 kcal mol⁻¹, showing that the thermodynamic dependence of 1,2-hydride shifts depends less on the substitution of the product carbocation and more on that of the reactant carbocation. Additional insights can be gained: if the E_min is lower when the migration occurs from the variable aryl, this suggests that the distortion will vary with the substituents more strongly (aligning better with thermodynamics) than if the substituents are located on a different ring, matching general chemical intuition. Through the thermodynamic-independent minimum preorganisation and reorganisation energies, E_min and E_eq can be instrumental in explaining whether the initial distortion is eased by the chemical modifications.

Plotting the same results with Gibbs energy also lead to an analysable plot (Fig. 4, triangles), although the dispersion of the data points was slightly higher, presumably owing to the additional non-trivial effects of the entropy and zero-point energy. The change in the parameters shows G_min (the Gibbs energy analogue of E_min) to be 1.9 kcal mol⁻¹ lower than E_min, G_eq (the Gibbs energy analogue of E_eq) 2.7 kcal mol⁻¹ lower than E_eq, and the resulting reorganisation energy (G_min − G_eq) 0.8 kcal mol⁻¹ higher than its potential energy analogue. These changes suggest that the inclusion of thermal terms, likely with the changes related to zero-point energy as the greatest contribution, significantly reduces the barriers at higher driving forces, and with the shift of G_eq – towards a more equal response – also conditions how the ΔG^‡(ΔG_r) function grows. However, the change in θ makes the activation energy much more responsive to changes in the thermodynamic driving force. On its own, θ_G accounts for a reduction of up to 2.5 kcal mol⁻¹ at G_eq with respect to θ_E, and a reduction of at least 1 kcal mol⁻¹ between around ΔG_r = −15 and 25 kcal mol⁻¹. Had θ remained the same, the stabilisation would have been lower than 1 kcal for ΔG_r ≥ 0 kcal mol⁻¹, and actually destabilising for barriers above ΔG_r = 15.7 kcal mol⁻¹. This exemplifies how the observed energy barriers arise from a complex interplay of the parameters, even when studying the effect of thermal corrections. The influence of entropy is known to explain patterns in selectivity,³⁵ reaction kinetics,³⁶ and reaction outcomes,³⁷ with consequences on the free energy surface including the appearance of entropic intermediates.³⁸ Here, entropy can be seen to change not just the absolute value of the activation energy itself, but also its responsiveness to reaction energy. The changes in the non-linear parameters can enable the analysis of this effect.

We then investigated the kinetic–thermodynamic response in the gold(I)-catalysed 5-exo-dig cyclisation of a set of 2-aryl-1,5-enynes C, affording vinylgold(I) species D (Scheme 2). Using a trimethylphosphine ligand, as a reasonable model for most experimentally used phosphines, the aim of showing how the global kinetic–thermodynamic relationship can also be realised in synthetically relevant metal-catalysed transformations. Varying the substitution on the phenyl ring with a total of 25 different examples allowed modulating the reactivity of the alkene, as well as changing the stability of the product. This resulted in a reaction energy range of about 30 kcal mol⁻¹, allowing the observation of the curvature in spite of the inherent complexity of the reaction and the inclusion of implicit solvation.

Scheme 2 Substrates studied in the gold(I)-mediated 5-exo-dig cyclisation of 2-aryl-1,5-enynes. Non-linear relationship between the activation and reaction potential energies (circles) and Gibbs energies (triangles) of 5-exo-dig cyclisations of 1,5-enynes listed in Scheme 2. R² = 0.9216 (ΔE), 0.8391 (ΔG).

These reactions showed a shallower curvature, as evidenced by its lower value of θ, at 0.0655 mol kcal⁻¹. In addition, the minimum preorganisation energy E_min was also particularly high, at 21.3 kcal mol⁻¹. This was expected for this reaction, as even with highly stabilising factors in the product, the TS must be reached by significant distortion of the alkyne and the alkene. This distortion is independent of the stabilisation of the final vinylgold(I) species and is not as strongly affected by the electronics as the reaction energy itself. The parent enyne is always intrinsically robust, as electronic factors do not radically alter the stability of the styrene fragment. However, due to the benzylic cation, the stability of the product is very strongly affected by the electronics. For that reason, the reverse barrier should be generally much more dependent on thermodynamics than the forward reaction. Thus, in very disfavoured reactions, the energetic cost to reach the TS is much lower. That intuitive behaviour is observed in the notably lower reorganisation energy for the process (E_min − E_eq) of −2.0 kcal mol⁻¹. This would mean that the barrier would be zero for any reaction with a reaction energy higher than +53.2 kcal mol⁻¹, were it to remain well-behaved and, in those, the product would return barrierlessly to the reactant configuration as the barrier would cease to exist. We can see that the parent enyne is intrinsically robust (it remains a stationary point at any driving force), at least with respect to modification on the ring, whereas the product carbocation is not intrinsically robust as it is not stable at all energy ranges (it would become a non-stationary geometry given a large enough destabilisation). Since the reactant and product of these reactions respond so differently to the same effect, consequently, the reaction is highly anisodromic, with an E_eq of 23.3 kcal mol⁻¹. The Gibbs energy analogue similarly led to a plot with a higher variability. However, there was a slight decrease in the minimum preorganisation, with G_min = 19.0 kcal mol⁻¹, accompanied by a slightly higher asymmetry in response (G_eq = 25.2 kcal mol⁻¹) and a marginally shallower curvature, θ_G = 0.0549 mol kcal⁻¹. Due to the greater scatter of the data and low local gradient, the significance of these small variations is harder to assess conclusively.

Wherever an apparent linear relationship is observed, leading to the use of the Leffler equation, the true physical origin of the observation can lead to conceptually very different interpretations of what happens in a given chemical system. Several of these examples are shown in Fig. 8, showing locally linear approximations may originate from different curves. Our model explains why rate–driving force relationships display the Brønsted slopes they do. The chemical insights gained from studying their origins are greater than from the simpler linear model because it enables the approach to rather complex questions. For instance, two non-crossing Leffler plots with different slopes are typically interpreted as representing distinct classes of reactions, each with its own set of thermodynamic-independent parameters. However, our model shows that such behaviour may also arise from a single reaction class sampled at different thermodynamic regions (Fig. 8a). Is the minimum preorganisation of the system always high, regardless of driving force (Fig. 8b)? Two apparently coincident linear plots may indeed belong to the same global non-linear energy relationship or happen to match an unrelated curve due to close asymptotic behaviour (Fig. 8c). These reinterpretations show that even when linear models like Leffler's appear valid, they reflect only a local view of a deeper structure. Our equation explains why they succeed where they do, providing a global context that transforms heuristics into physical understanding. Does the reaction ever become barrierless and perhaps only limited by the diffusion rate? Will extrapolating to more exergonic reactions allow a large improvement in selectivity between two competing barriers? As evidenced in the analysis of the reactions in this work, the magnitude of the thermodynamic-independent component in the preorganisation or reorganisation can be fully consistent with chemical intuition. Furthermore, in any case, the values of the parameters may explain patterns that are not immediately obvious, but are critical to explaining why an energy barrier has a particular energy. Understanding the origin of a given linear approximation can provide many more insights into a system of interest. Having examined systems with clearly visible curvature, we now turn to cases where the response appears deceptively linear and show how our model can still extract mechanistic insight.

Fig. 8 Conceptual map showing how classical rate–driving force relationships, particularly the Leffler equation, emerge as limiting cases of the global kinetic–thermodynamic model derived here. The figure illustrates how to causally interpret apparent linear relationships between the reaction rate and thermodynamic driving force within the broader, non-linear energy landscape governed by E_min, E_eq, and θ. (a) Two distinct linear plots may belong to the same curve, and so share all three parameters, or to two different ones. (b) Two linear plots that intersect may curve at similar rates, with only a minor difference in preorganisation, or respond very differently over a greater range of reaction energy. (c) Two apparently coincident linear plots may indeed belong to the same global non-linear energy relationship or happen to match an unrelated curve due to close asymptotic behaviour.

We considered the Beckmann rearrangement,³⁹ a classical reaction which would support a similar 1,2-migration to the hydride shift albeit concerted with the concomitant cleavage of a leaving group. In order to find several related curves, we calculated six different sets of reactions, each of them with a different migrating group. Within each set, the migrating group was kept constant while varying the adjacent substituent of the imine (Scheme 3). As the reaction involves the initial formation of a nitrilium species (F), these were used as the thermodynamic reference for the products.

Scheme 3 Substrates studied in the initial step of the Beckmann rearrangement, with a pair of numbers defining the substituents. Plots of six different reaction sets, each with a different migrating R group (labelled), with linear fitting shown. Each point is a reaction with a perturbation on the non-migrating R′ group. The leaving group is modelled as a water dimer. Observed Brønsted slopes and intercepts (intrinsic barriers) from a linear fitting.

For each migrating group, an essentially linear plot was obtained, suggesting that the value of θ was significantly lower in all cases (Scheme 3). However, each of the observed barriers at ΔE_r = 0 and Brønsted slopes were significantly different, the latter varying from 0.48 to 0.74. While the curvature is not directly observable, the general trend that higher intercepts display steeper gradients would be a result of the similarities of the reaction sets. Methoxymethyl substitution (E6-0 through E6-10) is the only example where several points (corresponding to E6-6, E6-9 and E6-10) appear significantly separated from the rest of the points, suggesting either the influence of non-perfect synchronisation in the process or the existence of markedly different stabilisation factors in some substrates.⁴⁰ As all linear plots constitute the local tangents to the corresponding curves at the average of the reaction energies of the data points, they cannot be directly compared without accounting for reaction energy. However, if one parameter is fixed, the other two become solvable with a linear plot. Studying what values of E_min, E_eq and θ could match the observed linear fits by keeping one constant allowed determination of the trifluoroethyl migrating group as an outlier, displaying a higher E_min than any other reaction class. Similarly, the only main difference with the methoxymethyl migrating group, and responsible for its much lower activation energies across the board, seems to arise from a lower reorganisation energy, and so its more isodromic energy response. This is precisely where the value of a physically derived model emerges: even where behaviour appears linear, the framework reveals the hidden structure beneath. This analysis shows how misleading apparent linearity can be, and how only a global model reveals what parameters actually control reactivity. In practical terms, this enables synthetic chemists to tune barriers through design principles derived from physical causality, not regression.

To the question of how one may increase or lower the rate and selectivity of a given reaction, the answer involves the interplay of the exergonicity of the reaction with respect to G_eq, the curvature of the non-linear ΔG^‡(ΔG_r) function and the proximity to its asymptotic behaviour. This contributes to understanding the patterns in general chemical intuition, where different reactions respond differently to modifications, while providing a measure of which aspects are unnecessary to modify. The non-linear equation, even when only implicitly observed due to local linearity, provides insights into how reactivity can be modulated. In a highly exergonic regime, while a simpler linear plot might suggest modifications could be futile, our equation would suggest that an increase in θ could be very effective, attempting to make the barrier more sensitive to the stabilising factors through changes such as steric environment or strain. This breaks away from the reaction energy-based approach with a linear model (and mostly intractable outliers), instead providing alternatives for synthetic chemists to improve even already very exergonic reactions and rationally lower their barriers. This approach can be of utility regardless of the exergonicity of the reaction. Its application in new methodology will be ultimately determined by the understanding of how different factors affect the parameters that guide the non-linear kinetic–thermodynamic relationship. The answers to these questions should lie closer to the chemical properties, as the meaning of θ, E_min and E_eq is inherently fundamentally distinct from heuristic interpolation. Full accommodation of asynchronous processes is ongoing in our laboratory and would further broaden the scope of applicability.

3 Conclusions

We have derived a universal non-linear relationship between the reaction energies and activation energies of a reaction set that not only captures rate–driving force relationships across thermodynamic regimes but also explains why classical models like Bell–Evans–Polanyi relationships and Marcus equation succeed or fail in specific limits. The derivation is grounded in a formulation stemming from the principle of microscopic reversibility, with an exponential dependence on the energies, and otherwise minimal assumptions about the reactions themselves. Consequently, the observed Brønsted slopes and intercepts of the linear approximations of a system constitute local descriptors, dependent on the underlying parameters and the studied thermodynamic range. The full non-linear understanding allows global comparisons between very different ranges of reaction energies and reaction families. We expect that factors that affect the sensitivity of energy barriers to the reaction energies to correlate strongly with the parameters (θ, E_min and E_eq) governing the non-linear relationship, and therefore only indirectly with the local slopes and intercepts. This deductively reached global non-linear model allows direct mechanistic comparisons between different reaction classes across wide ranges of reaction energies. In systems with visible curvature, such as 1,2-hydride shifts and gold(I)-catalysed cyclisations, the model enables the extraction of the three parameters and how energy barriers depend on their combined effects. Even where linear fitting may appear adequate, this is deceptively simple and the non-linear model clarifies how differences in parameters, not just thermodynamic range, govern barrier variation. This reframing offers a more unified and practical understanding of reactivity, particularly for identifying outlier behaviour and guiding new reaction discovery, by moving from local interpolation to derived reactivity parameters. While the model can reinterpret past observations, such as the non-linear rate–driving force behaviour in iron(IV) oxo-mediated reactions, its greater value lies in guiding new mechanistic design and uncovering reactivity principles that remain hidden under traditional frameworks. It remains open for future work to discern what chemical properties impact the values of θ, E_min and E_eq, together with how their rational modulation can lead to its adoption in organic and organometallic chemistry. Through the exploration of the kinetic–thermodynamic relationship with the new overarching model, we acquire a better perspective on the importance of driving force in reaction kinetics. This work invites to rethink conventional strategies and intuition on the effects of different types of perturbation on chemical reactions: what appears locally adequate may conceal global inconsistency, and this model equips chemists to see beyond the curve, to recognise structure where none was thought to exist.

Author contributions

E.G.-P. derived the equation and led the study. G.Q. supervised the project. Both authors jointly wrote the manuscript.

Note after first publication

This article replaces the version published on 3rd September 2025, where eqn (5) was appearing incorrectly, this has now been resolved.

Conflicts of interest

There are no conflicts to declare.

Data availability

All computed data are available within the paper and its SI. See DOI: https://doi.org/10.1039/d5sc04829j.

A collection of the computational data underlying this work has been made available in the ioChem-BD repository at https://doi.org/10.19061/iochem-bd-6-443.⁴¹

Acknowledgements

We thank the MPI für Kohlenforschung for funding. We gratefully acknowledge the storage resources at BSC provided to ioChem-BD. Open access was funded by Max Planck Society.

Notes and references

1. (a) R. A. Marcus, J. Chem. Phys., 1956, 24, 979–989; (b) R. A. Marcus, J. Chem. Phys., 1956, 24, 966–978; (c) R. A. Marcus and N. Sutin, Biochim. Biophys. Acta, 1985, 811, 265–322.
2. N. Agmon, J. Chem. Soc., Faraday Trans. 2, 1978, 74, 388–404.
3. (a) J. E. Leffler, Science, 1953, 117, 340–341; (b) G. S. Hammond, J. Am. Chem. Soc., 1955, 77, 334–338.
4. (a) J.-M. Savéant, Acc. Chem. Res., 1993, 26, 455–461; (b) D. Bím, M. Maldonado-Domínguez, L. Rulíšek and M. Srnec, Proc. Natl. Acad. Sci. U.S.A., 2018, 115, E10287–E10294.
5. (a) R. P. Bell, Proc. R. Soc. A, 1936, 154, 414–429; (b) H. Eyring and M. Polanyi, Trans. Faraday Soc., 1936, 32, 1333–1360.
6. (a) J. N. Brønsted and K. Pedersen, Z. Phys. Chem., 1924, 108, 185–235; (b) D. Usharani, W. Lai, C. Li, H. Chen, D. Danovich and S. Shaik, Chem. Soc. Rev., 2014, 43, 4968–4988.
7. M. Eigen, Angew. Chem., Int. Ed., 1964, 3, 1–72.
8. (a) W. D. Morris and J. M. Mayer, J. Am. Chem. Soc., 2017, 139, 10312–10319; (b) X. Chen and J. I. Brauman, J. Am. Chem. Soc., 2008, 130, 15038–15046; (c) J. J. A. Garwood, A. D. Chem and D. Nagib, J. Am. Chem. Soc., 2024, 146, 28034–28059; (d) Y. Sun, F. Liu, J. N. Sanders and K. N. Houk, J. Org. Chem., 2023, 88, 12668–12676; (e) T. V. Chciuk, W. R. Anderson Jr. and R. A. Flowers II, J. Am. Chem. Soc., 2018, 140, 15342–15352.
9. R. A. Marcus, J. Phys. Chem., 1968, 72, 891–899.
10. (a) S. Fletcher, J. Solid State Electrochem., 2010, 14, 705–739; (b) G. W. Koeppl and A. J. Kresge, J. Chem. Soc., Chem. Commun., 1973, 371–373; (c) E. Laborda, M. C. Henstridge and R. G. Compton, J. Electroanal. Chem., 2012, 667, 48–53; (d) E. Laborda, M. C. Henstridge and R. G. Compton, J. Electroanal. Chem., 2012, 681, 96–102.
11. G. W. Koeppl and A. J. Kresge, J. Chem. Soc., Chem. Commun., 1973, 32, 371–372.
12. J. P. Richard, G. Williams, A. C. O'Donoghu and T. L Amyes, J. Am. Chem. Soc., 2002, 124, 2957–2968.
13. (a) J. P. Richard, J. Org. Chem., 1994, 59, 25–29; (b) J. P. Richard, T. L. Amyes and M. M. Toteva, Acc. Chem. Res., 2001, 34, 981–988; (c) J. P. Richard and W. P. Jencks, J. Am. Chem. Soc., 1984, 106, 1396–1401.
14. (a) T. F. Markle, J. W. Darcy and J. M. Mayer, Sci. Adv., 2018, 4, eaat5776; (b) G. Qiu and R. R. Knowles, J. Am. Chem. Soc., 2019, 141, 2721–2730; (c) G. Qiu and R. R. Knowles, J. Am. Chem. Soc., 2019, 141, 16574–16578; (d) G. Qiu and P. R. Schreiner, ACS Cent. Sci., 2023, 9, 2129–2137.
15. (a) W. P. Jencks, S. R. Brant, J. R. Gandler, G. Fendrich and C. Nakamura, J. Am. Chem. Soc., 1982, 104, 7045–7051; (b) W. v. E. Doering and E. A. Barsa, J. Am. Chem. Soc., 2004, 126, 12353–12362; (c) B. K. Carpenter, Chem. Rev., 2013, 113, 7265–7286; (d) Y. Yoneda, S. J. Mora, J. Shee, B. L. Wadsworth, E. A. Arsenault, D.- Hait, G. Kodis, D. Gust, G. F. Moore, A. L. Moore, M. Head-Gordon, T. A. Moore and G. R. Fleming, J. Am. Chem. Soc., 2021, 143, 3104–3112.
16. These reactions have been referred to by some authors as symmetric. This term is also in use, however, for identity reactions in the narrower sense. To clearly differentiate these two related but distinct concepts, and for clarity distinguishing function symmetry from reaction symmetry, we propose the use of isodromic instead.
17. Marcus theory, assuming both reactants and products distort identically and quadratically, approximates local isodromic behaviour. This leads to considering the local slope around zero driving force to be fixed at 0.5, which often does not hold. See the SI for an example of the consequences with an anisodromic reaction.
18. M. Solà, J. Mestres, R. Carbó and M. Duran, J. Am. Chem. Soc., 1994, 116, 5909–5915.
19. See further details on the derivation in the SI.
20. H. Feng, K. Zhang and J. Wang, Chem. Sci., 2014, 5, 3761–3769.
21. (a) R. G. Pearson and R. L. Dillon, J. Am. Chem. Soc., 1953, 75, 2439–2443; (b) N. Agmon, J. Am. Chem. Soc., 1980, 102, 2164–2167.
22. Selected examples: (a) G. L. Closs and J. R. Miller, Science, 1988, 240, 440–447; (b) M. Kuss-Petermann and O. S. Wenger, Phys. Chem. Chem. Phys., 2016, 18, 18657–18664; (c) J. Wang, T. Ding, K. Gao, L. Wang, P. Zhou and K. Wu, Nat. Commun., 2021, 12, 6333.
23. (a) M. Mandal, C. E. Elwell, C. J. Bouchey, T. J. Zerk, B. W. Tolman and C. J. Cramer, J. Am. Chem. Soc., 2019, 141, 17236–17244; (b) M. K. Goetz and J. S. Anderson, J. Am. Chem. Soc., 2019, 141, 4051–4062; (c) Y. Yang, X. Zhang, L.-P. Zhong, J. Lan, X. Li, C.-C. Li and L. W. Chung, Nat. Commun., 2020, 11, 1850.
24. K. Ando, Y. Shimazu, N. Seki and H. Yamataka, J. Org. Chem., 2011, 76, 3937–3945.
25. (a) A. J. Kresge, Can. J. Chem., 1974, 52, 1897–1903; (b) F. G. Bordwell, W. J. Boyle Jr., J. A. Hautala and K. C. Yee, J. Am. Chem. Soc., 1969, 91, 4002–4003; (c) F. G. Bordwell and W. J. Boyle Jr., J. Am. Chem. Soc., 1972, 94, 3907–3911.
26. X. Wang, S. Peter, M. Kinne, M. Hofrichter and J. T. Groves, J. Am. Chem. Soc., 2012, 134, 12897–12900.
27. C. C. Lee and E. C. F. Ko, J. Org. Chem., 1975, 40, 2132–2133.
28. F. Neese, Wiley Interdiscip. Rev.:Comput. Mol. Sci., 2022, 12(1), e1606.
29. S. Grimme, J. Antony, S. Ehrlich and H. Krieg, J. Chem. Phys., 2010, 132, 154104.
30. Y.-S. Lin, G.-D. Li, S.-P. Mao and J.-D. Chai, J. Chem. Theory Comput., 2013, 9, 263–272.
31. F. Weigend and R. Ahlrichs, Phys. Chem. Chem. Phys., 2005, 7, 3297–3305.
32. (a) F. Neese, J. Comp. Chem., 2003, 24, 1740–1747; (b) G. L. Stoychev, A. Auer and F. Neese, J. Chem. Theory Comput., 2017, 13, 554–562.
33. See the SI for details and additional calculations.
34. As the curvature was shallow, the inclusion of a quadratic term was confirmed to be statistically significant (p < 7.10⁻⁴). Likewise, analysis of covariance partitioning the values in two groups, ΔE_r < 0 and ΔE_r ≥ 0, showed statistical significance for the different slopes of each group (p < 2.10-4), of 0.1821 and 0.4902, respectively.
35. D. W. Flaherty and E. Iglesia, J. Am. Chem. Soc., 2013, 135, 18586–18599.
36. M. Xu, S. Jing, B. Sun, C. He, J. N. H. Reek and C. Duan, Angew. Chem., Int. Ed., 2023, 135, e202310420.
37. (a) W. v. E. Doering, J. L. Ekmanis, K. D. Belfield, F.-G. Klärner and B. Krawczyk, J. Am. Chem. Soc., 2001, 123, 5532–5541; (b) M. E. Abbasov, B. M. Hudson, D. J. Tantillo and D. Romo, Chem. Sci., 2017, 8, 1511–1524.
38. D. J. Tantillo, J. Am. Chem. Soc., 2022, 144, 13996–14004.
39. E. Beckmann and B. Dtsch, Chem. Ges., 1886, 19, 988–993.
40. (a) C. Bernasconi, Acc. Chem. Res., 1987, 20, 301–308; (b) C. Bernasconi, Acc. Chem. Res., 1992, 25, 9–16; (c) Z. Zhong, T. S. Snowden, M. D. Best and E. V. Anslyn, J. Am. Chem. Soc., 2004, 126, 3488–3495; (d) F. Liu, Z. Yang, Y. Yu, Y. Mei and K. N. Houk, J. Am. Chem. Soc., 2017, 139, 16650–16656.
41. M. Álvarez-Moreno, C. de Graaf, N. Lopez, F. Maseras, J. M. Poblet and C. Bo, J. Chem. Inf. Model., 2015, 55, 95–103.

---