From Free-Energy Profiles to Activation Free Energies

Given a chemical reaction going from reactant (R) to the product (P) on a potential energy surface (PES) and a collective variable (CV) that discriminates between R and P, one can define a free-energy profile (FEP) as the logarithm of the marginal Boltzmann distribution of the CV. The FEP is not a true free energy, however, it is common to treat the FEP as the free-energy analog of the minimum energy path on the PES and to take the activation free energy, $\Delta F^\ddagger_\mathrm{RP}$, as the difference between the maximum of the FEP at the transition state and the minimum at R. We show that this approximation can result in large errors. Since the FEP depends on the CV, it is therefore not unique, and different, discriminating CVs can yield different activation free energies for the same reaction. We derive an exact expression for the activation free energy that avoids this ambiguity with respect to the choice of CV. We find $\Delta F^\ddagger_\mathrm{RP}$ to be a combination of the probability of the system being in the reactant state, the probability density at the transition state surface, and the thermal de~Broglie wavelength associated with the transition from R to P. We then evaluate the activation free energies based on our formalism for simple analytic models and realistic chemical systems. The analytic models show that the widespread FEP-based approximation applies only at low temperatures for CVs for which the effective mass of the associated pseudo-particle is small. Most chemical reactions of practical interest involve polyatomic molecules with complex, high-dimensional PES that cannot be treated analytically and pose the added challenge of choosing a good CV, typically through heuristics. We study the influence of the choice of CV and find that, while the reaction free energy is largely unaffected, $\Delta F^\ddagger_\mathrm{RP}$ is quite sensitive.


I. INTRODUCTION
Computer simulations of chemical systems are valuable for the explanation of their experimental counterparts. In the case of chemical reactions, 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 ∆F ‡ RP (the difference between free energies of transition state and reactants). Indeed, the computation of such free energy differences has a long history. [1][2][3][4][5][6][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 a) Electronic mail: rafagb@mit.edu 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 non-overlapping domains that correspond to the domains of R and P. It is practically impossible to find the optimal CV for a complex realistic system. 8 One must therefore base the choice of CV either on chemical intuition or on recently developed machine learning-based methods. [9][10][11][12][13] The free-energy profile (FEP) 14 (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][16][17][18][19][20][21][22] These algorithms usually directly yield the FEP.
Contrary to what the name implies, the FEP is not a true Helmholtz or Gibbs free energy. 23 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 . 23 The choice of the CV has a large influence on the FEP. In fact, the FEP has no meaning independent of the CV [23][24][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., ∆F ‡ RP ) derived from the FEP, which depends solely on the PES and does not account for particle masses, must be approximations. The rigorous formula for ∆F ‡ RP 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 ∆F ‡ RP than on ∆F RP .
The remainder of the article is organized as follows. In Section II we first derive an expression for the rate constant k R→P . Then, using the Eyring equation, we derive the connection between ∆F ‡ RP and k R→P . The physical interpretation of the components that constitute the correct activation free energy is discussed. In Section III we employ simple analytic models to assess the error incurred by the common practice of taking ∆F ‡ RP 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 ∆F ‡ RP to the choice of the CV. To emphasize the errors that can result from estimating ∆F ‡ RP directly from the FEP, we examine in Section V a numerical onedimensional 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 Section VII.

A. Description of the System
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 = (x 1 , x 2 , . . . , x 3N ) T denotes the column vector of Cartesian coordinates that specify the atomic config- is the potential energy surface (PES), and β ≡ 1/k B T . 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 non-overlapping 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.

B. Curvilinear Coordinates
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 q 1 ( x) = ξ( x). From the inverse transformation x = x(q) we obtaiṅ where [J] ij = ∂ xi ∂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,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 q i ) T = (∂q i /∂ x 1 , ∂q i /∂ x 2 , . . . , ∂q i /∂ x 3N ) is the 3N -dimensional mass-weighted gradient. Using eq. (3), we get from eq. (9) Note the distinction between ∇ x for the Cartesian gradient and ∇ x for the gradient with respect to massweighted coordinates.

C. Reaction Rate Constant
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, 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)−z TS changes sign). Since the forward and backward rates are equal, then either rate must equal ν/2. Thus, focusing on the forward rate, we have from eq. (11) The following alternative expression for the rate constant is frequently used: [29][30][31][32][33][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.

D. Frequency of Crossing the Dividing Surface
The frequency of crossing the dividing surface can be expressed formally as the time average of the frequency with which ξ( x) − z TS changes sign 35 : 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), (6), and (7) we geṫ where the second equality invokes the definition of the inverse Jacobian. Substitution of eq. (16) into eq. (15) and transformation to curvilinear coordinates yields ν = dq e −βU (q) dp e − β 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 2 , q 3 , . . . , q 3N be orthogonal to q 1 = ξ, which constraint is expressed by In general, the construction of the orthogonal set can be achieved in a variety of ways. 16 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 −1 and the (3N − 1)dimensional momentum vector p = (p 2 , p 3 , . . . , p 3N ) T . Likewise, we can simplify eq. (16) Plugging eqs. (19) and (20) into eq. (17), we get Performing the integration on p 1 gives Inserting the identity 1 = |∇ x ξ|(2πk B T ) −1/2 ∞ −∞ dp 1 e − β|∇ x ξ| 2 p 2 1 2 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 re-cast 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

E. Free Energy of Activation
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 Planck's 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 λ ξ ≡ h 2 /2πm ξ k B T . 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 23 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 k B T ln P(R) ≈ −A(z R,min ), according to eqs. (32) and (33). Under this condition the approximate formula agrees with the exact, except for the term −k B T 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 23 and where q R is the molecular partition function of R and T (the product of all Cartesian de Broglie wavelengths), we rewrite the exact expression as The second term on the right side of eq. (37) is the (negative of the) free energy of R. 23 Likewise, if we regard 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 case the curvilinear coordinates form a complete orthogonal set. Then we can rewrite eq. (37) as This form of ∆F ‡ RP 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. 36

III. IMPACT OF APPROXIMATING THE ACTIVATION FREE ENERGY
In order to gauge the error incurred by approximating the activation free energy ∆ F ‡ RP (eq. (34)) in comparison to the "exact" ∆F ‡ RP (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 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 neither on temperature (T ), nor 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 −1 . 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 multi-dimensional 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 B T • /h ≈ 6.0 10 12 s −1 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.

IV. THE INFLUENCE OF THE CHOICE OF CV
The validity of the formulas describing the activation free energy (eq. (30)), and the reaction free energy 23 depend 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 ∆F ‡ RP . For this purpose we employ the following model: a single particle of mass m moving in two dimensions on the PES  yield realistic free energies (∆F RP = −12.28 kJ/mol and ∆F ‡ RP = 16.06 kJ/mol, which is roughly the activation free energy of the internal rotation of butane 37 ). The ideal CV is ξ(x, y) = x and the dividing surface coincides with the line x = x max = −0.06725 (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. 35 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 e S , the unit vector parallel with the true separatrix S (i.e., e y ). In other words, a, and therefore ξ, are determined by the condition ∇ξ |∇ξ| · e S = 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 , where x max is the x-coordinate of the saddle point on the PES. When a = 1, then θ = 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 1 = ξ(x, y) and q 2 = (a − 1)x + ay. We obtain the FEP using eq. (32). Illustrative plots of ρ(z) and A(z) are shown in Fig. 3a-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 23 , 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. 3d 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 3e and 3f show the orthogonality measure for non-ideal 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 ∆F ‡ RP on θ in Fig. 4).
Using the numerically computed ρ(z), we calculate the reaction free energy and activation free energy, which are given, respectively, by eq. (42) and eq. (30), where we set z TS = z max . In Fig. 4 we plot ∆F RP , ∆F ‡ RP , D(z max ), and D s (z max ) as functions of θ. Fig. 4d 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 • to 90 • separate the minima well. As a consequence, ∆F RP is essentially independent of the choice in this range (see Fig. 4a). 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. 4b, 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 • , ∆F ‡ RP 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 ∆F ‡ RP and ∆F RP one must choose the CV with a great deal of care.

V. PITFALLS IN THE ESTIMATION OF THE ACTIVATION FREE ENERGY FROM THE FEP
To further illustrate the errors that one may incur by estimating ∆F ‡ RP 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.

A. One-dimensional Model
We consider a single particle of mass m moving in one dimension on the PES 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Å −2 and b = 1Å. The PES, plotted in Fig. 5a 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: I. Activation free energies (kJ mol −1 ) for one-dimensional model PES U (eq. (46)) with = 5 kJ/mol for selections of temperatures (Kelvin) and particle masses (amu). ∆ F ‡ 1 = A1(zTS) − A1(zR,min) Letters above columns specify following differences: a) A2(zTS)−A2(zR,min), b) A2(zmax)−A2(zR,min), c) A2(zTS)−A2(zP,min), and d) A2(zmax)−A2(zP,min). Numbers above columns specify particle masses. 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: ξ 1 (x) = x and ξ 2 (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. Figs. 5b and 5c show plots of the FEPs based on eqs. (49) and (50). We note the strong distortion of configuration space induced by ξ 2 (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 Tab. I, along with "exact" values ∆F ‡ RP 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. Tab. I bears out these expectations.
The dominant impression of Tab. 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 ξ 2 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 TABLE II. Activation free energies (kJ mol −1 ) for onedimensional model PES U1 (eq. (46)) with = 50 kJ/mol for selections of temperatures (Kelvin) and particle masses (amu). ∆ F ‡ 1 = A1(zTS) − A1(zR,min) Letters above columns specify following differences: a) A2(zmax) − A2(zR,min) and b) A2(zmax) − A2(zP,min). Numbers above columns specify particle masses. the backward reaction it increases markedly with T .
Examination of the exact data reveals the following general trends. At fixed T , ∆F ‡ RP increases with particle mass m; the higher T , the greater the increase. At fixed m, ∆F ‡ RP increases with T ; the greater 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 Tab. II lists approximate and exact free energies of activation. In this case the immediate impression from Tab. II is the greatly improved agreement between approximate and "exact" results. Though the symmetry is still lost by ξ 2 , 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 ∆F ‡ RP 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 ∆F ‡ RP 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 ∆F ‡ RP 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(NH 3 ) 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][39][40][41][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 8ring (8 silicon sites) window and into cavity B. Millan et al. 43 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 8-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 4 Si and 2 O atoms of the ring that remain nearly in the same plane.
Our primary purpose is to analyze the data of Millan et al. 43 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 44 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. 45 The weights were used to re-compute the FEP and compare it with the result of Millan et al. 43   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. 43 is excellent. The probability densities are normalized according to 4 −4 dz e −βA(z) = 1. Millan et al. 43 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. Tab. 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. 43 . 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. 43 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.

C. Chemically Realistic Model -Radical Cyclization
As a second chemical example, we consider the intramolecular cyclization of the 5-hexenyl radical (see Fig. 8), a radical clock reaction. 46 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, ∆F ‡ exp (300 K) = 42 ± 4 kJ/mol. 47 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. 50 We employed WTM-eABF 51-53 as enhanced sampling algorithm. The unbiased weights were recovered with the recently developed combination of eABF and MBAR. 54 Details of the simulation are given in the supplementary material. The FEP (Fig. 9a) 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. 9b) 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 Tab. 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.

