notebook.community

Edit and run

Definition: "The gentlest ascent dynamics (Weinan and Zhou, 2011)

Application: "Atomistic simulations of rare events using gentlest ascent dynamics" (Samanta and Weinan, 2011)

Background:

Transition states correspond to saddle points on the potential energy surface
We can characterize critical points (minima, maxima, saddle points) by computing properties of the matrix of second-order partial derivatives (the "Hessian" matrix) at each critical point
- Eigenvalues of the Hessian tell us whether the surface is concave up, concave down, or a mixture of both
- Index of saddlepoint: the fraction of negative eigenvalues of the hessian:
  - All negative: concave "down"
  - All positive: concave "up"
  - Mixed: inconsistent concavity

Notes on Weinan and Zhou, 2011

Steepest descent dynamics: Given an energy function $V$ on $\newcommand{\reals}{\mathbb{R}} \reals^n$, steepest descent dynamics are: $$ \dot{ \newcommand{\x}{\mathbf{x}} \newcommand{\v}{\mathbf{v}} \x} = - \nabla V (\x)$$
- If $\x(\cdot)$ is a solution, then $V(\x(t))$ is a decreasing function of $t$
- Stable fixed points: local minima of $V$
- Basins of attraction (aka potential wells of $V$) are separated by "separatrices" on which dynamics converges to saddle points
Goal:
- Opposite: climb out of a basin of attraction
Strategy:
- Naïve suggestion: flip sign (from $- \nabla V (\x)$ to $ \nabla V (\x)$)
  - Outcome: find local maxima
- We need a dynamics that converges to index-1 saddle points of $V$, of interest to noise-induced transition between metastable states: $$ \begin{align} \dot{\x} & = -\nabla V(\x) + 2 \frac{(\nabla V, \v)}{(\v,\v)} \v\\ \dot{\v} & = -\nabla^2 V (\x) \v + \frac{(\v,\nabla^2 V \v)}{(\v,\v)} \v\\ \end{align}$$

Notes on Samanta and Weinan 2011

Complex system dynamics frequently involve rare transition events between metastable states

Next steps and generalization:

Can we automatically construct dynamical systems that have the desired properties? (e.g. of hopping out of basins of attraction: create a number of test potentials and see how far a walker gets on each one)
What is a good set of design criteria for metadynamics systems?

Implementation

(Based on MATLAB implementation: http://web.math.princeton.edu/string/gad/)



In [55]:

    
import numpy as np
import numpy.random as npr
import pylab as pl
%matplotlib inline



In [56]:

    
dt = 0.01
T = 100
ndim=2

X = np.zeros((T,ndim))
X[0] = [0.45,-0.45]
direction = -np.array([0.87,0.07])
normalize = lambda vec : vec / np.sqrt(sum(vec**2))
direction = normalize(direction)
V = np.zeros((T,ndim))



In [57]:

    
from numpy import pi, sin, cos



In [58]:

    
def force(x):
    return np.array((-pi*cos(pi*x[0])*sin(pi*x[1]),
                     -pi*sin(pi*x[0])*cos(pi*x[1])))



In [59]:

    
force(X[0])









    Out[59]:





array([ 0.48540276, -0.48540276])



In [60]:

    
def hessian(x):
    return pi**2 * np.array([[-sin(pi*x[0])*sin(pi*x[1]),
                              cos(pi*x[0])*cos(pi*x[1])],
                             [cos(pi*x[0])*cos(pi*x[1]),
                              -sin(pi*x[0])*sin(pi*x[1])]])



In [61]:

    
hessian(X[0])









    Out[61]:





array([[ 9.62807799,  0.24152641],
       [ 0.24152641,  9.62807799]])



In [68]:

    
for b in range(10):
    #direction = npr.randn(2)
    #direction = normalize(direction)
    for i in range(T-1):
        F = force(X[i])# + npr.randn(ndim)*0.01
        H = hessian(X[i]).T# + npr.randn(ndim,ndim)*0.01
        c1 = direction.dot(F)
        X[i+1] = X[i] + F*dt - 2*dt*c1.dot(direction)
        direction = normalize(direction - dt*H.dot(direction))
    pl.plot(X[:,0],X[:,1])



In [69]:

    
c1









    Out[69]:





array([ -6.38299509e-14,   6.38299509e-14])



In [70]:

    
direction









    Out[70]:





array([[-0.70710678,  0.70710678],
       [ 0.70710678, -0.70710678]])



In [41]:

    
normalize(direction)









    Out[41]:





array([ 0.5741445 ,  0.81875399])

MD-GAD



In [ ]: