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)
tps= paths.TPSNetwork.from_states_all_to_all(vol_state.values())
tps_ensemble = tps.sampling_ensembles[0]
trajectories = tps_ensemble.split(trajectory)
print(trajectories)
Find initial samples¶
It can be difficult to get initial samples for an ensemble for various reasons. And in fact for this example it might be the most tricky part in setting up the path simulation.
In theory you could just pick a frame and extend forward and backward until an ensemble reports cannot append anymore and check if that sample is in the ensemble. While this works in many cases it can be rather inefficient. For us it is better to reuse the previously computed trajectories as a basis to generate samples.
Bottomline: The next command just does what you would try to do with the long trajectory for you: Use split and extend to generate the required set of samples to start.
If you are interested in details keep reading, otherwise execute the next cell.
To do so, we use a number of functions from Ensemble that will try to generate a sample from a list of initial trajectories in different ways:
ensemble.get_sample_from_trajectories: directly search for candidate sample among the given trajectoriesensemble.split_sample_from_trajectories: split the given trajectories in valid samplesensemble.extend_sample_from_trajectories: try to extend suitable sub-trajectories into valid samples
Instead of calling these functions directly we make use of a function from MoveScheme that uses its knowledge about the necessary ensembles to run these function for all necessary ensembles to create a suitable set of initial samples.
When doing this there are several aspects that we might want to take into account
- also use reversed copies and try these
- avoid reusing samples for better decorrelation.
- try to use smaller trajectories as these are usually better equilibrated
The function scheme.initial_conditions_from_trajectories can take care of all of these. We will run it with the default setting but also use the extend-complex and extend-minimum strategies besides split. That means we try different strategies in a specific order.
split: try to split the trajectory and look for valid samplesextend-complex: try to look for suitable sub-trajectories that can be extended into valid samples. Currently this works only forMinusInterfaceEnsemblesand it will extendA-X-Atrajectories into Minus Ensembles.extend-minimum: try to look for trajectories that cross an interfaceA-Xand try to extend these into a full ensemble. This one will work forTISEnsembleandMinusInterfaceEnsemble. We will use this last method to create minus ensembles for the last state found. Since we stop when we discover it, we will not have generated a sample of typeA-X-A. This might also happen if we only visit a state once.
Note, that the creation of sample might fail: Extending might just not reach the desired target state. By default, it will make 3 attempts per extendable subtrajectory. The whole process may take some time. Basically because we need to split a very long trajectory and to use the shortest possible ones we need to really search for all possible sub-trajectories.
total_sample_set = scheme.initial_conditions_from_trajectories(
trajectories)
Let's see, if we were successful. If so, the next function should not show any missing or extra ensembles.
print(scheme.initial_conditions_report(total_sample_set))
plt.figure(21, (12,5))
plt.subplot(121)
plt.title('High temperature')
plot.main()
for s in total_sample_set:
plot.add_trajectory(s.trajectory, line=True)
If it does, this means that our attempt to extend was probably not successful. So we can just try again until it does. To extend an existing sample_set you need to pass it to the function.
# loop until we are done
while scheme.check_initial_conditions(total_sample_set)[0]:
total_sample_set = scheme.initial_conditions_from_trajectories(
trajectory,
sample_set=total_sample_set, # this is new
strategies=[ # we omit the `split` strategy
'extend-complex',
'extend-minimal'],
engine=engine_high)
plt.figure(21, (12,5))
plt.subplot(121)
plt.title('High temperature')
plot.main()
for s in total_sample_set:
plot.add_trajectory(s.trajectory, line=True)
Done.
Equilibration¶
For the next steps we need the engine running at the desired (room) temperature. So we create a new engine. Since we might run on a local machine and do not have several instances of GPU available we will unload the existing engine and create a new one.
engine_high.unload_context() # delete the current openmm.Context object
engine_low.initialize(platform)
refresh_output('Low-Engine uses', refresh=False)
print('platform `%s`' % engine_low.platform)
print('temperature `%6.2f K' % (
float(engine_low.integrator.getGlobalVariableByName('kT')) / 0.0083144621))
In molecular dynamics, you need to equilibrate if you don't start with an equilibrium frame (e.g., if you start with solvent molecules on a grid, your system should equilibrate before you start taking statistics). Similarly, if you start with a set of paths which are far from the path ensemble equilibrium, you need to equilibrate. This could either be because your trajectories are not from the real dynamics (generated with metadynamics, high temperature, etc.) or because your trajectories are not likely representatives of the path ensemble (e.g., if you put transition trajectories into all interfaces).
As with MD, running equilibration can be the same process as running the total simulation. However, in path sampling, it doesn't have to be: we can equilibrate without replica exchange moves or path reversal moves, for example. In the example below, we create a MoveScheme that only includes shooting movers.
equil_scheme = paths.OneWayShootingMoveScheme(mstis, engine=engine_low)
Note, that some Schemes may change the actual required set of initial samples. For
the equilibration we apply shooting moves only to TIS ensembles and not to the
MinusInterfaceEnsembles.
Create the PathSampler object that can run the equilibration.
equilibration = paths.PathSampling(
storage=storage,
sample_set=total_sample_set,
move_scheme=equil_scheme,
)
And run lots of equilibration steps. This is not really necessary, but we generated from the wrong path distribution and this should give more likeli paths for the right temperature.
equilibration.run(500)
equilibrated_sset = equilibration.sample_set
Let's have a quick look how many frames in TIS ensembles are already generated by the default engine at 300K and how many were done using the high engine at 1000K.
engine_list = {}
unique_snapshots = set(sum(
[
samp.trajectory for samp in equilibrated_sset
if isinstance(samp.ensemble, paths.TISEnsemble)
],
[]))
for samp in equilibrated_sset:
if isinstance(samp.ensemble, paths.TISEnsemble):
for snap in samp.trajectory:
eng = snap.engine
engine_list[eng] = engine_list.get(eng, 0) + 1
for eng, counts in engine_list.items():
print('%20s : %6d' % (eng.name, counts))
Finally, save the list of samples for later without any reference to previous samples. This is necessary since usually a sample knows its origin sample and saving a sample as is, would trigger saving the whole history of all intermediate samples used. Since we are not interested in the history of equilibration, this saves a lot of storage.
storage.tag['sampleset'] = equilibrated_sset.copy_without_parents()
Let's have a final look at what we cover with all our initial samples.
for s in total_sample_set:
print(s)
for s in equilibrated_sset:
print(s)
plt.figure(21, (12,5))
plt.subplot(121)
plt.title('High temperature')
plot.main()
for s in total_sample_set:
plot.add_trajectory(s.trajectory, line=True)
plt.subplot(122)
plt.title('Equilibrated')
plot.main()
plot.zoom = 180/3.1415926
for s in equilibrated_sset:
plot.add_trajectory(s.trajectory, line=True)
A final check that we still have a valid sampleset
equilibrated_sset.sanity_check()
print('snapshots:', len(storage.snapshots))
print('trajectories:', len(storage.trajectories))
print('samples:', len(storage.samples))
print('filesize:', storage.file_size_str)
storage.close()
storage =paths.AnalysisStorage("ala_mstis_bootstrap.nc")
print("PathMovers:", len(storage.pathmovers))
print("Engines:", len(storage.engines))
print("Samples:", len(storage.samples))
print("Ensembles:", len(storage.ensembles))
print("SampleSets:", len(storage.samplesets))
print("Snapshots:", len(storage.snapshots))
print("Trajectories:", len(storage.trajectories))
print("Networks:", len(storage.networks))
scheme = storage.schemes[0]
scheme.move_summary(storage.steps)
scheme.move_summary(storage.steps, 'shooting')
tree = ops_vis.PathTree(
storage.steps[0:500],
ops_vis.ReplicaEvolution(
replica=0
)
)
f = open("myfile.svg","w")
f.write(tree.svg())
f.close()
SVG(tree.svg())
path_lengths = [len(step.active[7].trajectory) for step in storage.steps]
plt.hist(path_lengths, bins=40, alpha=0.5);
print("Maximum:", max(path_lengths), "("+str(max(path_lengths)*engine_low.snapshot_timestep)+")")
print("Average:", "{0:.2f}".format(np.mean(path_lengths)), "("+(np.mean(path_lengths)*engine_low.snapshot_timestep).format("%.3f")+")")
traj = storage.steps[-1].active[0].trajectory
len(traj)
print((phi(traj))*(180/np.pi))
print((psi(traj))*(180/np.pi))
(AD_mstis_1_setup.ipynb; AD_mstis_1_setup.py)
Run from bootstrap paths¶
Now we will use the initial trajectories we obtained from bootstrapping to run an MSTIS simulation. This will show both how objects can be regenerated from storage and how regenerated equivalent objects can be used in place of objects that weren't stored.
Tasks covered in this notebook:
- Loading OPS objects from storage
- Ways of assigning initial trajectories to initial samples
- Setting up a path sampling simulation with various move schemes
- Visualizing trajectories while the path sampling is running
%matplotlib inline
import openpathsampling as paths
import numpy as np
import math
# the openpathsampling OpenMM engine
import openpathsampling.engines.openmm as eng
Loading things from storage¶
First we'll reload some of the stuff we stored before. Of course, this starts with opening the file.
old_store = paths.AnalysisStorage("ala_mstis_bootstrap.nc")
A lot of information can be recovered from the old storage, and so we don't have the recreate it. However, we did not save our network, so we'll have to create a new one. Since the network creates the ensembles, that means we will have to translate the trajectories from the old ensembles to new ensembles.
print("PathMovers:", len(old_store.pathmovers))
print("Engines:", len(old_store.engines))
print("Samples:", len(old_store.samples))
print("Trajectories:", len(old_store.trajectories))
print("Ensembles:", len(old_store.ensembles))
print("SampleSets:", len(old_store.samplesets))
print("Snapshots:", len(old_store.snapshots))
print("Networks:", len(old_store.networks))
Loading from storage is very easy. Each store is a list. We take the 0th snapshot as a template (it doesn't actually matter which one) for the next storage we'll create. There's only one engine stored, so we take the only one.
# template = old_store.snapshots[0]
engine = old_store.engines['default']
mstis = old_store.networks[0]
sset = old_store.tag['sampleset']
initialize engine¶
if we do not select a platform the fastest possible will be chosen but we explicitly request to use the one in the config file
#platform = 'CUDA'
#engine.initialize(platform)
print('Engine uses platform `%s`' % engine.platform)
sset.sanity_check()
Running RETIS¶
Now we run the full calculation. Up to here, we haven't been storing any of our results. This time, we'll start a storage object, and we'll save the network we've created. Then we'll run a new PathSampling calculation object.
# logging creates ops_output.log file with details of what the calculation is doing
#import logging.config
#logging.config.fileConfig("logging.conf", disable_existing_loggers=False)
storage = paths.storage.Storage("ala_mstis_production.nc", "w")
storage.snapshots.save(old_store.snapshots[0]);
Before we can sample we still need to set the actual MoveScheme which determines the
set of moves to apply to our set of samples and effectively doing the steps in
replica (sampleset) space. We pick the default scheme for mstis and feed it with
the engine to be used.
scheme = paths.DefaultScheme(mstis, engine)
and finally generate the PathSampler object to conduct the simulation.
mstis_calc = paths.PathSampling(
storage=storage,
sample_set=sset,
move_scheme=scheme
)
mstis_calc.save_frequency = 10
#mstis_calc.save_frequency = 1
Now everything is ready: let's run the simulation! The first step takes a little since all necessary information, i.e. the engines, topologies, initial snapshots, ..., need to be stored. Then the monte carlo steps will be performed.
mstis_calc.run(10000)
print(len(storage.steps))
# commented out during development, so we can "run all" and then do more
storage.close()
(AD_mstis_2_run.ipynb; AD_mstis_2_run.py)
Analyzing the MSTIS simulation¶
Included in this notebook:
- Opening files for analysis
- Rates, fluxes, total crossing probabilities, and condition transition probabilities
- Per-ensemble properties such as path length distributions and interface crossing probabilities
- Move scheme analysis
- Replica exchange analysis
- Replica move history tree visualization
- Replaying the simulation
- MORE TO COME! Like free energy projections, path density plots, and more
from __future__ import print_function
%matplotlib inline
import matplotlib.pyplot as plt
import openpathsampling as paths
import numpy as np
The optimum way to use storage depends on whether you're doing production or analysis. For analysis, you should open the file as an AnalysisStorage object. This makes the analysis much faster.
%%time
storage = paths.AnalysisStorage("ala_mstis_production.nc")
print("PathMovers:", len(storage.pathmovers))
print("Engines:", len(storage.engines))
print("Samples:", len(storage.samples))
print("Ensembles:", len(storage.ensembles))
print("SampleSets:", len(storage.samplesets))
print("Snapshots:", len(storage.snapshots))
print("Trajectories:", len(storage.trajectories))
print("Networks:", len(storage.networks))
%%time
mstis = storage.networks[0]
%%time
for cv in storage.cvs:
print(cv.name, cv._store_dict)
Reaction rates¶
TIS methods are especially good at determining reaction rates, and OPS makes it extremely easy to obtain the rate from a TIS network.
Note that, although you can get the rate directly, it is very important to look at other results of the sampling (illustrated in this notebook and in notebooks referred to herein) in order to check the validity of the rates you obtain.
By default, the built-in analysis calculates histograms the maximum value of some order parameter and the pathlength of every sampled ensemble. You can add other things to this list as well, but you must always specify histogram parameters for these two. The pathlength is in units of frames.
mstis.hist_args['max_lambda'] = { 'bin_width' : 2, 'bin_range' : (0.0, 90) }
mstis.hist_args['pathlength'] = { 'bin_width' : 5, 'bin_range' : (0, 100) }
%%time
mstis.rate_matrix(storage.steps, force=True)
The self-rates (the rate of returning the to initial state) are undefined, and return not-a-number.
The rate is calcuated according to the formula:
$$k_{AB} = \phi_{A,0} P(B|\lambda_m) \prod_{i=0}^{m-1} P(\lambda_{i+1} | \lambda_i)$$
where $\phi_{A,0}$ is the flux from state A through its innermost interface, $P(B|\lambda_m)$ is the conditional transition probability (the probability that a path which crosses the interface at $\lambda_m$ ends in state B), and $\prod_{i=0}^{m-1} P(\lambda_{i+1} | \lambda_i)$ is the total crossing probability. We can look at each of these terms individually.
Total crossing probability¶
stateA = storage.volumes["A"]
stateB = storage.volumes["B"]
stateC = storage.volumes["C"]
tcp_AB = mstis.transitions[(stateA, stateB)].tcp
tcp_AC = mstis.transitions[(stateA, stateC)].tcp
tcp_BC = mstis.transitions[(stateB, stateC)].tcp
tcp_BA = mstis.transitions[(stateB, stateA)].tcp
tcp_CA = mstis.transitions[(stateC, stateA)].tcp
tcp_CB = mstis.transitions[(stateC, stateB)].tcp
plt.plot(tcp_AB.x, tcp_AB)
plt.plot(tcp_CA.x, tcp_CA)
plt.plot(tcp_BC.x, tcp_BC)
plt.plot(tcp_AC.x, tcp_AC) # same as tcp_AB in MSTIS
We normally look at these on a log scale:
plt.plot(tcp_AB.x, np.log(tcp_AB))
plt.plot(tcp_CA.x, np.log(tcp_CA))
plt.plot(tcp_BC.x, np.log(tcp_BC))
Flux¶
Here we also calculate the flux contribution to each transition. The flux is calculated based on
import pandas as pd
flux_matrix = pd.DataFrame(columns=mstis.states, index=mstis.states)
for state_pair in mstis.transitions:
transition = mstis.transitions[state_pair]
flux_matrix.set_value(state_pair[0], state_pair[1], transition._flux)
flux_matrix
Conditional transition probability¶
outer_ctp_matrix = pd.DataFrame(columns=mstis.states, index=mstis.states)
for state_pair in mstis.transitions:
transition = mstis.transitions[state_pair]
outer_ctp_matrix.set_value(state_pair[0], state_pair[1], transition.ctp[transition.ensembles[-1]])
outer_ctp_matrix
Path ensemble properties¶
hists_A = mstis.transitions[(stateA, stateB)].histograms
hists_B = mstis.transitions[(stateB, stateC)].histograms
hists_C = mstis.transitions[(stateC, stateB)].histograms
Interface crossing probabilities¶
We obtain the total crossing probability, shown above, by combining the individual crossing probabilities of
for hist in [hists_A, hists_B, hists_C]:
for ens in hist['max_lambda']:
normalized = hist['max_lambda'][ens].normalized()
plt.plot(normalized.x, normalized)
# add visualization of the sum
for hist in [hists_A, hists_B, hists_C]:
for ens in hist['max_lambda']:
reverse_cumulative = hist['max_lambda'][ens].reverse_cumulative()
plt.plot(reverse_cumulative.x, reverse_cumulative)
for hist in [hists_A, hists_B, hists_C]:
for ens in hist['max_lambda']:
reverse_cumulative = hist['max_lambda'][ens].reverse_cumulative()
plt.plot(reverse_cumulative.x, np.log(reverse_cumulative))
Path length histograms¶
for hist in [hists_A, hists_B, hists_C]:
for ens in hist['pathlength']:
normalized = hist['pathlength'][ens].normalized()
plt.plot(normalized.x, normalized)
for ens in hists_A['pathlength']:
normalized = hists_A['pathlength'][ens].normalized()
plt.plot(normalized.x, normalized)
Sampling properties¶
The properties we illustrated above were properties of the path ensembles. If your path ensembles are sufficiently well-sampled, these will never depend on how you sample them.
But to figure out whether you've done a good job of sampling, you often want to look at properties related to the sampling process. OPS also makes these very easy.
Move scheme analysis¶
scheme = storage.schemes[0]
scheme.move_summary(storage.steps)
scheme.move_summary(storage.steps, 'shooting')
scheme.move_summary(storage.steps, 'minus')
scheme.move_summary(storage.steps, 'repex')
scheme.move_summary(storage.steps, 'pathreversal')
Replica exchange sampling¶
See the notebook repex_networks.ipynb for more details on tools to study the convergence of replica exchange. However, a few simple examples are shown here. All of these are analyzed with a separate object, ReplicaNetwork.
repx_net = paths.ReplicaNetwork(scheme, storage.steps)
Replica exchange mixing matrix¶
repx_net.mixing_matrix()
Replica exchange graph¶
The mixing matrix tells a story of how well various interfaces are connected to other interfaces. The replica exchange graph is essentially a visualization of the mixing matrix (actually, of the transition matrix -- the mixing matrix is a symmetrized version of the transition matrix).
Note: We're still developing better layout tools to visualize these.
repxG = paths.ReplicaNetworkGraph(repx_net)
repxG.draw('spring')
Replica exchange flow¶
Replica flow is defined as TODO
Flow is designed for calculations where the replica exchange graph is linear, which ours clearly is not. However, we can define the flow over a subset of the interfaces.
Replica move history tree¶
import openpathsampling.visualize as vis
#reload(vis)
from IPython.display import SVG
tree = vis.PathTree(
[step for step in storage.steps if not isinstance(step.change, paths.EmptyMoveChange)],
vis.ReplicaEvolution(replica=3, accepted=False)
)
tree.options.css['width'] = 'inherit'
SVG(tree.svg())
decorrelated = tree.generator.decorrelated
print ("We have " + str(len(decorrelated)) + " decorrelated trajectories.")
Visualizing trajectories¶
Histogramming data (TODO)¶