VI. DISCUSSION AND CONNECTION TO PRIOR WORK
This study is not the first work to present expressions for the rate constant and activation free energy based on transition state theory. [29][30][31][32][33][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 ∆F ‡ RP 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) = ΩR dz ρ(z), depend only on ρ(z), the marginal probability density that the CV ξ(x) takes the value z. The third, λ ξ zTS , can be rewritten 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 −1 ξ zTS 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 36 the rate constant is expressed in terms of canonical partition functions for reac-tant 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 −1 ξ zTS is readily expressible explicitly in terms of the discrete masses. 56 In general, however, the CV-conditioned ensemble average must be computed.
The "gauge-independent geometric" free-energy profile, given by has been proposed 24,25 as an alternative to the "standard" FEP (eq. (32)). Since the geometric FEP at the transition point is related to ∆F ‡ RP according to it is also referred to as the "kinetic" free-energy profile. 57 On one hand, like A(z), A G (z) cannot alone provide ∆F ‡ RP . On the other, unlike A(z), A G (z) cannot alone furnish ∆F RP . The essential reason is that e −βA G (z) is generally not a probability density, whereas e −βA(z) always is.
We remark on an apparent inconsistency in the dimensions of terms in eq. (31), as noted in Ref. 57. We observe that the dimensions of ρ(z) are those of ξ −1 and the dimensions of λ ξ are those of ξ. The argument of the logarithm is therefore dimensionless, as it should be. Thus, there is no inconsistency. It appears only because of the tendency to overlook that the definition of the FEP includes an implicit scaling factor, which is unfortunately rarely, if ever, pointed out. The same remarks apply as well to the geometric FEP.

VII. CONCLUSION
Our applications of the exact formula for the activation free energy demonstrate how significant errors can arise when ∆F ‡ RP is approximated simply by the difference between the values of the FEP at the transition state and reactant.
The often employed procedure to obtain ∆F ‡ RP solely from the FEP (by taking the difference between the values at the transition state and reactant (eq. (34))) is an approximation. If ρ(z) is strongly peaked in the vicinity of the minimum of R (i.e, at low temperature and small effective mass m ξ ), then eq. (34) may be satisfactory (see Section V C). However, it is especially questionable when the temperature is high, m ξ is large, and the barrier of the PES between R and P is low (see Section V B).
The exact formula for ∆F ‡ RP (eq. (30)) assumes implicitly that the CV is good (i.e., it is orthogonal to the separatrix). According to our study of the two-dimensional model PES with a systematically variable CV, as the CV becomes less good the reliability of ∆F ‡ RP decreases markedly, while that of the reaction free energy ∆F RP is only slightly affected. We conclude that one must choose the CV with considerable caution in order to achieve the same accuracy for both kinetic and thermodynamic properties.
The exact formulas for ∆F ‡ RP (eq. (30)) and ∆F RP (eq. (42)) depend only on CV-conditioned ensemble averages, which are readily available from enhanced-sampling simulations via reweighting techniques. 45,54,[58][59][60][61][62][63] Therefore, it should be more convenient to use these formulas than to resort to alternative special sampling strategies such as infrequent metadynamics. 64,65 In light of the results of the present study and those of our prior work 23 , we recommend less reliance on the FEP alone and more on the exact formulas, which can be easily evaluated from data provided by commonly employed advanced-sampling algorithms. The exact formulas are more reliable and can be clearly related to experimental data. In this regard we agree with Ref. 34 that use of the FEP alone should be discouraged, except we think that ∆F ‡ RP is a better touchstone for comparison between theory and experiment than the rate constant itself.

SUPPLEMENTARY MATERIAL
The supplementary material contains the following: (1) proof that eq. (14) yields the frequency of crossing the dividing surface; (2) proof of the equivalency of eqs. (12) and (13)

Conflict of Interest
The authors have no conflicts of interest to disclose.

DATA AVAILABILITY STATEMENT
The data that support the findings of this study are available from the corresponding author upon reasonable request. Starting with eq. (14) of the article, we apply the chain rule of differentiation to obtain We now utilize the property of the Dirac distribution 66 where f (t i ) = 0 and (df /dt) t=ti = 0, to recast eq. (A2) as where N τ is the number of zeroes of ξ(t) − z TS on the interval [0, τ ], which is equal to the number of times ξ(t)− z TS changes sign during the interval. Hence, N τ /τ is just the frequency of crossing the dividing surface.
Appendix B: Proof of the Equivalency of Eqs. (12) and (13) We assume that in general the CV is "good" in that it distinguishes properly between R and P (i.e., no configuration of R has the same value of the CV as a configuration of P). Moreover, implicit in eq. (13) is the assumption that the value of z TS − ξ(x) is positive for configurations of R and negative for those of P. It follows that where the Heaviside function is 1 in the domain of R, where its argument is positive.
Assuming the system to be ergodic, we can replace the time average in eq. (S2) with the ensemble average (B2) We note that since ν is the frequency of crossings due to both forward and reverse reactions, the absolute value of the rate of change of the CV is necessary to prevent cancellations of forward and reverse contributions. Further, because the system is taken to be in thermodynamic equilibrium, the forward and reverse reactions occur with the same frequency. Hence, we can simply count reactions in one direction, say forward from R to P, whereξ > 0. Then the absolute value ofξ becomes unnecessary. The sign of the velocity is enforced by introducing a Heaviside function. Thus, we have This is precisely the numerator of the expression in eq. (13) in the article. Dividing eq. (S7) by (S5) finally gives Appendix C: Analytical One-Dimensional Models of Section III Here we treat a one-dimensional system consisting of a single particle of mass m moving on a PES with two minima separated by a maximum. We consider two model PESs, comparing the approximate and "exact" free energies of activation.
We take the PES of the first to be a square well specified piecewise by Assuming that δ L, B > R , and B > P , and taking the CV to be ξ(x) = x, we derive the probability density Hence the probability of observing R is According to eq. (30), we have for the "exact" free energy of activation where we use the relation λ ξ zTS = h 2 /2πmk B T . For the second model we consider a double-well PES having minima of R at x R,min and P at x P,min , separated by a maximum of B at the transition state. We approximate this PES about the minima by the harmonic-oscillator (HO) approximation (e.g., U (x) ≈ U HO (x) = R + k/2(x − x R,min ) 2 , where the force con- ). We again take the CV to be ξ(x) = x. Thus, the probability of observing R is where we approximate the probability density in the domain of R by ρ HO (ξ) = Z −1 e −βUHO(x) . Using eq. (30), we obtain the "exact" activation free energy where we again invoke the relation < λ ξ > zTS = h 2 /2πmk B T . According to eq. (34), the approximate activation free energy for both models is given by: for both models. Comparing eq. (C7) with eq. (C4) and with eq. (C6), we see that the difference between "exact" and approximate activation free energies is, respectively yields m = a/ (1 − a), which, when substituted back into eq. (D7), gives e S(θ) = (a − 1)e x + ae y

