Methods for the study of rare events

Contents

Metastability

We consider the overdamped Langevin equation

dx = - grad V(x) dt + (2kBT)(1/2)dw

where V is a potential, kB is Boltzmann's constant, T is the temperature, and w is the n-dimensional Brownian motion. This equation is largely used in molecular dynamics and material science as a simplified model. On the one hand, it is convenient for theoretical study and numerical experiments. On the other hand, it allows to explain many phenomena.

If the temperature is low, such a system spends the most of the time in the neighborhoods of the potential minima and performs rare transitions between these neighborhoods. We are interested in finding the most likely transition paths for such a system.

Minimum Energy Paths

As the temperature approaches zero, the most likely transitions between two given minima follow so-called Minimum Energy Paths or MEP's, i.e., paths, either parallel or antiparallel to the gradient of the potential.
However, if the potential has Morse index 2 or higher saddles, there are families of MEP's continuously transformed one into another. The shaded region in the Figure on the right consists of MEP's. Each of these MEP's have the same probability according to the Wentzel-Freidlin theory of Large Deviation. Therefore, the transition paths are not quite predictable in such a situation.

The String Method as a Dynamical System

A number of methods for finding MEP's have been developed. Two of them are
the String Method (W. E, W. Ren, and E. Vanden-Eijnden, 2002, 2007) and
the Nudged Elastic Band method (H. Jonsson, G. Mills, and K.W. Jacobsen, 1998).
These methods are different numerical procedures, but they correspond to the same curve evolution in the continuous-time continuous-space setting.

It has been assumed that any curve evolved according to the string method equation converges to an MEP. We have realized that it is not so. The large-time behavior of a curve evolved according to the string method equation can be much more complicated. The curve evolution according to the string method equation

&phit = - grad V(&phi) + &lambda&tau

(&tau is the unit tangent vector to the curve) is equivalent to its evolution according to the gradient flow

&phit = - grad V(&phi).

One can see this from the fact that the normal velocity is the same for these two equations.

The gradient flow is more convenient for the analysis in the continuous-time continuous-space setting.
We have shown that the &omega-limit set of a curve evolving according to the gradient flow consists of critical points and heteroclinic trajectories, i.e., MEP's. If the &omega-limit set of a curve is a curve, i.e., a single MEP, the curve necessarily converges to it.
However, it is not necessarily a single MEP if the potential has a Morse index 2 or higher saddle. It can be of any dimension up to the dimension of the space. In this case, the curve can either endlessly wander around its limit set without convergence, or it can converge to it, i.e., fill it. The movies below demonstrate these phenomena.
Movie: A curve filling a 2D region.
Its &omega-limit set is the 2D region.
The curve converge to it. 		Movie: A revolving curve in 3D. Its &omega-limit set is the surface of revolution. No convergence. Movie: The initial curve does not pass through the Morse index 2 saddle. Its &omega-limit set is 2D and reminiscent of an umbrella. No convergence.

Conclusion

Curve evolution according to the string method equation can account for physical ambiguity.
The details and the exact criterion of when a curve does not converge to a single MEP are in the reference.

References

  • The String Method as a Dynamical System,
          Maria Cameron, Robert Kohn, and Eric Vanden-Eijnden, Journal of Nonlinear Science, in revision

The MaxFlux functional

The MaxFlux functional is given by

E(&phi) = &int 0Le &beta V(&phi(s))ds, where &beta &equiv (kBT)-1.

Its minimizer is the path along which the reactive flux is maximal at a given finite temperature. The MaxFlux path is always smooth; it "flies through" the unimportant local minima but passes close to the highest saddle, the one which is crucial for the estimation of the transition rate.

The MaxFlux functional has been around in the chemical physics community for almost 30 years. However, it has not been widely used. The reason is that in the case of most interest where the temperature is low, the dimension of the space is large, and the potential is stiff, its minimization is a difficult problem. In the limit T &rarr 0 the functional becomes singular. Moreover, its original derivation (Berkowitz, 1983) is heuristic and hard to justify.

