Included in this notebook:
NOTE: This notebook uses our old analysis approach. See the "new analysis" appendix notebook for details on how to customize analysis.
In [1]:
from __future__ import print_function
In [2]:
# If our large test file is available, use it. Otherwise, use file generated
# from toy_mstis_2_run.ipynb. This is so the notebook can be used in testing.
import os
test_file = "../toy_mstis_1k_OPS1.nc"
filename = test_file if os.path.isfile(test_file) else "mstis.nc"
print("Using file `"+ filename + "` for analysis")
In [3]:
%matplotlib inline
import matplotlib.pyplot as plt
import openpathsampling as paths
import numpy as np
In [4]:
%%time
storage = paths.Storage(filename, mode='r')
In [5]:
# the following works with the old file we use in testing; the better way is:
# mstis = storage.networks['mstis'] # when objects are named, use the name
mstis = storage.networks.first
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.
In [6]:
mstis.hist_args['max_lambda'] = { 'bin_width' : 0.05, 'bin_range' : (0.0, 0.5) }
mstis.hist_args['pathlength'] = { 'bin_width' : 5, 'bin_range' : (0, 150) }
In [7]:
%%time
mstis.rate_matrix(storage.steps)
Out[7]:
The self-rates (the rate of returning the to initial state) are undefined, and return not-a-number.
The rate is calculated 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.
In [8]:
stateA = storage.volumes["A"]
stateB = storage.volumes["B"]
stateC = storage.volumes["C"]
In [9]:
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, '-r')
plt.plot(tcp_CA.x, tcp_CA, '-k')
plt.plot(tcp_BC.x, tcp_BC, '-b')
plt.plot(tcp_AC.x, tcp_AC, '-g') # same as tcp_AB in MSTIS
Out[9]:
We normally look at these on a log scale:
In [10]:
plt.plot(tcp_AB.x, np.log(tcp_AB), '-r')
plt.plot(tcp_CA.x, np.log(tcp_CA), '-k')
plt.plot(tcp_BC.x, np.log(tcp_BC), '-b')
plt.xlim(0.0, 1.0)
Out[10]:
Now, in case you want to know the total crossing probabability at each interface (for example, to use as a bias in an SRTIS calculation):
In [11]:
# TODO: MOVE THESE TO A METHOD INSIDE THE CODE; MAKE THEM WORK WITH NEW ANALYSIS
In [12]:
import pandas as pd
def crossing_probability_table(transition):
tcp = transition.tcp
interface_lambdas = transition.interfaces.lambdas
values = [tcp(x) for x in interface_lambdas]
return pd.Series(values, index=interface_lambdas, name=transition.name)
In [13]:
def outer_crossing_probability(transition):
tcp = transition.tcp
interface_outer_lambda = transition.interfaces.lambdas[-1]
return tcp(interface_outer_lambda)
In [14]:
crossing_probability_table(mstis.from_state[stateA])
Out[14]:
In [15]:
outer_crossing_probability(mstis.from_state[stateA])
Out[15]:
In [16]:
tcp_AB(mstis.from_state[stateA].interfaces.lambdas[-1])
Out[16]:
In [17]:
tcp_A = mstis.from_state[stateA].tcp
In [18]:
flux_dict = {
(transition.stateA, transition.interfaces[0]): transition._flux
for transition in mstis.transitions.values()
}
flux_dict
Out[18]:
In [19]:
paths.analysis.tis.flux_matrix_pd(flux_dict)
Out[19]:
In [20]:
state_names = [s.name for s in mstis.states]
outer_ctp_matrix = pd.DataFrame(columns=state_names, index=state_names)
for state_pair in mstis.transitions:
transition = mstis.transitions[state_pair]
outer_ctp_matrix.at[state_pair[0].name, state_pair[1].name] = transition.ctp[transition.ensembles[-1]]
outer_ctp_matrix
Out[20]:
In [21]:
state_pair_names = [(t[0].name, t[1].name) for t in mstis.transitions]
ctp_by_interface = pd.DataFrame(index=state_pair_names)
for state_pair in mstis.transitions:
transition = mstis.transitions[state_pair]
for ensemble_i in range(len(transition.ensembles)):
state_pair_name = (transition.stateA.name, transition.stateB.name)
ctp_by_interface.at[state_pair_name, ensemble_i] = transition.conditional_transition_probability(
storage.steps,
transition.ensembles[ensemble_i]
)
ctp_by_interface
Out[21]:
In [22]:
hists_A = mstis.transitions[(stateA, stateB)].histograms
hists_B = mstis.transitions[(stateB, stateC)].histograms
hists_C = mstis.transitions[(stateC, stateB)].histograms
In [23]:
hists = {'A': hists_A, 'B': hists_B, 'C': hists_C}
plot_style = {'A': '-r', 'B': '-b', 'C': '-k'}
In [24]:
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)
In [25]:
# add visualization of the sum
In [26]:
for hist_type in hists:
hist = hists[hist_type]
for ens in hist['max_lambda']:
reverse_cumulative = hist['max_lambda'][ens].reverse_cumulative()
plt.plot(reverse_cumulative.x, reverse_cumulative, plot_style[hist_type])
plt.xlim(0.0, 1.0)
Out[26]:
In [27]:
for hist_type in hists:
hist = hists[hist_type]
for ens in hist['max_lambda']:
reverse_cumulative = hist['max_lambda'][ens].reverse_cumulative()
plt.plot(reverse_cumulative.x, np.log(reverse_cumulative), plot_style[hist_type])
plt.xlim(0.0, 1.0)
Out[27]:
In [28]:
for hist in [hists_A, hists_B, hists_C]:
for ens in hist['pathlength']:
normalized = hist['pathlength'][ens].normalized()
plt.plot(normalized.x, normalized)
In [29]:
for ens in hists_A['pathlength']:
normalized = hists_A['pathlength'][ens].normalized()
plt.plot(normalized.x, normalized)
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.
In [ ]:
In [30]:
scheme = storage.schemes[0]
In [31]:
scheme.move_summary(storage.steps)
In [32]:
scheme.move_summary(movers='shooting')
In [33]:
scheme.move_summary(movers='minus')
In [34]:
scheme.move_summary(movers='repex')
In [35]:
scheme.move_summary(movers='pathreversal')
In [36]:
repx_net = paths.ReplicaNetwork(scheme, storage.steps)
In [37]:
repx_net.mixing_matrix()
Out[37]:
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.
In [38]:
repxG = paths.ReplicaNetworkGraph(repx_net)
repxG.draw('spring')
In [ ]:
In [39]:
import openpathsampling.visualize as vis
from IPython.display import SVG
In [40]:
tree = vis.PathTree(
storage.steps[0:500],
vis.ReplicaEvolution(replica=2, accepted=True)
)
SVG(tree.svg())
Out[40]:
In [41]:
decorrelated = tree.generator.decorrelated
print("We have " + str(len(decorrelated)) + " decorrelated trajectories.")
In [42]:
# we use the %run magic because this isn't in a package
%run ../resources/toy_plot_helpers.py
background = ToyPlot()
background.contour_range = np.arange(-1.5, 1.0, 0.1)
background.add_pes(storage.engines[0].pes)
In [43]:
xval = paths.FunctionCV("xval", lambda snap : snap.xyz[0][0])
yval = paths.FunctionCV("yval", lambda snap : snap.xyz[0][1])
visualizer = paths.StepVisualizer2D(mstis, xval, yval, [-1.0, 1.0], [-1.0, 1.0])
visualizer.background = background.plot()
In [44]:
visualizer.draw_samples(list(tree.samples))
Out[44]:
In [45]:
#! skip
# The skip directive tells our test runner not to run this cell
import time
max_step = 10
for step in storage.steps[0:max_step]:
visualizer.draw_ipynb(step)
time.sleep(0.1)
In [ ]: