The Standard Potential Energy Surface

 The parameters for the standard PES employed in the present study are

$$\begin{eqnarray} \kappa&=&101 \ \ \,\qquad\qquad\qquad\sigma=1\\ E_2 &=&3.79688 ... ...ad\qquad l_{max}=10 \\ \Delta\nu&=&0.01 \ \qquad\qquad\qquad b=0.5\end{eqnarray}$$

(5.9)

(5.10)

(5.11)

(5.12)

These values produce a total of $1.11\times 10^{9}$ minima on the PES. Some of the thermodynamic properties of the system are shown in Figure 5.2. These were calculated using the following definitions: the microcanonical temperature, T, is given by $1/kT=\left(\partial\ln\Omega /\partial E\right)_{N,V}$,the canonical internal energy by $U=-\left(\partial\ln Z/\partial\beta\right)_{N,V}$,and the canonical heat capacity by $C_v=\left(\partial U/\partial T\right)_{N,V}$.In each case analytic derivatives of the appropriate partition functions were employed.

  

Figure 5.2: Thermodynamics of our standard PES: (a) caloric curves for the microcanonical (solid line) and canonical (dashed line) ensembles, (b) canonical heat capacity, and (c) equilibrium probabilities of occupation of level l, P_l^eq, for the microcanonical ensemble.

It can be seen from Figure 5.2 that the system shows the finite-size analogue of a first-order phase transition; there are prominent features in the caloric curves and a large peak in the heat capacity. The transition has a two-state character involving the global minimum and minima in the highest energy region of the PES, with intermediate minima never being significantly populated at equilibrium (Figure 5.2c). The energy at which the global minimum has an equilibrium probability of 0.5, E_f, can be used as a definition of the transition energy. For our standard' PES, E_f=38.2745.

Also of note is the Van der Waals loop that occurs in the microcanonical caloric curve (Figure 5.2a), which is a feature unique to finite systems[134]. Such loops have been commonly observed for clusters[135,136] and have also been noted in protein simulations[258]. In small clusters the physical cause of this region of negative heat capacity is the absence of phase separation, which is prevented by the relatively high energetic cost of an interface between the solid-like and liquid-like regions[176].

  

Figure 5.3: (a) Variation of the folding time, t_f, with energy for our standard PES. t_f is defined as the time for P₁ to reach 0.8, from an initial state P₁₀(0)=1. The vertical dashed line denotes the energy at which P₁^eq=0.8 and so is the upper limit for which t_f is defined. (Above this energy P₁(t) is always less than 0.8). (b) The variation of the probability of being in the global minimum, P₁, as a function of time and energy. The ten lines shown are for t=100 to 1000 in intervals of 100 and are labelled by the appropriate value of t.

The microcanonical dynamics of the system were examined by initially populating the highest energy state, i.e. setting P₁₀(0)=1, and observing the relaxation towards equilibrium at constant energy. The corresponding canonical results are omitted for brevity. We consider in particular the time required for the global minimum to develop an 80% probability of occupation. By analogy to relaxation in proteins we call this the `folding' time, t<sub>f</sub>; i.e. P<sub>1</sub>(t<sub>f</sub>)=0.8. The dependence of t<sub>f</sub> on E for the standard PES is shown in Figure 5.3a. There is a clear minimum in t<sub>f</sub> which occurs at an energy, E<sub>max</sub>, of 26.41 (=0.69E<sub>f</sub>). The folding time increases rapidly at low energy because it becomes more difficult to overcome the barriers on the PES, and increases at higher energies because the thermodynamic driving force towards the global minimum is diminished. Similar maxima have been seen in the folding rates of model proteins[247,259,260,261,262] in the `crystallization' rate of a (KCl)₃₂ cluster[91], and in the nucleation rate of glasses[263]; these maxima occur for the same reasons as for our model PES. Maxima in the folding rates have also been observed in experiments on proteins which are rapid folders[264,265,266], but the interpretation of these results is complicated by the fact that real proteins can undergo low temperature denaturation[267].

Another view of the relaxation dynamics is provided by Figure 5.3b which shows the increase in the population of the ground state with time; there is again an energy window for which the global minimum can become significantly populated. It is interesting to note that the energy at which the maximum in the population of the global minimum occurs decreases somewhat as the time increases. This indicates that to optimize the population of the global minimum one should not quench the system to a low energy and let it evolve, but rather use a cooling schedule whereby the energy is decreased as a function of time, as in constant thermodynamic speed annealing[268].

To further understand the dynamics of the system it is helpful to use Landau functions. In the microcanonical ensemble the presence of a first-order phase transition is indicated by two maxima in the Landau entropy separated by a well. We choose to call the minimum in the Landau entropy an entropy bottleneck, since it represents a constricted region of phase space which has a lower density of states than the regions that it connects. The Landau profiles also enable us to understand the influence of the thermodynamics on the dynamics. If there is one stable state, relaxation to that state is likely to be rapid unless other kinetic factors intervene. If there are two stable states, the rate of transition from the metastable to the stable state decreases as the entropy bottleneck becomes narrower or the free energy barrier increases. The importance of the Landau free energy in understanding protein folding has been recognized by Bryngelson , who considered the dynamics expected for a number of different Landau free energy profiles[245].

  

