Can any local geometric features of the potential predict kinetic features of a configuration?

Related work: Weinan E and colleagues identify the index-1 saddle points of the potential with transition states in their work on "gentlest ascent dynamics." They can then construct a simple dynamical system (given a potential function that is guaranteed to have the index-1 saddle points of the potential function as its set of invariant points.

does this work for more complicated potentials?
- alanine dipeptide, etc. -- the progression of test problems in MSMbuilder
are there other features of interest?
- distance to nearest saddlepoint / minimum?
how can we recover kinetic features if we use altered dynamics?
- transition path sampling: use this new dynamics to find the metastable states efficiently, then compute determine their transition rates separately

Concrete steps:

Plot saddlepoint index of $U(x)$ vs. kinetic properties of $x$ (e.g. stationarity, committor function)
- committor function!



In [1]:

    
from __future__ import print_function
import numpy as np
from matplotlib import pyplot as plt
from matplotlib.colors import LogNorm

from msmbuilder.decomposition import tICA
from msmbuilder.example_datasets import fetch_met_enkephalin
from msmbuilder.featurizer import AtomPairsFeaturizer
from sklearn.pipeline import Pipeline

%matplotlib inline



In [2]:

    
dataset = fetch_met_enkephalin()
print(dataset.DESCR)









    



The dataset consists of ten ~50 ns molecular dynamics (MD) simulation
trajectories of the 5 residue Met-enkaphalin peptide. The aggregate
sampling is 499.58 ns. Simulations were performed starting from the 1st
model in the 1PLX PDB file, solvated with 832 TIP3P water molecules using
OpenMM 6.0. The coordinates (protein only -- the water was stripped)
are saved every 5 picoseconds. Each of the ten trajectories is roughly
50 ns long and contains about 10,000 snapshots.

Forcefield: amber99sb-ildn; water: tip3p; nonbonded method: PME; cutoffs:
1nm; bonds to hydrogen were constrained; integrator: langevin dynamics;
temperature: 300K; friction coefficient: 1.0/ps; pressure control: Monte
Carlo barostat (interval of 25 steps); timestep 2 fs.

The dataset is available on figshare at

http://dx.doi.org/10.6084/m9.figshare.1026324



In [3]:

    
from itertools import combinations
heavy_atoms = dataset.trajectories[0].topology.select_atom_indices('heavy')
heavy_pairs = list(combinations(heavy_atoms, 2))



In [4]:

    
traj = dataset.trajectories



In [5]:

    
feat = AtomPairsFeaturizer(heavy_pairs)
featurized = feat.fit_transform(traj)



In [6]:

    
feat_stack = np.vstack(featurized)



In [7]:

    
import msmbuilder



In [ ]:

    
msmbuilder.featurizer.



In [15]:

    
from simtk import openmm
from simtk.openmm import app



In [16]:

    
b = openmm.AmoebaBondForce()



In [18]:

    
forcefield = app.ForceField('amber99sb.xml', 'tip3p.xml')



In [21]:

    
t = traj[0][0]



In [26]:

    
t.topology.to_openmm()









    Out[26]:





<simtk.openmm.app.topology.Topology at 0x1412853d0>



In [31]:

    
f = forcefield.createSystem(traj[0][0].topology.to_openmm())



In [ ]:



In [34]:

    
forces = f.getForces()
forces









    Out[34]:





[<simtk.openmm.openmm.HarmonicBondForce; proxy of <Swig Object of type 'OpenMM::HarmonicBondForce *' at 0x14121ca20> >,
 <simtk.openmm.openmm.HarmonicAngleForce; proxy of <Swig Object of type 'OpenMM::HarmonicAngleForce *' at 0x14121c510> >,
 <simtk.openmm.openmm.PeriodicTorsionForce; proxy of <Swig Object of type 'OpenMM::PeriodicTorsionForce *' at 0x14121c780> >,
 <simtk.openmm.openmm.NonbondedForce; proxy of <Swig Object of type 'OpenMM::NonbondedForce *' at 0x14121c4b0> >,
 <simtk.openmm.openmm.CMMotionRemover; proxy of <Swig Object of type 'OpenMM::CMMotionRemover *' at 0x14121c600> >]



In [37]:

    
f0 = forces[0]
f0









    Out[37]:





<simtk.openmm.openmm.HarmonicBondForce; proxy of <Swig Object of type 'OpenMM::HarmonicBondForce *' at 0x14121ca20> >



In [ ]: