molPX Di-Ala example¶
Guillermo Pérez-Hernández guille.perez@fu-berlin.de
In this notebook we will be using a trajectory of Di-Ala-peptide to easily identify conformations in the Ramachandran plot.
[4]:
from os.path import exists
import molpx
from matplotlib import pylab as plt
%matplotlib ipympl
import pyemma
import numpy as np
Start from files on disk¶
[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
Featurize to Ramachandran \((\phi,\psi)\)-pairs with PyEMMA¶
[6]:
feat = pyemma.coordinates.featurizer(top)
feat.add_backbone_torsions()
src = pyemma.coordinates.source(MD_trajfiles, features=feat)
Y = src.get_output()
Visualize a FES and the molecular structures behind it¶
Execute the following cell and click either on the FES or on the slidebar
[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
Visualize trajectories, FES and molecular structures¶
[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
Paths samples along the different projections (=axis)¶
[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,
)
[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]
[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])
/home/guille/miniconda3/lib/python3.6/site-packages/ipykernel_launcher.py:4: RuntimeWarning: divide by zero encountered in log
after removing the cwd from sys.path.
Let’s do TICA and try to look a the correlations in a TICA analysis¶
[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()
01-11-17 21:01:44 pyemma.coordinates.data.featurization.featurizer.MDFeaturizer[7] WARNING The 1D arrays input for add_distances() have been sorted, and index duplicates have been eliminated.
Check the output of describe() to see the actual order of the features
[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
[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
[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,
)
[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]
[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
/home/guille/miniconda3/lib/python3.6/site-packages/ipykernel_launcher.py:3: RuntimeWarning: divide by zero encountered in log
This is separate from the ipykernel package so we can avoid doing imports until
WARNING:root:DEPRECATED: Please use 'center' method
[ ]: