We develop computational tools for spectral analysis of stochastic networks representing energy landscapes of atomic and molecular clusters. Physical meaning and some properties of eigenvalues, left and right eigenvectors, and eigencurrents are discussed. We propose an approach to compute a collection of eigenpairs and corresponding eigencurrents describing the most important relaxation processes taking place in the system on its way to the equilibrium. It is suitable for large and complex stochastic networks where pairwise transition rates, given by the Arrhenius law, vary by orders of magnitude. The proposed methodology is applied to the network representing the Lennard-Jones-38 cluster created by Wales's group. Its energy landscape has a double funnel structure with a deep and narrow face-centered cubic funnel and a shallower and wider icosahedral funnel. However, the complete spectrum of the generator matrix of the Lennard-Jones-38 network has no appreciable spectral gap separating the eigenvalue corresponding to the escape from the icosahedral funnel. We provide a detailed description of the escape process from the icosahedral funnel using the eigencurrent and demonstrate a superexponential growth of the corresponding eigenvalue. The proposed spectral approach is compared to the methodology of the Transition Path Theory. Finally, we discuss whether the Lennard-Jones-38 cluster is metastable from the points of view of a mathematician and a chemical physicist, and make a connection with experimental works.

Modeling by means of stochastic networks has become widespread in chemical physics. Efficient methods for the conversion of energy landscapes into stochastic networks have been developed. Wales's group constructed a large number of networks representing energy landscapes of atomic and molecular clusters and proteins.1–3 A number of researchers developed efficient techniques for constructing continuous-time Markov chains representing coarse grained dynamics of biomolecules and protein folding called Markov State Models.4–11 

Stochastic networks are attractive for analysis. On one hand, they are more mathematically tractable than the original continuous systems. On the other hand, they are designed to preserve important features of the underlying systems. Nevertheless, analysis of large and complex networks still presents a challenge. The numbers of states (vertices) and edges can be of the order of 10n, n = 3, 4, 5, 6, …, and the pairwise transition rates may vary by tens of orders of magnitude. The study of such networks motivates development of new techniques and leads to a new level of understanding of complex systems.

In this work, we focus on developing a methodology for computing eigenvalues, eigenvectors, and eigencurrents for networks representing energy landscapes whose pairwise transition rates are given by the Arrhenius law. The spectral decomposition of the generator matrix of a stochastic network exposes all of the relaxation processes taking place in the system on its way to the equilibrium. The concept of metastability and its relationship with the spectral structure has received a lot of attention in the scientific literature.12–19 The assumption of the existence of a spectral gap separating a few smallest eigenvalues allows us to simplify the dynamics of a complex stochastic system by combining its states into quasi-invariant subsets and considering transitions between them. As the temperature tends to zero, the eigenvectors can be approximated by the indicator functions of those quasi-invariant subsets, while the eigenvalues give the exit rates from them.12,13

We have recently shown that, under certain nonrestrictive assumptions, the zero-temperature asymptotics for eigenvalues and eigenvectors of the generator matrices of such networks can be found recursively in the increasing order of the absolute values of the eigenvalues.20 This recursion can be realized in the form of an efficient algorithm which, in a nutshell, is a procedure for removing edges one-by-one in a certain order from the minimum spanning tree for the given network.20 Each asymptotic eigenpair is associated with a particular state of the stochastic network. This state has the smallest potential value out of those where the asymptotic eigenvector is nonzero.

Here, we upgrade the zero temperature asymptotic analysis with a continuation technique for computing a collection of eigenpairs of interest at a finite temperature range. It is based on the Rayleigh Quotient iteration supplemented with some precautions. Each eigenpair defines a relaxation process from a perturbed probability distribution where the perturbation, given by the left eigenvector, decays uniformly throughout the network, and the decay rate is given by the eigenvalue. The quantitative description of the relaxation process for each eigenpair is given by the corresponding eigencurrent.21 We provide its physical interpretation and prove some of its properties. Eigencurrents can be readily computed from the corresponding finite temperature eigenvectors. At each temperature, the collection of states can be partitioned into those where the eigencurrent is emitted and those where it is absorbed. The collection of edges connecting emitting and absorbing states can be viewed as the true transition state for the relaxation process. We will call it the emission-absorption cut.

It is instructive to compare the key concepts of the spectral analysis with those of the Transition Path Theory (TPT).22–25 We will discuss similarities and differences between the escape process quantified by the left and right eigenvectors, the eigenvalue, and the eigencurrent and the Transition Path Process21,25 described by the committor, the transition rate, and the reactive current.

We apply the proposed methodology to the Lennard-Jones-38 (LJ38) network with 71 887 states (local minima) and 119 853 edges (transition states) created by Wales's group and available on the web.26 The LJ38 cluster has been profoundly studied by a number of researchers and used as a test problem for various newly developed techniques.20,21,27–33 The LJ38 cluster is remarkable due to the fact that it is the smallest cluster whose global minimum is achieved at a configuration based on other than icosahedral packing.34 Its energy landscape has a double-funnel structure.2,3,27 The deep and narrow face-centered cubic (fcc) funnel ends at the global minimum, the face-centered cubic truncated octahedron with the point group Oh (Fig. 1, left). We will denote it by FCC. The wider and shallower icosahedral funnel is crowned by the second lowest minimum, a fivefold rotationally symmetric incomplete icosahedron with the point group C5v (Fig. 1, right). We will denote it by ICO. At temperatures close to zero, the system spends most of the time at the fcc funnel, while at higher temperatures, the icosahedral funnel with the greater configurational entropy becomes more preferable.2,28,35 The LJ38 stochastic network is a feasible model for the low-temperature dynamics of the LJ38 cluster.

FIG. 1.

Left: the global minimum of the LJ38 cluster, the face-centered cubic truncated octahedron with the point group Oh. Right: the second lowest minimum of the LJ38 cluster, an incomplete icosahedron with the fivefold rotational symmetry (the point group C5v).

FIG. 1.

Left: the global minimum of the LJ38 cluster, the face-centered cubic truncated octahedron with the point group Oh. Right: the second lowest minimum of the LJ38 cluster, an incomplete icosahedron with the fivefold rotational symmetry (the point group C5v).

Close modal

The double-funnel structure of the energy landscape might make us expect that the smallest nonzero eigenvalue of the generator matrix of the LJ38 network corresponds to the escape process from the icosahedral funnel. However, we have found that the eigenvalue corresponding to this escape process is only 245th smallest in the absolute value at low temperatures.20 Moreover, there is no appreciable spectral gap separating it from the rest. There are a few notable gaps separating eigenvalues corresponding to the escape processes from some isolated high-lying liquid-like states over extremely high potential barriers. The rest of the spectrum is pretty dense, and gaps appear in it only at temperatures very close to zero. Roughly speaking, the spectrum is polluted by a large number of eigenvalues corresponding to the escapes from isolated high-lying liquid-like states, while these states are dynamically accessible from icosahedral and fcc states, only if the waiting time is extremely long. If we bound the observation time36 which is equivalent to imposing a cap on the potential energy, the expected spectral gap separating the eigenvalue, responsible for the escape from the icosahedral funnel, appears.

We have computed a collection of finite temperature eigenpairs for the LJ38 network corresponding to the escape processes from the largest quasi-invariant sets. The eigenvalue corresponding to the escape process from the icosahedral funnel exhibits a superexponential growth. A quantitative description of this escape process at different values of temperature is provided by the corresponding eigencurrent. The emission distribution of the eigencurrent together with the distribution of the eigencurrent in the emission-absorption cut explain the observed superexponential temperature dependence of the escape rate from the icosahedral funnel.

Finally, we discuss the issue of metastability of LJ38. The absence of the spectral gap renders this question nontrivial from the mathematical point of view. On the other hand, the set of icosahedral configurations is certainly metastable from the point of view of a chemical physicist. We undertake an attempt to make a connection between our results obtained for the LJ38 network and the experimental results where rare gas clusters such as argon, krypton, and xenon have been studied using electron diffraction37–43 or x-rays.44 

The rest of the paper is organized as follows. Section II is devoted to a physical interpretation and some properties of eigenvalues, left and right eigenvectors, and eigencurrents of the generator matrices of networks with detailed balance. In Sec. III, we give an overview of zero-temperature asymptotic estimates for spectra of stochastic networks. Numerical algorithms are presented in Sec. IV. Section V is devoted to the application to the LJ38 network. The analyses of the LJ38 network by means of the proposed spectral approach and the Transition Path Theory21–25 are compared in Sec. VI. Metastability issues and a connection with experimental works are discussed in Sec. VII. Conclusions are summarized in Sec. VIII.

In this section, we discuss what can we learn from the spectral decomposition of the generator matrix describing the dynamics of a stochastic network.

In this section, we provide some background and introduce some notations. Let G(S, E, L) be a network, where S and E are its sets of states and edges, respectively. We assume that the number of states is finite and denote it by N. The dynamics of this network is described by the generator matrix

$L =\lbrace L_{ij}\rbrace _{i,j=1}^N$
L={Lij}i,j=1N⁠. If states i and j (ij) are connected by an edge, then Lij is the pairwise transition rate from i to j, otherwise Lij = 0. The diagonal entries Lii are defined so that the row sums of the matrix L are zeros. This fact implies that L has a zero eigenvalue λ0 = 0 whose right eigenvector can be chosen to be e ≔ [1, 1, …, 1]T, i.e., Le = 0. The corresponding left eigenvector can be chosen to be the equilibrium probability distribution π = [π1, π2, …, πN],

\begin{equation*}\pi ^TL=0,\quad \sum _{i=1}^N\pi _i= 1.\end{equation*}
πTL=0,i=1Nπi=1.

The generator matrices of networks representing energy landscapes possess the detailed balance property, i.e., πiLij = πjLji, which means that on average, there is the same numbers of transitions from i to j and from j to i per unit time. The detailed balance implies that the matrix L can be decomposed as

\begin{equation*}L=P^{-1}Q,\end{equation*}
L=P1Q,

where P is the diagonal matrix

\begin{equation*}P={\rm diag}\lbrace \pi _1,\pi _2,\ldots ,\pi _N\rbrace ,\end{equation*}
P= diag {π1,π2,...,πN},

and Q is symmetric. Hence, L is similar to a symmetric matrix

\begin{equation}L_{sym}:=P^{1/2}LP^{-1/2}=P^{-1/2}QP^{-1/2}.\end{equation}
Lsym:=P1/2LP1/2=P1/2QP1/2.
(1)

Using Eq. (1) and the fact that Lii = −∑jiLij where Lij ⩾ 0, one can show that eigenvalues of L are real and nonpositive. Furthermore, we assume that the network is connected. Hence, L is irreducible and the zero eigenvalue has algebraic multiplicity one. We will denote the nonzero eigenvalues of L by −λk, k = 1, 2, …, N − 1, and order them so that

\begin{equation*}0<\lambda _1\le \lambda _2\le \ldots \le \lambda _{N-1}.\end{equation*}
0<λ1λ2...λN1.

The eigendecomposition of a generator matrix L leads to a nice representation of the time evolution of the probability distribution. The matrix L can be written as

\begin{equation}L=\Phi \Lambda \Phi ^TP,\end{equation}
L=ΦΛΦTP,
(2)

where Φ = [ϕ0, ϕ1, …, ϕN − 1] is a matrix of right eigenvectors, normalized so that P1/2Φ is orthogonal. Note that P1/2Φ is the matrix of eigenvectors of Lsym. Therefore,

\begin{align}\Vert P^{1/2}\phi ^k\Vert ^2&=\sum _{j=1}^N\pi _j|\phi ^k_j|^2=1,\end{align}
P1/2ϕk2=j=1Nπj|ϕjk|2=1,
(3)
\begin{align}\sum _{j=1}^N\pi _j\phi ^k_j\phi ^l_j& = 0,\quad k\ne l,\end{align}
j=1Nπjϕjkϕjl=0,kl,
(4)

and in particular, ϕ0 = e and Pϕ0 = π.

The probability distribution evolves according to the Forward Kolmogorov (also known as the Fokker-Planck) equation

\begin{equation*}\frac{dp}{dt}=L^Tp,\quad p(0)=p_0,\end{equation*}
dpdt=LTp,p(0)=p0,

where p0 is the initial distribution. Using Eqs. (2)–(4) we get

\begin{equation}p(t) = \sum _{k=0}^{N-1}c_ke^{-\lambda _kt}P\phi ^k, \ \text{ where } \ c_k=(\phi ^k)^Tp_0.\end{equation}
p(t)=k=0N1ckeλktPϕk,whereck=(ϕk)Tp0.
(5)

Note that

$c_0=e^Tp_0=\sum _{j=0}^N(p_0)_j=1$
c0=eTp0=j=0N(p0)j=1⁠. Equation (5) shows that the k-th eigen-component of p(t) remains significant only on the time interval
$O(\lambda _k^{-1})$
O(λk1)
. Eventually, all eigen-components, except for the zeroth, decay. Hence, p(t) → π as t → ∞.

For k > 0, the left and right kth eigenvectors of L, Pϕk, and ϕk, respectively, can be understood from a recipe for preparing the initial probability distribution so that only the coefficients c0 and ck in Eq. (5) are nonzero. Imagine n ≫ 1 particles distributed in the stochastic network according to the equilibrium distribution π. Since

$(\phi ^0)^TP\phi ^k=\sum _{j=1}^N\pi _j\phi _j^k\break =0$
(ϕ0)TPϕk=j=1Nπjϕjk=0 for any k > 0, the set of states S can be divided into two parts,

\begin{align}S_{+}^k &:= \lbrace i\in S : \phi ^k_i\ge 0\rbrace \nonumber\\[-5pt]\hspace*{-10pc} {\rm and} \hspace*{5.9pc} \\[-5pt]S^k_{-}&:= \lbrace i\in S : \phi ^k_i<0\rbrace . \nonumber\end{align}
S+k:={iS:ϕik0} and Sk:={iS:ϕik<0}.
(6)

In order create a component of the initial distribution parallel to the left eigenvector Pϕk, we pick some α satisfying

$\alpha \phi ^k_j\break \le 1$
αϕjk1 for all
$j\in S^k_{-}$
jSk
, remove
$\alpha n \pi _j|\phi ^k_j|$
αnπj|ϕjk|
particles from each state
$j\in S^k_{-}$
jSk
and distribute them in
$S^k_{+}$
S+k
so that each state
$j\in S^k_{+}$
jS+k
gets
$\alpha n\pi _j\phi ^k_j$
αnπjϕjk
particles. The resulting distribution is

\begin{equation*}p_0 = \pi +\alpha P\phi ^k.\end{equation*}
p0=π+αPϕk.

Therefore, the left eigenvector αPϕk is a perturbation of the equilibrium distribution decaying uniformly across the network with the rate given by the corresponding eigenvalue λk. (The keyword here is “uniformly.”) The corresponding right eigenvector αϕk for each state jS gives the factor by which it is over- or underpopulated with respect to the equilibrium distribution π, as

\begin{equation*}\alpha \phi ^k_j \equiv \frac{ (p_0)_j - \pi _j}{\pi _j}.\end{equation*}
αϕjk(p0)jπjπj.

Therefore, ϕk is, in essence, a fuzzy signed indicator function of the perturbation Pϕk.

Probability currents are important tools for the quantitative description of transition processes in the system.45 For example, the reactive current is one of the key concepts of the Transition Path Theory.21–24 In the context of spectral analysis, Vanden-Eijnden proposed to consider eigencurrents.21,46 While eigenvectors determine the perturbations to the equilibrium distribution decaying with the rates given by the corresponding eigenvalues, eigencurrents give a quantitative description of the escape process from these perturbed distributions.

Eigencurrents are defined as follows. The time derivative of the ith component of the probability distribution p(t) is

\begin{equation}\frac{dp_i}{dt}=\sum _{j=1}^NL_{ji}p_j=\sum _{j\ne i}(L_{ji}p_j-L_{ij}p_i).\end{equation}
dpidt=j=1NLjipj=ji(LjipjLijpi).
(7)

Plugging in expressions for pi and pj from Eq. (5) into Eq. (7) and using the detailed balance property πiLij = πjLji, we obtain

\begin{equation}\frac{dp_i}{dt}=\sum _{k=0}^{N-1}c_ke^{-\lambda _kt}\sum _{j\ne i}\pi _iL_{ij}[\phi ^k_j-\phi ^k_i].\end{equation}
dpidt=k=0N1ckeλktjiπiLij[ϕjkϕik].
(8)

The collection of numbers

\begin{equation}F^k_{ij}:=\pi _iL_{ij}e^{-\lambda _kt}[\phi ^k_i-\phi ^k_j]\end{equation}
Fijk:=πiLijeλkt[ϕikϕjk]
(9)

is called the eigencurrent associated with the kth eigenpair. (Here, we have switched the sign and incorporated the factor

$e^{-\lambda _kt}$
eλkt into the definition in comparison with the one in Ref. 21 in order to make its physical sense more transparent.) In terms of the eigencurrent
$F^k_{ij}$
Fijk
, Eq. (8) can be rewritten as

\begin{equation}\frac{dp_i}{dt}=-\sum _{k=0}^{N-1}c_k\sum _{j\ne i}F^k_{ij}.\end{equation}
dpidt=k=0N1ckjiFijk.
(10)

Hence, the eigencurrent

$F^k_{ij}$
Fijk is the net probability current of the kth perturbation Pϕk along the edge (i, j) per unit time. In other words, if the system is originally distributed according to p0 = π + αPϕk, then the current
$\alpha F^k_{ij}$
αFijk
gives the difference of the average numbers of transitions from i to j and from j to i per unit time.

The eigencurrent

$F^k_{ij}$
Fijk can be compared to the reactive current
$F^{R}_{ij}$
FijR
in Ref. 21 (which is the same as the quantity
$f_{ij}^{AB}-f_{ji}^{AB}$
fijABfjiAB
in Ref. 24). The reactive current
$F^R_{ij}$
FijR
is conserved at every node of the network24 except for the specially designated subsets of source states A and sink states B. Contrary to it, the eigencurrent is either emitted or absorbed in all states i where
$\phi ^k_i\ne 0$
ϕik0
. Indeed, for any iS we have

\begin{align}\sum _{j\ne i}F^k_{ij}=\, &\pi _ie^{-\lambda _kt}\left[\phi _i^k\sum _{j\ne i}L_{ij}-\sum _{j\ne i}L_{ij}\phi _j^k\right] \nonumber \\=\, &\pi _ie^{-\lambda _kt}\left[-\sum _{j=1}^NL_{ij}\phi _j^k\right]=e^{-\lambda _kt}\lambda _k\pi _i\phi _i^k.\end{align}
jiFijk=πieλktϕikjiLijjiLijϕjk=πieλktj=1NLijϕjk=eλktλkπiϕik.
(11)

Therefore, every state i emits or absorbs (depending on the sign of

$\phi ^k_i$
ϕik⁠)
$e^{-\lambda _kt}\lambda _k\pi _i|\phi _i^k|$
eλktλkπi|ϕik|
units of the eigencurrent
$F^k_{ij}$
Fijk
per unit time.

Let us partition the set of states S into

$S^k_{+}$
S+k and
$S^k_{-}$
Sk
(see Eq. (6)). The corresponding cut-set (also known as cut) consists of all edges (i, j) where
$\phi _i^k\ge 0$
ϕik0
and
$\phi _j^k<0$
ϕjk<0
. We will call this cut the emission-absorption cut as it separates the states where the eigencurrent
$F^k_{ij}$
Fijk
is not absorbed from those where it is absorbed. It can be compared to the isocommittor cut21 corresponding to the committor value q = 0.5.

Now we will show that for every fixed time t, the eigencurrent over the emission-absorption cut is maximal among all possible cuts of the network, and it is equal to the total eigencurrent emitted by the states

$S^k_{+}$
S+k per unit time at time t, i.e.,

\begin{align}&\max _{S^\prime ,S^{\prime \prime }\ : S=S^\prime \dot{\cup }S^{\prime \prime }}\sum _{i\in S^\prime ,j\in S^{\prime \prime }}F^k_{ij}\nonumber \\&\quad =\sum _{i\in S^k_{+},j\in S^k_{-}}F^k_{ij}=e^{-\lambda _kt}\lambda _k\sum _{i\in S^k_{+}}\pi \phi ^k_i.\end{align}
maxS,S:S=ṠSiS,jSFijk=iS+k,jSkFijk=eλktλkiS+kπϕik.
(12)