This work consists of three parts.

  1. We have derived the MaxFlux functional in the framework of the Transition Path Theory (Weinan E and Eric Vanden-Eijnden) using the lower bound approach. Further, we have shown that the Euler-Lagrange equation for the this functional can be obtained by a short geometric argument. Both of these arguments are valid at any finite temperature. However, the MaxFlux functional can be used for estimating of the transition rate only if the temperature is low. We have obtained an estimate of the transition rate (the exponent and the prefactor) consistent with the classic result of Kramers and the estimate by the large deviation theory.
  2. We have developed an efficient way to find a minimizer of this functional. We solve a discretized and preconditioned version of the Euler-Lagrange equation for a uniformly parametrized path using a combination of a quasi-Newton Broyden's method (see "Numerical Optimization" by Nocedal & Wright) and a temperature-relaxation strategy. We have shown that the resulting path is necessarily uniformly parametrized, and we have to do neither reparametrization nor recomputation of the second derivative operator during the Broyden's iterations. This feature may give an advantage to our MaxFlux method over the string method. If the potential is stiff, dimension is high, and the temperature is low, reparametrization destroys convergence.
  3. We have applied our computational procedure to the problem of rearrangement of the Lennard-Jones cluster of 38 atoms (LJ38). We have chosen this problem because
    • on one hand, it is difficult due to high dimensionality (38 &sdot 3 &sdot # of points in the path) and stiffness of the Lennard-Jones potential;
    • on the other hand, its potential landscape is well-studied by D. Wales et al.
    The two lowest minima of (LJ38) are:
    Face-centered-cubic truncated octahedron with the point group Oh
    (the lowest minimum).     		Incomplete icosahedron with the point group C5v
    (the second lowest minimum).

Application to LJ38

One of the problems arising here is the choice of the endpoints for the transition path. (There are totally 38! &asymp 5.2 &sdot 1044 possible pairs of endpoints.) We have found a collection of reasonable endpoints as follows. We define an artificial pairwise potential of interaction between the atoms in the original configuration and the positions at target configuration. This potential needs to be

  • non-stiff,
  • prohibiting the atoms from penetrating each other, and
  • its minima are only achieved if all of the atoms have achieved different positions in the target configuration.
In result, the flow according to this potential is easy-to-run and leads to some reasonable correspondence between the original and the target configuration. One can readily imagine this procedure in 2D. Suppose we have a collection of sponge balls on a horizontal surface with a collection of smooth sinks such that each sink can hold exactly one ball. The centers of the sinks are the centers of the atoms in the target configuration of atoms. At time zero the arrangement of the balls corresponds to the original configuration of atoms. We let the sponge balls roll. This process stops if each of the balls falls into one of the sinks. Thus, we can find a correspondence between the balls and the sinks, i.e., the atoms in the original and the target configurations.

Changing the artificial potential we have found several different pairs of endpoints. The potential energy along the transition paths for a sequence of temperatures for two pairs of endpoints are shown in Figures below.
Movie: Click on the image to watch a movie. Movie: Click on the image to watch a movie.

The MaxFlux paths are smooth and do not pass exactly through the critical points. In order to find the highest saddles for the transitions shown above we used the string method. The found MEP's corresponding the the MaxFlux paths at kBT = 6 are shown below.
Movie: Click on the snapshots to watch a movie.
Movie: Click on the snapshots to watch a movie.

References

  • Revisiting the MaxFlux functional with application to Lennard-Jones-38,
          Maria Cameron and Eric Vanden-Eijnden, J. Chem. Phys., about to be submitted
  • The MaxFlux functional: Derivation, Numerics, and Application to 38-Lennard-Jones cluster,
          Maria Cameron, Robert Kohn, and Eric Vanden-Eijnden, J. Comp. Phys., in preparation