Weinan E. Principles of Multiscale Modeling

Подождите немного. Документ загружается.

Chapter 10

Rare Events

We have already discussed various examples of problems with multiple time scales.

In this chapter, we will tu rn to another important class of such problems: Rare events.

These are the kind of events that happen rarely, compared with the relaxation time scale

of the system, but when they do happen, they happen rather quickly and they have

important consequences. Typically, a small amount of noise is present in the system,

and these rare events are driven by this small amount of noise. In order for the event to

happen, the system has to wait for the diﬀerent components of the noise to work together

to bring the system over some energy barrier, or to go through a sequence of correlated

changes. Chemical reactions, conformation changes of bio-molecules, nucleation events,

and in some cases extreme events that lead to material failure or system failure, are all

examples of rare events.

It should be emphasized that in a way, the rare events we are interested in are not

really rare. For example, conformation changes of biomolecules are rare on the time scale

of molecular vibration (which is typically on the order of femto-seconds, 10

−15

s), but

they are not rare on the time scale of our daily lives, which is often measured in minutes,

hours or days. After all, all biologoical processes are driven by such events.

An important concept is th at of metastability. Roughly speaking, metastable systems

are systems which are stable over a very long time scale compared with the relaxation

time scale of the system, but ultimately they undergo some noise-induced transition.

It is helpful to distinguish several diﬀerent situations:

1. Gradient and non-gradient systems. A gradient system has an underlying potential

energy and an associated energy landscape. The diﬀerent stable states are the

427

428 CHAPTER 10. RARE EVENTS

Figure 10.1: Time series from molecular dynamics simulation for the conformation

changes of alanie dipeptide. Shown here is the time series of a dihedral angle.

local minima of the energy. They are separated by some barriers or sequence of

barriers as well as other intermediate stable states. For this reason, transition events

between these locally stable states are also referred to as barrier-crossing events.

The d ynamics of the system is such that it spends most of its time vibrating around

diﬀerent local energy minima, with occasional transitions between the diﬀerent local

minima. For a given transition, we call the local minimum b efore the transition the

initial state and the local minimum after the transition the ﬁnal state or the end

state.

The mechanism of transition in n on-gradient systems can be very diﬀerent from

that of the gradient system. Here we will limit ourselves to gradient systems.

2. Energetic barriers and entropic barriers. Entropic barriers refer to the cases where

the energy landscape is rather ﬂat, bu t the system has to wonder through a very

large set of conﬁgurations in order to reach the ﬁnal state. Picture this: You are

in a very large room with no lights, and you are supposed to ﬁnd the small exit

which is the only way of getting out of the room.

3. Systems with simple or complex energy landscap es. This notion is only a relative

one. There are two important factors to consider. One is the length scale over which

the energy landscape varies. The second is the size of the barriers compared with

429

thermal noise. If the energy landscape is rather smooth and the individual barriers

are much larger than the typical amplitude of the noise (which is characterized

by k

T in physical systems), we refer to it as a system with simple or smooth

energy landscape, or simply a simple system. This is the case when the system has

a few large barriers. As we will discuss below, in t his case, the relevant objects

for the transition events are the saddle points which serve as the transition states,

and minimum energy p aths (MEP) which are the most probable paths. The other

extreme is the case when the landscape is rugged with many saddle points and

barriers, and most barriers are of comparable size to the thermal noise. We refer

to this class of systems as systems with complex or rough energy landscapes. In

this situation, the individual saddle points are sometimes not enough to act as the

obstacle for the transition. The real obstacle to the transition is caused by the

combined eﬀect of many individual barriers. There may not be a well-deﬁned most

probable path. Instead, many paths in the transition path ensemble contribute to

the transition.

Figure 10.2: An example of a rough energy landscape and one particular transition path

between the two metastable states.

Figure 10.2 illustrates the situation with rough energy landscapes in a temperature

regime where the thermal energy is comparable to the small scale energy ﬂuctuations in

430 CHAPTER 10. RARE EVENTS

the landscape. The two big blue regions are the metastable sets. One exp ect s that tran-

