In [1]:
%matplotlib inline
import matplotlib.pyplot as plt
import openpathsampling as paths
import numpy as np
import pandas as pd
pd.options.display.max_rows = 10
In [2]:
storage = paths.Storage("committor_results.nc", "r")
In [3]:
phi = storage.cvs['phi']
psi = storage.cvs['psi']
In [4]:
%%time
C_7eq = storage.volumes['C_7eq']
alpha_R = storage.volumes['alpha_R']
experiments = storage.tag['experiments']
The experiments object is a list of tuples (snapshot, final_state). Each snapshot is an OPS snapshot object (a point in phase space), and the final_state is either the C_7eq object or the alpha_R object.
In [5]:
%%time
committor_analyzer = paths.ShootingPointAnalysis.from_individual_runs(experiments)
Before going further, let's talk a little bit about the implementation of the ShootingPointAnalysis object. The main thing to understand is that the purpose of that object is to histogram according to configuration. The first snapshot encountered is kept as a representative of that configuration.
So whereas there are 10000 snapshots in experiments (containing the full data, including velocities), there are only 1000 entries in the committor_analyzer (because, in this data set, I ran 1000 snapshots with 10 shots each.)
In [6]:
committor_analyzer.to_pandas()
Out[6]:
You can also pass it a function that takes a snapshot and returns a (hashable) value. That value will be used for the index. These collective variables return numpy arrays, so we need to cast the 1D array to a float.
In [7]:
psi_hash = lambda x : float(psi(x))
committor_analyzer.to_pandas(label_function=psi_hash)
Out[7]:
You can also directly obtain the committor as a dictionary of (representative) snapshot to committor value. The committor here is defines as the probability of ending in a given state, so you must give the state.
In [8]:
committor = committor_analyzer.committor(alpha_R)
In [9]:
# show the first 10 values
{k: committor[k] for k in committor.keys()[:10]}
Out[9]:
In [10]:
hist1D, bins = committor_analyzer.committor_histogram(psi_hash, alpha_R, bins=20)
In [11]:
bin_widths = [bins[i+1]-bins[i] for i in range(len(bins)-1)]
plt.bar(left=bins[:-1], height=hist1D, width=bin_widths, log=True);
In [12]:
ramachandran_hash = lambda x : (float(phi(x)), float(psi(x)))
hist2D, bins_phi, bins_psi = committor_analyzer.committor_histogram(ramachandran_hash, alpha_R, bins=20)
In [13]:
# not the best, since it doesn't distinguish NaNs, but that's just a matter of plotting
plt.pcolor(bins_phi, bins_psi, hist2D.T, cmap="winter")
plt.clim(0.0, 1.0)
plt.colorbar();
The information committor_results.nc should be everything you could want, including initial velocities for every system. In principle, you'll mainly access that information using collective variables (see documentation on using MDTraj to create OPS collective variables). However, you may decide to access that information directly, so here's how you do that.
In [14]:
# let's take the first shooting point snapshot
# experiments[N][0] gives shooting snapshot for experiment N
snapshot = experiments[0][0]
OpenMM-based objects come with units. So snapshot.coordinates is a unitted value. This can be annoying in analysis, so we have a convenience snapshot.xyz to get the version without units.
In [15]:
snapshot.coordinates
Out[15]:
In [16]:
snapshot.xyz
Out[16]:
For velocities, we don't have the convenience function, but if you want to remove units from velocities you can do so with velocity / velocity.unit.
In [17]:
snapshot.velocities
Out[17]:
In [18]:
snapshot.velocities / snapshot.velocities.unit
Out[18]:
Note that snapshots include coordinates and velocities. We have several sets of initial velocities for each initial snapshot. Taking the second shooting snapshot and comparing coordinates and velocities:
In [19]:
snapshot2 = experiments[1][0]
In [20]:
np.all(snapshot.coordinates == snapshot2.coordinates)
Out[20]:
In [21]:
np.any(snapshot.velocities == snapshot2.velocities)
Out[21]: