text:
# Large well-relaxed models of vitreous silica, coordination numbers and entropy

## I Introduction

If liquid silica (SiO2) is cooled below its melting temperature, it usually does not crystallize, but stays a supercooled liquid for an extended period of time. Upon further cooling, to below a certain temperature known as the glass transition temperature $`T_g`$, it becomes highly viscous and shows many properties of ordinary solids, but stays disordered and does not evolve to the thermodynamically stable crystalline phase, at least on earthly time scales. The study of silica in its disordered phase (vitreous silica) is hampered by a lack of knowledge on its microscopic structure. Experimental techniques alone, such as neutron and X-ray scattering, which can uniquely resolve the structure of a crystal, are less successful in this case because of the lack of a repeating unit cell. Typically, only information averaged over many atoms, such as the radial distribution function, can be obtained from experiments.

An alternative method to study vitreous silica is through computer simulations. The usual computational approach to generate vitreous silica structures is to simulate a quench from the melt within the framework of molecular dynamics (MD). This resembles closely the experimental process used to prepare vitreous silica, except that the typical computational cooling rates are about ten orders of magnitude faster than experimental ones. In agreement with experiment, these MD simulations result in structures of vitreous silica in which almost all silicon atoms are bonded to four oxygen atoms and almost all oxygen atoms are bonded to two silicon atoms, without any sign of long-range order vollmayr1996a . This in turn is in agreement with the continuous random network (CRN) model, proposed by Zachariasen zachariasen1932a .

The high cooling rates used in MD simulations result in structures with a strain higher than observed in experiments. This shows up, for instance, as an anomalously high density of coordination defects and a larger spread in bond lengths and bond angles. In this paper we present an alternative approach to generate models of vitreous silica. In contrast to MD, we do not attempt to simulate the details of the dynamics of the melt-quenching process. Instead, we use a Monte Carlo (MC) scheme with an artificial dynamics consisting of bond transpositions. Our scheme is similar in spirit to the one used by Tu et al. in their study of silica and silica sub-oxides tu1998a ; tu2000a , which in turn is based on the algorithm of Wooten, Winer and Weaire (WWW) for the generation of four-fold coordinated CRNs wooten1985a ; djordjevic1995a .

The outline of this paper is as follows. We first describe the MC scheme used by Tu et al. We then move on to describe a number of optimizations to this algorithm. These optimizations allow us to generate large and well-relaxed silica networks containing up to 300,000 atoms. The properties of these networks are discussed and compared to MD-prepared networks and experiment in Section VI. We end with a number of conclusions in Section VIII.

## II Bond-switching algorithm for silica

In the approach followed by Tu and co-workers, a vitreous silica network consists of $`N_s`$ silicon atoms and $`N_o=2N_s`$ oxygen atoms. The total number of atoms in the network thus equals $`N=N_s+N_o`$. At all times, an explicit list of bonds is maintained which uniquely determines the connectivity or topology of the network. In the list of bonds, each silicon atom is bonded to four oxygen atoms, and each oxygen atom is bonded to two silicon atoms. The total number of bonds in the list thus equals $`4N_s`$.

To calculate energy and forces the following potential is used:

$`E_{TT}`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{i}{}}k_b(b_ib_0)^2`$ (1)
$`+`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{i,j}{}}k_\theta (\mathrm{cos}\theta _{ij}\mathrm{cos}\theta _0)^2,`$
where the first summation runs over the list of bonds, and the second summation over all pairs of bonds that share one atom; $`b_i`$ represents the length of the $`i`$-th bond, and $`\theta _{ij}`$ is the angle between bonds $`i`$ and $`j`$ (which share one atom). The potential parameters for the two-body Si-O interaction are given by $`k_b=27.0`$ eV and $`b_0=1.60`$ Å. For the three-body O-Si-O interaction common values are $`k_\theta =4.32`$ eV and $`\mathrm{cos}\theta _0=1/3`$; for Si-O-Si interactions these values are $`k_\theta =0.75`$ eV and $`\mathrm{cos}\theta _0=1`$. Since the number of bonds in the network grows linearly with the number of atoms it contains, each energy evaluation is an $`O(N)`$ operation.

The network is evolved by making explicit changes in the list of bonds, each time followed by a local energy minimization with respect to the geometric coordinates of all atoms (using for example the method of steepest descent). The changes in the list of bonds consist of bond transpositions as illustrated in Fig. 1. Each bond transposition is accepted with the Metropolis acceptance probability metropolis1953a :

$$P_a=\mathrm{min}[1,\mathrm{exp}\left(\frac{E_bE_f}{k_BT}\right)],$$
(2)
where $`k_B`$ is the Boltzmann constant, $`T`$ is the temperature, and $`E_b`$ and $`E_f`$ are the (minimized) energies of the network before and after the proposed bond transposition.

## III Optimized bond-switching algorithm for silica

We have optimized the bond-switching algorithm of Tu et al. for the generation of large and well relaxed silica networks. These optimizations are similar in spirit to the optimizations used in the scalable WWW algorithm barkema2000a ; vink2001c .

The first optimization is aimed at reducing the CPU time spent on rejected bond transpositions. After a bond transposition in the original algorithm, the energy of the network is always completely minimized. After the minimization, the bond transposition is either accepted or rejected based on the Metropolis probability.

In contrast, we reformulate the Metropolis algorithm. We first determine a threshold energy $`E_t=Ek_BT\mathrm{ln}(1r)`$, where $`r`$ is a random number uniformly drawn from the interval $`[0,1`$. We then proceed with the minimization procedure. During minimization, the converged energy is continuously estimated. Assuming that the total energy is harmonic around the minimum with curvature $`c`$, the energy $`E(\stackrel{}{R}_0)`$ at the (unknown) minimum $`\stackrel{}{R}_0`$ can be estimated from a nearby configuration $`\stackrel{}{R}`$ as $`E(\stackrel{}{R}_0)E(\stackrel{}{R})|\stackrel{}{F}|^2/c`$, with $`|\stackrel{}{F}|`$ the magnitude of the total force at $`\stackrel{}{R}`$. Minimization is aborted as soon as it becomes clear that the threshold energy will not be reached and the bond transposition will be rejected. This leads to a large reduction in the number of force evaluations associated with rejected bond transpositions.

The second optimization is aimed at exploiting the local nature of the bond transposition depicted in Fig. 1. Immediately after a bond transposition, only a small cluster of atoms in the network will experience a significant force. This cluster consists of the atoms directly involved in the bond transposition, marked $`\{S1,S2,O1,O2,O3\}`$ in Fig. 1, and of nearby atoms (typically extending to the fourth neighbor shell around these five transposition atoms). The number of atoms in such a cluster is about 400. It therefore suffices to calculate the force locally (only for the 400 or so atoms inside the cluster) rather than globally (for all the atoms in the network). Calculating the force on a cluster of atoms is an $`O(1)`$ operation, which means that it is independent of the total number of atoms in the network. Local force calculations are therefore much cheaper than global $`O(N)`$ force calculations.

Already after a converged local minimization, it is often clear that the threshold energy will not be reached, and the bond transposition can be rejected. If this is not the case, a few global minimization steps are required additionally, usually resulting in an accepted bond transposition. The combination of these two improvements reduces the computational effort from $`O(N)`$ operations per attempted bond transposition to $`O(1)`$ operations per attempted bond transposition, plus a few $`O(N)`$ operations per accepted bond transposition. Since common acceptance probabilities are well below 1%, the speed-up for large systems is about two orders of magnitude.

## IV Initial configurations

The optimized bond-switching algorithm for silica requires an initial network, which should be free of coordination defects and crystalline regions. As noted in Ref. thorpe1986a, , four-fold coordinated CRNs decorated with one oxygen atom on each bond already provide structures that compare reasonably well to vitreous silica. As initial network we therefore use a periodic, four-fold coordinated CRN generated as described in Refs. barkema2000a ; vink2001c . This CRN serves as the silicon backbone. Next, one oxygen atom is placed on the center of each Si-Si bond to obtain a properly coordinated silica network. The simulation volume is then scaled to obtain the experimental density of silica $`\rho =2.20`$ g cm-3 mazurin1983a .

## V BKS sampling

A computationally more expensive, but also more realistic potential for silica is the BKS potential beest1990a :

$$E_{BKS}=\underset{i<j}{}\left(\frac{q_iq_j}{r_{ij}}+A_{ij}\mathrm{exp}\left(B_{ij}r_{ij}\right)\frac{C_{ij}}{r_{ij}^6}\right),$$
(3)
where the sum runs over all pairs of atoms in the system, $`r_{ij}`$ is the distance between atoms $`i`$ and $`j`$, and $`q_i`$ is the charge of atom $`i`$. The atomic charges and the values of the potential parameters $`A_{ij}`$, $`B_{ij}`$ and $`C_{ij}`$ are also listed in Ref. beest1990a, . We use the fast multipole code developed at Duke University rankin1995a  to compute the Coulomb and van der Waals interactions. The exponential term is truncated at $`5.50`$ Å and then shifted, as described in Ref. vollmayr1996a, .

The bond-switching algorithm described above generates networks that are local energy minima according to the potential of Tu et al. The algorithm can be changed such that these networks will be sampled from the Boltzmann distribution according to the BKS potential at temperature $`T`$. To this end, the acceptance probability of Eq. (2) is replaced by:

$`P`$ $`=`$ $`\mathrm{min}[1,\mathrm{exp}\left({\displaystyle \frac{\mathrm{\Delta }E_{TT}}{k_BT^{}}}\right)]`$ (4)
$`\times \mathrm{min}[1,\mathrm{exp}\left({\displaystyle \frac{\mathrm{\Delta }E_{TT}}{k_BT^{}}}{\displaystyle \frac{\mathrm{\Delta }E_{BKS}}{k_BT}}\right)],`$
where $`\mathrm{\Delta }E_{TT}`$ and $`\mathrm{\Delta }E_{BKS}`$ are the energy differences between the networks before and after the bond transposition, calculated with the potential of Tu et al. and the BKS potential, respectively. The algorithm satisfies detailed balance, also when the fictitious temperatures $`T`$ and $`T^{}`$ are different. In the expression above, the first factor biases moves according to the potential of Tu et al. at temperature $`T^{}`$, after which the second factor corrects for the difference between this potential at temperature $`T^{}`$ and the BKS potential at temperature $`T`$.

With the above probability, bond transpositions that are accepted with the potential of Tu et al. undergo an additional accept/reject decision based on the BKS potential. By increasing the temperature $`T^{}`$ we can increase the number of moves that will be accepted under the potential of Tu et al. In the limit $`T^{}\mathrm{}`$, all such moves are accepted. In this limit, the accept/reject decision is based solely on the BKS potential and $`T`$. While correct, this situation is undesirable because BKS energy evaluations are expensive. In practice, a lower $`T^{}`$ is used and the temperatures $`\{T,T^{}\}`$ are tuned to maximize the decrease of the BKS energy of the network as a function of simulation time.

## VI Results

We have used the optimized bond-switching algorithm described in Section III to generate vitreous silica networks containing 3000 atoms (‘3k’), 60,000 atoms (‘60k’) and 300,000 atoms (‘300k’). For the 3000-atom and the 60,000-atom networks the bond-switching algorithm was used in conjunction with BKS sampling at temperatures $`k_BT=0.30`$ eV and $`k_BT^{}=0.70`$ eV (see Section V). For the 300,000-atom network, BKS sampling proved to be too computationally demanding and was not implemented; for this network the Metropolis probability of Eq. (2) was used with $`k_BT=0.15`$ eV and the potential of Tu et al.

All three networks were evolved with approximately 10 attempted bond transpositions per atom from their starting configurations. The energy of the resulting networks was then minimized in a single quench using the method of steepest descent and the BKS potential without volume optimization. We observed that if volume optimization were used with the BKS potential, the density of the networks becomes unphysically high, typically in the range 2.27-2.37 g cm-3. Similar high densities are also reported in a number of MD studies vollmayr1996a ; horbach1999a .

In Fig. 2 we show the radial distribution functions (RDFs) $`g_{SiSi}(r)`$, $`g_{SiO}(r)`$, and $`g_{OO}(r)`$ for model ‘60k’. The RDFs are density normalized such that  $`lim_r\mathrm{}g_{\alpha \beta }(r)=1`$. The RDFs for models ‘3k’ and ‘300k’ are similar (not shown).

The Si-O RDF was used to extract a cutoff for the Si-O bond length at $`r_c=1.80`$ Å, where $`r_c`$ is taken at the minimum between the first and the second peak. Based on this $`r_c`$ we find structural properties as reported in Table 1. For comparison, we also report results obtained with MD using the same BKS potential. Compared to MD, we observe that our networks are better relaxed, as is evident from the smaller variations in the O-Si-O and Si-O-Si bond angles. Moreover, our networks are nearly free of coordination defects.

A more stringent test is to compare the properties of our networks to experimental data. The crucial quantity to consider is the total correlation function $`T(r)`$ as obtained in X-ray or neutron scattering experiments. In these experiments, the quantity that can be measured directly is the total interference function $`Qi(Q)`$; the correlation function $`T(r)`$ is obtained via a Fourier transform of $`Qi(Q)`$. In practice, the resolution in $`T(r)`$ is determined by the maximum attainable inverse wavelength $`Q_{max}`$ in the experiment. As discussed in Ref. wright1994a, , peaks in $`T(r)`$ are significantly broadened because $`Q_{max}`$ is not infinite in real experiments. In addition to Fourier broadening, there is thermal broadening because experiments are typically carried out at room temperature. For a meaningful comparison it is essential to take these factors into account.

To capture the effect of thermal broadening, all three networks are thermalized over approximately 1 ps with MD using the BKS potential in the NVE-ensemble at $`T300`$ K. The structural properties of the thermalized networks are summarized and compared to experiment in Table 2. We observe good agreement on the Si-O bond length, the variation in the bond length and the mean bond angles. Compared to experiment, our networks slightly overestimate the variation in the O-Si-O bond angle. This may indicate that experimental vitreous silica is more relaxed than our networks. Experimental estimates for the variation in the Si-O-Si bond angle range from 12.8 degrees grandinetti2002a  to 36 degrees mozzi1969a . Compared to these data our networks coincide with the higher value. The coordination numbers reported in Table 2 indicate that the vast majority of atoms is still properly coordinated after thermalization.

Next, Fourier broadening is taken into account. In the absence of Fourier broadening, the total correlation function of a network is simply a weighted sum of the RDFs:

$`T(r)`$ $`=`$ $`{\displaystyle \frac{r}{w_{SiSi}+w_{SiO}+w_{OO}}}[w_{SiSi}\times g_{SiSi}(r)`$ (5)
$`+w_{SiO}\times g_{SiO}(r)+w_{OO}\times g_{OO}(r)],`$
in which each RDF is normalized such that $`lim_r\mathrm{}T(r)=r`$. We focus on neutron scattering, in which case the weights are given by: $`w_{SiSi}=0.1722`$, $`w_{SiO}=2\times 0.4817`$ and $`w_{OO}=4\times 0.3370`$, see Ref. wright1994a, . These weights include the neutron scattering lengths of silicon and oxygen atoms; the factors of 2 and 4 account for the chemical composition of silica (two oxygen atoms for every silicon atom).

To obtain the Fourier-broadened correlation function $`T_b(r)`$, the correlation function of Eq. (5) is transformed to obtain the interference function $`Qi(Q)`$:

$$Qi(Q)=_{r=0}^{r=r_{max}}\left[T(r)r\right]\mathrm{sin}(Qr)𝑑r,$$
(6)
where $`r_{max}`$ is half the edge of the cubic simulation volume of the model. Next, we transform the interference function back to obtain $`T_b(r)`$:

$$T_b(r)=r+\frac{2}{\pi }_{Q=0}^{Q=Q_{max}}M(Q)Qi(Q)\mathrm{sin}(Qr)𝑑Q,$$
(7)
using the experimentally relevant $`Q_{max}`$. The function $`M(Q)`$ is an (optional) modification function commonly used to reduce Fourier artifacts. A popular choice is the Lorch modification function lorch1969a  given by:

$$M(Q)=\frac{Q_{max}}{\pi Q}\mathrm{sin}\left(\frac{\pi Q}{Q_{max}}\right).$$
(8)
The neutron $`T(r)`$ of vitreous silica was accurately measured by Grimley, Wright and co-workers in 1991 grimley1990a ; wright1991a  using high-energy neutrons with $`Q_{max}=45.2`$ Å-1. For each of the three thermalized networks we determined the Fourier-broadened correlation function $`T_b(r)`$ using the same value for $`Q_{max}`$ and the Lorch modification function. The experimental correlation function and the correlation function of the thermalized ’60k’ network are compared to each other in Fig. 3. We observe excellent agreement. As a more quantitative measure of the agreement between model and experiment we consider the $`R_\chi `$ factor wright1994a  defined as:

$$R_\chi ^2=\frac{\left[T_{exp}(r)T_{model}(r)\right]^2𝑑r}{T_{exp}^2(r)𝑑r}.$$
(9)
The discrepancy between the experimental $`T(r)`$ and $`T_b(r)`$ of the thermalized ‘60k’ network, over the range $`1.0r10.0`$ Å, equals $`R_\chi =5.1\%`$. For the thermalized networks ‘3k’ and ‘300k’ we obtain $`R_\chi `$ factors of 4.7% and 4.9%, respectively, over the same range.

In the top frame of Fig. 4 we compare the total interference function of the thermalized ‘60k’ network to experiment. There is excellent agreement on the overall peak positions and also on the damping of the interference function for large $`Q`$. We observe some discrepancy at $`Q15`$ Å-1, where the small peak visible in the experimental data is not reproduced by our network. The discrepancy seems to be temperature-related: the lower frame of Fig. 4 shows that the network at zero temperature does reproduce the peak correctly.

We observe in Fig. 3 that the total correlation function is significantly broadened by Fourier and thermal effects, in particular the first peak at $`r1.6`$ Å. We have quantified the broadening of this peak by fitting it to the form:

$$P(r)=\frac{1}{\sqrt{2\pi }}\frac{\eta }{r\sigma }\mathrm{exp}\left[\frac{(rr_0)^2}{2\sigma ^2}\right],$$
(10)
with fit parameters $`\{r_0,\sigma ,\eta \}`$. This follows the same procedure used in Ref. grimley1990a,  to obtain bond lengths and coordination numbers from experimentally obtained correlation functions. The parameters $`r_0`$ and $`\sigma `$ provide estimates for the mean Si-O bond length and its variation. The parameter $`\eta `$ can be used to extract the mean coordination number $`z`$ of silicon atoms:

$$z=4\pi \eta \rho _O\frac{w_{SiSi}+w_{SiO}+w_{OO}}{w_{SiO}},$$
(11)
with $`\rho _O`$ the number density of oxygen atoms, $`w_{\alpha \beta }`$ the neutron weights and $`T(r)`$ normalized such that $`lim_r\mathrm{}T(r)=r`$.

The effects of Fourier and thermal broadening on the parameters $`\{r_0,\sigma ,z\}`$ are reported in Table 3. As expected, the broadening is significant: the combined effect of Fourier and thermal broadening will boost $`\sigma `$ by a factor of over five.

It is also interesting to compare the silicon coordination number obtained directly from the atomic coordinates of a network, to that obtained from a fit to the correlation function of that network. To this end the data of the networks at 300 K without Fourier broadening in Table 3 are compared to Table 2. The fitting procedure accurately predicts the Si-O bond length and its variation $`\sigma `$, but not the coordination number which we know to be 4.000. This result is important because it shows that in the ideal situation free of Fourier broadening, fits based on Eq. (10) will systematically underestimate the coordination number. The cause of this underestimation is due to the non-gaussian nature of the peak in $`rT(r)`$: fitting with a more elaborate function, for instance with a sum of two gaussians, yields the correct coordination number. The table also shows that fits with Eq. (10) applied to Fourier-broadened data tend to overestimate the coordination numbers.

## VII Configurational entropy

Recently, we developed a method to determine the configurational entropy of a network vink2002a . The method requires only the atomic coordinates and a list of bonds between particles, for a single well-relaxed configuration. In the original paper, the method was applied to a silica network consisting of 3000 atoms. Starting with the atomic coordinates of this network, the list of bonds was constructed using $`r_c=1.80`$ Å for the Si-O bond-length cutoff. Each Si-O-Si bridge was replaced by a single Si-Si bond. We then applied our method to the resulting silicon backbone network to obtain for the configurational entropy $`s=0.88`$ $`k_B`$ per silicon atom.

The accuracy of the method is sensitive to the size of the network. As mentioned in the original paper, the configurational entropy of vitreous silica reported above is most likely an underestimate due to the limited size of that network. The large networks generated in this work allow us to quantify these finite size effects. They also allow us to determine the minimum size of a network that is required to accurately measure the configurational entropy.

To this end, we use our method to determine the entropy of networks ‘3k’, ‘60k’ and ‘300k’ in two different ways: first, by using the entire simulation volume of the network; secondly, by using only half the volume of the simulation cell. The first procedure yields entropies of 0.83, 0.99 and 1.04 $`k_B`$ per silicon atom, respectively, while the second procedure yields 0.78, 0.97 and 1.03 $`k_B`$ per silicon atom, respectively. This illustrates that finite-size effects are significant for the smallest network, but small for the large ones. Most of the difference in entropy between the three networks is due to the varying degree of relaxation in these networks.

## VIII Conclusions

In summary, we have presented a Monte Carlo-based approach for the generation of well-relaxed networks of vitreous silica, based on earlier work of Tu et al. With this method, networks containing 3000, 60,000 and 300,000 atoms are generated. Compared to networks generated with molecular dynamics, our networks have smaller bond-angle variations and are nearly defect-free, indicating they are better relaxed.

The total correlation function $`T(r)`$ of our networks is in excellent agreement with neutron scattering experiments (with $`R_\chi `$ factors around 5%), provided that thermal effects and the maximum experimental inverse wavelength $`Q_{max}`$ are included in the comparison. We also determined the silicon coordination number of our networks by fitting a gaussian to $`rT(r)`$, as is commonly done in experiments. We observe that the results of this procedure are biased to lower numbers by the non-gaussian nature of the peaks, and to higher numbers due to finite value of $`Q_{max}`$. The configurational entropy of vitreous silica was determined to be 0.99 $`k_B`$ and 1.04 $`k_B`$ per silicon atom, for the networks containing 60,000 and 300,000 atoms, respectively.

Upon request, we will communicate the atomic coordinates of our networks.


format: text

text:
# NONEXTENSIVE INFORMATION ENTROPY FOR STOCHASTIC NETWORKS

## Abstract

Nature is full of random networks of complex topology describing such apparently disparate systems as biological, economical or informatical ones. Their most characteristic feature is the apparent scale-free character of interconnections between nodes. Using an information theory approach, we show that maximalization of information entropy leads to a wide spectrum of possible types of distributions including, in the case of nonextensive information entropy, the power-like scale-free distributions characteristic of complex systems.

PACS numbers: 89.75.-k 89.70.+c 24.60.-k
Keywords: Complex networks, Information theory, Nonextensive statistics

Random networks have recently found applications in the description of complex systems in different, apparently very disparate branches of modern science such as, for example, molecular biology, sociology, economy and computer science . For example, living organisms form huge genetic networks the nodes of which are proteins and links represent the corresponding chemical interactions . A similarly big network is formed by the nervous system the nodes of which are connected by axions . Comparable complexity show networks existing in the sociological systems in which nodes are countries, organizations or single persons whereas links characterize their mutual interactions , in the world of finances and computer networks (with World Wide Web being the most known example where nodes are HTML documents connected via hiper-links URL ). For most recent reviews of random networks see .

Analysis of different random networks clearly indicate that the probability $`P(k)`$ of joining a given node with other nodes is described by the power law $`P(k)k^\gamma `$ . For example, the most convincing analyses of computer networks with over $`800`$ milion nodes  lead to a power-like distribution of $`P(k)`$ with exponent equal $`\gamma 2.1÷2.45`$. This contradicts the existing models of random networks  predicting instead exponential distributions: $`P(k)\mathrm{exp}(k)`$. The most popular model (ER) dealing with a fixed number of nodes $`N`$ was proposed in , where the Poisson distribution was advocated to be used for probability that a given node has $`k`$ links (with the mean number of connections being $`\lambda _0`$) 111Notice that in the limit of large $`k`$ distribution (1) can be approximated by $`P(k)=(2\pi k)^{1/2}\left(\lambda _0/k\right)^k\mathrm{exp}\left(k\lambda _0\right)`$, whereas for large values of $`\lambda _0`$ it becomes a gaussian distribution, $`P(k)=\left(2\pi \lambda _0\right)^{1/2}\mathrm{exp}\left[\left(k\lambda _0\right)^2/2\lambda _0\right]`$.,

$$P(k)=\frac{\lambda _0^ke^{\lambda _0}}{k!}.$$
(1)
However, in order to get the observed power-like form of $`P(k)`$ one has to allow for growing $`N`$ and replace the democratic law of attachement a new link used in deriving (1) by a preferential one. This means that distribution $`P(k)`$ is determined by the dynamics of the growth of network . Starting from a small number $`m_0`$ of nodes, adding in each time step a new node with $`mm_0`$ possible connections and assuming that this new node joins the already existing nodes with equal, $`k_i`$-independent, probability $`\mathrm{\Pi }(k_i)=\frac{1}{m_0+t1}`$, the evolution (growth) of network is described by the following equation :

$$\frac{k_i}{t}=m\mathrm{\Pi }(k_i)=\frac{m}{m_0+t1}$$
(2)
leading for long times $`t`$ to exponential stationary distribution:

$$P(k)=\frac{1}{m}\mathrm{exp}\left(\frac{k}{m}\right).$$
(3)
On the other hand, assuming that probability $`\mathrm{\Pi }(k_i)`$ is selective, for example that $`\mathrm{\Pi }(k_i)=\frac{k_i}{_{j=1}^{m_0+t1}k_j}=\frac{k_i}{2mt}`$, one gets instead asymptotically a simple power law for $`P(k)`$:

$$P(k)=\frac{2m^2t}{m_0+t}k^3k^3.$$
(4)
Apparently such distribution with universal exponent $`\gamma =3`$ shows up in different situations (and under different names). As Pareto distribution  it describes the growth of the wealth of persons living in stable economical systems, as Zipf’s law  it is applied in linguistic and it also describes the distribution of the citations of the scientific works .

In the limit of large $`t`$, i.e., when stationary state already develops, this problem can be also studied from the information theory point of view. In it one asks the following question : what is the informational content of data represented by distributions $`P(k)`$? In other words, what is the minimal number of parameters needed to reproduce a given shape of $`P(k)`$? The question asked in such approach is: suppose that we know only that a network we are interested in, which we would like to describe by $`P(k)`$, leads to some mean value of $`k`$, i.e., we know that:

$$\underset{k=1}{\overset{\mathrm{}}{}}P(k)=\mathrm{\hspace{0.17em}1}\mathrm{and}k=\underset{k=1}{\overset{\mathrm{}}{}}kP(k)=\lambda _0=const,$$
(5)
what would be then the most probable and least biased form of $`P(k)`$ in such a case (i.e., describing the existing data and given entirely in terms of $`\lambda _0`$)? To answer this question one maximalizes the corresponding information entropy associated with probability distribution $`P(k)`$ under constraints (5). The usual form of such entropy is Shannon entropy ,

$$S=\underset{k=1}{\overset{\mathrm{}}{}}P(k)\mathrm{ln}P(k).$$
(6)
The conditions (5) representing our a priori knowledge of the problem lead to the exponential probability distribution

$$P(k)=\frac{1}{\lambda _0}\mathrm{exp}\left(\frac{k}{\lambda _0}\right),$$
(7)
closely resembling eq.(3) 222It is interesing to mention that additional knowledge that all entities represented by $`k`$ are indistinguishable, which results in the necessity of introducing in the correspodning sumations the additional weight factor $`1/k!`$, changes the above distribution to Poisson distribution of the ER model mentioned above, cf., eq.(1)..

There are, however, many systems with properties preventing the use of Shannon type of information entropy and calling for its generalization. For example, they posses some intrinsic fluctuations resulting in the whole spectrum of parametr $`\lambda `$, $`P(\lambda )`$, with $`\lambda _0=\lambda `$ being only its mean value  or they develop some correlations introducing memory effects, cf.  (in statistical physics such situation leads to the necessity of departing from the use of the usual Boltzmann-Gibbs statistics in favour of some sort of generalized one ). Out of many possible generalizations we shall use in this note the nonextensive Tsallis entropy defined as :

$$S_q=\frac{1}{1q}\left\{1\underset{k=1}{\overset{\mathrm{}}{}}\left[P(k)\right]^q\right\}.$$
(8)
It can be regarded as a minimal (i.e., one parameter) extension of Shannon entropy (6), to which it reduces when $`q1`$. Parameter $`q`$ describes therefore summarily all effects preventing the use of Shannon entropy mentioned above.

Using $`S_q`$ as a measure of information about our system, i.e., maximalizing $`S_q`$ with constraints (equivalent to (5) above):

$$\underset{k=1}{\overset{\mathrm{}}{}}P(k)=\mathrm{\hspace{0.17em}1}\mathrm{and}k_q=\frac{_{k=1}^{\mathrm{}}k[P(k)]^q}{_{k=1}^{\mathrm{}}[P(k)]^q}=\lambda _0=const,$$
(9)
one obtains as result a power-like distribution of the form:

$$P(k)=P_q(k)=C\left[1(1q)\frac{k}{\lambda _0}\right]^{\frac{q}{1q}},$$
(10)
where $`C=1/_{k=1}^{\mathrm{}}[1(1q)k/\lambda _0]^{q/(1q)}=1/\lambda _0`$ is normalization333It should be stressed that maximalization of entropy provides us in this case only with the shape of distribution $`P(k)`$, eq. (10), and gives no information on the particular values of parameters $`\lambda _0`$ and $`q`$. Only knowledge of moments $`k`$ and $`Var(k)`$ of $`P_q(k)`$, as given by eq. (11), allows for determining these two parameters.. In this case

$$k=\frac{\lambda _0}{(2q)}\mathrm{and}Var(k)=\frac{\lambda _0^2}{(32q)(2q)^2}.$$
(11)
Notice that for $`q1`$ this distribution becomes exponential, as in eq.(7). On the other hand, for large values of $`k`$, $`k>>\lambda _0/(q1)`$, it becomes a power-like distribution of the form

$$P_q(k)k^\gamma \mathrm{with}\gamma =\frac{q}{q1},$$
(12)
i.e., our distribution becomes in this limit a scale-free one. It is easy to check that if we demand that $`k<\mathrm{}`$ then $`q<2`$. It is interesting to note at this point that $`\gamma =3`$ in eq.(4) corresponds precisely to $`q=3/2`$ at which variation $`Var(k)`$ diverges.

In Fig. 1 we show (as example) distribution of the number of connections in the WWW network  with $`325729`$ nodes and the mean values of connections equal $`k=5.46`$ fitted by $`P_q(k)`$ as given by eq.(10) with parameters $`\lambda _0=1.91`$ and $`q=1.65`$ 444Notice that although in fitting procedure both parameters were varied independently, they are connected via distribution moments $`k`$ and $`Var(k)`$, cf. (11).. Notice that eq.(10) describes the whole range of $`k`$ whereas the purely power-like distribution $`k^\gamma `$ with $`\gamma =q/(q1)=2.54`$ occurs only for large values of $`k`$. In the spirit of information theory this result can therefore be interpreted in the following way: $`(a)`$ the system forming network described in Fig. 1 posseses some features (mentioned above) preventing the use of Shannon information entropy and $`(b)`$ data represented by $`P(k)`$ can be quite adequately described in terms of only two parameters: the Tsallis entropy parameter $`q`$ and the mean number of links $`k_q=\lambda _0`$, i.e., their informational capacity is rather limited.

The question now is: what is the physical meaning of the $`q`$ parameter in the context of stochastic networks? There is a long list of possibilities in what concerns of the origin of $`q1`$ to be found in the literature dealing with nonextensivity . Out of these we shall only mention two: fluctuations and correlations. In  it was demonstrated that $`q`$ reflects fluctuations of the parameter $`\lambda _0`$ in exponential distribution (1) above. In fact, it turns out that $`(q1)/q=\pm Var(1/\lambda )/1/\lambda ^2`$. As is known from other branches of physics where Tsallis statistics is applied , the appearance of $`q`$ can also be caused by some correlations existing in the system under consideration. Such correlations (resulting, for example, from preferential attachements and ”rich-get-richer” phenomenon ) seem to play a decisive role in the description of stochastic networks. Therefore when choosing vertices with connectivity $`k`$, to which a new vertex is going to be connected, we shall assume that it will do so with probability that depends on the connectivity $`k`$. To illustrate this point let us introduce in the evolution equation

$$\frac{dP(k)}{dk}=\frac{1}{\lambda (k)}P(k),$$
(13)
parameter $`\lambda =\lambda (k)`$ given by a simple linear function of $`k`$:

$$\lambda (k)=\frac{\left[\lambda _0+(q1)k\right]}{q}.$$
(14)
It is easy to see that in this case one gets immediately $`P_q(k)`$ in the form of eq. (10). Notice that for $`q1`$ (i.e., for $`\lambda \lambda _0`$) one recovers the exponential distribution (7) 555Eq. (13) can be derived by dividing master equation, $`P(k)/t=cP(k)`$, by the ”growth of network” $`k/t`$. One gets then evolution equation, $`P(k)/k=cP(k)t/k`$, which for the linear dependence of the grow of network assumed here, $`k/t=a+bk`$, leads to eq. (13) with $`\lambda (k)=(a+bk)/c`$. Our example considered here corresponds to the choice: $`c=1/\lambda _0`$, $`a=1/q`$ and $`b=(11/q)/\lambda _0`$. Notice that $`1/q`$ plays role of weight with which we select the constant and linear terms in the equation describing the growth of network. Notice also that this kind of the growth of network, i.e., its dependence on $`k`$ (cf. ) corresponds to selective probability $`\mathrm{\Pi }(k)`$, which leads to power-like distribution (4)..

It must be stressed, however, that the information theory approach leads in a natural way (via maximalization of the respective information entropy) only to equilibrium (or stationary) distributions $`P_q(k)`$, whereas in models describing evolving complex networks  the functional form of $`P(k)`$ is determined by the growth equation $`k/t`$. Introducing more complicated network evolution than the one presented above when deriving eq.(4), for example allowing for the occurence of local events in the form of internal edges and rewirings, one gets (see eq. (111) of )

$$P(k)\left[k+\kappa (p,r,m)\right]^{\gamma (p,r,m)},$$
(15)
where $`p`$ is the probability that one is adding $`m`$ new edges to the system, $`r`$ is probability that one is rewiring $`m`$ edges and $`1pr`$ is probabibility that one is adding a new node to the system. The $`\kappa `$ and $`\gamma `$ in eq.(15) are given by :

$$\kappa (p,r,m)=A(p,r,m)+1\mathrm{and}\gamma (p,r,m)=B(p,r,m)+1,$$
(16)
where, in turn,

$`A(p,r,m)`$ $`=`$ $`(pr)\left[{\displaystyle \frac{2m(1r)}{1pr}}+1\right],`$
$`B(p,r,m)`$ $`=`$ $`{\displaystyle \frac{2m(1r)+1pr}{m}}.`$ (17)
One can now write formally eq.(15) in the form resembling eq.(10), i.e., as

$$P(k)\left[\kappa (p,r,m)\right]^{\gamma (p,r,m)}\left[\mathrm{\hspace{0.17em}1}(1q)\frac{k}{(1q)\kappa (p,r,m)}\right]^{\frac{q}{1q}}$$
(18)
and identify $`\gamma (p,r,m)=q/(q1)`$ and $`\lambda _0=(1q)\kappa (p,r,m)`$, or

$$q=1+\frac{1}{B}\mathrm{and}\lambda _0=\frac{A+1}{B}.$$
(19)
However, it must be noticed that, at least in the example considered here, the limit $`q=1`$ cannot be achieved because the quantity $`B`$ above is finite (for all reasonable values of parameters ). It means then that formula (10) is more general and captures (by means of parameters $`q`$ and $`\lambda _0`$) some additional feature of complex networks, not present in its simple formulation (as, for example, given by eq.(15)).

To summarize: We have demonstrated that, in order to apply the information theory approach to analysis of stochastic networks one has to use the nonextensive Tsallis information entropy $`S_q`$  leading to distribution $`P_q(k)`$ as given by eq.(10). As shown in Fig. 1, such distribution provides satisfactory description of data on number of links in random networks in the whole range of variable $`k`$ by means of only two parameters: the mean value of $`k`$ and the parameter characterizing the type of information entropy to be chosen, $`q`$. In this way one describes such disparate situations as the exponential model ER  (for $`q=1`$) and the scale-free, power-like models  (with $`q=\gamma /(\gamma 1)`$). For the value of nonextensivity parametr $`q=3/2`$, for which variance of our system is divergent, one obtains the exponent $`\gamma =3`$, which seems to be limited value observed in analyses of diverse systems displaying complex topology. Although only one example has been shown here in Fig. 1, it is obvious that one can just as easily also fit other, similar results discussed in the literature (cf., for example, ) 666For the most recent application of Tsallis statistics to investigation of Internet traffic problems see , for the previous one see .. The other point is the possible systematics of the $`q`$ parameter emerging from such a search, but this problem is outside of the scope of our presentation.

We conclude by saying that from the information theory point of view eq. (10) could be used to fit different data providing a pair of numbers $`(q,\lambda _0)`$ for each example. All competing models could be then checked for their ability to correctly reproduce these $`(q,\lambda _0)`$ and all models reproducing them correctly should be regarded as equally good from the point of view of distribution $`P(k)`$ because, according to the philosophy of infromation theory approach, they apparently contain the same amount of information existing in data which have been used. To distinguish between them further one would have to use some additional information contained in other network measures like, for example, clustering coefficient, distance between nodes, cycles or graph spectra.


text:
# Energy Landscapes and Rare Events

## 1. Introduction

Many problems in biology, chemistry and material science can be formulated as the study of the energy or free energy landscape of the underlying system. Well-known examples of such problems include the conformational changes of macromolecules, chemical reactions and nucleation in condensed systems. Very often the dimension of the state space is very large, and the energy landscape exhibits a hierarchy of structures and scales. These problems are becoming a major challenge in their respective scientific disciplines and are beginning to receive attention from the mathematics community. In this article, we report recent work in this direction. For a detailed account, we refer to .

We begin with a simple example. Plotted in Figure 1 is the solution of the stochastic differential equation

$$dx(t)=_xV(x(t))dt+\sqrt{\epsilon }dW(t)$$
(1.1)
where the potential

$$V(x)=\frac{1}{4}(1x^2)^2$$
(1.2)
and $`dW(t)`$ is Gaussian white noise, $`\epsilon =0.06,x(0)=1`$. Without the random perturbation, the solution would be $`x(t)x(0)=1`$. Indeed the deterministic part of dynamics in (1.1) does nothing but taking the system to local equilibrium states. With the random perturbation, the solution, over long time, exhibits completely different behavior. It fluctuates around the two local minima of $`V,x=1`$ and 1, with sudden transitions between these two states. The time scale of the transition, $`t_M`$ is much larger than the time scale of the fluctuation around the local minima, $`t_R`$. For this reason, we refer to $`x=1`$ and 1 as the metastable states.

Figure 1. Time series of the solution to the stochastic differential equation (1.1), with $`\epsilon =0.06`$.

Obviously the transition between the metastable states is of more interest than the local fluctuation around them. The transition time is much larger since it requires the system to overcome the energy barrier between the two states. This is only possible because of the noise. When $`\epsilon `$ is small, a huge noise term is required to accomplish this. For this reason, such events are very rare, and this is the origin of the disparity between the time scales $`t_M`$ and $`t_R`$.

This simple example illustrates one of the major difficulties in modeling such systems, namely the disparity of the time scales. It does not, however, illustrate the other major difficulty, namely, the large dimension of the state space and the complexity of the energy landscapes. Indeed for typical systems of interests the energy landscape can be very complex. There can be a huge number of local minima in the state space. The usual concept of hopping over barriers via saddle points may not apply (see ).

In applications, the noise comes typically from thermal noise. In this case, we should note that even though the potential energy landscapes might be rough and contain small scale features, the system itself experiences a much smoother landscape, the free energy landscape, since some of the small scale features on the potential energy landscape are smoothed out by the thermal noise.

Our objective in modeling such systems are the following:

1. Find the transition mechanism between the metastable states.
2. Find the transition rates.
3. Reduce the original dynamics to the dynamics of a Markov chain on the metastable states.

Our discussion will be centered around the following model problems:

$$\gamma \dot{x}(t)=V(x(t))+\sqrt{\epsilon }\dot{W}(t)$$
(1.3)
or

$$m\ddot{x}(t)+\gamma \dot{x}(t)=V(x(t))+\sqrt{m\epsilon }\dot{W}(t)$$
(1.4)
$`\epsilon `$ is related to the temperature of the system by $`\epsilon =2\gamma k_BT`$ where $`k_B`$ is the Boltzmann constant. We refer to (1.3) as type-I gradient flow and (1.4) as type-II gradient flow.

Before proceeding further, let us remark that there is a very well-developed theory, the large-deviation theory, or the Wentzell-Freidlin theory , that deals precisely with questions of the type that we discussed above. However as was explained in , this theory is not best suited for numerical purpose. Therefore we will seek an alternative theoretical framework that is more useful for numerical computations.

## 2. Transition state theory

Transition state theory (TST)  has been the classical framework for addressing the questions we are interested in. It assumes the existence and explicit knowledge of a reaction coordinate, denoted by $`q`$, that connects the two metastable states. In addition it assumes that along the reaction coordinate there exists a well-defined transition state, which is typically the saddle point configuration, say at $`q=0`$, and the two regions $`\{q<0\}`$ and $`\{q>0\}`$ defines the two metastable regions $`A`$ and $`B`$. For these reasons, transition state theory is restricted to cases when the system is simple and the energy landscape is smooth, i.e. the energy barriers are larger than the thermal energy $`k_BT`$.

Knowing the transition state, TST calculates the transition rates by placing particles at the transition state, and measuring the flux that goes into the two regions. For example, the transition rate from $`A`$ to $`B`$ is given approximately by

$$k_{AB}=\frac{1}{Z_0}\dot{q}(t)\delta _\mathrm{\Gamma }(q(t))\theta (q(t))𝑑\mu _A(q(0))$$
(2.1)
where

$$Z_0=𝑑\mu _A(q(0))$$
(2.2)
Here $`\delta _\mathrm{\Gamma }`$ is the surface delta function at $`q=0`$, $`\theta `$ is the Heaviside function, $`\mu _A`$ is the Lebesgue measure restricted to $`A`$. For a system with a single particle of mass $`m`$ and potential $`V`$, this gives

$$k_{AB}=\frac{\omega _0}{2\pi }e^{\frac{\delta E}{k_BT}}$$
(2.3)
where $`\delta E`$ is the energy barrier at the transition state, $`\omega _0=\left(\frac{V^{\prime \prime }(x_A)}{m}\right)^{\frac{1}{2}}`$, $`x_A`$ denotes the location of the local minimum inside $`A`$. Formulas such as (2.3) are the origin of the Arrhenius law for chemical reaction rates and Boltzmann factor for hopping rates in kinetic Monte Carlo models.

## 3. Reduction to Markov chains on graphs

For simplicity, we will discuss mainly type-I gradient flows (1.3). The Fokker-Planck equation can be expressed as

$$\frac{p}{t}(x,t)=\left(p_s(x)\left(\frac{p(x,t)}{p_s(x)}\right)\right)$$
(3.1)
where $`p_s`$ is the equilibrium distribution

$$p_s(x)=\frac{1}{Z}e^{\frac{V(x)}{k_BT}}$$

$`Z`$ is the normalization constant $`Z=_{R^n}e^{\frac{V(x)}{k_BT}}𝑑x`$. The states of the Markov chain consist of the sets $`\{B_j\}_{j=1}^J`$, where $`\{B_j\}_{j=1}^J`$ satisfies

1. The $`B_j`$’s are mutually disjoint
2. $$_Bp_s(x)𝑑x=1+o(k_BT)$$
(3.2)

where $`B=_{j=1}^JB_j`$. An illustration of the collection the $`\{B_j\}_{j=1}^J`$ is given in Figure 2. $`\{B_j\}`$ depends on $`T`$. As $`T`$ decreases, the $`B_j`$’s exhibits a hierarchical structure.

Figure 2. Illustration of the collection of metastable sets $`\{B_j\}`$ at different energies.

Having defined the states of the Markov chain, we next compute the transition rates between neighboring states. Denote by $`A`$ and $`B`$ two such neighboring states. We would like to compute the transition rate from $`A`$ to $`B`$. Without loss of generality, we may assume $`J=2`$. Let $`B_1,B_2`$ be the metastable region containing $`A`$ and $`B`$ respectively, and let $`n_j(t)=_{B_j}p(x,t)𝑑x`$, $`N_j=_{B_j}p_s(x)𝑑x`$, $`j=1,2`$. Applying Laplace’s method to (3.1) we get

$$\dot{n}_j(t)=\frac{\epsilon }{\kappa }\left(\frac{n_2(t)}{N_2}\frac{n_1(t)}{N_1}\right)+\text{higher order terms}$$
(3.3)
where

$$\kappa =_0^1𝑑\alpha \left(_{S^0(\alpha )}p_s(x)𝑑x\right)^1|\phi ^0(\alpha )|$$
(3.4)
$`\{\phi ^0(\alpha ),0\alpha 1\}`$ is a so-called minimal energy path, to be defined below, $`\{S^0(\alpha )\}`$ is the family of hyperplanes normal to $`\phi ^0`$.

The minimal energy path (MEP) is defined as follows. If $`V`$ is smooth, then $`\phi ^0`$ is a MEP if

$$(V)^{}(\phi ^0(\alpha ))=0$$
(3.5)
for all $`\alpha [0,1]`$, i.e. $`V`$ restricted to $`\phi ^0`$ is parallel to $`\phi ^0`$. In general there is not a unique $`\phi ^0`$ that is particularly significant, but rather a collection (a tube or several tubes) of paths contribute to the transition rates. However, one can define a MEP self-consistently via the equation

$$\phi ^0(\alpha )=\frac{1}{Z(\alpha )}_{S^0(\alpha )}xe^{\frac{V(x)}{k_BT}}𝑑x$$
(3.6)
where $`Z(\alpha )=_{S^0(\alpha )}e^{\frac{V(x)}{k_BT}}𝑑x`$. In the case when $`V`$ has two scales: $`V=\overline{V}+\delta V`$, $`|\delta V|O(k_BT)`$, and $`\overline{V}`$ is smooth, then $`\phi ^0`$ can be defined as the MEP of $`\overline{V}`$.

For type-I gradient systems, if the two metastable sets are separated by a single saddle point, then the MEP is the unstable manifold associated with the saddle point. MEP for type-II gradient systems is less trivial. In this case (3.3)-(3.6) have to be modified . Consider the simple example

