# The Dynamics of Structural Transitions in Sodium Chloride Clusters

## I Introduction

Understanding the relationship between the potential energy surface (PES), or energy landscape, and the dynamics of a complex system is a major research effort in the chemical physics community. One particular focus has been the dynamics of relaxation from an out-of equilibrium starting configuration down the energy landscape to the state of lowest free energy, which is often also the global minimum of the PES.

The possible relaxation processes involved can be roughly divided into two types. The first is relaxation from a high-energy disordered state to a low-energy ordered state, and examples include the folding of a protein from the denatured state to the native structure, and the formation of a crystal from the liquid.

The second kind of relaxation process, which is the focus of the work here, is relaxation from a low-energy, but metastable, state to the most stable state. This second situation often arises from the first; the initial relaxation from a disordered state can lead to the population of a number of low-energy kinetically accessible configurations. The time scales for this second relaxation process can be particularly long because of the large free energy barriers that can separate the states.

Some proteins provide an instance of this second type of relaxation. Often, as well as a rapid direct path from the denatured state to the native structure, there is also a slower path which passes through a low-energy kinetic trap. As this trapping is a potential problem for protein function, the cell has developed its own biochemical machinery to circumvent it. For example, it has been suggested that the chaperonin, GroEL, aids protein folding by unfolding those protein molecules which get stuck in a trapped state.

There are also a growing number of examples of this second type of relaxation involving clusters. For Lennard-Jones (LJ) clusters there are a number of sizes for which the global minimum is non-icosahedral. For example, for LJ<sub>38</sub> the global minimum is a face-centred-cubic truncated octahedron, but relaxation down the PES almost always leads to a low-energy icosahedral minimum. Similarly, for LJ<sub>75</sub> the icosahedral states act as a trap preventing relaxation to the Marks decahedral global minimum. A similar competition between face-centred-cubic or decahedral and icosahedral structures has recently been observed for metal clusters.

For these clusters unbiased global optimization is difficult because the icosahedral states act as an effective trap. More generally, kinetic traps are one of the major problems for a global optimization algorithm. Therefore, much research has focussed on decreasing the ‘life-time’ of such traps. For example, some methods simulate a non-Boltzmann ensemble that involves increased fluctuations, thus making barrier crossing more likely. Other algorithms transform the energy landscape in a way that increases the temperature range where the global minimum is populated, thus allowing one to choose conditions where the free energy barriers relative to the thermal energy are lower.

Recently, clear examples of trapping associated with structural transitions have emerged from experiments on NaCl clusters. These clusters have only one energetically favourable morphology: the magic numbers that appear in mass spectra correspond to cuboidal fragments of the bulk crystal (rocksalt) lattice, hence the term nanocrystals. Indirect structural information comes from the experiments of Jarrold and coworkers which probe the mobility of size-selected cluster ions. For most (NaCl)<sub>N</sub>Cl<sup>-</sup> with $`N>30`$, multiple isomers were detected which were assigned as nanocrystals with different cuboidal shapes. The populations in the different isomers were not initially equilibrated, but slowly evolved, allowing rates and activation energies for the structural transitions between the nanocrystals to be obtained.

In a previous paper we identified the mechanisms of these structural transitions by extensively searching the low-energy regions of the PES of one of these clusters, (NaCl)<sub>35</sub>Cl<sup>-</sup>, in order to obtain paths linking the different cuboidal morphologies. The key process in these transitions is a highly cooperative rearrangement in which two parts of the nanocrystal slip past one another on a $`\{110\}`$ plane in a $``$11̄0$``$ direction.

Here we continue our examination of the structural transitions by investigating the dynamics of (NaCl)<sub>35</sub>Cl<sup>-</sup>. Given the long time scales for which the clusters reside in metastable forms, it is not feasible to probe the transitions with conventional dynamics simulations. Instead, we use a master equation that describes the probability flow between the minima on the PES. This method has the advantage that we can relate the dynamics to the topography of the PES . In this paper we are particularly concerned with obtaining activation energies for the structural transitions: firstly, in order to compare with experiment, and secondly, to understand how the activation energy for a process that involves a series of rearrangements and a large number of possible paths is related to the features of the energy landscape. In section II we outline our methods, and then in Section III, after a brief examination of the topography of the PES and the thermodynamics, we present our results for the dynamics of the structural transitions.

## II Methods

### A Potential

The potential that we use to describe the sodium chloride clusters is the Tosi-Fumi parameterization of the Coulomb plus Born-Mayer form:

$$E=\underset{i<j}{}\left(\frac{q_iq_j}{r_{ij}}+A_{ij}e^{r_{ij}/\rho }\right),$$

where $`q_i`$ is the charge on ion $`i`$, $`r_{ij}`$ is the distance between ions $`i`$ and $`j`$ and $`A_{ij}`$ and $`\rho `$ are parameters. Although simple, this potential provides a reasonable description of the interactions. For example, in a previous study we compared the global minima of (NaCl)<sub>N</sub>Cl<sup>-</sup> ($`N35`$) for this potential with those for a more complex potential derived by Welch et al. which also includes terms due to polarization. Most of the global minima were the same for both potentials. Given some of the other approximations we use in this study, the small advantages gained by using the Welch potential do not warrant the considerable additional computational expense.

We should also note that a well-known problem associated with the above family of potentials for the alkali halides is that they never predict the CsCl structure to be the most stable. This problem arises because the potentials do not allow the properties of an ion to be dependent on the local ionic environment. This deficiency should not greatly affect the relative energies of the low-lying (NaCl)<sub>35</sub>Cl<sup>-</sup> minima because they all have the same rock-salt structure, but it may affect the barriers for rearrangements where some ions experience a different local environment at the transition state.

### B Searching the potential energy surface

The samples of 3518 minima and 4893 transition states that we use here were obtained in our previous study on the mechanisms of the structural transitions for (NaCl)<sub>35</sub>Cl<sup>-</sup>. This sampling was performed by repeatedly stepping across the PES from minimum to minimum via transition states, thus giving a connected set of stationary points. We biased this search to probe the low-energy regions of the PES either by using a Metropolis criterion to decide whether to accept a step to a new minimum or by systematically performing transition state searches from the lower energy minima in the sample. Thus, although our samples of stationary points only constitute a tiny fraction of the total number on the PES, we have a good representation of the low-energy regions that are relevant to the structural transitions of (NaCl)<sub>35</sub>Cl<sup>-</sup>.

### C Thermodynamics

We used two methods to probe the thermodynamics: first, conventional Monte Carlo simulations and second, the superposition method. The latter is a technique to obtain the thermodynamics of a system from a sample of minima. It is based on the idea that all of configuration space can be divided up into the basins of attraction surrounding each minimum. The density of states or partition function can then be written as a sum over all the minima on the PES, e.g. $`Z=_iZ_i`$, where $`Z_i`$ is the partition function of minimum $`i`$.

The limitations of the superposition method are that the $`Z_i`$ are not known exactly and that, for all but the smallest systems, the total number of minima on the PES is too large for us to characterize them all. However, the harmonic expression for $`Z_i`$ leads to a reasonable description of the thermodynamics. Furthermore, anharmonic forms are available which allow the thermodynamics of larger clusters to be reproduced accurately.

The incompleteness of the sample can be overcome by weighting the contributions from the minima in a representative sample. However, this approach is not necessary in the present study since we are interested in low temperature behaviour where the number of thermodynamically relevant minima is still relatively small. Furthermore, in this temperature regime the superposition method has the advantage that it is unaffected by large free energy barriers between low-energy minima which can hinder the determination of equilibrium thermodynamic properties by conventional simulation.

Here we use the harmonic form of the superposition method, because we later use the harmonic approximation to derive rate constants (reliable anharmonic expressions for the rate constants are not so readily available). The partition function is then

$$Z=\underset{i}{}\frac{n_ie^{\beta E_i}}{(\beta h\overline{\nu }_i)^\kappa },$$
(1)   
where $`\beta =1/kT`$, $`E_i`$ is the energy of minimum $`i`$, $`\overline{\nu }_i`$ is the geometric mean vibrational frequency of $`i`$, $`\kappa =3N6`$ is the number of vibrational degrees of freedom and $`n_i`$ is the number of permutational isomers of $`i`$. $`n_i`$ is given by $`2N!/h_i`$, where $`h_i`$ is the order of the point group of $`i`$. From this expression thermodynamic quantities such as the heat capacity, $`C_v`$, can be obtained by the application of standard thermodynamic formulae. The superposition method also allows us to examine the contributions of particular regions of configuration space to the the thermodynamics. For example, the probability that the system is in region $`A`$ is given by

$$p_A=\underset{iA}{}\frac{Z_i}{Z},$$
(2)   
where the sum runs over those minima that are part of region A.

### D Dynamics

The master equation approach is increasingly being used to describe the inter-minimum dynamics on a multi-dimensional PES with applications to, for example, clusters, glasses, proteins and idealized model landscapes. The master equation is defined in terms of $`𝐏(𝐭)=\{P_i(t)\}`$, the vector whose components are the ensemble-average probabilities that the system is associated with each minimum at time $`t`$:

$$\frac{dP_i(t)}{dt}=\underset{ji}{\overset{n_{\mathrm{min}}}{}}[k_{ij}P_j(t)k_{ji}P_i(t)],$$
(3)   
where $`k_{ij}`$ is the rate constant for transitions from minimum $`j`$ to minimum $`i`$. Defining the matrix

$$W_{ij}=k_{ij}\delta _{ij}\underset{m=1}{\overset{n_{\mathrm{min}}}{}}k_{mi}$$
(4)   
allows Equation (3) to be written in matrix form: $`d𝐏(t)/dt=\mathrm{𝐖𝐏}(t)`$.

If the transition matrix $`𝐖`$ cannot be decomposed into block form, it has a single zero eigenvalue whose corresponding eigenvector is the equilibrium probability distribution, $`𝐏^{eq}`$. As a physically reasonable definition for the rate constants must obey detailed balance at equilibrium, i.e. $`W_{ij}P_j^{\mathrm{eq}}=W_{ji}P_i^{\mathrm{eq}}`$, the solution of the master equation can be expanded in terms of a complete set of eigenfunctions of the symmetric matrix, $`\stackrel{~}{𝐖}`$, defined by $`\stackrel{~}{W}_{ij}=(P_j^{\mathrm{eq}}/P_i^{\mathrm{eq}})^{1/2}W_{ij}`$. The solution is

$$P_i(t)=\sqrt{P_i^{\mathrm{eq}}}\underset{j=1}{\overset{n_{\mathrm{min}}}{}}\stackrel{~}{u}_i^{(j)}e^{\lambda _jt}\left[\underset{m=1}{\overset{n_{\mathrm{min}}}{}}\stackrel{~}{u}_m^{(j)}\frac{P_m(0)}{\sqrt{P_m^{\mathrm{eq}}}}\right],$$
(5)   
where $`\stackrel{~}{u}_i^{(j)}`$ is component $`i`$ of the $`j^{\mathrm{th}}`$ eigenvector of $`\stackrel{~}{𝐖}`$ and $`\lambda _j`$ is the $`j^{\mathrm{th}}`$ eigenvalue.

The eigenvalues of $`𝐖`$ and $`\stackrel{~}{𝐖}`$ are identical and the eigenvectors are related by $`u_i^{(j)}=\stackrel{~}{u}_i^{(j)}\sqrt{P_i^{eq}}`$. Except for the zero eigenvalue, all $`\lambda _j`$ are negative. Therefore, as $`t\mathrm{}`$ the contribution of these modes decays exponentially to zero and $`𝐏𝐏^{eq}`$.

To apply Equation (5) we must first diagonalize $`\stackrel{~}{𝐖}`$. The computer time required for this procedure scales as the cube of the size of the matrix and the memory requirements scale as the square. Therefore, it is advantageous for the matrix $`\stackrel{~}{𝐖}`$ to be as small as possible. For this reason we recursively removed those minima that are only connected to one other minimum; these ‘dead-end’ minima do not contribute directly to the probability flow between different regions of the PES. After pruning our samples have 1624 minima and 2639 transition states. To test the effect of this pruning, we performed some calculations using both the full and the pruned samples. The effect on the dynamics of the structural transitions was negligible.

As the temperature of a system is decreased the spread of eigenvalues can increase rapidly. When the ratio of the largest to smallest eigenvalues is of the order of the machine precision of the computer, the accuracy of the extreme eigenvalues can become degraded by rounding errors. We encountered these problems below 275K. Without pruning the samples these numerical difficulties are more pronounced.

We model the rate constants, which are needed as input to Equation (3), using RRKM theory in the harmonic approximation. Therefore, the rate constant for a transition from minimum $`i`$ to minimum $`j`$ via a particular transition state (denoted by $``$) is given by

$$k_{ij}^{}(T)=\frac{h_i}{h_{ij}^{}}\frac{\overline{\nu }_i^\kappa }{\overline{\nu }_{ij}^{(\kappa 1)}}\mathrm{exp}((E_{ij}^{}E_i)/kT).$$
(6)   
## III Results

### A Topography of the (NaCl)<sub>35</sub>Cl<sup>-</sup> PES

In our previous study of (NaCl)<sub>35</sub>Cl<sup>-</sup> we found that the low-energy minima all had rock-salt structures. The different minima have four basic shapes: an incomplete $`5\times 5\times 3`$ cuboid, a $`6\times 4\times 3`$ cuboid with a single vacancy, a $`8\times 3\times 3`$ cuboid with a single vacancy and an incomplete $`5\times 4\times 4`$ cuboid. The lowest-energy minimum for each of these forms is shown in Figure 1.

In the experiments on (NaCl)<sub>35</sub>Cl<sup>-</sup> by Jarrold and coworkers the three peaks that were resolved in the arrival time distribution were assigned on the basis of calculated mobilities as $`5\times 5\times 3`$, $`5\times 4\times 4`$ and $`8\times 3\times 3`$ nanocrystals. However, when the $`6\times 4\times 3`$ nanocrystal is also considered, better agreement between the calculated and observed mobilities can be obtained by assigning the three experimental peaks to the $`6\times 4\times 3`$, $`5\times 5\times 3`$ and $`8\times 3\times 3`$ nanocrystals in order of increasing drift time. This reassignment is also in better agreement with the thermodynamics since the clusters convert to (what is now assigned as) the $`5\times 5\times 3`$ nanocrystal as time progresses, indicating that this structure has the lowest free energy. In our calculations a $`5\times 5\times 3`$ isomer is the global potential energy minimum, and the $`5\times 5\times 3`$ nanocrystal is always more stable than the other nanocrystals (See Section III B).

Disconnectivity graphs provide a way of visualizing an energy landscape that is particularly useful for obtaining insight into dynamics, and have previously been applied to a number of protein models and clusters. The graphs are constructed by performing a ‘superbasin’ analysis at a series of energies. This analysis involves grouping minima into disjoint sets, called superbasins, whose members are connected by pathways that never exceed the specified energy. At each energy a superbasin is represented by a node, and lines join nodes in one level to their daughter nodes in the level below. Every line terminates at a local minimum. The graphs therefore present a visual representation of the hierarchy of barriers between minima.

The disconnectivity graph for (NaCl)<sub>35</sub>Cl<sup>-</sup> is shown in Figure 2. The barriers between minima with the same cuboidal form are generally lower than those between minima that have different shapes. Therefore, the disconnectivity graph splits into funnels corresponding to each cuboidal morphology. (A funnel is a set of downhill pathways that converge on a single low-energy minimum or a set of closely-related low-energy minima. In a disconnectivity graph an ideal funnel is represented by a single tall stem with lines branching directly from it, indicating the progressive exclusion of minima as the energy is decreased.) The separation is least clear for the $`5\times 4\times 4`$ minima because of the large number of different ways that the nine vacant sites can be arranged. For example, these vacancies are organized very differently in the two lowest-energy $`5\times 4\times 4`$ minima (Figure 1), and in fact the barrier between minimum O and the low-energy $`6\times 4\times 3`$ isomers is lower than the barrier between O and minimum L. Therefore, the minima associated with O form a sub-funnel that splits off from the $`6\times 4\times 3`$ funnel, rather than being directly connected to the main $`5\times 4\times 4`$ funnel.

The disconnectivity graph shows that the barriers between the $`5\times 5\times 3`$, $`6\times 4\times 3`$ and $`5\times 4\times 4`$ nanocrystals are of similar magnitude, while the $`8\times 3\times 3`$ minima are separated from the rest by a considerably larger barrier. The values of some of the barrier heights are given in Table I.

The disconnectivity graph is also helpful for interpreting the (NaCl)<sub>35</sub>Cl<sup>-</sup> dynamics observed in experiments. In the formation process it is likely that a high-energy configuration is initially generated. The cluster then relaxes to one of the low-energy nanocrystals. Simulations for potassium chloride clusters indicate that this relaxation is particularly rapid for alkali halides because of the low barriers (relative to the energy difference between the minima) for downhill pathways. However, there is a separation of time scales between this initial relaxation and the conversion of the metastable nanocrystals to the one with the lowest free energy. the large barriers between the different cuboids make them efficient traps.

### B Thermodynamics of (NaCl)<sub>35</sub>Cl<sup>-</sup>

Some thermodynamic properties of (NaCl)<sub>35</sub>Cl<sup>-</sup> are shown in Figures 3 and 4. The caloric curve shows a feature at $``$700K which indicates melting (Figure 3a); the melting temperature is depressed relative to the bulk due to the cluster’s finite size. The effects of the barriers between nanocrystals are apparent in the MC simulations. The radius of gyration, $`R_g`$, provides a means of differentiating the nanocrystals. It can be seen from the plot of $`R_g`$ for simulations started in the lowest-energy minima of the four cuboidal forms that each simulation is stuck in the starting structure (Figure 3b) up to temperatures close to melting, implying that there are large free energy barriers between the nanocrystals.

These free energy barriers prevent an easy determination of the relative stabilities of the different nanocrystals by conventional simulations. Therefore, we use the superposition method to examine this question. First we assign the fifty lowest-energy minima to one of the cuboidal forms by visual inspection of each structure. Then, using these sets as definitions of the nanocrystals in Equation (2), we calculate the equilibrium probabilities of the cluster being in the different cuboidal morphologies as a function of temperature. It can be seen from Figure 4a that the $`5\times 5\times 3`$ nanocrystal is most stable up until melting. The $`6\times 4\times 3`$ nanocrystal also has a significant probability of being occupied. However, the probabilities for the $`8\times 3\times 3`$ and $`5\times 4\times 4`$ nanocrystals are always small. The onset of the melting transition is indicated by the rise of $`p_{\mathrm{rest}}`$ in Figure 4a and by the peak in the heat capacity (Figure 4b). However, this transition is much too broad and the heat capacity peak occurs at too high a temperature because the incompleteness of our sample of minima leads to an underestimation of the partition function for the liquid-like minima. Given these expected failings of the superposition method at high temperature when the partition functions of the minima are not reweighted, we restrict our dynamics calculations to temperatures below 600K.

Although the probabilities of being in the different morphologies show little variation at low temperature, there are significant changes in the occupation probabilities of specific minima. For example, the small low temperature peak in the heat capacity is a result of a redistribution of probability amongst the low-energy $`5\times 5\times 3`$ minima; the third lowest-energy minimum becomes most populated. It is also interesting to note that the second lowest-energy $`5\times 4\times 4`$ minimum (O) becomes more stable than the lowest-energy $`5\times 4\times 4`$ minimum (L) for temperatures above approximately 220K. Both these changes are driven by differences in vibrational entropy.

### C Dynamics of (NaCl)<sub>35</sub>Cl<sup>-</sup>

Some examples of the interfunnel dynamics that we find on solution of the master equation are depicted in Figure 5. The time scales involved are much longer than those accessible by conventional simulations.

The dynamics of relaxation to equilibrium depend significantly on the starting configuration. When the lowest-energy $`6\times 4\times 3`$ minimum is the initial configuration there is a small transient population in $`5\times 4\times 4`$ minima before the system adopts a $`5\times 5\times 3`$ structure. This is consistent with the lowest-energy pathway that was found between the two nanocrystals; it passes through some intermediate $`5\times 4\times 4`$ minima.

The relaxation to equilibrium is much slower when the cluster starts from the lowest-energy $`8\times 3\times 3`$ minimum. This is a result of the large barrier to escape from this funnel (Figure 2). The probability flow out of the $`8\times 3\times 3`$ funnel leads to a simultaneous rise in the occupation probabilities of both the $`5\times 5\times 3`$ and $`6\times 4\times 3`$ nanocrystals towards their equilibrium values, even though the lowest-energy pathway out of the $`8\times 3\times 3`$ funnel directly connects it to the low-energy $`6\times 4\times 3`$ minima. This occurs because the time scale for interconversion of these latter two nanocrystals is much shorter than that for escape from the $`8\times 3\times 3`$ funnel. However, we do find evidence that the cluster first passes through the $`6\times 4\times 3`$ minima if we examine the probabilities on a log-scale. At times shorter than that required for local equilibrium between the $`5\times 5\times 3`$ and $`6\times 4\times 3`$ minima the occupation probability for the $`6\times 4\times 3`$ minima is larger.

As the two lowest-energy $`5\times 4\times 4`$ minima are well-separated in configuration space we considered relaxation from both these minima. In both cases, there is a large probability flow into the $`6\times 4\times 3`$ minima, which is then transferred to the $`5\times 5\times 3`$ funnel on the same time scale as when initiated in the $`6\times 4\times 3`$ funnel. However, the time scale for the build-up of population in the $`6\times 4\times 3`$ minima depends on the initial configuration. Probability flows more rapidly and directly into the $`6\times 4\times 3`$ minima when initiated from minimum O, reflecting the low barriers between these minima (Figure 2). For relaxation from minimum L there are two active pathways, leading to an increase in the population of both the $`5\times 5\times 3`$ and the $`6\times 4\times 3`$ minima. The direct path into the $`5\times 5\times 3`$ funnel has the lower barrier (Table I) but is long (96.7Å), and so has a smaller rate than the path into the $`6\times 4\times 3`$ funnel, which has a slightly higher barrier (by 0.05 eV). The small shoulder in the occupation probability of the $`5\times 5\times 3`$ minima occurs in the time range when the occupation probability of the $`5\times 4\times 4`$ minima has reached a value close to zero (thus reducing the contribution from the direct path) and when the probability flow out of the $`5\times 4\times 4`$ minima is only just beginning.

The combination of our thermodynamics and dynamics results for (NaCl)<sub>35</sub>Cl<sup>-</sup> enable us to explain why a peak associated with the $`5\times 4\times 4`$ cuboids was not observed experimentally. The $`5\times 4\times 4`$ minima have a shorter lifetime than the other cuboidal forms and have a low equilibrium occupation probability.

Another way to analyse the dynamics is to examine how local equilibration progresses towards the point where global equilibrium has been obtained. To accomplish this we define two minima to be in local equilibrium at the time when

$$\frac{\left|P_i(t)P_j^{\mathrm{eq}}P_j(t)P_i^{\mathrm{eq}}\right|}{\sqrt{P_i(t)P_j(t)P_i^{\mathrm{eq}}P_j^{\mathrm{eq}}}}ϵ$$
(7)   
is obeyed for all later times. In the present work we set $`ϵ=0.01`$, i.e. the two minima are within 1% of equilibrium. Using this definition we can construct equilibration graphs, in which nodes occur when two (groups of) minima come into local equilibrium.

We show an example of an equilibration graph in Figure 6. Equilibration first occurs at the bottom of the $`6\times 4\times 3`$ and $`5\times 5\times 3`$ funnels between those minima that are connected by low barrier paths, then progresses to minima with the same cuboidal shape but which are separated by larger barriers, and finally occurs between the funnels. The order of interfunnel equilibrium agrees with the time scales that we observe in the time evolution of the occupation probabilities of the nanocrystals (Figure 5). Minimum O, then minimum L, come into equilibrium with the $`6\times 4\times 3`$ funnel. Then, the $`5\times 5\times 3`$ and $`6\times 4\times 3`$ funnels reach local equilibrium. Finally, the $`8\times 3\times 3`$ funnel reaches equilibrium with the rest of the PES. As one of the major determinants of the time scale required for local equilibrium is the height of the barriers between minima, it is unsurprising that the equilibration graph reflects the structure of the disconnectivity graph (Figure 2).

In the experiments on (NaCl)<sub>35</sub>Cl<sup>-</sup> rate constants and activation energies were obtained for the conversion of the $`6\times 4\times 3`$ and $`8\times 3\times 3`$ nanocrystals into the $`5\times 5\times 3`$ nanocrystal. It would, therefore, be useful if we could extract rate constants for the different interfunnel processes from the master equation dynamics.

For a two-state system, where $`AB`$ and $`k_+`$ and $`k_{}`$ are forward and reverse rate constants, respectively, it can be shown that

$$\mathrm{ln}\left[\frac{P_A(t)P_A^{\mathrm{eq}}}{P_A(0)P_A^{\mathrm{eq}}}\right]=(k_++k_{})t,$$
(8)   
and the equivalent expression for B are obeyed. This is a standard result for a first-order reaction. This expression will also hold for the rate of passage between two funnels in our multi-state system if the interfunnel dynamics are the only processes affecting the occupation probabilities of the relevant funnels, and if the interfunnel dynamics cause the occupation probabilities of the two funnels to converge to their equilibrium values.

In Figure 7 we test the above expression by applying it to the interconversion of $`6\times 4\times 3`$ and $`5\times 5\times 3`$ nanocrystals. The two lines in the graph converge to the same plateau value, before both falling off beyond 0.001s. This plateau corresponds to the time range for which the interfunnel passage dominates the evolution of the probabilities for the two funnels. At shorter times, when the occupation probabilities for the two funnels are still close to their initial values, there are many other contributing processes. At longer times the probabilities are both very close to their equilibrium values, and the slower equilibration with the $`8\times 3\times 3`$ funnel dominates the probability evolution. From the plateau in Figure 7 we obtain $`k^++k^{}=5320\mathrm{s}^1`$. The individual rate constants can be obtained by using the detailed balance relation: $`k^+P_A^{\mathrm{eq}}=k^{}P_B^{\mathrm{eq}}`$.

The application of Equation (8) to escape from the $`8\times 3\times 3`$ funnel also leads to a range of $`t`$ where there is a well-defined plateau. However, this approach works less well for the interconversion of $`5\times 4\times 4`$ and $`6\times 4\times 3`$ minima (Figures 5a and b). This is because the assumption that no other processes contribute to the probability evolution of the two funnels is obeyed less well, and because the occupation probabilities do not converge to their equilibrium values, but near to the equilibrium values that would obtain if the $`5\times 5\times 3`$ minima were excluded. Nevertheless, approximate values of $`k^+`$ and $`k^{}`$ can be obtained.

Diagonalization of the matrix $`\stackrel{~}{𝐖}`$ produces a set of eigenvalues that give the time scales for a set of characteristic probability flows. The dynamical processes to which the eigenvalues correspond can be identified by examining the eigenvectors. Flow occurs between those minima for which the corresponding components of the eigenvector have opposite sign (Equation (5)). This observation forms the basis for Kunz and Berry’s net-flow index which quantifies the contribution of an eigenvector $`i`$ to flow out of a funnel $`A`$. The index is defined by

$$f_i^\mathrm{A}=\underset{j\mathrm{A}}{}\stackrel{~}{u}_j^{(i)}\sqrt{P_j^{\mathrm{eq}}}.$$
(9)   
The index allows the interfunnel modes to be identified; the values of $`f^A`$ and $`f^B`$ for these modes will be large and of opposite sign. For example, at $`T`$=400K the mode with the most $`5\times 5\times 36\times 4\times 3`$ character has $`f^{5\times 5\times 3}=0.339`$, $`f^{6\times 4\times 3}=0.331`$ and $`\lambda =5275\mathrm{s}^1`$. The eigenvalue is in good agreement with the sum of interfunnel rates obtained using Equation (8).

The extraction of the interfunnel rate in this manner is hindered by the fact that the eigenvalues of $`\stackrel{~}{𝐖}`$ cannot cross as a function of temperature. Instead, there are avoided crossings and mixing of modes. For example, the small difference between the two values for the sum of the $`5\times 5\times 36\times 4\times 3`$ interfunnel rate constants that we obtained above is probably due to mixing. The eigenvector with which the interfunnel mode mixes gains some $`5\times 5\times 36\times 4\times 3`$ character; it has $`f^{5\times 5\times 3}=0.014`$, $`f^{6\times 4\times 3}=0.017`$ and $`\lambda =6388\mathrm{s}^1`$.

By calculating and diagonalizing $`\stackrel{~}{𝐖}`$ at a series of temperatures we can examine the temperature dependence of the interfunnel rate constants. For most of the interfunnel processes, the rate constants that we obtain fit well to the Arrhenius form, $`k=A\mathrm{exp}(E_a/kT)`$, where $`E_a`$ is the activation energy and $`A`$ is the prefactor (Figure 8).

In our previous study of the interconversion mechanisms in (NaCl)<sub>35</sub>Cl<sup>-</sup> we estimated the activation energies in order to compare with the experimental values. Using the analogy to a simple one-step reaction, we equated the activation energy with the difference in energy between the highest-energy transition state on the lowest-energy path between the relevant nanocrystals and the lowest-energy minimum of the starting nanocrystal. However, it was not clear how well this analogy would work. In Table I we compare these estimates with the activation energies obtained from our master equation results. There is good agreement, confirming the utility of the approximate approach. Similar agreement has also been found for the interfunnel dynamics of a 38-atom Lennard-Jones cluster.

A simple explanation for this correspondence can be given. We first label the minima along the lowest-barrier path between the two funnels in ascending order and define $`l`$ so that the highest-energy transition state on this path lies between minima $`l1`$ and $`l`$. If the minima behind the highest-energy transition state (i.e. 1 to $`l1`$) are in local equilibrium, then the occupation probability of the minimum $`l1`$ is given by

$$\frac{p_{l1}}{p_1}\mathrm{exp}\left((E_{l1}E_1)/kT\right)$$
(10)   
Then, if the rate of interfunnel flow, $`k^+p_A`$, is equated to the rate of passage between minima $`l1`$ and $`l`$,

$`k^+p_A`$ $``$ $`p_{l1}\mathrm{exp}((E_{l1,l}^{}E_{l1})/kT)`$ (11)   
$``$ $`p_1\mathrm{exp}((E_{l1,l}^{}E_1)/kT),`$ (12)   
Therefore, if the occupation probability for funnel A is dominated by the occupation probability for the lowest-energy minimum in the funnel, i.e. $`p_Ap_1`$ for all $`T`$, then the activation energy is equal to the energy difference between the highest-energy transition state on the lowest-barrier path and the energy of the lowest-energy minimum in the starting funnel.

We should note that in the above derivation the interfunnel probability flow is assumed to all pass through a single transition state. However, if there is competition between two paths, one with a low barrier and a small prefactor and one with a larger barrier and a large prefactor, we expect that the low-barrier path would dominate at low temperature and the high-barrier path at high temperature. This behaviour would give rise to a interfunnel rate constant with a positive curvature in an Arrhenius plot. However, the lines are either straight or have a small amount of negative curvature.

The lack of positive curvature, and the agreement between the estimated and the observed activation energies, probably indicates that the interfunnel probability flow is dominated by paths which pass through the highest-energy transition state on the lowest barrier path. It is interesting that, on a PES with so many minima and transition states, a single transition state can have such a large influence on the dynamics.

At low enough temperature, $`d(p_1/p_A)/dT0`$, and so the interfunnel barrier height can be measured with respect to the lowest-energy minimum in the starting funnel (Equation (12)). However, as the occupation probabilities of other minima in the funnel becomes significant relative to that for the lowest-energy minimum in the funnel, the ratio $`p_1/p_A`$ decreases, thus giving the lines in the Arrhenius plot their slight negative curvature. In other words, the apparent activation energy decreases with increasing temperature, because the barrier height should be measured with respect to some kind of average minimum energy for the funnel, perhaps $`E_A=_{iA}p_iE_i`$. The negative curvature is most pronounced when the occupation probabilities in a funnel change considerably. For example, the $`5\times 4\times 4(L)6\times 4\times 3`$ rate constant has the most curvature because minimum L has a particularly low vibrational entropy leading to population of other minima within that funnel (Figure 2).

In Table I we also compare our activation energies to the two experimental values. Our values are too large by 0.24 and 0.49 eV. There are a number of possible sources of error. Firstly, the samples of minima and transition state provide only an incomplete characterization of the PES and so it is possible that the nanocrystals are connected by undiscovered lower-barrier paths. However, we believe that our sample of minima is a good representation of the low-energy regions of the PES and consider it improbable that undiscovered pathways could account for all of the discrepancy. Secondly, in our calculations the input rate constants $`k_{ij}^{}`$ were calculated on the basis of the harmonic approximation. Although this is likely to have a significant effect on the absolute values of the interfunnel rate constants and prefactors (the latter are too large compared to the experiment values), it should not have such a significant effect on the activation energies. Instead, we consider the most likely source of the discrepancy between theory and experiment to be inaccuracies in the potential. When polarization is included by using the Welch potential, the estimated barriers become closer to the experimental values (the discrepancies are then 0.16 and 0.34 eV) but significant differences still remain. For better agreement we may need a potential that allows the properties of the ions to depend on the local environment. Unfortunately, although such potentials have been developed for a number of systems, one does not yet exist for sodium chloride.

## IV Conclusion

The (NaCl)<sub>35</sub>Cl<sup>-</sup> potential energy surface has a multiple-funnel topography. Structurally, the different funnels correspond to rocksalt nanocrystals with different cuboidal forms. The large potential energy barriers between the funnels causes the time scales for escape from metastable nanocrystals to be far longer than those accessible by conventional dynamics simulations. Therefore, we examined the interfunnel dynamics by applying the master equation approach to a database of (NaCl)<sub>35</sub>Cl<sup>-</sup> minima and transition states. The slowest rate constant we obtained was $`k_{5\times 5\times 38\times 3\times 3}=1.46\times 10^8\mathrm{s}^1`$ at $`T=275`$K.

Using a net flow index we were able to identify the eigenvalues of the transition matrix, $`𝐖`$, which correspond to interfunnel probability flow. Thus, we were able to obtain rate constants and activation energies for the interconversion of the different nanocrystals. One particularly interesting finding is that the activation energies correspond fairly closely to the potential energy differences between the highest-energy transition state on the lowest-energy path between two nanocrystals and the lowest-energy minimum of the starting nanocrystal. This is the result one might expect by a simple extrapolation from the dynamics of a simple molecular reaction. However, it holds despite the multi-step, and potentially multi-path, nature of the interfunnel dynamics. The question of whether this result is generally true for interfunnel dynamics involving large potential energy barriers or reflects some of the particulars of the (NaCl)<sub>35</sub>Cl<sup>-</sup> system is an interesting subject for further research. We already know that this simplification holds for the 38-atom Lennard-Jones cluster.

###### Acknowledgements.

J.P.K.D. is the Sir Alan Wilson Research Fellow at Emmanuel College, Cambridge. D.J.W. is grateful to the Royal Society for financial support. We would like to thank Mark Miller for helpful discussions.