Exploring the Mechanism of Catalysis with the Unified Reaction Valley Approach (URVA)—A Review

Abstract

The unified reaction valley approach (URVA) differs from mainstream mechanistic studies, as it describes a chemical reaction via the reaction path and the surrounding reaction valley on the potential energy surface from the van der Waals region to the transition state and far out into the exit channel, where the products are located. The key feature of URVA is the focus on the curving of the reaction path. Moving along the reaction path, any electronic structure change of the reacting molecules is registered by a change in their normal vibrational modes and their coupling with the path, which recovers the curvature of the reaction path. This leads to a unique curvature profile for each chemical reaction with curvature minima reflecting minimal change and curvature maxima, the location of important chemical events such as bond breaking/forming, charge polarization and transfer, rehybridization, etc. A unique decomposition of the path curvature into internal coordinate components provides comprehensive insights into the origins of the chemical changes taking place. After presenting the theoretical background of URVA, we discuss its application to four diverse catalytic processes: (i) the Rh catalyzed methanol carbonylation—the Monsanto process; (ii) the Sharpless epoxidation of allylic alcohols—transition to heterogenous catalysis; (iii) Au(I) assisted [3,3]-sigmatropic rearrangement of allyl acetate; and (iv) the Bacillus subtilis chorismate mutase catalyzed Claisen rearrangement—and show how URVA leads to a new protocol for fine-tuning of existing catalysts and the design of new efficient and eco-friendly catalysts. At the end of this article the pURVA software is introduced. The overall goal of this article is to introduce to the chemical community a new protocol for fine-tuning existing catalytic reactions while aiding in the design of modern and environmentally friendly catalysts.

1. Introduction

More than 180 years ago the Swedish chemist Berzelius defined a catalyst as a substance that increases the rate of a chemical reaction [1,2]. Today, more than 90% of all chemical processes involve at least at one point a catalytic transformation [3,4,5,6]. Prominent examples are the Ru/Os-catalyzed hydrogenation and oxidation reactions [7], the Ru/Mo-catalyzed olefin metathesis [8], and Pd-catalyzed cross-coupling reactions [9] which were awarded with the Nobel Prize in Chemistry in 2001, 2005, and 2010, respectively. It is interesting to note that 20 Chemistry Nobel Prizes have been associated with discoveries related to catalysts or catalysis [10]. According to a current report by Grand View Research, Inc., the global catalyst market size is expected to reach USD 34.1 billion by 2025 [11]. This market growth is related mainly to rising demand for polymers from automotive, electrical, and electronic industries. The growing area of nanocatalysis is an excellent example for so-called business to business (B2B) markets, leading to a new generation of patented nanocatalysts [12]. Another promising area is related to single-molecule catalysis monitored by fluorescence microscopy [13,14]. Biocatalysis is another expanding field of catalysis, which is inspired by enzymatic processes developed by nature [15,16]. Biocatalysts have become an essential tool in industrial biotechnology aimed at producing sustainably cheaper chemical products in chemical and pharmaceutical industries in an eco-friendly way [17].

These developments go hand-in-hand with the urgent need to optimize existing catalysts and to design new ones, a challenging endeavor for computational chemistry. A lot of effort has been put into the modeling of catalysis and their catalytic activities, ranging from the design and optimization of catalysts based on quantum chemical reaction modeling [18,19,20,21,22,23,24] and automated in silico catalyst design [25] to machine learning based procedures [26].

Modeling the mechanism of chemical reactions, in particular, complex catalytic cycles, at the atomic level is a formidable task, caused by high dimensionality of the potential energy surface (PES) which has to be mapped and analyzed for this purpose. Over the past few years, different methods have been explored to construct the PES, even for larger systems aiming toward chemical accuracy [27,28,29,30]; recently machine learning and artificial intelligence have been used [29,31,32,33]. Related approaches have been based on the fact that not all locations on the PES are of equal importance for a chemical reaction; examples include procedures mapping parts of the PES along reaction pathways [34,35,36,37,38,39], algorithms for automatically finding reaction paths [39,40,41,42,43], and methods working in a reduced dimensional space [44]. Numerous efforts have been devoted to the analysis of the information collected while tracing reaction paths on PESs. Some representative examples include: (i) The analysis of local reactivity descriptors, e.g., the Fukui function, and how they infer reactivity trends along the reaction path [45,46]. (ii) The reaction force and reaction force constant formalism, which is aimed at extracting information from higher derivatives of the energy profile taken along the reaction path [47,48,49], wherein one has to consider that the energy is a cumulative property, and as such the mechanistic information obtained is of more holistic nature. (iii) Monitoring changes of the topological features of the electron density [50,51] or the electron localization function [52] along the reaction path [53,54,55]. (iv) The recently suggested exploration of the PES with immersive virtual reality [56].

A complementary way of exploring the reaction mechanism is pursued in molecular dynamics (MD) simulations that allow an atomic visualization into the detail of molecular reactions on a femtosecond time scale [57,58,59,60]. MD simulations can be performed on a previously constructed PES or alternatively with a direct methodology calculating the trajectories on the fly [61,62,63]. In contrast to classical MD relying on Newton’s equation, ab initio MD is based on the Schrödinger equation. It offers a more realistic simulation of complex molecular systems and processes from first principles [64,65]. Recently, MD has also been paired with machine learning [32,66,67].

The unified reaction valley approach (URVA) [68,69,70,71,72] offers a different opportunity of gaining detailed mechanistic insights into chemical reaction mechanisms, perfectly complementing the above mentioned methodologies. URVA monitors all electronic changes in the reaction complex (RC; i.e., the union of the reacting molecules), from the pre-chemical to the post-chemical stages of the reaction based on the curving of the reaction path (RP), which the RC traces out on PES during the reaction. With regard to catalysis, URVA specifies and quantifies different catalytic activities for each reaction of a catalytic cycle and clarifies how catalysis can be related to (un)known electronic structure principles. We will present in this work the underlying concepts of URVA, followed by the application of URVA to a variety of catalytic reactions. We will conclude our article by discussing how, based on our findings, new concepts for rational catalyst design can emerge, including aspects of both chemo- and bio-catalysis, thereby serving as an incubator for enhanced interdisciplinary efforts in the design of the next generation of catalysts. In the Appendix A, the new pURVA software will be introduced.

2. The Unified Reaction Valley Approach (URVA)

2.1. Background: The Reaction Path Hamiltonian

In 1980, Miller, Handy, and Adams introduced in their seminal paper the reaction path Hamiltonian (RPH) that describes a chemical reaction in a specific region of the PES, the so-called reaction valley stretching around the RP connecting reactants, TS, and products [73]. In the same year, similar work of Kato and Morokuma [74] was published. The underlying idea was (i) the reaction valley is the location where all important mechanistic and dynamic steps of a chemical reaction take place; and (ii) the RC can be described via a large amplitude motion along the RP and the vibrations perpendicular to it, which has its roots in the early work of Hofacker [75], Marcus [76,77,78], and Hougen [79].

The RPH is a classical Hamiltonian H describing the kinetic energy T and potential energy V of the RC when moving along the RP. As the reaction path one generally chooses the steepest descent path in mass-weighted Cartesian coordinates, leading from the TS into reactant and product channels; the so-called intrinsic reaction coordinate (IRC) path, first introduced by Eliason and Hirschfelder [80], and formalized by Fukui [81].

$$H[s,p_s,\left\{Q_{\mu }\right\},P_{\mu }]=T[s,p_s,\left\{Q_{\mu }\right\},P_{\mu }]+V(s,Q_{\mu })$$

(1)

Considering a RC of N atoms and $N_{vib}=3N-L$ internal coordinates (L = 6 for non-linear and 5 for linear molecules), the $N_{vib}$ dimensional configurational space of H is partitioned into the one-dimensional reaction parameter s (i.e., the arc length of the RP) and its conjugated momentum $p_s$ describing the motion along the RP, and a set of normal coordinates $Q_{\mu }$ and conjugated momenta $P_{\mu }$ ($\mu$ = 2, ⋯, $N_{vib})$ forming the reaction valley; i.e., the valley walls perpendicular to the $\mu$ = 1 direction along the RP.

The potential $V[s,Q_{\mu }]$ is approximated at each point s by the potential $V_0\left(s\right)$ at s plus the potential for the harmonic displacements perpendicular to the path, defining the curvature of the valley; i.e., its steepness

$$V[s,Q_{\mu }]=V_0\left(s\right)+1/2\sum _{\mu =2}^{N_{vib}}{\left(k_{\mu }^g\right)}^2\left(s\right)Q_{\mu }^2\left(s\right)$$

(2)

with $k_{\mu }^g\left(s\right)$ being the generalized force constant. It has to be noted that instead of a standard normal mode analysis performed at a stationary point, a generalized normal mode analysis [73,82] is necessary, which projects out at each path point s in addition to the rotations and translations of the reaction complex also the translational motion in internal coordinates along the RP from the mass-weighted Hessian matrix ${\tilde{\mathbf{f}}}^x\left(s\right)$ with a suitable projection operator $\mathbf{P}\left(s\right)$ [83]. This leads to the generalized Hessian ${\tilde{\mathbf{f}}}^g\left(s\right)$

$${\tilde{\mathbf{f}}}^g\left(s\right)=\left(\mathbf{I}-\mathbf{P}\left(s\right)\right){\tilde{\mathbf{f}}}^x\left(s\right)\left(\mathbf{I}-\mathbf{P}\left(s\right)\right)$$

(3)

The generalized mass-weighted generalized normal vibrational modes ${\tilde{\ell }}^g\left(s\right)$ defined in the ($N_{vib}-1$) dimensional reaction valley space and corresponding frequencies ${\omega }_{\mu }^g\left(s\right)$ can then be determined by solving a generalized Wilson equation [82,83]

$${\tilde{\mathbf{f}}}^g\left(s\right){\tilde{\ell }}_{\mu }^g\left(s\right)={\left(2\pi c{\omega }_{\mu }^g\left(s\right)\right)}^2{\tilde{\ell }}_{\mu }^g\left(s\right)$$

(4)

where constant c is the speed of light.

The kinetic energy $T[s,p_s,\left\{Q_{\mu }\right\},P_{\mu }]$ is composed of two terms; the first term corresponds to the kinetic energy part of the movement along the RP and the second term represents the kinetic energy part associated with a movement orthogonal to the path direction.

$$T[s,p_s,\left\{Q_{\mu }\right\},P_{\mu }]=\frac{1}{2}\frac{{[p_s-{\sum }_{\mu ,\nu =2}^{N_{vib}}{B}_{\mu \nu }\left(s\right)Q_{\mu }\left(s\right)P_{\nu }\left(s\right)]}^2}{{[1+{\sum }_{\mu =2}^{N_{vib}}{B}_{\mu s}\left(s\right)Q_{\mu }\left(s\right)]}^2}+\frac{1}{2}\sum _{\mu =2}^{N_{vib}}{P}_{\mu }{\left(s\right)}^2$$

(5)

Central to the first part of the kinetic energy are the so-called coupling coefficients $B_{\mu \nu }\left(s\right)$ and $B_{\mu s}\left(s\right)$. The coupling term $B_{\mu \nu }\left(s\right)$ reflects the energy transfer between modes $\mu$ and $\nu$ induced by the motion along the RP. As this motion involves a rotation of the transverse vibrational modes about the RP, they are often referred to as Coriolis couplings. The $B_{\mu \nu }\left(s\right)$ terms are given by the dot product between the normal mode vector ${\tilde{\ell }}_{\mu }^g\left(s\right)$ and the change of normal mode vector ${\tilde{\ell }}_{\nu }^g\left(s\right)$ with regard to s and vice versa [73]:

$$B_{\mu \nu }\left(s\right)={\tilde{\ell }}_{\nu }^g{\left(s\right)}^{†}\frac{d{\tilde{\ell }}_{\mu }^g\left(s\right)}{ds}=-B_{\nu \mu }\left(s\right)$$

(6)

The coupling term $B_{\mu s}\left(s\right)$ reflects the translational–vibrational energy transfer between the translational motion along the RP and vibrational mode $\mu$ [73] induced by the curving of the RP. Since the RP is a curved line in $N_{vib}$ dimensional space, its direction and curvature can be derived with the Frenet–Serret formalism [84]. The reaction path direction at a path point s is given by the unit vector $\eta \left(s\right)$:

$$\mathbf{\eta }\left(s\right)=\frac{d\tilde{\mathbf{x}}\left(s\right)}{d\mathit{s}}=-\frac{\tilde{\mathbf{g}}\left(\tilde{\mathbf{x}}\left(s\right)\right)}{\parallel \tilde{\mathbf{g}}\left(\tilde{\mathbf{x}}\left(s\right)\right)\parallel }$$

(7)

where the derivative of the mass-weighted reaction coordinate $\tilde{\mathbf{x}}\left(s\right)$ with regard to s is the normalized mass-weighted gradient vector $\tilde{\mathbf{g}}\left(s\right)\equiv \tilde{\mathbf{g}}\left(\tilde{\mathbf{x}}\left(s\right)\right)={\mathbf{M}}^{1/2}\mathbf{g}\left(s\right)$ and $\mathbf{M}$ is a diagonal matrix of atomic masses. The curvature vector $\mathit{\kappa }\left(s\right)$ is given by [68,83]

$$\mathit{\kappa }\left(s\right)=\frac{d^2\tilde{\mathbf{x}}\left(s\right)}{d{\mathit{s}}^2}=\frac{d\mathit{\eta }\left(s\right)}{d\mathit{s}}=\frac{-1}{\parallel \tilde{\mathbf{g}}\left(s\right)\parallel }\left({\tilde{\mathbf{f}}}^x\left(s\right)\eta \left(s\right)-\left[{\left(\mathit{\eta }\left(s\right)\right)}^{†}{\tilde{\mathbf{f}}}^x\left(s\right)\mathit{\eta }\left(s\right)\right]\mathit{\eta }\left(s\right)\right)$$

(8)

The $B_{\mu s}\left(s\right)$ can then be calculated as the dot product between the reaction path vector $\mathit{\eta }\left(s\right)$ and the change of normal mode vector ${\tilde{\ell }}_{\mu }^g\left(s\right)$ equivalent to the dot product of ${\tilde{\ell }}_{\mu }^g\left(s\right)$ and the change of $\mathit{\eta }\left(s\right)$ with regard to s, which corresponds to the reaction path curvature $\mathit{\kappa }\left(s\right)$

$$B_{\mu s}\left(s\right)=\mathit{\eta }{\left(s\right)}^{†}\frac{d{\tilde{\ell }}_{\mu }^g\left(s\right)}{ds}=-{\tilde{\ell }}_{\mu }^g{\left(s\right)}^{†}\frac{d{\mathit{\eta }}_{\mu }\left(s\right)}{ds}={\tilde{\ell }}_{\mu }^g{\left(s\right)}^{†}\mathit{\kappa }\left(s\right)$$

(9)

As $B_{\mu s}\left(s\right)$ becomes larger, the mixing between the RP and vibrational mode $\mu$ increases and the translational–vibrational energy transfer becomes more pronounced. This can be directly monitored by the scalar curvature $\kappa \left(s\right)=\left|\right|\mathit{\kappa }\left(s\right)\left|\right|$ which is defined as [73]

$$\kappa \left(s\right)={\left[\mathit{\kappa }\left(s\right)\mathit{\kappa }\left(s\right)\right]}^{1/2}={\left[\sum _{\mu =2}^{N_{vib}}{B}_{\mu s}^2\left(s\right)\right]}^{1/2}$$

(10)

Because of this relationship, the $B_{\mu s}\left(s\right)$ coefficients are called curvature couplings.

The main focus of the RPH was and still is to be used as a tool for the calculation of the dynamics of a chemical reaction, and in particular, the calculation of rate constants and tunneling coefficients [85,86]; or, as a valuable resource for laser spectroscopists working in the field of vibrationally driven reactions, which includes enhancement of reaction rates, manipulation of energy disposal, and promotion of a certain product channel by mode selective excitation [87,88,89,90,91,92]. However, the depth of mechanistic information provided by the RPH was not fully exploited in a systematic way, until Kraka, Cremer, and co-workers started to transform the RPH approach into an advanced mechanistic tool, known as the unified reaction valley approach (URVA) [68,69,70,93], which will be discussed in the following.

2.2. Basic Methodology of URVA

During the course of a chemical reaction, the RC changes its electronic structure. This is directly registered by the vibrational modes, which are sensitive to even the smallest electronic structure changes. The change in the vibrations leads to a change in the coupling between valley and path motions as is described by $B_{\mu s}\left(s\right)$ coefficients, which altogether define the scalar reaction path curvature, as shown in Equation (10). Therefore, URVA’s main focus is on the scalar curvature $\kappa \left(s\right)$. Monitoring the $\kappa \left(s\right)$ along s leads to a unique curvature profile for each chemical reaction, with curvature maxima K and minima M as schematically shown in Figure 1a. The curvature maxima define the locations of electronic structure change, such as charge transfer and charge polarization, rehybridization, bond cleavage, and formation, and also change in the optimal orientation of the reactants for the reactive collision. Each curvature peak corresponding to a certain chemical event is flanked by two curvature minima, e.g., locations with minimal change of the RC, which can be classified as “hidden intermediates” [71,72,94,95]. Accordingly, Kraka and Cremer [69,71] defined a reaction phase as the reaction path region from one curvature minimum to the next, characterized by a curvature maximum. Different chemical reactions possess different curvature patterns with different numbers of reaction phases, which can be used as their fingerprints.

Further insight into the reaction mechanism, e.g., answering the question of which bond is broken/formed during a certain reaction phase, requires a decomposition of $\kappa \left(s\right)$. In the original RPH of Miller, Handy, and Adams, a decomposition into $B_{\mu s}\left(s\right)$ coefficients was performed. However, for a mechanistic analysis, the coupling coefficients $B_{\mu s}\left(s\right)$ are of limited use because they refer to generalized normal vibrational modes, as shown in Equation (6). Normal vibrational modes are generally delocalized over part or even the entirety of a molecule because of mode–mode coupling [82]. Therefore, the first milestone of the URVA analysis was to introduce the curvature decomposition into local curvature coupling coefficients $A_{ns}^{\kappa }$ based on local vibrational modes ${\mathbf{a}}_n$ that are associated with the internal coordinates q${}_n\left(s\right)$ used to describe the RC [68]. The local vibrational mode theory, originally introduced by Konkoli and Cremer [96,97,98,99,100,101,102] derives local vibrational modes ${\mathbf{a}}_n$ via mass-decoupled Euler Lagrange equations, which then can be used to define the local curvature coupling coefficients $A_{ns}^{\kappa }$ for the local decomposition of the reaction path curvature [68,69]

$$\begin{array}{ccc}\hfill A_{ns}^{\kappa }\left(s\right)& =& \mathit{\kappa }{\left(s\right)}^{†}\frac{{\mathbf{a}}_n^g\left(s\right)}{\parallel {\mathbf{a}}_n^g\left(s\right)\parallel }\hfill \end{array}$$

(11)

In this way, each curvature peak can be associated with a particular chemical event, as sketched in Figure 1b. Often, only a few vibrational modes at a given path position s contribute to the curving of the reaction path, which makes the analysis of larger chemical reactions feasible [68,95,103,104,105,106,107,108,109,110]. It is noteworthy that the sign of the component denotes whether it supports the chemical change (positive sign) or resists (negative sign).

The curvature decomposition into local modes ${\mathbf{a}}_n^g\left(s\right)$ via Equation (11) has been successfully applied to a number of organic reactions, such as the hydrogenation of methyl radical [68,105], FH addition to ethylene [95], 1,2-H-shift in methylchlorocarbene [103], the Diels–Alder reaction [104,106], the cycloaddition of 1,3 dipoles and ethylene or acetylene [107,108], and other pericyclic reactions [109]. However, it fails if the description of the reaction path is hampered by path instabilities (reflected by the occurrence of imaginary reaction valley frequencies). These can be caused by methodological shortcomings (e.g., in the case of bond-breaking or bond-forming processes leading to multi-reference character, which is often encountered for catalytic reactions involving transition metals) or can have a chemical origin (e.g., the reaction path bifurcates as an intrinsic feature of the PES). Any path instability precludes the description of electronic structure changes in terms of local mode curvature coupling coefficients $A_{ns}^{\kappa }$ being based on local modes ${\mathbf{a}}_n^g\left(s\right)$. Therefore, the decomposition of the reaction path curvature $\kappa \left(s\right)$ in terms of geometrically-based local modes ${\mathbf{u}}_n$ (dubbed internal coordinate component analysis) was developed, thereby offering a robust reaction path analysis which is no longer sensitive to path instabilities [110,111,112,113]. This can be considered as the second milestone of URVA that lead to a breakthrough in particular for the investigation catalysis reactions, which will be demonstrated below.

For each internal coordinate $q_n$, a unit column vector ${\mathbf{u}}_n$ can be defined [101,111] through its local mass $m_n^q=G_{n,n}^{-1}$ [96,97] and Wilson’s B-matrix formalism [82], which connects internal coordinate $q_n$ to the Cartesian coordinates $\mathbf{x}$ via ${\mathbf{b}}_n=d{q}_n/d\mathbf{x}$,

$$\begin{array}{c}\hfill {\mathbf{u}}_n=\frac{{\mathbf{M}}^{-1/2}{\mathbf{b}}_n^{†}}{\parallel {\mathbf{M}}^{-1/2}{\mathbf{b}}_n^{†}\parallel }=G_{n,n}^{-1/2}\left({\mathbf{M}}^{-1/2}{\mathbf{b}}_n^{†}\right)\end{array}$$

(12)

where $G_{n,n}={\mathbf{b}}_n{\mathbf{M}}^{-1}{\mathbf{b}}_n^{†}$. With the help of ${\mathbf{u}}_n$, Equation (7) can be rewritten in the mass-weighted internal coordinate ${\tilde{q}}_n={\left(m_n^q\right)}^{1/2}{q}_n$ [111],

$$\begin{array}{ccc}\hfill {\eta }_n^q\left(s\right)=\frac{d{\tilde{q}}_n\left(s\right)}{d\mathit{s}}& =& G_{n,n}^{-1/2}\frac{d{q}_n\left(s\right)}{d\tilde{\mathbf{x}}\left(s\right)}\frac{d\tilde{\mathbf{x}}\left(s\right)}{d\mathit{s}}=G_{n,n}^{-1/2}{\mathbf{b}}_n\left(s\right){\mathbf{M}}^{-1/2}\mathit{\eta }\left(s\right)={\mathbf{u}}_n^{†}\left(s\right)\mathit{\eta }\left(s\right)\hfill \end{array}$$

(13)

leading to a decomposition of the reaction path direction $\eta \left(s\right)$ into internal coordinate components.

Similarly, Equation (8) may also be rewritten in ${\tilde{q}}_n$ using Equation (13) [111]

$$\begin{array}{ccc}\hfill \frac{d^2{\tilde{q}}_n\left(s\right)}{d{s}^2}& =& \frac{d{\eta }_n^q\left(s\right)}{ds}=\frac{d}{ds}\left({\mathbf{u}}_n^{†}\left(s\right)\mathbf{\eta }\left(s\right)\right)={\mathbf{u}}_n^{†}\left(s\right)\mathit{\kappa }\left(s\right)+\frac{d{\mathbf{u}}_n^{†}\left(s\right)}{ds}\mathbf{\eta }\left(s\right)\hfill \end{array}$$

(14)

$$\begin{array}{cc}=& \parallel \mathit{\kappa }\left(s\right)\parallel cos{\beta }_n\left(s\right)+\parallel \frac{d{\mathbf{u}}_n\left(s\right)}{ds}\parallel cos{\gamma }_n\left(s\right)={\kappa }_n^q\left(s\right)+{\eta }_n^u\left(s\right)\hfill \end{array}$$

(15)

The first term on the right hand side of Equation (15) is the projection of u${}_n\left(s\right)$ onto the curvature vector $\mathit{\kappa }\left(s\right)$ leading to the amplitude $\parallel {\mathbf{u}}_n\left(s\right)\parallel cos{\beta }_n\left(s\right)=cos{\beta }_n\left(s\right)={\kappa }_n^q\left(s\right)$, as each ${\mathbf{u}}_n$ is a unit vector describing the local motion, which is driven by the internal coordinate $q_n\left(s\right)$. The amplitude is scaled by the scalar curvature $\parallel \mathit{\kappa }\left(s\right)\parallel$ corresponding to the length of the curvature vector in $N_{vib}$-dimensional space. The mixed second order term ${\eta }_n^u\left(s\right)$ is determined by the change in the direction of ${\mathbf{u}}_n\left(s\right)$ with s and the tangent vector. The vector derivative du${}_n\left(s\right)$/ds that is orthogonal to ${\mathbf{u}}_n\left(s\right)$ is projected onto vector $\mathit{\eta }$ and therefore does not have any information on the curvature [111]. Hence, the curvature contribution of the internal coordinate $q_n$ is defined by the first term only; i.e.,

$$\begin{array}{c}\hfill {\kappa }_n^q\left(s\right)={\mathbf{u}}_n^{†}\left(s\right)\kappa \left(s\right)\end{array}$$

(16)

It has to be noted that the decomposition into any coordinate, e.g., Cremer–Pople puckering coordinates [114], is possible, as long as the corresponding B-matrix vector can be derived [106,109].

Other important features of URVA include the advanced path following procedure of Hratchian and Kraka [115], which allows a chemical reaction to be tracked far into entrance and exit channels. This is where decisive events often happen [106], particularly in the case of catalysis, as will be discussed below. URVA is not restricted to the IRC. It can be applied to any representative path. For example, in the case of barrier-less reactions without a TS for which an IRC does not exist, a representative path based on Newton trajectories can be used [94,116,117]. This led to the surprising result that barrier-less reactions such as the chelotropic reaction between methylene and ethene [94] or the spontaneous cycloaddition between vinylidene and acetylene reaction often possess a complex reaction mechanism [117]. It is important to note that the use of mass-weighted coordinates allows for a one-to-one comparison of reaction and curvature profiles in different media, i.e., gas phase versus solution or enzyme environment described by a QM/MM (quantum mechanics/molecular mechanics) approach [118]; and of reactions with different substituents. For example, the URVA study of 20 different 1,3 dipolar cycloadditions elucidated different substituent effects and the exchange of reaction partners in a compact and concise form [107,108].

The URVA study involving a multitude of catalyzed reactions in comparison with their non-catalyzed counterparts led to a number of interesting mechanistic insights which will be discussed in the following representative examples: (i) Rh-catalyzed methanol carbonylation; (ii) Sharpless epoxidation of of allylic alcohols; (iii) Au(I) assisted [3,3]-sigmatropic rearrangement of allyl acetate; and (iv) Bacillus subtilis chorismate mutase catalyzed Claisen rearrangement.

3. Computational Methods

The calculations of the methanol carbonylation catalyzed by dicarbonyldiiodorhodium [Rh(CO)${}_2$I${}_2$]${}^{-}$ were performed with the B3LYP density functional [119,120] utilizing Pople’s 6-31G(d,p) basis set [121,122] and the Stuttgart–Dresden (SDD) effective core potential [123] with the corresponding basis set for the Rh and I atoms [124,125]. The calculations of the Sharpless epoxidation reaction of allylic alcohols catalyzed by [Ti(tartrate)(OR)${}_2$]${}_2$ were carried out with the B3LYP/6-31G(d,p) level of theory and utilizing the SDD effective core potential with a corresponding basis set for the Ti atom [126]. The calculation of the non-catalyzed [3,3]–sigmatropic rearrangement of allyl acetate was performed with B3LYP and Pople’s 6-31+G(d,p) basis set [127]. The Au(I) catalyzed reactions were performed at the same level of theory utilizing using the SDD effective core potential with the corresponding basis set for the Au atom [123,124]. The corresponding energetics were calculated at the DLPNO-CCSD(T) level of theory [128,129] for the DFT geometries, using for Au the same basis set as for the DFT calculations, and the cc-pVQZ Dunning basis set for all other atoms [130,131,132]. The IRC calculations of the Claisen rearrangement of chorismate to prephenate in the gas phase (R1) were performed using the B3LYP/6-31+G(d,p) level of theory. The calculations in implicit aqueous solution (R2) were performed at the same level of theory as in the gas phase using the polarizable continuum solvent model (PCM) of Tomasi [133]. The calculations in explicit water, reaction R3, were performed with a hybrid of QM/MM (quantum mechanics/molecular mechanics) as provided by the ONIOM methodology [134]. The QM part was described at the B3LYP/6-31+G(d,p) level of theory; water was modeled with TIP3P [135]. The enzyme reactions R4–R6 were also performed with QM/MM; for the MM part the AMBER General Force Field for organic molecules was utilized [136,137]. Technical details for the QM/MM calculations are provided in the Supplementary Materials.

All IRC calculations were performed with Gaussian09 [138] applying the improved reaction path following procedure of Hratchian and Kraka, enabling the tracking of a chemical reaction far out into the entrance and exit channels [115]. URVA and the local mode analysis were carried out with the COLOGNE20 program [139]. DLPNO-CCSD(T) calculations were executed with the ORCA 4.0.1 program [128]. Natural bond orbital (NBO) charges were calculated with the program NBO 6 [140]. A standalone URVA program written in python (pURVA), is described in the Appendix A.

4. Results

4.1. Rh Catalyzed Methanol Carbonylation—The Monsanto Process

Acetic acid is an important commodity chemical used in a broad range of applications [141]. In 2014, the worldwide demand was 12.1 million tons; projections estimate at least 17 million tons by 2022 [142]. Originally, acetic acid was produced by aerobic fermentation of ethanol, which is still the major process for the production of vinegar [143]. In 1970, Monsanto commercialized an improved homogeneous methanol carbonylation process using a methyl-iodide-promoted Rh catalyst [144,145]. Compared to ethanol fermentation, homogeneous Rh catalyzed methanol carbonylation is an efficient route that exhibits high productivity and yields. However given the reaction conditions (pressure: 30–60 bar; temperature: 150–200 ${}^{\circ }$C) there is room for improvement [141,146,147]. As shown in the center of Figure 2, the catalytic cycle involving the original Monsanto catalyst [RhI${}_2$(CO)${}_2$]${}^{-}$ involves four steps: step 1, oxidative addition; step 2, methyl group migration; step 3, carbon monoxide insertion; and step 4, reductive elimination of product and catalyst regeneration [144,145]. The corresponding curvature profiles and their decomposition into internal coordinate components provide direct access to the most important chemical events taking place in each catalysis step, which will be discussed in the following.

Step 1 is characterized by a longer preparation phase 1 aligning the reaction partners. In phase 2, a smaller curvature peak develops caused by the beginning of the formation of the new Rh–C bond (blue line, resisting) and the formation of the new Rh–I bond (red line, supportive). In phase 3 with a more pronounced curvature peak, both Rh–I and Rh–C bond formation continue and C–I bond cleavage of the CH${}_3$I reactant starts (all supportive). Three observations are noteworthy which are typical for catalysis reactions: (i) several chemical events take place at the same time; (ii) they are often accompanied by a change of the coordination number of the metal atom—in this case the Rh coordination number changes from 4 to 6; and (iii) bond formation/breakage can stretch over several phases—here, Rh–I and Rh–C bond formation is finalized in phase 5, i.e., after the TS, where the reaction is already on a downhill path, and therefore does not require further energy to proceed. In step 2 migration of the methyl group from Rh to the one of the CO ligands takes place. This process is initiated by the bending of the equatorial iodine atom (I–Rh–C angle, green line) dominating phases 2 and 3. In phase 4 the breakage of the Rh–C bond (blue line) starts with a small contribution from the new C–C bond formation (orange line). The actual finalization of the bond breaking/forming process occurs simultaneously with a reorganization of the Rh–CO bond (purple line) after the TS. During this process, Rh changes its coordination number from 6 to 5. The following step (step 3) is characterized by CO insertion, a barrier-less process in which Rh restores its coordination number from 5 to 6. The final step (step 4) involves breakage of the Rh–C bond between Rh and the acetic acid group, migration of one iodine atom to the acetic acid group, i.e., breakage of one Rh–I bond, and formation of the new C-I bond. As depicted by the curvature decomposition, this complex process starts in phase 2 with a reorientation of the methyl hydrogens, as reflected by the torsional angle RhCCH (purple line, shown representatively for one of the methyl hydrogens) helping to initiate the breaking of the Rh–C bond (blue line) and the formation of the new C–I bond (green line, first resisting) which takes the lead in phase 4 before the TS. The breakage of the Rh–I bond (read line) contributes only marginally to the curvature peaks before the TS. All three bond breaking/forming processes are finalized synchronously in phase 4 after the TS, i.e., without further energy demand. Reaction movies for each steps can be found in the Supplementary Materials, Videos S1–S4.

The next step of the analysis is to combine mechanistic and energy information. As obvious from the data in

Table 1. Activation energy $\mathsf{\Delta }{\mathbf{E}}^{‡}$, reaction energy $\mathsf{\Delta }{\mathbf{E}}_{\mathbf{R}}$, activation enthalpy $\mathsf{\Delta }{\mathbf{H}}^{‡}$, and reaction enthalpy $\mathsf{\Delta }{\mathbf{H}}_{\mathbf{R}}$ for each step of the [Rh(CO)${}_2$I${}_2$]${}^{-}$ catalyzed methanol carbonylation ${}^a$. Energy values $\mathsf{\Delta }\mathbf{E}\left(\mathbf{m}\right)$ denote the energy consumed in phase m up to the TS. ${}^c$

Reaction Step	$\mathsf{\Delta }{\mathbf{E}}^{‡}$	$\mathsf{\Delta }{\mathbf{E}}_{\mathbf{R}}$	$\mathsf{\Delta }{\mathbf{H}}^{‡}$	$\mathsf{\Delta }{\mathbf{H}}_{\mathbf{R}}$	$\mathsf{\Delta }\mathbf{E}\left(1\right)$	$\mathsf{\Delta }\mathbf{E}\left(2\right)$	$\mathsf{\Delta }\mathbf{E}\left(3\right)$	$\mathsf{\Delta }\mathbf{E}\left(4\right)$
1	45.0	2.4	43.5	2.1	13.2	12.0	17.1	2.7
2	20.0	−9.3	19.1	−8.5	0.0	0.9	3.1	16.0
3 ${}^b$	-	−8.0	-	−7.5	-	-	-	-
4	29.9	0.5	25.0	0.1	0.5	6.5	12.7	6.2

${}^a$ B3LYP/6-31G(d,p)/SDD level of theory. All values (kcal/mol) are given relative to the first point in the entrance channel of the corresponding reaction step. ${}^b$ Downhill reaction without TS. ${}^c$$\mathsf{\Delta }\mathbf{E}\left(\mathbf{m}\right)$ of phase m is the energy difference between the energy of phase m and energy of phase m-1, $\mathsf{\Delta }\mathbf{E}\left(\mathbf{4}\right)$ is an energy difference between energy at the TS and the energy of phase 3 of the corresponding reaction step.

showing the energetics of each reaction step, step 1 is most demanding with an activation energy of 45.0 kcal/mol, whereas steps 2 and 4 have activation energies of 20.0 and 29.9 kcal/mol, whereas 3 proceeds without a barrier. Clearly, the high barrier of step 1 is the reason for the harsh reaction conditions. Further insight is needed to clarify which chemical events of this reaction step lead to such a high barrier. This can be gained by analyzing how much energy is consumed in each reaction phase up to the TS. As depicted in Table 1 approach of the reactants in step 1 requires 13.2 kcal/mol, and formation of the Rh–C bond 12.0 kcal/mol, whereas the cleavage of the C–I bond is, with 17.1 kcal/mol, the most energy demanding process. In step 2 the breakage of the Rh–C bond is with 16.0 kcal/mol most expensive, and in step 4 the breakage of the Rh–C bond with 12.7 kcal/mol. These findings provide valuable information for catalyst modification. We are currently investigating a selenium functionalized bis(phosphino)methane rhodium catalyst [Rh(CO)Cl(PH${}_2$CH${}_2$PH${}_2$(Se))] [148] which could help via chelating the reactant to shorten the long approach phase and via polarization of the C–I bond supporting C–I bond breakage in the first step, lowering in this way the activation energy of step 1—an important move toward milder reaction conditions.

4.2. Sharpless Epoxidation of Allylic Alcohols—Transition to Heterogenous Catalysis

In 1980, Katsuki and Sharpless reported the enantioselective synthesis of 2,3-epoxyalcohols in the presence of hydroperoxide, (mostly t-butylhydroperoxide), a titanium (IV)-tartrate-diester catalyst, and a secondary allylic alcohol [149]. Given the commercial availability and low cost of the reagents, the high enantioselectivity of the reaction, and the versatility of the resulting chiral epoxides, this method has become an essential tool in organic synthesis [150,151,152,153,154]. As sketched in Figure 3a the Ti-tartrate -diester catalyst [Ti(IV)(tartrate) (OR)${}_2$] forms a chiral dimer that controls the stereochemical outcome of the epoxidation process. The allylic alcohol binds to one of the Ti centers together with the hydrogen peroxide oxidant, as shown in Figure 3b, so that the metal center can support the breakage of the peroxide O–O bond and guide the subsequent transfer of the coordinated hydrogen peroxide oxygen atom to the ethylene unit of the allylic alcohol to form the two new CO epoxide bonds. After epoxide formation, Ti coordinates the cleaved OR part of the hydrogen peroxide.