$$\{\begin{array}{cc}\dot{q}\hfill & =p\hfill \\ \dot{p}\hfill & =\frac{V}{q}(q)p+\sqrt{\epsilon }\dot{W}\hfill \end{array}$$

where $`V(q)=\frac{1}{4}(1q^2)^2`$. It has two local equilibrium states $`A=(1,0)`$ and $`B=(1,0)`$. The MEP that connects $`A`$ and $`B`$ is plotted in Figure 3. It is not a smooth curve. The velocity is reversed at the saddle point. To verify that the MEP does reflect the true behavior of the transition path, we also plot the transition path obtained from direct simulation of the stochastic differential equation.

Figure 3. MEP for type-II gradient systems. The red line is the MEP in phase space, the green line is the transition path computed from solving the stochastic differential equation.

MEP is a very important concept since it defines the “most probable” transition path from which transition rates can be computed via equation (3.3). However we should emphasize that from a numerical point of view our task is not that of a conventional optimization or control problem, since there is not an objective function that we can easily work with. Instead our aim is to perform importance sampling in path space to sample the paths that contribute significantly to the switching.

Finally, if there exists a MEP that connects two metastable sets without going through a third one, then we connect these two metastable sets by a link. In this way, we form a graph. The original dynamics is then reduced to a Markov chain on this graph.

## 4. Previous numerical techniques

A variety of numerical techniques have been developed, most prominently in chemistry, but also in biology and material science, for computing MEPs and sometimes transition rates. Among the most well-known techniques in the chemistry literature are the nudged elastic band method and the transition path sampling technique. The nudged elastic band method (NEB)  aims at computing the MEP defined by (3.5). It represents the MEP by a discrete chain of states. These states evolve according to the potential forces of the system. To prevent the states from all falling to the two local equilibrium states, a spring force is applied to neighboring states to penalize the non-uniformity in the distribution of the states along the chain. This by itself may cause convergence to a path which is not a MEP. Hence a nudging technique is used, namely only the normal component of the potential force and the tangential component of the spring force is applied. NEB is a very effective method for small systems with relatively smooth energy landscapes. It has two main drawbacks. One is that it is highly inefficient and may not even be applicable to systems with rough energy landscapes. The other is the choice of the elastic constant. A large elastic constant requires a small time step in the evolution of the states. A small elastic constant will not achieve the desired uniformity of the states and hence will not give the required accuracy for the energy barrier.

A second important technique is transition path sampling (TPS) . This method aims at complex systems with rough energy landscapes by developing a Monte Carlo technique that samples the path space. Its efficiency hinges on the ability to produce new accepted paths from old ones.

Other techniques include the ridge method, blue mooth sampling, etc. . Often these methods require knowing beforehand the reaction coordinate.

Elber et. al propose to minimize the Onsager-Machlup action as a way of finding the most probable path for macromolecular systems . The Onsager-Machlup action is the same as the Wentzell-Freidlin action. From a numerical point of view, there are certain difficulties associated with minimizing this action functional. These issues are discussed in .

## 5. The string method

The basic idea in the string method, developed in , is to represent transition paths by their intrinsic parameterization in order to efficiently evolve and sample paths in path space. It has two versions. The zero temperature version is designed for smooth energy landscapes. The finite temperature version is designed for rough energy landscapes in which case thermal noise acts to smooth out the small scale features.

The simplest example of a zero-temperature string method is to evolve curves in path space by the gradient flow

$$\phi _t(\alpha ,t)=(V)^{}(\phi (\alpha ,t))+r(\alpha ,t)\widehat{\tau }(\alpha ,t)$$
(5.1)
Here $`\widehat{\tau }`$ is the tangent vector of the curve $`\{\phi (,t)\}`$, $`(V)^{}(\phi )`$ denotes the component of $`V`$ normal to $`\widehat{\tau }`$, $`r`$ is the Lagrange multiplier that enforces certain specific parameterization of the curves. For example if we require equal arclength parameterization, then we need $`\frac{}{\alpha }|\phi _\alpha |=0`$, i.e.

$$r(\alpha ,t)=\alpha _0^1V(\phi (\alpha ^{},t))\widehat{\tau }(\alpha ^{},t)𝑑\alpha ^{}_0^\alpha V(\phi (\alpha ^{},t))\widehat{\tau }(\alpha ^{},t)𝑑\alpha ^{}$$
(5.2)
We call such curves with intrinsic parameterization strings.

In practice the strings are discretized into a collection of points. These points move according to the normal component of the potential force. A reparameterization step is applied once in a while to enforce the proper parameterization of the strings.

The finite temperature string method is designed for systems with rough energy landscapes, particularly the case when the potential can be expressed in the form

$$V(x)=\overline{V}(x)+\delta V(x)$$

where $`\overline{V}`$ is smooth and $`|\delta V|O(k_BT)`$. In this case we would like to compute the MEP of $`\overline{V}`$ without first computing $`\overline{V}`$ explicitly. This is achieved by creating an ensemble of a special type. Our computational object will be a string connecting the two metastable sets, together with a family of probability measures on the hyperplanes normal to the string. Consider the stochastic equation

$$\phi _t^\omega (\alpha ,t)=𝐏_\alpha ^0(V(\phi ^\omega (\alpha ,t)))+r^0(\alpha ,t)\widehat{\tau }^0(\alpha ,t)+𝐏_\alpha ^0\eta ^\omega (\alpha ,t)$$
(5.3)
where $`\eta ^\omega `$ is Gaussian noise with mean 0 and correlation

$$𝐄\eta ^\omega (\alpha ,t)\eta ^\omega (\alpha ^{},\tau )=\{\begin{array}{cc}2k_BT\delta (t\tau ),\hfill & \text{if }\alpha =\alpha ^{}\hfill \\ 0,\hfill & \text{if }\alpha \alpha ^{}\hfill \end{array}$$

The projection operator $`𝐏_\alpha ^0`$ is defined by projecting to the hyperplane normal to the string $`\{\phi ^0(\alpha ,t),0\alpha 1\}`$, where

$$\phi ^0(\alpha ,t)=𝐄\phi ^\omega (\alpha ,t)$$

$`\widehat{\tau }^0(\alpha ,t)`$ is the tangent vector of $`\phi ^0`$ at $`\alpha `$, $`r^0`$ is the Lagrange multiplier that enforces proper parameterization of $`\phi ^0`$.

Theorem 1.1.

1. The statistical steady state of (5.3) satisfies:

$$\phi ^0(\alpha )=\frac{1}{Z(\alpha )}_{S^0(\alpha )}xe^{\frac{V(x)}{k_BT}}𝑑x$$
(5.4)
where $`S^0(\alpha )`$ is the hyperplane normal to $`\phi ^0`$ at $`\alpha `$, $`Z(\alpha )=_{S^0(\alpha )}e^{\frac{V(x)}{k_BT}}𝑑x`$.
2. The stationary distribution of (5.3) is given by the family of distributions

$$\mu _\alpha (x)=\frac{1}{Z(\alpha )}e^{\frac{V(x)}{k_BT}}\delta _{S^0(\alpha )}(x)$$
(5.5)

Knowing $`\{\phi ^0\}`$ and $`\{\mu _\alpha \}`$, the transition rates and free energy landscapes can be computed, see .

The finite temperature string method is applied to the perturbed Mueller potential

$$V(x)=V_m(x)+\delta V(x)$$

where $`V_m(x)`$ is the so-called Mueller potential (see ), $`\delta V`$ is a random perturbation. The results of the (finite temperature) string method is shown in Figure 4. Also plotted is the MEP of $`V_m`$ (since $`\overline{V}=V_m`$ is explicitly known for this particular example) as well as the fluctuations around $`\phi ^0`$.

Figure 4. Effective MEP and local fluctuations for the perturbed mueller potential. The red curve is the MEP for $`\overline{V}=V_m`$. The black curve is the MEP computed from the finite temperature string method. The green curves show the size of the fluctuations.

## 6. Concluding remarks

There are several important topics that we did not cover in this brief report. These include the effect of dynamics, non-gradient systems, and acceleration techniques. These are discussed in . Also found in these references are applications of the ideas discussed here to thermal activated reversal of magnetic thin films, models of martensitic transformations, and the formation of $`C_{60}`$ from 60 carbon atoms. The last example is a case when the barrier is entropic. Even though the potential energy is mainly going downhill, the free energy has barriers because of entropic effects. Such examples are found frequently in biopolymers.

From a numerical point of view, our main idea for overcoming the difficulty caused by the disparity of the times scales is to reformulate the problem as a boundary value problem instead of initial value problem, since we have some knowledge of the initial and final state of the system. Compared with other existing methods that assume explicit knowledge of a reaction coordinate, our method finds the reaction coordinate self-consistently during the computation.

The topic discussed here is relatively new in applied mathematics, but it is of paramount importance in science and particularly computational science. Progress in this area will likely have a fundamental impact in many areas of applications.

Acknowledgment. We are grateful to Bob Kohn for many stimulating discussions. This work is partially supported by NSF via grant DMS01-30107.