John W. R.
Morgan
^{a},
Dhagash
Mehta
^{b} and
David J.
Wales
*^{c}
^{a}Department of Chemical Engineering, University of Michigan, Ann Arbor, MI 48109-2136, USA. E-mail: jwrm@umich.edu
^{b}Department of Applied and Computational Mathematics and Statistics, University of Notre Dame, Notre Dame, IN 46556, USA. E-mail: dmehta@nd.edu
^{c}University Chemical Laboratories, Lensfield Road, Cambridge CB2 1EW, UK. E-mail: dw34@cam.ac.uk
First published on 7th September 2017
A database of minima and transition states corresponds to a network where the minima represent nodes and the transition states correspond to edges between the pairs of minima they connect via steepest-descent paths. Here we construct networks for small clusters bound by the Morse potential for a selection of physically relevant parameters, in two and three dimensions. The properties of these unweighted and undirected networks are analysed to examine two features: whether they are small-world, where the shortest path between nodes involves only a small number or edges; and whether they are scale-free, having a degree distribution that follows a power law. Small-world character is present, but statistical tests show that a power law is not a good fit, so the networks are not scale-free. These results for clusters are compared with the corresponding properties for the molecular and atomic structural glass formers ortho-terphenyl and binary Lennard-Jones. These glassy systems do not show small-world properties, suggesting that such behaviour is linked to the structure-seeking landscapes of the Morse clusters.
The degree of a node is the number of edges connected to it. The degree distribution is then defined as the number of nodes with given degrees.^{5} A path is a sequence of nodes connected by edges, with the length of the path being the number of edges it contains. The shortest path length gives the minimum number of edges between a pair of nodes, and the average shortest path length is the shortest path length averaged over all pairs of nodes. The clustering coefficient is the ratio of the number of connections between neighbours of a node to the number of such connections that could exist.^{6}
Two useful references are provided by models based on a lattice graph and a random graph. Random graphs have a small clustering coefficient and a slowly growing average shortest path length.^{5,7} In contrast, lattice graphs, which are composed of a regular array of nodes with edges only between nearest-neighbours in space, have a comparatively large average shortest path length and a larger clustering coefficient.^{6} Small-world networks were introduced by Watts and Strogatz^{6} and defined as networks that show a high degree of local clustering behaviour, similar to a lattice graph, but also a short path length, even between distant nodes, as exhibited by random graphs. The name comes from experiments performed by Milgram with letter passing, which demonstrated that most citizens of the USA are separated from each other by a surprisingly small number of social contacts.^{8} This result is now known in terms of “six degrees of separation” from the estimated average path length.^{9} Watts and Strogatz showed that networks from a wide variety of areas show small-world behaviour, including the neural network of C. elegans, the power grid of the Western USA, and the collaborations of film actors (linked to the concept of the Bacon number).^{10}
A scale-free network has a degree distribution with a power law tail, as defined by Barabási and Albert.^{11,12} This distribution implies there are a small number of nodes with a very high degree, called hubs, and more nodes with a smaller degree, acting as local hubs, down to nodes with only a few connections, in a hierarchical fashion. Power laws were fitted for a variety of available networks, including the World Wide Web and the film actor collaboration graph, and a preferential attachment model was suggested to explain how the behaviour arises. Following this work, power laws were fitted to many other networks, such as the interactions between solar magnetic loops leading to a solar flare,^{13} and aftershocks of earthquakes.^{14} However, few of these studies used robust statistical methods to determine whether or not a power law was a good fit to the data, and many of the results have been called into question by Clauset et al.^{15}
Doye and Massen^{16,17} studied Lennard-Jones (LJ) clusters and concluded that the networks were both small-world and scale-free. Small-world behaviour is believed to have implications for self-assembly, as it suggests that the global minimum can be reached from anywhere on the PES by a relatively short transition pathway. One would therefore expect that systems with a funnel-like PES, characteristic of structure-seekers, may display small-world behaviour.^{18} Scale-free behaviour was more surprising, as there is no obvious analogy to the addition of new nodes in the preferential attachment model, the most straightforward route to a power law distribution. However, noticing a correlation between minima with a high degree and a low potential energy, it was suggested that due to the larger basins of attraction for low energy minima, any other minimum was more likely to be connected to a low energy minimum than a high energy one, thus establishing an analogy to preferential attachment. This possibility was further investigated using Apollonian packings,^{19} concluding that the contacts between discs in the 2D Apollonian packing form a scale-free network with a spatial distribution that may resemble the catchment basins of a PES.
Other energy landscapes have been studied with a view to analysing the same properties. Protein folding networks have been considered by Rao and Caflisch.^{20} They used molecular dynamics to generate snapshots of the structure and then derived a conformation based on the secondary structure each residue belonged to. These snapshots do not precisely correspond to minima on the PES, so the formation of the network is not equivalent to the LJ results considered above, but the idea is similar. A power law was fitted to the degree distribution with the conclusion was that the network was scale-free, drawing attention to the hubs present in the network.
Bowman and Pande studied a Markov state model for another protein folding landscape.^{21} They noted that native states were hubs and suggested that this was characteristic of a scale-free network, but without consideration of the degree distribution. Similarly, Chakraborty et al. considered an RNA folding network, noting that it may be scale-free due to the presence of hubs, but without providing the degree distribution.^{22} It is, however, possible to construct networks with hubs that are not scale-free, the simplest example being the star network,^{5} in which there is one central node to which all others are attached and there no other edges. Therefore, although scale-free networks do contain hubs, the presence of hubs is not in itself enough evidence to conclude that a network is scale-free. Recently, Mehta et al.^{23} constructed stationary point databases for the Thomson problem^{24} up to 150 particles. They demonstrated that the networks show small-world behaviour, and that the presence of hubs suggested the possibility of the networks being scale-free, but acknowledged that further statistical testing was required. Recently, an analysis of networks of minima for machine learning problems has also been initiated.^{25}
Potential energy surfaces for glassy systems have a large number of competing low-lying amorphous minima with similar energies, separated by high barriers.^{1} Here we consider two examples: bulk ortho-terphenyl (OTP), which is often represented by a course-grained model consisting of a LJ site at the centre of each benzene ring;^{26} and bulk binary Lennard-Jones (BLJ), in which there are two different types of particles, interacting via LJ potentials with different parameters.^{27,28} If the small-world behaviour observed for atomic clusters is the result of a funnelled landscape, then we anticipate that glassy systems will be different.
The Morse potential^{29} is a widely used pairwise representation for modelling atomic clusters.^{30–37} The functional form is
V_{M}(R) = εe^{ρ(1−R/Re)}(e^{ρ(1−R/Re)} − 2), | (1) |
The Morse potential is pairwise additive so the total energy of the system is a sum over the interaction between all pairs of particles. This approximation neglects many-body interactions, which can be important in some systems, especially when the particle interaction range is large.^{44} The pair approximation has been successfully used to predict experimental properties, in colloids for example, when the range is small.^{39}
In the present work we explore the PES for small atomic clusters bound by the Morse potential with a range of values for ρ, namely 3, 6, 10, 14 and 30. Our databases are almost complete, meaning that they were expanded until no new minima were found, but the search methods used do not rigorously guarantee that all the minima have been located. Disconnectivity graphs are used to visualise the potential energy landscape.^{45,46} Previously, global minima and some low-lying minima have been described,^{31–37,42,47} but attempts to generate relatively complete databases have been made only for a few specific cases.^{30,43}
For each database, a network was constructed by considering the minima as nodes with edges between any two minima directly connected by a transition state. The barrier height for the transition state was ignored and the path through the transition state can be followed in either direction, so the network was treated as undirected and unweighted. The network was analysed using the Python NetworkX package.^{68} A further check on the completeness of the database was to ensure that all minima were connected to each other by one or more transition states, i.e. that the network was connected. Finally, transition states connecting the same pair of minima as an earlier pathway were removed, to leave at most one transition state connecting any pair. Transition states that connected permutation-inversion isomers of the same minimum, which correspond to degenerate rearrangements,^{69} were also removed.
We chose to study clusters with ρ = 3, 6, 10, 14 and 30, in both two and three dimensions. These range parameters correspond to previous values used for the Morse potential and give a good spread of physically relevant pair potentials. The number of particles was chosen to produce databases with between several hundred and a few thousand minima. Selecting these sizes meant we could be confident of generating essentially complete databases in a reasonable time, improving the reliability of statistical tests. The symbol M^{xD}_{N} is used to refer to a Morse cluster in x dimensions with N particles.
We are interested in whether the degree distribution follows a power law,^{11} such that
p_{k} = Ck^{−α}, |
lnp_{k} = lnC − αlnk, | (2) |
(3) |
Typically, a power law is only obeyed in the tail of the distribution, so a cutoff k_{min} above which the power law applies needs to be determined. Clauset et al. also have a method for determining the best value of k_{min}.^{71} A fit is calculated for all possible values of k_{min} and the value which gives the smallest Kolmogorov–Smirnov (KS) statistic^{72} is chosen.
A goodness-of-fit test is used to check whether the proposed power law is a plausible fit for the measured data.^{15} Data sets are generated from the power law and compared to the measured data. The p value is the fraction of generated data sets that are a worse fit than the measured data. The fit is rejected if the p value is less than 0.05. If p > 0.05, it does not necessarily indicate that the power law is a good fit, merely that we do not have sufficient evidence to reject it.
Other distributions must be considered that could be a better fit to the data, such as the log-normal distribution. If we cannot say that a power law is a better fit than the log-normal, we should not conclude that the data follows a power law.^{15} A likelihood ratio test is used to compare the fit of two distributions, indicating which distribution is more likely, and a p value is used to judge whether the result is significant.^{73} A small p value (<0.05) means there is sufficient evidence to accept the result of the tests.
(4) |
The local clustering coefficient can be extended to a global clustering coefficient for the network in two ways. The average clustering coefficient is the mean of the local clustering coefficients:^{6}
(5) |
(6) |
These two measures differ in the weighting they give to nodes of small degree. Such nodes have a small number of possible neighbour pairs, so they have a small denominator in their local clustering coefficient. They therefore dominate the sum in eqn (5). In networks with a large number of small degree nodes, the average clustering coefficient does not depend much on the larger degree nodes. In PES networks, the large degree nodes tend to be those with a low potential energy, which are considered more important for the properties of the network. We therefore prefer the transitivity to the average clustering coefficient.
To determine whether or not the transitivity should be considered large or small, it is useful to compare it to the value for an Erdös-Rényi random network^{75} with the same number of nodes and edges. In this model, every edge has an equal probability of being present, independent of the existence of any other edge. Therefore the probability of a given triad being a triangle, which is the probability of the final edge being present, is the probability of any edge existing. The maximum number of edges is the number of possible pairs of nodes, n(n − 1)/2, so given the number of edges (m) and nodes (n), the transitivity is:
(7) |
(8) |
Dimensions | ρ | N | Minima | Transition states |
---|---|---|---|---|
2 | 3 | 27 | 1135 | 17006 |
2 | 6 | 14 | 805 | 13098 |
2 | 10 | 14 | 840 | 19965 |
2 | 14 | 14 | 843 | 20552 |
2 | 30 | 13 | 358 | 6966 |
3 | 3 | 30 | 663 | 6678 |
3 | 6 | 13 | 1478 | 25173 |
3 | 10 | 12 | 2258 | 41441 |
3 | 14 | 12 | 2980 | 60615 |
3 | 30 | 11 | 1127 | 13752 |
Fig. 1 The number or minima (black) and transition states (red) for all the two-dimensional (top) and three-dimensional (bottom) clusters considered. |
Some general patterns are clear: 3D clusters have more minima for the same range and number of particles than 2D clusters, except at very long-range with ρ = 3. The well-known approximate exponential growth of the number of minima and transition states with system size is present.^{4,79} A shorter-range potential produces an energy landscape that is locally rough, having more minima, but globally smooth with a less overall funnelled structure.^{42} To avoid the databases growing too large to be confident of finding the complete network, it was necessary to restrict the number of particles in the cluster as the range decreased. The difference in the database size is more apparent in three dimensions than two dimensions. This effect is evident in the disconnectivity graphs in Fig. 2, where the funnel for a ρ = 6 landscape is apparent, suggesting structure-seeking behaviour, but there is no such structure for ρ = 30, where the landscape is more frustrated in the region of the global minimum.
Fig. 2 Disconnectivity graphs for M^{3D}_{13} with ρ = 6 (left) and M^{3D}_{11} with ρ = 30 (right). The difference in the global structure of the landscape is clear. |
Fig. 3 Plots of the average shortest path length against the number of minima for all the 2D clusters (top) and all the 3D clusters (bottom). Note the logarithmic scale on the horizontal axis. |
The transitivity for a small-world network is expected to be larger than for the equivalent random graph with the same number of minima and edges.^{6} The transitivity scaled by the value for a random graph with the same number of nodes and edges is plotted for all the clusters in Fig. 4. The expected small-world behaviour is again apparent, with the transitivity larger than for a random graph in all the clusters, except for three-dimensional systems with very few minima, where finite size effects dominate. There is no significant difference between the values for different ranges, except for ρ = 3, which appears to have a slightly higher transitivity. The long range of the potential produces a smoother landscape, which makes it more likely that minima in distant parts of configuration space are directly connected. Short-range clusters do not display the structure-seeking properties of long-range clusters, since their landscapes are not funnelled, and there are relatively large barriers separating some minima with similar energies. Since the present analysis accounts only for the existence of connections, regardless of barrier heights, we find that the small-world properties of the network are similar for short- and long-range clusters.
Insight into the structure of the landscape and network can be obtained by inspecting a plot of the degree of the nodes against their potential energy, shown in Fig. 5 for M^{3D}_{13} with ρ = 6 and the M^{3D}_{11} with ρ = 30. The funnelled landscape for the ρ = 6 cluster has one node with a significantly lower potential energy, as seen in the disconnectivity graph (Fig. 2), which is also the node with the highest degree. Over a third of the minima in the network are directly connected to the global minimum. At higher energies, although there are many minima with different degrees at similar energies, the maximum degree for a given energy decreases as the energy increases. These observations can be related to funnels on the landscape: the global minimum has a large basin of attraction, which occupies a significant volume in configuration space, giving a large boundary surface containing many transition states. High energy minima have small basins of attraction with a small boundary surface and fewer transition states. In comparison, the global minimum for M^{3D}_{11} with ρ = 30 is not the minimum with the highest degree, and the most highly connected minimum has transition states linking it to less than a quarter of the other minima. The general trend for higher energy minima to have a lower degree is still present. The prominent vertical stripes are a result of the short-ranged potential. Since the potential energy is primarily determined by the number of nearest-neighbour contacts, an integer, it is approximately discretised.^{80–84}
Fig. 5 Degree of each minimum as a function of potential energy for M^{3D}_{13} with ρ = 6 (top) and M^{3D}_{11} with ρ = 30 (bottom). |
Dimensions | ρ | N | α | k _{min} | n _{tail} | p |
---|---|---|---|---|---|---|
2 | 3 | 27 | 2.32 ± 0.32 | 23 ± 12 | 255 | 0.00 |
2 | 6 | 14 | 3.50 ± 0.16 | 26 ± 4 | 302 | 0.00 |
2 | 10 | 14 | 2.55 ± 0.14 | 18 ± 7 | 765 | 0.00 |
2 | 14 | 14 | 2.03 ± 0.21 | 10 ± 5 | 840 | 0.00 |
2 | 30 | 13 | 2.90 ± 0.15 | 17 ± 4 | 332 | 0.00 |
3 | 3 | 30 | 2.79 ± 0.17 | 10 ± 2 | 219 | 0.18 |
3 | 6 | 13 | 2.66 ± 0.16 | 19 ± 8 | 698 | 0.02 |
3 | 10 | 12 | 2.56 ± 0.09 | 28 ± 3 | 498 | 0.03 |
3 | 14 | 12 | 2.37 ± 0.12 | 65 ± 11 | 174 | 0.02 |
3 | 30 | 11 | 3.30 ± 0.25 | 48 ± 8 | 47 | 0.16 |
3 | LJ | 14 | 2.82 ± 0.10 | 37 ± 8 | 513 | 0.07 |
The cumulative degree distribution is plotted in Fig. 6 for M^{3D}_{13} clusters with ρ = 6 (top) and M^{3D}_{11} with ρ = 30 (bottom), along with the best power law fit calculated by the MLE method. The low degree end does not follow a power law, because the number of minima for a cluster is finite. Above k_{min}, the distribution appears to be approximately linear. However, the goodness-of-fit test for the ρ = 6 cluster gives p < 0.05, indicating the power law is not a good fit for the data.
Fig. 6 Plots of the normalised cumulative degree distribution for M^{3D}_{13} with ρ = 6 (top) and M^{3D}_{11} with ρ = 30 (bottom). The calculated best fit is indicated by the black dashed line. |
In fact, according to the suggested criteria, only two clusters have a p value compatible of a power law fit: M^{3D}_{30} with ρ = 3 and M^{3D}_{11} with ρ = 30. Both of these clusters have a large value of k_{min} compared to the number of minima in the cluster, giving a small value of n_{tail}. Therefore the fit should be considered unreliable, especially as it is clearly the case that these clusters do not generally display scale-free behaviour.
Testing whether the incomplete OTP database is likely to be representative of the full network was assessed by taking increasingly large subsets of the database and looking at the convergence of the network properties. Here the appropriate network properties are the ones that are global, but do not directly depend on the size of the network. The average shortest path length is therefore unsuitable, as it is expected to increase with the network size. Three properties were chosen: the transitivity, the average clustering coefficient (the mean of the local clustering coefficients) and the assortativity.
All three properties are only defined for a connected network, but an arbitrary subset of the database will not generally be connected because of the way it was constructed. Hence transition states were added to the network one by one, in the order located by the original search, and the largest connected component selected. Whenever the largest connected component increased in size, the network properties were calculated. Graphs of the three properties against the number of transition states in the connected component are shown in Fig. 8. There are fluctuations in the values, but all three seem to be approaching a limit. However, caution is needed in interpreting these results. Sampling biases in generating the database could potentially lead to limiting values of these properties without truly being a representative sample of the network.
To test for small-world properties, the values of the average shortest path length and transitivity can be compared to those of a random graph. These ratios are plotted in Fig. 9 for increasingly large connected subgraphs of the network. The transitivity and average shortest path length increase more quickly than the values for the equivalent random graph suggesting that the network (or at least this subset of the full PES) exhibits lattice graph properties, rather than small-world behaviour. The network is mildly assortative, although the value is close to zero, indicating a slight preference for minima to connect to other minima with similar degree. Assortative networks tend to have cores of highly connected nodes with a periphery of low degree nodes.^{5} This pattern could arise in a glassy landscape, with well connected groups of structurally closely related minima, and higher energy minima forming the periphery between the numerous groups.
Scale-free networks are small-world, so this network is also not scale-free. We found p = 0.0 for the power law goodness-of-fit test.
The database was expanded until it contained 703827 minima and 360836 transition states. However, the largest connected component consisted of only 12905 minima and 14363 transition states.
Our analysis followed the steps described above for OTP. The convergence of the average clustering coefficient, transitivity and assortativity are shown in Fig. 10, and the ratios of the average shortest path length and transitivity to those for the equivalent random graph are shown in Fig. 11. The results are qualitatively similar to OTP, although the convergence of the average clustering coefficient and transitivity is less clear due to the smaller sample size. The value of the assortativity is negative, indicating slight disassortativity compared to the slight assortativity of the OTP system. Park and Newman suggest that disassortativity can arise due to the exclusion of alternative and self-connections,^{88} so this result is unsurprising, although it is unclear why the opposite behaviour to OTP is observed. The growing ratios of average shortest path length and transitivity to the equivalent random graph suggest the same conclusion: that the system is showing lattice graph behaviour rather than small-world behaviour.
The Morse networks were compared with glassy systems, which have very different energy landscapes characterised by numerous low energy amorphous structures separated by relatively high barriers. The two glassy networks did not show small-world behaviour, which suggests that the small-world properties of the LJ and Morse networks are a feature of their single-funnel potential energy landscapes. However, further investigation is required. Larger clusters with more minima and transition states will make statistical tests for scale-free behaviour more reliable.
Using the convergence of network properties to test whether a sample is representative opens the possibility of studying systems where it is impractical to find the complete network. While it remains computationally intractable to find the complete network for glassy systems, as the number of minima and transition states is far too large, we believe the approach here will allow further insight for such systems. It is certainly possible to generate and analyse larger samples, and if databases generated by different schemes converge on the same properties we can have greater confidence about whether the samples are truly representative. Analysing the network properties of a representative sample will provide deeper understanding of the behaviour of bulk glasses. The LJ_{38} cluster, which has a double-funnel landscape,^{89} could produce interesting results if the two funnels are considered separately and together. Studying a fast-folding protein network would also be insightful. On the basis of the results here, we anticipate that the folding network is not really scale-free, but due to its known structure-seeking nature, small-world character is likely.
This journal is © the Owner Societies 2017 |