(The symbol

$\dot{\cup }$
̇ denotes the disjoint union.) Since Eq. (9) implies that
$F^k_{ij}=-F_{ji}^k$
Fijk=Fjik
, for any subset SS we have

\begin{equation*}\sum _{i\in S^\prime ,j\in S^\prime }F^k_{ij}=0.\end{equation*}
iS,jSFijk=0.

Therefore, the eigencurrent over the cut (S, S) is

\begin{align*}\sum _{i\in S^\prime , j\in S^{\prime \prime }}F_{ij}^k &= \sum _{i\in S^\prime , j\in S}F_{ij}^k - \sum _{i\in S^\prime , j\in S^\prime }F_{ij}^k \nonumber \\&= e^{-\lambda _kt}\lambda _k\sum _{i\in S^\prime }\pi _i\phi ^k_i,\end{align*}
iS,jSFijk=iS,jSFijkiS,jSFijk=eλktλkiSπiϕik,

i.e., it is the total eigencurrent emitted in S per unit time. The maximum of the last sum is achieved if

$S^\prime =S^k_{+}$
S=S+k⁠, i.e., if S consists of all non-absorbing states.

The cut-set of the emission-absorption cut is the true transition state of the relaxation process from the perturbed distribution π + αPϕk.

If the considered network is relatively small and the pairwise transition rates are of the same order of magnitude, eigenvalues and eigenvectors of its generator matrix can be readily computed using standard methods.47,48 For large networks with complex underlying potential energy landscapes and pairwise transition rates varying by orders of magnitude, computing a collection of eigenpairs corresponding to small eigenvalues is a challenging problem. We approach this problem by taking advantage of the specific form of the pairwise transition rates (see Eq. (13) below) and employing theoretical results of Wentzel49,50 and Bovier and collaborators.12,13

In 1970s, Wentzel considered the spectral problem for stochastic networks without the detailed balance but with pairwise transition rates of the form

$e^{-U_{ij}/2\epsilon ^2}$
eUij/2ε2⁠, where ε was a small parameter49,50 (see also Ref. 51, Chap. 6). Such networks come, e.g., from continuous systems evolving according to non-gradient stochastic differential equations with small white noise.51 Wentzel obtained asymptotic estimates for eigenvalues up to the logarithmic order of the form
$\lambda _k\asymp e^{-\Delta _k/\epsilon }$
λkeΔk/ε
where the numbers Δk are the solutions of certain optimization problems over the so called W-graphs.49–51 

In 2000s, Bovier and collaborators developed sharp estimates for low lying eigenvalues and the corresponding eigenvectors for stochastic networks with detailed balance using the conceptual apparatus of the classic potential theory.12,13 (Low lying eigenvalues are those small eigenvalues that are separated from each other and the rest of the spectrum by notable spectral gaps.) They also extended their results for continuous systems.14,15 Their sharp estimates are obtained under the assumption of the existence of the hierarchy of sets of so called metastable points. The estimates for eigenvalues are expressed in terms of expected exit times from certain subsets of states (see below) while right eigenvectors are approximated by capacitors (also known as committors) for certain pairs of source and sink subsets of states.

The sharp estimates of Bovier et al.12,13 are of great theoretical value. However, in order to implement them in practice, one actually needs to construct the hierarchy of sets of metastable points. A recipe for doing this for networks representing energy landscapes was proposed in Ref. 20. The nonzero off-diagonal entries of the generator matrix L for such networks are of the form29 

\begin{equation}L_{ij} = \frac{a_{ij}}{a_i}e^{-(V_{ij}-V_i)/T},\end{equation}
Lij=aijaie(VijVi)/T,
(13)

where ai and aij are the prefactors for the local minimum i and the saddle (i, j) separating the local minima i and j, respectively, Vi and Vij are the values of the potential at i and (i, j), respectively, and T is the temperature. This recipe involves the construction of the set of the optimal W-graphs and also relies on estimates of Boviers et al. for eigenvectors12,13 and Wentzel's estimates for eigenvalues.49–51 

In the graph theory, a tree is a connected graph with no cycles; a forest is a collection of trees; a spanning tree for a given undirected graph is a tree with the same set of vertices; and a minimum spanning tree for a weighted undirected graph is a spanning tree such that the sum of the weights of its edges is minimal possible.52 A W-graph with k sinks is a directed forest such that in each of its k connected components (subtrees) all directed paths lead to a unique vertex called sink. The optimal W-graph with one sink is the minimum spanning tree where the cost of each edge (i, j) is Vij, and the directions of the edges are chosen so that all directed paths lead to the state with the smallest value of the potential. We will call it sink

$s^{\ast }_0$
s0*⁠.

We have proven20 that for the networks representing energy landscapes, the optimal W-graphs have a nested property which allows us to construct them recursively from top to bottom by solving a simple optimization problem at each step. At step k, we add one sink

$s^{\ast }_k$
sk* and remove one edge
$(p^{\ast }_k,q^{\ast }_k)$
(pk*,qk*)
. Then the asymptotics for the kth eigenvalue is
$\lambda _k\asymp e^{-\Delta _k/T}$
λkeΔk/T
, where
$\Delta _k:=V_{p^{\ast }_kq^{\ast }_k}-V_{s^{\ast }_k}$
Δk:=Vpk*qk*Vsk*
. The asymptotic eigenvector ϕk is the indicator function of the set of states of the connected component of the optimal W-graph at step k containing the sink
$s^{\ast }_k$
sk*
. The sinks fit the definition of metastable points12,13 if the temperature is low enough.

In this section, we propose a two-step approach for computing eigenvalues and eigenvectors. Step one is the algorithm introduced in Ref. 20 for finding zero-temperature asymptotics. Here, we will give its brief description. Step two is computing eigenpairs of interest at finite temperatures using the asymptotic eigenvectors as initial guesses.

For simplicity, we make the following genericness assumption: all numbers Vi, Vij and the differences VpqVi are different. The cases where this assumption is violated will be addressed in future work. Note that the genericness assumption implies that the minimum spanning tree is unique.

The necessary preprocessing for the algorithm20 is computation of the minimum spanning tree

$\mathcal {T}^{\ast }$
T* for the given network, where the cost of the edge (i, j) is Vij, the value of the potential at the saddle (i, j). In other words, one needs to extract a collection of N − 1 edges such that all N states are connected by these edges, and the sum of the costs of these edges is minimal possible. The rest of the edges are removed from the network. There are several greedy algorithms for performing this task.52 We have used Kruskal's algorithm.53 Note that all edges constituting the minimum spanning tree are such that for any pair of states i and j, the maximal cost Vpq along the unique path w*(i, j) in
$\mathcal {T}^{\ast }$
T*
connecting i and j is minimal possible.52 

Central to the algorithm20 are the barrier function u(i) and the escape function v(i), iS. For a given set of sinks

$s^{\ast }_j$
sj*⁠, j = 0, 1, …, k,

\begin{align}u(i) &= \min _{0\le j\le k}\max _{(p,q)\in w^{\ast }(i,s^{\ast }_j)}V_{pq},\end{align}
u(i)=min0jkmax(p,q)w*(i,sj*)Vpq,
(14)
\begin{align}v(i) &= u(i)-V_i,\end{align}
v(i)=u(i)Vi,
(15)

where

$w^{\ast }(i,s^{\ast }_j)$
w*(i,sj*) is the unique path in the minimum spanning tree
$\mathcal {T}^{\ast }$
T*
connecting states i and
$s^{\ast }_j$
sj*
. Hence, u(i) is the height of the lowest possible highest saddle separating state i from the set of sinks
$\lbrace s^{\ast }_j\rbrace _{j=0}^k$
{sj*}j=0k
, while v(i) is the minimal possible maximal potential barrier to overcome in order to escape from i to one of the sinks.

The algorithm20 proceeds as follows. The initial sink is

$s^{\ast }_0=\arg \min _{i\in S} V_i$
s0*=argminiSVi⁠. The functions u and v are computed with respect to this sink. The initial forest
$\mathcal {T}_0^{\ast }$
T0*
is the minimum spanning tree
$\mathcal {T}^{\ast }$
T*
. Set S0 = C0 = S. Then for k = 1, 2, …

  1. Pick the state with the largest value of the escape function v and make it the next sink

    $s^{\ast }_{k}$
    sk*⁠.

  2. Find the edge with the maximal value Vpq in the path connecting the new sink

    $s^{\ast }_{k}$
    sk* with one of the existing sinks
    $s^{\ast }_j$
    sj*
    , 0 ⩽ jk − 1, denote it by
    $(p^{\ast }_k,q^{\ast }_k)$
    (pk*,qk*)
    and remove it, i.e., set
    $\mathcal {T}^{\ast }_{k}=\mathcal {T}^{\ast }_{k-1}\backslash (p^{\ast }_k,q^{\ast }_k)$
    Tk*=Tk1*(pk*,qk*)
    .

  3. Set

    $\Delta _k:=V_{p^{\ast }_kq^{\ast }_k}-V_{s^{\ast }_{k}}$
    Δk:=Vpk*qk*Vsk*⁠. Set Sk to be the set of states in the connected component of the forest
    $\mathcal {T}^{\ast }_{k}$
    Tk*
    containing the new sink
    $s^{\ast }_{k}$
    sk*
    .

  4. Update the values of u and v in the connected component of the forest

    $\mathcal {T}^{\ast }_{k}$
    Tk* containing the new sink
    $s^{\ast }_{k}$
    sk*
    . Denote by Ck the set of states where the values of u and v have actually changed.

The first step of this algorithm for a chain-of-states example is illustrated in Fig. 2. An application of this algorithm to the network in Fig. 3 associated with the Lennard-Jones cluster of 7 atoms (LJ7)54,55 is worked out in Fig. 4.

FIG. 2.

A chain-of-states example. In this case, the minimum spanning tree is the graph itself. The sink

$s^{\ast }_0$
s0* is the state with the minimal value of the potential. The sink
$s_1^{\ast }$
s1*
is the state for which the escape function v computed with respect to
$s^{\ast }_0$
s0*
is maximal. The cutting edge
$(p^{\ast }_1,q^{\ast }_1)$
(p1*,q1*)
is the edge with the maximal potential value in the unique path connecting
$s^{\ast }_0$
s0*
and
$s^{\ast }_1$
s1*
. Remove the edge
$(p^{\ast }_1,q^{\ast }_1)$
(p1*,q1*)
. The quasi-invariant set S1 is the set of states in the connected component of the resulting forest containing the new sink
$s^{\ast }_1$
s1*
. The Freidlin's cycle C1 is the subset of S1 consisting of states separated from
$s_1^{\ast }$
s1*
by lower potential barriers than
$V_{p^{\ast }_1q^{\ast }_1}$
Vp1*q1*
.

FIG. 2.

A chain-of-states example. In this case, the minimum spanning tree is the graph itself. The sink

$s^{\ast }_0$
s0* is the state with the minimal value of the potential. The sink
$s_1^{\ast }$
s1*
is the state for which the escape function v computed with respect to
$s^{\ast }_0$
s0*
is maximal. The cutting edge
$(p^{\ast }_1,q^{\ast }_1)$
(p1*,q1*)
is the edge with the maximal potential value in the unique path connecting
$s^{\ast }_0$
s0*
and
$s^{\ast }_1$
s1*
. Remove the edge
$(p^{\ast }_1,q^{\ast }_1)$
(p1*,q1*)
. The quasi-invariant set S1 is the set of states in the connected component of the resulting forest containing the new sink
$s^{\ast }_1$
s1*
. The Freidlin's cycle C1 is the subset of S1 consisting of states separated from
$s_1^{\ast }$
s1*
by lower potential barriers than
$V_{p^{\ast }_1q^{\ast }_1}$
Vp1*q1*
.

Close modal
FIG. 4.

Initialization (k = 0) and three iterations (k = 1, 2, 3) of the algorithm for computing zero-temperature asymptotics for eigenvalues and eigenvectors applied to the stochastic network in Fig. 3 associated with LJ7.

FIG. 4.

Initialization (k = 0) and three iterations (k = 1, 2, 3) of the algorithm for computing zero-temperature asymptotics for eigenvalues and eigenvectors applied to the stochastic network in Fig. 3 associated with LJ7.

Close modal
FIG. 3.

A network associated with the energy landscape of the LJ7.54,55 The values of Vi and Vij, i, j ∈ {1, 2, 3, 4} are given with respect to the global potential minimum (the pentagonal bipyramid, V = −16.505). The edges constituting the minimum spanning tree are marked red.

FIG. 3.

A network associated with the energy landscape of the LJ7.54,55 The values of Vi and Vij, i, j ∈ {1, 2, 3, 4} are given with respect to the global potential minimum (the pentagonal bipyramid, V = −16.505). The edges constituting the minimum spanning tree are marked red.

Close modal

To summarize, for each k = 1, 2, …, N − 1 we get

  • the quasi-invariant set Sk;

  • the sink

    $s^{\ast }_k$
    sk*⁠, which is the deepest minimum in the quasi-invariant set Sk;

  • the cutting edge

    $(p^{\ast }_k,q^{\ast }_k)$
    (pk*,qk*)⁠, which corresponds to the lowest possible highest saddle separating the sink
    $s^{\ast }_k$
    sk*
    from the previously chosen sinks
    $s^{\ast }_j$
    sj*
    , j = 0, 1, …, k − 1;

  • the number

    $\Delta _k=V_{p^{\ast }_kq^{\ast }_k}-V_{s^{\ast }_{k}}$
    Δk=Vpk*qk*Vsk*⁠;

  • Freidlin's cycle33,51,56,57CkSk which is the set of states separated from the sink

    $s^{\ast }_k$
    sk* by smaller potential barriers than
    $V_{p^{\ast }_kq^{\ast }_k}$
    Vpk*qk*
    . It is the largest Freidlin's cycle20 containing the new sink and not containing any state with a smaller potential value. Its significance will be explained in Sec. VII.

The asymptotics for the eigenvalues and the right eigenvectors are given by

\begin{align}\lambda _k&\asymp e^{-\Delta _k/T},\end{align}
λkeΔk/T,
(16)
\begin{align}\phi ^k_i&\approx \left\lbrace \begin{array}{@{}l@{\quad }l@{}}1,&i\in S_k\\[4pt]0,& i\notin S_k.\end{array}\right.\end{align}
ϕik1,iSk0,iSk.
(17)

Note the ordering: Δ1 < Δ2 < …. Hence, the eigenpairs are computed in the increasing order of the corresponding eigenvalues. Therefore, one can stop the for-cycle whenever all eigenpairs of interest are computed.

Finally, we remark on what happens if the genericness assumption is violated. If Δk = Δk+1 for some k then two cases are possible. If the corresponding cutting edges belong to different connected components of the forest

$\mathcal {T}^{\ast }_k$
Tk*⁠, then no error occurs. If the cutting edges belong to the same connected component, then the asymptotic eigenvectors ϕk and ϕk+1 might be determined wrongly, but the error does not propagate for larger k in the for-cycle. In other words, the possible effect from the violation of the genericness assumption is local.

Difficulties in computing eigenpairs at finite (but still low) temperatures are associated with the facts that (i) all eigenvalues tend to zero as T → 0, and the ones of interest can be extremely small at the desired range of temperatures; (ii) the eigenvalues can cross as the temperature increases; (ii) the right eigenvectors can be nearly parallel. Experimenting with a few continuation techniques we have found out that the Rayleigh Quotient iteration (see, e.g., Refs. 47,48) with the asymptotic estimates (17) for the eigenvectors as initial guesses is a feasible and robust approach as soon as some precautions are taken.

The initial approximations given by Eq. (17) are sufficiently accurate so that the Rayleigh Quotient iteration converges rapidly (cubically). However, occasionally, convergence to a wrong eigenpair may occur. For example, this can happen for eigenpairs number k and l such that (i) the sets Sk and Sl are both large, (ii) SkSl, (iii) the number of states in Sl\Sk is small in comparison with the number of states in Sk, and (iv) the potential values at sinks

$s_k^{\ast }$
sk* and
$s^{\ast }_l$
sl*
are close. Fortunately, the convergence to a wrong eigenpair is easy to detect using the orthogonality relationship (4).

Now we provide some technical details of computing finite temperature eigenpairs. Instead of L, we work with the matrix Lsym given by Eq. (1). Its eigenvectors ψk are orthonormal and relate to the right eigenvectors ϕk of L via ψk = P1/2ϕk. To compute the kth eigenpair of Lsym at temperature T, we set the initial guess for the eigenvector ψk as

\begin{equation*}\left(\psi _i^k\right)_0(T)=\left\lbrace \begin{array}{@{}l@{\quad }l@{}}\pi ^{1/2}_i(T),&i\in S_k,\\[4pt]0,&i\notin S_k,\end{array}\right.\end{equation*}
ψik0(T)=πi1/2(T),iSk,0,iSk,

run the Rayleigh Quotient iteration and obtain an eigenpair (λ, ξ). Then we check whether (λ, ξ) is a correct eigenpair by calculating the dot product ξTk)0. Since eigenvectors of Lsym corresponding to distinct eigenvalues are orthogonal, this inner product is close to zero if (λ, ξ) is a wrong eigenpair. Otherwise, it is close to one. If (λ, ξ) is the desired eigenpair, we set λk(T) = λ and ϕk(T) = P−1/2(T)ξ. If (λ, ξ) is a wrong eigenpair, we repeat the Rayleigh Quotient iteration with a more accurate initial guess for the eigenvector (and/or eigenvalue) obtained from already computed eigenpairs for the same k but different values of temperature.

The LJ38 dataset created by Wales's group26 contains 100 000 local minima and 138 888 transition states. Its largest connected component consists of 71 887 local minima (states) and 119 853 transition states (edges) and contains the two deepest minima FCC and ICO. We will analyze this connected component and refer to it as the LJ38 network. Local minima, other than FCC and ICO, will be referred to by their indices in the list.26 Here are the indices of some important minima: FCC has index 1, ICO has index 7, the third lowest minimum has index 16, the lowest possible highest saddle separating ICO and FCC is the transition state between minima 342 and 354.

Following Ref. 29 we define pairwise transition rates Lij using the harmonic approximation. If there is a transition state in the database26 separating local minima i and j, then

\begin{equation}L_{ij}(T)=\frac{O_i\bar{\nu }_i^{\kappa }}{O_{ij}\bar{\nu }_i^{\kappa -1}}e^{-(V_{ij}-V_i)/T},\end{equation}
Lij(T)=Oiν¯iκOijν¯iκ1e(VijVi)/T,
(18)

where Oi, Vi, and

$\bar{\nu }_i$
ν¯i are, respectively, the point group order, the value of the potential energy, and the geometric mean vibrational frequency for minimum i; Oij, Vij, and
$\bar{\nu }_{ij}$
ν¯ij
are the same parameters for the saddle (i, j); and κ = 3 × 38 − 6 = 108 is the number of vibrational degrees of freedom. Otherwise, if there is no such a transition state (i.e., there is no edge (i, j)), Lij = 0. The equilibrium probability distribution is given by

\begin{equation}\pi _i(T)=\frac{1}{Z(T)}\frac{e^{-V_i/T}}{O_i\bar{\nu }_i^{\kappa }},\quad Z(T)=\sum _{i\in S}\frac{e^{-V_i/T}}{O_i\bar{\nu }_i^{\kappa }}.\end{equation}
πi(T)=1Z(T)eVi/TOiν¯iκ,Z(T)=iSeVi/TOiν¯iκ.
(19)

To choose the range of temperatures at which we will investigate the LJ38 network, we use its critical temperature of the solid-solid phase transition as a reference. It has been established2,28,35 that the solid-solid phase transition, when fcc structures give place to icosahedral packings, occurs in the LJ38 cluster at T = 0.12. The outer layer starts to melt, i.e., icosahedral states become less populated than liquid-like states, at T = 0.18. These critical temperatures were obtained in the continuous setting using either Monte Carlo simulations35 or the superposition method for determination of the density of states and an anharmonic approximation for the pairwise transition rates.28,58 However, if the continuous LJ38 cluster is modeled by the described above network with the pairwise rates given by Eq. (18), the critical temperatures are shifted upward. One can determine them, e.g., by plotting the formally computed2,28 heat capacity at constant volume: the solid-solid transition occurs at T = 0.16, while the population of the liquid-like states passes the population of the icosahedral ones at T = 0.27. The critical temperature 0.16 for this network model with harmonic pairwise rates (Eq. (18)) was also obtained in Ref. 21 as the point where the transition rates FCC → ICO and ICO → FCC cross. This discrepancy in critical temperatures illustrates and quantifies the difference between the continuous LJ38 cluster and its network approximation.

Originally, the results involving the zero-temperature asymptotics for the eigenvalues and eigenvectors of the LJ38 network were presented in Ref. 20. Here, we highlight the main one of them. The results associated with the finite temperature continuation are new.

1. Is there a spectral gap?

The numbers Δk determining the zero-temperature asymptotics for eigenvalues of the generator matrix L of the LJ38 network by

$\lambda _k\asymp e^{-\Delta _k/T}$
λkeΔk/T have been originally presented in Ref. 20. Here, we reproduce them for reader's convenience (Fig. 5). It is surprising at the first glance that the numbers Δk are arranged pretty densely. A few largest Δk's separated by visible gaps correspond to high-lying isolated liquid-like local minima 76 525 (k = 1), 91 143 (k = 2), 16 497 (k = 3), etc., of no particular interest. ICO, the second lowest minimum, corresponds to k = 245. The numbers Δk around k = 245 are zoomed in a separate window in Fig. 5. The difference between Δ245 and Δ246 is only 0.0036. Therefore, the corresponding eigenvalues are separated by an appreciable gap only for extremely low temperatures, at least T < 0.0036.

FIG. 5.

The numbers Δk, k = 1, 2, …, 71 886, defining the zero-temperature asymptotics for the eigenvalues of the LJ38 network by λk ≈ − Δk/T. The white spot on the graph is blown up in the separate window. The escape process from the quasi-invariant subset of states where the deepest minimum is the second lowest minimum ICO corresponds to k = 245. Reproduced with permission from M. K. Cameron, “Computing the asymptotic spectrum for networks representing energy landscapes using the minimum spanning tree,” Networks Heterogeneous Media 9, 383 (2014). Copyright 2014 American Institute of Mathematical Sciences.

FIG. 5.

The numbers Δk, k = 1, 2, …, 71 886, defining the zero-temperature asymptotics for the eigenvalues of the LJ38 network by λk ≈ − Δk/T. The white spot on the graph is blown up in the separate window. The escape process from the quasi-invariant subset of states where the deepest minimum is the second lowest minimum ICO corresponds to k = 245. Reproduced with permission from M. K. Cameron, “Computing the asymptotic spectrum for networks representing energy landscapes using the minimum spanning tree,” Networks Heterogeneous Media 9, 383 (2014). Copyright 2014 American Institute of Mathematical Sciences.

Close modal

The reason for absence of significant spectral gaps lies in the fact that many liquid-like states in the LJ38 network are (i) separated from the icosahedral and/or the fcc funnels by extremely high saddles and (ii) isolated, i.e., correspond to small sets Sk and Ck. The majority of those sets Sk contain less than 5 states. The values of the potential at the saddles

$V_{p^{\ast }q^{\ast }}$
Vp*q* and the local minima
$V_{s^{\ast }_k}$
Vsk*
corresponding to the first 1500 sinks are shown in Fig. 6.

FIG. 6.

The values of the potential at the saddles

$V_{p^{\ast }q^{\ast }}$
Vp*q* and the local minima
$V_{s^{\ast }_k}$
Vsk*
corresponding to the first 1500 sinks. Among them, only two sinks correspond to icosahedral minima: ICO (k = 245) and minimum 264 (k = 1379).

FIG. 6.

The values of the potential at the saddles

$V_{p^{\ast }q^{\ast }}$
Vp*q* and the local minima
$V_{s^{\ast }_k}$
Vsk*
corresponding to the first 1500 sinks. Among them, only two sinks correspond to icosahedral minima: ICO (k = 245) and minimum 264 (k = 1379).

Close modal

Actually, if, following Ref. 28, we call a state liquid-like if its energy exceeds −171.64 (or 2.29 with respect to FCC), the first two states with the smallest positive sink index that are not liquid-like are ICO (k = 245) and state 264 (k = 1379).

Roughly speaking, those liquid-like states, separated by extremely high saddles, pollute the interesting part of the spectrum, which corresponds to the states that are dynamically accessible from icosahedral or fcc states within reasonable waiting times. This fact motivates us to restrict the observation time.36 This is equivalent to the introduction of a cap Vmax  on the potential energy. Then only the subnetwork G(S, E), where the sets of states and edges are given by

\begin{align*}S^\prime &=\lbrace i\in S \,\vert\, \max _{(p,q)\in w^{\ast }(i,{\rm FCC)}}V_{pq}< V_{\max }\rbrace ,\\E^\prime &=\lbrace (i,j)\in E \,\vert\, i\in S^\prime , j\in S^\prime , V_{ij}< V_{\max }\rbrace ,\end{align*}
S={iS|max(p,q)w*(i, FCC )Vpq<Vmax},E={(i,j)E|iS,jS,Vij<Vmax},

is considered. For Vmax  = 6 with respect to FCC (the lowest possible highest barrier separating ICO and FCC is27 4.219), the resulting network contains 30 520 states and 71 750 edges. In this truncated network, Δ1, corresponding to the sink

$s^{\ast }_1 =$
s1*= ICO, is separated by a gap of approximately 0.19 from Δ2 corresponding to the sink
$s^{\ast }_2 = 223$
s2*=223
(minimum 223) and the cutting edge with potential energy
$V_{p^{\ast }_2q^{\ast }_2}=5.840$
Vp2*q2*=5.840
. If we pick Vmax  = 5.5, the gap between Δ1 and Δ2 will increase to 1.00, while Vmax  = 5 will make this gap 1.05. More details can be found in Ref. 20.

2. Finite temperature eigenvalues, eigenvectors, and eigencurrents

The eigenvalues corresponding to the eigenpairs with the largest sets Sk are plotted in Fig. 7. These eigenpairs are identified by20k = 245 (sink ICO (7), |S(ICO)| = 56 290), k = 6910 (sink 3, |S(3)| = 4252), k = 7482 (sink 4, |S(4)| = 1316), k = 4143 (sink 5215, |S(5215)| = 990), k = 11 750 (sink 3551, |S(3551)| = 680), and k = 11 961 (sink 2052, |S(2051)| = 1758). The linear least squares fits for these eigenvalues and the numbers Δk are

\begin{equation*}\begin{array}{lll}\rm ICO: & \Delta = 3.543,&\lambda \approx 1.417\times 10^5\cdot e^{-3.570/T},\\3:& \Delta = 1.408,&\lambda \approx 4.811\times 10^4\cdot e^{-1.485/T},\\4:& \Delta = 1.356,&\lambda \approx 8.222\times 10^3\cdot e^{-1.360/T},\\5215:& \Delta = 1.753,&\lambda \approx 4.253\times 10^3\cdot e^{-1.706/T},\\3551:& \Delta = 1.039,&\lambda \approx 1.122\times 10^4\cdot e^{-1.037/T},\\2052:& \Delta = 1.027,&\lambda \approx 2.029\times 10^3\cdot e^{-1.007/T}. \end{array}\end{equation*}
ICO :Δ=3.543,λ1.417×105·e3.570/T,3:Δ=1.408,λ4.811×104·e1.485/T,4:Δ=1.356,λ8.222×103·e1.360/T,5215:Δ=1.753,λ4.253×103·e1.706/T,3551:Δ=1.039,λ1.122×104·e1.037/T,2052:Δ=1.027,λ2.029×103·e1.007/T.
FIG. 7.

Eigenvalues corresponding to a collection of large sets Sk are plotted versus 1/T where T is the temperature. The sinks for these sets are ICO (minimum 7) and minima 3, 4, 5215, 3551, and 2052. The lines are obtained by the least squares fit for the low temperature data with linear functions.

FIG. 7.

Eigenvalues corresponding to a collection of large sets Sk are plotted versus 1/T where T is the temperature. The sinks for these sets are ICO (minimum 7) and minima 3, 4, 5215, 3551, and 2052. The lines are obtained by the least squares fit for the low temperature data with linear functions.

Close modal

