# Exploring canyons in glassy energy landscapes using metadynamics

Amruthesh Thirumalaiswamy  Robert A. Riggleman<EMAIL_ADDRESS>John C.
Crocker<EMAIL_ADDRESS>Department of Chemical and Biomolecular
Engineering, University of Pennsylvania, Philadelphia, Pennsylvania

###### Abstract

The complex physics of glass forming systems is controlled by the structure of
the low energy portions of their potential energy landscapes. Here, we report
that a modified metadynamics algorithm efficiently explores and samples low
energy regions of such high-dimensional landscapes. In the energy landscape
for a model foam, metadynamics finds and descends meandering ‘canyons’ in the
landscape, which contain dense clusters of energy minima along their floors.
Similar canyon structures in the energy landscapes of two model glass
formers—hard sphere fluids and the Kob-Andersen glass—allow metadynamics to
reach low energies. In the hard sphere system, fluid configurations are found
to form continuous regions that cover the canyon floors, but only up to a
volume fraction close to that predicted for kinetic arrest. For the Kob-
Andersen glass former, metadynamics reaches the canyons’ ends with modest
computational effort; with the lowest energies found approaching the predicted
Kauzmann limit.

††preprint: APS/123-QED

##  I Introduction

Many outstanding problems in understanding glasses have been related to their
potential energy landscapes [1], hypersurfaces describing a system’s total
potential energy spanning the high-dimensional configuration space formed by
all the particles’ spatial coordinates [2]. Navigating such high-dimensional
spaces and mapping these vast landscapes has been a major challenge. The
heterogeneous nature of glassy dynamics [3] makes the use of collective
descriptors of system dynamics ineffective. Previous work has revealed a broad
distribution of energy basin hypervolumes [4, 5], highly tortuous steepest
descent paths [6, 7], and fractal clustering of local minima [7]. Rapid
quenching from a very high temperature (equivalent to following a steepest
descent path from a random configuration) yields local minima, or inherent
structures (ISes) [2], weighted by their associated basin hypervolume [4, 5].
The ensemble of such quenched configurations [8] having the largest basin
hypervolumes are distinctly different from the lowest energy inherent
structures [9, 10] that control glass transitions [11, 12, 13]. Conventional
methods for exploring low energy states, such as molecular dynamics and Monte
Carlo, naturally become inefficient at the low temperatures where the
corresponding physical system undergoes kinetic arrest [11, 9]. Moreover,
methods like eigenvector-following [14], and techniques for exhaustively
enumerating inherent structures [15] tend to be costly [16, 17] in large
systems.

Here, we modify a high-dimensional metadynamics [18] algorithm by Yip [19],
pairing it with with a fast energy minimizer (FIRE) [20]. The result uses
localized bias potentials to cross energy barriers enabling the efficient
exploration of glassy potential energy landscapes. We term our approach the
Metadynamics-Inspired Multi-fractal Sampling Explorer (MIMSE). Applying this
algorithm to a model of foams or soft glassy materials (SGMs) [7], we discover
that it finds and rapidly descends canyon-like structures in the landscape
that contain many inherent structures along the canyons’ floors, but very few
along the canyons’ smooth walls. Analyzing small ensembles of biased
trajectories reveals that these canyons have well-defined widths, which become
narrower at lower energies. To apply this approach to a popular model of dense
hard-sphere (HS) fluids [9], we examine the corresponding extended energy
landscape for compressible or soft spheres. The algorithm again finds and
descends canyons, finding HS configurations (where the soft sphere potential
energy is zero) that cover the canyon floor. These configurations form a
continuous, connected domain up to but not above a maximum volume fraction,
$\phi_{max}\approx 0.669$, near the value predicted for the HS fluid to
undergo kinetic arrest [9]. Last, when applied to a simple model of an atomic
glass former, the Kob-Andersen (KA) model [21], our algorithm again finds and
rapidly descends canyons, reaching their ends with modest computational
effort. The canyons’ ends contain low-dimensional and fractal clusters of
inherent structures. Exploring many canyons yields a distribution of canyon
bottom energies that approaches the predicted Kauzmann energy limit [22].

##  II The MIMSE Algorithm

MIMSE is an athermal, metadynamics based approach [18] directly applied to the
$3N$-dimensional potential energy landscape of $N$ particles, and uses
sequentially added bias potentials to overcome energy barriers [19]. We use
the Fast Inertial Relaxation Engine (FIRE) [20] to relax the system
configuration on the biased energy landscape; tests using steepest descent
gave similar results but at higher computational cost. The bias potential
employed is a $3N$-dimensional, bell-shaped, smooth quartic function with the
form

$U(\mathbf{r})=\begin{cases}\mathcal{U}_{0}{\left(1-{\left(\frac{\lVert{\underaccent{\bar}{\mathbf{r}}-\underaccent{\bar}{\mathbf{r}}_{\text{m}}}\rVert}{\mathcal{U}_{\sigma}}\right)}^{2}\right)}^{2},&\text{if
$\lVert{\underaccent{\bar}{\mathbf{r}}-\underaccent{\bar}{\mathbf{r}}_{\text{m}}}\rVert<\mathcal{U}_{\sigma}$}\\\
0,&\text{otherwise}.\end{cases}$ (1)

Here, $\underaccent{\bar}{\mathbf{r}}$ represents the $3N$ dimensional system
configuration, and $\underaccent{\bar}{\mathbf{r}}_{\text{m}}$ the IS minimum
location around which the bias is centered. The bias parameters
$\mathcal{U}_{\sigma}$ and $\mathcal{U}_{0}$, represent the $3N$-dimensional
radial extent and energetic height of the bias respectively. To ensure that
the center of mass remains at rest, we modify the corresponding bias forces to
ensure the total force on the system is zero. To ensure stable configuration
dynamics, we use models with energy landscapes that are continuous and
differentiable.

Algorithmic flow starts at a local minimum of the landscape, where a bias is
added centered on that location, forming a local maximum on the biased
landscape. To move away from this unstable maximum, the system is given a
small, random $3N$-dimensional displacement (corrected so as to preserve the
position of the system’s center of mass). This configuration is then relaxed
to the nearest minimum on the biased landscape using FIRE, and the process
repeated by adding a new bias at each subsequent minimum. We distinguish among
the minima encountered and primarily store and analyze the subset which lie
outside the bias potentials, which correspond to physical inherent structures
(ISes) of the unbiased energy landscape.

In conventional metadynamics, the bias radius is smaller than the separation
between neighboring ISes, leading to basin filling and exhaustive enumeration
of minima [18, 19]. Here, we use larger bias radii, such that when a bias is
added, many ISes within $\mathcal{U}_{\sigma}$ of
$\underaccent{\bar}{\mathbf{r}}_{\text{m}}$ are also covered by that bias,
effectively forcing the system over their associated nearby energy barriers.
Values of $\mathcal{U}_{\sigma}$ and $\mathcal{U}_{0}$ that lead to efficient
landscape descent in energy are found by manual ‘tuning’ using a divide and
conquer-like approach, see Appendix A. For computational efficiency, a
$3N$-dimensional neighbor list is maintained to track the different biases
affecting the system at its current $3N$-dimensional position. In practice, we
find that reasonably ‘old’ biases can often be retired without difficulty,
keeping computational and memory costs manageable. Additional details of the
algorithm are provided in Appendix A.

##  III SGM Energy Landscape

(a)(a)

(a)

(a)

(a)

(b)

Figure 1: Exploring a soft-glassy material landscape with MIMSE finds low-
energy ISes efficiently. (a) The left panel shows the probability distribution
of energies for an ensemble of $1000$ different ‘quenched’ initial
configurations, formed by FIRE relaxation of random configurations. Applying
MIMSE, the energy of the configurations decreases rapidly, depending on bias
radius, $\mathcal{U}_{\sigma}$. The right axis rescales the energy relative to
the mean and standard deviation of the quenched energies $U_{q}$, also known
as the z-score, $\mathcal{Z}_{U}$. (b) The descending energy trajectory
appears roughly logarithmic in time. Symbols indicate three inherent
structures that are analysed in detail in Fig. 2. (Inset) The number of biases
that are required to descend a given energy range (dashed lines in Fig. 1b)
shows a strong dependence on $\mathcal{U}_{\sigma}$, which can be approximated
by a power-law over a narrow range.

We first demonstrate our approach on a model of foam or soft glassy matter
(SGM) [23, 24, 25], whose energy landscape was characterized in an earlier
study [7]. The system consists of a highly polydisperse collection of $N\sim
350$ compressible ‘soft’ spheres interacting via a purely repulsive harmonic
interactions when they overlap, at a volume fraction, $\phi=0.75$, above the
threshold for solidity [7]. The distribution of radii ($R$) resembles that of
bubbles in a ripening foam, described by a Weibull distribution with
polydispersity $\Delta=(\langle R^{2}\rangle-{\langle
R\rangle}^{2})^{1/2}/\langle R\rangle=0.64$, see Appendix B for model details.
To create initial configurations of this model, we select random points in
configuration space and relax them to their nearest energy minimum using FIRE.
This ensemble of ‘quenched’ inherent structures has a roughly Gaussian
distribution of energies, see Fig. 1a (left panel).

Our key finding is that for certain bias parameters, MIMSE rapidly descends
down the potential energy landscapes of glassy systems. Typical results for
the SGM model landscape are shown in Fig. 1, reaching energies that are much
lower than the energies of quenched configurations after only a few hundred
biases. These proceed to lower energy lower logarithmically with additional
effort, reaching energies which, compared to the quenched energy distribution,
are $\mathcal{Z}_{U}\approx-5$ standard deviations below the mean of the
initial quenched inherent structures, Fig. 1(right scale). Notably, such
significant descent in energy is only found to occur in a window of bias
radii, $2\lesssim\mathcal{U}_{\sigma}\lesssim 5$, and the number of biases
required to descend by a given energy shows a strong dependence on the bias
radius, Fig. 1b(inset). Larger biases descend the landscape far faster than
smaller ones, up to a maximum radius. The energetic height of the bias
$\mathcal{U}_{0}$ appears less critical, descent down the energy landscape is
observed in a wide window,
$10^{-3}\lesssim\mathcal{U}_{0}/\mathcal{U}_{\sigma}\lesssim 10^{3}$.

Our first task is to understand what features of the SGM landscape enable our
algorithm to descend to such low energies, and to explain the narrow window of
bias radii where it occurs. To this end, we formed ensembles of $1100$ FIRE
trajectories, each descending the same bias around the same starting IS
minimum, but initialized with different random $3N$-dimensional displacements,
and examined the results as a function of bias radius. We find that, after
initially moving ‘radially’ away from the center in random directions due to
the strong bias force, these paths then drift in the ‘angular’ $(N-1)$
dimensions according to the gradient of the underlying landscape. Finally,
these paths terminate in $1100$ different energy minima near the edge of the
bias.

To descend the energy landscape, the FIRE minimizer must be following
gradients in the underlying landscape; perhaps the window of bias radii is due
to suitable gradients only being present on corresponding lengthscales? To
test this idea, we first compute the average ‘drift’ of the ensemble of
trajectories, such as would be caused by an underlying gradient. Specifically
we first compute a 3N-dimensional unit vector $\hat{u}$ pointing from the
initial IS to the new energy minimum. Next, we compute the Euclidean length of
the average of these unit vectors, $\|\langle\hat{u}\rangle\|$, akin to the
center of mass of the new minima, as a function of bias radius,
$\mathcal{U}_{\sigma}$. The results for three starting ISes at different
depths in the landscape are shown in Fig. 2a. While the average drift shows a
broad peak for the $\mathcal{U}_{\sigma}$ where energy descent is observed, a
similar amount of drift is seen for much smaller and larger radii, failing to
explain the lack of energy descent in those cases.

(a)(a)

(a)

(a)

(a)

(b)

Figure 2: Trajectory ensembles reveal landscape gradients and IS clustering.
Ensembles of $1100$ FIRE trajectories descending the same hyperspherical bias
potential after different random initial displacements, move radially before
being affected by the underlying landscape and terminating near the bias edge
(inset schematics). The final relative displacements are converted to unit
vectors, $\hat{u}$. (a) The norm of the average of the $1100$ unit vectors
$\|\langle\hat{u}\rangle\|$, reveals the average drift due to an energy
gradient, for three different depths in the landscape:
$\mathcal{Z}_{U}\approx-2.8$ (squares), $-3.8$ (circles), $-4.7$ (triangles),
showing a broad peak around $\mathcal{U}_{\sigma}\approx 3$. (b) Performing
single-linkage hierarchical cluster analysis [26] on the same $\hat{u}$
ensembles yields a single cluster only for $\mathcal{U}_{\sigma}\gtrsim 2$,
with an upper limit that depends on $\mathcal{Z}_{U}$. Landscape energy
descent only occurs for bias radii and energies where a single cluster is
found, as explained in the text.

To gather clues about barrier crossing in the underlying landscape, we applied
a hierarchical single-linkage clustering algorithm [26, 27] to the unit vector
ensembles analyzed above. The idea is that if a subset of the trajectories are
crossing a large barrier, they will form a separate cluster from the others.
Remarkably, we now find a one-to-one correspondence between the bias radii
that give efficient energy descent and the $\hat{u}$ vectors that form a
single cluster, for ISes at the three different energies considered, Fig. 2b.
The maximum $\mathcal{U}_{\sigma}$ that forms a single cluster gets
progressively smaller at lower energy domains in the landscape, decreasing
from $7$ to $6$ to $5$ units. The correspondence between bias radius and
energy descent is confirmed by another observation: biases radii in the range
$5-7$ rapidly descend the landscape at first, but then get ‘stuck’, plateauing
at intermediate energies, see Appendix B, Fig. 6. Analysis of the same vectors
using complete linkage clustering fails to resolve these clusters well,
suggesting that they are extended and intertwined, rather than compact and
well-separated.

These findings indicates how the algorithm works: it is simply following
clusters of ISes like a trail of ‘bread-crumbs’ to low energy portions of the
landscape. Finding multiple clusters with too large biases indicates the
crossing of barriers into other domains, losing the trail. Finding many
clusters with too small biases indicates the failure to cross many small-scale
barriers, and the algorithm remains trapped within IS clusters at nearly the
same location. Lastly, a fractal scaling analysis, Appendix C, of the
configuration space path (or trail) followed with properly sized biases
reveals the path is fractal on the longest lengthscales, with a fractal
dimension close to that of a random walk.

(a)(a)

(a)

(a)

(a)

(b)

Figure 3:  Schematic ‘aerial view’ of the energy landscape showing ‘canyons’
(white region) with tortuous bottoms (dashed line) descending to lower energy
(left to right). (a) For the SGM system, inherent structures (ISes) (black
dots) are clustered along the canyon bottom and MIMSE skips from IS to IS
(solid red path), when biases (large red circles) are sequentially applied
(left to right). Red open circles represent biased minima; Unbiased minima
(black arrows) are stored as output. (b) For the HS glass, the bottom of the
canyon floor is covered by a slender connected domain having zero energy
(blue), as if flooded by a ‘stream’ of water, and MIMSE skips between
configurations where the bias and the zero energy domain intersect (small red
circles). Some ISes in (a) may have zero energy (blue) but may be either
blocked by other ISes or too rare to sample.

The above analyses can be summarized schematically in Fig. 3a. The landscape
consists of meandering canyons or tubes (as suggested by earlier studies [6,
7]), with the ISes forming extended dense clusters along the canyon floors.
These dense clusters act as traps for MIMSE when the bias radius is too small.
For optimal bias radii, as shown in Fig. 3a, the biased FIRE path bounces off
the canyon walls (finding no ISes there), follows the landscape gradient and
ends at a new inherent structure further down the canyon. Sometimes this path
extends out of the biased region completely, ending in a physical IS. Too
large biases push the configuration into adjacent canyons apparently
preventing energy descent. Since the largest $\mathcal{U}_{\sigma}$ that
descends the landscape decreases at lower energies, it appears that the canyon
tapers as it descends to lower energy. Given that the canyon is high-
dimensional, the reduction in its hypervolume due to even a slight radial
decrease is enormous.

##  IV Hard Sphere fluid

To determine if the energy landscapes of glass formers also contain canyon
structures that resemble those in the SGM system, we apply MIMSE to two well-
studied models. The first is a model of hard-sphere (HS) fluids [9, 10, 28],
having a size distribution given by $\textbf{P}(R)\sim
R^{-3},R\in[R_{\text{min}},R_{\text{max}}],R_{\text{min}}/R_{\text{max}}=0.4492$,
and $N=1000$ particles. Such a polydisperse system has been shown to
consistently avoid crystallization, instead undergoing kinetic arrest into a
‘hard sphere glass’ as the volume fraction is increased [9, 10].

Importantly, the hard sphere potential segments configuration space into two
domains, a physical one with zero energy and a forbidden domain with one or
more particle overlaps. To form a continuous energy landscape, we consider a
‘soft-sphere’ extension to the HS model, with harmonic repulsion. See Appendix
B for details. This model thus resembles the SGM model, but has a different
particle radius distribution and volume fraction. We then consider any
configurations found with energy per particle less than a small tolerance
($U/N<U_{tol}=10^{-16}$) to be physical HS configurations [9]. For such states
with zero energy, we use a different order parameter to follow the progress of
our algorithm across the landscape, corresponding to the mean coordination or
contact number $\langle z\rangle$. Since pair contacts containing
unconstrained or ‘rattler’ particles can be trivially relaxed, $\langle
z\rangle$ is a measure of the force bearing contacts. Its value is confirmed
to be insensitive to the choice of $U_{tol}$, see Appendix B, Fig. 7.

(a)(a)

(a)

(a)

(a)

(b)

Figure 4: MIMSE descends both the soft-sphere and HS portions of the HS glass
landscape. (a) MIMSE efficiently descends down canyons, reducing the energy
(left panel) of the soft-sphere system, ending in a nearly isostatic, ‘jammed’
HS state. Further application of biases (right panel), shows that the average
coordination number of these zero-energy configurations, $\langle z\rangle$
drops significantly, indicating structural relaxation, before reaching a
steady state value. (b) Above a volume fraction, $\phi_{max}$, MIMSE returns
finite energies, when initialized with a quenched configuration (blue circles)
or with a zero energy state (red triangles). The dashed lines (blue and red)
indicate the transitions observed using the the two different initialization
schemes.

We perform MIMSE relaxation experiments varying the volume fraction over a
small range, $0.66<\phi<0.67$. To provide initial ISes, we FIRE relax random
configurations. Since this range is above the ‘jamming’ or random close packed
volume fraction for this model, $\phi_{J}\simeq 0.65$, this initial quench
step consistently yields ISes having finite energy, residing on the soft
sphere portion of the landscape. Applied to these quenched states, we find
MIMSE again rapidly descends the energy landscape, Fig. 4a, when the bias
radius is in a small range ($3.5\lesssim\mathcal{U}_{\sigma}\lesssim 4.5$),
indicating that there are canyons similar to those in the SGM landscape.

For the lower range of volume fractions we study, $\phi\lesssim 0.667$ MIMSE
is able to consistently reach states with zero energy, corresponding to hard
sphere configurations, Fig. 4a. The first HS configurations found by MIMSE are
isostatic (or nearly so), with $\langle z\rangle\approx 6$. Further
application of biases produces a continuous series of HS states, with $\langle
z\rangle$ dropping by roughly three orders of magnitude after $\sim 10^{3}$
biases, Fig. 4a. Notably, once it finds its first zero-energy state, MIMSE
finds a new, unbiased zero-energy state after every bias addition, and no
finite energy ISes. The $g(r)$ of these configurations on the $\langle
z\rangle$ plateau resembles the jammed case for large $r$ suggesting that the
algorithm is sampling the landscape of the marginal glass adjacent to the
jamming line [9, 29].

The behavior of the MIMSE algorithm while exploring such HS states provides
insights about their distribution in configuration space, summarized in Fig.
3b. Starting from an initial quenched configuration, MIMSE descends a finite
energy canyon until it reaches a portion of the landscape where the canyon
floor is covered by a slender, connected domain of HS configurations, as if
the canyon floor were flooded by a level set having zero energy. In this case,
for every added bias, the canyon walls funnel all FIRE paths to points on the
intersection of the edge of the bias and the edge of a zero energy domain, see
Fig. 3b. FIRE then halts when the total energy per particle first drops below
$U_{tol}$. The fact that the hyperspherical edge of every bias intersects a
zero energy domain leads us to conclude that the domain is continuous (at
least on the bias length scale used to search the configuration space).

