From free-energy profiles to activation free energies

Computer simulations of chemical systems are valuable for the explanation of their experimental counterparts. In the case of chemical reactions, the quantities of primary interest are equilibrium constants and reaction rate constants or quantities directly related to these, i.e., the reaction free energy ΔF_RP (difference between free energies of products and reactants) and the activation free energy $ΔFRP‡$ (the difference between free energies of transition state and reactants). Indeed, the computation of such free energy differences has a long history.^1–7 

The kinetics of a chemical reaction can be modeled as a transition from a reactant well (R) on the potential energy surface (PES) to a product well (P). The two local minima are separated by a potential energy barrier that must be overcome as the atomic configuration changes and the reaction progresses. The total configuration space is partitioned into (hyper) volumes corresponding to R and P by a dividing (hyper) surface, the separatrix. The atomic rearrangement is described by a collective variable (CV) (or reaction coordinate), which is a function of some subset of Cartesian coordinates that gives the degree of reaction progress (e.g., 0 at R and 1 at P). In order to describe a reaction well, one needs to choose a “good” CV, i.e., one that distinguishes properly between configurations of R and P. The CV is chosen so that it has two nonoverlapping domains that correspond to the domains of R and P. It is practically impossible to find the optimal CV for a complex realistic system.⁸ One must, therefore, base the choice of CV either on chemical intuition or on recently developed machine learning-based methods.^9–13 

The free-energy profile (FEP)¹⁴ (also referred to as the potential of mean force) is defined, up to a scaling constant, as the logarithm of the marginal Boltzmann distribution of the CV (Fig. 1). The FEP is determined in practice by molecular dynamics (MD) or Monte Carlo simulations. Because R and P are often separated by high potential energy barriers that are not overcome on simulation timescales, special simulation techniques, such as importance-sampling algorithms, must often be employed to sample configuration space properly.^15–22 These algorithms usually directly yield the FEP.

Contrary to what its name implies, the FEP is not a true Helmholtz or Gibbs free energy.²³ Treating the FEP as if it were a free-energy analog of the minimum energy path is pervasive in the field and rarely acknowledged explicitly as the approximation that it is. Differences in the FEP between local extrema are then misinterpreted as reaction and activation free energies (see red highlight in Fig. 1). We have recently shown that this misconception leads to significant errors in reaction free energies, ΔF_RP.²³ The choice of the CV has a large influence on the FEP. In fact, the FEP has no meaning independent of the CV^23–25 and the structure of the FEP (e.g., the breadth and depth of local extrema or even their existence) depends on the CV. Thus, a treatment that relies solely on the shape of the FEP yields CV-dependent activation free energies. Moreover, kinetic quantities (e.g., $ΔFRP‡$⁠) derived from the FEP, which depends solely on the PES and does not account for particle masses, must be approximations. The rigorous formula for $ΔFRP‡$ derived here (see green highlight of Fig. 1 and Sec. II E) is independent of the precise mathematical form of the CV, as long as it discriminates between R and P. We show below that a poor choice of CV has an even bigger impact on $ΔFRP‡$ than on ΔF_RP.

The remainder of this article is organized as follows: In Sec. II, we first derive an expression for the rate constant k_R→P. Then, using the Eyring equation, we derive the connection between $ΔFRP‡$ and k_R→P. The physical interpretation of the components that constitute the correct activation free energy is discussed. In Sec. III, we employ simple analytic models to assess the error incurred by the common practice of taking $ΔFRP‡$ to be the difference between the values of the FEP at the maximum (transition state) and the minimum at R. Section IV is devoted to an analysis of the sensitivity of ΔF_RP and $ΔFRP‡$ to the choice of the CV. To emphasize the errors that can result from estimating $ΔFRP‡$ directly from the FEP, we examine in Sec. V a numerical one-dimensional model and two realistic chemical processes. Section VI consists of a summary of our findings and a discussion of open questions on the computation of the activation free energy. Our conclusions are summarized in Sec. VII.

The interconversion of R and P is represented by the chemical reaction

State α(=R,P) is defined by the region of configuration space it occupies, designated by Ω_α. Thus, we define the configuration integral associated with the state α by

Here, $x=x1,x2,…,x3NT$ denotes the column vector of Cartesian coordinates that specify the atomic configuration; $dx=∏i=13Ndxi$ is the 3N-dimensional volume element; U(x) is the potential energy surface (PES), and β ≡ 1/k_BT. Only those configurations x that belong to Ω_α contribute to Z_α, which is the effective volume of configuration space occupied by state α. We assume that Ω_R and Ω_P constitute the whole configuration space available to the system and they are separated by a (3N − 1)-dimensional dividing (hyper) surface, normally taken to contain the ridge of the barrier of the PES between the minima corresponding to R and P.

The course of the reaction can be monitored by a (scalar) CV (or reaction coordinate), ξ(x), which is a function of a subset of the atomic coordinates that gives a measure of the progress of the reaction. The CV is chosen such that Ω_R and Ω_P correspond to nonoverlapping domains of the CV. Ideally, the gradient of ξ(x) should be normal to the dividing surface, on which the CV assumes a particular value z_TS. In this case, the CV discriminates properly between R and P.

It is convenient to introduce mass-weighted coordinates

where M stands for the 3N × 3N diagonal matrix of atomic masses. In terms of mass-weighted coordinates, the Hamiltonian is

where $p̃i=x̃̇i$ is the momentum conjugate to the coordinate $x̃i$⁠. Henceforth, we employ the condensed notation of the second line of Eq. (4), where $p̃$ stands for the column vector of momenta.

The treatment of the reaction rate is facilitated by employment of a special set of coordinates, one of which is the CV. Hence, we transform from mass-weighted coordinates to a complete set of curvilinear coordinates, $q=q(x̃)$⁠, of which we take $q1(x̃)=ξ(x̃)$⁠. From the inverse transformation $x̃=x̃(q)$⁠, we obtain

where $[J]ij=∂x̃i∂qj$ is an element of the Jacobian. The momentum conjugate to q is

where

the mass matrix in curvilinear coordinates, is also referred to as the mass-metric tensor (see, for example, Refs. 2 and 26–28). In general, M_q is a full matrix. The Hamiltonian is given in curvilinear coordinates by

From Eq. (7), we deduce the following expression for the effective inverse mass matrix:

where we employ $[J−1]ik=∂qi∂x̃k$ and $(∇x̃qi)T=(∂qi/∂x̃1,∂qi/∂x̃2,…,∂qi/∂x̃3N)$ is the 3N-dimensional mass-weighted gradient. Using Eq. (3), we get the following from Eq. (9):

Note the distinction between ∇_x for the Cartesian gradient and $∇x̃$ for the gradient with respect to mass-weighted coordinates.

We assume the system to be in thermodynamic equilibrium. Then, the rate of the forward reaction equals the rate of the backward reaction,

where k_R→P and k_P→R are the forward and backward rate constants, respectively, and $P(R)$ and $P(P)$ are the respective probabilities of observing R and P. The rate can also be expressed in terms of the frequency ν of crossing the dividing surface in either the forward or backward direction [i.e., of the number of times per unit time that $ξ(x̃)−zTS$ changes sign]. Since the forward and backward rates are equal, either rate must equal ν/2. Thus, focusing on the forward rate, we have from Eq. (11) that

The following alternative expression for the rate constant is frequently used:^29–34 

Here $p,q$ denotes the ensemble average over all of phase space, δ the Dirac delta function, Θ the Heaviside function, and $ξ̇$ the time derivative of the CV. The equivalency of the two expressions is proven in the supplementary material.

The frequency of crossing the dividing surface can be expressed formally as the time average of the frequency with which $ξ(x̃)−zTS$ changes sign,³⁵ 

A proof of this expression is provided in the supplementary material. Assuming the system to be ergodic, we can recast the time average as an ensemble average,

where $H$ is given by Eq. (4). We next transform from mass-weighted to curvilinear coordinates. From Eqs. (5)–(7), we get

where the second equality invokes the definition of the inverse Jacobian. Substitution of Eq. (16) into Eq. (15) and transformation to curvilinear coordinates yield

To simplify this expression, we exploit the freedom afforded by curvilinear coordinates. While the “first” is chosen to be the CV, the remaining 3N − 1 are as yet unspecified. Hence, we require that q₂, q₃, …, q_3N be orthogonal to q₁ = ξ, which constraint is expressed by

In general, the construction of the orthogonal set can be achieved in a variety of ways.¹⁶ 

Invoking Eq. (18), we can express the kinetic energy as

where in analogy to Eq. (9) we define the (3N − 1) × (3N − 1) inverse mass matrix M′⁻¹ and the (3N − 1)-dimensional momentum vector $p′=(p2,p3,…,p3N)T$⁠. Likewise, we can simplify Eq. (16) as follows:

Plugging Eqs. (19) and (20) into Eq. (17), we get

Performing the integration on p₁ gives

Inserting the identity $1=|∇x̃ξ|(2πkBT)−1/2∫−∞∞dp1e−β|∇x̃ξ|2p122$ into Eq. (22) at the place indicated, we obtain

Transforming back to Cartesian coordinates yields

where $$ indicates the ensemble average over configuration space. Using the fact that

is the normalized probability density of observing an atomic configuration x such that ξ(x) = z, we can recast Eq. (24) as

where $zTS$ signifies an average over the dividing surface. The second line of Eq. (26) follows from Eq. (3); the third line implicitly defines m_ξ, which we interpret as the effective mass of the pseudo-particle associated with the coordinate ξ(x),

which is the 1,1 element of the inverse mass-metric tensor [see Eq. (10)].^2,26–28 Finally, combining Eqs. (12) and (26), we obtain

Eyring’s equation relates the rate constant to a free energy of activation by defining a modified equilibrium constant for the formation of activated complex from reactant R (see, for example, Ref. 36). In the present notation, the equation is

where h is the Planck constant. We use the symbol F for the Helmholtz free energy in order to distinguish it from the free-energy profile denoted by A [see Eq. (32)]. We solve Eq. (29) for the activation free energy and combine the result with Eq. (28) to get

where $λξ≡h2/2πmξkBT$⁠. We interpret λ_ξ as the de Broglie thermal wavelength of the pseudo-particle associated with the CV.

By expanding the logarithm in Eq. (30), we can recast the “exact” expression for the activation free energy as

The second line of Eq. (31) depends on the definition of the free-energy profile (FEP),^16,23

and on the relation²³ 

A frequently employed procedure is to set the activation free energy equal to the difference between the maximum of the FEP at z_TS and the minimum at z_R,min,

We place a tilde on this formula to distinguish it from the “exact” one in Eq. (30). Thus, $ΔF̃RP‡$ can be viewed as an approximation. For example, if the density is strongly peaked about z_R,min, then $kBT⁡lnP(R)≈−A(zR,min)$⁠, according to Eqs. (32) and (33). Under this condition, the approximate formula agrees with the exact, except for the term $−kBT⁡lnλξzTS$⁠. Therefore, the influence of distortions of the coordinate system induced by ξ(x) is ignored by $ΔF̃RP‡$⁠, as is the influence of mass [see Eq. (27)].

An alternative recasting of the exact formula for the activation free energy, Eq. (30), is instructive. Invoking the relations²³ 

and

where q_R is the molecular partition function of R and $Λ≡∏i=13Nh2/2πmikBT$ (the product of all Cartesian de Broglie wavelengths), we rewrite the exact expression as

The second term on the right-hand side of Eq. (37) is the negative of the free energy of R.²³ Likewise, if we regard $q‡≡Zρ(zTS)λξzTSΛ$ as the effective partition function with z fixed at z_TS, then the first term is the free energy of the constrained system. That q^‡ has the stated character can be demonstrated explicitly in the case the curvilinear coordinates form a complete orthogonal set. Then, we can rewrite Eq. (37) as

This form of $ΔFRP‡$ is very intuitive: The activation free energy is the difference between the free energy of the system constrained to the dividing surface, F^‡, and the free energy of the reactant, F_R. Moreover, it is noteworthy that Eq. (38) assumes the same form as the corresponding expression derived by conventional transition state theory.³⁶ 