Figure 5.4: (a) Landau entropy profiles for our standard PES using the order parameter l. The seven lines shown are for E=20 to 80 in intervals of 10 and are labelled by the appropriate value of E. (b) Depth of the Landau entropy bottleneck relative to the Landau entropy maximum at large values of l (`folding') and relative to the Landau entropy maximum at the global minimum (`unfolding'). The dashed line gives the magnitude of the Landau entropy difference between the two stable states.

Application of a Landau-type analysis using a restricted thermodynamic potential depends on finding a suitable order parameter which can differentiate between the states of interest. For many systems this can be a difficult task, but for our model PES l provides a natural order parameter since it is a measure of the distance from the global minimum. The Landau entropy is then defined as $S_L(l)=S(E)+k\ln P_l^{eq}(E)$ and the Landau free energy by $A_L(l)=A(T)-kT\ln P_l^{eq}(T)$. Landau entropy profiles for our standard PES are given in Figure 5.4 for a range of energies. At low energy there is a single maximum corresponding to the global minimum, and relaxation to the global minimum is thermodynamically `downhill'. Consequently at E<sub>max</sub>, the probability density `flows' smoothly down the PES with the intermediate states having significant transient populations (Figure 5.5).

The equilibration trees in Figure 5.5 map out the times at which different levels come in to approximate equilibrium with each other (as defined by equation 5.3). At E_max levels 1 and 2 first reach local equilibrium (Figure 5.5e), i.e. equilibrium first occurs between the levels corresponding to the Landau entropy maximum. Next, level 3 reaches equilibrium with levels 1 and 2, and so on to higher values of l: the local equilibrium propagates to higher energy away from the Landau entropy maximum. Similar behaviour occurs at lower energies, except that the barriers on the PES have a larger retarding effect on the dynamics. Therefore, the equilibration tree at E=15 (Figure 5.5c) shows the same structure as at E_max except that the time scales are around a thousand times slower. In the microcanonical ensemble the kinetic energy increases as the potential energy decreases. Consequently, at low energies the transition rate over the barriers is much slower at the top of the PES, and the system accelerates as it descends the PES. Therefore, transient populations in the intermediate states are less pronounced at E=15 than at E_max. In the canonical ensemble this effect is absent because the kinetic energy is independent of the position in configuration space.

  

Figure 5.5: (a)-(c) Time evolution of the probabilities, P_l, and (d)-(f) equilibration trees for our standard PES at an energy of (a) and (d) 15; (b) and (e) 26.41; and (c) and (f) 36. Initially, P₁₀=1. The energy 26.41 corresponds to the minimum in the folding' time of Figure 5.3a. Each line in the equilibration tree represents a single level or a set of levels that are in equilibrium with each other. The lines join at times when states first come into local equilibrium as defined by equation 5.33. The vertical position of these nodes is given by the average value of l for the levels that are in equilibrium with each other. The horizontal positions of the ends of lines on the left-hand side of the graph, which correspond to individual levels, are arbitrary.

For E>28.238, there are two Landau entropy maxima, which correspond to the global minimum and the high potential energy `unfolded' states. The entropy bottleneck occurs because at lower potential energy the loss in entropy due to the reduced number of minima is not completely compensated by the greater vibrational entropy that results from the larger kinetic energy. As a result of the entropy bottleneck, the equilibration trees have a qualitatively different form from those at lower energies; equilibration initially occurs within the secondary Landau entropy maximum corresponding to the higher lying levels and then propagates down the PES (Figure 5.5f). It is interesting to note that this change in direction of equilibration appears at approximately the same energy as the Landau entropy bottleneck. For E=36 there is a separation of time scales between equilibration within the secondary Landau entropy maximum, and equilibration between the maxima. This is reflected in the time evolution of the probability distribution: the probabilities of occupying higher potential energy levels rapidly attain approximate local equilibrium values which then slowly decay as the global minimum becomes more populated (Figure 5.5c).

Similar behaviour is seen for the Landau free energy in the canonical ensemble. In fact, the temperature at which A_L develops a double well can be derived analytically (assuming that $\Delta \nu$ is small). Two free energy minima are observed if  

$$kT\gt{\Delta E \over \ln g + \kappa \Delta \nu}.$$ (5.13)

This corresponds to kT>0.242 for our standard PES. As relaxation to the global minimum is impeded by a free energy barrier, equation 5.13 implies that large $\Delta E$ and small g and $\Delta \nu$ assist relaxation to the global minimum. These are the same conditions as those that lead to a large value of T_f.

Landau free energy profiles have previously been calculated for lattice models of proteins[233,234,247,250,258]. They show a similar behaviour to that found for our model PES, with a single minimum at low temperature, and two stable states at higher temperatures corresponding to the folded and unfolded protein[233,234,247]. For clusters there can be a range of temperature around the melting transition for which Landau functions reveal two stable states corresponding to the solid-like and liquid-like states[112,133]. Similarly it has been shown that the nucleation of a crystalline phase from the bulk Lennard-Jones liquid involves a large free energy barrier[269,270].