1.1. Interactive exploration

In this part, we can interactively explore the structures found on the pca landscape. For this, klick the –> Live Code button above and the cells below will run automatically. Click on the PCA plot to see the corresponding 3d structures.

import matplotlib
%matplotlib widget
import numpy as np
import nglview as nv
import pandas as pd
import mdtraj as md
import sklearn as sk
import matplotlib.pyplot as plt

def pickle_load(path):
    """Load object from pickle file, given by path"""
    import pickle
    with open(path, "rb") as f:
        object = pickle.load(f)
    return object


def link_ngl_wdgt_to_ax_pos(ax, pos, ngl_widget):
    """source: https://github.com/gph82/nglview_notebooks"""
    from matplotlib.widgets import AxesWidget
    from scipy.spatial import cKDTree
    r"""
    Initial idea for this function comes from @arose, the rest is @gph82 and @clonker
    """

    kdtree = cKDTree(pos)
    assert ngl_widget.trajectory_0.n_frames == pos.shape[0]
    x, y = pos.T

    lineh = ax.axhline(ax.get_ybound()[0], c="black", ls='--')
    linev = ax.axvline(ax.get_xbound()[0], c="black", ls='--')
    dot, = ax.plot(pos[0,0],pos[0,1], 'o', c='red', ms=7)

    ngl_widget.isClick = False

    def onclick(event):
        linev.set_xdata((event.xdata, event.xdata))
        lineh.set_ydata((event.ydata, event.ydata))
        data = [event.xdata, event.ydata]
        _, index = kdtree.query(x=data, k=1)
        dot.set_xdata((x[index]))
        dot.set_ydata((y[index]))
        ngl_widget.isClick = True
        ngl_widget.frame = index

    def my_observer(change):
        r"""Here comes the code that you want to execute
        """
        ngl_widget.isClick = False
        _idx = change["new"]
        try:
            dot.set_xdata((x[_idx]))
            dot.set_ydata((y[_idx]))
        except IndexError as e:
            dot.set_xdata((x[0]))
            dot.set_ydata((y[0]))
            print("caught index error with index %s (new=%s, old=%s)" % (_idx, change["new"], change["old"]))

    # Connect axes to widget
    axes_widget = AxesWidget(ax)
    axes_widget.connect_event('button_release_event', onclick)

    # Connect widget to axes
    ngl_widget.observe(my_observer, "frame", "change")
    

folder = "test_report/assets/22/"
topology = f"{folder}mc_gas.prmtop"
trajectory = f"{folder}ed6dd3148ef9b069_short-traj.netcdf"
dihedral_file = f"{folder}ed6dd3148ef9b069_dihedrals_short.dat"
dpca_file = f"{folder}ed6dd3148ef9b069_dPCA.dat"
weight_file = f"{folder}ed6dd3148ef9b069_dPCA_weights_MC_short.dat"

traj = md.load(trajectory, top=topology)
# Plot dihedral PCA
dihedrals = pickle_load(dihedral_file)
dpca = pickle_load(dpca_file)
weights = pickle_load(weight_file)
pca_dihedrals = dpca.transform(dihedrals)

---------------------------------------------------------------------------
FileNotFoundError                         Traceback (most recent call last)
Input In [2], in <cell line: 1>()
----> 1 traj = md.load(trajectory, top=topology)
      2 # Plot dihedral PCA
      3 dihedrals = pickle_load(dihedral_file)

File /opt/hostedtoolcache/Python/3.8.13/x64/lib/python3.8/site-packages/mdtraj/core/trajectory.py:396, in load(filename_or_filenames, discard_overlapping_frames, **kwargs)
    394 topkwargs.pop("stride", None)
    395 topkwargs.pop("start", None)
--> 396 kwargs["top"] = _parse_topology(kwargs.get("top", filename_or_filenames[0]), **topkwargs)
    398 #get the right loader
    399 try:
    400     #loader = _LoaderRegistry[extension][0]

File /opt/hostedtoolcache/Python/3.8.13/x64/lib/python3.8/site-packages/mdtraj/core/trajectory.py:177, in _parse_topology(top, **kwargs)
    175     topology = _traj.topology
    176 elif isinstance(top, (string_types, os.PathLike)) and (ext in ['.prmtop', '.parm7', '.prm7']):
--> 177     topology = load_prmtop(top, **kwargs)
    178 elif isinstance(top, (string_types, os.PathLike)) and (ext in ['.psf']):
    179     topology = load_psf(top, **kwargs)

File /opt/hostedtoolcache/Python/3.8.13/x64/lib/python3.8/site-packages/mdtraj/formats/prmtop.py:123, in load_prmtop(filename, **kwargs)
    120 raw_data   = {}
    121 ignoring = False
--> 123 with open(filename, 'r') as f:
    124     for line in f:
    125         if line[0] == '%':

FileNotFoundError: [Errno 2] No such file or directory: 'test_report/assets/22/mc_gas.prmtop'

v = nv.show_mdtraj(traj)
v.clear_representations()
v.add_representation('licorice')
fig, ax = plt.subplots()
link_ngl_wdgt_to_ax_pos(ax, pca_dihedrals, v)
ax = plt.scatter(pca_dihedrals[:,0], pca_dihedrals[:,1])
v