Reaction rate theory

Already in 1940, Hendrik Kramers understood well the mechanism of climbing the reaction barrier in solution as governed by the thermal (Brownian) motions of the solvent. The above mentioned rare event problem did not play a role in the description he used, as the solvent effect was not described by explicit molecules, but rather by Gaussian random fluctuations $f(t)$ and by a linear friction force $- \gamma m \frac{\mathrm{d}x}{\mathrm{d}t}$ working on a reaction coordinate $x$ with a mass $m$. A one-dimensional asymmetric double-well potential $U(x)$ represents the reactant and product well and the separating barrier along by the reaction coordinate, which enter with the solvent forces Newton's equation of motion of $x$ in the form of a Langevin equation:

$$m \frac{\mathrm{d}^2 x}{\mathrm{d}t^2} = -\frac{\mathrm{d}U(... ...\mathrm{d}x} - \gamma m \frac{\mathrm{d}x}{\mathrm{d}t} + f(t)$$ (11)

Here, $\gamma$ is a constant friction rate, and the fluctuating force $f(t)$ denotes Gaussian white noise with zero mean, which obeys the fluctuation-dissipation relation,

$$\left< f(t) \right> = 0 \quad ,$$

$$\left< f(t) f(t^\prime)\right> = 2m \gamma k_BT \delta(t-t^\prime) \quad .$$ (12)

The delta function in 2.11 ensures instantaneous dissipation, which is known as the Markov property. From the time evolution of the probability density $p(x,v,t)$, expressed in a so-called Fokker-Planck equation, he studied the reactant flux over the barrier for low and high friction $\gamma$ (also referred to as damping, coupling or viscosity, and related to the diffusion coefficient via Einstein's relation: $\gamma = k_BT/D$). For moderate-to-strong friction he came to the cornerstone result in rate theory for the reaction rate constant[4]:

$$k_R = \frac{1}{\omega_b} \left( - \frac{\gamma}{2} + \sqrt... ...ight) \Bigg\{ \frac{\omega_R}{2\pi} \exp{[-\beta E_b]} \Bigg\}$$ (13)

where $\omega_R$ and $\omega_b$ the frequencies that confine the parabolic reactant well and barrier top respectively, as they are both of the form $U(x)=U(x^\prime) - \frac{1}{2}m\omega^\prime(x-x^\prime)^2$. In this expression (2.12), the part between the curly brackets gives the transition state theory (TST) result for the reaction rate, so that the expression in front of the curly brackets approximates the correction to the TST result, known as the transmission coefficient $\kappa=k_R/k_\mathrm{TST}$, for the moderate-to-strong solvent friction regime.

Transition state theory[5] is based on two assumptions, namely 1) thermodynamic equilibrium must prevail throughout the entire system for all degrees of freedom and 2) any transition pathway starting from the reactant state which crosses the dividing transition state surface will end up in the product state. For strong (over-damped) friction, $\gamma \gg \omega_b$, the expression for the rate constant simplifies to,

$$k_R = \frac{\omega_b}{\gamma} \frac{\omega_R}{2 \pi } \exp{[-\beta E_b]}$$ (14)

so that the correction factor to TST, $\kappa= \frac{\omega_b}{\gamma}$ becomes inversely proportional to the solvent friction, and $\kappa \to 0$ as $\gamma \to \infty$. From equation 2.12, we can deduce that the transmission coefficient is always smaller than unity for any friction $\gamma$ larger than zero, and that the transition state theory rate $k_\mathrm{TST}$ thus provides an upperbound for the true rate. Kramers simulations showed that the second TST assumption does not hold for a reaction which is strongly coupled to the thermal motions of a solvent, because the strong solvent friction forces the reaction coordinate into something like a random walk at the barrier top so that many transition pathways recross and do not end up in the product state. The result is referred to as spatial-diffusion-controlled rate.

In the low solvent friction limit, Kramers found

$$k_R = p \gamma \frac{I(E_b^+)}{k_BT} \Bigg\{\frac{\omega_R}{2\pi}\exp{[-\beta E_b^+]}\Bigg\}$$ (15)

where $I(E_b^+)$ is the action at the barrier top in the forward direction, $E_b^+$ the barrier top energy with respect to the reactant state and $p$ the probability to thermalize in the product well. This probability equals $p=\frac{1}{2}$ for a symmetric double well potential. The transmission coefficient (again the terms in front of the curly brackets) in the small friction regime is thus proportional to the friction and it goes to zero as $\gamma \to 0$. This friction regime describes the energy-diffusion-controlled rate, because the reaction coordinate moving over the barrier top will recross many times from the product state to the reactant state and back before it loses enough energy to the solvent to thermalize in either stable state. At some intermediate friction, there will be a crossover between the spatial-diffusion regime and the energy-diffusion regime, known as Kramers turnover, which is highly system dependent.

