Fundamentals

Consider a unimolecular isomerization reaction:

$$A ^{\stackrel{\textstyle k_{AB}}{\textstyle \rightharpoonup}}_{\stackrel{\textstyle\leftharpoondown}{\textstyle k_{BA}}} B$$ (246)

$k_{AB}$ is the forward rate constant for the reaction, and $k_{BA}$ is the backward rate constant. They are required to express the concentration balances for species $A$ and $B$ in the traditional way:

$$\frac{d C_A}{d t} = -k_{AB}C_A + k_{BA}C_B$$ (247)

$$\frac{d C_B}{d t} = k_{AB}C_A - k_{BA}C_B$$ (248)

At equilibrium, the time derivatives vanish, and $C_A \rightarrow \left<C_A\right>$, and

$$K = \frac{\left<C_A\right>}{\left<C_B\right>} = \frac{k_{BA}}{k_{AB}}$$ (249)

We now ask, how does an initial perturbation in the concentration of $A$, written $\Delta A$, decay with time, assuming $C_A + C_B = C$ = const.?

$$\frac{d \left(C_A + \Delta C_A\right)}{d t} = -k_{AB}\left(C_A + \Delta C_A\right) + k_{BA}\left(C_B - \Delta C_A\right)$$ (250)

$$\frac{d C_A}{d t} = -k_{AB}\Delta C_A - k_{BA}\Delta C_A = -\left(k_{AB}+k_{BA}\right)\Delta C_A$$ (251)

We can easily solve this to yield:

$$\Delta C_A\left(t\right) = \Delta C_A\left(0\right)\exp\left... ...A}\right)t\right] \equiv \Delta C_A\left(0\right)e^{-t/\tau_R}$$ (252)

which defines the ``decay time'':

$$\tau_R = \left(k_{AB} + k_{BA}\right)^{-1}$$ (253)

$$= k_{AB}^{-1}\left(1+\frac{\left<C_A\right>}{\left<C_B\right>}\right)^{-1}$$ (254)

$$= \frac{\left<C_B\right>}{k_{AB}}$$ (255)

(The last equality assumes the concentrations are normalized such that $C_A + C_B = 1$; that is, we can consider $C_A$ as the probability of observing ``state'' $A$.)

Now we make a connection to the microscopic. Imagine that the two states labelled $A$ and $B$ are separated from each other along some general reaction coordinate $q$ which has a large free energy barrier at position $q^*$:

This is a free energy barrier because the system presumably has many more degrees of freedom other than (which may or may not contribute to) $q$. Now, we enter the realm of linear response theory (See Appendix C in F&S), and ask, what is the behavior or the system with a finitely small, static perturbation toward state $A$? We express this as a perturbation Hamiltonian:

$$\mathscr{H} = \mathscr{H}_0 -\epsilon g_A\left(q-q^*\right)$$ (256)

where $\mathscr{H}_0$ is the reference Hamiltonian of the unperturbed system. Here, $g_A$ is an indicator function which is 1 if we are in state $A$ and 0 otherwise:

$$g_A\left(x\right) = \left\{\begin{array}{ll} 1 & \mbox{$x<0$}\\ 0 & \mbox{$x>0$} \end{array}\right.$$ (257)

Notice that the peturbation lowers the energy by a little bit ($\epsilon$) when the system chooses to be in state $A$. Here, we think of $\epsilon$ as a switch we can flip in order to perturb the system. We measure the static pertubation as the difference between the average concentration of $A$ in the perturbed state to that in the unperturbed state:

$$\Delta C_A = \left<C_A\right>_\epsilon - \left<C_A\right>_0 = \left<g_A\right>_\epsilon - \left<g_A\right>_0$$ (258)

Notice that because of our choice in magnitude for $g_A$, we can consider its ensemble average, $\left<g_A\right>$, to be the probability of observing the system in state A, which is, due to our normalization of concentration, equivalent to the concentration of A.

Our first step is to find the linear response of $\Delta C_A$ to $\epsilon$, defined as

$$\left(\frac{\partial\Delta C_A}{\partial\epsilon}\right)_{\e... ...partial\left<g_A\right>}{\partial\epsilon}\right)_{\epsilon=0}$$ (259)

So let us first compute $\left<g_A\right>_\epsilon$ in the canonical ensemble:

$$\left<g_A\right>_\epsilon = \frac{ \mbox{\begin{minipage}{9c... ...\right)\right]}_{Z_\epsilon} \end{displaymath}\end{minipage}}}$$ (260)

Here, we have defined $Z_\epsilon$ as the unnormalized canonical partition function, merely to make the following gymnastics a little more transparent. Differentiating with respect to $\epsilon$:

$$\frac{\partial\left<g_A\right>_\epsilon}{\partial\epsilon} = \frac{1}{Z_\epsilon}\int d{\bf q}^Nd{\bf p}^N \exp\left(-\beta\mathscr{H}_0+\beta\epsilon g_A\right)\beta g_A^2$$ (261)

$$- \frac{1}{Z_\epsilon^2}\left[\int d{\bf q}^Nd{\bf p}^N \exp\left(-\beta\mathscr{H}_0+\beta\epsilon g_A\right)g_A\right]$$

$$\times \left[\int d{\bf q}^Nd{\bf p}^N \exp\left(-\beta\mathscr{H}_0+\beta\epsilon g_A\right)\beta g_A\right]$$

$$= \beta\left<g_A^2\right>_\epsilon - \beta\left<g_A\right>^2_\epsilon$$ (262)

$$\rightarrow \left(\frac{\partial\Delta C_A}{\partial\epsilon}\right)_{\epsilon=0} = \beta\left(\left<g_A^2\right>_0 - \left<g_A\right>_0^2\right)$$ (263)

Now, we haven't said anything in detail about the structure of $g_A$, but in fact, to perform this analysis, we require that it be continuous through the origin. But, we can make it vary from 1 to 0 across as narrow a region of $q$ around the origin as we like. Then it is true that, at all values of $q$ except perhaps right at the barrier position $q^*$,

$$g_A^2 = g_A \mbox{so } \left<g_A^2\right> = \left<g_A\right>$$ (264)

This implies

$$\left(\frac{\partial\Delta C_A}{\partial\epsilon}\right)_{\e... ...ight>_0\right)\right] = \beta\left<C_A\right>\left<C_B\right>$$ (265)

The next step is to introduce a clock that begins ticking the moment we set $\epsilon = 0$. As long as $\epsilon$ was small, we can compute the decay of the static perturbation $\Delta C_A\left(t\right)$ to first order in $\epsilon$ as²

$$\Delta C_A\left(t\right) = \beta\epsilon \frac{ \mbox{ \begi... ... p}^Ne^{-\beta\mathscr{H}_0} \end{displaymath}\end{minipage}}}$$ (269)

$iL_0$ is the Liouville operator, which we first encountered in the Liouville operator formalism for deriving MD integrators (Sec. 4.1.2), and $e^{iL_0t}$ is the classical propagator:

$$e^{iL_0t}f\left(0\right) = f\left(t\right)$$ (270)

This yields:

$$\Delta C_A\left(t\right) = \beta\epsilon\left<\Delta g_A\lef... ...right)\right>, \Delta g_A \equiv g_A - \left<g_A\right>$$ (271)

Now, Eq. 266, when integrated implies

$$\Delta C_A\left(0\right) = \beta\epsilon\left<C_A\right>\left<C_B\right>$$ (272)

which we can use to eliminate $\epsilon$ in Eq. 272, yielding

$$\Delta C_A\left(t\right) = \Delta C_A\left(0\right) \frac{\l... ...ta g_A\left(t\right)\right>}{\left<C_A\right>\left<C_B\right>}$$ (273)

Now, recalling Eq. 253, we can make a connection to the macroscopic:

$$e^{-t/\tau_R} = \frac{\left<\Delta g_A\left(0\right)\Delta g_A\left(t\right)\right>}{\left<C_A\right>\left<C_B\right>}$$ (274)

This result states that the decay behavior is governed by the autocorrelation function of concentration fluctuations. Eq. 275 is itself a remarkable result of modern nonequilibrium statistical mechanics, attributable to Onsager.

We have implicitly assumed that the system spends all of its time in either $A$ or $B$; the system spends virtually no time at the barrier crossing itself. Eq. 275 is therefore valid as long as the decay time, $\tau_R$, is much greater than the ``barrier crossing time''; usually, this is true (we assume).

Now, we differentiate Eq. 275 with respect to time:

$$-\frac{1}{\tau_R}e^{-t/\tau_R} = \frac{\left<g_A\left(0\righ... ...t{g}_A\left(t\right)\right>}{\left<C_A\right>\left<C_B\right>}$$ (275)

where $\dot{g_A} \equiv d g_A / dt$. Note that the $\Delta$'s are gone thanks to the fact that $d\left<g_A\right>/dt = 0$.

For time much less than $\tau_R$, but still much greater than the barrier crossing time, $(t \ll \tau_R)$,³

$$-\left(\tau_R\right)^{-1} = \frac{\left<g_A\left(0\right)\do... ...ht)g_A\left(t\right)\right>}{\left<C_A\right>\left<C_B\right>}$$ (278)

Recalling from way back from Eq. 256 how $\tau_R$ and the rate constant $k_{AB}$ are related, we find finally that

$$k_{AB}\left(t\right) = \frac{\left<\dot{g}_A\left(0\right)g_A\left(t\right)\right>}{\left<C_A\right>}$$ (279)

So: A rate constant is the time-derivative of a concentration autocorrelation function.

The idea of path sampling Monte Carlo is that, because $k_{AB}$ is an ensemble average, we can compute it using MC, if we determine first what ensemble we need to sample.