# Evolution of the Potential Energy Surface with Size for Lennard-Jones Clusters

## I Introduction

The relationship between dynamics and the potential energy surface (PES), or potential energy ‘landscape’, holds the key to understanding a wide range of phenomena, including protein folding, global optimization and the properties of glasses. For example, what are the topographical features of the PES that enable a protein to fold to its native state? What sort of PES renders global optimization tractable, and how can it be transformed to make global optimization easier? Which features of the PES produce good glass formers, and what distinguishes ‘strong’ from ‘fragile’ liquids?

Stillinger and Weber first showed the utility of partitioning configuration space into basins of attraction surrounding each local minimum. The thermodynamics of the system can then be formulated in terms of the properties of these local minima. However, the minima alone contain no dynamic information. Therefore, to obtain insight into the dynamics it is necessary to locate the transition states which connect the local minima. This information can then be incorporated into a master equation that describes the flow of probability between minima as the system evolves towards equilibrium. As well as providing a picture of the relationship between the dynamics and the PES, this approach can probe time scales much longer than those accessible in conventional simulations.

We also wish to relate the thermodynamics and dynamics to general features of the PES in a qualitative manner. Disconnectivity graphs are a recently-developed tool that provide a helpful visualization of the PES. These graphs have now been calculated for a number of polypeptides and various clusters. They show which of the minima in a sample are connected by pathways lying below a series of energy thresholds, and so provide a picture of the hierarchy of energy barriers in a system.

In this paper we present disconnectivity graphs for a series of Lennard-Jones (LJ) clusters. Our aims are to gain more insight into the potential energy landscapes of these clusters and to explore systematically how the graphs evolve with size. LJ clusters provide an ideal test case because their dynamics and thermodynamics have been thoroughly studied and they have been used in the development of global optimization algorithms.

As the size of the cluster increases the number of minima and transition states is expected to grow at least exponentially, even when permutational isomers are not included. Furthermore, cluster properties do not vary smoothly with the number of atoms in this size regime, and we have chosen to examine sizes which should illustrate particularly interesting features.

Due to the exponential increase in the number of stationary points with size it is difficult for a disconnectivity graph to retain a global picture of the PES. We have therefore explored the use of monotonic sequence basins as the fundamental topographical unit rather than minima (Section III B). Finally, in Section III C we use disconnectivity graphs to probe the transformed energy landscape that is searched by the Monte Carlo minimization or basin-hopping algorithm.

## II Methods

The atoms interact via a Lennard-Jones potential:

$$E_c=4ϵ\underset{i<j}{}\left[\left(\frac{\sigma }{r_{ij}}\right)^{12}\left(\frac{\sigma }{r_{ij}}\right)^6\right],$$
(1)   
where $`ϵ`$ is the pair well depth and $`2^{1/6}\sigma `$ is the equilibrium pair separation.

To examine the topography of the PES, we need to locate its minima and the network of transition states and pathways that connect them. Transition states can be found efficiently using eigenvector-following, in which the energy is maximized along one direction and simultaneously minimized in all the others. The minima connected to a transition state are defined by the two steepest-descent paths which begin parallel and antiparallel to the unique Hessian eigenvector whose corresponding eigenvalue is negative. The steepest-descent paths were calculated using a method which employs analytic second derivatives.

The samples of minima and transition states were generated using an approach which has also been described before. The elementary step in this process is to perform a transition state search after stepping away from a minimum along one of the Hessian eigenvectors and then, if this search finds a transition state, to generate the corresponding pathway. If the transition state is connected to the minimum from which the search was started, we add it to our database along with the connected minimum. Occasionally we find that the transition state is not connected to the starting minimum. In this case we only add it to the database if it is connected to a minimum that is already present. At each step we start from the lowest-energy minimum from which fewer than a specified number, 2$`n_{\mathrm{ev}}`$, of transition state searches have been performed. Searches are carried out in positive and negative directions along each eigenvector in order of increasing eigenvalue. The value chosen for $`n_{\mathrm{ev}}`$ determines how thoroughly the PES is searched in the vicinity of a given minimum.

By repeating this process until no new minima are found we can obtain a nearly exhaustive catalogue of the minima. This approach was used for the 13-atom cluster, LJ<sub>13</sub>, but for the larger clusters it is impossible to find all the minima. We therefore stopped searching once we were confident that we had obtained an accurate representation of the low-energy regions of the PES. For LJ<sub>75</sub> an alternative strategy was required due to the double-funnel character of the landscape; this approach is described in the following section.

## III Results

We now characterize the PES’s of LJ<sub>13</sub>, LJ<sub>19</sub>, LJ<sub>31</sub>, LJ<sub>38</sub>, LJ<sub>55</sub> and LJ<sub>75</sub>. To explain this selection of cluster sizes, and as a background to the interpretation of the disconnectivity graphs, we briefly review some of the structural properties of small LJ clusters. Comparisons of particularly stable sequences of clusters indicate that structures based on Mackay icosahedra (e.g. 13.1 and 55.1 in Figure 1) are dominant up to $`N1600`$. Therefore, as can be seen from Figure 2, most of the clusters in the size range that we are considering have a global minimum based on icosahedral packing. The sizes at which complete Mackay icosahedra are geometrically possible ($`N`$=13, 55,$`\mathrm{}`$) are particularly stable compared to other sizes, leading to particularly low-energy global minima (Figure 2). For these clusters, there is a large energy gap between the two lowest-energy minima ($`2.85ϵ`$ and $`2.64ϵ`$ for LJ<sub>13</sub> and LJ<sub>55</sub>, respectively) and we expect the PES to have a single deep funnel.

The term ‘funnel’ was introduced in the protein folding literature to describe the situation where a collection of downhill pathways converge on the native structure of the protein. As a result of this PES topography the protein is ‘funnelled’ towards the native structure on relaxation. We use the term here in a similar manner except that, instead of converging on the native structure of the protein, the funnel should converge on a single low-energy minimum, or possibly to a collection of structurally-related low-energy minima.

Between $`N`$=13 and 55 the energy first increases relative to an energy function interpolated between the Mackay icosahedra as an overlayer is added to the 13-atom icosahedron, and then decreases as the overlayer approaches completion at $`N`$=55 (Figure 2). Of course, there are smaller features superimposed on this broad maximum. At certain sizes particularly stable overlayers with no low-coordinate atoms can be formed, giving rise to subsidiary minima in Figure 2. For example, at $`N`$=19 a six-atom cap on the 13-atom icosahedron leads to the relatively stable double icosahedron shown in Figure 1. We would expect LJ<sub>19</sub> to exhibit a single-funnel PES but to relax less efficiently to the global minimum than LJ<sub>13</sub> or LJ<sub>55</sub>.

We also chose to study a size, $`N`$=31, from near the top of the broad maximum in Figure 2, where there are a large number of low-energy minima associated with the various ways that the atoms of the overlayer on the 13-atom icosahedron can be arranged. Moreover, for this size the lowest-energy decahedral cluster is only a little higher in energy than the icosahedral structures.

For most sizes the global minimum has icosahedral character, but there are a few clusters for which this is not the case. For $`N`$=38 the global minimum is the face-centred-cubic (fcc) truncated octahedron (38.1 in Figure 1) and for $`75N77`$ the global minima are based on the Marks decahedron (75.1 in Figure 1). We chose to study $`N`$=38 and 75 where the global minimum is structurally very different from the second lowest-energy minimum, which in each case is icosahedral. The two lowest-energy minima are therefore rather distant in configuration space, and they are separated by large potential energy and free energy barriers. Therefore, we expect these surfaces to have two funnels: one which leads to the low-energy icosahedral minima and one which leads to the global minimum.

Details of the samples of minima and the transition states that we found for each of the clusters are given in Table I. For LJ<sub>13</sub>, we were able to gauge how close the samples were to convergence by following the number of minima and transition states as $`n_{\mathrm{ev}}`$ was increased. The number of minima seemed close to reaching an asymptote, but the number of transition states was still increasing at $`n_{\mathrm{ev}}`$=15.

### A Disconnectivity graphs

The conceptual basis for the disconnectivity graph approach is as follows. At a given total energy, $`E`$, minima can be grouped into disjoint sets, called superbasins, whose members are mutually accessible at that energy. In other words, each pair of minima in a superbasin are connected directly or through other minima by a path whose energy never exceeds $`E`$, but would require more energy to reach a minimum in another superbasin. At low energy there is just one superbasin—that containing the global minimum. At successively higher energies, more superbasins come into play as new minima are reached. At still higher energies, the superbasins coalesce as higher barriers are overcome, until finally there is just one containing all the minima (provided there are no infinite barriers).

To construct a disconnectivity graph, the superbasin analysis is performed at a series of energies. 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. The horizontal position of a node has no significance, and is chosen for clarity. Every line terminates at a local minimum.

The disconnectivity graphs produced using this procedure are shown in Figure 3. For $`N`$=13 it is possible to include all the minima that we have found in the disconnectivity graph, which therefore provides a practically complete global picture of the PES. The graph has the form expected for an ideal funnel: there is a single stem, representing the superbasin of the global minimum, with branches sprouting directly from it at each level, indicating the progressive exclusion of minima as the energy is decreased.

The form of this graph implies that all the minima are directly connected to the superbasin associated with the global minimum. In fact, 99% of the minima are within two rearrangements of the global minimum and none are further than three steps away. Furthermore, there are 911 structurally distinct transition states connecting 535 minima directly to the global minimum. If the reaction path degeneracy is taken into account, there are $`\mathrm{108\hspace{0.17em}967}`$ transition states (some of which are permutational isomers) connected to each permutational isomer of the global minimum.

If we examine the distance of a minimum from the global minimum as a function of its potential energy, for LJ<sub>13</sub> the resulting plot (Figure 4a) has the form we would expect of a funnel: the distance from the global minimum increases as the potential energy of the minimum increases. It is also worth noting that the slope of the funnel is steeper than for any of the other clusters. Hence it is no surprise that relaxation down the PES to the global minimum is relatively easy for this cluster. However, relaxation is hindered somewhat by the fairly large barriers (1–2$`ϵ`$) that exist for escaping from some of the minima into the superbasin of the global minimum (Figure 3a).

For all the other clusters it is not computationally feasible to obtain a nearly complete set of minima. Moreover, if we attempt to represent all the minima of our samples on a disconnectivity graph, the density of lines simply becomes too great. Instead we only represent those branches that lead to a specified number of the lowest-energy minima, both for clarity and because our samples of minima and transition states are likely to be most complete for the low-energy regions of the PES. We should note that the minima that are not represented can still contribute to the appearance of the graph if they mediate low barrier paths between minima that are included. However, this approach does have the consequence that, as the size increases, the graphs increasingly focus on a smaller and smaller proportion of the surface.

We can see this effect if we examine the graph of LJ<sub>55</sub> for which we also expect the PES to exhibit a single funnel. Although there is more fine structure, the form is similar to that of LJ<sub>13</sub>. However, all the minima represented in the disconnectivity graph have relatively ordered structures and so, unlike the graph for LJ<sub>13</sub>, this representation does not tell us whether there is a funnel leading from the disordered liquid-like minima to the global minimum Mackay icosahedron. Instead, it only shows that the low-energy region of the PES associated with structures based on the Mackay icosahedron is funnel-like. The graph probably only represents the bottom of a larger funnel leading down from the huge number of disordered minima. The plot of the distance from the global minimum as a function of potential energy shows a glimpse of this larger funnel (Figure 4e).

It is easy to understand the differences between the graphs for LJ<sub>13</sub> and LJ<sub>55</sub>. The energy gap between the global minimum and the disordered liquid-like minima is an extensive quantity, whereas the energy gap between two ordered structures, where a few atoms have changed position, is related to the change in the number of nearest-neighbour contacts and is independent of size. Furthermore, the number of possible ways of arranging surface atoms and vacancies in an ordered structure becomes increasingly large. Therefore, the number of ordered minima which lie below the liquid-like band of minima increases rapidly with the number of atoms.

The fine structure of the LJ<sub>55</sub> disconnectivity graph reveals some interesting features. The minima separate into bands related to the number of defects present in the Mackay icosahedron. For example, all the minima in the first band above the global minimum are Mackay icosahedra with a missing vertex and an atom on the surface. The eleven minima in this band correspond to the eleven possible sites for this atom that are unrelated by symmetry. Of these minima the four lowest in energy have the atom located over the centre of one of the faces of the Mackay icosahedron, and the other seven have the atom off-centre. In the disconnectivity graph these minima split into four groups corresponding to the four symmetry-unrelated faces on which the adatom can be located. The splitting occurs because the barriers for rearrangements in which the adatom passes between faces are larger than the barriers for changing the position of the adatom on a face. Therefore, in the disconnectivity graph the minima first become connected to the other (one or two) minima with the adatom on the same face before they are connected to the stem associated with the global minimum.

The second band of minima consists of Mackay icosahedra with two missing vertices and two surface atoms. The lower-lying minima in this band have the two adatoms in contact, either on the same face or bridging an edge. The minima with the two adatoms unpaired give rise to a repeated motif, an example of which is illustrated in the inset. This feature is repeated over ten times on the left-hand side of the graph with the top node always at $`273.0ϵ`$. The minima split into these sets because of the lower barriers for an adatom moving between sites on the same face. If the two faces with adatoms are unrelated by symmetry, there are sixteen distinct ways of arranging the atoms on the two faces. The lowest-energy minimum of these sixteen has the two adatoms in the central sites. Slightly higher in energy are the six minima with one adatom central and one off-centre. Higher still are the nine minima with both adatoms off-centre.

The features noted above are reflected in the dynamics of LJ<sub>55</sub>. Just below the melting point the cluster passes between the Mackay icosahedron and states with one and two defects, residing in each state for periods of the order of nanoseconds (using potential parameters appropriate to argon). While in one of these states, the surface atoms and vacancies diffuse across the surface. This time scale separation between diffusive motion and the formation and annihilation of defects is a consequence of the higher barriers for the latter processes and the distances over which a surface atom and vacancy must diffuse to recombine. The latter is reflected in the wide range of distances that minima in this second band are away from the global minima.

The two minima (55.13 and 55.149) in Figure 3e that must overcome the highest barriers to reach the superbasin containing the global minimum do not have icosahedral structures. Both have $`D_{5h}`$ point group symmetry and are depicted in Figure 1. The lower-energy of the two is actually the thirteenth lowest-energy minimum and was recently discovered by Wolf and Landman. The other minimum has previously been noted by Wales. Both can be converted into a Mackay icosahedron by a single cooperative rearrangement in which parts of the structure twist around the fivefold axis. As a significant number of nearest-neighbour contacts are broken in these rearrangements, the barriers are higher than those for the localized rearrangements by which the defective Mackay icosahedra interconvert.

The disconnectivity graph for LJ<sub>19</sub> shows that the PES is again funnel-like (Figure 3b), as expected. The branches for most of the minima connect directly to the superbasin containing the global minimum. Although our graph only shows branches leading to the lowest 250 minima, Figure 4b reveals that the funnel continues up to higher energies.

The PES of LJ<sub>19</sub> has been analysed previously in several studies. The profiles of the downhill pathways to the global minimum were described as ‘sawtooth-like’ rather than ‘staircase-like’ because the barrier heights are relatively large compared to the energy differences between the minima. On this basis Ball et al. concluded that LJ<sub>19</sub> has topographical features typical of a glass-former. The disconnectivity graph also shows that some of the downhill barriers are quite large. However, as the global minimum is at the bottom of a funnel the barriers only slow down the rate of relaxation towards the double icosahedron, rather than preventing it. Lowering the energy does not take the system away from the global minimum, but rather towards it (Figure 4b). The appearance of $`N=19`$ as a magic number in mass spectra of rare gas clusters confirms that the global minimum is kinetically accessible.

The disconnectivity graph for LJ<sub>31</sub> (Figure 3c) is fundamentally different from those we have considered so far. The energetic bias towards the global minimum is smaller (Figure 4a). In fact there are a number of minima with energies close to that of the global minimum which are separated from it by fairly large energy barriers. This situation is partly the result of competition between two distinct types of overlayer. In the first type, the anti-Mackay overlayer, atoms are added to the faces and vertices of the underlying 13-atom icosahedron (giving rise, for example, to the double icosahedron, 19.1). In the second type, the Mackay overlayer atoms are added to the edges and vertices. The completion of the Mackay overlayer leads to the next Mackay icosahedron. LJ<sub>31</sub> is the first size for which a cluster with the Mackay overlayer is the global minimum. It can be seen from Figure 1 that minimum 31.1 is a fragment of the 55-atom Mackay icosahedron. The second lowest-energy minimum, 31.2, has an anti-Mackay overlayer.

There are also some low-energy decahedral minima for LJ<sub>31</sub>. For example, minimum 31.9 (Figure 1) is separated by a large energy barrier from the global minimum, because not only must there be a change in morphology from decahedral to icosahedral, but the cluster must also change shape. Some of the decahedra with more spherical shapes are connected to the superbasin associated with the global minimum by smaller barriers.

We can deduce something of the relaxation behaviour of LJ<sub>31</sub> from its disconnectivity graph. Once the cluster has reached a low-energy configuration, presumably by rapid descent of a funnel from the liquid, subsequent relaxation towards the global minimum may be considerably slower. There is little energetic bias at the bottom of the PES to guide the system towards the global minimum and the barriers for interconversion of the low-energy minima can be relatively large. Therefore, it is not surprising that the time required to find the global minimum using the basin-hopping global optimization algorithm shows a maximum at $`N`$=31.

The effects of competing structures that we noted for LJ<sub>31</sub> appear in a more extreme form for LJ<sub>38</sub> and LJ<sub>75</sub>. For both these clusters the disconnectivity graph clearly separates the low-energy minima into two main groups, namely those associated with the global minimum and those with icosahedral structure. These two groups of minima are separated by a large energy barrier, so the graph splits into two stems at high energy which lead down to two structurally distinct sets of low-energy minima. This splitting is characteristic of a multiple funnel PES. The separation is particularly dramatic for LJ<sub>75</sub>, where the decahedral to icosahedral barrier is over $`3ϵ`$ larger than any of the other barriers between the 250 lowest-energy minima. Although the disconnectivity graphs clearly show only the bottom of the PES, these two groups of minima can be associated with separate funnels, and give rise to distinct thermodynamic states. The double funnel structure is also apparent from the plots in Figure 4d and f. Passing from the global minimum to the icosahedral funnel, the energy first increases as the primary funnel is ascended and then decreases during the descent into the second funnel.

From the relative numbers of minima associated with each funnel it is clear that the icosahedral funnel is much wider. This is one of the reasons why relaxation is much more likely to lead to the icosahedral minima. There are also thermodynamic effects: the region containing the icosahedral minima becomes lowest in free energy at low temperature ($`T`$$``$$`0.12ϵk^1`$ for LJ<sub>38</sub> and $`T`$$``$$`0.09ϵk^1`$ for LJ<sub>75</sub>) because of the entropy that arises from the large number of low-energy icosahedral minima. Between this transition temperature and the melting point there is therefore a thermodynamic driving force towards the icosahedral structures. Moreover, for LJ<sub>38</sub> we have shown that the free energy barrier for entering the icosahedral region of configuration space is lower than for entering the fcc funnel.

Once the cluster enters the icosahedral funnel it is likely to be trapped because of the large energy and free energy barriers to passing between the two funnels. The energy barriers are $`4.22ϵ`$ and $`3.54ϵ`$ for LJ<sub>38</sub> and $`8.69ϵ`$ and $`7.48ϵ`$ for LJ<sub>75</sub>. At higher temperatures the entropy of the intermediate states reduces the free energy barriers from these zero temperature limits. For the above reasons LJ<sub>38</sub> and LJ<sub>75</sub> provide relatively difficult test cases for any unbiased global optimization algorithm. It is only relatively recently that algorithms have begun to find the truncated octahedron, and only one unbiased method and one method that involves seeding have reported finding the Marks decahedron.

The difficulty in finding the global minimum for these clusters is illustrated by statistics for the basin-hopping algorithm. The time required to find the global minimum for LJ<sub>38</sub> is a maximum with respect to neighbouring sizes and for LJ<sub>75-77</sub> the time is so long that good statistics for the first passage time have not yet been obtained. Global optimization is an order of magnitude more difficult for LJ<sub>75</sub> than LJ<sub>38</sub> for a combination of reasons. First, LJ<sub>75</sub> has a much larger search space and a much greater number of minima than LJ<sub>38</sub>. Second, the temperature at which the global free energy minimum becomes associated with the global potential energy minimum lies further below the melting point ($`T_\mathrm{m}`$$``$$`0.17ϵk^1`$ for LJ<sub>38</sub> and $`T_\mathrm{m}`$$``$$`0.29ϵk^1`$ for LJ<sub>75</sub>). Third, the path between the global minimum and the lowest-energy icosahedral structure is also longer, more complicated and higher in energy than for LJ<sub>38</sub>. The transition involves not only a change in morphology, but also a change in shape—the Marks decahedron is oblate whereas the icosahedral structures are prolate. The pathway therefore involves both cooperative rearrangements, where the structure twists around a quasi-fivefold axis, and surface diffusion steps. The lowest-energy pathway that we found between the two lowest-energy minima is depicted in Figure 5a. It passes through sixty-five transition states. There are significantly shorter paths between the same two minima (the shortest is $`41.5\sigma `$), but these all involve higher barriers. There is also a significant difference in the character of the pathways between the minima at the bottom of the two funnels for LJ<sub>38</sub> and LJ<sub>75</sub>. All the structures along the LJ<sub>75</sub> pathway are ordered and solid-like, whereas at its highest points the LJ<sub>38</sub> pathway passes through disordered structures.

For LJ<sub>75</sub> the highest-energy minimum on the lowest-energy interfunnel path lies at position $`\mathrm{10\hspace{0.17em}909}`$ in terms of an energy ranking of the minima in our sample, where 1 is the global minimum. Many more low-lying minima could have been found if we had searched the icosahedral region of configuration space more intensively. Therefore, if we had simply performed consecutive transition state searches from minima in order of their energetic rank it is unlikely that a pathway between the two funnels could have been found—rather the search would have most likely been stuck in one of the funnels cataloguing the multitude of ordered structures with that morphology. To find a pathway connecting the two lowest-energy minima we therefore had to bias the search to probe intermediate structures. We first constructed a series of decahedral minima which were decreasingly prolate and then increasingly oblate. We then started searching the PES around these structures with the aim of finding pathways connecting them both to the Marks decahedron and to the low energy icosahedral structures.

Once we had found a pathway between the two structures, we only performed transition state searches from those minima that were connected to either of the two lowest-energy minima by a pathway that had no transition states higher than the highest-energy transition state on the lowest-energy path between structures 75.1 and 75.2. Moreover, this search concentrated on those minima in the set that had an intermediate value of the bond-order parameter, $`Q_6`$. By this procedure, increasingly low-energy paths were found between the two funnels. There is, of course, no guarantee that we have found the lowest-energy pathway, but we doubt if a significantly lower one exists.

From Figure 5b we can see how the bond-order parameter can separate the minima into various groups. The group with $`Q_6`$$``$0.02 correspond to icosahedral minima with an anti-Mackay overlayer (such as 75.2), whereas $`Q_6`$=0.306 for the Marks decahedron. The decahedral minima based upon 75.1 have similar values, these being generally lower for structures that are less oblate. There are also bands of minima at intermediate values of $`Q_6`$. For example, icosahedral structures with a Mackay overlayer have $`Q_6`$$``$0.15. Some of the other minima with intermediate values have in part motifs similar to 55.13, and are connected to icosahedral or decahedral minima by rearrangements involving concerted twists around quasi-fivefold axes.

This distribution of $`Q_6`$ values provides us with another way of visualizing the double-funnel structure of the LJ<sub>75</sub> PES. Figure 5c maps out the energy of the transition state at the top of the lowest-barrier pathway from each minimum to either 75.1 or 75.2, whichever is lower, as a function of $`Q_6`$ for the minimum. The plot shows that the barrier separating a decahedral minimum from the Marks decahedron increases as the value of $`Q_6`$ deviates further from the value for the global minimum. This trend continues until $`Q_6`$$``$0.16. Below this value some of the minima have lower barriers for paths to minimum 75.2, and the energies generally decrease as the value of $`Q_6`$ approaches that of 75.2 from above. The icosahedral minima with a Mackay overlayer stand out in the plot as they have both an intermediate value of $`Q_6`$ and low barriers connecting them to minimum 75.2. The transformations between the two types of icosahedral overlayer are usually achieved by a concerted twisting of the overlayer around one of the vertices of the underlying icosahedron.

### B Coarse-graining the PES

As the size of the cluster increases our disconnectivity graphs focus on an increasingly small proportion of the whole PES to avoid being swamped by the rapidly increasing number of minima. However, it would be desirable to retain a more global picture of the PES. To do so, the disconnectivity graphs would need to be based not on the barriers between minima, but between larger topographical units. For example, Levy and Becker ‘diluted’ their sample of hexapeptide conformations by removing conformations that were energetically and structurally similar. However, this approach requires a meaningful measure of similarity, which is probably harder to devise for a cluster than a molecule with a bonded framework.

Here, we explore the use of monotonic sequences, a concept introduced by Berry and coworkers, to produce a more coarse-grained picture of the PES. Monotonic sequences are series of connected minima where the potential energy decreases with every step. The collection of sequences leading to a particular minimum defines a ‘basin’. To avoid confusion with the various other ‘basins’ that have been defined, we will always refer to such a set as a monotonic sequence basin (MSB). The MSB leading to the global minimum is termed ‘primary’, and is separated from neighbouring ‘secondary’ MSB’s by transition states lying on a ‘primary divide’, and so on. It is important to realise that, above such a divide, it is possible for a minimum to belong to more than one MSB through different monotonic sequences.

For the division of the PES provided by a monotonic sequence analysis to be useful in an analysis of the dynamics, transitions between minima in an MSB must be more rapid than transitions between different MSB’s. Kunz and Berry found evidence for this in a simplified model of LJ<sub>19</sub>. However, there is no guarantee that this separation will always hold, because MSB’s are defined by the connections between minima without reference to the size of the barriers between them.

Disconnectivity graphs that only include the minima at the bottom of each MSB can be produced by excluding all the minima directly connected to a lower-energy minimum. In the resulting graphs, branches are joined by a node at the energy of the lowest transition state on the divide between the MSB’s. For the clusters we consider here, the number of MSB’s is small enough (Table I) for all of them to be represented on the graph. Two examples are shown in Figure 6.

One implication of the method we use to generate our samples of stationary points is that higher-energy minima that were not used as starting points for transition state searches are unlikely to lie at the bottom of an MSB. These minima lie at the higher end of a pathway that was found in a search from a lower-energy minimum. Therefore, most of these minima are directly connected to a lower-energy minimum, and so cannot be at the bottom of an MSB. Occasionally, the pathways do not connect back to the starting minimum but to another unsearched minimum in the sample. Only in these rare instances can an unsearched minimum be at the bottom of an MSB in our sample. Therefore, although the disconnectivity graphs based upon MSB’s provide a more global picture of the PES than the graphs in the previous section, they are still limited by the incompleteness of our samples. The graphs only provide a reliable picture of the PES around the $`n_{\mathrm{search}}`$ lowest-energy minima, where $`n_{\mathrm{search}}`$ is the number of minima from which transition state searches have been performed.

For LJ<sub>13</sub> there is only one MSB, reflecting its ideal funnel character and the remarkable degree of connectivity. For LJ<sub>19</sub> there are only a few MSB’s, and these are all directly connected to the primary MSB, again reflecting the single funnel character of the PES that we noted earlier. The double-funnel character of the LJ<sub>38</sub> PES is still apparent in the MSB disconnectivity graph, but now we get a better impression of the overall shape of the PES (Figure 6a). There is a wide, gently sloping funnel down from the higher-energy minima towards the low-energy icosahedral funnel, whereas the funnel down to the global minimum is much narrower. 2292 of the minima lie on monotonic sequences to the lowest-energy icosahedral minimum, whereas only 518 lie on sequences leading to the global minimum. Only 12 of the minima lie on sequences to both, showing that there is little overlap between the two MSB’s.

For LJ<sub>55</sub> the MSB analysis leads to a remarkable simplification of the disconnectivity graph (Figure 6b). The single-defect minima produce just one MSB, and the fine structure of the two-defect minima collapses onto the band of MSB’s that branch off at $`273ϵ`$ and $`272ϵ`$. The remaining branches are mostly three-defect minima with the three surface atoms close together. The graph clearly shows the single-funnel character of the PES.

### C Transforming the PES

One approach to global optimization is to transform the PES to a form for which it is easier to find the global minimum. One efficient approach of this type is the Monte Carlo minimization or basin-hopping algorithm. This unbiased method succeeded in finding the global minimum for all the clusters considered in this paper. In this approach the transformed potential energy $`\stackrel{~}{E}_c`$ is given by

$$\stackrel{~}{E}_c(𝐗)=\mathrm{min}\left\{E_c(𝐗)\right\}.$$
(2)   
Hence, the potential energy at any point in configuration space is assigned to that of the local minimum obtained by a minimization starting from that point, and the PES is mapped onto a set of interpenetrating staircases with steps corresponding to the basins of attraction surrounding a particular local minimum. Figure 8 illustrates the transformation for a simple one-dimensional example. Note that the transformation does not change the identity of the global minimum, nor the relative energies of any of the minima, but it does remove the downhill barriers between directly connected minima. This latter change has a significant effect on the dynamics and thermodynamics.

Due to our PES search strategy virtually all the unsearched higher-energy minima are directly connected to lower-energy minima. As the transformation removes the downhill barriers, in the disconnectivity graphs these unsearched minima would be directly connected to the stem associated with the superbasin which contains the relevant lower-energy minimum. This connectivity makes the higher-energy regions of the PES look funnel-like irrespective of their actual character. Therefore, in the graphs we only depict branches which lead to the lowest $`n_{\mathrm{search}}`$ steps.

The two examples shown in Figure 7 illustrate the effects of the staircase transformation on the disconnectivity graphs. It removes the barriers to progress down a funnel, and the disconnectivity graph for LJ<sub>19</sub> is therefore transformed into that for an ideal funnel. The long, dangling branches that are indicative of large barriers have disappeared, and so relaxation is now easier. However, barriers to interfunnel passage remain because the major component of such barriers usually arises from the high energy minima that the system has to pass through to go between funnels. (Figure 5a). For example, the transformation reduces the potential energy barriers for interfunnel passage by only $`0.68ϵ`$ for LJ<sub>38</sub> and $`0.86ϵ`$ for LJ<sub>75</sub>. Therefore, the splitting of the LJ<sub>38</sub> disconnectivity graph into two funnels is still clear, and perhaps even more obvious because many of the other barriers have been removed (Figure 7b). As for LJ<sub>19</sub>, the icosahedral region of configuration space now looks much more funnel-like.

The retention of the energy barriers for interfunnel passage in LJ<sub>38</sub> and LJ<sub>75</sub> means that the global minima of these clusters are still hard to find even on the transformed PES. In fact, the success of the basin-hopping algorithm for these clusters lies in the changes to the thermodynamics caused by the transformation. The thermodynamic transitions are broadened, thus providing a temperature window where both the states at the top of the paths for interfunnel passage and the states at the bottom of the two funnels have a significant probability of being occupied, making passage between funnels easier.

On the transformed PES any pathway that is monotonically decreasing in energy (albeit a stepped rather than a smooth one) must end at the step corresponding to the bottom of an MSB. Therefore, an MSB on the transformed PES is analogous to a basin of attraction surrounding a minimum on the original PES. Figure 8 illustrates this point: on the transformed PES the two MSB’s are like ‘stepped’ minima. However, the analogy breaks down in one respect. A basin of attraction surrounding a minimum is defined as the set of points from which steepest-descent paths lead down to that minimum. The use of the steepest-descent path ensures that each point in configuration space is mapped uniquely to a minimum even when a point is higher in energy than the lowest transition state connected to that minimum. However, on a step on the transformed PES the gradient is zero and so no steepest-descent direction is defined. Therefore, the mapping of a point in configuration space to a MSB bottom may not be unique if it lies above the divide between two MSB’s, and so MSB’s, unlike basins of attraction, can overlap.

In Figure 7c we show the disconnectivity graph for the transformed PES of LJ<sub>38</sub> with only the branches leading to the bottoms of MSB’s displayed. It has a similar structure to Figure 7b but the density of branches is greatly reduced.

## IV Conclusions

In this paper we have shown that disconnectivity graphs can provide a valuable tool for visual representation of a PES, and that they can provide clear physical insights into structure, dynamics and thermodynamics.

The exponential increase in the number of minima with system size means that disconnectivity graphs with branches representing individual minima must, as the size increases, inevitably focus on low-energy regions of the PES which represent a decreasing proportion of the whole configuration space. By taking monotonic sequences basins as the basic topographical unit, we have extended the cluster size for which a disconnectivity graph can represent the full sample of minima. However, these graphs still cannot provide a truly global picture of the PES because of the incompleteness of our samples of stationary points.

We have also used disconnectivity graphs to probe the transformed PES that is searched by the Monte Carlo minimization or ‘basin-hopping’ global optimization algorithm. The transformation removes downhill barriers, making relaxation down a funnel much easier. For example, the algorithm can find the global minimum of LJ<sub>55</sub>, which has a single-funnel PES, in on average fewer than 150 steps when started from a random geometry. However, the barriers between funnels remain, and so the success of basin-hopping for double-funnel PES’s must be explained in terms of the different thermodynamics for the transformed landscape.

###### 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 and M.A.M. to the Engineering and Physical Sciences Research Council for financial support. We would also like to thank Dr David Manolopoulos for his suggestion that disconnectivity graphs could be used to probe the transformed energy landscape used by the basin-hopping algorithm.