In order to gauge the error incurred by approximating the activation free energy $ΔF̃RP‡$ [Eq. (34)] in comparison to the “exact” $ΔFRP‡$ [Eq. (30)], we study the behavior of two analytically treatable models. Each consists of a single particle of mass m moving in one dimension. The PESs are meant to represent a system with two minima, which are approximated either by square wells (SW) or by parabolic (harmonic oscillator) wells (HO). Their detailed treatment is presented in the supplementary material. We take the difference between approximate and “exact” activation free energy as a correction term, which we derive to be

In Eq. (39), L_R denotes the width of the reactant square well. In Eq. (40), k is the force constant of the harmonic well.

We note that $ΔF̃‡$ does not depend on particle mass (m) (as it is only derived from a marginal Boltzmann distribution), and in the one-dimensional case, it depends neither on temperature (T) nor on the parameters of the PES (L_R and k). Thus, we regard the difference as a correction of $ΔF̃‡$ that accounts for the influence of these parameters. Though the corrections for the two models exhibit different dependencies on the parameters, they can nevertheless be correlated. We note directly, for example, that both corrections increase at the same rate with increasing m. Further, both increase with increasing T, although corr_HO increases more rapidly. Concerning the PES parameters, we observe that corr_SW increases with increasing L_R, whereas corr_HO increases with increasing k⁻¹. This is expected, since as k decreases, the harmonic potential broadens, allowing the particle to move in an effectively larger domain of R, just as an increase in L_R does.

The one-dimensional HO model can be roughly correlated with realistic multidimensional systems. We observe that $νR=k/m/2π$ is the frequency of oscillation of the particle about the minimum x_R,min. Hence, we can recast the correction given by Eq. (40) as

For reactions carried out around room temperature T_◦ = 300 K, a reference frequency ν_◦ = k_BT_◦/h ≈ 6.0 × 10¹² s⁻¹ can be defined. Thus, for molecular vibrations around this frequency, the correction is negligible. In the typical case, where the masses of constituent atoms (e.g., H, C, and O) are small, and the bonds are stiff, ν_R → ν_◦ and the correction is small. On the other hand, for reactions involving more massive atoms and “soft” degrees of freedom, ν_R < ν_◦ and we expect substantial corrections.

The validity of the formulas describing the activation free energy [Eq. (30)], and the reaction free energy,²³ 

depends on the assumption that the CV distinguishes properly between R and P, as defined by the dividing surface $S$⁠. Thus, knowledge of $S$ is crucial to the proper choice of CV. For low-dimensional model systems, the choice is generally clear, but for realistic multidimensional systems, one usually has little or no information about $S$ and must base their choice on heuristics and chemical intuition. Such intuitive CVs can lead to significant errors.

In this section, we systematically explore the influence of the choice of the CV on ΔF_RP and $ΔFRP‡$⁠. For this purpose, we employ the following model: a single particle of mass m moving in two dimensions on the PES given by

a contour plot of which is shown in Fig. 2. The particle coordinates x and y are given in units of Å and the energy in units of kJ/mol. The parameters ϵ, b, and c are taken to be 25 kJ mol⁻¹ Å⁻⁴, 2 Ų, and 0.25 ų, respectively. The parameter ϵ effectively controls the height of the barrier of the PES between R and P; c controls the difference between the minima of R and P. The values are chosen to yield realistic free energies (ΔF_RP = −12.28 kJ/mol and $ΔFRP‡=16.06$ kJ/mol, which is roughly the activation free energy of the internal rotation of butane³⁷). The ideal CV is ξ(x, y) = x and the dividing surface coincides with the line x = x_max = −0.067 25 (see Fig. 2). Clearly, ∇ξ · ∇U vanishes on $S$⁠, which is the constraint that should be obeyed by a CV that properly discriminates between R and P.³⁵ 

To vary the choice of the CV systematically, we define the CV by

