Multiple state TIS on Alanine Dipeptide
This examples uses a 6-state model of alanine dipeptide and shows how to use multiple state TIS. This is based on previous work by Du and Bolhuis.
This uses the state definition from [1]. 6 states named A,B,C,D,E and F
[1] W.-N. Du, K. A. Marino, and P. G. Bolhuis, “Multiple state transition interface sampling of alanine dipeptide in explicit solvent,” J. Chem. Phys., vol. 135, no. 14, p. 145102, 2011.
Import and general setup¶
First we tell the ipython notebook to put figures inside the notebook
%matplotlib inline
%config InlineBackend.figure_format = 'svg'
And import lots of modules
from __future__ import print_function
# standard packages
import numpy as np
import mdtraj as md
import pandas as pd
import math
import random
import time
# helpers for phi-psi plotting
import alatools as ala
import matplotlib.pyplot as plt
# OpenMM
from simtk.openmm import app
import simtk.unit as unit
# OpenMMTools
import openmmtools as omt
# OpenPathSampling
import openpathsampling as paths
from openpathsampling.tools import refresh_output
# Visualization of PathTrees
import openpathsampling.visualize as ops_vis
#from openpathsampling.visualize import PathTree
from IPython.display import SVG
# the openpathsampling OpenMM engine
import openpathsampling.engines.openmm as eng
Set simulation options and create a simulator object¶
Create an AlanineOpenMMSimulator for demonstration purposes
We will need a openmm.System
object and an openmm.Integrator
object.
To learn more about OpenMM, read the OpenMM documentation. The code we use here is based on output from the convenient web-based OpenMM builder.
We first create a snapshot using an existing PDB file. The contained topology is necessary for the OpenMM system object to compute all forces.
pdb_file = "../resources/AD_initial_frame.pdb"
template = eng.snapshot_from_pdb(pdb_file)
1. the force field¶
using AMBER96 as in the original paper with Tip3P water.
forcefield = app.ForceField('amber96.xml', 'tip3p.xml')
2. the system object¶
and we use the template to create the necessary OpenMM Topology
object. Note that an openmm topology is (a little bit) different from an mdtraj.Topology
so we need to convert.
pdb = app.PDBFile(pdb_file)
system = forcefield.createSystem(
pdb.topology,
nonbondedMethod=app.PME,
nonbondedCutoff=1.0*unit.nanometers,
constraints=app.HBonds,
rigidWater=True,
ewaldErrorTolerance=0.0005
)
3. the integrator¶
integrator_low = omt.integrators.VVVRIntegrator(
temperature=300 * unit.kelvin, # temperature
timestep=2.0 * unit.femtoseconds # integration step size
)
integrator_low.setConstraintTolerance(0.00001)
4. the platform¶
we let OpenMM decide to use the fastest platform available. You can ask an OpenMM engine what platform it is using with the .platform
attribute. You can actually override this behaviour when you initialize the engine and pass either a string or a platform object.
print(eng.Engine.available_platforms())
5. OpenMM properties¶
There are lots of options. For speed we pick mixed precision on GPUs. Most GPU programming libraries have difficulty with double precision. For our purposes this should be fine and it is faster than double precision. These are the options passed
platform = 'CPU'
if platform == 'OpenCL':
openmm_properties = {
'OpenCLPrecision': 'mixed',
'OpenCLDeviceIndex': "2" # pick the correct index if you have several instances
}
elif platform == 'CUDA':
openmm_properties = {'CudaPrecision': 'mixed'}
elif platform == 'CPU':
openmm_properties = {}
else:
openmm_properties = {}
6. OPS options¶
An engine in OpenPathSampling also has general options that are the same for any engine, like the number of steps until a frame is stored or the maximal number of steps until an engine will automatically stop.
engine_low_options = {
'n_frames_max': 5000,
'n_steps_per_frame': 10
}
Finally build the main engine from the parts
engine_low = eng.Engine(
template.topology,
system,
integrator_low,
openmm_properties=openmm_properties,
options=engine_low_options
)
engine_low.name = 'default'
For the exploration of state space and getting an initial trajectory we also want to simulate at a very high temperature. Here 750K. This is totally unphysical, but all we want is to generate conformations that are not completely wrong in the sense that they could exist at low temperatures.
Set a high temperature simulation for exploration.
For OpenMM all we need to change is the replace the integrator with a new temperature and half the stepsize. We will later pick only every second frame
integrator_high = omt.integrators.VVVRIntegrator(
temperature=1000 * unit.kelvin, # temperature
timestep=2.0 * unit.femtoseconds # integration step size
)
integrator_high.setConstraintTolerance(0.00001)
engine_high_options = {
'n_frames_max': 5000,
'n_steps_per_frame': 10 # twice as many steps with half stepsize
}
An clone the engine using a new integrator
engine_high = engine_low.from_new_options(
integrator=integrator_high,
options=engine_high_options)
engine_high.name = 'high'
For now we use the default engine at
engine_high.initialize(platform)
print('High-Engine uses')
print('platform `%s`' % engine_high.platform)
print('temperature `%6.2f K' % (
float(engine_high.integrator.getGlobalVariableByName('kT')) / 0.0083144621))
Create the storage¶
We open a new empty storage
storage_file = 'ala_mstis_bootstrap.nc'
storage = paths.Storage(storage_file, 'w', template=template)
And store both engines in it. Since these are named you can load these later using storage.engines[name]
.
storage.save(engine_low);
storage.save(engine_high);
And store a template for convenience.
storage.tag['template'] = template
State Definitions¶
We define states A-F. The specifications are taken from the paper [^1].
The list of state names
states = ['A', 'B', 'C', 'D', 'E', 'F']
If you just want to play around you might want to just simulate the first 4 states, which is much faster. You can do so by uncommenting the line below
states = ['A', 'B', 'C', 'D']
Define the centers.
state_centers = {
'A' : [-150, 150],
'B' : [-70, 135],
'C' : [-150, -65],
'D' : [-70, -50],
'E' : [50, -100],
'F' : [40, 65]
}
And the radii of the interfaces
interface_levels = {
'A' : [20, 45, 65, 80],
'B' : [20, 45, 65, 75],
'C' : [20, 45, 60],
'D' : [20, 45, 60],
'E' : [20, 45, 65, 80],
'F' : [20, 45, 65, 80],
}
And also save these information for later convenience
storage.tag['states'] = states
storage.tag['state_centers'] = state_centers
storage.tag['interface_levels'] = interface_levels
Order Parameters¶
this generates an order parameter (callable) object named psi (so if we call psi(trajectory)
we get a list of the values of psi for each frame in the trajectory). This particular order parameter uses mdtraj's compute_dihedrals function, with the atoms in psi_atoms.
The .with_diskcache
will tell the CVs to also create some space to save the values of this CV in the same storage where the CV gets saved. Usually you want to save the direct CVs with diskcache while CVs that build upon these CVs should be fine. The reason is this:
Loading large snapshots from disk is expensive, while compute something from snapshots in memory is fast. So for analysis re-computing values is fine as long as we do not need to reload the snapshots into memory. In our case this means that the direct CVs (CVs that require accessing snapshots properties) are phi
and psi
since everything else is based upon this, while the distance to the state center opA
to opF
is an indirect CV and should not be stored.
psi_atoms = [6,8,14,16]
psi = paths.MDTrajFunctionCV(
name="psi",
f=md.compute_dihedrals,
topology=template.topology,
indices=[psi_atoms]
).with_diskcache()
phi_atoms = [4,6,8,14]
phi = paths.MDTrajFunctionCV(
name="phi",
f=md.compute_dihedrals,
topology=template.topology,
indices=[phi_atoms]
).with_diskcache()
storage.save([psi, phi]);
Define a function that defines a distance in periodic $\phi,\psi$-space.
def circle_degree(snapshot, center, phi, psi):
import numpy
p = numpy.array([phi(snapshot), psi(snapshot)]) / numpy.pi * 180.0
delta = numpy.abs(center - p)
delta = numpy.where(delta > 180.0, delta - 360.0, delta)
return numpy.hypot(delta[0], delta[1])
Create CVs for all states by using the circle_degree
function with different centers
cv_state = dict()
for state in state_centers:
op = paths.FunctionCV(
name = 'op' + state,
f=circle_degree,
center=np.array(state_centers[state]),
phi=phi,
psi=psi
)
cv_state[state] = op
Volumes¶
Volume define regions in state space using given CVs. We define here the regions around the state centers. Their boundaries correspond to the interfaces as used for TIS. Crossing an interface thus corresponds to leaving an interface volume into the next larger one.
interface_sets = {}
for state, levels in interface_levels.items():
print(levels)
interface_sets[state] = \
paths.VolumeInterfaceSet(cv_state[state], 0.0, levels)
Create Volume
objects for all states. In our setup we chose the lowest interface volume for the state definition .
vol_state = {}
for state, levels in interface_levels.items():
# vol_state[state] = interface_sets[state][0]
vol_state[state] = paths.VolumeInterfaceSet(cv_state[state], 0.0, 10)[0]
vol_state[state].name = state
Visualize in Phi/Psi space¶
We use a littl helper function specific for phi/psi plots to illustrate what is happening. This code will not affect the generation of data but will help in visualizing the system
# reload(ala)
plot = ala.TwoCVSpherePlot(
cvs=(phi, psi),
states=[vol_state[vol] for vol in states],
state_centers=[state_centers[vol] for vol in states],
interface_levels=[interface_levels[vol] for vol in states]
)
plot.new()
# plot the phi/psi plot and show all states and interfaces
plot.main()
# add the initial template
plot.add_snapshot(template, label='initial')
Set up the MSTIS network¶
We want to simulate using multistate transition interface sampling and to alleviate
the effort of setting up all states, interfaces, appropriate pathensembles and
a scheme on how to generate samples we use some helper construct like in the other
examples. In our case the MSTISNetwork
comes first and contains the definitions of
all states, their interfaces as well as associated CVs.
To make use of the optional Multi-State Outer Interface MSOuterInterface
we need to define it. It will allow reversing samples that connect two states and overcomes the issue that in all present moves the initial state needs to remain the same.
ms_outers = paths.MSOuterTISInterface.from_lambdas(
{ifaces: max(ifaces.lambdas) + 4
for state, ifaces in interface_sets.items()}
)
mstis = paths.MSTISNetwork([
(vol_state[state], interface_sets[state]) # core, interface set
for state in states],
ms_outers=ms_outers
)
And save the network
storage.save(paths.Trajectory([template]))
storage.tag['network'] = mstis
Set up the MoveScheme
¶
Next, we need to specify the actual way in which we want to simulate the Network. This could be done by
constructing a suitable PathMover yourself or employing the power of a MoveScheme that will generate all nevessary movers for you. The MoveScheme also effectively knows which initial samples in which ensembles you need. Beside the ones indirectly defined by the network the movescheme might not need some of these because certain moves will not be used. It could also happen that a movescheme require an extra
ensemble to work for special moves that are not directly known from the network definition.
In our case we will use the DefaultScheme
that uses RepEx, Reversal, Shooting and MinusMoves.
scheme = paths.DefaultScheme(mstis)
Initial trajectories¶
Next, we want to generate initial pathways for all ensembles, e.g. The TISEnsembles
and MinusInterfaceEnsembles
. There are several ways to do so. A very easy approach is to explore the state space at a high temperature and wait until all states have been visited. Splitting the long trajectory appropriately would give us at least one valid trajectory for each TIS ensemble. In a second step we can use over sub parts of the long trajectory and extend these into minus ensembles.
So let's try this: We generate an ensemble that is True
if a trajectory is at least completely outside one state A-F. It will report false once all states have been hit, which is what we want. We will ues this as a stopping condition to generate a first long trajectory at high temperature.
state_volumes = list(vol_state.values())
visit_all_states = paths.join_ensembles([paths.AllOutXEnsemble(state) for state in state_volumes])
# trajectory = engine_high.generate(template, [hit_all_states.can_append])
This will take some time, even when running at a high temperature. Since this type of ensemble is quite common, and since we'd also like get progress updates while the trajectory is running, we created a specialized ensemble that wraps around the idea sketched out above. This ensemble gives a trajectory that will have visited all states, and reports progress while it is running.
visit_all_states = paths.VisitAllStatesEnsemble(state_volumes)
trajectory = engine_high.generate(template, [visit_all_states.can_append])
Output the long trajectory and see if we hit all states
plot.new()
plot.main()
plot.add_trajectory(trajectory)