At higher volume fractions $\phi\gtrsim 0.667$, MIMSE descends the canyon in
the energy landscape, but never reaches a zero energy or HS domain, Fig. 4b.
Instead the canyon floor contains a cluster of finite energy inherent
structures as in the SGM case (which was explored at a volume fraction well
above its $\phi_{J}$). In a sense, it is as if the entire landscape has been
raised and the river drained away, Fig. 3a. In the opposite limit, the jamming
volume fraction $\phi_{J}$ corresponds to lowering the landscape so far as to
flood the entire canyon, such that the HS domain covers even the highest ISes
on the landscape.

These findings beg the question: is the lack of HS states at high $\phi$ a
property of the landscape itself, or is MIMSE simply failing to descend far
enough along the canyon to reach them? To address this question, we repeated
our experiment using zero-energy initial configurations. Above $\phi\approx
0.666$ these were created using an adapted affine compression scheme [30, 31]
with MIMSE relaxation. When applying MIMSE to these configurations at high
volume fractions, they quickly jump out of the HS domain, finding soft sphere
ISes and suggesting that the HS domain is small compared to our bias radius.
In practice, most runs with volume fractions $\phi\gtrsim 0.669$ show mixed
behavior, switching intermittently between zero and finite energy states,
revealing HS domains of finite size somewhat larger than our bias radius. Runs
at such high volume fractions return many biased minima (as opposed to there
being zero biased minima at lower phi), suggesting a qualitative change in the
landscape structure.

Overall, this suggests that at low volume fractions, HS states form a
connected network of domains consisting of thin channels, which break up into
clusters of ISes around $\phi_{max}\simeq 0.669$. This change in the domain
connectivity would presumably correlate with changes in the relaxation
dynamics of the HS fluid. Indeed, the value of $\phi_{max}$ we find is close
to the predicted value for kinetic arrest $\phi_{0}=0.672$, obtained by
extrapolating relaxation times of the equilibrium HS fluid to a divergence
[9]. It seems likely that this breakup occurs at somewhat different values of
$\phi_{max}$ in different portions of the landscape.

##  V Kob-Andersen glass former

Lastly, to determine if the canyons we find are also a feature of the energy
landscapes of atomic glasses, we consider the Kob-Andersen (KA) model. This
consists of a binary mixture having a total of $N=256$ particles with
$x_{A}=0.8$ at a number density of $\rho_{total}=1.2$. To ensure force
continuity, we employ a quadratically-smoothed, truncated form of the Lennard-
Jones potential [32]:

$V(\mathbf{r}_{ij})=\begin{cases}k\epsilon_{ij}\left({\left(\frac{\sigma_{ij}}{\lVert{\mathbf{r}_{ij}}\rVert}\right)}^{1}2-{\left(\frac{\sigma_{ij}}{\lVert{\mathbf{r}_{ij}}\rVert}\right)}^{6}\right)\\\
{}+\nu(\lVert{\mathbf{r}_{ij}}\rVert),&\text{if
$\lVert\mathbf{r}_{ij}\rVert<\mathbf{r}_{c}$}\\\
0,&\text{otherwise},\end{cases}$ (2)

where $\epsilon_{AA}=1.0$, $\epsilon_{BB}=0.5$ and $\epsilon_{AB}=1.5$;
$\sigma_{AA}=1.0$, $\sigma_{BB}=0.88$ and $\sigma_{AB}=0.8$, and
$\nu(\lVert{\mathbf{r}_{ij}}\rVert)$ is a smoothing function (see Appendix B
for details). As before, we first choose random points in configuration space
and FIRE relax them to their first energy minimum to form a quenched ensemble
of ISes (Fig. 5a, left panel).

(a)(a)(a)[22][22] [22][22]

(a)

(a)

(a)

(b)

(a)

(c)

Figure 5: MIMSE descends the landscape of the Kob-Andersen glass former. (a)
The left panel shows the probability distribution of energies for an ensemble
of $1000$ ‘quenched’ initial configurations, formed by FIRE relaxation of
random configurations. In the middle panel we find that MIMSE successfully
reaches low energy configurations for bias radii, $\mathcal{U}_{\sigma}\approx
1$. Unlike the prior systems, here the algorithm reaches a point where it
abruptly finds no new ISes, even after a large number of additional biases are
applied. This suggests that it has reached the bottom of a metabasin and
completely covered the ISes there. The lowest energy IS or ‘basin-bottom’
energy, $U_{bb}$, found, when computed for many initial configurations, yield
a distribution shown in the right panel. (b) The cumulative density function
(CDF) [33] of the pairwise separations of the ISes near the bottom energy
indicates that ISes form a low-dimensional cluster, explained in the text. (c)
Replotting the distribution of basin bottom energies shows that it peaks at
the mode-coupling energy (gray line) and closely approaches the Kauzmann limit
(black line) [22].

Despite the likely differences between the KA and soft sphere energy
landscapes, we find that MIMSE again efficiently descends to low energy
portions of the KA landscape, Fig. 5a. This occurs in a narrower window of
bias radius, $0.9<\mathcal{U}_{\sigma}<1.3$. The upper bound value for
$\mathcal{U}_{\sigma}$ is much smaller than the length of the descent path,
confirming a narrow canyon-like structure. The KA system does require more
biases for each new IS sampled, typically $\sim 8-22$ versus $\sim 1.4$ in the
SGM case, indicating a different structure for the IS clusters. Earlier
studies using simulated annealing of the KA model [11] provide a relationship
between IS energy and fictive equilibrium temperature. Characteristic
temperatures such as the mode-coupling and Kauzmann temperatures can thus be
mapped to their corresponding IS energies [22]. Benchmarking our algorithm
against a Langevin dynamics implementation of simulated annealing finds a
$\approx 100$-fold improvement in computational speed to reach configurations
near the mode-coupling energy.

Unexpectedly, after MIMSE finds roughly $100$ Kob-Andersen ISes, the number of
biases that must be added to generate each new IS appears to increase
drastically, soon precluding the finding of any new ISes. Further, the path
starts to cross back on itself (exploring configurations covered by previous
biases), suggestive of basin-filling behavior. This suggests that the path
being followed by MIMSE effectively ends after it reaches the bottom of a
canyon or ‘meta-basin’. To rule out the possibility that the algorithm is
getting stuck due to canyon tapering, we verified that smaller biases did not
descend to lower energy.

MIMSE being unable to find new ISes after reaching the bottom of a meta-basin
suggests that the ISes there form a cluster whose hypervolume is small enough
to be covered completely by a small number of biases, and that the cluster is
surrounded by a large hypervolume that is nearly devoid of ISes. We determined
the effective dimensionality of the IS cluster at the canyon bottom by
computing its correlation dimension [33]; that is, the scaling exponent of the
cumulative distribution function of the $3N$-dimensional Euclidian distances
between all pairs of ISes. Were the cluster a fractal, the correlation
dimension would report the fractal dimension, $D_{f}$. The power-law exponent
seen for small separations Fig. 5b is consistent with the cluster being
fractal, $D_{f}\approx 1.35$; the roll-over at longer scales being due to the
cluster’s finite-size [33]. The cluster is clearly a low-dimensional object,
explaining the ease of covering it with biases.

The meta-basin architecture of the KA landscape means that for every random
starting configuration, MIMSE returns a single lowest energy minimum for the
corresponding meta-basin. Due to sparse sampling this lowest value is only an
upper-bound for the true minimum of the metabasin, and its difference from the
true minimum is presumably sensitive to details of our algorithm. Sampling
many meta-basins shows that these lowest energies are roughly Gaussian
distributed, with a standard deviation comparable to that of the quenched ISes
on this landscape, Fig. 5a. A structural analysis of these different low
energy configurations using common neighbor analysis [34, 35] show no sign of
crystallinity. An earlier study using ‘basin-hopping’ minimization [36]
yielded similar minimum energy values for amorphous structures in the same
model [37]. We note that the mean energy of these bottoms is a characteristic
energy of the landscape itself, and corresponds roughly to this model’s mode-
coupling energy, $U_{MC}\approx-7.00$, [22] see Fig. 5c. The predicted value
of the Kauzmann energy, $U_{K}$, is about $\sim 3$ standard deviations below
the mean, $U_{K}\approx-7.08$ [22]. Sampling several hundred meta-basins
yields ISes very close to but not below the Kauzmann energy, consistent with
this empirical prediction.

##  VI Discussion

We have found that our algorithm provides surprising insights into the large-
scale geometry and structure of three glassy landscapes, in particular the
routes leading to the landscapes’ low energy ISes. A cluster analysis or
trajectory ensembles shows a structural bifurcation between small and large
lengthscales. It finds that barriers between minima are common on short
length-scales, indicating that ISes are clustered densely together; while
these clusters in turn are separated by barriers on much larger lengthscales.
In between these lengthscales, there is a window where motion is essentially
unrestricted by barriers, corresponding to the smooth walled canyons we
report. The dense IS clusters along the floors of these canyons resemble the
well-known rugged and barrier-filled landscapes of glasses; our work reveals
that these rugged and low-dimensionality domains are surrounded by large,
high-dimensional domains devoid of minima and barriers, and which contain
energy gradients that lead directly to lower energies in the landscape. In the
the hard-sphere fluid, the breakup of the zero-energy domain at a threshold
volume fraction suggests an intuitive geometrical explanation for a
corresponding kinetic arrest of the fluid. In the Kob-Andersen glass former,
the algorithm is able to consistently reach the ends of the canyons. This may
be due to the smaller configuration space, a result of having both
indistinguishable and fewer particles. The fact that the lowest canyon
energies observed approach the Kauzmann limit confirms that by avoiding being
trapped, MIMSE can act like a global optimizer similar to basin hopping [36].
The presence of similar canyon-like features in these seemingly different
systems, might explain the broad range of qualitatively similar physical
properties and phenomenon in these systems [13, 38, 29].

Future work will apply MIMSE to other glassy landscapes, including bonded
model systems. Moreover, we have not yet studied the (presumed) branching of
these canyons nor the barriers between different canyons. While informative,
MIMSE does not obey detailed balance, nor does it return samples corresponding
to any canonical ensemble. We anticipate that future work will fruitfully
hybridize MIMSE with other methods, such as swap [10, 39] or ghost-particle
[40] Monte Carlo, and tempering [41] to enable barrier crossing while
obtaining more canonical sampling. Recently the loss landscape of deep neural
networks has been shown to have similar features [42] to that of a soft glassy
matter system, perhaps MIMSE will enable useful exploration of such landscapes
as well.

###### Acknowledgements.

We are grateful for useful conversations with Andrea Liu and Talid Sinno. This
work was supported by NSF-DMR 1609525 and 1720530 and computational resources
provided by XSEDE through TG-DMR150034.

##  Appendix A Details of MIMSE

A detailed description of the algorithmic procedure is provided below.

  1. 1.

A random energetic minima sampled via FIRE from a random configuration
(quenched ensemble) is selected as the starting point.

  2. 2.

A bias potential is added on to the system centered around the selected
minimum in 3N dimensional configuration space. (Note: to ensure the bias
potential is center of mass preserving during evaluation of the bias, we
subtract the mass weighted sum of forces on all particles.)

  3. 3.

The system is given a small, positional displacement in a random direction in
the configuration space, scaled such that at least one particle exceeds the
system force tolerance, $F_{tol}$. This displacement is adjusted to preserve
the center of mass location. The system is then quenched on the biased energy
landscape (${U}_{\text{tot}}={U}+{U}_{\text{i, bias}}$, where $i$ signifies
the bias number)

  4. 4.

Biases are added in succession until a new unbiased minimum is reached. This
is verified by ensuring that ${U}_{\text{bias}}/N<U_{tol}$.

  5. 5.

The above process is repeated until a desired amount of the energy landscape
is sampled.

To undertake the above process with utmost efficiency and accuracy for our
systems of interest with a highly rugged energy landscapes, we use a finite-
ranged bias potential viz. a symmetric quartic function, as described in Eq.
1.

Parameters $\mathcal{U}_{\sigma}$ and $\mathcal{U}_{0}$ represent the $3N$
dimensional Euclidean extent and energetic height of the bias respectively.
Throughout our simulation, the center of mass of the system remains at rest,
due to the bias force exerting no force on the center of mass. This is
enforced by subtracting the mass weighted total force from individual bias
force calculations. Despite this correction, our bias remains isotropic in the
$3N-1$ angular dimensions.

We select a pair of parameters that efficiently samples lower energy states
effectively and does so with a small number of biases per every new (unbiased)
minima. We use a divide and conquer like parameter search procedure to
determine a suitable pair of parameters for our algorithm. Our initial guess
for both parameters are based on the potential energy of interaction. For
example, in the SGM system, we search for an optimal $\mathcal{U}_{0}$ between
$\sim\epsilon/2$ and $\sim N^{2}\epsilon/2$. Meanwhile, a search between
$\sim\langle R\rangle$ and $\sim\sqrt{N}\langle R\rangle$. The selection of
parameters is further verified by studying the gradient and clusters formed by
trajectories emerging from a starting IS using different
$\mathcal{U}_{\sigma}$ values (see Fig. 2).

To keep track of the biases, we use a modified extended $3N$-dimensional
neighbor list that keeps track of all biases in high dimensional space around
the system configuration within a cut-off distance $\mathcal{U}_{c1}$.
Further, we have a long range secondary cut-off length $\mathcal{U}_{c2}$
which defines our neighbor search boundary; biases are retired when the system
is farther than $\mathcal{U}_{c2}$ from a particular bias. Thus, in effect we
maintain neighbors (biases) within $\mathcal{U}_{\sigma}+\mathcal{U}_{c1}$,
and update the neighbor list by searching for biases within $\mathcal{U}_{c2}$
from the system position. The update is done when the system moves an
Euclidean displacement more than $\mathcal{U}_{c1}$, whose value is chosen for
computational efficiency. Parameter $\mathcal{U}_{c2}$ is chosen so that the
ultimate trajectory of the system is not affected. It may be noted that both
of these parameters are system and bias size dependent, and need to be tuned
for each system to balance the efficiency of the method with not retiring
biases prematurely.

Our metadynamics-based algorithm differs from that developed by Yip in three
key ways: (i) We use a smooth, truncated bias potential 1 instead of the
traditional Gaussian bias, that enables the use of a high dimensional neighbor
list. This allows us to significantly reduce the computational overhead
associated with computing the total bias potential, and also facilitates the
retirement of biases beyond a certain distance. (ii) Further, we make sure
that the center of mass of the system remains at rest, throughout the
simulation. We enforce this on the bias potential and the random displacement
after bias addition, which appears to lead to more efficient descent down the
landscape. (iii) Last, we use an optimized minimzer - FIRE [20], that helps
boost the computational efficiency while exploring such glassy high-
dimensional landscapes.

These modifications and additions allow us to overcome the challenges posed by
our high dimensional system. It must be noted that each system responds
differently to the biases in a way that depends on the characteristic energy
and length scales of that particular system.

##  Appendix B  Model Details

Figure 6:  Energy descent versus bias radius $\mathcal{U}_{\sigma}$ confirms
that SGM landscape canyons are tapered. The optimal choice of
$\mathcal{U}_{\sigma}=3$ descends to significantly lower energy states.
Smaller biases also descend, but more slowly. Larger biases descend rapidly at
first, but plateau at an energy when they have become too big for the
narrowing canyon.

SGM Model: The Soft Glassy material system is a set of $N=342$ particles
interacting via

$V(\mathbf{r}_{ij})=\begin{cases}\frac{\epsilon}{\alpha}{\left(1-\frac{\lVert{\mathbf{r}_{ij}}\rVert}{R_{i}+R_{j}}\right)}^{\alpha},&\text{if
$\lVert\mathbf{r}_{ij}\rVert<R_{i}+R_{j}$}\\\ 0,&\text{otherwise},\end{cases}$
(3)

with $\alpha=2$. The particles have a radii distribution resembling a Weibull
distribution, $\textbf{P}(R)\sim(k/\lambda)(x/\lambda)^{k-1}$, where
$k=1.66,\lambda=0.78$. The exact distribution snapshot is taken from a model
foam simulation [7]. To ensure meaningful descent, we further tweak our
algorithm so that the biases only act on non-rattler particles (those having a
coordination $z>d+1=4$). We use tolerance values of $U_{tol}\sim 10^{-12}$ and
$F_{tol}\sim 10^{7}$ for simulations for this system. As noted in the main
text, for slightly too large values of $\mathcal{U}_{\sigma}$ the algorithm
gets stuck, plateauing at intermediate energy, Fig. 6. This confirms the
interpretation that the canyon is tapered as it descends to lower energies.

Figure 7:  We use contact coordination, $\langle z\rangle$, as an order
parameter for hard sphere configurations (see Fig. 4a). We verified that this
order parameter is physically sensible by finding that its value is
insensitive to the choice of energy tolerance, $U_{tol}$.

Cluster Analysis: We use a single-linkage clustering analysis [26] to analyze
the linkage between the different $\hat{u}$ obtained through the $1100$
($>3N$) different trajectories. We use an angle of $\theta_{0}\sim
1.31\text{rad}$ as the linkage threshold to classify any two points as part of
the same cluster. We see that for other angles greater and smaller than
$\theta_{0}$, the ISes form a single or large number of clusters. Sweeping
through the angles, we find that at $\theta_{0}$, the system shows a clear
trend for the number of clusters detected (see Fig. 2b). Meanwhile, we ignore
any outliers detected as small or singleton clusters in the analysis using the
following criteria: $n_{cluster}\leq 0.25\times N_{points}/N_{cluster}$ [27].
It may also be noted that the bias potential acts on all particles unlike a
typical SGM MIMSE run.

HS glass: To study this system, we use a well-studied [9, 10, 28] model of
hard-sphere (HS) glasses, having a size distribution given by:
$\textbf{P}(R)\sim
R^{-3},R\in[R_{\text{min}},R_{\text{max}}],R_{\text{min}}/R_{\text{max}}=0.4492$,
and $N=1000$ particles; and $R_{\text{max}}$ is taken as the reference $=1$.

To compare our results with earlier work [9], we use a similar tolerance value
of $U_{tol}=10^{-16}$. We further use a $F_{tol}\sim 10^{-7}$ to determine
positional displacements in MIMSE. MIMSE works efficiently on the HS portion
and the soft sphere extension of the landscape. All results (Fig. 4) discussed
and shown are normalized with respect to the $\phi_{ref}=0.66$ system. We use
$\mathcal{U}_{\sigma,ref}=4$ and $\mathcal{U}_{0,ref}=0.1$ or $1$ as reference
bias dimensions and scale them with changes in $\phi$. The scaling relation
used to determine bias parameters at $\phi$ are:
$\mathcal{U}_{0}=\frac{\mathcal{U}_{0},ref}{U_{IS,ref}}U_{IS}$, and
$\mathcal{U}_{\sigma}=\frac{\mathcal{U}_{\sigma}}{\phi^{(1/3)}}\phi_{ref}^{(1/3)}$;
where $U_{IS}$ represents the IS energy obtained through random quenching at
$\phi$. To study the changing landscape with $\phi$, we use MIMSE on this
landscape starting from $2$ kinds of starting points. When starting from a
random quenched IS on the soft-sphere portion of the landscape, we calculate
the average $U/N$ of the unbiased ISes found over $100$ consecutive applied
biases, starting after evolving the system by applying $300$ biases.
Meanwhile, when a compression generated HS IS is used as starting point, we
calculate the average $U/N$ of the unbiased ISes found over $90$ consecutive
applied biases, starting after evolving the system by applying $10$ biases.
This operational measure serves as an indicator of a transition in the density
and distribution of the hypervolume of HS domains in the landscape.

We use the average coordination number $\langle z\rangle$ as an order
parameter to study the HS state exploration facilitated by MIMSE. We find that
the effective $\langle z\rangle$ values remain independent of $U_{tol}$, see
Fig. 7, suggesting that it is a reliable feature of the configuration and not
the algorithm.

To create HS states at high volume fraction we combine MIMSE with an affine
compression and relaxation scheme by Weeks [31]. Specifically, we run the
compression scheme with FIRE as a minimizer until the system jams and is
unable to reach higher volume fractions. We then run MIMSE on that jammed
configuration, reducing its contact coordination significantly, as in Fig. 4a.
Rerunning the compression scheme then allows a still higher volume fraction to
be reached before the system again jams. This alternating compression and
unjamming (using MIMSE) process was repeated to attain HS states up to
$\phi=0.670$ Fig. 4b.

(a)(a)

(a)

(a)

(a)

(b)

Figure 8:  Fractal scaling analysis of the configuration space path on the SGM
landscape: (a) As explained in the text, departures from a single power-law
indicate a multi-fractal scaling with a range of fractal dimensions $2.2-2.7$.
(b) The probability distribution of the number of biases that must be applied
between two consecutive ISes appears insensitive to $\mathcal{U}_{\sigma}$,
suggesting a self-similar rather than uniform distribution of ISes along the
MIMSE descent path.

KA glass former: The Kob-Andersen (KA) model glass is a binary system with a
total of $N=256$ particles with number fraction $x_{A}=0.8$ at a number
density of $\rho_{total}=1.2$ [11]. To ensure force continuity, we employ a
quadratically-smoothed, truncated form of the Lennard-Jones potential [32]. We
use a tolerance value of $U_{tol}\sim 10^{-12}$ and $F_{tol}\sim 10^{-7}$ to
define IS minima in this model.

Simulated Annealing: To compare MIMSE’s ability to explore and generate low-
energy states, we employ Langevin Dynamics based simulated annealing. The
system is first equilibrated at $k_{B}T=1.0-1.5$, and then systematically
cooled using a linear cooling ramp. The theraml configurations are
periodically quenched using FIRE to obtain corresponding ISes for these
thermal states. Our results reproduce those in a previous study [11].
Comparing the computational cost to generate low energy ISes using the two
different methods shows $\approx 100$ fold speed-up when using MIMSE at final
energy values of $U/N\approx-7.00$.

##  Appendix C Configuration path fractal analysis

We analyzed the geometry of the high-dimensional path followed by MIMSE on the
SGM landscape, using methods from an earlier study [7]. For each pair of ISes
on the path both their Euclidean distance, $\Delta r^{2}$, and the contour
length $\Delta s$ of the path between them are first computed, then pooled
together, partitioned in $\Delta s$ and averaged. This yields the scaling of
Euclidean distance with contour length, analogous to a mean-squared
displacement in lower-dimensions. For a fractal curve, these quantities obey a
power-law relationship, $\Delta r^{2}\sim\Delta s^{(2/D_{f})}$, where $D_{f}$
is the fractal dimension. Fig. 8a shows the results for the configuration
space path on the SGM landscape for different bias parameters, which resembles
a crossover between two power-laws, indicating only rough self-similarity. At
the longest scales, the configuration path scaling suggests
$D_{f}^{MSD}\approx 2.2$. This value is slightly higher than that for a random
walk $D_{f}=2$, suggesting that the canyon MIMSE is following is slightly more
confined than a random walk.

The small scale tortuosity of the configuration path explains the power-law
dependence of the descent rate on $\mathcal{U}_{\sigma}$, Fig. 1b (inset). If
MIMSE essentially follows the the same path with different-sized biases, more
smaller biases would be required to ‘cover’ a tortuous path. This is analogous
to another method for computing fractal dimension, which consists of counting
the number of boxes required to cover a curve as a function of the box size.
We suppose that the scaling exponent seen in Fig. 1b(inset) effectively
provides such a ‘box-counting’ measure of the path fractal dimension,
$D_{f}^{bias}\approx 2.8$ which is in fact consistent with $D_{f}^{MSD}\approx
2.7$ at that length-scale.

As the algorithm proceeds we can track the number of new biases $n$ required
to produce each unbiased IS. Unexpectedly, both $\langle n\rangle\sim 1.4$ and
the distribution $P(n)$ are similar and insensitive to the bias radius
$\mathcal{U}_{\sigma}$, Fig. 8b. If the ISes were distributed uniformly along
the path, larger biases would be more likely to find biased minima, and
$\langle n\rangle$ would increase with $\mathcal{U}_{\sigma}$. This
insensitivity can be naturally explained by structure of the IS cluster being
self-similar on these length scales.

## References

  * Wales  _et al._ [2007] D. J. Wales, J. P. K. Doye, M. A. Miller, P. N. Mortenson, and T. R. Walsh,  _Energy Landscapes: From Clusters to Biomolecules_, Vol. 115 (John Wiley & Sons, Ltd, 2007) pp. 1–111.
  * Stillinger and Weber [1982] F. H. Stillinger and T. A. Weber, Hidden structure in liquids, Physical Review A 25, 978 (1982).
  * Berthier [2011] L. Berthier, Dynamic Heterogeneity in Amorphous Materials, Physics 4 (2011).
  * Massen and Doye [2007] C. P. Massen and J. P. Doye, Power-law distributions for the areas of the basins of attraction on a potential energy landscape, Phys. Rev. E 75, 10.1103/PhysRevE.75.037101 (2007), arXiv:0509185 [cond-mat] .
  * Xu  _et al._ [2011] N. Xu, D. Frenkel, and A. J. Liu, Direct Determination of the Size of Basins of Attraction of Jammed Solids, Physical Review Letters 106, 1 (2011), arXiv:1101.5879 .
  * Ashwin  _et al._ [2012] S. S. Ashwin, J. Blawzdziewicz, C. S. O’Hern, and M. D. Shattuck, Calculations of the structure of basin volumes for mechanically stable packings, Phys. Rev. E 85, 061307 (2012), arXiv:1112.4234 .
  * Hwang  _et al._ [2016] H. Hwang, R. Riggleman, and J. Crocker, Understanding soft glassy materials using an energy landscape approach, Nature Materials 15, 1031 (2016).
  * Liu and Nagel [2010] A. J. Liu and S. R. Nagel, The jamming transition and the marginally jammed solid, Annual Review of Condensed Matter Physics 1, 347 (2010), https://doi.org/10.1146/annurev-conmatphys-070909-104045 .
  * Berthier  _et al._ [2016] L. Berthier, D. Coslovich, A. Ninarello, and M. Ozawa, Equilibrium Sampling of Hard Spheres up to the Jamming Density and Beyond, Phys. Rev. Lett. 116, 1 (2016), arXiv:1511.06182 .
  * Ozawa  _et al._ [2017] M. Ozawa, L. Berthier, and D. Coslovich, Exploring the jamming transition over a wide range of critical densities, SciPost Physics 3, 1 (2017), arXiv:1705.10156 .
  * Sastry  _et al._ [1998] S. Sastry, P. G. Debenedetti, and F. H. Stillinger, Signatures of distinct dynamical regimes in the energy landscape of a glass-forming liquid, Nature 393, 554 (1998).
  * Charbonneau  _et al._ [2014] P. Charbonneau, J. Kurchan, G. Parisi, P. Urbani, and F. Zamponi, Fractal free energy landscapes in structural glasses, Nature Communications 5, 1 (2014), arXiv:1404.6809 .
  * Royall  _et al._ [2018] C. P. Royall, F. Turci, S. Tatsumi, J. Russo, and J. Robinson, The race to the bottom: approaching the ideal glass?, Journal of Physics: Condensed Matter 30, 363001 (2018), arXiv:1711.04739 .
  * Cerjan and Miller [1981] C. J. Cerjan and W. H. Miller, On finding transition states, The Journal of Chemical Physics 75, 2800 (1981), https://doi.org/10.1063/1.442352 .
  * Wales [2018] D. J. Wales, Exploring Energy Landscapes, Annual Review of Physical Chemistry 69, 401 (2018), pMID: 29677468, https://doi.org/10.1146/annurev-physchem-050317-021219 .
  * Munro and Wales [1999] L. J. Munro and D. J. Wales, Defect migration in crystalline silicon, Phys. Rev. B 59, 3969 (1999).
  * Middleton  _et al._ [2001] T. F. Middleton, J. Hernández-Rojas, P. N. Mortenson, and D. J. Wales, Crystals of binary lennard-jones solids, Phys. Rev. B 64, 184201 (2001).
  * Laio and Parrinello [2002] A. Laio and M. Parrinello, Escaping free-energy minima, Proceedings of the National Academy of Sciences 99, 12562 (2002), https://www.pnas.org/doi/pdf/10.1073/pnas.202427399 .
  * Kushima  _et al._ [2009] A. Kushima, X. Lin, J. Li, J. Eapen, J. C. Mauro, X. Qian, P. Diep, and S. Yip, Computing the viscosity of supercooled liquids, The Journal of Chemical Physics 130, 10.1063/1.3139006 (2009), https://doi.org/10.1063/1.3139006 .
  * Bitzek  _et al._ [2006] E. Bitzek, P. Koskinen, F. Gähler, M. Moseler, and P. Gumbsch, Structural Relaxation Made Simple, Phys. Rev. Lett. 97, 1 (2006).
  * Kob and Andersen [1995] W. Kob and H. C. Andersen, Testing mode-coupling theory for a supercooled binary Lennard-Jones mixture I: The van Hove correlation function, Phys. Rev. E 51, 4626 (1995).
  * Angell  _et al._ [1999] C. A. Angell, B. E. Richards, and V. Velikov, Simple glass-forming liquids: their definition, fragilities, and landscape excitation profiles, Journal of Physics: Condensed Matter 11, A75 (1999).
  * Sollich  _et al._ [1997] P. Sollich, F. Lequeux, P. Hébraud, and M. E. Cates, Rheology of Soft Glassy Materials, Phys. Rev. Lett. 78, 2020 (1997).
  * Sollich [1998] P. Sollich, Rheological constitutive equation for a model of soft glassy materials, Phys. Rev. E 58, 738 (1998).
  * Sollich and Cates [2012] P. Sollich and M. E. Cates, Thermodynamic interpretation of soft glassy rheology models, Phys. Rev. E 85, 031127 (2012).
  * Gower and Ross [1969] J. C. Gower and G. J. S. Ross, Minimum Spanning Trees and Single Linkage Cluster Analysis, Journal of the Royal Statistical Society. Series C (Applied Statistics) 18, 54 (1969).
  * Christy  _et al._ [2015] A. Christy, G. M. Gandhi, and S. Vaithyasubramanian, Cluster Based Outlier Detection Algorithm for Healthcare Data, Procedia Computer Science 50, 209 (2015), big Data, Cloud and Computing Challenges.
  * Berthier  _et al._ [2017] L. Berthier, P. Charbonneau, D. Coslovich, A. Ninarello, M. Ozawa, and S. Yaida, Configurational entropy measurements in extremely supercooled liquids that break the glass ceiling, Proceedings of the National Academy of Sciences of the United States of America 114, 11356 (2017), https://www.pnas.org/doi/pdf/10.1073/pnas.1706860114 .
  * Ediger  _et al._ [2021] M. D. Ediger, M. Gruebele, V. Lubchenko, and P. G. Wolynes, Glass Dynamics Deep in the Energy Landscape, The Journal of Physical Chemistry B 125, 9052 (2021), pMID: 34357766, https://doi.org/10.1021/acs.jpcb.1c01739 .
  * Xu  _et al._ [2005] N. Xu, J. Blawzdziewicz, and C. S. O’Hern, Random close packing revisited: Ways to pack frictionless disks, Phys. Rev. E 71, 061306 (2005).
  * Desmond and Weeks [2009] K. W. Desmond and E. R. Weeks, Random close packing of disks and spheres in confined geometries, Phys. Rev. E 80, 051305 (2009), arXiv:0903.0864 .
  * Stoddard and Ford [1973] S. D. Stoddard and J. Ford, Numerical Experiments on the Stochastic Behavior of a Lennard-Jones Gas System, Phys. Rev. A 8, 1504 (1973).
  * Grassberger and Procaccia [1983] P. Grassberger and I. Procaccia, Characterization of Strange Attractors, Phys. Rev. Lett. 50, 346 (1983).
  * Stukowski [2012] A. Stukowski, Structure identification methods for atomistic simulations of crystalline materials, Modelling and Simulation in Materials Science and Engineering 20, 045021 (2012).
  * Stukowski [2009] A. Stukowski, Visualization and analysis of atomistic simulation data with OVITO–the open visualization tool, Modelling and Simulation in Materials Science and Engineering 18, 015012 (2009).
  * Wales and Doye [1997] D. J. Wales and J. P. K. Doye, Global optimization by basin-hopping and the lowest energy structures of lennard-jones clusters containing up to 110 atoms, The Journal of Physical Chemistry A 101, 5111 (1997), https://doi.org/10.1021/jp970984n .
  * Middleton and Wales [2001] T. F. Middleton and D. J. Wales, Energy landscapes of some model glass formers, Phys. Rev. B 64, 024205 (2001).
  * Debenedetti and Stillinger [2001] P. G. Debenedetti and F. H. Stillinger, Supercooled liquids and the glass transition, Nature 410, 259 (2001).
  * Parmar  _et al._ [2020] A. D. Parmar, B. Guiselin, and L. Berthier, Stable glassy configurations of the Kob-Andersen model using swap Monte Carlo, The Journal of Chemical Physics 153, 134505 (2020), arXiv:https://doi.org/10.1063/5.0020208 .
  * Turci  _et al._ [2017] F. Turci, C. P. Royall, and T. Speck, Nonequilibrium phase transition in an atomistic glassformer: The connection to thermodynamics, Phys. Rev. X 7, 031028 (2017).
  * Hansmann [1997] U. H. Hansmann, Parallel tempering algorithm for conformational studies of biological molecules, Chemical Physics Letters 281, 140 (1997).
  * Chen  _et al._ [2022] G. Chen, C. K. Qu, and P. Gong, Anomalous diffusion dynamics of learning in deep neural networks, Neural Networks 149, 18 (2022).