Methods

  A master equation can be used to describe the dynamics on a multi-dimensional PES[254,255], thus giving the time evolution of the probability distribution, ${\mathbf{P}}(t)=\{P_i(t)\}$,where P_i is the probability of the system being in state i. The master equation has the form  

$${dP_i\over dt}=\sum_j w_{ij} P_j,$$ (5.1)

where $w_{ij}=W_{ij}-\delta_{ij}\sum_{k} W_{ki}$ and W_ij is the rate of transition from state j to state i. The diagonal elements of the matrix, w_ii, give the total rate of transition out of state i.

If the transition matrix w is neither decomposable nor splitting, the matrix has a single zero eigenvalue[255], whose eigenvector corresponds to the equilibrium probability distribution, ${\mathbf{P}}^{eq}$. The master equation describes the evolution of a system from an initial distribution, ${\mathbf{P}}(0)$,towards ${\mathbf{P}}^{eq}$;at infinite time the probability distribution must be equal to ${\mathbf{P}}^{eq}$.The dynamics and thermodynamics must be consistent in this limit.

The transition matrix must satisfy detailed balance to be physically reasonable, i.e. W_ijP^eq_j=W_jiP^eq_i. As a consequence, the solution of the master equation can be expanded in a complete set of eigenfunctions of the symmetric matrix, $\mathbf{w'}$, defined as $w'_{ij}=\sqrt{P^{eq}_j/P^{eq}_i} w_{ij}$. The result is[254]

$$P_i(t)=\sqrt{P^{eq}_i}\sum_{j,k} u_i^j e^{\lambda_j t} u_k^j {P_k(0)\over \sqrt{P^{eq}_k}},$$ (5.2)

where u_i^j is the $i^{\mathrm{th}}$ element of the $j^{\mathrm{th}}$ eigenvector of $\mathbf{w'}$ and $\lambda_j$ is the $j^{\mathrm{th}}$ eigenvalue of $\mathbf{w'}$. Apart from the zero eigenvalue, the $\lambda_j$ are all negative. In all our calculations, we have used the above analytical solution. An alternative method of solving the master equation is to integrate equation 5.1 numerically.

The concept of detailed balance also allows us to define when two states come into local equilibrium; i.e. $W_{ij}P_j(t)\approx W_{ji}P_i(t)$. The precise condition that we use is  

$${\left\vert P_i(t) P_j^{eq} - P_j(t) P_i^{eq}\right\vert\over \sqrt{P_i(t) P_j(t) P_i^{eq} P_j^{eq}}}\leq 0.01,$$ (5.3)

i.e. the two states are within 1% of equilibrium. The equilibration time for two states is then defined as the time at which the above inequality is first satisfied. This allows us to construct equilibration trees mapping out how the system proceeds to the state where the whole system is in equilibrium[256,257].

  

Figure 5.1: (a) One-dimensional cross-section through our standard PES showing a reaction pathway from the top of the PES to the global minimum. Integer values of l correspond to minima and half-integer values to transition states. (b) Schematic depiction of a PES with g=2 and l_max=6 showing the dramatic increase in the number of minima with energy.

The standard PES that we consider is depicted in Figure 5.1. It consists of a single funnel in which the number of minima increases rapidly with energy (Figure 5.1b). To simplify the calculations the minima have been grouped into l_max levels. The minima in each level are assumed to have identical properties, and to be always in equilibrium with each other, thus allowing us to consider each level as a single state in the master equation. This framework is equivalent to assuming that the barriers between minima in the same level are zero. Level l=1 is the global minimum, and the number of minima increases geometrically as the PES is ascended. There are g times more minima in each subsequent level. Therefore the number of minima in level l, n_l, is g^l-1 and the total number of minima on the PES is (g^l_max-1)/(g-1). We do not consider permutational isomers explicitly since they do not affect the relaxation dynamics. For a real PES, there are $\mathcal{O}(N!)$ permutational isomers of each minimum. Therefore, there will be many funnels leading down to the different permutational isomers of the global minimum. However, as each funnel is identical they do not need to be considered separately.

We assume that the minima in level l are only connected to minima in levels $l\pm1$. The transition states connecting minima in adjacent levels are all assumed to lie an energy b above the higher minimum. Each minimum in level l is assumed to be connected to $\sigma$ minima in level l-1 and hence to $g\sigma$ minima in level l+1. From this information, and using RRKM theory[166] within the harmonic approximation, the microcanonical rates for transitions between levels are

$$\begin{eqnarray} W_{l+1,l}&=&g \sigma \left( {E-E_{l+1}-b\over E-E_l}\right)^{\k... ...over{\overline\nu_{l-1,l}}^{\kappa-1}} \quad\quad\hbox{downhill},\end{eqnarray}$$

(5.4)

where E_l is the energy of level l, $\kappa$ is the number of internal degrees of freedom, and $\overline\nu_l$ and $\overline\nu_{l,l+1}$ are the geometric mean vibrational frequencies of the minima in level l, and the transition states between levels l and l+1, respectively. All other elements of the matrix W are assumed to be 0 in the present treatment. Similarly, the canonical transition rates are given by

$$\begin{eqnarray} W_{l+1,l}&=&g \sigma {{\overline\nu_l}^\kappa\over{\overline\nu... ...exp\left(-\beta b)\right) \qquad\qquad\qquad\quad\hbox{downhill},\end{eqnarray}$$

(5.5)

where $\beta$ is the inverse temperature. We have performed calculations both in the microcanonical and canonical ensembles, but since most of the properties probed show a weak ensemble dependence the majority of the results reported here are microcanonical. The microcanonical ensemble is appropriate to describe an isolated cluster in vacuo, but the canonical ensemble is probably more appropriate for proteins where the solvent can act as a heat bath.

The thermodynamics of the system can be described using a superposition method (Chapter 3), whereby the total energy density of states, $\Omega(E)$, is constructed by summing the density of states for all the energetically accessible minima on the PES. Applying this method to the model PES within the harmonic approximation gives

$$\Omega(E)=\sum_{l=1}^{l_{max}}{n_l (E-E_l)^{\kappa-1} \over\Gamma(\kappa){\overline\nu_l}^\kappa},$$ (5.6)

where $\Gamma$ is the Gamma function. It follows that the microcanonical equilibrium probabilities are given by

$$P_l^{eq}(E) ={n_l (E-E_l)^{\kappa-1}\over \Gamma(\kappa){\overline\nu_l}^\kappa\Omega(E)}.$$ (5.7)

The canonical partition function and equilibrium probabilities can be constructed in a similar manner:

$$Z(T)=\sum_{l=1}^{l_{max}} {n_l e^{-\beta E_l}\over (\beta h... ...T)={n_l e^{-\beta E_l}\over (\beta h \overline\nu_l)^\kappa Z}.$$ (5.8)

The global minimum defines the energy zero, i.e. E₁=0. Beyond l=2 the potential energy is assumed to increase linearly with l. Therefore, $E_l=E_2+(l-2)\Delta E$ for $l\ge 2$, where $\Delta E$ is a measure of the potential energy gradient of the funnel. For a typical cluster PES, the mean vibrational frequency is smaller for minima of higher potential energy[153]--the stabilization of the liquid-like phase at high temperatures is due to both the large number of minima and the greater vibrational entropy. Therefore, we use $\overline\nu_l=1-(l-1)\Delta\nu.$This defines the unit of time as the vibrational period of the global minimum. The mean vibrational frequency of a transition state is assumed to be the geometric mean of the vibrational frequencies of the two minima it connects, i.e. $\overline\nu_{l+1,l}=\sqrt{\nu_l \nu_{l+1}}$