Kramers reaction rate theory predicts that the solvent dynamics will always decrease the reaction rate, or that in the most favorable situation no barrier recrossings take place and the rate constant corresponds with that of Rice-Ramsperger-Kassel-Marcus (RRKM) rate theory, which gives the TST rate as a function of the collision rate for independent polyatomic molecules. However, for general realistic applications, Kramers' theory (and RRKM theory) fails when the time scale of barrier crossing (which is of course much faster than $1/k$) is in the same order or even slower than the time scale of the correlations in the random solvent fluctuations. Moreover in RRKM theory, the reacting polyatomic molecules provide there own (infinitely fast) energy sinks, whereas in realistic applications the coupling of the reaction coordinate with these intramolecular modes can be weaker than the coupling (friction) with the solvent. In certain cases, the solvent dynamics can therefore increase the reaction rate, as it provides an extra energy dissipation source after barrier crossing. An important improvement on Kramers formulation is the generalization to non-Markovian solvents which introduces memory effects in the solvent fluctuations expressed by equation 2.11 (Grote-Hynes theory[6]). An excellent review of improvements on Kramers theory of barrier crossings in many-particle systems is found in ref kramersreview90.

From the phenomenological point of view, the Langevin equation and Fokker-Planck equation methods to describe barrier crossings in the condensed phase have been (and still are) very important for our understanding of solvent effects on the rate constant. A disadvantage of these methods is that they depend heavily on parameters such as the potential felt by the reaction coordinate $U(x)$, the coupling (friction) between the reactants and the solvent $\gamma$ and the dissipative fluctuating force $f(t)$, which are not known for general applications. Moreover, we are also interested in the behavior of the solvent environment during a chemical reaction, which requires the inclusion of explicit molecules via molecular dynamics. We will therefore make use of an alternative microscopic expression for the reaction rate derived from the macroscopic equations (2.1) for which in equilibrium $k_R\left<c_R\right> = k_P\left<c_P\right>$. We introduce a microscopic reaction coordinate $x$ and the Heavyside step-function $\Theta(x-x^\ddagger)$ which indicates whether the system finds itself in the reactant state or in the product state,

$$\Theta(x-x^\ddagger) = \left\{ \begin{array}{ll} 0 & \textrm... ...e of the transition state } x^\ddagger \\ \end{array} \right.$$

The decay of a deviation at time zero, $\Delta\Theta(x(0)-x^\ddagger)$, from the equilibrium $\left<\Theta(x-x^\ddagger)\right>$ can be written as[8]:

$$\Delta\Theta(x-x^\ddagger) = \Theta(x-x^\ddagger) - \left<\Theta(x-x^\ddagger)\right>$$ (16)

$$\Delta\Theta(x(t)-x^\ddagger) = \Delta\Theta(x(0)-x^\ddagger) \exp[-(k_R+k_P)t]$$ (17)

Applying again Onsager's hypothesis and taking the time derivative, we obtain equation 2.17 for the forward reaction rate constant[9,10,11,8]:

$$k_R(t) = \frac{\left< \delta(x(0)-x^\ddagger) \dot{x}(0) \Th... ...x(t)-x^\ddagger) \right>} {\Big<\Theta(-x(0)+x^\ddagger)\Big>}$$ (18)

The rate constant is thus expressed as a time correlation function of the velocity along the reaction coordinate ( $\dot{x}=\mathrm{d}x/ \mathrm{d}t$) at the transition state at time zero, $\delta(x(0)-x^\ddagger) \dot{x}(0)$, with the probability that the system finds itself in the product state at time $t$, normalized by the fraction of transitions that started in the reactant state. A typical appearance of the rate constant as a function of time is shown in figure 2.1. For very short correlation times, $t\approx 0$, the system, crossing the transition state with a positive velocity $\dot{x}(0)$, will find itself in the product state ( $\Theta(x(t)-x^\ddagger)=1$), because it has not yet had the time to recross back. For $t$ equals zero, we thus recover the transition state theory result for the rate constant, $k_\mathrm{TST}$. For larger correlation times, recrossings do take place, and $k(t)$ decreases. For very long times, in the order of the time scale of typical microscopic relaxation times $\tau_\mathrm{relax}$, $k(t)$ decays as $k\exp{[-kt]}$, as the system will have thermalized in one of the stable states and no further barrier recrossings take place until the system is reactivated. Since for this time scale $\exp{[-kt]}\approx1$, $k(t)$ appears to reach a plateau value, which equals the macroscopic rate constant $k$.

Figure 2.1: Reaction flux rate constant as a function of time, as expressed by equation 2.17. Figure adapted from ref chandlerlecture.

In the shape of equation 2.17, apparently we still have to perform one prohibitively long simulation and wait for every occasion that the barrier top is reached to set our stopwatch to $t=0$ and add a measurement of $\Theta(x(t)-x^\ddagger)$ to our statistics, in order to obtain $k$. For practical calculations we therefore rewrite equation 2.17 by multiplying numerator and denominator with
$\left< \delta(x(0)-x^\ddagger) \dot{x}(0) \Theta(\dot{x}(0)) \right> \left< \delta(x(0)-x^\ddagger) \right>$, giving: 2pt

$$k_R(t) = \frac{\left< \delta(x(0)-x^\ddagger) \dot{x}(0) \Theta(x(t)-x^\dd... ...ht>} {\Big< \delta(x(0)-x^\ddagger) \dot{x}(0) \Theta(\dot{x}(0)) \Big>} \times$$

$$\times \frac{\left< \delta(x(0)-x^\ddagger) \dot{x}(0) \Theta(\do... ...ac{\left< \delta(x(0)-x^\ddagger) \right>} {\Big<\Theta(-x(0)+x^\ddagger)\Big>}$$ (19)

The last term at right-hand-side of equation 2.18 is simply the probability to find the system at the transition state relative to the probability of finding the system in the reactant state. The term in the middle gives the average forward crossing speed at the transition state, which equals together with the last term the transition state theory rate constant, $k_\mathrm{TST}$. The first term at right-hand-side is again the transmission function, which corrects upon $k_\mathrm{TST}$ for long enough correlation times $t$, by giving the fraction of transitions that end up in the product state, irrespective of the crossing direction, relative to the fraction that crosses in the forward direction. Or using Chandlers notation[10]:

$$k_R = \kappa \frac{1}{2}\left< \left\vert \dot{x} \right\ver... ...-\infty}^{x^\ddagger} \mathrm{d}x \exp{[-\beta \Delta G(x)]} }$$ (20)

The separated terms can in principle all be computed using molecular dynamics. The transmission function $\kappa$, expressed by the time correlation function in the first term of eqation 2.18, is usually obtained by performing a dynamics run with the reaction coordinate constrained at the transition state $x^\ddagger$ to obtain a large number of starting configurations with a positive crossing velocity, which are then used to initiate unconstrained MD runs to evaluate $\Theta(x(t)-x^\ddagger)$. The average crossing velocity (i.e. the second term) can be obtained from the initial constrained run. The last term, the probability to reach the transition state from the reactant state, can be obtained using a method to calculate free energy differences, such as umbrella sampling[13,14,15] or thermodynamic integration in combination with the method of constraint[16]. In the umbrella sampling method, the system is biased by addition of an (umbrella) potential $U(x)$, which approximately cancels the free energy barrier. In the MD run, the system can now move barrier-free back and forth from the reactant state to the product state. Accumulating the probability distribution $P(x)$ gives the desired free energy profile, $\Delta G(x) = -U(x) + k_BT \ln{P(x)}$. Using the method of constraint, one performs a number of constrained MD runs at different constraint values of the reaction coordinate. In each constrained run, one accumulates the force required to keep the reaction coordinate fixed. The obtained average constraint forces can be related to the thermodynamic force, $\mathrm{d}G(x)/ \mathrm{d}x$, along the reaction coordinate[17], so that integration gives the desired free energy profile.

Summarizing, we see that chemical reactions are activated processes, and in particular that barrier crossings are, on the time scale of the thermal motions, rare but very important events. Performing the reaction in solution affects the reaction rate constant in two ways: 1) the static energy barrier is changed, due to different energies of solvation for the reactant state and the transition state (on the same ground the equilibrium constant $K=k_R/k_P$ can be changed by different energies of solvation for the reactant state and the product state) and 2) the solvent dynamics can alter the rate by providing an extra energy sink after barrier crossing (increasing the rate constant) or by forcing the reactants into recrossings at the barrier top region (decreasing the rate constant). By separating the rate constant in terms of the probability to reach the barrier top, the crossing velocity and the transmission coefficient, solvent effects on the reaction rate constant can be studied with molecular dynamics techniques.