The exponents of the least squares are close to their asymptotic values. Note that the eigenvalues corresponding to sinks 3 and 4 cross at T = 0.07. Also note that the eigenvalue corresponding to sink ICO grows faster than predicted by the Arrhenius law. This is consistent with the graph of the transition rate from icosahedral to fcc states in Ref. 2. This phenomenon is due to the fact that weights of some states other than ICO (e.g., 16 and 8) in the left eigenvector P(T(ICO) grow dramatically as the temperature increases. As a result, the effective potential at the starting position effectively increases. On the other hand, there is an abundance of pathways in the LJ38 network connecting icosahedral and fcc states that pass through higher than the lowest possible one but not very high saddles.21 This, in turn, leads to a dramatic broadening of the eigencurrent distribution in the emission-absorption cut (Fig. 8). Nearly a power law distribution is observed. Fig. 8 can be compared to Fig. 9 in Ref. 21 where the reaction flux distribution in the isocommittor cut corresponding to the committor value q = 0.5 is shown.

FIG. 8.

The distribution of the eigencurrent in the emission-absorption cut at different temperatures is shown. The empirical cumulative distribution function is shown in the inset.

FIG. 8.

The distribution of the eigencurrent in the emission-absorption cut at different temperatures is shown. The empirical cumulative distribution function is shown in the inset.

Close modal

The asymptotic eigenvector ϕ(ICO)(0) corresponding to the sink ICO and its finite temperature counterparts ϕ(ICO)(0.10) and ϕ(ICO)(0.18) at T = 0.10 and T = 0.18, respectively, are shown in Fig. 9 in the form of disconnectivity graphs.61,62 There is a close agreement between ϕ(ICO)(0) and ϕ(ICO)(0.10), while the components of ϕ(ICO)(0.18) corresponding to the fcc states grow in the negative direction. This reflects the fact that icosahedral states are more populated than fcc states in the equilibrium distribution at T = 0.18. The majority of the components of ϕ(ICO) remain in the shown interval [−2, 2]. At the same time, a few components corresponding to high-lying states grow large in the absolute value. We could not show them without making the colormap non-informative. The ranges of values of ϕ(ICO)(0.10) and ϕ(ICO)(0.18) are [− 9408.0, 577.1] and [− 116.8, 874.0], respectively. At T = 0.10, 104 out of 71 887 components of ϕ(ICO) are greater than 2, and 60 components are less than −2, while at T = 0.18, there are 498 components greater than 2, and 323 components less than −2. The emission-absorption cut changes as the temperature grows. In particular, one can spot a state in Fig. 9 separated by a high barrier from both ICO and FCC, that switches sign as the temperature increases from 0.10 to 0.18.

FIG. 9.

The asymptotic eigenvector (left) and the eigenvectors at T = 0.10 (middle) and at T = 0.18 (right) corresponding to the sink ICO are visualized in the form of disconnectivity graphs. The shown states are selected as follows. First, 100 states with the smallest potential energy are selected. Then those that are separated by barriers less than 0.2 are lumped together. The states are colored according to

$\phi ^{({\rm ICO})}_j(T)$
ϕj( ICO )(T)⁠, T = 0, 0.10, and 0.18, respectively. The states are arranged along the x-axis in the decreasing order according to
$\phi ^{({\rm ICO})}_j(T=0.18)$
ϕj( ICO )(T=0.18)
.

FIG. 9.

The asymptotic eigenvector (left) and the eigenvectors at T = 0.10 (middle) and at T = 0.18 (right) corresponding to the sink ICO are visualized in the form of disconnectivity graphs. The shown states are selected as follows. First, 100 states with the smallest potential energy are selected. Then those that are separated by barriers less than 0.2 are lumped together. The states are colored according to

$\phi ^{({\rm ICO})}_j(T)$
ϕj( ICO )(T)⁠, T = 0, 0.10, and 0.18, respectively. The states are arranged along the x-axis in the decreasing order according to
$\phi ^{({\rm ICO})}_j(T=0.18)$
ϕj( ICO )(T=0.18)
.

Close modal

The relaxation process corresponding to the kth eigenpair is quantified by the kth eigencurrent defined by Eq. (9). The eigencurrents for the sink ICO at temperatures T = 0.10 and T = 0.18 are depicted in Figs. 10(a) and 10(b), respectively. The arrows representing the maximal eigencurrent, i.e.,

$F_{\max }^{({\rm ICO})}(T):=\max _{ij}F^{({\rm ICO})}_{ij}(T)$
Fmax( ICO )(T):=maxijFij( ICO )(T)⁠, have the unit (the maximal) width. Only the edges (i, j) with
$F_{ij}^{({\rm ICO})}\ge 0.05\cdot F_{\max }^{({\rm ICO})}(T)$
Fij( ICO )0.05·Fmax( ICO )(T)
are shown, and their widths are proportional to
$F_{ij}^{({\rm ICO})}$
Fij( ICO )
.

FIG. 10.

The eigencurrents corresponding to sink ICO at temperatures T = 0.10 (a), and T = 0.18 (b). The widths of the arrows are proportional to

$F^{({\rm ICO})}_{ij}$
Fij( ICO )⁠. The states which emit or absorb at least 5% of the eigencurrent are shown by maroon or bark blue solid circles, respectively, and the areas of these circles are proportional to the percentages of emitted/absorbed eigencurrent. Other states are represented by empty circles. The dashed lines indicate the emission-absorption cuts.

FIG. 10.

The eigencurrents corresponding to sink ICO at temperatures T = 0.10 (a), and T = 0.18 (b). The widths of the arrows are proportional to

$F^{({\rm ICO})}_{ij}$
Fij( ICO )⁠. The states which emit or absorb at least 5% of the eigencurrent are shown by maroon or bark blue solid circles, respectively, and the areas of these circles are proportional to the percentages of emitted/absorbed eigencurrent. Other states are represented by empty circles. The dashed lines indicate the emission-absorption cuts.

Close modal

In Fig. 10(a) for T = 0.10, the sequence of states

\begin{equation*}354\rightarrow 3590\rightarrow 3185\rightarrow 958\end{equation*}
35435903185958

stands out by thicker arrows connecting them. This sequence is a part of the asymptotic zero-temperature path (also known as the MinMax path) connecting ICO and FCC:33 {ICO, 9, 8, 248, 284, 312, 313, 371, 372, 374, 375, 376, 17895, 3213, 350, 351, 352, 342, 354, 3590, 3185, 958, 2831, 607, 3, 4, FCC}. This path can be traced almost fully in Fig. 10(a). The two highest saddles in the MinMax path are V342, 354 = 4.219 and V958, 2831 = 4.186. The height of the saddle between minima 958 and 1, corresponding to another thick arrow across the emission-absorption cut, is26 V958, 1 = 4.272.

At T = 0.18, the eigencurrent is considerably more spread out, especially to the right from the emission-absorption cut (see Fig. 10(b)). The number of pathways carrying notable amount of the eigencurrent from the emission-absorption cut to FCC dramatically increases, and this cut itself shifts away from FCC. The fraction of the eigencurrent carried by the sequence 354 → 3590 → 3185 → 958 notably drops in comparison with it at T = 0.10.

The collection of saddles corresponding to the edges of the emission-absorption cut is the natural transition state of the relaxation process. As shown in Fig. 8, the distribution of current in this cut is wide, nearly a heavy tailed power law. This means that the relaxation process dramatically diversifies as the temperature increases.

Figs. 10(a) and 10(b) also show the distributions of the emission and the absorption of the eigencurrent. At the range of temperatures 0.09 ⩽ T ⩽ 0.18, minimum 16 carries the maximal weight in the perturbation given by the left eigenvector Pϕ(ICO). The weight of minimum 8 considerably increases. At the same time, the weight of ICO, the deepest minimum among the emitting states, drops due to its small prefactor (which, in turn, is due to its high order of the point group (OICO = 10)). As the temperature increases, more and more states emit appreciable amounts of the eigencurrent. On the other hand, at this temperature range, FCC remains the dominant absorbing state: it absorbs over 99% of the emitted eigencurrent even at T = 0.18.

In Ref. 21, the LJ38 network was analyzed using tools of the TPT.22–25 Throughout this paper, we have indicated some points of comparison between the proposed spectral approach and TPT. Now we will summarize similarities and differences between these two approaches.

The basic steps of the analysis of a transition process taking place in the stochastic network G(S, E, L) by means of TPT are the following. We assume detailed balance for simplicity.

  1. Pick two disjoint subsets of states: a source set AS and a sink set BS.

  2. Solve the committor equation
    \begin{equation}\sum _{j\in S\backslash (A\cup B)}L_{ij}q_j=0,\quad q(A)=0,\quad q(B)=1.\end{equation}
    jS(AB)Lijqj=0,q(A)=0,q(B)=1.
    (20)
    For each state i, the value of the committor qi is the probability that the random walk starting at i will reach first B rather than A.
  3. Calculate the reactive current
    \begin{equation*}F^{R}_{ij}:=\pi _iL_{ij}(q_j-q_i)\end{equation*}
    FijR:=πiLij(qjqi)
    and analyze it. For example, one can consider its distribution in the isocommittor cuts and visualize the reactive tube with their aid. (For a given q ∈ [0, 1) the isocommittor cut21C(q) is the collection of edges (i, j) such that qiq and qj > q.)
  4. Calculate the transition rate νR, i.e., the average number of transitions from A to B per unit time, e.g., by
    \begin{equation*}\nu _R=\sum _{(i,j)\in C(q)}F^{R}_{ij},\end{equation*}
    νR=(i,j)C(q)FijR,
    where C(q) is an arbitrary isocommittor cut. Also calculate the rates kA, B and kB, A which are the inverse average times of going back to B after hitting A and going back to A after hitting B, respectively. They are given by
    \begin{equation*}k_{A,B}=\frac{\nu _R}{\sum _{i\in S}\pi _i(1-q_i)},\quad k_{A,B}=\frac{\nu _R}{\sum _{i\in S}\pi _iq_i}.\end{equation*}
    kA,B=νRiSπi(1qi),kA,B=νRiSπiqi.

In TPT, the probability distribution in the network is assumed to be equilibrium, i.e., π. The transition process from A to B is stationary. The committor, the reactive current, and the transition rates are time-independent. Contrary to this, the eigenpair (λk, ϕk) describes the time-dependent relaxation process starting from the non-equilibrium distribution π + αPϕk, where the perturbation αPϕk decays uniformly throughout the network with the rate λk.

Nevertheless, it is instructive to compare the right eigenvector ϕk to the committor, the eigencurrent

\begin{equation*}F^k_{ij}=\pi _iL_{ij} \big(\phi ^k_j-\phi ^k_i\big)e^{-\lambda _k/T}\end{equation*}
Fijk=πiLijϕjkϕikeλk/T

to the reactive current

$F^{R}_{ij}$
FijR⁠, and the transition rate kA, B to λk. We will do it on the example of the LJ38 network.

As the temperature tends to zero, the right eigenvector ϕk tends12,13 to the committor qAB with the source set

$A=\lbrace s^{\ast }_0,s^{\ast }_1,\ldots ,s^{\ast }_{k-1}\rbrace$
A={s0*,s1*,...,sk1*} and the sink set
$B=\lbrace s^{\ast }_k\rbrace$
B={sk*}
, which, in turn, tends to the indicator function of the set Sk. In Ref. 21, the source and sink sets were chosen to be A ≡ ICO and B ≡ FCC. Therefore, if we narrow our consideration to the subset of states
$S\backslash \bigcup _{j=1}^{244}S_j$
Sj=1244Sj
(recall that ICO is the sink
$s^{\ast }_{245}$
s245*
), we can expect that the right eigenvector ϕ(ICO) is close to 1 − q(ICO, FCC) at low temperatures. Fig. 9 visualizing the eigenvector ϕ(ICO) at T → 0, T = 0.10, and T = 0.18 can be compared to Fig. 4 in Ref. 21 depicting the committor at T = 0.06, 0.12, and 0.18.

The committor and the right eigenvector play similar roles in the definitions of the corresponding currents. The net average numbers of transitions per unit time along the edge (i, j) in the Transition Path Process21,25 in TPT and the relaxation process are proportional to πiLij(qjqi) and

$\pi _iL_{ij}(\phi ^k_j-\phi ^k_i)$
πiLij(ϕjkϕik)⁠, respectively. However, there are important differences between the right eigenvector and the committor. While the committor indicates the probability to reach first B rather than A starting from the given state, the right eigenvector indicates how the given state is over- or underpopulated relative to the equilibrium distribution. The committor takes values in the interval [0, 1], while the range of values of the components of the right eigenvector ϕ(ICO), normalized as described in Sec. II A, is temperature dependent. Some components acquire values much larger than 1 or much less than −1.

Both the reactive current

$F^R_{ij}$
FijR and the eigencurrent
$F^{({\rm ICO})}_{ij}$
Fij( ICO )
describe the net probability flow for the corresponding processes. The reactive current is time-independent, while the eigencurrent uniformly decays with time at the rate λ(ICO). The reactive current is emitted by the source states A, absorbed by the sink states B, and conserved at all other states.24 There is no reactive current between any two source states, and there is no reactive current between any two sink states. Contrary to this, the eigencurrent is emitted at all states where ϕk > 0, absorbed at all states where ϕk < 0, and there is a nonzero eigencurrent along any edge (i, j) as long as
$\phi ^k_i\ne \phi ^k_j$
ϕikϕjk
. The total flux of the reactive current is the same through any cut of the network separating the sets A and B. The total flux of the eigencurrent is maximal through the emission-absorption cut, which is the cut separating the states with
$\phi ^k_j\ge 0$
ϕjk0
and
$\phi ^k_j<0$
ϕjk<0
. Note if
$\phi ^k_j=0$
ϕjk=0
at some state j, then the emission-absorption cut is not the only cut separating the states with ϕk > 0 and ϕk < 0. In this case, the flux of the eigencurrent is maximal through any cut separating the states with ϕk > 0 and ϕk < 0.

Despite the differences, the reactive current can serve as a good approximation to the quantity

$\pi _iL_{ij}(\phi ^k_j-\phi ^k_i)$
πiLij(ϕjkϕik) if almost all eigencurrent is emitted and absorbed by a few states. For example, this is the case for the eigencurrent
$F^{({\rm ICO})}_{ij}$
Fij( ICO )
in the LJ38 network at low temperatures, say T ⩽ 0.12. The distributions of the eigencurrent
$F^{({\rm ICO})}_{ij}$
Fij( ICO )
in the emission-absorption cut and the reactive current
$F^R_{ij}$
FijR
in the isocommittor cut corresponding to q = 0.5 are similar (compare Fig. 8 and 9 in Ref. 21).

At low temperatures, both the transition rate kICO, FCC from ICO to FCC and the escape rate λ(ICO) from the set S(ICO) are well-approximated by the Arrhenius law with the same parameters. The least squares fits to the Arrhenius law obtained here and in Ref. 21 are close,

\begin{align*}k_{{\rm ICO,FCC}}&=9.81\times 10^4\cdot e^{-3.525/T},\\\lambda ({\rm ICO})&=1.417\times 10^5\cdot e^{-3.570/T}.\end{align*}
k ICO , FCC =9.81×104·e3.525/T,λ( ICO )=1.417×105·e3.570/T.

The theoretical asymptotic value of the exponent is Δ = 3.543. As the temperature increases, the graphs of kICO, FCC and λ(ICO) diverge (compare Fig. 7 with Fig. 5 in Ref. 21). λ(ICO) grows faster than kICO, FCC due to the fact that the emission of the eigencurrent is distributed. The fraction emitted by ICO drops in favor of higher lying states as the temperature increases.

Finally, we would like to make a few remarks about technical difficulties in applying the TPT approach and the proposed spectral approach. The major computational challenge in the TPT lies in solving the committor equation (20) which is a large linear system. However, for networks representing energy landscapes, the matrix of the system can be made symmetric. In addition, it is positive-definite by construction. Hence, the powerful Conjugate Gradient (CG) method with the incomplete Cholesky preconditioning47 can be used to solve it. The main technical difficulty in the proposed spectral approach is solving the large linear system in the Rayleigh Quotient iteration. The system is symmetric but indefinite, and its numerical solution is harder than the one of the committor equation. Note that the computation of the zero-temperature asymptotics in the spectral approach is robust as soon as the genericness assumption mostly holds. Its violations might lead only to local effects.

A state of an isolated physical system is called metastable, if it is not its ground state, but the system can remain in it for a long time. This concept is clear when the system has just two states, but even then some reference time is needed in order to determine what period of time is long. For example, it can be the observation time. It is a nontrivial task to extend the concept of metastability to complex systems. In a number of works,12,13,16,17 the concept of metastability was associated with the presence of a spectral gap separating a group of small eigenvalues from the rest of the spectrum.

In this section, we will discuss whether LJ38 is metastable from two points of view: of a mathematician and of a chemical physicist.

A detailed discussion on the metastability of the LJ38 network from the point of view of a mathematician can be found in Ref. 20. Here, we highlight its conclusions. A formal definition of metastability for the case of Markov chains with detailed balance was introduced by Bovier, Eckhoff, Gayrard, and Klein.12,13 They have defined the set of metastable points

$\mathcal {M}\subset S$
MS⁠, each of which, in essence, is a representative of a quasi-invariant subset of states. Then the stochastic network is metastable with respect to the set
$\mathcal {M}$
M
, if the minimal expected time to reach any metastable point from another metastable point is much larger than the maximal expected time to reach any metastable point from any state not belonging to
$\mathcal {M}$
M
. The absence of the spectral gap in the full LJ38 network renders it not metastable in the sense of Bovier et al.12,13 with respect to a set of metastable points one of which is ICO, unless the temperature is extremely small, i.e., T < 0.0036. However, capping the potential we can make ICO a metastable point for the considered range of temperatures 0 < T ⩽ 0.18. Huisinga, Meyn, and Schuette proposed another definition of metastability in the context of continuous systems where the time reversibility (the detailed balance) was not assumed.18,19 Their definition links metastability with ergodicity. In their sense, the Freidlin's cycle C(ICO) is metastable for the whole considered range of temperatures 0 < T ⩽ 0.18.

Now we assume the point of view of a chemical physicist. The Lennard-Jones-38 network is a model for the low temperature dynamics of rare gas clusters of 38 atoms such as argon, krypton, and xenon. Rare gas clusters have been studied via electron diffraction37–43 and extended x-ray absorption fine structure spectroscopy44 experiments. Mass spectra revealing increased stability of clusters of “magic” numbers of atoms corresponding to completion (13, 55, 147, 309, …) or partial completion (19, 25, …) of icosahedral layers were obtained.38–41 The number 38 is not in these lists. The comparison of the experimentally determined mean radial distribution function with the simulated one allowed Kakar et al.44 to infer the structure of clusters for different sizes. Smaller clusters tend to be icosahedral, while larger ones exhibit signs of the fcc packing. The structural transition from icosahedral packing to fcc occurs around cluster sizes from44N = 200 to37,40N = 800 atoms. It has been hypothesized by van der Waal59 and confirmed experimentally by Kovalenko et al.42,43 that the structural transition from icosahedral packing to fcc occurs via the build up of faulty face-centered cubic layers around the icosahedral core rather than via cluster rearrangement. Simulated diffraction patterns for faulty fcc clusters were found to be in a better agreement with the experimental ones than the patterns simulated from other possible packings.

These experimental results suggest that the structure of the clusters of 38 atoms generated in experiments is icosahedral. Otherwise, a peak at 38 atoms in the mass spectra would be present. Before the structural transition to the global minimum (FCC) has chance to occur, these clusters are likely to acquire more atoms. Therefore, the Lennard-Jones-38 clusters are more likely to be found in the icosahedral state that is metastable at low temperatures.

In order to make our results comparable with experimental ones, we divide the collection of states into icosahedral and fcc. This partition roughly corresponds to the division into emitting and absorbing states for the eigencurrent FICO. Using parameters appropriate for argon2 we obtain the escape rates listed in Table I. T, * and r, * denote the temperature and the escape rate in the reduced units used throughout the text. T, K, is the temperature in Kelvins. r, s−1, is the escape rate in inverse seconds. The calculated rates are rather small but comparable with the time scales of argon clusters dynamics.60 Therefore, if there would be set up an experiment where argon clusters of 38 atoms would be selected and isolated, structural rearrangements from icosahedral packing to fcc could be observed. Then the measured rates of such structural rearrangements could be compared with the calculated rates in Table I.

Table I.

Escape rates from icosahedral states to fcc in the LJ38 network for parameters appropriate for argon. The symbol “*” denotes the reduced units used throughout the text.

T, *T (K)r, *r (s−1)
0.06 7.2 2.0 × 10−21 9.5 × 10−10 
0.08 9.6 5.9 × 10−15 2.7 × 10−3 
0.10 12.0 4.3 × 10−11 2.0 × 101 
0.12 14.4 1.8 × 10−8 8.1 × 103 
0.14 16.8 1.6 × 10−6 7.3 × 105 
0.16 19.2 5.7 × 10−5 2.6 × 107 
0.18 21.6 1.1 × 10−3 4.9 × 108 
T, *T (K)r, *r (s−1)
0.06 7.2 2.0 × 10−21 9.5 × 10−10 
0.08 9.6 5.9 × 10−15 2.7 × 10−3 
0.10 12.0 4.3 × 10−11 2.0 × 101 
0.12 14.4 1.8 × 10−8 8.1 × 103 
0.14 16.8 1.6 × 10−6 7.3 × 105 
0.16 19.2 5.7 × 10−5 2.6 × 107 
0.18 21.6 1.1 × 10−3 4.9 × 108 

We have proposed an approach for computing eigenvalues, eigenvectors, and eigencurrents of stochastic networks representing energy landscapes. Step one of this approach is finding zero-temperature asymptotics for eigenvalues and eigenvectors using the earlier introduced20 efficient graph-cutting algorithm. Step two is the finite temperature continuation using the Rayleigh Quotient iteration supplemented with some precautions.

The proposed methodology has been applied to the LJ38 network26 with 71 887 states and 119 853 edges. The results turned out to be stereotype-breaking due to the absence of any significant spectral gaps. This happens due to the presence of a large number of high-lying liquid-like local minima separated by extremely high (perhaps artificial) potential barriers from the icosahedral and fcc funnels. However, the one-to-one correspondence between local minima and eigenvalues obtained by our algorithm renders the task of identifying important eigenvalues trivial.

We have explained the physical meaning of the eigencurrent and proved some of its properties. We demonstrated how it can be used for quantification and visualization of relaxation processes on the example of the LJ38 network. Each eigencurrent partitions the set of local minima into emitting and absorbing. The emission-absorption cut is the true transition state for the relaxation process described by the given eigencurrent.

We have compared the results of the analysis of the LJ38 network by means of the proposed spectral approach and by means of TPT. The outcomes of these two analyses are shown to be consistent.

Finally, we have discussed the question of metastability of LJ38 from the points of view of a mathematician and a chemical physicist. A connection with the experimental results for rare gas clusters has been made.

I am grateful to Professor E. Vanden-Eijnden for inspiring discussion on spectra and eigencurrents and useful suggestions about the present paper. I thank Professor D. Wales for providing the data of the LJ38 network and for valuable comments. This work is partially supported by DARPA YFA Grant No. N66001-12-1-4220 and National Science Foundation (NSF) Grant No. 1217118.

1.
Wales group web site
, http://www-wales.ch.cam.ac.uk.
2.
D. J.
Wales
,
Energy Landscapes: Applications to Clusters, Biomolecules and Glasses
(
Cambridge University Press
,
2003
).
3.
D. J.
Wales
,
Int. Rev. Chem. Phys.
25
,
237
(
2006
).
4.
C.
Schuette
, Habilitation thesis,
ZIB
, Berlin, Germany,
1999
, see http://publications.mi.fu-berlin.de/89/1/SC-99-18.pdf.
5.
W. C.
Swope
,
J. W.
Pitera
, and
F.
Suits
,
J. Phys. Chem. B
108
,
6582
(
2004
).
6.
S.
Elmer
,
S.
Park
, and
V. S.
Pande
,
J. Chem. Phys.
123
,
114902
(
2005
).
7.
J. D.
Chodera
,
N.
Singhal
,
V. S.
Pande
,
K. A.
Dill
, and
W. C.
Swope
,
J. Chem. Phys.
126
,
155101
(
2007
).
8.
F.
Noe
,
I.
Horenko
,
C.
Schuette
, and
J. C.
Smith
,
J. Chem. Phys.
126
,
155102
(
2007
).
9.
F.
Noe
,
C.
Schuette
,
E.
Vanden-Eijnden
,
L.
Reich
, and
T. R.
Weikl
,
Proc. Natl. Acad. Sci. U.S.A.
106
,
19011
(
2009
).
10.
J.-H.
Prinz
,
H.
Wu
,
M.
Sarich
,
B.
Keller
,
M.
Senne
,
M.
Held
,
J. D.
Chodera
,
C.
Schuette
, and
F.
No
,
J. Chem. Phys.
134
,
174105
(
2011
).
11.
C.
Schuette
,
F.
Noe
,
J.
Lu
,
M.
Sarich
, and
E.
Vanden-Eijnden
,
J. Chem. Phys.
134
,
204105
(
2011
).
12.
A.
Bovier
,
M.
Eckhoff
,
V.
Gayrard
, and
M.
Klein
,
Commun. Math. Phys.
228
,
219
(
2002
).
13.
A.
Bovier
,
Methods of Contemporary Mathematical Statistical Physics
(
Springer-Verlag
,
Berlin
,
2009
), pp.
177
222
.
14.
A.
Bovier
,
M.
Eckhoff
,
V.
Gayrard
, and
M.
Klein
,
J. Eur. Math. Soc.
6
,
399
(
2004
).
15.
A.
Bovier
,
M.
Eckhoff
,
V.
Gayrard
, and
M.
Klein
,
J. Eur. Math. Soc.
7
,
69
(
2005
).
16.
B.
Gaveau
and
L. S.
Schulman
,
J. Math. Phys.
39
,
1517
(
1998
).
17.
O.
Cepas
and
J.
Kurchan
,
Eur. Phys. J. B
2
,
221
(
1998
).
18.
C.
Schuette
,
W.
Huisinga
, and
S.
Meyn
,
Iutam Symposium on Nonlinear Stochastic Dynamics
(
Springer Science+Business Media
,
B.V.
,
2003
), pp.
71
82
.
19.
W.
Huisinga
,
S.
Meyn
, and
C.
Schuette
,
Ann. Appl. Probab.
14
,
419
(
2004
).
20.
M. K.
Cameron
, “
Computing the asymptotic spectrum for networks representing energy landscapes using the minimum spanning tree
,”
Networks Heterogeneous Media
9
,
383
(
2014
); e-print arXiv:1402.2869.
21.
M. K.
Cameron
and
E.
Vanden-Eijnden
,
J. Stat. Phys.
156
,
427
(
2014
).
22.
W.
E
and
E.
Vanden-Eijnden
,
J. Stat. Phys.
123
,
503
(
2006
).
23.
W.
E
and
E.
Vanden-Eijnden
,
Annu. Rev. Phys. Chem.
61
,
391
(
2010
).
24.
P.
Metzner
,
C.
Schuette
, and
E.
Vanden-Eijnden
,
SIAM Multiscale Model. Simul.
7
,
1192
(
2009
).
25.
J.
Lu
and
J.
Nolen
, “
Reactive trajectories and the transition path process
,”
Probab. Theory Relat. Fields
(published online
2014
).
26.
See http://www-wales.ch.cam.ac.uk/examples/PATHSAMPLE/ for the database for the Lennard-Jones-38 network.
27.
J. P. K.
Doye
,
M. A.
Miller
, and
D. J.
Wales
,
J. Chem. Phys.
110
,
6896
(
1999
).
28.
J. P. K.
Doye
and
D. J.
Wales
,
J. Chem. Phys.
109
,
8143
(
1998
).
29.
30.
M.
Picciani
,
M.
Athenes
,
J.
Kurchan
, and
J.
Taileur
,
J. Chem. Phys.
135
,
034108
(
2011
).
31.
J. P.
Neirotti
,
F.
Calvo
,
D. L.
Freeman
, and
J. D.
Doll
,
J. Chem. Phys.
112
,
10340
(
2000
).
32.
J. C.
Hamilton
,
D. J.
Siegel
,
B. P.
Uberuaga
, and
A. F.
Voter
, see http://www-personal.umich.edu/~djsiege/Energy_Storage_Lab/Publications_files/LJ38_v14.pdf for the preprint.
33.
M. K.
Cameron
,
J. Stat. Phys.
152
,
493
(
2013
).
34.
D. J.
Wales
and
J. P. K.
Doye
,
J. Phys. Chem. A
101
,
5111
(
1997
).
35.
V. A.
Mandelshtam
and
P. A.
Frantsuzov
,
J. Chem. Phys.
124
,
204511
(
2006
).
36.
D. J.
Wales
and
P.
Salamon
,
Proc. Natl. Acad. Sci. U.S.A.
111
,
617
(
2013
).
37.
J.
Farges
,
M. F.
de Feraudy
,
B.
Raoult
, and
G.
Torchet
,
Surf. Sci.
106
,
95
(
1981
).
38.
O.
Echt
,
K.
Sattler
, and
E.
Recknagel
,
Phys. Rev. Lett.
47
,
1121
(
1981
).
39.
I. A.
Harris
,
R. S.
Kidwell
, and
J. A.
Northby
,
Phys. Rev. Lett.
53
,
2390
(
1984
).
40.
I. A.
Harris
,
K. A.
Norman
,
R. V.
Mulkern
, and
J. A.
Northby
,
Chem. Phys. Lett.
130
,
316
(
1986
).
41.
O.
Echt
,
O.
Kandler
,
T.
Leisner
,
W.
Mlechle
, and
E.
Recknagel
,
J. Chem. Soc. Faraday Trans.
86
,
2411
(
1990
).
42.
S. I.
Kovalenko
,
D. D.
Solnyshkin
,
E. T.
Verkhovtseva
, and
V. V.
Eremenko
,
Chem. Phys. Lett.
250
,
309
(
1996
).
43.
S. I.
Kovalenko
,
D. D.
Solnyshkin
, and
E. T.
Verkhovtseva
,
Low Temp. Phys.
26
,
207
(
2000
).
44.
S.
Kakar
,
O.
Bjoerneholm
,
J.
Wiegelt
,
A. R. B.
de Castro
,
L.
Troeger
,
R.
Frahm
, and
T.
Moeller
,
Phys. Rev. Lett.
78
,
1675
(
1997
).
45.
J.
Kurchan
,
Six Out of Equilibrium Lectures
(
les Houches Summer School
,
2009
); e-print arXiv:0901.1271.
46.
E.
Vanden-Eijnden
,
An Introduction to Markov State Models and their Application to Long Timescale Molecular Simulation
(
Springer
,
Dordrecht
,
2014
), pp.
91
100
.
47.
J. W.
Demmel
,
Applied Numerical Linear Algebra
(
SIAM
,
1997
).
48.
L. N.
Trefethen
and
D.
Bau
,
Numerical Linear Algebra
(
SIAM
,
1997
).
49.
A. D.
Wentzel
,
Dokl. Akad. Nauk SSSR
202
,
19
(
1972
).
50.
A. D.
Wentzel
,
Sov. Math. Dokl.
13
,
65
(
1972
).
51.
M. I.
Freidlin
and
A. D.
Wentzell
,
Random Perturbations of Dynamical Systems
, 3rd ed. (
Springer-Verlag
,
2012
).
52.
R. K.
Ahuja
,
T. L.
Magnanti
, and
J. B.
Orlin
,
Network Flows: Theory, Algorithms, and Applications
(
Prentice Hall
,
New Jersey
,
1993
).
53.
J. B.
Kruskal
,
Proc. Am. Math. Soc.
7
,
48
(
1956
).
54.
C. J.
Tsai
and
K. D.
Jordan
,
J. Phys. Chem.
97
,
11227
(
1993
).
55.
M. A.
Miller
and
D. J.
Wales
,
J. Chem. Phys.
107
,
8568
(
1997
).
56.
M. I.
Freidlin
,
Sov. Math. Dokl.
18
,
1114
(
1977
).
57.
58.
J. P. K.
Doye
and
D. J.
Wales
,
J. Chem. Phys.
102
,
9659
(
1995
).
59.
W.
van der Waal
,
Phys. Rev. Lett.
76
,
1083
(
1996
).
60.
H.
Haberland
,
Clusters of Atoms and Molecules: Theory, Experiment, and Clusters of Atoms
(
Springer-Verlag
,
Berlin
,
1994
), pp.
374
395
.
61.
O. M.
Becker
and
M.
Karplus
,
J. Chem. Phys.
106
,
1495
(
1997
).
62.
D. J.
Wales
,
M. A.
Miller
, and
T. R.
Walsh
,
Nature (London)
394
,
758
(
1998
).