where a is restricted to the interval [0, 1]. We determine the value of a by specifying the angle θ between ∇ξ and $eS$⁠, the unit vector parallel with the true separatrix $S$ (i.e., e_y). In other words, a and, therefore, ξ, are determined by the condition $∇ξ|∇ξ|⋅eS=cos⁡θ$⁠. (Details of the calculation are provided in the supplementary material.) Corresponding to a given θ (i.e., a given choice of the CV) is a “trial” separatrix S(θ), which is a line having the equation y = −a(x − x_max)/(1 − a), where x_max is the x-coordinate of the saddle point on the PES. When a = 1, θ = 90°. In this limit, S(90°) coincides with $S$⁠. As a decreases from 1 to 0, S(θ) rotates counterclockwise about the point (x_max, 0). The trial separatrix S(45°) is shown in Fig. 2. In the limit a = 0, ∇ξ = e_y, θ = 0. Hence, S(0°) is normal to $S$⁠, which makes ξ(x, y) = y the worst possible choice of the CV.

For a given θ, we calculate the probability density ρ(z) using Eq. (25). As shown in the supplementary material, the evaluation of the required double integrals is facilitated by transforming from Cartesian to orthogonal coordinates q₁ = ξ(x, y) and q₂ = (a − 1)x + ay. We obtain the FEP using Eq. (32). Illustrative plots of ρ(z) and A(z) are shown in Figs. 3(a)–3(c) for three CV choices. The local maximum of the FEP, z_max, defines the domains of R and P. We note, however, that the FEPs for θ < 32° lack any such local maximum. We henceforth ignore these choices, as the CV cannot distinguish R from P at all.

As a measure of the quality of the chosen CV, we adopt a modification of the procedure introduced previously,²³ which was to monitor the quantity $D(z)=∇ξ(x)⋅∇U(x)z$⁠. We note that D(z_TS) is exactly zero on $S$ for the ideal CV (i.e., the one that discriminates perfectly between R and P). However, away from $S$⁠, or in case the choice of CV is not ideal, D(z) is difficult to interpret, because it depends so strongly on the local gradient of the PES. To ameliorate this defect, we propose a scaled, dimensionless orthogonality measure defined by

where we replace the gradients of U and ξ with their corresponding unit vectors. Thus, D_s(z_TS) is zero on $S$ for the ideal CV, where the gradients of U and ξ are perpendicular, and unity where they are parallel.

One can see in Fig. 3(d) that for the ideal CV ξ(x, y) = x, D and D_s have very sharp roots at z_TS, indicating that the CV is orthogonal to the separatrix. Because of the symmetry of the PES, the two measures have two additional roots located at the minima of reactant and product. D_s does not actually reach zero on account of the finite numerical resolution of our computation. However, the sharp minima are still visible. Figures 3(e) and 3(f) show the orthogonality measure for nonideal CVs. The shape of the D-measures changes drastically. Most significantly, the sharp root or minimum at the maximum of the FEP turns into a local maximum for both D and D_s, which is an unmistakable sign that results for these CVs cannot be trusted (see the dependence of the $ΔFRP‡$ on θ in Fig. 4).

Using the numerically computed ρ(z), we calculate the reaction free energy and activation free energy, which are given, respectively, by Eqs. (42) and (30), where we set z_TS = z_max. In Fig. 4, we plot ΔF_RP, $ΔFRP‡$⁠, D(z_max), and D_s(z_max) as functions of θ. Figure 4(d) shows clearly how sensitive D_s(z_max) is to the choice of CV. At θ = 90° D_s(z_max) vanishes, since the chosen CV coincides with the ideal one. But, as θ decreases, D_s(z_max) rises sharply over a narrow interval of about 10°. That is, for large θ, ∇U and ∇ξ are almost orthogonal, whereas with decreasing θ, they become nearly parallel. The fall off of D_s(z_max) as θ decreases from about 45° is due to the interference of force vectors that are almost isotropically distributed, and result in essentially randomized alignment of the force and CV gradient vectors.

Since ρ(z) is strongly peaked around the minima of R and P (see Fig. 3), the choices of CV in the range of 45°–90° separate the minima well. As a consequence, ΔF_RP is essentially independent of the choice in this range [see Fig. 4(a)]. In other words, over this range of choices, one obtains an accurate value of the reaction free energy. Only for θ < 45°, where the CV begins to fail to discriminate between R and P, does the error in ΔF_RP set in rapidly.

As seen in Fig. 4(b), the activation free energy is dramatically more sensitive than ΔF_RP to the choice of CV. It deviates from the correct value by more than “chemical accuracy” (1 kcal/mol) at θ ≈ 60°. For θ < 40°, $ΔFRP‡$ even becomes negative. If this were correct, the rate of reaction would decrease with increasing temperature. This apparent sensitivity can be reasoned as follows: All points on the true separatrix have very low likelihood. A trial separatrix with θ < 90° includes more likely configurations and, therefore, overestimates ρ(z_TS). Since the true ρ(z_TS) is very small, the relative error is large. For large probabilities, e.g., $P(R)$⁠, the same absolute error would incur a much smaller relative error. The relative error in the density directly translates to an absolute error in the activation free energy because of the logarithm of ρ(z_TS) [see Eq. (31)].

The fact that ΔF_RP is largely unaffected by the choice of the CV explains why CVs based purely on chemical intuition can yield reaction free energies comparable with experiment. However, ΔF_RP is expected to become somewhat more sensitive to the choice of CV for more complex PES. Compared with the reaction free energy, the activation free energy is generally more sensitive. Hence, to achieve the same accuracy for $ΔFRP‡$ and ΔF_RP, one must choose the CV with a great deal of care.

To further illustrate the errors that one may incur by estimating $ΔFRP‡$ directly from the FEP alone [i.e., by invoking Eq. (34)], we consider first a simple one-dimensional model that can be treated for the most part analytically and then models of two real chemical processes.

We consider a single particle of mass m moving in one dimension on the PES given by

where ϵ, which controls the steepness of the potential barrier, has units of kJ/mol. The parameter a, which controls the width of the barrier, is set to 1 Å⁻² and b = 1 Å. The PES, plotted in Fig. 5(a) for the case ϵ = 5 kJ/mol, has two equal minima separated by a maximum at x = 0. Because U(x) diverges as x approaches −5 or 5, the particle is confined to the domain −5 < x < 5. R and P correspond, respectively, to the domains −5 < x < 0 and 0 < x < 5. The symmetry of the PES dictates that $P(R)=P(P)=0.5$⁠. Therefore, from Eq. (12), we get

where ν is the crossing frequency. From Eqs. (29) and (47), we deduce the following expression:

We compute ν by molecular dynamics (MD) simulation, as detailed in the supplementary material. MD simulations were carried out at five temperatures in the range of 100–1000 K and for five different particle masses in the range of 1–100 amu.

We consider two CVs: ξ₁(x) = x and ξ₂(x) = 1/(x + 5). Using Eq. (32), we obtain the corresponding FEPs,

Setting ϵ = 5 kJ/mol ensures that even the most massive particle considered crosses the dividing surface at the lowest temperature during the 10 ns time interval of the MD simulation. Figures 5(b) and 5(c) show plots of the FEPs based on Eqs. (49) and (50). We note the strong distortion of configuration space induced by ξ₂(x). The domains of R and P are reversed, the minima are not equal, and the maximum of the barrier between R and P does not occur precisely at z = 0.2, the inverse of the position of the maximum of the barrier of the PES at x = 0.

Approximate activation free energies obtained according to Eq. (34) are listed in Table I, along with “exact” values $ΔFRP‡$ obtained from Eq. (30), which yields exactly the same result for both CVs, and from Eq. (48) via MD. The excellent agreement between the values obtained from Eqs. (30) and (48) is gratifying. According to Eq. (49), $ΔF̃1‡$ should be independent of both temperature and particle mass. Likewise, $ΔF̃2‡$ should depend on temperature, but we note that by definition $ΔF̃2‡$ is independent of mass. Table I bears out these expectations.

The dominant impression of Table I is the severe lack of agreement between approximate and exact activation free energies. The impact of the loss of the symmetry of the PES by ξ₂ is particularly evident. Since z_TS ≈ z_max, the results in columns a and b, which correspond to the forward reaction, agree quite well, as do those of columns c and d for the backward reaction. However, the magnitudes of the forward and backward activation free energies differ greatly. Even more noteworthy is the contrary dependence of the activation free energy on temperature. For the forward reaction, it decreases with T, whereas for the backward reaction it increases markedly with T.

Examination of the exact data reveals the following general trends. At fixed T, $ΔFRP‡$ increases with particle mass m; the higher the value of T, the greater the increase. At fixed m, $ΔFRP‡$ increases with T; the greater the value of m, the greater the increase. Those are the same trends observed for the analytical models in Sec. III.

To see the influence of the parameter ϵ, we set ϵ = 50 kJ/mol. Unbiased molecular dynamics simulations were not performed for this choice of ϵ as no barrier crossings would be observed within the previously employed simulation time. Figure S2 of the supplementary material displays plots of the PES and FEPs and Table II lists approximate and exact free energies of activation. In this case, the immediate impression from Table II is the greatly improved agreement between approximate and “exact” results. Though the symmetry is still lost by ξ₂, the distortion is relatively less severe, so that forward and backward activation energies differ less. The contrary dependence of forward and reverse activation energy on T persists, but it is relatively weaker.

The trends in $ΔFRP‡$ noted above for the case ϵ = 5 kJ/mol hold for ϵ = 50 kJ/mol, but the observed variations are relatively smaller. For example, whereas the change in $ΔFRP‡$ for ϵ = 5 kJ/mol at T = 300 K is about 80% over the range of particle mass considered, it is only 13% for ϵ = 50 kJ/mol. A similar observation holds for variations of $ΔFRP‡$ with T at fixed m.

We stress that since both CVs perfectly distinguish between R and P, the computed “exact” activation free energy is identical for either, even though the CVs are very dissimilar.

We consider the realistic three-dimensional model system pictured in Fig. 6(a): a $[Cu(NH3)2]+$-complex migrating between cavities (A and B) in chabazite, a mixed crystal of the family of zeolites. This process is of importance in the deactivation of nitrogen oxides where copper-exchanged zeolites are used as catalysts.^38–42 The migration can be regarded as a “chemical reaction,” in which the Cu-complex in cavity A or B is the “reactant” or “product,” respectively. The reaction consists of the complex diffusing out of cavity A through the eight-ring (eight silicon sites) window and into cavity B. Millan et al.⁴³ have simulated this system by means of ab initio MD combined with umbrella sampling (for details, see Ref. 43). The CV they employ, which is depicted in Fig. 6(b), is defined with respect to the eight-ring window that separates the cavities. It is the projection of the vector position of the Cu atom onto the normal to the “average” plane of the central four Si and two O atoms of the ring that remain nearly in the same plane.

Our primary purpose is to analyze the data of Millan et al.⁴³ in order to determine the exact values of the reaction free energy and activation free energy for the migration reaction described above. We are especially interested in the effect of mass on the activation free energy. The authors of Ref. 43 supplied the coordinates of the trajectories and the bias used for the umbrella sampling for every frame. We implemented the CV in pyTorch⁴⁴ to gain easy access to ∇ξ and, consequently, $mξ−1$ [see Eq. (27)], through the automatic differentiation in Torch. We computed the weights of every frame with an in-house implementation of MBAR.⁴⁵ The weights were used to recompute the FEP and compare it with the result of Millan et al.,⁴³ as well as to compute the conditional ensemble average of $mξ−1$ needed for the calculation of $λξzTS$ [see Eq. (30)].

The FEPs are plotted in Fig. 7, which shows that the agreement of our FEP with that of Millan et al.⁴³ is excellent. The probability densities are normalized according to $∫−44dze−βA(z)=1$⁠. Millan et al.⁴³ take the maximum of A(z), located at z = 0.35 Å, to be the position of the TS. According to the definition of the CV, the TS should be at z = 0.0 Å. We computed exact and approximate reaction and activation free energies for both choices of the TS. Table III shows very clearly the large influence of mass on the activation free energy. Further, the approximate free energies (⁠$ΔF̃AB$ and $ΔF̃AB‡$⁠) obtained by us agree well with those of Millan et al.⁴³ The precise choice of z_TS has little effect on the activation free energies, because the FEP is quite flat around z = 0.