Computation of the Probability Density
The marginal probability density is given by To facilitate the evaluation of the double integrals, we transform from Cartesian coordinates to the orthogonal coordinates defined by That ∇q 1 · ∇q 2 = 0 is manifest. The inverse transformation is where d = a 2 + (1 − a) 2 . Hence, the Jacobian is J = ∂x/∂q 1 ∂x/∂q 2 ∂y/∂q 1 ∂y/∂q 2 From eq. (D10) we have ρ(z) = Z −1 dq 1 dq 2 |J| e −βU (x,y) δ(q 1 − z) where the Cartesian coordinates that are the arguments of the PES are given in terms of q 1 = ξ and q 2 by eq. (D14). Using eq. (D15) and the definition of d, we get |J| = 1/d. Hence, ρ(z) = dq 2 e −βU (x,y) dq 1 dq 2 e −βU (x,y) For each simulation corresponding to a given temperature and particle mass, ten independent Langevin dynamics simulations were carried out with a friction constant of 1 ps −1 , a time step of 1 fs, and a total time of 10 ns. The system was propagated using the velocity Verlet algorithm. Crossing frequencies and activation free energies were computed for each simulation independently, only the final values were used for averages and estimation of the standard deviation.

Use of the Heaviside Function
Approximating the time derivative by the forward finite-difference formula, we rewrite eq. (A1) as where the summation on i is over consecutive time steps of length ∆t, N f is the number of steps of the MD simulation, and τ = (N f − 1)∆t is the duration of the MD trajectory. That the expression in eq. (E1), which is straightforward to implement, yields a proper count may be seen as follows. If, at time step i, ξ(t i ) and ξ(t i−1 ) are both greater than or less than z TS , the contribution is zero, since z TS is not crossed during the step. If, on the other hand, ξ(t i−1 ) < z TS and ξ(t i ) > z TS , or ξ(t i−1 ) > z TS and ξ(t i ) < z TS , then the contribution is 1, as z TS is crossed in one direction or the other during the step. We note that ∆t must be sufficiently small that highly frequent crossings are not inadvertently missed.

Use of the Dirac Delta Function and the Atomic Velocities
We consider here an alternative approach to the computation of ν. We begin by recasting eq. (A2) as Here we express the rate of change of the CV aṡ where (∇ x ξ) T = (∂ξ/∂x 1 , ∂ξ/∂x 2 , . . . , ∂ξ/∂x 3N ) is the 3N-dimensional gradient. Discretizing the integration on time, we rewrite eq. (E2) as In practice, of course, the duration of the MD simulation, and therefore the number of time steps, are finite. It is very unlikely that during a finite simulation ξ(t i ) is ever exactly equal to z TS . Hence, almost every configuration x(t) gives zero contribution. To circumvent this problem we introduce a continuous function to represent the delta function approximately. We begin by defining the continuous approximation to the Heaviside function where α is a positive real number having dimension reciprocal length. (Observe that we can formally express the true Heaviside "function" by Θ(x) = lim α→∞ Θ(x; α).) The corresponding delta function is given by the derivative δ(x; α) = dΘ(x; α) dx = α e −αx (1 + e −αx ) 2 . (E6) Note that δ(x; α) satisfies exactly the relation ∞ −∞ dx δ(x; α) = 1 .
We can now rewrite eq. (E4) as where δ[ξ(t i ) − z TS ; α] is approximated by eq. (E6). The index α emphasizes the dependence of the crossing frequency on this methodological parameter.
The quality of this approximation depends on the choice of α. On one hand, if α is too small, MD frames that are far away from the dividing surface contribute significantly, thus yielding too large ν. On the other hand, if α is too large, all frames are weighted so lightly that ν is too small. The influence of α is shown graphically in Fig. 10. Note that we plot the magnitude of the difference so that small differences are visible on the logarithmic scale. One can clearly see that extreme choices of α can lead to errors in ν of up to a factor of 100. As predicted, ν α tends to zero for very large α, which is reflected by the curves approaching 1 (10 0 ). The optimal choice of α seems to be between 10 2 and 10 3 . Even choices in the range 10 1 to 10 4 generally yield ν α 's which are around 0.