We performed an URVA study of the Sharpless epoxidation reaction for a number of different epoxy alcohols to clarify the question of whether the curvature profile for such a complex system still discloses the basic reaction mechanism, i.e., uncovers the most important features of the reaction mechanism. In the following we will discuss the reaction for but-2-en-1-ol, methyl-hydroperoxide, and a dimeric Ti(IV) model catalyst, which was derived from the crystal structure of original Sharpless catalyst [154] with the simplifications shown in Figure 3b. A reaction movie can be found in the Supplementary Materials, Video S5.

As shown in Figure 4, we found seven distinct reaction phases. The calculated activation energy is 20.8 kcal/mol, which is almost 20 kcal/mol lower than that for the non-catalyzed reaction [155]. The reaction is exothermic by 42.2 kcal/mol; see Figure 4a. The corresponding curvature diagram shown in Figure 4b provides the mechanistic details. The reaction starts with the preparation phases 1 and 2, which are characterized by allyl group rotations and the orienting of the migrating oxygen atom O1 of the peroxide into a perfect position above the double bond, (see reaction Video S5, Supplementary Materials). Phase 3, the most energy consuming phase requiring 14.2 kcal/mol, is dominated by the breakage of the peroxide bond (red line). The formation of the new epoxide bonds starts (O1–C2, green line; O1–C1 blue line, both resisting) and there are also smaller contributions from TiO1 (orange line); i.e., the Ti–O bond to be broken and from TiO2 (olive line) the new Ti–O bond to be formed during the epoxidation process. In phase 4 including the TS, the TiO1 bond takes over the lead; it lengthens and O1 approaches the double bond. It is noteworthy that the curvature around the TS is small, confirming a general finding that in many reactions, including catalytic reactions, important chemical events take place before or after TS. (A prominent example is the Diels Alder reaction. The curvature profile of this reaction shows only one large curvature peak far out in the exit channel corresponding to the formation of the two new CC bonds—a signature of symmetry allowed reactions [104,106,111]). The synchronous finalization of the new C–O bonds occurs in phase 5, after the TS, i.e., without further energy request. This process is accompanied by the finalization of O–O bond breakage, which stretches into phase 6, the cleavage of TiO1, and the formation of TiO2 bonds. Phase 7 is devoted to allyl group rotation and reorientation of the expoxide, (see reaction Video S5, Supplementary Materials). Clearly, the catalytic activity is based on the topology of the dimeric Ti catalyst framework, channeling the hydroperoxide oxidant into a central position for the attack of the CC double bond of the allyl alcohol, which is coordinated in a way that allyl group rotations are feasible for optimal orientation.

One can say that the dimeric Ti catalyst mimics a surface typical of the heterogenous catalysis, thereby facilitating a stereospecific collision of the reaction partners. The peroxide O–O bond glides over the Ti atom. During this process, the metal atom polarizes the oxygen atoms, facilitating O–O bond breakage. Another important feature of the Sharpless reaction is that both new C–O epoxide bonds are synchronously finalized after the TS, i.e., without energy consumption; i.e., catalyst optimization should predominantly focus on improving O–O breakage. This example clearly shows that the curvature diagram unravels all mechanistic details, even for complex reactions, confirming that the chemical events are mostly characterized by a few significant parameters at a time.

4.3. Au(I) Assisted [3,3]-Sigmatropic Rearrangement of Allyl Acetate

Gold catalysis offers a large spectrum of organic transformations under mild conditions with high yields [156,157,158,159,160,161,162,163,164,165]. In particular, cationic gold complexes efficiently activate intra or intermolecularly, CC, CN, or CO multiple bonds for a nucleophilic attack [112,159,166,167,168], whereas Au(III) catalysts can be considered as hard Lewis acids with a higher oxophilic character, being useful when d-holes are needed [169,170]; Au(I) species are soft Lewis acids with high carbophilic character and high affinity to coordinate with carbon $\pi$-systems [166,167]. Au(I) phosphine complexes are primarily applied. However, Au(I) complexes with heterocyclic carbenes (NHC) ligands [171,172] are an interesting alternative; NHCs are easier to synthesize and to handle than their phosphine counterparts, because they are less likely to dissociate during a chemical reaction from the attached metal atom; they are stronger $\sigma$-donors than most electron rich phosphines, and they display a broader range of versatility and regio-selectivity [173,174,175,176].

Gold-catalyzed sigmatropic rearrangements have been successfully applied in gold catalysis over the past decade [177], such as the gold catalyzed Claisen rearrangement of propargylic esters, which undergo intramolecular 1,2-ester or 1,3-ester migration [168]. We elucidated the mechanistic differences of the non-catalyzed and the Au(I)-NHC catalyzed [3,3]-sigmatropic rearrangement of allyl acetate via a comprehensive URVA study, including the investigation of the influence of solvation, counterions, bulky, and electron withdrawing/donating substituents as well as the exchange of the Au(I)-NHC with a Au(I)-phosphine catalyst [110].

In agreement with previous theoretical studies [178], we could confirm that the non-catalyzed rearrangement reaction proceeds in one step (DLPNO-CCSD(T): activation energy $\mathsf{\Delta }$E${}^{\#}$ of 45.9 kcal/mol and a reaction energy $\mathsf{\Delta }$E${}_R$ of 0.7 kcal/mol) [110]; see Figure 5. A barrier of this height prevents the non-catalyzed reaction from being useful in synthesis. In contrast, the catalyzed reaction proceeds in two steps via a Au(I)-acetonium intermediat (DLPNO-CCSD(T): step 1, $\mathsf{\Delta }$E${}^{\#}$ of 18.2 kcal/mol and $\mathsf{\Delta }$E${}_R$ of 12.5 kcal/mol; step 2, $\mathsf{\Delta }$E${}^{\#}$ of 2.3 kcal/mol and $\mathsf{\Delta }$E${}_R$ of -9.2 kcal/mol) [110].

In the following we will focus on one important aspect of this study; namely, how the analysis of the reaction path curvature leads to a detailed understanding of the mechanistic differences between the gas phase reaction and its catalyzed counterpart. Reaction movies following the evolution of the RC along the reaction path are provided for both the non-catalyzed and the two steps of the catalyzed reaction in the Supplementary Materials, Videos S6–S8. More details about environmental and substituent effects can be found in [110].

Figure 6a shows the curvature diagram of the non-catalyzed reaction. The pre-chemical phases 1 and 2 are characterized by the rotation of methyl groups and the acetate unit and phase 3 by a pseudo-rotation of the 6-membered ring reducing the distance between the C${}_b$ and O${}_b$ atoms, which are supposed to form a new chemical bond (see Video S6, Supplementary Materials). C${}_b$O${}_b$ bond formation (red line) starts already in the pre-chemical phases 2–4. The essential chemical event in phase 5 is the cleavage of the C${}_a$O${}_a$ bond dominating curvature peak K3 (blue line). Phase 6, containing the TS, is characterized by a smaller curvature peak K4. It is interesting to note that the RC adapts a boat form at the TS with C${}_a$O${}_a$ and C${}_b$O${}_b$ distances close to 2.1 Å. Both the C${}_a$O${}_a$ and the C${}_b$O${}_b$ bonds resist a further change in this phase (see Video S6, Supplementary Materials), indicating the possibility to break up the reaction in this area into two steps with lower activation energies, e.g., with a catalyst that transforms the TS of the non-catalyzed reaction into a stable intermediate, as depicted in Figure 5. In phase 7, leading to the curvature peak K5, the formation of the new C${}_b$O${}_b$ bond is finalized. Phases 8–10 are characterized by six-membered ring pseudo-rotations and a rotation of acetate methyl-group. In summary, the most important curvature peaks comprise: K3 before the TS (breakage of the C${}_a$O${}_a$ bond), K4 close to the TS (boat form of the RC with equal C${}_a$O${}_a$ and C${}_b$O${}_b$ bond lengths), and K5 after the TS (finalization of C${}_b$O${}_b$ bond formation).

Figure 6b,c shows the corresponding curvature diagrams of step 1 and step 2 of the Au(I)-NHC catalyzed reaction, revealing the major difference between the non-catalyzed and catalyzed reaction (see also Videos S7 and S8, Supplementary Materials). In step 1 (Figure 6b), first, the new C${}_b$O${}_b$ bond is formed with the C${}_a$O${}_a$ bond to be broken still being intact. During this process, the $\pi$-bonded gold complex transforms into an intermediate $\sigma$-bonded Au(I) complex, strongly resembling the TS of the non-catalyzed reaction. The curvature peaks K1 and K2 before the TS are dominated by the loosening of the Au(I)-C${}_b$ interaction (green line), which is less energy demanding than breaking a CO bond. The formation of the Au(I))-$\sigma$-complex is endothermic by 12.5 kcal/mol (DLPNO-CCSD(T)) [110]; therefore, in step 2 (see Figure 6c) the regeneration of a stable Au(I)-$\pi$-complex takes place, which includes the almost effortless breakage of the C${}_a$O${}_a$ bond.

In essence, the URVA analysis clearly shows how the Au(I) catalyst breaks up the non-catalyzed rearrangement into two steps by easily switching between Au(I)-$\pi$ and Au(I)-$\sigma$ complexation. These findings form the basis for a more general route map: (i) study the curvature diagram of the non-catalyzed reaction; (ii) identify energy consuming events before the TS and hidden intermediates, which could be the starting point for splitting the reaction into several steps and (iii) finding a catalyst which can (a) transform the hidden intermediates into real intermediates, and (b) can move energy-consuming events (i.e., bond breakage) into the exit channel.

4.4. Bacillus Subtilis Chorismate Mutase Catalyzed Claisen Rearrangement

In recent years, one has put a lot of emphasis on substituting traditional chemical processes relying on heavy metal catalysts (often leading to toxic byproducts) with nature’s powerful way of producing environmentally pure compounds under eco-friendly conditions using enzymes as catalysts. This new trend goes hand–in–hand with the development of new technologies in enzyme engineering [179,180,181,182,183]. If chemists would be able to use catalysts as effectively as nature uses enzymes, most of the sustainability problems in connection with chemistry could be solved. If we succeeded in mimicking aspects of natural processes while implementing features of traditional catalysts, we could take the precision and selectivity of catalysis to an entirely new level. One important prerequisite for this ambitious goal is the detailed understanding of the mechanism of enzyme catalysis at the atomic level, which is still in its infancy and often refers to more approximate methods and/or descriptions [184,185,186,187,188,189,190,191,192].

Our first successful URVA QM/MM enzyme study was the Bacillus subtilis chorismate mutase (BsCM) catalysed Claisen rearrangement of chorismate to prephenate [118]. This rearrangement forms an important part of the shikimate pathway controlling the synthesis of aromatic amino acids in the cells of fungi, bacteria and plants [193,194,195]. It is also interesting from a mechanistic point of view; being one of the few pericyclic processes taking place in a biological system, it provides a rare opportunity for studying how nature promotes such rearrangements so successfully.

Although it has been the object of numerous experimental and theoretical studies (for a comprehensive overview, see [118]), the actual enzyme mechanism lowering the reaction barrier has been still under debate suggesting the following alternatives: (i) the enzyme changes the rearrangement mechanism, (ii) space confinement of the chorismate in the active site plays a key role, or (iii) the enzyme lowers the barrier by stabilizing the TS via specific H-bonding, in particular between the migrating ether oxygen and Arg90. In order to shed light into this debate, we compared the mechanism in the gas phase (reaction R1), in solution (reactions R2 and R3) and in the enzyme (reactions R4–R6), as depicted in Figure 7. Reaction movies for reactions R1–R5 showing the evolution of the RC along the reaction path are provided in the Supplementary Materials, Videos S9–S14. In Figure 8 the QM/MM-optimized geometry of the BsCM-chorismate complex is shown.

Figure 9 shows the energy profiles for reactions R1–R6. The activation energy of the chorismate rearrangement decreases successively from the gas phase to the solution for the reaction in the BsCM; $\mathsf{\Delta }{E}^{‡}\left(\mathbf{R}\mathbf{1}\right)$ = 24.4 kcal/mol; $\mathsf{\Delta }{E}^{‡}\left(\mathbf{R}\mathbf{2}\right)$ = 19.7 kcal/mol; $\mathsf{\Delta }{E}^{‡}\left(\mathbf{R}\mathbf{3}\right)$ = 17.3 kcal/mol; $\mathsf{\Delta }{E}^{‡}\left(\mathbf{R}\mathbf{4}\right)$ = 11.6 kcal/mol; $\mathsf{\Delta }{E}^{‡}\left(\mathbf{R}\mathbf{5}\right)$ = 10.4 kcal/mol. Additionally, as expected, modification of the enzyme increases the activation energy, $\mathsf{\Delta }{E}^{‡}\left(\mathbf{R}\mathbf{6}\right)$ = 14.9 kcal/mol, (values from [118]). While these findings confirm the biological activity of BsCM, they do not provide any mechanistic details, because the energy is a cumulative property, showing the overall effects. Mechanistic details can be disclosed by the investigation of the corresponding curvature diagrams.

In Figure 10 the curvature diagrams of gas phase reaction R1 and enzyme reaction R5 are compared. This comparison clearly reveals that the reaction mechanism is the same for both media. There are two pronounced curvature peaks K1 in phase 3 and K3 in phase 5, separated by a smaller curvature enhancement K2 located in phase 4 containing the TS. The actual chemical process of C6O8 bond cleavage and C1C10 formation follows a similar pattern also for reactions R2–R4 and R6 [118] (not shown here), in line with the suggestion of several authors [196,197] that the mechanism of the chorismate rearrangement is the same in gas phase, solution, and in the enzyme. K1 located in phase 3 is the most pronounced peak before the TS. This peak is dominated by the C6O8 bond to be broken (red line) accompanied by a small resisting contribution of C1C10 (blue line), the new bond to be formed. Phase 4 is characterized by negative contributions from C6O8 and C1C10; i.e., a resistance to further change. The finalization of C1C10 bond formation and C6O8 bond cleavage occur after the TS in phase 5, characterized by curvature peak K3. This confirms the suggestion of several authors [198,199,200], that the actual bond breaking and bond forming processes are asynchronous. Noteworthy is that again the TS does not play any significant mechanistic role, as is reflected by the small curvature enhancement K2. In contrast to the Au(I)-assisted [3,3]-sigmatropic rearrangement of allyl acetate changing the bond forming/breaking process, the mechanistic differences between R1 and R5 leading to the different activation energies are related to the pre-chemical phases. While R1 is characterized by long pre-chemical phases 1 and 2 with curvature peak K4, in R5 the chemical phase 3 starts right away; i.e., while docked into the enzyme, the enolpyruvyl side chain of the chorismate is already in an optimal position to react with the cyclohexadienyl part. This clearly proves that space confinement is the key feature of the BsCM catalysis.

It has been suggested that intramolecular charge transfer from the enolpyruvyl part to the cylohexadienyl unit helps accelerating the rearrangement by decreasing the repulsion between these two negatively charged groups, which have to approach each other to form prephenate [196]. Figure 11a,b shows the total NBO charge of the enolpyruvyl group (in red color) and that of the 4-hydroxy-1,5-cyclohexdiene-1-carboxylate group (in blue color) as a function of the reaction parameter s, for the gas phase reaction R1 and the enzyme reaction R5, respectively. In the gas phase, charge transfer from enolpyruvyl to cyclohexadienly starts at K1, goes through a minimum at M4 where both parts are equally charged. After M4, a small charge transfer occurs in the opposite direction; see Figure 11a. In the enzyme reaction R5 the actual charge transfer from enolpyruvyl to cyclohexadienly starts after the TS, see Figure 11b. This clearly shows that charge transfer from the enolpyruvyl to the cyclohexdienyl group cannot contribute to the lowering of the energy barrier.

Several authors emphasized that BsCM could lower the barrier by selective H-bonding with the TS [201,202,203]. Most of these studies assume that the strength of a H-bond can be assessed from the H-bond distance. However, the shorter bond is not always the stronger bond [204,205,206]. Our local mode analysis has led to a unique and quantitative measure of bond strength based on vibrational spectroscopy [118,205,206,207,208,209,210,211,212,213,214,215,216,217,218,219,220,221,222,223,224,225,226,227,228,229,230,231,232] utilizing local mode force constants and related bond strength orders (BSO)s, which can be derived from local mode force constants via a generalized Badger rule [101,233] and two reference compounds with known bond strength orders. For H-bonds usually FH, n = 1 and F${}_2$H${}^{-}$, n = 0.5 are taken [213]. Applying the local mode analysis to the reactant, product, and TS of reaction R5 we could quantify that the enzyme does not lower the barrier via specific H-bonding in the TS, refuting these suggestions.

In Figure 12 the intermolecular H-bonds between the BsCM active site pocket and the reactant, TS, and product of reaction R5 are shown. According to Figure 12 there are 5 H-bonds (labelled with blue numbers) between the enzyme and the substrate which do not change during the reaction, showing that the enzyme confines the chorismate in an optimal position for the rearrangement already in the entrance channel. This complies with our finding that in the enzyme there is no longer a pre-chemical phase. The carboxylate group of the enolpyruvyl part of the reactant forms two strong H-bonds with Arg7 (BSO values of 0.332 and 0.333, H-bonds 1 and 2, respectively), one H-bond with Arg90 (BSO value of 0.329, H-bond 3), and one H-bond with Glu78 (BSO value of 0.278, H-bond 4), with comparable values for TS and product, and also comparable with the H-bonds of some reference systems shown in Figure 12. For H-bond 5 between the ether oxygen O8 and Arg90, which has been declared in several studies as the major source for the enzyme-TS stabilization [201,202,203,234] we find only a slight increase of the BSO for the TS, (BSO value of 0.304 compared with 0.290 for chorismate and 0.301 for prephenate).

In summary, the analysis of the reaction path curvature indisputably confirms that the actual mechanism of the chorismate rearrangement is the same in the gas phase, solution and enzyme. The process of CO bond cleavage starts before the TS and the new CC bond formation is finalized after the TS. There are significant differences in the pre-chemical phases, which are a result of the different environments. The pre-chemical phases become shorter in aqueous solution and disappear for the enzyme reaction; i.e., the chemical process of CO bond cleavage starts directly in the entrance channel. These findings are supported by the local mode analysis revealing that the inter-molecular H-bond network between BsCM and chorismate keeps the chorismate in an optimal position for the rearrangement process.

5. Conclusions

URVA is a powerful quantum chemical tool for a comprehensive and detailed study of chemical reaction mechanisms. URVA records all chemical events as the RC follows the energy valley from the entrance channel through the TS into the exit channel leading to the final products. All chemical events can be identified by analyzing the curvature of the reaction path, as shown in this review. This makes URVA feasible also for large molecular systems with many degrees of freedom and for complex reaction mechanisms, as found in catalysis. The unique curvature profile, i.e., the fingerprint of the chemical reaction, helps finding in the most efficient way, the mechanistic needles in the haystack.

The four reactions presented in this work provide a flavor of how URVA discloses the different facets of catalysis, and how based on the URVA analysis new guidelines for catalyst fine-tuning and rational catalyst design can emerge.

• The Rh catalyzed methanol carbonylation is an example for a coordination-sphere-driven catalysis catalyzing a chemical reaction by changing the coordination sphere of the transition metal to facilitate bond forming/bond breaking processes. The Rh coordination number changes from 4 to 6 in step 1, from 6 to 5 in step 2 from 5 to 6 in step 3 and back to 4 in step 4. Step 1 with an activation energy of 45.0 kcal/mol is the cause of the harsh reaction conditions. URVA identifies the approach of the reactants and the cleavage of the C–I bond as the most energy demanding processes of this step. These findings provide valuable information for catalyst modification aiming at milder reaction conditions. The long approach phase can be shortened by chelating the reactant and C–I bond breakage can be supported via polarization of the C–I bond, which will lower the overall activation energy. Work is in progress along these lines.
• In the Sharpless epoxidation of allylic alcohols the dimeric Ti catalyst mimics a surface typical of heterogenous catalysis, thus facilitating a stereospecific collision of the reaction partners, where the O–O bond of the oxidizing peroxide glides over the Ti atom. During this process the metal atom polarizes the peroxide oxygen atoms facilitating O–O bond breakage. Another important feature of the Sharpless reaction revealed by URVA is that both new C–O epoxide bonds are synchronously finalized after the TS, i.e., without energy consumption. Therefore, catalyst optimization should predominantly focus on further support of O–O breakage.
• The URVA analysis of Au(I) assisted [3,3]-sigmatropic rearrangement of allyl acetate shows how the Au(I) catalyst breaks up the non-catalyzed rearrangement into two energy saving steps by switching between Au(I)-$\pi$ and Au(I)-$\sigma$ complexation. The unfavorably high activation energy of non-catalyzed reaction is caused the fact that the migrating CO bond is broken before the TS. In contrast, the $\pi$-acidic cationic Au(I) catalyst forms a Au(I)-$\pi$-complex via the ethylene unit in the first step, supports the formation of the new CO bond while conserving the CO bond to be broken, and transforms at the end of this step into an intermediate Au(I)-$\sigma$-complex resembling the TS of the non-catalyzed reaction. In the second step, the Au(I)-$\sigma$-complex transforms back into a more stable Au(I)-$\pi$-complex including the energy conserving breakage of the migrating CO bond.
• The Bacillus subtilis chorismate mutase catalyzed Claisen rearrangement is an example for space-confinement-driven catalysis, perfectly designed by nature. URVA indisputably proves that the actual mechanism of the chorismate rearrangement is the same in the gas phase, solution and in the enzyme. The process of CO bond cleavage starts before the TS and the new CC bond formation is finalized after the TS. There are subtitle differences in the pre-chemical phases, which are a result of the different environments. The pre-chemical phases become shorter in aqueous solution and disappear for the reaction in the enzyme, where the chemical process of CO bond cleavage starts directly in the entrance channel. The local mode analysis reveals that the inter-molecular H-bond network between chorismate and BsCM does not change during the whole rearrangement, which eliminates suggestions that the enzyme lowers the barrier by stabilizing the TS through specific H-bonding.

The comparison of non-catalyzed and catalyzed reactions has led to a first roadmap for catalyst design comprising: (i) Analysis of the curvature diagram of the non-catalyzed reaction. (ii) Identification of the energy consuming events before the TS and of potential hidden intermediates as key locations for splitting the reaction into several energy conserving steps. (iii) Identification of a catalyst which can transform the hidden intermediates into real intermediates and/or can move energy consuming events into the exit channel. We hope that this review will attract the attention of the computational community and that in the future URVA will be applied as a mechanistic tool to a broad range of chemical reactions, in particular catalytic reactions.