Since Millan et al.⁴³ used the same CV for all of the systems they simulated, the correction of the activation free energy should be about the same for all. Therefore, the correction should not affect the ordering of the barriers $(ΔF̃AB‡)$ they determined approximately. However, we would expect any comparison with experimental activation barriers to depend strongly on the difference between the approximate and exact treatments.

As a second chemical example, we consider the intramolecular cyclization of the 5-hexenyl radical (see Fig. 8), a radical clock reaction.⁴⁶ The forward reaction involves the formation of a new single bond and the conversion of a C–C double bond to a single bond. Carbon single bonds are usually stiff and have high activation barriers, as reflected in the experimental activation free energy for the cyclization, $ΔFexp‡(300K)=42±4$ kJ/mol.⁴⁷ Hence, we expect the approximate relation in Eq. (34) to hold.

As CV, we choose the distance between the two carbon atoms (C1 and C5) that form a new bond, ξ = d(C1 − C5). The associated mass m_ξ is constant and equal to the reduced mass of the two carbon atoms (i.e., 6 amu).

The system was simulated at 300 K by means of ab initio MD at the ωB97M-V/def2-TZVP^48,49 level of theory and solvated in benzene with the COSMO continuum solvation model.⁵⁰ We employed WTM-eABF^51–53 as enhanced sampling algorithm. The unbiased weights were recovered with the recently developed combination of eABF and MBAR.⁵⁴ Details of the simulation are given in the supplementary material.

The FEP [Fig. 9(a)] shows one deep minimum for P (methylcyclopentane radical) and three shallow minima for R (5-hexenyl radical). We take all configurations with z > 2.2 Å to belong to R.

The scaled orthogonality measure D_s [Eq. (45), Fig. 9(b)] is lower than 0.25 for almost the entire range of z values, rising sharply only at the ends of the simulated range. The plot of D_s shows a clear local minimum near the local maximum of the FEP, indicating that it is a good CV.

In Table IV, we can see that the exact reaction and activation free energies obtained from Eqs. (42) and (30), respectively, agree well with the approximate ones.

Hence, this example confirms that Eq. (34) does hold in cases of high barriers, low temperatures, light CVs, and narrow wells about the minima of R and P on the PES.

This study is not the first work to present expressions for the rate constant and activation free energy based on transition state theory.^29–34,55 However, previous work often lacks a stepwise derivation of their expression for the rate constant. Further, Refs. 55 and 34, which also present equations for the activation free energy, still include local differences of the FEP in their final expressions, which can thus be interpreted as corrections to the approximate treatment. Because of complex notation, it is difficult to verify whether their expressions are equivalent to our Eq. (30). It is perhaps due to the complexity of the equations and lack of physical interpretability that their expressions have not been widely adopted. Therefore, we are motivated to present a meticulous and straightforward derivation of the exact formula [Eq. (30)] for the activation free energy $ΔFRP‡$ for the two-state process from a reactant R to a product P in a novel form. The formula involves three key quantities having clear physical interpretations. Two of these, ρ(z_TS) and $P(R)=∫ΩRdzρ(z)$⁠, depend only on ρ(z), the marginal probability density that the CV ξ(x) takes the value z. The third, $λξzTS$⁠, can be rewritten as $h2/2πkBTmξ−1zTS$ to indicate explicitly the dependence on the effective mass of the pseudo-particle associated with the CV. The three clearly defined terms also facilitate implementation.

The presence of the factor $mξ−1zTS$ in the exact formula for k_R→P [Eq. (28)] shows that knowledge of ρ(z) [or alternatively A(z)] alone is insufficient to determine the rate constant k_R→P. We note that in the “conventional” transition state theory,³⁶ the rate constant is expressed in terms of canonical partition functions for reactant and activated complex [minus that associated with the CV (reaction coordinate)] and the discrete masses of the atoms enter into them. In the present treatment, the effective mass m_ξ depends not only on the discrete masses of atoms but also on the gradient of the CV [see Eq. (27)]. If the CV is linear in the Cartesian coordinates, then $mξ−1zTS$ is readily expressible explicitly in terms of the discrete masses.⁵⁶ In general, however, the CV-conditioned ensemble average must be computed.

The “gauge-independent geometric” free-energy profile, given by