In [4]:
from os.path import exists
import molpx
from matplotlib import pylab as plt
%matplotlib ipympl
import pyemma
import numpy as np
In [5]:
top = molpx._molpxdir(join='notebooks/data/ala2.pdb')
# What data do we have?
if exists('/group/ag_cmb/scratch/gph82/Di-Ala-nbdata/ala2.dcd'):
MD_trajfiles = ['/group/ag_cmb/scratch/gph82/Di-Ala-nbdata/ala2.dcd'] #long trajectory
elif exists('/home/guille/ala2.dcd'):
MD_trajfiles = ['/home/guille/ala2.dcd'] # extra for Stralsund
else:
MD_trajfiles = [molpx._molpxdir(join='notebooks/data/ala2.mini.xtc')] #short trajectory
In [6]:
feat = pyemma.coordinates.featurizer(top)
feat.add_backbone_torsions()
src = pyemma.coordinates.source(MD_trajfiles, features=feat)
Y = src.get_output()
In [7]:
mpx_widget_box = molpx.visualize.FES(MD_trajfiles,
top,
Y,
#proj_idxs=[1],
nbins=50,
proj_labels=['$\phi$',
'$\psi$'],
atom_selection="symbol != H",
#n_overlays=5,
#sticky=True,
#color_list='random'
)
mpx_widget_box
In [8]:
from molpx import visualize, _linkutils
from imp import reload
reload(visualize)
reload(_linkutils)
mpl_wdg_box = molpx.visualize.traj(MD_trajfiles,
top,
Y,
plot_FES = True,
#dt = dt*1e-6, tunits='ms',
max_frames=10000,
proj_idxs=[0, 1],
panel_height=2,
proj_labels=['$\phi$', '$\psi$']
)
mpl_wdg_box
In [10]:
paths_dict, idata = molpx.generate.projection_paths(MD_trajfiles,
top,
Y,
n_points=50,
proj_idxs=[0,1],
n_projs=3,
proj_dim = 3,
verbose=False,
)
In [11]:
# Choose the coordinate and the type of path
coord = 1
path_type = 'min_rmsd'
#path_type = 'min_disp'
igeom = paths_dict[coord][path_type]["geom"]
ipath = paths_dict[coord][path_type]["proj"]
# Choose the proj_idxs for the path and the FES
# to be shown
proj_idxs = [0,1]
In [12]:
plt.ioff() # Turn of interactive plotting
plt.figure(figsize=(4,4))
h, (x,y) = np.histogramdd(np.vstack(Y)[:,proj_idxs], bins=50)
plt.contourf(x[:-1], y[:-1], -np.log(h.T), alpha=.50)
plt.ion()
linked_NGL_wdg, linked_ax_widget = molpx.visualize.sample(ipath[:,proj_idxs],
igeom,
plt.gca(),
clear_lines=True,
n_smooth = 2,
plot_path=True,
#radius=True,
)
linked_NGL_wdg._set_size('4in', '4in')
from ipywidgets import HBox
HBox([linked_NGL_wdg, linked_ax_widget.canvas])
In [13]:
feat = pyemma.coordinates.featurizer(top)
#feat.add_backbone_torsions(cossin=True)
feat.add_distances(feat.topology.select('symbol != H'))
src = pyemma.coordinates.source(MD_trajfiles, features=feat)
tica = pyemma.coordinates.tica(src, lag=np.int(src.trajectory_lengths()/3000))
Y_tica = tica.get_output()
In [15]:
mpx_wdg_box = molpx.visualize.FES(MD_trajfiles,
top,
Y_tica,
n_overlays=5,
atom_selection='backbone',
#sticky=True,
#color_list='rand'
)
mpx_wdg_box
In [17]:
mpx_wdg_box = molpx.visualize.traj(MD_trajfiles,
top,
Y_tica,
plot_FES = True,
#dt = dt*1e-6, tunits='ms',
max_frames=10000,
#proj_idxs=[0,1],
panel_height=2,
projection=tica
)
mpx_wdg_box
In [20]:
paths_dict, idata = molpx.generate.projection_paths(MD_trajfiles,
top,
Y_tica,
n_points=50,
proj_idxs=[0,1],
n_projs=2,
proj_dim = 2,
verbose=False,
)
In [21]:
# Choose the coordinate and the type of path
coord = 0
path_type = 'min_rmsd'
#path_type = 'min_disp'
igeom = paths_dict[coord][path_type]["geom"]
ipath = paths_dict[coord][path_type]["proj"]
# Choose the proj_idxs for the path and the FES
# to be shown
proj_idxs = [0,1]
In [22]:
plt.figure(figsize=(4,4))
h, (x,y) = np.histogramdd(np.vstack(Y_tica)[:,proj_idxs], bins=50)
plt.contourf(x[:-1], y[:-1], -np.log(h.T), alpha=.50)
linked_wdg, axes_widget = molpx.visualize.sample(ipath[:,proj_idxs],
igeom,
plt.gca(),
clear_lines=True,
n_smooth = 1,
plot_path=True,
)
# You can even choose to add the correlations a posteriori
molpx.visualize.correlations(tica, widget=linked_wdg, proj_idxs=0)
linked_wdg.center_view()
linked_wdg
In [ ]: