Included in this notebook:
In [1]:
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.
In [2]:
%%time
storage = paths.AnalysisStorage("ala_mstis_production.nc")
In [3]:
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))
In [4]:
%%time
mstis = storage.networks[0]
In [5]:
%%time
for cv in storage.cvs:
print(cv.name, cv._store_dict)
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' : 2, 'bin_range' : (0.0, 90) }
mstis.hist_args['pathlength'] = { 'bin_width' : 5, 'bin_range' : (0, 100) }
In [7]:
%%time
mstis.rate_matrix(storage.steps, force=True)
Out[7]:
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.
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)
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
Out[9]:
We normally look at these on a log scale:
In [10]:
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))
Out[10]:
In [11]:
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
Out[11]:
In [12]:
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
Out[12]:
In [13]:
hists_A = mstis.transitions[(stateA, stateB)].histograms
hists_B = mstis.transitions[(stateB, stateC)].histograms
hists_C = mstis.transitions[(stateC, stateB)].histograms
In [14]:
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 [15]:
# add visualization of the sum
In [16]:
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)
In [17]:
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))
In [18]:
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 [19]:
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 [20]:
scheme = storage.schemes[0]
In [21]:
scheme.move_summary(storage.steps)
In [22]:
scheme.move_summary(storage.steps, 'shooting')
In [23]:
scheme.move_summary(storage.steps, 'minus')
In [24]:
scheme.move_summary(storage.steps, 'repex')
In [25]:
scheme.move_summary(storage.steps, 'pathreversal')
In [26]:
repx_net = paths.ReplicaNetwork(scheme, storage.steps)
In [27]:
repx_net.mixing_matrix()
Out[27]:
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 [28]:
repxG = paths.ReplicaNetworkGraph(repx_net)
repxG.draw('spring')
In [29]:
import openpathsampling.visualize as vis
#reload(vis)
from IPython.display import SVG
In [30]:
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())
Out[30]:
In [32]:
decorrelated = tree.generator.decorrelated
print ("We have " + str(len(decorrelated)) + " decorrelated trajectories.")
In [ ]:
In [ ]: