## New Research In

### Physical Sciences

### Social Sciences

#### Featured Portals

#### Articles by Topic

### Biological Sciences

#### Featured Portals

#### Articles by Topic

- Agricultural Sciences
- Anthropology
- Applied Biological Sciences
- Biochemistry
- Biophysics and Computational Biology
- Cell Biology
- Developmental Biology
- Ecology
- Environmental Sciences
- Evolution
- Genetics
- Immunology and Inflammation
- Medical Sciences
- Microbiology
- Neuroscience
- Pharmacology
- Physiology
- Plant Biology
- Population Biology
- Psychological and Cognitive Sciences
- Sustainability Science
- Systems Biology

# A geometrical approach to computing free-energy landscapes from short-ranged potentials

Edited by José N. Onuchic, Rice University, Houston, TX, and approved November 7, 2012 (received for review July 10, 2012)

## Abstract

Particles interacting with short-ranged potentials have attracted increasing interest, partly for their ability to model mesoscale systems such as colloids interacting via DNA or depletion. We consider the free-energy landscape of such systems as the range of the potential goes to zero. In this limit, the landscape is entirely defined by geometrical manifolds, plus a single control parameter. These manifolds are fundamental objects that do not depend on the details of the interaction potential and provide the starting point from which any quantity characterizing the system—equilibrium or nonequilibrium—can be computed for arbitrary potentials. To consider dynamical quantities we compute the asymptotic limit of the Fokker–Planck equation and show that it becomes restricted to the low-dimensional manifolds connected by “sticky” boundary conditions. To illustrate our theory, we compute the low-dimensional manifolds for identical particles, providing a complete description of the lowest-energy parts of the landscape including floppy modes with up to 2 internal degrees of freedom. The results can be directly tested on colloidal clusters. This limit is a unique approach for understanding energy landscapes, and our hope is that it can also provide insight into finite-range potentials.

The dynamics on free-energy landscapes are a ubiquitous paradigm for characterizing molecular and mesoscopic systems, from atomic clusters, to protein folding, to colloidal clusters (1⇓⇓–4). The predominant strategy for understanding the dynamics on an energy landscape has focused on the stationary points of the energy, the local minima and the transition states, and seeks the dynamical paths that connect these to each other, whereas more recent models generalize to metastable states connected by paths as a Markov state model (5). These techniques have proved to be extremely powerful, giving innumerable insights into the behavior of complex systems (6⇓⇓⇓⇓⇓⇓–13). On the other hand, a major issue has been the difficulty of finding the transition paths, connecting local minima or metastable states to each other, especially given a complex energy landscape in a high-dimensional space. A variety of creative methods have been developed in recent years for efficiently finding transition paths (14⇓⇓⇓⇓⇓⇓⇓⇓⇓–24) but for a given system, there is no guarantee that all such paths have been found.

Here we present a different point of view for understanding an energy landscape that occurs when the range over which particles interact is much smaller than their size. Such is the case in certain mesoscale systems, for example, for molecules (25, 26), or for colloids interacting via depletion (4) or coated with polymers or cDNA strands (27⇓⇓–30). We show that in this limit, the free-energy landscape is described entirely by geometry, plus a single control parameter *κ* that is a function of the temperature, depth, and curvature of the original potential. This limit is related to the sticky sphere limit of a square-well potential (31), which has been used to investigate thermodynamic properties of hard sticky spheres (32⇓–34). The landscape can be thought of as a polygon living in a high-dimensional space, whose corners (zero-dimensional manifolds) are joined to each other by lines (one-dimensional manifolds) that in turn form the boundaries of faces (two-dimensional manifolds), and so on. These manifolds are fixed functions of the particles’ geometries, independent of the details of the original interaction potential from which the limit was taken.

Once the geometrical manifolds composing the landscape are computed, nonequilibrium quantities characterizing the dynamics can be calculated by solving the Fokker–Planck equation or its adjoint on these manifolds. We show that in the short-ranged limit these equations acquire an effective boundary condition at the boundary of every *p*-dimensional manifold in the polygon. This makes the kinetics computationally tractable, because the stiff modes of a narrow potential become a set of boundary conditions.

The geometrical nature of the energy landscape does not mitigate its high dimensionality, but at low temperatures (high *κ*) both the free energies and the kinetics are dominated by the lowest-dimensional manifolds. This means that the description of the free-energy landscape and kinetics for short-ranged potentials reduces to a problem in discrete and computational geometry—to find all of the low-dimensional manifolds for a given set of interacting particles.

As an illustration of the framework, we characterize the geometrical landscape for particles with identical potentials and demonstrate how these solutions lead to a complete description of the energy landscape and the kinetics of this system. This solution describes the geometrical limit of small atomic clusters as well as a direct prediction for colloidal clusters interacting with depletion forces (4, 6, 35, 36), where the predictions could be tested experimentally. The solution also provides a framework for understanding and extending simulations on clusters with finite range potentials (1). Our calculation of the energy landscape builds on the enumeration of all finite sphere packings of *n* particles with at least contacts (37, 38). With these as the starting point, we compute every one- and two-dimensional manifold of motions that contains a zero-dimensional manifold at its boundary, from which we can extract statistical quantities such as the relative entropies of the different types of motions. Then, we solve Fokker–Planck equations on these manifolds to obtain transition rates between the lowest energy states, the zero-dimensional manifolds.

## Geometrical Landscape

We begin by showing how the geometrical free-energy landscape arises as a distinguished limit of particles interacting with arbitrary potentials of finite range. We consider a point in configuration space, with the potential energy given as a sum of potentials concentrated near different geometrical boundaries as

The functions , () represent the geometrical boundaries via their level sets , and is the potential energy near each boundary. For concreteness, let us suppose this is a model for *n* spheres with centers at , , so the configuration is , and we take to be the excess bond distance, as , where *d* is the particles’ diameter. Then is the pairwise interaction potential that we assume has minimum at and is negligible beyond some cutoff . For ease of exposition, the interaction potential is assumed to be identical between each pair of particles, but this is not a necessary restriction for the geometrical landscape to apply. For the total potential , the parameter *ε* characterizes the range of the potential, whereas is proportional to the depth. The geometrical limit requires taking as in a manner that we specify momentarily.

We consider particles that evolve according to the overdamped Langevin dynamics (39) at temperature *T*,

where *γ* is the friction coefficient, is the diffusion coefficient, , is the Boltzmann constant, and is a -dimensional white noise. The equilibrium probability for this system is the Gibbs distribution (40),

where is the normalizing constant.

The geometrical free-energy landscape occurs when the range . The relationship between the depth and the range is critical to obtaining an interesting limit. If only the range goes to zero, then particles are bonded for shorter and shorter times so that in the limit they simply behave like hard spheres. On the other hand, if the depth goes to too quickly, then the particles simply stick together and unbind only on exponentially long timescales (41). The interesting limit occurs when particles stick to each other but unbind on accessible timescales; for this we must choose so that the Boltzmann factor for two particles to be bonded approaches a finite, nonzero constant, , where we define Boltzmann factors nondimensionally by scaling by the diameter *d*. Evaluating the integral using Laplace’s method then implies

We call the constant *κ* the sticky parameter. Note that *κ* is a function of both the potential and the temperature. The constant for hard spheres.

In the geometrical limit, the probability measure becomes concentrated at the exact locations in configuration space where a bond forms, i.e., on the level sets and all possible multiway intersections. Thus, the limiting probability distribution will be a weighted sum of delta functions, each defined on a manifold corresponding to a different set of bond constraints. The weight of each delta function depends on the number of bonds and a geometrical factor associated with the entropy of the configuration. This gives a geometrical picture of the energy landscape: When *κ* is large, the occupation probabilities will be dominated by configurations where the number of bonds *m* is large. For identical particles, with , the maximum number of contacts is (37, 38): These are rigid structures that have no internal degrees of freedom, so they correspond to zero-dimensional manifolds or “points”. The next lowest configurations in potential energy are manifolds with , which are obtained from rigid structures by breaking one bond. These have one internal degree of freedom so are one-dimensional manifolds or “lines”. The lines form the boundaries of two-dimensional manifolds or “faces”, when another bond is broken, and continuing up in dimensionality we obtain the entire energy landscape as the union of manifolds of different dimensions.

Fig. 1 shows a schematic contrasting this limiting case with the traditional picture of an energy landscape. The traditional picture is of an undulating surface, with local minima connected through saddle points, whose heights provide an activation barrier that determines the transition rates between basins. In contrast, in the geometrical limit, the local minima are infinitely narrow and deep, separated by long, relatively flat spaces in between—the landscape can be thought of as a golf course punctuated with deep trenches and very deep holes. Kinetics on this landscape are determined partly by an activation barrier—the time it takes to climb out of the hole—and partly by diffusion, or the time to cross the horizontal barrier.

Fig. 1 also shows an explicit one-dimensional manifold taken from the landscape for , an example we return to throughout the text. There are two ground states, the polytetrahedron and the octahedron (38), and the manifold is the set of deformations corresponding to the transition path between these when a single bond is broken.

To explicitly calculate the equilibrium probabilities of the different states in the geometrical landscape we consider a configuration with *m* constraints or equivalently bonds broken. The constraints are written as an ordered multiindex so the manifold of configurations satisfying such constraints is

We write to mean the region where no constraints are active and let be the full space of accessible configurations. The limiting partition function associated with these constraints is

where is the neighborhood surrounding the manifold where the potential associated with the constraints is active:

This splits configuration space near each manifold into two parts—the fast variables changing rapidly along directions associated with the constraints (sometimes called the vibrational modes) and the slow variables that are the unconstrained configuration.

To compute the integral in Eq. **6** we need a parameterization of the manifolds associated with the constraints. It is convenient to parameterize the fast directions by the constraint variables themselves, . We choose the additional variables so that on ; i.e., the variables *y*, are orthogonal on the manifold. As discussed in *SI Appendix*, it is possible to find such a parameterization locally as long as the manifold associated with the constraint variables is regular—i.e., the coordinate transformation for the constraint variables must be smooth and invertible. This happens when the Hessian of the potential energy (say at )

has *m* nonzero eigenvalues.

We can now evaluate

where with , is the metric tensor associated with the transformation with components , where , and is its determinant. Let the metric tensor separate into blocks as

The first set of indexes runs from and describes the fast variables, perpendicular to the manifold, whereas the second set runs from and describes the slow variables, along the constraint manifold. Let be the determinant of a particular block. It follows from the definition of the metric that , where the s are the nonzero eigenvalues of , and the condition of orthogonality gives . Evaluating the integral in Eq. **9** over the fast variables using Laplace’s method, and letting , shows the limit is

where

is a geometrical factor (representing the “vibrational” degrees of freedom) that depends only on the level sets of the constraints, and is the metric on manifold inherited from the ambient space by restriction. The integral Eq. **11** is simply the volume integral of over .

The manifold contains 6 degrees of freedom representing translation and rotation of the cluster, and the partition function integral can be further simplified by integrating over the subspace spanned by these motions. This introduces a factor in the integral, the square root of the moment of inertia tensor (40). If we let be the quotient space formed by identifying points if *x* can be mapped to *z* by a combination of translations and rotations, i.e., , where is the special Euclidean group, then we can write

where is the metric on , and we have dropped constants (such as free volume) that are the same for all configurations. For convenience later on, we define as the part of the partition function that is independent of *κ*. *SI Appendix* has a detailed discussion for how to construct , which is critical for using the formalism developed here for practical calculations. Fig. 1 (*Bottom*) is an example of such a quotient manifold, where each point on the manifold represents the six-dimensional space of clusters in configuration space that are related by rotations or translations.

The particles in Fig. 1 are different colors to identify the different transitions that occur when moving around configuration space. However, when (as imagined here) the particles are identical, permuting the colors of any particular structure yields a geometrically isomorphic structure on a separate part of the quotient space. The free energy must account for the number of distinct manifolds that are geometrically isomorphic to . When , this is , where *σ* is the symmetry number, i.e., the number of particle permutations that are equivalent to a rotation, and if the structure is chiral and otherwise (42). For , we count the multiplicities by counting how many times a mode occurs from the perspective of each zero-dimensional “corner” of the mode and dividing by the total number of corners. (This is a combinatorial argument; it is equivalent to considering the molecular symmetry group for nonrigid molecules as in ref. 1, section 3.4.) For example, Fig. 1 has corners from two different ground states, the polytetrahedron and the octahedron, which occur with multiplicities respectively. For each polytetrahedron, there is a line coming out of it that is isomorphic to this transition, and for the octahedron there are distinct lines. (The numbers , are indicated on the arrows connecting red circles to blue circles in Fig. 2, where the transition under consideration is mode 7.) Therefore, the total multiplicity of the line is . Consider a transition connecting a polytetrahedron to a distinct copy of itself, say mode 5. Here there are such lines connected to each polytetrahedron, so the multiplicity of the line is . In general, each *p*-dimensional manifold contains a total of corners from nonisomorphic ground states, each ground state having multiplicity , and such that each single ground state is connected to distinct manifolds isomorphic to , so the multiplicity is .

Putting this all together, the total partition function of all structures isomorphic to a given constraint manifold is , and the free energy of these structures is . We can separate this free energy into , using the definition of in Eq. **13**. The first term () depends on the temperature, bond energy, and width of the potential, whereas the second term with is entirely geometrical and essentially is the entropy of structures corresponding to the constraint set . As an example, we have plotted along the polytetrahedral–octahedral transition in Fig. 1 (*Bottom*). This varies smoothly along the manifold as , vary. The endpoints of the manifold are the ground states, where the free energy changes discontinuously because the number of bonds has changed—this causes a jump in the energetic part via *m* and in the entropic part via .

With these results in hand, we can now compare the total entropies of floppy manifolds of different dimensions, to understand the temperature range in which the different manifolds occur. Let the total geometrical partition function of manifolds of dimension *p* be

Here is independent of the temperature and potential, whereas *Z* combines everything to obtain the entire landscape. Note that lower-dimensional manifolds have more bonds and thus are favored in the partition function when the temperature (or equivalently *κ*) is small. As temperature increases, *κ* shrinks and higher-dimensional manifolds become more highly populated. Eventually clusters will fall apart into single particles. The temperature dependence of how clusters fall apart is encoded in the relative sizes of the s. The temperature where the landscape transitions from having more *p*-dimensional structures than -dimensional structures is found by solving for κ, which gives roughly .

## Free-Energy Landscape for Identical Particles

To illustrate our asymptotic calculations with a concrete example we have computed the geometrical manifolds up to dimension for identical particles with diameter . To do this, we begin with the set of clusters with bonds derived in Arkus et al. (38) For every rigid cluster we break each single bond in turn and move along the internal degree of freedom until we form another bond. This is the set of one-dimensional manifolds. For the 2D manifolds, we break each pair of bonds from the rigid clusters in turn and move along the internal degrees of freedom to compute the boundaries, corners, and interior of the 2D manifolds (details in *SI Appendix*). This algorithm ensures we have every floppy manifold that can eventually access one of the rigid clusters in our list only by forming bonds, but never breaking them. Our analysis makes three assumptions that we believe to be true, but await rigorous proof: First, we are assuming that the list of clusters in Arkus et al. (38) is the complete set of rigid (zero-dimensional) clusters; this is true as long as there are no rigid clusters with bonds or less, a condition that was not checked. (The existence of such pathological examples is not ruled out purely by rank constraints on the Jacobian as there could be singular structures; see e.g., the examples in ref. 43.) Second, in the calculations of the entropy of the 2D floppy manifolds we assume that the manifolds are topologically equivalent to a disk. The fact that our parameterization algorithm works is evidence for this claim, although we have not proved this rigorously. Third, we assume that all floppy manifolds can eventually access a rigid mode and are not, for example, circles. We mention these caveats because although we are confident that they do not apply in the low *n* examples described here, it is possible that potentially significant exceptions arise at higher *n*.

The landscape for is shown in Fig. 2. (Quantitative summaries are given in the *SI Appendix*). There are two ground states, denoted by the red circles, each with contacts; the area of the red circles is proportional to the probability of each state, with the polytetrahedral ground state times more likely than the octahedral. The light blue circles denote the 5 topologically unique structures that are missing a single bond in the ground states—such structures correspond to a one-dimensional manifold, with continuous deformation along the direction of the missing bonds. The yellow circles denote the 13 unique structures that are missing two bonds from the ground state. Again, the area of the circles is proportional to the occurrence probability of these modes. These structures correspond to 2D manifolds, with continuous deformations allowed along both of the directions between the missing bonds. The connections between the different modes are denoted by arrows in Fig. 2, with structures missing (say) two bonds generally arising from breaking a single bond from several different pathways.

Each element of a mode with two bonds broken can be mapped to a polygon in , and these parameterizations are also shown in Fig. 2. We have chosen one parameterization (mode 18) to illustrate in detail in Fig. 3. The interior of the manifold represents structures with 10 bonds, with each point representing a different set of coordinates for the particles. The edges correspond to structures with 11 bonds, whereas the corners are structures with a full bonds. Of the four corners, three correspond to polytetrahedra with different permutations of the particles, and the final one is an octahedron. The one-dimensional edges connecting the corners are possible transition paths that can be followed by breaking only one bond.

In general, the number of corners varies among modes, with no a priori way to determine this without solving the full geometry problem. It ranges from three to six for and from three to seven for . Many two-dimensional modes contain several permutations of a given zero-dimensional mode as a corner.

Table 1 summarizes the partition function data. The number of different modes increases combinatorially with *n*, as do the geometrical partition functions . Strikingly, the ratios , remain virtually constant as *n* increases. This implies that the temperature dependence of the landscape is independent of *n*. Indeed, Fig. 4 shows the relative probabilities of zero-,one-, and two-dimensional modes as the temperature varies, for , with the yield of a *p*-dimensional mode given by . As an illustration, we have chosen , , so that . Because the ratios are nearly the same, these graphs are essentially indistinguishable for different numbers of particles. Moreover, the critical temperature for transitioning from mostly zero-dimensional structures to one-dimensional structures () is quite close to that for transitioning from one-dimensional structures to two-dimensional structures (). If this remains true as *p* increases, it would imply that clusters melt explosively at some critical temperature, rather than incrementally: Clusters would occupy either mostly the zero-dimensional modes or a gaseous, no- or few-bond-state, but not the chain-like floppy configurations in between.

## Kinetics on the Geometrical Landscape

We now consider kinetics on the geometrical landscape. The concentration of equilibrium probabilities on manifolds with varying dimensions also applies to nonequilibrium quantities, such as transition rates, first-passage times, the evolution of probability density, etc. Such quantities can be computed from the time-dependent transition probabilities, which for dynamics given by Eq. **2** are obtained from the corresponding forward or backward Fokker–Planck (FP) equations (44). We show that in the geometrical limit, the FP equation asymptotically approaches a hierarchy of FP equations, one on each manifold of each dimension. These equations are coupled to each other, with equations on manifolds with dimension *p* serving as boundary conditions for the equations on manifolds with dimension .

The idea behind the derivation is quite natural: If a potential is deep and narrow, then it equilibrates much more rapidly in the directions along the bonds than along a cluster’s internal degrees of freedom. Therefore, the probability density near a cluster with *p* bonds broken is approximately , where *y* parameterizes the internal degrees of freedom of the cluster. The “constant” evolves slowly in the transverse directions according to the Fokker–Planck dynamics (which include an “effective” potential that arises from the curvature of the manifold), and it also changes due to the flux of probability out of the -dimensional manifolds for which it forms part of their boundary. This gives a hierarchy of coupled FP equations.

To see how this comes about in detail, we examine solutions to the Fokker–Planck equation for the evolution of the probability density . Given a parameterization of configuration space with metric tensor *g*, the nondimensional Fokker–Planck equation corresponding to Eq. **2** is

with boundary conditions at each level set

where is the outward normal to the boundary. We have nondimensionalized lengths by *d*, times by , and energy by . Away from all boundaries there is no force, so in the limiting probability *p* evolves only by diffusion as

Now consider the evolution near a manifold , a “wall”. The dynamics in the directions orthogonal to the wall, where bond distances are changing, are much faster than those along it, so near the wall the probability density will rapidly approach a multiple of the equilibrium probability. Parameterizing the region near the wall as as in the previous section, we obtain

where is the correction to the leading-order formula. This satisfies the matching condition that be asymptotically continuous. This ansatz can also be derived from a consistent asymptotic expansion of Eq. **15** after the change of variables .

Substituting Eq. **17** into the Fokker–Planck Eq. **15** gives

We want to integrate out the fast variables so we need to separate these from the slow ones. This is most conveniently done using the metric *g* with block decomposition (Eq. **10**). Making the change of variables shows that Eq. **18** can be written as

where we abbreviate .

We now integrate Eq. **19** over the fast variables and keep the leading-order parts. The terms vanish in the limit because . The terms require evaluating an integral similar to Eq. **9**, which converts to the factor at each point on the manifold. The first term is the most interesting. This is the divergence in the fast variables, and although *p* does not depend on these, the unknown might and contributes at . However, we can use the divergence theorem to replace this term with an integral of the flux through the *p*-dimensional fast-variable boundary, at . We then introduce a second matching condition that requires the flux to be continuous, so we can replace it with the sum of fluxes from each -dimensional manifold that has as part of its boundary and evaluate the limit of the boundary integral for these matched fluxes. [Specifically, we require that , where is the normal to level set and is the normal to level set , and with a similar expression for .]

Finally, we integrate over the space of rotations and translations of each point, assuming is constant on orbits. The divergence in these directions will disappear by Stokes’ theorem, and the remaining directions provide dynamics on the quotient space. Combining with the previous calculations yields

where the fluxes have leading-order part

Here combines the sticky parameter and the geometric factor at the wall, and the metric tensor is the quotient metric obtained from the metric on , which in turn is inherited from the original metric *g* in the ambient space by restriction: . The sum is over *β* such that is part of the boundary of the -dimensional manifold and where is an outward normal vector to at . We define to be the differential operators on the quotient manifold .

Eq. **20** has a more intuitive interpretation as the evolution of the probability along the wall. It is clear from the derivation that the total probability density (with respect to the wall coordinates) of being on the wall is . This satisfies

The dynamics along the wall are therefore a combination of diffusion, plus drift due to an effective potential , plus a flux in and out of the wall.

The effective potential is entropic in nature and comes from the changing wall curvature, which makes the potential look wider in some places than in others. A particle will spend more time in the wide places than in the narrow ones, and because it reaches equilibrium much more quickly in the transverse directions than in the along-wall directions, it looks like there is an effective force pushing it to the wider areas. This is the same equation one obtains by letting the depth of the potential become infinite without changing the width (41), but with the addition of the flux in and out of the wall.

This flux term is illustrated schematically in Fig. 5 for a simple case where the configuration space has a single constraint, —a wall at the horizontal axis. The probability integrated over a small box (red) with length near the wall changes in a time increment due to two processes: probability fluxing along the wall in the *x*-coordinate, which contributes a change of , and the flux of probability from the wall to the interior, which contributes a change of . Equating with gives the effective boundary condition.

To summarize, the limiting FP equation and boundary conditions are Eq. **16**, plus Eq. **20** on every . Substituting for the time derivatives shows that the boundary conditions are second order, and this is why conditions are needed on every boundary and not only on those with codimension 1. We call this set of equations the “sticky” equations because they are the Fokker–Planck equations for sticky Brownian motion (45), a stochastic process that has a probability atom on the boundary of its domain.

## Transition Rates for Sticky Brownian Motion

To transition between the rigid configurations in the geometrical limit, a cluster must break one (or more) bonds and then diffuse across the line segment (or face, etc.) until it hits the other endpoint. The time it takes to do this depends on the length of the line and on how the entropy of the configuration varies along the line, and we can find this by solving an equation on the full line segment (face, etc.). In our asymptotic limit there are no meaningful transition “states”—rather, the entire line segment can be thought of as a transition state. Once the energetic barrier of breaking a bond has been overcome, transitions are dominated by diffusion.

We now consider how to compute transition rates using the sticky equations (Eq. **20**), supposing we have a stochastic process whose probability evolution is well approximated by these. (Note that we have not actually constructed a process that satisfies this in the strong sense. However, we began with a process satisfying Eq. **2** and showed the sticky equations describe its probability evolution asymptotically so we use these to compute rates.) Consider the transition rate between sets , , and for simplicity let us focus only on the case where these are both (disjoint) subsets of the zero-dimensional manifolds. For example, we may be interested in the transition rate between the octahedron and the polytetrahedron (introduced in Fig. 1), in which case *A* would contain all points in the quotient space representing the octahedron and *B* all those representing the polytetrahedron. These rates can be computed using transition path theory, which provides a mathematical framework for computing transition rates directly from the Fokker–Planck equations. We simply state the facts that are relevant to our example and refer the reader to other resources for more details (46⇓–48).

The committor function is the probability, starting from point *x*, of reaching set *B* before set *A*. This solves the stationary backward Fokker–Planck equations with boundary conditions , , plus any other boundary conditions remaining from the equations. As shown in *SI Appendix*, the forward sticky equations (Eq. **20**) are self-adjoint (with respect to the invariant measure) so these are also the backward sticky equations.

A reactive trajectory is a segment of the path that hits *B* before *A* going forward in time and *A* before *B* going backward in time. The probability current of reactive trajectories is a vector field that, when integrated over a surface element, gives the net flux of reactive trajectories through it. Because our process is time reversible, this current is (47)

where is the equilibrium probability measure. From Eq. **13** we find this is

where is the singular measure on ; i.e., it satisfies . Finally, the transition rate is calculated by integrating this flux over a surface *S* dividing the two states, giving

where is the normal pointing from the side containing *A* into the side containing *B*.

Computing this exactly requires solving the backward equations on a high-dimensional space and integrating over a high-dimensional surface—a computationally infeasible proposition. However, when the sticky parameter *κ* is large, most of the probability is concentrated on the lowest-dimensional manifolds so we expect these to contribute the most to the rates. Therefore, let us expand the equations in powers of . We suppose all variables have an asymptotic expansion as

and similarly for *q*, *μ*, *ρ*, *J*, etc. Expanding *ρ* shows that to first order it is a measure on points, to second order it is a measure on points and lines, etc. Measures on points will not contribute to the rate because the dividing surface *S* can be chosen to avoid them, so and the leading-order part of the rate is , computed from .

Expanding the backward sticky equations in powers of gives a set of equations for ,

with boundary conditions , , and on all other zero-dimensional manifolds.

To solve these equations we can first find the solution on the zero- and one-dimensional manifolds and then use this as a boundary condition for the solution on the two-dimensional manifolds, which becomes in turn a boundary condition for the solution on the three-dimensional manifolds, etc. The leading-order rate requires only the solution for . If we enumerate the lines connecting a point to a point and use an arc-length parameterization for the *k*th line whose total length is , then this is given analytically as , where

is the normalization factor. On all other lines .

Any dividing surface *S* hits each line at a single point, so the leading-order rate (in dimensional units) is

where the sum is over all connecting lines. This is asymptotically equivalent to the rate one would obtain simply by restricting the full committor function and invariant measure to the set of zero- and one-dimensional manifolds.

### Transition Rates for Hard Spheres.

We used the formalism to compute the leading-order rates for hard spheres with diameter . To do this we took the set of one-dimensional solutions computed as part of the free-energy landscape and computed the factor from Eq. **29** on each manifold. Summing over all of the one-dimensional manifolds that connect a zero-dimensional manifold numbered *a* to a zero-dimensional mode numbered *b* gives the transition rate between *a* and *b.* We also include transitions between different ground states belonging to the same mode (e.g., ), by multiplying the previous calculation by 2 because transitions can go in either direction along the line.

Fig. 6 shows the network of zero-dimensional states and the reaction rates between each state. The numbers reported are the dimensionless, purely geometrical parts of the rates and should be multiplied by to give the dimensional rate. These rates, when multiplied by the total time of a simulation or experiment, give the average number of transitions one would expect to observe, so they are equal for both directions and because our system is time reversible. To obtain the rate relative to a particular state, i.e., the rate at which one leaves state *a* to visit state *b* next, given that the last state visited was state *a*, one should divide by the so-called reactive probability of *a* (49). To leading order, this is equivalent to dividing by the equilibrium probability of *a*.

### Simulations.

We have verified our results by performing Brownian dynamics simulations of Eq. **2** with a short-range Morse potential. The results agree very well with our calculations of both free-energy and transition rates. Fig. 7 shows a comparison of the simulated probabilities vs. theoretical probabilities of each mode for (see *SI Appendix* for ), for particles interacting with a Morse potential with range parameter , a range of ∼5% of the particle diameter. The agreement is nearly perfect. Fig. 7 also compares the number of each type of transition we saw in the simulations to that predicted from theory. The theory slightly overpredicts the total number of transitions—this is what one should expect from the geometrical picture, as these leading-order rates neglect the time the system spends in the floppy manifolds, which would tend to slow it down.

These results are encouraging partly because they are evidence that we have executed these calculations correctly, but also because they suggest the asymptotic limit may apply for experimental systems, such as Meng et al. (4) where the potential reportedly had roughly this width.

### Comparison with Other Numerical Approaches to the Free-Energy Landscape.

We have compared our results with those obtained by a numerical study that directly searched an energy landscape of a short-ranged potential (a Morse potential, with range roughly 0.05 particle diameter) for local minima and transition states (42). The numerical method found fewer local minima than there are zero-dimensional modes and fewer transition states than one-dimensional modes, suggesting that the asymptotic theory has the roughest landscape. By computing adjacency matrices for the numerical states with a cutoff bond distance of 1 + 1 , we verify that for each local minimum corresponds to a unique zero-dimensional mode, and each transition state lies on a unique one-dimensional manifold.

We identify the point on the one-dimensional manifold that is closest to each transition state. The transition states are very close to the local maxima of the vibrational factor (see *SI Appendix* for plot), consistent with it being a saddle point of the potential energy. We believe the small discrepancy in location can be attributed to the finite width of the Morse potential used in the numerical procedure.

The missing zero-dimensional modes occur for manifolds that are very close to each other in the quotient space metric (separation ), which is the case for modes () and modes (). For these only one local minimum was found for the entire group. We hypothesize that because the separation is within the range of the potential, the zero-dimensional modes merge to form a single local minimum.

The missing one-dimensional modes often, but not always, correspond to self–self transitions—these transitions do not matter when particles are identical; however, they will account for transitions between different states when the particles are not all the same (e.g., ref. 50). For example, for the numerical procedure identifies 45 transition states, compared with our 75 one-dimensional manifolds. Of the missing manifolds, 16 are self–self transitions, 9 are nonself–nonself transitions within the group , and 5 are nonself–nonself transitions for endpoints not both in the group.

## Discussion/Conclusions

We have developed a unique framework for understanding energy landscapes when particles interact with a short-ranged potential. We show that in the limit as the range goes to zero and the depth goes to , the energy landscape becomes entirely governed by geometry, with a single parameter *κ* encapsulating details about the potential and temperature. When *κ* is large, only the lowest-dimensional geometrical manifolds contribute significantly to the landscape and this makes a computational approach tractable. To illustrate the limit, we have computed the set of low-dimensional manifolds for hard, spherical particles. This solution to a nontrivial problem in statistical mechanics can be used to compute equilibrium or nonequilibrium quantities for any potential whose range is short enough.

We were able to calculate this set of low-dimensional manifolds because we began with the set of rigid clusters and made the conjecture that all floppy modes can be accessed from these by breaking bond constraints. Solving for the complete set of rigid clusters is a difficult problem in discrete geometry that has been done only for (38, 51), but with current computational power and novel approaches (51, 52) one can anticipate reaching larger *n*. Very large *n* will eventually require making approximations to the geometry problem. We speculate that as *n* increases, structures with extra bonds, as well as “singular” structures whose Jacobians have extra zero eigenvalues, will come to dominate the landscape—these have not yet been considered in our asymptotic framework but they are observed with high probability in experiments (4).

We have compared our results to those obtained by numerically searching the free-energy landscape of a short-ranged Morse potential for local minima and transition states. Our method finds more minima and transition regions of the potential energy than the numerical search procedure, and this points to a potentially useful extension of our theory—one can imagine starting with the limiting geometrical manifolds and following these in some way as the range of the potential is increased, to obtain a low-dimensional approximation to the free-energy landscape for finite-width potentials, such as Lennard–Jones or Van der Waals clusters. This method would overcome a major issue with numerical searches, which is that there is no way to ensure that all important parts of the landscape have been found—we claim that our manifolds are the complete set of low-energy states so they will remain so under small enough perturbations. In addition, this would provide a way to deal with the increasing ruggedness of energy landscapes with short-ranged potentials, which are a challenge for numerical methods—we start with the most rugged landscape and would need only to smooth it.

## Acknowledgments

We thank David Wales for generously providing data and Vinothan Manoharan and Eric Vanden-Eijnden for helpful discussions. This research was funded by the National Science Foundation through the Harvard Materials Research Science and Engineering Center (DMR-0820484), the Division of Mathematical Sciences (DMS-0907985), and the Kavli Institute for Bionano Science and Techology at Harvard University. M.P.B. is an investigator of the Simons Foundation.

## Footnotes

- ↵
^{1}To whom correspondence should be addressed. E-mail: holmes{at}cims.nyu.edu.

Author contributions: M.H.-C. and M.P.B. designed research; M.H.-C. and M.P.B. performed research; S.J.G. contributed new reagents/analytic tools; and M.H.-C. and M.P.B. wrote the paper.

The authors declare no conflict of interest.

This article is a PNAS Direct Submission.

See Author Summary on page 12 (volume 110, number 1).

This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1211720110/-/DCSupplemental.

## References

- ↵
- Wales DJ

- ↵
- Stillinger F,
- Weber T

- ↵
- Beberg AL,
- Ensign DL,
- Jayachandran G,
- Khaliq S,
- Pande VS

*IEEE International Symposium on Parallel & Distributed Processing*, May 23–29, 2009, IEEE Computer Society. - ↵
- Meng G,
- Arkus N,
- Brenner MP,
- Manoharan VN

- ↵
- Bowman GR,
- Pande VS

- ↵
- ↵
- Stillinger FH

- ↵
- ↵
- Liwo A,
- Lee J,
- Ripoll DR,
- Pillardy J,
- Scheraga HA

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- Weinan E,
- Ren W,
- Vanden-Eijnden E

- ↵
- ↵
- ↵
- Abrams CF,
- Vanden-Eijnden E

- ↵
- Wales DJ,
- Scheraga HA

- ↵
- ↵
- ↵
- ↵
- Noé F,
- Schütte C,
- Vanden-Eijnden E,
- Reich L,
- Weikl TR

- ↵
- ↵
- ↵
- ↵
- ↵
- Macfarlane R,
- et al.

- ↵
- Rogers WB,
- Crocker JC

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- Allen MP,
- Tildesley DJ

- ↵
- Landau LD,
- Lifshitz EM

*Statistical Physics*, Course of Theoretical Physics (Pergamon, Oxford), Vol 5. - ↵
- ↵
- ↵
- ↵
- Gardiner C

*Stochastic Methods: A Handbook for the Natural and Social Sciences*(Springer, Berlin), 4th Ed. - ↵
- Ikeda N,
- Watanabe S

- ↵
- ↵
- ↵
- ↵
- ↵
- Hormoz S,
- Brenner M

- ↵
- ↵
- Sommese A,
- Wampler C

## Citation Manager Formats

## Sign up for Article Alerts

## Jump to section

## You May Also be Interested in

### More Articles of This Classification

### Related Content

- No related articles found.

### Cited by...

- Colloidal matter: Packing, geometry, and entropy
- Size limits of self-assembled colloidal structures made using specific interactions
- Mesoscale molecular network formation in amorphous organic materials
- Hydrodynamics selects the pathway for displacive transformations in DNA-linked colloidal crystallites