Describing a potential energy surface in terms of local minima and the transition states that connect them provides a conceptual and computational framework for understanding and predicting observable properties. Visualizing the potential energy landscape using disconnectivity graphs supplies a graphical connection between different structure-seeking systems, which can relax efficiently to a particular morphology. Landscapes involving competing morphologies support multiple potential energy funnels, which may exhibit characteristic heat capacity features and relaxation time scales. These connections between the organization of the potential energy landscape and structure, dynamics and thermodynamics are common to all the examples presented, ranging from atomic and molecular clusters to biomolecules and soft and condensed matter. Further connections between motifs in the energy landscape and the interparticle forces can be developed using symmetry considerations and results from catastrophe theory.
Where are the genes for Paulingite?
Title of a lecture by A. L. Mackay, Lund, 20 March 1984 .
The present contribution will address this famous question using the framework of energy landscapes . What we seek is an insightful way to understand and predict properties of molecular, mesoscopic and condensed matter systems that lies somewhere between the bare foundations of the Schrödinger equation and a direct investigation of every specific problem. The energy landscapes scheme provides one possible approach to this task, which has revealed connections between a remarkably diverse range of problems.
Our starting point is a potential energy surface having implicitly assumed that the Born–Oppenheimer approximation  can be employed. For a system containing N atoms, the potential energy is then a function of the 3N Cartesian coordinates, and the potential energy surface is a 3N+1-dimensional object, where the extra dimension corresponds to the value of the function. For a finite system in field-free space, internal coordinates can be used to eliminate overall translations and rotations, which reduces the dimension by six. The precise details of how the potential energy is actually calculated are generally system-dependent, and may also hinge upon the questions we are trying to answer. For example, if making or breaking covalent bonds is important then an explicit treatment of the electronic structure will probably be required. Accuracy is also essential if we aim to reproduce observable properties that are sensitively dependent on details of the energy landscape, such as tunnelling splittings [4,5]. However, empirical representations of the interparticle forces may be more appropriate if the aim is to understand trends in a range of different systems. Similarly, coarse-grained models are important for describing mesoscopic systems, particularly to shed light on how the formation of structures such as shells and spirals depends on the forces between building blocks [6,7].
It is certainly appealing to think about low-dimensional surfaces in terms of familiar geographical concepts, such as valleys and hills. This was the description adopted by James Clark Maxwell in his article for Philosophical Magazine ‘On Hills and Dales’ in 1870 . Unfortunately, 3N+1-dimensional surfaces usually cannot be visualized directly, and low-dimensional projections can be misleading. For example, projection onto one or two order parameters to obtain a free energy surface will not preserve barriers faithfully if the corresponding pathway does not involve a large enough change in the order parameter [9–16].
An alternative approach, which does preserve the lowest potential or free energy barriers between local minima, can be realized using disconnectivity graphs [17–20]. This framework is described in §3. To conclude this introduction it may be helpful to focus on how particular observable properties of an atomic, molecular or condensed matter system are encoded in the potential energy surface (PES), V (X). Here, X is the 3N-dimensional vector of atomic coordinates. If we think of a structure in terms of observable isomers, or distinct atomic frameworks with well-defined characteristic properties, then we can usually make a connection with local minima of the PES. These are stationary points of V (X), where the gradient vanishes, i.e. ∇V (X)=0, and any infinitesimal displacement of internal coordinates raises the energy. Thermodynamic properties can be defined locally for particular isomers, or as global ensemble averages. In either case, we require a density of states or a canonical partition function, which can be expressed in terms of V (X). The classical dynamical properties of the system are determined by forces derived from the gradient of V (X). Hence, we possess all the theoretical tools required, in principle, to extract the observables of interest that are associated with a particular PES. In fact, a more coarse-grained approach, described in §2, provides a way to unify our understanding of systems ranging from atomic and molecular clusters, through biomolecules, to soft and condensed matter. For example, the ability to self-organize reliably into a particular structure, or the emergence of glassy properties, depends on the organization of the potential energy landscape according to universal principles.
2. Potential energy landscapes
In §1, we have already associated stable isomers with local minima of the PES. Using the superposition approach [2,21–28], thermodynamic properties can be obtained from a partition function written as a sum over the basins of attraction of local minima j at temperature T: 2.1Here, the contribution from minimum j is obtained by restricting the corresponding integral to the basin of attraction for j [2,29,30], which is the region of configuration space for which steepest-decent paths lead to j. This approach is particularly efficient for cases of broken ergodicity , when alternative low-lying minima are separated by barriers that are large compared with kBT, where kB is Boltzmann's constant.
Knowledge of pathways between local minima can be used to improve approximations to the superposition formula in equation (2.1)  and to address dynamics. The rearrangements between local minima can be analysed by calculating transition states of the PES, which are here defined according to the geometrical criterion of Murrell & Laidler  as stationary points of V (X) that have a single imaginary normal mode frequency. Each transition state is associated with two downhill steepest-descent paths corresponding to positive and negative displacements along the reaction coordinate associated with the unique imaginary normal mode frequency. Following these paths enables us to define the connectivity of the PES in terms of the local minima connected by each transition state, and approximate rate constants can be evaluated for these transitions using statistical rate theory [33–36].
Analysing structure, thermodynamics and kinetics in terms of local minima and transition states of the PES provides a coarse-grained approach that serves both as a conceptual and a computational framework. In effect, we have reduced the exploration of the landscape to geometry optimization, namely the characterization of local minima and the pathways that connect them via transition states. All the applications presented below were made possible by efficient geometry optimization. Details of the minimization [37,38] and transition state searching [39–46] algorithms employed can be found elsewhere. Computer programs that implement these algorithms can be downloaded for use under the Gnu General Public Licence. This software includes tools for global optimization  via basin-hopping [48–50], location of transition states and characterization of pathways , and methods for expanding connected stationary point databases  using discrete path sampling [53,54], combined with kinetic analysis [55,56]. Here, we will focus on the insight that can be gained from visualizing the network [14,57] defined by a database of local minima and transition states, and on the emergence of characteristic thermodynamic and kinetic properties as a consequence of particular motifs.
3. Visualizing the landscape
(a) Disconnectivity graphs
A connected database of local minima and transition states not only defines a kinetic transition network [14,16,57], but also provides the information required to visualize the landscape using disconnectivity graphs . Constructing such graphs for a wide range of atomic and molecular clusters, biomolecules and condensed matter systems has enabled us to understand how observable properties are governed by the organization of the landscape [2,58,59].
To construct a potential energy disconnectivity graph from a connected set of local minima, we first choose an energy spacing, ΔV , and then determine how the minima are partitioned into subsets (superbasins ) at threshold energies V 0,V 0+ΔV,V 0+2ΔV,…. The subsets are groups of minima, which can interconvert via transition states that lie below the potential energy threshold. At sufficiently high energy, all the minima can interconvert without exceeding the threshold (unless there are infinite barriers) and there is a single superbasin. As we lower the threshold energy, the superbasins progressively split apart, and the disconnectivity graph follows this splitting by connecting subgroups to the parent superbasin at the threshold energy above. Hence, this scheme explicitly embraces hierarchical organization of the energy landscape, in keeping with the theme of Mackay . Eventually, the superbasins split into individual local minima, which are represented by points at the corresponding potential energy on a vertical scale, connected to the parent superbasin at the closest superbasin energy above.
The disconnectivity graph construction is perhaps best understood by considering some examples. A wide range of atomic and molecular clusters have been found to favour icosahedral packing, with complete icosahedra defined by Mackay's scheme  having special stability. The sizes at which such stability occurs are often referred to as ‘magic numbers’, since they exhibit enhanced abundances in molecular beam experiments. Disconnectivity graphs for the first two clusters in the series are illustrated in figures 1 and 2. These graphs were calculated for particles interacting via the isotropic, pairwise additive Lennard-Jones potential . In each case, the spacing in the horizontal direction was simply chosen to clarify the structure as much as possible, and highlight the special position of the global potential energy minimum. When permutation–inversion isomers of the global minimum are collected together , the Mackay icosahedra appear well separated in energy from the next-lowest structures, which are icosahedra with single vacancies in the outer shell and capping atoms on the surface . Disconnectivity graphs are technically tree graphs , and the organization of the landscape for LJ13 and LJ55 has been informally described as a ‘palm tree’ in previous work . This is the form that we associate with efficient relaxation to the global minimum, since there is a strong energetic driving force towards this structure, and relatively low downhill barriers [18,58,59]. Alternatively, there is a well-defined free energy minimum over a wide range of temperature, which is kinetically accessible under the same conditions. Hence, we can identify such landscapes with ‘funnelling’ properties [65–69] and convergent kinetic pathways . The palm tree structure is also associated with ‘minimal frustration’ [66,71], since there are no low-lying competing morphologies separated from the global minimum by high barriers. Geometrical and visual analysis of potential energy landscapes has enabled us to generalize these ideas from the field of protein folding and unify our understanding of efficient structure-seeking systems throughout molecular science. Some further examples are presented in §3b, and contrasted with frustrated and glassy landscapes in §3c.
Mackay icosahedra are also expected to appear as global minima for atomic clusters bound by long- and medium-range isotropic potentials for larger sizes. However, the strain  involved in packing 20 distorted tetrahedra will eventually lead to competition from structures where the nearest-neighbour distances are closer to the optimal separation. The size at which this crossover might occur has been investigated in a number of studies [73–85]. There are also indications that competition from alternative structures may affect the observed thermodynamic properties. Although LJ147 was found to exhibit fairly straightforward melting behaviour [86,87], describable as the finite analogue of a bulk phase transition, the behaviour of LJ309 is more complicated .
(b) A unified view of self-organizing systems
The palm tree potential energy landscape described above for LJ13 and LJ55 is associated with efficient relaxation to a well-defined morphology in a wide variety of different systems. For example, figure 3 illustrates a free energy disconnectivity graph [19,20] for a peptide containing 16 amino acids, which forms a β-hairpin in solution [90–93], and in the larger B1 domain of protein G . At this temperature, β-hairpin structures correspond to a well-defined free energy minimum, and folding is expected to occur in a microsecond or less [89,95–100].
For C60, rearrangements between different fullerene isomers are associated with relatively large barriers [18,101–103]. However, there is still a strong energetic bias towards the icosahedral buckminsterfullerene structure. The corresponding disconnectivity graph has longer branches because of the high downhill barriers (figure 4), and looks more like a willow tree [18,101,103]. Efficient relaxation to the global minimum now requires a relatively high temperature , as expected from experiment [104,105].
Recent work has explored the energy landscapes of mesoscopic systems, exploiting an angle-axis coordinate system to describe the interactions between rigid building blocks [6,7,58,106,107]. A generic anisotropic potential , combined with two axial repulsive sites, has proved to be a remarkably versatile representation, supporting spheroidal shells, tubular, helical and even head–tail structures for suitable parametrizations [6,7]. Figure 5 shows the disconnectivity graph obtained for 32 particles, which corresponds to a triangulation number  of T=3 for a virus capsid shell. The global minimum has icosahedral symmetry, and the structure corresponds to the face dual of the truncated icosahedral geometry for buckminsterfullerene in figure 4. Efficient self-assembly of an icosahedral shell is predicted for this system over a wide range of temperature, starting from any high-energy 32-particle structure .
(c) Competing morphologies and frustration
Energy landscapes that correspond to efficient relaxation into a well-defined morphology are greatly outnumbered by landscapes that exhibit multiple low-lying competing morphologies. Figure 6 shows the disconnectivity graph associated with the LJ38 cluster, which exhibits two main features corresponding to incomplete icosahedral and truncated octahedral minima, separated by a relatively high barrier . This landscape provides an interesting challenge for global optimization algorithms as well as sampling schemes designed to address thermodynamics and kinetics [2,49,50,53,54,110–114]. The global potential energy minimum is the truncated octahedron, but the configuration space corresponding to incomplete icosahedra is favoured by entropy. Hence, this cluster exhibits the finite system analogue of a first-order phase transition, with an associated feature in the heat capacity. The global minimum is harder to locate than for clusters of comparable size where the lowest minimum is based on icosahedral packing, because the icosahedral minima act as a kinetic trap for LJ38. The trap results in two distinct time scales for relaxation to the global minimum [2,115]. The fastest process corresponds to direct relaxation to the truncated octahedron, and the slower process corresponds to trajectories that visit low-lying icosahedra first, and must then escape. These time scales are evident in figure 7, where plateaux occur in the occupation probability at an intermediate time scale (around 109 reduced time units) and for true equilibrium beyond about 1013 reduced time units. Such features are directly related to the problem of ‘broken ergodicity’, where the time scale for the system to reach global equilibrium becomes comparable to the observation or simulation time scale.
For systems with more numerous competing morphologies, separate heat capacity features and relaxation time scales may be harder to resolve. A wide range of examples for such multi-funnel landscapes have been identified in previous work, and results for a structural glass-former are illustrated in the following section. One further example will be described first, namely the (H2O)20 cluster illustrated in figure 8. For this system, well-defined regions of the graph can be disconnected by cutting single edges . Hence, this landscape has distinct signatures of hierarchical organization, and the overall appearance has been likened to a banyan tree, where new trunks are formed by branches that dangle to ground level. The global minimum for this cluster is now relatively difficult to locate, and the structure of the landscape resembles disconnectivity graphs constructed for fractal topologies . The greater complexity of this landscape is associated with barriers that appear on different energy scales. For the (H2O)20 cluster, this range of barrier heights results from the interplay between translational and orientational degrees of freedom for the individual water molecules. An equivalent view is that minima with different hydrogen-bonding patterns, but essentially the same oxygen framework, have energies that span a very wide range of energy. Hence it is possible to identify structures that have practically the same oxygen arrangement as the global minimum, but lie much higher in energy.
(d) Contrasting landscapes for glassy and crystalline behaviour
To conclude this brief overview of visualizations and phenomenology within the energy landscape approach, we contrast two bulk systems modelled using periodic boundary conditions. Figure 9 shows a disconnectivity graph for low-lying minima associated with the diamond cubic crystal for silicon. This landscape corresponds to a modified Stillinger–Weber (SW) silicon potential  under constant volume conditions, where the supercell was optimized for the crystal . The global minimum is readily located for this system, and there are no low-lying competing morphologies. The lowest defect structures are clearly separated from the global minimum, both for the empirical SW potential and when the electronic structure is accounted for explicitly [43,44].
A very different picture is obtained for the binary Lennard-Jones system in figure 10. This binary representation was first introduced to model the metallic glass Ni0.8P0.2 , and it has proved popular for simulations because it does not crystallize readily on a molecular dynamics time scale [118,121–129]. In fact it has been claimed that ‘This model is suitable for such problems because of the lack of a crystalline state’, although crystalline phases have subsequently been characterized [130,131]. The mixture corresponding to the graph in figure 10 contains two different types of atom, A and B, with an A : B ratio of 80 : 20. The Lennard-Jones parameters are σAA=1, σAB=0.8, σBB=0.88, ϵAA=1, ϵAB=1.5 and ϵBB=0.5 . The graph was constructed from minima sampled from the supercooled region of configuration space over a locally ergodic time interval for a supercell containing 60 particles at a number density of 1.3 and kBT/ϵAA=0.96 . At temperatures around the glass transition for this system, the barriers between different amorphous structures in distinct funnels are typically of order 30kBT. The disconnectivity graph in this region of configuration space is therefore highly frustrated [66,71], and provides a useful reference for a ‘rough’ or ‘rugged’ landscape. This picture has been confirmed in subsequent calculations [132,133], and analysis of rearrangements in terms of cage-breaking events has revealed  a further level of hierarchical structure and a possible definition of ‘metabasins’ [134,135].
4. Organization of the energy landscape
The above examples illustrate how observables such as structure, dynamics and thermodynamics for a particular system are determined by the organization of the energy landscape. Extracting numerical values for quantities such as free energy barriers and rate constants generally requires additional calculations. However, the appearance of alternative funnels and competing morphologies immediately suggests the existence of multiple time scales for relaxation and possible features in the heat capacity. At this level of description, we can therefore explain how various properties of interest are encoded in the energy landscape.
An equally important problem concerns the relation between the energy landscape and the underlying interatomic or intermolecular forces . Here, we cannot expect to find a simple connection between the interparticle interactions and properties such as the ability to self-organize. Since clusters of different sizes composed of identical atoms can themselves exhibit a wide range of thermodynamic and dynamical behaviour, we know that such properties are not a universal function of the interparticle potential. However, there are some additional theoretical tools that may prove to be useful in uncovering further general principles, namely catastrophe theory and symmetry.
(a) Families of landscapes: catastrophe theory
As noted above, the organization of a single potential energy surface can be usefully described in terms of its stationary points. The organization of families of potential energy surfaces as a function of parameters in the interparticle forces is determined by stationary points that have additional zero eigenvalues of the Hessian (second derivative) matrix , beyond those required by conservation laws. Catastrophe theory  provides canonical forms for the expansion of the potential around such non-Morse points. These predictions have been tested for atomic clusters as a function of the range of the interatomic potential , and have also been applied to analyse first-order phase transitions .
For families of landscapes related by variation of a single parameter, such as the strength of an applied field or the range of interparticle forces, the geometry is expected to be described by a fold catastrophe. In the limit of short pathways, where the system is close to a non-Morse point, a universal connection between the barrier height, ΔV , the path length, Δs, and the smallest non-zero Hessian eigenvalue, λ, is predicted, namely  4.1Alternatively, for a pathway with twofold symmetry, the geometry for a single control parameter corresponds to a symmetrical cusp catastrophe and the predicted relationship is 4.2where λ is the smallest non-zero Hessian eigenvalue of the minimum, and −2λ is the smallest non-zero Hessian eigenvalue of the corresponding transition state .
The connection between path lengths, barrier heights and curvatures has been used to explain the locally rougher but globally flatter appearance observed for landscapes corresponding to short-ranged potentials . It also provides a basis for Hammond's postulate, that the transition state most closely resembles the higher in energy of the products and reactants . Trends such as these, which are realized across a range of different systems, may therefore be a consequence of the universality embodied by catastrophe theory.
The appearance of symmetry has often been connected with notions of structural stability. Pierre Curie suggested that ‘the symmetry characteristic of a phenomenon is the maximal symmetry compatible with the existence of the phenomenon' , while D'Arcy Thompson commented that ‘the perfection of mathematical beauty is such… that whatsoever is most beautiful and regular is also found to be most useful and excellent’ . Looking through the point groups of global minima deposited in the Cambridge Cluster Database  certainly suggests that higher symmetry is prevalent among these structures. However, it is not difficult to think of counter-examples, where a given system does not exhibit the highest symmetry imaginable.
While it can be proved that certain high-symmetry ‘balanced’ geometries are stationary points of some index, a rigorous connection to energetics is lacking. However, it is possible to derive a connection that is expected to hold in an average sense [2,144,145]. Here, we consider the total energy as a sum over terms from a many-body expansion, involving single atom, pairwise and three-body terms, etc. For a system with geometrical symmetries, some of these terms must appear in degenerate sets. Approximate symmetry can be allowed for by considering energy bins with a finite width. Suppose that all the energetic contributions are drawn from the same distribution, but contributions that are equivalent owing to symmetry are only drawn once and replicated appropriately. The mean value for the total energy would be the same, but the variance will be larger when there are duplicate contributions. For example, if every term appeared twice, the variance would be doubled. Hence structures with more symmetry seem more likely to have either particularly low or particularly high energy. High-symmetry structures are therefore likely to appear at both ends of the energy range for local minima of a given composition. In fact, the contributions to the total energy are not independent random variables. Nevertheless, there is a significant body of numerical evidence to support this ‘principle of maximum symmetry’ , as well as connections to designability  and minimal frustration [147,148].
Analysis of potential energy surfaces in terms of their stationary points, particularly local minima and transition states, provides a powerful framework for addressing the properties of clusters, biomolecules and soft and condensed matter. For example, visualizing the connectivity of stationary point databases, which constitute kinetic transition networks [14,16,57], enables us to connect the self-organizing properties of a diverse range of systems. In this sense, disconnectivity graphs [17,18] provide a distinct ‘phenotype’ for a zeroth-order prediction of observable properties, such as heat capacity features, and competing time scales for relaxation. The associated ‘genes’ could then correspond to the interparticle forces, as well as the number of particles of each type and external parameters, such as temperature and pressure. These forces in turn depend on the underlying Schrödinger equation or the classical Hamiltonian for the system in question. At each level of description, we can infer more about the likely observable properties. However, to make connections between generic behaviour, we need to identify common features before reaching the full detail of a complete numerical analysis of individual systems. The energy landscape approach outlined above represents one possible way of addressing this issue. It has provided insight into the behaviour of a wide range of different problems in molecular and condensed matter science. It has also enabled us to begin the closely related task of understanding how distinct motifs in the energy landscape are determined by the underlying interparticle forces.
I am very grateful to everyone who has contributed to this research effort, and to the EPSRC and BBSRC for funding.
One contribution of 14 to a Theme Issue ‘Beyond crystals: the dialectic of materials and information’.
- This journal is © 2012 The Royal Society