For chemical reactions that occur via the rearrangement of atoms from a configuration about one minimum (reactant, R) of the potential energy surface (PES) to a configuration about another minimum (product, P), an exact relation between the Helmholtz reaction free energy (ΔFRP) and the free-energy profile (FEP) can be derived. Since the FEP assumes a form similar to that of the PES along the minimum energy path between R and P, there is an unfortunate tendency to regard the FEP as the “free-energy” analog of the minimum energy path and consequently to equate ΔFRP to the difference between the values of the FEP at the minima corresponding to R and P. Analytic treatments of one- and two-dimensional models are presented that show how this mistaken idea leads to errors. In effect, treating the FEP by analogy with the minimum energy path neglects the role of entropy. The FEP is a function of a collective variable (CV), which must be chosen to describe the course of the rearrangement consistently with the exact relation between ΔFRP and the FEP. For large systems of common interest, the PES is often so complex that a straightforward way of choosing a CV is lacking. Consequently, one is forced to make an educated guess. A criterion for judging the quality of the guess is proposed and applied to a two-dimensional model.
I. INTRODUCTION
We are concerned in this article with the theoretical treatment of molecular processes (e.g., unimolecular isomerization) that involve the rearrangement of the atoms from one configuration (reactant, R) to another (product, P). The rearrangement can be represented as a chemical reaction
where R and P are treated as distinct molecules. The configurations of the molecules are localized about minima of the potential energy surface (PES) and these minima are separated by a barrier that must be surmounted as the atomic configuration changes from that of R to that of P. We focus on the property that is typically used to characterize the reaction: the change in free energy that accompanies the conversion of R to P or the Helmholtz reaction free energy, ΔFRP, which can be expressed analytically in terms of the partition functions of R and P.
Thus, in principle, the straightforward route to ΔFRP is via the partition functions, although, in practice, the calculation of reliable partition functions is very difficult. In some situations, one may wish to employ economically taxing numerical methods in order to explore precisely only the domains of the PES in the vicinity of the minima about R and P. In such cases, one may then invoke the ideal-gas, rigid-rotor, harmonic-oscillator (IGRRHO) approximation1 for the partition functions from which one obtains an estimate of ΔFRP.2,3
Since direct calculation of accurate partition functions is generally infeasible, alternative approaches for the determination of ΔFRP have been proposed. One of these relates ΔFRP to the so-called free-energy profile (FEP),4 A(ξ), which is a function of a collective variable (CV), ξ(x) (or reaction coordinate), that describes the course of the rearrangement. Other frequently used terms for the FEP are the free energy surface, which is often used when more than one CV is employed, and the potential of mean force. The FEP can be expressed in terms of the probability density ρ(ξ) of observing the system in a configuration corresponding to a particular value of the CV. The CV is ideally chosen so that R and P correspond to non-overlapping domains of the CV.
In principle, ρ(ξ) can be computed by means of molecular dynamics (MD) for the given PES determined by either force-field approximations or quantum mechanics. However, the complex, multidimensional systems of the greatest interest usually have a high potential barrier separating R and P. In this case, the standard MD procedure fails to properly explore the whole configuration space. To overcome this barrier, one usually resorts to advanced sampling algorithms. A number of computational techniques based on various definitions of the FEP have been developed and applied to complex biomolecular reactions. These include, but are not limited to,5,6 the following: imposition of a pre-determined biasing potential (i.e., umbrella sampling)7–10 coupled with the post hoc removal of the bias,11–14 use of adaptive biasing potentials that converge to the inverse of the sought FEP,15–17 and use of adaptive biasing forces.18–27
The apparent success of these methods is reflected in the vast number of publications that make use of them. In many instances, however, the proper relation between ΔFRP and the FEP is confused. One pervasive tendency is to interpret the FEP by analogy to the minimum energy path, i.e., the minimum potential energy along the reaction path between R and P (see Fig. 1). Although we could cite a number of articles that exhibit such a confusion, we prefer to mention a few that employ the relation correctly.21,28–30
Schematic FEP [A(ξ)] with two minima separated by a maximum. Domain to the left of the maximum corresponds to the reactant (R) and that to the right to the product (P). Incorrect estimation of reaction free energy highlighted in red and correct relation in green.
Schematic FEP [A(ξ)] with two minima separated by a maximum. Domain to the left of the maximum corresponds to the reactant (R) and that to the right to the product (P). Incorrect estimation of reaction free energy highlighted in red and correct relation in green.
The principal purpose of this article is to clearly demonstrate the connection between ΔFRP and the FEP, especially the quantitative relation between the two. In Sec. II, we develop the theory. To illustrate the types of errors that occur, we examine in Sec. III one- and two-dimensional models that can be handled analytically. Section IV summarizes our conclusions.
II. THEORY
A. Chemical equilibrium and reaction free energy
We assume that the system consists of a dilute solution of species R and P undergoing the reaction represented by Eq. (1). The condition for chemical equilibrium is
where μα is the chemical potential of species α (=R, P). Since we assume the solution to be dilute, we can write1
where T is the absolute temperature, kB is the Boltzmann constant, qα is the canonical molecular partition function, and Nα is the number of species α. In the present theory, we take the species to be defined by the region of the configuration space it “occupies,” designated by Ωα. Hence, in the classical limit, we have
Here, denotes the configuration of atoms in Cartesian coordinates and dx is the corresponding 3N-dimensional volume element; denotes the vector of momenta conjugate to the coordinates, h is Planck’s constant, β ≡ 1/kBT, is the Hamiltonian
M is the diagonal mass matrix, and U is the potential energy surface (PES).
We assume that ΩR and ΩP constitute the whole 3N-dimensional configuration space. ΩR is separated from ΩP by a (3N − 1)-dimensional dividing surface. In principle, the dividing surface can be chosen arbitrarily. However, it is normally taken to contain the ridge of the potential barrier of the PES that separates the minima corresponding to R and P, in which case it is termed the separatrix . From the definition, we conclude that on , the gradient of the potential energy, , and the surface normal are orthogonal to each other.
Inserting the expression for into Eq. (4) and performing the integration on momenta, we obtain
where the second line implicitly defines the configuration integral Zα of molecule α,
and is the standard de Broglie thermal wavelength of atom i. Only those configurations belonging to Ωα contribute to Zα, which is the effective 3N-dimensional volume of the configuration space occupied by α. Since ΩR and ΩP constitute the whole configuration space available to the system, the total volume is given by
The fraction of the total volume of the configuration space occupied by α,
is the probability of observing an α molecule (i.e., the probability that a molecule chosen at random from the solution in equilibrium is an α).
which means that the ratio of the number density [α] of the product to that of the reactant, fixed at equilibrium, is equal to the ratio of the respective configuration integrals, as well as the ratio of the respective probabilities of observation. The Helmholtz free energy of a single molecule of α is
For pedagogical reasons, we choose here the symbol F to distinguish the free energy from the FEP, which we denote by A. Solving Eq. (11) for qα and substituting the expression into Eq. (10) gives
where
is the free-energy difference between the product and the reactant (i.e., the reaction free energy).
B. Free-energy profile and its connection with reaction free energy
The atomic rearrangement can be described by a collective variable (CV), denoted by ξ(x), which is a function of some subset of the Cartesian coordinates. The CV is chosen so that it has two non-overlapping domains that correlate with the domains of R and P, as defined by the dividing surface. In this circumstance, we say that the CV discriminates between R and P. Variation in the CV generates the FEP that corresponds to the rearrangement. Various definitions of the FEP have been proposed. Here, we take the FEP to be given by18,31
where
We note that the symbol ξ, without the argument x, stands for a particular value of the CV. The second line of Eq. (15) shows that ρ(ξ) is the normalized probability density (in the variable ξ) of observing the molecule in a configuration x with ξ(x) = ξ.
Equation (14) can be recast as
which, upon integration on ξ over the domain of molecule α, gives
which is the sought relation between the FEP, A(ξ), and the reaction free energy, ΔFRP.
This correct relation is sometimes overlooked in favor of an intuitive notion that ΔFRP is simply equal to the difference between the two minima of A(ξ) (see Fig. 1). This is true if ρ(ξ) is so strongly peaked about the minima that
where A(ξα,min) is the minimum value. In Sec. III, we explore the limits of this approximation (i.e., treating the FEP by analogy to the PES and equating ΔFRP to the difference between the minima of the FEP).
We note that all relations in this section for the case of a single CV can be generalized to the case of multiple CVs.
C. Criterion for a well discriminating CV
The validity of the key formula [Eq. (18)] depends on the assumption that the CV distinguishes correctly between R and P, as defined by the dividing surface. Hence, knowledge of the separatrix is critical to the choice of the CV. However, for large, complex systems, we usually have no specific information about and must, therefore, rely on our chemical intuition to guess the CV. Such intuitive CVs may work, but, as shown in Sec. III C, they can lead to significant errors.
We recall the requirement that a “good” CV should satisfy: its gradient should be normal to . This condition of orthogonality is31
However, our lack of knowledge of renders this relation alone useless. To proceed, more or less blindly, we propose to monitor the mean value of along the CV, which can be expressed as
where denotes the ensemble average over all configurations with ξ(x) = ξ. We choose the absolute value of the scalar product, instead of the scalar product itself, in order to prevent the cancellation of negative and positive values, which tends to occur near the extrema of the FEP. It can be shown19,32–34 that the slope of the FEP is given by
In the case that the second term of the right member of Eq. (22) vanishes for all values of ξ, e.g., if the CV is a linear combination of Cartesian coordinates, all extrema of the FEP [and of ρ(ξ)] coincide with roots of .
We remark on a noteworthy contrast between [i.e., the first term of the right member of Eq. (22)] and . Unless ∇ξ is constant, the zeroes of do not coincide with those of , whereas the zeroes of always coincide with those of . Now, in order for to vanish, the absolute value |∇U · ∇ξ| must be zero at every point x for fixed ξ, that is, ∇ξ must be orthogonal to ∇U everywhere on the isosurface corresponding to ξ(x) = ξ (we assume that ∇U is not zero everywhere on that surface). The same is true for . As the scaling of |∇U · ∇ξ| by 1/|∇ξ|2 does not alter the roots of , we omit it from the definition of D(ξ).
D(ξ) vanishes on if ∇ξ is orthogonal to everywhere or, equivalently, if ∇ξ is everywhere parallel to the surface normal. If orthogonality is unfulfilled, then D(ξ) has at best a minimum in the vicinity of an extremum of A(ξ). However, it may happen that on account of the symmetry of the PES, ∇ξ · ∇U = 0 everywhere on an isosurface other than . In such an instance, the vanishing of D is a spurious indication of . Distinguishing anomalous zeroes of D from the one corresponding to the true will still require one’s chemical intuition.
D. The IGRRHO approximation
In case the PES is known in regions of the configuration space around the minima of R and P, one can often obtain a good estimate of the reaction free energy using the IGRRHO approximation to compute the partition functions qR and qP [see Eq. (12)]. The procedure is summarized as follows: Considering one of the minima, we expand the PES about the minimum configuration x0 in Taylor’s series, neglecting terms of degree higher than quadratic. In terms of mass-weighted coordinate displacements , we obtain
where U(0) is the potential energy at the minimum and U″(0) is the Hessian matrix {i.e., }. Since M−1/2U″(0)M−1/2 is real symmetric, it can be diagonalized by a unitary transformation. Thus, we define the normal coordinates , converting Eq. (23) to
Taking W to be the required unitary transformation, we obtain
where Qi are the normal modes and are the corresponding eigenvalues. The Hamiltonian is, therefore, separable into 3N independent terms. Five or six of the eigenvalues are zero; they correspond to the three translations of the whole molecule and to the two (for linear molecules) or three (for non-linear molecules) rotational modes. The remaining (3N − 5) or (3N − 6) modes correspond to intramolecular vibrations. The partition function can then be factored as
The translational partition functions cancel in the quotient qP/qR; the rotational partition functions, which depend on the principal moments of inertia of the rigid rotor, are generally different for R and P. The same is true for the vibrational partition functions, given in the classical limit by qvib = ∏ikBT/ℏωi, where ωi denotes the fundamental radial frequency of the ith vibrational mode.
The IGRRHO approximation is commonly invoked when quantum mechanics is used to determine the PES in vacuum or an implicit solvent.35–38
III. PITFALLS IN THE INTERPRETATION OF THE FEP
We turn now to specific applications of the theory that illustrate potential traps in the interpretation of the FEP. To render the presentation lucid, we employ one- and two-dimensional models that can be handled largely analytically. Where required, numerical analysis is handled by NumPy.39 All graphs are generated with Matplotlib.40 The Jupyter41 Notebook for all computations is given in the supplementary material.
A. Asymmetric one-dimensional PES with equally deep minima
We consider first a single particle that moves in one dimension on the asymmetric model PES
where x is the dimensionless coordinate and the positive parameter ϵ scales the energy. The PES, plotted in Fig. 2, has equal minima at xR,min = 2.084 43 and xP,min = 6.000 00, which are separated by a barrier with a maximum at xmax = 3.479 45. The domains x < xmax and x > xmax correspond to R and P, respectively.
Since the FEP is merely shifted vertically by a constant with respect to the PES, it has equal minima at ξ = xR,min and ξ = xP,min. Hence, were we to treat the FEP in the same way as the PES, we would conclude that the reaction free energy ΔFRP = A1(ξ = xP,min) − A1(ξ = xR,min) = 0, according to the procedure indicated in Fig. 1 in red highlight. However, the disparity between the widths of the potential wells associated with R and P (see Fig. 2) suggests that P is observed more frequently than R [i.e., ]. Indeed, numerical integration of ρ(ξ) according to Eq. (17) yields and , if, for convenience, we set ϵ = kBT. From Eq. (18), we get .
B. Symmetric one-dimensional PES with equal minima
We look again at the one-dimensional motion of a single particle. The model PES
which is plotted in Fig. 3(a), has two equal minima separated by a barrier with a maximum at xmax = 5. Since the PES diverges as x approaches 0 or 10, the particle is confined to the domain 0 < x < 10. R and P, respectively, correspond to the domains 0 < x < 5 and 5 < x < 10.
(a) Plot of the symmetric model PES U2(x)/ϵ [see Eq. (29)]; vertical dashed line indicates maximum at xmax = 5. (b) Plot of the FEP A2(ξ)/ϵ [see Eq. (33)]. For convenience, we set ϵ = kBT. Horizontal dashed lines indicate minima and vertical orange line indicates maximum ξmax = 1/4.786 = 0.2089. (c) Plot of D(ξ) [see Eq. (34)].
(a) Plot of the symmetric model PES U2(x)/ϵ [see Eq. (29)]; vertical dashed line indicates maximum at xmax = 5. (b) Plot of the FEP A2(ξ)/ϵ [see Eq. (33)]. For convenience, we set ϵ = kBT. Horizontal dashed lines indicate minima and vertical orange line indicates maximum ξmax = 1/4.786 = 0.2089. (c) Plot of D(ξ) [see Eq. (34)].
We now take the CV to be
The form of ξ(x), which appears to invert, or distort, the Cartesian space, may seem unlikely to be applicable to a realistic system. Nevertheless, such “distorting” CVs are indeed used (see, for example, the “coordination” CV in the program package PLUMED42). From Eq. (15), we have
Transforming the variable of integration to q = x−1, we can rewrite Eq. (31) as
Substituting this expression for ρ(ξ) into Eq. (14), we obtain the FEP
The form of this expression already suggests how the FEP viewed as a function of ξ = 1/x is a distortion of the PES viewed as a function of x. The plot of A2(ξ) in Fig. 3(b) shows the distortion explicitly. The domains of R and P are reversed in the FEP, which is strongly asymmetric. The sides of the potential well associated with P are much steeper than those of the potential well associated with R. Furthermore, the minima are not equal and the maximum of the barrier that separates R and P does not occur at precisely the inverse of the position xmax = 5 of the maximum of the barrier of the PES but rather at the slightly larger value ξmax = 0.2089, as indicated by the orange vertical line in Fig. 3(b). All of these elements of distortion evident in the FEP reflect the strongly distorting power of the CV, which effectively compresses large values of x while it stretches small values of x.
We note that if we were to regard the FEP in the same manner as the PES, we would estimate the reaction free energy from Fig. 3(b) to be . However, it is clear from Fig. 3(a) that by symmetry. Indeed, the correct treatment as delineated in Eq. (18) gives .
For complex, multi-dimensional systems, we generally lack knowledge of the (3N − 1)-dimensional dividing surface (separatrix) of the 3N-dimensional PES. To handle such situations, one usually makes an educated guess for the CV, uses it to compute the FEP, and takes the maximum of the FEP to correspond to the transition state (TS). The danger inherent in this procedure is illustrated for the present one-dimensional model by Fig. 3(b), in which the gray and orange vertical lines indicate a discrepancy between the maxima of the FEP and the PES. Thus, even correct use of the relation between the reaction free energy and the FEP [Eq. (18)] leads to an error in ΔFRP, since the domains of integration Ωα are slightly in error.
C. Choice of a good CV for multi-dimensional systems
To further illustrate the difficulties encountered in finding a “good” CV for a multi-dimensional system (i.e., one that discriminates properly between R and P), we consider a single atom moving on the model two-dimensional PES
where
We note that U3 is symmetric under reflection about the line x = 5. The contour plot in Fig. 4 shows that U3 has two minima. We take the portion of the x–y plane above the red line (y = 2.2), which has a broad, shallow minimum, to constitute ΩP and the portion below the red line, which has a narrow, deep minimum, to constitute ΩR. Thus, the red line is the dividing surface.
Contour plot of the PES U3(x, y) [see Eq. (35)]. Numbers on contours give energy in units of ϵ. The red line is the dividing surface between R (below line) and P (above line). The green line indicates the minimum energy path.
Contour plot of the PES U3(x, y) [see Eq. (35)]. Numbers on contours give energy in units of ϵ. The red line is the dividing surface between R (below line) and P (above line). The green line indicates the minimum energy path.
As the dimensionality of a system increases, so does the difficulty of choosing a CV that differentiates between R and P. One usually has some prior information (or chemical intuition) upon which to base a choice of the CV. However, for the purpose of demonstrating problems with the choice, which already arise for the present two-dimensional model, we suppose that nothing is known of the PES. We are forced to guess the CV blindly. Thus, we assume a CV and compute the probability density using Eq. (15). The FEP A3(ξ) is then given by Eq. (14). We take the relative maximum of A3(ξ) [or relative minimum of ρ(ξ)] to correspond to the dividing surface. Plots of A3(ξ) for three choices of the CV are displayed in Fig. 5. Plots of D(ξ) are also shown. For didactic purposes, we also show plots of . Since ξ(x) is a linear combination of Cartesian coordinates in all three cases, ∇ξ is constant everywhere, and Eq. (22) reduces to . Hence, vanishes at the extrema of A3(ξ), in contrast with D(ξ), which is zero only for isosurfaces on which ∇ξ is normal everywhere to ∇U3 (see Sec. II C). For the sake of simplicity of notation, we set T so that ϵ = kBT.
FEPs for three choices of the CV applied to the PES U3(x, y): (a) ξ(x) = y; (b) ξ(x) = x; and (c) ξ(x) = x + y. Position of maximum on the FEP is indicated by vertical dotted line. Probability that ξ is less than or greater than this maximum indicated by s.
FEPs for three choices of the CV applied to the PES U3(x, y): (a) ξ(x) = y; (b) ξ(x) = x; and (c) ξ(x) = x + y. Position of maximum on the FEP is indicated by vertical dotted line. Probability that ξ is less than or greater than this maximum indicated by s.
For the purpose of comparing the other two “blind” choices, we include the best guess for the CV, namely, ξ(x) = y [see Fig. 5(a)], which increases in the direction ∇ξ = ey. From Eqs. (15), (35), and (36), we find that ρ(ξ) ∝ g(ξ), and from Eq. (14), we get A3 = −kBT ln g(ξ) + const. Setting and solving the resulting non-linear equation, we find a relative maximum at ξmax = 2.2, which confirms our best guess. We also find relative minima at ξ ≈ 1.4 and ξ ≈ 6.1. Using Eqs. (21), (35), and (36), we obtain . The plots in Fig. 5(a) show that D(ξ) vanishes at the extrema of A3(ξ). One of these, namely, the relative maximum at ξmax = 2.2, corresponds to the orthogonality condition on . The other two are spurious, deriving from the symmetry of the PES, that is, there are surfaces at y ≈ 1.4 and y ≈ 6.1 on which is zero.
Next, we consider the CV ξ(x) = x, which increases parallel with the true dividing surface. In this case, ρ(ξ) ∝ f(ξ) such that A3(ξ) = kBT(ξ − 5)2/4 + const. Hence, the plot of A3(ξ) is a parabola with minimum at ξ = 5 [see Fig. 5(b)]. That A3(ξ) has no relative maximum shows immediately that this CV discriminates not at all between R and P. Indeed, it is the worst possible guess. For D(ξ), we get D(ξ) = ϵ|ξ − 5|/2, which tells us that D(ξ) vanishes at ξ = 5, even though this is not the true dividing surface. This anomaly is traceable to the symmetry of the PES. On the surface ξ = 5, everywhere. Therefore, over the entire surface ∇ξ · ∇U3 = 0.
We next examine the CV that is a compromise between the best and worst: ξ(x) = x + y, which increases in the direction of ∇ξ = ex + ey. The plot of A3(ξ) in Fig. 5(c) exhibits a weak relative maximum at ξmax = 7.45. For ξ < 7.45, the greater contribution arises mainly from ΩR; for ξ > 7.45, it comes mostly from ΩP. Hence, this CV distinguishes between R and P better than the worst CV. As the plot in Fig. 5(c) shows, D(ξ) far exceeds zero at ξmax = 7.45 corresponding to the relative maximum of A3(ξ). As D(ξ) has no roots, we conclude that the “compromised” CV is better than the worst but far from the best. However, unlike D(ξ), has three roots, one at each of the local extrema. Therefore, it is evidently not a good measure of the quality of the CV.
It is noteworthy that for the best CV, D(ξ) is strongly peaked. For bad CVs, D(ξ) is smooth and non-zero everywhere except for the worst CV.
D. Application of the IGRRHO approximation
To illustrate the utility of the IGRRHO approximation for the estimation of the reaction free energy, we apply it to the models considered in Secs. III A and III B. These one-dimensional models have only a single vibrational degree of freedom. Therefore, as demonstrated in the supplementary material, the expression for the reaction free energy in the HO approximation is
where kR and kP are the force constants for R and P, respectively. Table I compares the results of various approximations with the “exact” results [Eq. (12)]. Figures S1 and S2 display the PES with parabolas fitted to the minima.
It is apparent that though the difference between the FEPs at the minima disagrees markedly with the “exact” value, the HO approximation agrees quite well with the exact. We note that the precise agreement of the exact and HO values for U2 is due to the symmetry of the PES. As Fig. S2 shows, the shape of the minima of U2 is not as well fitted by the parabola as the shape of the minima of U1 (Fig. S1).
E. Confusion between free energy diagrams and free energy profiles
A common practice is to portray the course of a chemical reaction from reactants (R) to products (P) via the transition state (TS) by a diagram, as shown in Fig. 6(a). The horizontal lines indicate the free energy, as computed, for example, by the IGRRHO (see, e.g., Refs. 43 and 44). In particular, in the case of a multi-step reaction, the intermediates and transition states are often connected by a solid curve, as shown in Fig. 6(b). This gives the impression that the curve represents a one-dimensional PES, which it does not. The curve itself is merely an artistic flourish. Only the horizontal lines at the local extrema have physical significance, namely, the free energy of the state. The difference between the values indicated by these lines at the minima is the reaction free energy.
(a) Free energy diagram of a reaction, (b) same as (a) with added spline for visual purpose only, and (c) free energy profile that is fundamentally different from the curve in (b).
(a) Free energy diagram of a reaction, (b) same as (a) with added spline for visual purpose only, and (c) free energy profile that is fundamentally different from the curve in (b).
It is easy to confuse a free energy diagram [Fig. 6(b)] with the FEP, illustrated in Fig. 6(c). The FEP, A(ξ), is related to the probability density ρ(ξ) by Eq. (14) and its gradient is the mean force acting along the CV, ξ(x). The FEP is not, however, a free energy, as defined in Eq. (11) and as indicated by the horizontal lines in Fig. 6(b). As noted in Sec. II B, if the density ρ(ξ) is sharply peaked about the minima of R and P, then the difference between the FEP at the minima may be a good approximation to ΔFRP [see Eq. (19)]. However, whether ρ(ξ) is sufficiently strongly localized about the minima must be ascertained in each particular instance.
Another widespread practice adding to the confusion, as depicted in Fig. 1, is the denotation of the FEP by ΔA, or even ΔG, instead of A(ξ). This suggests that the FEP itself corresponds somehow to a free energy difference, which it does not. Although many articles in the literature employ this misleading notation, we prefer not to cite them, for obvious reasons.
IV. CONCLUSION
For chemical reactions that consist of the rearrangement of atoms from a configuration about one minimum (reactant) of the PES to a configuration about another minimum (product), we establish a precise relation between the reaction free energy (ΔFRP) and the so-called free-energy profile, A(ξ), which is a function of a collective variable ξ(x) that describes the rearrangement. A properly chosen CV discriminates between R and P, permitting a valid estimate of ΔFRP.
Because the FEP is similar to the PES along the minimum energy path between R and P (i.e., two relative minima separated by a relative maximum), one is tempted to view the FEP as the “free-energy” analog of the minimum energy path and consequently to equate ΔFRP to the difference between the values of the FEP at the minima corresponding to R and P. Using simple one-dimensional models, we demonstrate how this misguided notion leads to errors. For example, the mere fact that the values of the FEP at the minima of R and P are equal does not mean that ΔFRP = 0. This can happen when the PES is asymmetric so that the system is more likely to be found in the configurational space of R than in that of P or vice versa. Thus, treating the FEP by analogy with the minimum energy path neglects entropic effects.
A valid but non-linear CV applied to a highly symmetric PES (e.g., one for which the correct value of the reaction free energy is ΔFRP = 0) may yield a highly asymmetric FEP (e.g., one with unequal values at the minima corresponding to R and P). Thus, taking directly the difference between the values of the FEP at the minima of R and P obviously yields an erroneous ΔFRP ≠ 0. Even though the CV strongly distorts the PES, it nevertheless yields an FEP that can give the correct ΔFRP when properly treated according to the correct relation.
For large, complex systems, we usually lack knowledge of the ideal CV. As a consequence, we are forced to hazard a guess for the CV based on our chemical intuition. We propose a criterion to judge the quality of the guess: monitor the quantity D(ξ) along the CV. It must vanish at the dividing surface (i.e., it must exhibit a cusp near the relative maximum of the FEP). Note that in the case of distorting CVs, the cusps of D(ξ) and the local maximum of the FEP need not coincide. The better the guessed CV is, the closer to zero the D(ξ) is near ξmax. A bad CV will cause D(ξ) to remain positive and be smooth (i.e., its derivative is continuous). For symmetrical PESs, it can happen that D(ξ) also vanishes elsewhere. Such anomalies do not impair its utility.
In situations where the PES is known in regions of the configuration space about the minima of R and P, application of the IGRRHO approximation can yield a reasonably good estimate of ΔFRP. Strengths and limitations of the IGRRHO are discussed elsewhere.45–48
SUPPLEMENTARY MATERIAL
See the supplementary material for (i) a PDF with the details of the results presented in Table I and (ii) a Jupyter Notebook used to generate Figs. 2–5, provided for the interactive use of the reader.
ACKNOWLEDGMENTS
Financial support was provided by the “Deutsche Forschungsgemeinschaft” (DFG, German Research Foundation) (Grant No. SFB 1309-325871075) “Chemical Biology of Epigenetic Modifications.” C.O. acknowledges further support as Max Planck Fellow at the MPI-FKF Stuttgart.
AUTHOR DECLARATIONS
Conflict of Interest
The authors have no conflicts of interest to disclose.
DATA AVAILABILITY
The data that support the findings of this study are available within the article and its supplementary material.