sition between the two regions should proceed by wondering in the light blue/yellowish

region, i.e. most transition paths between the two metastable sets are restricted to a thin

tube in the light blue/yellowish region. One thing that we will do in this chapter is to

make such statements more precise.

Below we will ﬁrst discuss the theoretical background for describing transition between

metastable states. We will then discuss numerical algorithms that have been developed

to compute the relevant objects, e.g. transition states, transition pathways and transition

tubes, for such transitions. These algorithms require knowing something about the ﬁnal

states of the transition. Finally, we discuss some examples of t echniques that have been

proposed for directly accelerating the d ynamics. These techniques do not require knowing

anything about the ﬁnal states.

The setup is as follows. We are given two metastable states of the system, which are

respectively th e initial and ﬁnal states of the transition. For a system with a smooth

energy landscape, these metastable states are t he local minima of the potential en ergy,

or sometimes the free energy. In chemistry, one typically refers the initial state as the

reactant state, and the ﬁnal state as the product state. For the same reason, transition

events are also referred to as reactions.

Next, we introduce the kind of dynamics that we will consider, following [18]. From

a physical viewpoint, it is of interest to study

1. The Hamiltonian dynamics, with no explicit noise







x = M

−1

p = −∇U

(10.0.1)

where M is the mass matrix. Here we have used U (instead of V ) to denote the

potential energy of the system. In this case, the role of noise is played by the

chaotic dynamics of the system. For example, integrable system does not belong

to the class of systems that we are interested in.

2. Langevin dynamics







x = M

−1

p = −∇U − γp +

√

2γk

T M

1/2

(10.0.2)

431

where γ is a friction coeﬃcient, which can also be a tensor although we will only

consider the case when it is a scalar;

w is the standard Gaussian white noise.

3. Overdamped dynamics:

x = −∇V +

2γk

w (10.0.3)

All three examples can be cast into the form:

z = −K∇V (z) +

√

2εσ

w(t) (10.0.4)

where K ∈ R

× R

and σ ∈ R

× R

are matrices related by

(K + K

) = σσ

= a (10.0.5)

For the Hamiltonian dynamics, we have z = (x, p)

, V (z) =

hp, M

−1

pi+ U(x), σ = 0,

and

K =

0 −I

I 0

(10.0.6)

For the Langevin dynamics, we have

K =

0 −I

I γM

σ =

0 0

0 (γM)

1/2

, (10.0.7)

For the overdamped dynamics, we h ave z = x, K = γ

−1

I , and σ = γ

−1/2

I. In all three

cases, we have ε = k

T . In the following, we will focus on (10.0.4). We will assume that

the potential energy function V is a Morse function, namely, V is twice diﬀerentiable and

has non-degenerate critical points

det H

6= 0 (10.0.8)

where H

= ∇∇V (z

) is the Hessian of V at the critical point z

. It is relatively

straightforward to check that the distribution given by

dµ

(z) = Z

−1

−V (z)/ε

dz where Z =

Ω

−V (z)/ε

dz (10.0.9)

is invariant under the dynamics deﬁned by (10.0.4). We will also assume that the dy-

namics is ergodic, i.e. t he invariant distribution is unique.

432 CHAPTER 10. RARE EVENTS

10.1 Theoretical background

10.1.1 Metastable states and reduction to Markov chains

In the following, we will speak about metastable states or metastable sets. These are

sets of conﬁgurations that are stable during a certain time scale, often some relaxation

time scale. If initiated in that set, the system will remain in it during the speciﬁed time

scale. This time scale can be the relaxation time scale of a potential well, or a set of

neighboring wells. The important point is that the notion of metastability is a relative

one, relative to certain reference time scales.

In a simple system, each local minimum of the potential is a metastable set, if the

reference time scale is taken to be the relaxation time scale of the potential well for that

local minimum. The potential well, or the basin of attraction is the set of points in

the conﬁguration space that are attracted to this particular local minimum under the

gradient ﬂow:

= −∇V (x) (10.1.1)

where V is the potential of the system. The diﬀerent basins of attraction are separated by

the separatrices. The separatrices themselves make up an invariant set for the gradient-

ﬂow dynamics (10.1.1). The stable equilibrium points on this invariant set are the saddle

points of V .

For complex systems, metastable sets are less well-deﬁned since we must deal with a

set of wells. But their intuitive meaning is quite clear. In practice, they can be deﬁned

with the help of some coarse-grained variables or order parameters. For example, for a

molecule, we can deﬁne the metastable sets by restricting the values of a set of torsion

angles.

More precise deﬁnition of the metastable sets can be found using the spectral theory

of the transfer operators or the generators associated with the underlying stochastic dy-

namics. See for example [42, 36]. Roughly speaking, metastability can be related to gaps

in the spectrum of the relevant generators of the Markov process and metastable sets can

be deﬁned through the eigenfunctions associated with the leading eigenvalues. Indeed,

it can be shown that these eigenfunctions are approximately piecewise constant, and the

subsets on which the eigenfunctions are approximately constant are the metastable sets.

Going back to the global picture, our objective is not to keep track of the detailed

dynamics of the system, but rather to capture statistically the sequence of hops or tran-

10.1. THEORETICAL BACKGROUND 433

sitions between diﬀerent local minima or metastable states. This means that eﬀectively,

the dynamics of the system is modeled by a Markov chain: the local minima or the

metastable states are the states of the chain and the hopping rates are the transition

rates between diﬀerent metastable states. This idea has been used in two diﬀerent ways:

1. It has been used as the basis for accelerating molecular dynamics. One example is

the hyperd ynamics to be discussed below. Its purpose is not to capture the detailed

dynamics but only the eﬀective dynamics of the Markov chain.

2. It can be used to design more accurate kinetic Monte Carlo schemes. In stan-

dard kinetic Monte Carlo methods, the possible states and the transition rates are

pre-determined, often empirically. This can be avoided by ﬁnding the states and

transition rates “on-the-ﬂy” from a more detailed model. S ee for example [60].

10.1.2 Transition state theory

The classical theory that describes barrier-crossing events in gradient systems is the

transition state theory (TST), developed by Eyring, Polanyi, Wigner, Kramers, et al

[31]. There are two main components in the transition state theory, the notion of transi-

tion states and an expression for the transition rate. Transition states are the dynamic

bottlenecks for the transition. They are important for the following reasons:

1. With very high probability, transition paths have to pass through a very small

neighborbood of the transition states. Exceptional paths have exponentially small

relative probability.

2. Once the transition state is passed, the system can relax to the new stable state

via its own intrinsic dynamics, on the relaxation time scale.

For simple systems, given two neighboring local minima, the transition state is nothing

but the saddle point that separates the two minima.

The transition rate, i.e. the average number of transitions per unit time interval, can

be computed from the probability ﬂux of particles that pass through the neighborhood

of the transition state. As an example, let us consider the canonical ensemble of a

Hamiltonian system. For simplicity, we will focus on an example of a one-dimensional

system [31]. Let V be a one-dimensional potential with one local minimum at x = −1 and

434 CHAPTER 10. RARE EVENTS

another at x = 1 separated by a saddle point at x = 0. We are interested in computing

the transition rate from the local minimum at x = −1 to the local minimum at x = 1 in

a canonical ensemble with temperature T . According to the transition state theory, the

transition rate is given by the ratio of two quantities: The ﬁrst is the ﬂux of particles

that tranverse the transition state region from the side of the negative real axis to the

side of the positive real axis. The second is the equilibrium distribution of particles on

the side of the negative real axis, which is the basin of attraction of the initial state at

x = −1. This second quantity is given by:

x<0

−βH(x,p)

dxdp (10.1.2)

where H(x, p) =

+ V (x) is the Hamiltonian of the system, Z =

−βH(x,p)

dxdp is

the normalization constant, β = (k

T )

−1

. The ﬁ rst quantity is given by

x=0

pθ(p)e

−βH(x,p)

dp (10.1.3)

where θ is the Heaviside function: θ(z) = 1 if z > 0 and θ(z) = 0 if z < 0. To leading

order, we have

∼

Zβ

−βV (0)

(10.1.4)

If we approximate V further by a quadratic function centered at x = −1, V (x) ∼

V (−1) + 1/2V

′′

(−1)(x + 1)

, we obtain:

∼

2π

′′

(−1)

−βV (−1)

(10.1.5)

The ratio of the two quantities gives the rate we are interested in:

2π

′′

(−1)

−βδE

(10.1.6)

where δE = V (0) − V (−1) is the energy barrier.

Note that this rate depends on the mass of the particle and the second derivative of

the potential at the initial state, but it does not depend on the second derivative of the

potential at the saddle point.

(10.1.6) can obviously be extended to high dimension [31].

TST has been extended in a number of directions. The most important generalization

is to replace the region near the saddle point by a dividing surface, see for example [63].

10.1. THEORETICAL BACKGROUND 435

In this case, one can still compute the probability ﬂux through the surface and use it

as an estimate for the transition rate. Obviously this is only an upper bound for the

transition rate, since it does not take into account the fact that particles can recross the

dividing su rface and return to the reactant state. The correct rate is only obtained when

the success rate, often called the transmission coeﬃcient, is taken into account [26]. This

is conceptually appealing. But in practice, it is often diﬃcult to obtain the transmission

coeﬃcient.

Some of more mathematically oriented issues on transition state theory are discussed

in [54].

10.1.3 Large deviation theory

Large deviation theory is a rigorous mathematical theory for characterizing rare events

in general Markov processes [57, 27]. The special case for random perturbations of

dynamical systems was considered by Freidlin and Wentzell [27]. This theory is often

referred to as the Wentzell-Freidlin theory.

Consider the stochastically perturbed dynamical system:

x = b(x) +

√

w, x ∈ R

(10.1.7)

Fix two points A and B in R

and a time parameter T > 0. Let φ : [−T, T ] → R

a continuous and diﬀerentiable path such th at φ(−T ) = A, φ(T ) = B. The Wentzell-

Freidlin theory asserts roughly that

Prob.{x stays in a neighborhood of φ over th e interval [−T, T ]} ∼ e

−S(φ)/ε

(10.1.8)

where the action functional S is deﬁned by:

(φ) =

−T

φ(t) − b(φ(t))|

dt (10.1.9)

The precise statement and its proof can be found in [27]. For our purpose, it is

more instructive to understand intuitively (and formally) the origin of this result. From

(10.1.7), we have

w = (

x − b(x))/

√

ε (10.1.10)

Roughly sp eaking, the probability density for the white noise is e

−

. This gives

(10.1.8) and (10.1.9).

436 CHAPTER 10. RARE EVENTS

With (10.1.8), ﬁnding the path with maximum probability is turned into the problem

of ﬁnding the path with minimum action:

inf

(φ) (10.1.11)

subject to the constraint that φ(−T ) = A, φ(T ) = B. Consider the case when b(x) =

−∇V (x). Assume that A and B are two neighb oring local minima of V separated by

the saddle point C. Then we have:

Lemma (Wentzell-Freidlin [27]):

1. The minimum of (10.1.11) is given by:

∗

= 2(V (C) − V (A)) (10.1.12)

2. Consider the paths

(s) = ∇V (φ

(s)), φ

(−∞) = A, φ

(∞) = C (10.1.13)

(s) = −∇V (φ

(s)), φ

(−∞) = C, φ

(∞) = B (10.1.14)

then

∗

= S

∞

(φ

) + S

∞

(φ

) = S

∞

(φ

) (10.1.15)

The most probable transition path is then the combination of φ

and φ

: φ

goes

against the original dynamics and therefore requires the action of the noise. φ

simply

follows the original dynamics and therefore does not require the help from the noise.

It is not diﬃcult to convince oneself that the minimum in T in (10.1.11) is attained

when T = ∞. To see why the minimization problem in (10.1.11) is solved by the path

deﬁned above, note that

∞

[φ

] = 2(V (C) − V (A)), S

∞

[φ

] = 0. (10.1.16)

In addition, for any path φ connecting A and a point on the separatrix that separates