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.

## I. INTRODUCTION

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 10^{n}, *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 Process^{21,25} described by the committor, the transition rate, and the reactive current.

We apply the proposed methodology to the Lennard-Jones-38 (LJ_{38}) 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 LJ_{38} 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 LJ_{38} 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 *O*_{h} (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 *C*_{5v} (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 LJ_{38} stochastic network is a feasible model for the low-temperature dynamics of the LJ_{38} cluster.

The double-funnel structure of the energy landscape might make us expect that the smallest nonzero eigenvalue of the generator matrix of the LJ_{38} 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 time^{36} 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 LJ_{38} 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 LJ_{38}. 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 LJ_{38} network and the experimental results where rare gas clusters such as argon, krypton, and xenon have been studied using electron diffraction^{37–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 LJ_{38} network. The analyses of the LJ_{38} network by means of the proposed spectral approach and the Transition Path Theory^{21–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.

## II. SIGNIFICANCE OF SPECTRAL ANALYSIS

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

### A. The eigenstructure of networks with detailed balance

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

*i*and

*j*(

*i*≠

*j*) are connected by an edge, then

*L*

_{ij}is the pairwise transition rate from

*i*to

*j*, otherwise

*L*

_{ij}= 0. The diagonal entries

*L*

_{ii}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}],

The generator matrices of networks representing energy landscapes possess the detailed balance property, i.e., π_{i}*L*_{ij} = π_{j}*L*_{ji}, 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

where *P* is the diagonal matrix

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

Using Eq. (1) and the fact that *L*_{ii} = −∑_{j ≠ i}*L*_{ij} where *L*_{ij} ⩾ 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

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

where Φ = [ϕ^{0}, ϕ^{1}, …, ϕ^{N − 1}] is a matrix of right eigenvectors, normalized so that *P*^{1/2}Φ is orthogonal. Note that *P*^{1/2}Φ is the matrix of eigenvectors of *L*_{sym}. Therefore,

and in particular, ϕ^{0} = *e* and *P*ϕ^{0} = π.

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

where *p*_{0} is the initial distribution. Using Eqs. (2)–(4) we get

Note that

*k*-th eigen-component of

*p*(

*t*) remains significant only on the time interval

*p*(

*t*) → π as

*t*→ ∞.

### B. Interpretation of eigenvalues and eigenvectors

For *k* > 0, the left and right *k*th 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 *c*_{0} and *c*_{k} in Eq. (5) are nonzero. Imagine *n* ≫ 1 particles distributed in the stochastic network according to the equilibrium distribution π. Since

*k*> 0, the set of states

*S*can be divided into two parts,

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

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 *j* ∈ *S* gives the factor by which it is over- or underpopulated with respect to the equilibrium distribution π, as

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

### C. Eigencurrents

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 *i*th component of the probability distribution *p*(*t*) is

Plugging in expressions for *p*_{i} and *p*_{j} from Eq. (5) into Eq. (7) and using the detailed balance property π_{i}*L*_{ij} = π_{j}*L*_{ji}, we obtain

The collection of numbers

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

Hence, the eigencurrent

*k*th perturbation

*P*ϕ

^{k}along the edge (

*i*,

*j*) per unit time. In other words, if the system is originally distributed according to

*p*

_{0}= π + α

*P*ϕ

_{k}, then the current

*i*to

*j*and from

*j*to

*i*per unit time.

The eigencurrent

^{24}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

*i*∈

*S*we have

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

Let us partition the set of states *S* into

*i*,

*j*) where

^{21}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

*t*, i.e.,

(The symbol

*S*

^{′}⊂

*S*we have

Therefore, the eigencurrent over the cut (*S*^{′}, *S*^{″}) is

i.e., it is the total eigencurrent emitted in *S*^{′} per unit time. The maximum of the last sum is achieved 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}.

## III. ASYMPTOTIC ESTIMATES FOR EIGENVALUES AND EIGENVECTORS

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 Wentzel^{49,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

^{49,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

_{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 form^{29}

where *a*_{i} and *a*_{ij} are the prefactors for the local minimum *i* and the saddle (*i*, *j*) separating the local minima *i* and *j*, respectively, *V*_{i} and *V*_{ij} 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 eigenvectors^{12,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 *V*_{ij}, 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

We have proven^{20} 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

*k*th eigenvalue is

_{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

^{12,13}if the temperature is low enough.

## IV. COMPUTING EIGENVALUES AND EIGENVECTORS

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.

### A. Computing asymptotic eigenpairs

For simplicity, we make the following genericness assumption: all numbers *V*_{i}, *V*_{ij} and the differences *V*_{pq} − *V*_{i} 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 algorithm^{20} is computation of the minimum spanning tree

*i*,

*j*) is

*V*

_{ij}, 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

*V*

_{pq}along the unique path

*w**(

*i*,

*j*) in

*i*and

*j*is minimal possible.

^{52}

Central to the algorithm^{20} are the barrier function *u*(*i*) and the escape function *v*(*i*), *i* ∈ *S*. For a given set of sinks

*j*= 0, 1, …,

*k*,

where

*i*and

*u*(

*i*) is the height of the lowest possible highest saddle separating state

*i*from the set of sinks

*v*(

*i*) is the minimal possible maximal potential barrier to overcome in order to escape from

*i*to one of the sinks.

The algorithm^{20} proceeds as follows. The initial sink is

*u*and

*v*are computed with respect to this sink. The initial forest

*S*

_{0}=

*C*

_{0}=

*S*. Then for

*k*= 1, 2, …

Pick the state with the largest value of the escape function

*v*and make it the next sink$s^{\ast }_{k}$$sk*$.Find the edge with the maximal value

*V*_{pq}in the path connecting the new sink$s^{\ast }_{k}$$sk*$ with one of the existing sinks$s^{\ast }_j$$sj*$, 0 ⩽*j*⩽*k*− 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*=Tk\u22121*\u2216(pk*,qk*)$.Set

$\Delta _k:=V_{p^{\ast }_kq^{\ast }_k}-V_{s^{\ast }_{k}}$$\Delta k:=Vpk*qk*\u2212Vsk*$. Set*S*_{k}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*$.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*C*_{k}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 (LJ_{7})^{54,55} is worked out in Fig. 4.

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

the quasi-invariant set

*S*_{k};the sink

$s^{\ast }_k$$sk*$, which is the deepest minimum in the quasi-invariant set*S*_{k};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}}$$\Delta k=Vpk*qk*\u2212Vsk*$;Freidlin's cycle

^{33,51,56,57}*C*_{k}⊆*S*_{k}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 cycle^{20}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

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

^{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.

### B. Finite temperature continuation

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 *S*_{k} and *S*_{l} are both large, (ii) *S*_{k} ⊂ *S*_{l}, (iii) the number of states in *S*_{l}\*S*_{k} is small in comparison with the number of states in *S*_{k}, and (iv) the potential values at sinks

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

run the Rayleigh Quotient iteration and obtain an eigenpair (λ, ξ). Then we check whether (λ, ξ) is a correct eigenpair by calculating the dot product ξ^{T}(ψ^{k})_{0}. Since eigenvectors of *L*_{sym} 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.

## V. APPLICATION TO THE LJ_{38} NETWORK

The LJ_{38} dataset created by Wales's group^{26} 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 LJ_{38} 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.

### A. Settings and thermodynamics

Following Ref. 29 we define pairwise transition rates *L*_{ij} using the harmonic approximation. If there is a transition state in the database^{26} separating local minima *i* and *j*, then

where *O*_{i}, *V*_{i}, and

*i*;

*O*

_{ij},

*V*

_{ij}, and

*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*)),

*L*

_{ij}= 0. The equilibrium probability distribution is given by

To choose the range of temperatures at which we will investigate the LJ_{38} network, we use its critical temperature of the solid-solid phase transition as a reference. It has been established^{2,28,35} that the solid-solid phase transition, when fcc structures give place to icosahedral packings, occurs in the LJ_{38} 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 simulations^{35} 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 LJ_{38} 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 computed^{2,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 LJ_{38} cluster and its network approximation.

### B. Results

Originally, the results involving the zero-temperature asymptotics for the eigenvalues and eigenvectors of the LJ_{38} 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 LJ_{38} network by

_{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.

The reason for absence of significant spectral gaps lies in the fact that many liquid-like states in the LJ_{38} 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 *S*_{k} and *C*_{k}. The majority of those sets *S*_{k} contain less than 5 states. The values of the potential at the saddles

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 *V*_{max } on the potential energy. Then only the subnetwork *G*(*S*^{′}, *E*^{′}), where the sets of states and edges are given by

is considered. For *V*_{max } = 6 with respect to FCC (the lowest possible highest barrier separating ICO and FCC is^{27} 4.219), the resulting network contains 30 520 states and 71 750 edges. In this truncated network, Δ_{1}, corresponding to the sink

_{2}corresponding to the sink

*V*

_{max }= 5.5, the gap between Δ

_{1}and Δ

_{2}will increase to 1.00, while

*V*

_{max }= 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 *S*_{k} are plotted in Fig. 7. These eigenpairs are identified by^{20} *k* = 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

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 LJ_{38} 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.

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.

The relaxation process corresponding to the *k*th eigenpair is quantified by the *k*th 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.,

*i*,

*j*) with

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

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 *V*_{342, 354} = 4.219 and *V*_{958, 2831} = 4.186. The height of the saddle between minima 958 and 1, corresponding to another thick arrow across the emission-absorption cut, is^{26} *V*_{958, 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 (*O*_{ICO} = 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.

## VI. THE SPECTRAL ANALYSIS VERSUS TPT

In Ref. 21, the LJ_{38} 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.

Pick two disjoint subsets of states: a source set

*A*⊂*S*and a sink set*B*⊂*S*.- Solve the committor equationFor each state(20)\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}$\u2211j\u2208S\u2216(A\u222aB)Lijqj=0,q(A)=0,q(B)=1.$
*i*, the value of the committor*q*_{i}is the probability that the random walk starting at*i*will reach first*B*rather than*A*. - Calculate the reactive currentand analyze it. For example, one can consider its distribution in the isocommittor cuts and visualize the reactive tube with their aid. (For a given\begin{equation*}F^{R}_{ij}:=\pi _iL_{ij}(q_j-q_i)\end{equation*}$FijR:=\pi iLij(qj\u2212qi)$
*q*∈ [0, 1) the isocommittor cut^{21}*C*(*q*) is the collection of edges (*i*,*j*) such that*q*_{i}⩽*q*and*q*_{j}>*q*.) - Calculate the transition rate ν
_{R}, i.e., the average number of transitions from*A*to*B*per unit time, e.g., bywhere\begin{equation*}\nu _R=\sum _{(i,j)\in C(q)}F^{R}_{ij},\end{equation*}$\nu R=\u2211(i,j)\u2208C(q)FijR,$*C*(*q*) is an arbitrary isocommittor cut. Also calculate the rates*k*_{A, B}and*k*_{B, 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=\nu R\u2211i\u2208S\pi i(1\u2212qi),kA,B=\nu R\u2211i\u2208S\pi 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

to the reactive current

*k*

_{A, B}to λ

_{k}. We will do it on the example of the LJ

_{38}network.

As the temperature tends to zero, the right eigenvector ϕ^{k} tends^{12,13} to the committor *q*^{AB} with the source set

*S*

_{k}. 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

^{(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 Process^{21,25} in TPT and the relaxation process are proportional to π_{i}*L*_{ij}(*q*_{j} − *q*_{i}) and

*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

*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

*A*and

*B*. The total flux of the eigencurrent is maximal through the emission-absorption cut, which is the cut separating the states with

*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

_{38}network at low temperatures, say

*T*⩽ 0.12. The distributions of the eigencurrent

*q*= 0.5 are similar (compare Fig. 8 and 9 in Ref. 21).

At low temperatures, both the transition rate *k*_{ICO, 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,

The theoretical asymptotic value of the exponent is Δ = 3.543. As the temperature increases, the graphs of *k*_{ICO, FCC} and λ(ICO) diverge (compare Fig. 7 with Fig. 5 in Ref. 21). λ(ICO) grows faster than *k*_{ICO, 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 preconditioning^{47} 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.

## VII. Is LJ_{38} METASTABLE?

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 LJ_{38} is metastable from two points of view: of a mathematician and of a chemical physicist.

A detailed discussion on the metastability of the LJ_{38} 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

_{38}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 diffraction^{37–43} and extended x-ray absorption fine structure spectroscopy^{44} 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 from^{44} *N* = 200 to^{37,40} *N* = 800 atoms. It has been hypothesized by van der Waal^{59} 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 *F*^{ICO}. Using parameters appropriate for argon^{2} 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.

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 × 10^{1} |

0.12 | 14.4 | 1.8 × 10^{−8} | 8.1 × 10^{3} |

0.14 | 16.8 | 1.6 × 10^{−6} | 7.3 × 10^{5} |

0.16 | 19.2 | 5.7 × 10^{−5} | 2.6 × 10^{7} |

0.18 | 21.6 | 1.1 × 10^{−3} | 4.9 × 10^{8} |

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 × 10^{1} |

0.12 | 14.4 | 1.8 × 10^{−8} | 8.1 × 10^{3} |

0.14 | 16.8 | 1.6 × 10^{−6} | 7.3 × 10^{5} |

0.16 | 19.2 | 5.7 × 10^{−5} | 2.6 × 10^{7} |

0.18 | 21.6 | 1.1 × 10^{−3} | 4.9 × 10^{8} |

## VIII. CONCLUSION

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 introduced^{20} 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 LJ_{38} network^{26} 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 LJ_{38} 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 LJ_{38} 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 LJ_{38} 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.

## ACKNOWLEDGMENTS

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 LJ_{38} 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.