Working with trajectory ensembles#

Welcome

Welcome to the MD section of the EncoderMap tutorial. All EncoderMap tutorials are provided as jupyter notebooks, that you can run locally, on binderhub, or even on google colab.

Run this notebook on Google Colab:

Find the documentation of EncoderMap:

https://ag-peter.github.io/encodermap

Goals:

In this tutorial you will learn:

What CVs are.
How EncoderMaps’ new SingleTraj class loads MD data.
How a SingleTraj can be associated with CVs.
How to construct a TrajEnsemble from multiple SingleTrajs.
How you can write your own Features that will put CVs into SingleTraj and TrajEnsemble.
How to save your SingleTrajs and TrajEnsembles along with their CVs.

For Google colab only:#

If you’re on Google colab, please uncomment these lines and install EncoderMap.

[1]:

# !wget https://gist.githubusercontent.com/kevinsawade/deda578a3c6f26640ae905a3557e4ed1/raw/b7403a37710cb881839186da96d4d117e50abf36/install_encodermap_google_colab.sh
# !sudo bash install_encodermap_google_colab.sh

Primer#

Collective Variables#

Collective Variables (CVs) are often used to simplify and filter the xyz-Data of MD simulations. They are often employed to make sense of a complex protein (or more general molecular) system by breaking it down into just a few well-defined descriptors. They are similar to reaction coordinates, which are 1-dimensional variables along a reaction pathway but can be higher dimensional. When we think about a receptor-ligand system, the distance between the two species is often used as a reaction coordinate. A set of this distance CV and the relative rotation between receptor and ligand can add much more information and help understand the system. It becomes apparent that clever selection of CVs is an important task that many scientists in different fields face day to day.

With the tools presented in this notebook you will be able to work in a fluent and natural way with MD trajectories and their associated CV data. In the scope of this work we define CVs as:

A collection of values that is in one of its dimensions aligned with the frames/timesteps of the underlying simulation data.

Example CVs#

Here’s a list of widely used CVs:

Distances: This category of CVs can contain a 1-diemsnional scalar value describing the distance between two species in a receptor-ligand system. A 1-dimensional CV of the end-to-end distance of a protein can describe the proteins folding state in an approximate and generalizing manner. Distance CVs can also be higher-dimensional. The pairwise distances between \(\mathrm{C_\alpha}\) atoms of a protein with \(n\) residues can be captured as a \(n \times n\) matrix. A so-called hollow matrix, where all diagonal-elements are zero which is also symmetric. Most often a vector of length \(\binom{n}{2}\) is used to describe the pairwise distances in a protein. The distance between a protein and a membrane/interface would also fall into this category.
Angles: Angular CVs lie in a periodic space of \((-\pi, \pi]\) or \((0, 2\pi]\), or \((0, 360]\). The well-known Ramachandran plot uses a protein’s \(\phi\) and \(\psi\) angles to define regions of \(psi\)-\(phi\) combinations which correlate to different secondary structure motifs. Besides the backbone angles other angles can become important in understanding a protein’s conformations. Lysine is often modified via acetylation, phosphorylation, methylation, ubiquitylation and it’s sidechain angles (\(\chi_1\) to \(\chi_5\)) can be important descriptors in such a system. Pseudo-Dihedral angles lie also in this category.
Integer/binary values. If a protein has well-defined states (folded and unfolder) a binary value describing these states could also be a useful CV.
Values from other calculations and hyperparameters. An example for this category could be the temperature at which the simulation was carried out, when simulations were conducted at multiple temperatures. If phase space sub states were obtained by either using Markov-Chain models, or by using some sort of clustering algorithm (GROMOS), the membership to such cluster could also present CVs.
Positional values: Maybe even the full position of an atom or a group of atoms could be an important descriptor for a system.
etc.

We can view the raw xyz-data of trajectories as a high-dimensional array with a shape of (n_frames, n_atoms, 3). The last axis correxponds to the cartesian coordinates. Indexing this axis will give us the x, y, and z coordinates, respectively.

[2]:

# todo: delete
from __future__ import annotations
import plotly.graph_objects as go
import numpy as np
from pathlib import Path

import encodermap as em

%load_ext autoreload
%autoreload 2

/home/kevin/git/encoder_map_private/encodermap/__init__.py:194: GPUsAreDisabledWarning: EncoderMap disables the GPU per default because most tensorflow code runs with a higher compatibility when the GPU is disabled. If you want to enable GPUs manually, set the environment variable 'ENCODERMAP_ENABLE_GPU' to 'True' before importing EncoderMap. To do this in python you can run:

import os; os.environ['ENCODERMAP_ENABLE_GPU'] = 'True'

before importing encodermap.
  _warnings.warn(

[3]:

traj = em.load_project('linear_dimers', traj=0)

print(f"First frame, first atom, x-coordinate:\n{traj.xyz[0, 0, 0]=}")
print(f"First 2 frames, first atom, y-coordinates:\n{traj.xyz[:2, 0, 1]=}")
print(f"All frames, all atoms. Only z-coordinates:\n{traj.xyz[..., -1]}")

/home/kevin/git/encoder_map_private/encodermap/kondata.py:477: UserWarning: EncoderMap's repository at https://sawade.io/encodermap_data/linear_dimers is unreachable.
  warnings.warn(f"EncoderMap's repository at {url} is unreachable.")

First frame, first atom, x-coordinate:
traj.xyz[0, 0, 0]=6.1210003
First 2 frames, first atom, y-coordinates:
traj.xyz[:2, 0, 1]=array([6.42     , 6.3940005], dtype=float32)
All frames, all atoms. Only z-coordinates:
[[1.1060001  1.092      1.028      ... 5.7790003  5.7790003  5.881     ]
 [1.056      1.0610001  0.96400005 ... 5.76       5.689      5.8650002 ]
 [1.062      1.0090001  1.144      ... 5.8500004  5.9550004  5.8240004 ]
 ...
 [0.12100001 0.031      0.12400001 ... 4.783      4.846      4.7850003 ]
 [0.116      0.081      0.125      ... 4.617      4.6980004  4.5950003 ]
 [0.06900001 0.06       0.063      ... 4.5290003  4.486      4.6260004 ]]

The 3-dimensions of this data can be used for plotting. However, we don’t gain much from these plots.

[4]:

em.plot.plot_raw_data(
    traj,
    frame_slice=slice(0, 5001, 1000),
    atom_slice=slice(0, 1500, 150),
)

Adding more context to the raw xyz-data (like size, color, bonds) we can better understand the data.

[5]:

em.plot.plot_ball_and_stick(traj)

/home/kevin/git/encoder_map_private/encodermap/trajinfo/info_single.py:1691: UserWarning:

Dropping CVs from trajectory. Slicing CVs with this method is currently not possible. Raise an issue if you want to have this feature added.

We can get more information by using additional libraries like NGLView which also render secondary structure. However, we can’t really understand how the protein moves from one folding state to another. Hit the “Play” button and have a look at the first 200 frames of the simulation. Can you spot the LEU56-SER57-ASP58-TYR59 residues transition from a turn into an \(\alpha\)-helix?

[6]:

import nglview as nv
nv.show_mdtraj(traj.traj[::100])

That’s when collective variables can come in handy. They break down and combine different xyz data to help us understand the grand picture of protein folding.

A Ramachandran plot is a very general description of secondary structure motifs. The test protein does indeed contain \(alpha\)-helices and \(beta\)-sheets. We can’t really understand the conformation of the protein.

[7]:

em.plot.plot_ramachandran(traj, subsample=1000)

A DSSP plot doesn’t condense the \(phi\) and \(psi\) angles into a single plot. Here we can assign secondary structure motifs to certain residues at certain times of the simulation.

[8]:

em.plot.plot_dssp(
    traj,
    simplified=False,
    subsample=slice(0, 500),
)

Distances can also be used as useful collective variables. The end-to-end distance can give helpful insights into a proteins general conformation (collapsed vs. extended).

[9]:

em.plot.plot_end2end(traj)

The new `SingleTraj` class#

The SingleTraj class is meant as a container to hold a trajectory’s xyz coordinates, its topology, and its CVs. This class builds the backbone of the TrajEnsemble class which will be looked at later.

In the background of the SingleTraj class the MDTraj package (mdtraj/mdtraj) works its magic, however this class offers more benefits:

The SingleTraj class keeps track of CVs while indexing and slicing.
The trajectories’ xyz data is only loaded when it is actually needed. Most manipuilation can thus be done before creating huge objects in memory.

The SingleTraj class can be initialized in many ways:

From a trajectory file (.xtc, .dcd, .lammpstrj) and a topology file (.gro, .pdb).

[13]:

output_dir = Path(output_dir)
traj_from_xtc_file = em.SingleTraj(output_dir / "01.xtc", top=output_dir / "01.pdb")
print(traj_from_xtc_file)

encodermap.SingleTraj object. Currently not in memory. Basename is '01'. Not containing any CVs.

From an url of the pdb database.

[14]:

traj_from_url = em.SingleTraj('https://files.rcsb.org/view/1YUF.pdb')
print(traj_from_url)

encodermap.SingleTraj object. Currently not in memory. Basename is '1YUF'. Not containing any CVs.

Or by just by providing a pdb code to the from_pdb_id() constructor of the SingleTraj class.

[15]:

traj_from_pdb_id = em.SingleTraj.from_pdb_id('1YUG')
print(traj_from_pdb_id)

encodermap.SingleTraj object. Currently not in memory. Basename is '1YUG'. Not containing any CVs. Common string is '1YUG'.

You can also directly load the trajectory from an EncoderMap project. In this case, the trajectory already has CVs.

[16]:

# with the load_project() method
traj = em.load_project('linear_dimers', traj=0)
print(traj)

encodermap.SingleTraj object. Currently not in memory. Basename is 'trajs'. CV central_angles with shape (1, 5001, 454) loaded. CV central_cartesians with shape (1, 5001, 456, 3) loaded. CV central_dihedrals with shape (1, 5001, 453) loaded. CV central_distances with shape (1, 5001, 455) loaded. CV lowd with shape (1, 5001, 2) loaded. CV side_dihedrals with shape (1, 5001, 322) loaded.

``SingleTraj``s keep track of the files they have been instantiated from.

[17]:

print(
    f"basename  = {traj_from_xtc_file.basename:<70}{traj_from_url.basename:<70}\n"
    f"traj_file = {traj_from_xtc_file.traj_file:<70}{traj_from_url.traj_file:<70}\n"
    f"top_file  = {traj_from_xtc_file.top_file:<70}{traj_from_url.top_file:<70}\n"
    f"extension = {traj_from_xtc_file.extension:<70}{traj_from_url.extension:<70}\n"
    f"n_frames  = {traj_from_xtc_file.n_frames:<70}{traj_from_url.n_frames:<70}\n"
    f"n_atoms   = {traj_from_xtc_file.n_atoms:<70}{traj_from_url.n_atoms:<70}"
)

basename  = 01                                                                    1YUF
traj_file = /home/kevin/git/encoder_map_private/tests/data/linear_dimers/01.xtc   https://files.rcsb.org/view/1YUF.pdb
top_file  = /home/kevin/git/encoder_map_private/tests/data/linear_dimers/01.pdb   https://files.rcsb.org/view/1YUF.pdb
extension = .xtc                                                                  .pdb
n_frames  = 5001                                                                  16
n_atoms   = 1521                                                                  720

Difference between ``traj``, ``trajectory`` and ``top``, ``topology``

traj and top always give mdtraj.Trajectory and mdtraj.Topology, respectively. They are loaded on demand and return the corresponding mdtraj object.

trajectory and topology can be False and represent the current backend of the SingleTraj object.

When instantiated, SingleTrajd do not load the data from disk. Only when requested.

[18]:

print(traj_from_xtc_file.topology)
print(traj_from_xtc_file.top)
print(traj_from_xtc_file.topology)

False
<mdtraj.Topology with 1 chains, 152 residues, 1521 atoms, 1534 bonds>
False

[19]:

print(traj_from_xtc_file.trajectory)
print(traj_from_xtc_file.traj)
print(traj_from_xtc_file.trajectory)

False
<mdtraj.Trajectory with 5001 frames, 1521 atoms, 152 residues, and unitcells>
False

Loading can be forced

Using the load() method, the trajectory will be loaded. From that point forwards, the xyz data is kept in RAM and can be accessed. The unload() function does the reverse and frees up the RAM, but the xyz data can be loaded again. If your RAM is large enough you would not need the unload() function, but it is there nonetheless.

[20]:

traj_from_xtc_file.load_traj()
print(traj_from_xtc_file.topology)
traj_from_xtc_file.unload()
print(traj_from_xtc_file.topology)

<mdtraj.Topology with 1 chains, 152 residues, 1521 atoms, 1534 bonds>
False

Inside a context manager, the SingleTraj is always loaded and unloaded afterwards.

[21]:

with traj_from_xtc_file as t:
    print(t.topology)
print(traj_from_xtc_file.topology)

<mdtraj.Topology with 1 chains, 152 residues, 1521 atoms, 1534 bonds>
False

If nglview is set up, you can take a look at the trajectory with this code:

[22]:

view = traj_from_xtc_file.show_traj()
view

Some methods and attributes are duplicated from mdtraj. This allows us to call some mdtraj functions on the SingleTraj object directly.

[23]:

selection = traj_from_xtc_file.select('name CA')
print(selection[:5])
dssp = md.compute_dssp(traj_from_xtc_file)
print(dssp[0, :5].tolist())

[ 4 13 25 34 51]
['C', 'E', 'E', 'E', 'E']

[24]:

md.compute_center_of_mass(traj_from_xtc_file)[0]

[24]:

array([7.60928699, 7.31795958, 3.59687811])

By indexing a SingleTraj, you get another instance of SingleTraj containing only one frame.

[25]:

frame = traj_from_xtc_file[0]
print(frame)

encodermap.SingleTraj object. Data currently in memory. Basename is '01'. At indices (0,). Not containing any CVs.

The indexing of SingleTrajs is postponed, meaning, that for certain indexing types (integers, sequences of integers, slices without the step argument) the trajectory data is not loaded from disk. Only, when loaded, are the indexes applied one after the other.

This cell is fast. No file operations are carried out here.

[26]:

traj_from_xtc_file.unload()
frame = traj_from_xtc_file[:10][[0, 1]]
print(frame)
print(frame.trajectory)
print(frame.n_frames)

encodermap.SingleTraj object. Currently not in memory. Basename is '01'. At indices (slice(None, 10, None), [0, 1]). Not containing any CVs.
False
2

This cell is a bit slower. Here, the trajectory is loaded from file.

[27]:

frame.load_traj()
print(frame)

encodermap.SingleTraj object. Data currently in memory. Basename is '01'. At indices (slice(None, 10, None), [0, 1]). Not containing any CVs.

If the index is larger than the number of frames, the SingleTraj class will throw an exception.

[28]:

traj_from_xtc_file.unload()
frame = traj_from_xtc_file[10_000]
print(frame)

---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
Cell In[28], line 2
      1 traj_from_xtc_file.unload()
----> 2 frame = traj_from_xtc_file[10_000]
      3 print(frame)

File ~/git/encoder_map_private/encodermap/trajinfo/info_single.py:1986, in SingleTraj.__getitem__(self, key)
   1984 if isinstance(key, (int, np.integer)):
   1985     if key > self.n_frames:
-> 1986         raise IndexError(
   1987             f"Index {key} out of range for traj with "
   1988             f"{self.n_frames} frames."
   1989         )
   1990 if isinstance(key, (list, np.ndarray)):
   1991     if any([k > self.n_frames for k in key]):

IndexError: Index 10000 out of range for traj with 5001 frames.

SingleTrajs keep the number of their frames organized. If a SingleTraj is loaded from disk the frames are represented by a list of integers:

[29]:

traj_from_xtc_file.id

[29]:

array([   0,    1,    2, ..., 4998, 4999, 5000])

ndexing these frames with a [::2] slice will give the frames [0, 2, 4, ..., n].

[30]:

traj_from_xtc_file[::2].id

[30]:

array([   0,    2,    4, ..., 4996, 4998, 5000])

This begs the question what happens if we index these frames again? Which frame numbers will be chosen

The answer is: There are two methods.

The normal indexing [[0, 1, 2]] does not care for the original frame indices. It will select whichever frame is at positions 0, 1, and 2, regardless of what their original id was.

[31]:

traj_from_xtc_file[::2][[0, 1, 2, 3]].id

[31]:

array([0, 2, 4, 6])

[32]:

traj_from_xtc_file[::5][[0, 1, 2, 5, 6]].id

[32]:

array([ 0,  5, 10, 25, 30])

The other method uses the fsel[] selector. It expects the original frame indices.

[33]:

traj_from_xtc_file[::2].fsel[[0, 2, 4, 6]].id

[33]:

array([0, 2, 4, 6])

[34]:

traj_from_xtc_file[::5].fsel[[0, 5, 10, 25, 30]].id

[34]:

array([ 0,  5, 10, 25, 30])

This also means, that the fsel selector can fail for frames that are not present in a SingleTraj.

This means, that some frames are not accessible, when we do some weird slicing:

[35]:

traj_from_xtc_file[::5].fsel[3]

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[35], line 1
----> 1 traj_from_xtc_file[::5].fsel[3]

File ~/git/encoder_map_private/encodermap/trajinfo/info_single.py:199, in SingleTrajFsel.__getitem__(self, item)
    194         raise ValueError(
    195             f"The `fsel[]` method of `SingleTraj` takes {CanBeIndex} types, "
    196             f"but {type(item)} was provided."
    197         )
    198 if len(idx) == 0:
--> 199     raise ValueError(
    200         f"No frames with frame index {item} in trajectory {self.other} "
    201         f"with frames: {self.other._frames}"
    202     )
    203 return self.other[idx]

ValueError: No frames with frame index 3 in trajectory encodermap.SingleTraj object. Currently not in memory. Basename is '01'. At indices (slice(None, None, 5),). Not containing any CVs. with frames: [   0    5   10 ... 4990 4995 5000]

You can also give a numpy array, a list or even a slice into the indexing.

Indexing with without a stop ([::5]) will put the trajectory into memory.

[36]:

traj_from_xtc_file.unload()
subsample = traj_from_xtc_file[::2]
print(traj_from_xtc_file.n_frames)
print(subsample)
print(subsample.n_frames)

5001
encodermap.SingleTraj object. Currently not in memory. Basename is '01'. At indices (slice(None, None, 2),). Not containing any CVs.
2501

[37]:

traj_from_xtc_file.unload()
subsample = traj_from_xtc_file[[0, 1, 5, 6]]
print(subsample)
print(subsample.n_frames)

encodermap.SingleTraj object. Currently not in memory. Basename is '01'. At indices ([0, 1, 5, 6],). Not containing any CVs.
4

Note: Andvanced slicing can result in trajectories with 0 frames in them, or possibly reverse the time axis. Use this feature only if you are sure about what you are doing.

[38]:

traj_from_xtc_file.unload()
subsample = traj_from_xtc_file[5:46:3]
print(subsample)
print(subsample.n_frames)

encodermap.SingleTraj object. Currently not in memory. Basename is '01'. At indices (slice(5, 46, 3),). Not containing any CVs.
14

The HDF5 file format (ending wiht the .h5 extension) allows us to directly extract frames and accelerate loading. We can save a .h5 formatted file from an SingleTraj class by calling:

[39]:

with tempfile.TemporaryDirectory() as d:
    d = Path(d)
    traj_from_xtc_file.save(d / 'my_traj.h5', overwrite=True)

Loading of .h5 files is similar to all files:

[40]:

with tempfile.TemporaryDirectory() as d:
    d = Path(d)
    traj_from_xtc_file.save(d / 'my_traj.h5', overwrite=True)
    traj_h5 = em.SingleTraj(d / 'my_traj.h5')

    print(traj_h5[0])

encodermap.SingleTraj object. Currently not in memory. Basename is 'my_traj'. At indices (0,). Not containing any CVs.

The SingleTraj class can be used as an iterator in two ways:

Iterate over the frames with for frame in traj.
Iterate over frame_number and frame with for frame_no, frame in traj.iterframes().

[41]:

for frame in traj[::2][:5]:
    print(frame.index)

(None, slice(None, None, 2), slice(None, 5, None), 0)
(None, slice(None, None, 2), slice(None, 5, None), 1)
(None, slice(None, None, 2), slice(None, 5, None), 2)
(None, slice(None, None, 2), slice(None, 5, None), 3)
(None, slice(None, None, 2), slice(None, 5, None), 4)

[42]:

for frame_num, frame in traj[::2][:5].iterframes():
    print(frame_num)

The `TrajEnsemble` class contains multiple `SingleTraj`s#

The trajectory ensemble is everything you’ve always wanted, and more. Really, it is. Trajectory ensembles unlock fundamental ideas in statistical mechanics, including connections between equilibrium and non-equilibrium phenomena.

Daniel M. Zuckerman (https://statisticalbiophysicsblog.org/?p=92#more-92)

EncoderMap tries to implement the idea of a trajectory ensemble with the TrajectoryEnsemble class. A container for multiple SingleTrajs.

A trajectory ensemble can be created by providing it multiple trajectory files:

[43]:

output_dir = Path(em.get_from_kondata("pASP_pGLU", mk_parentdir=True, silence_overwrite_message=True))

trajs_from_files = em.load(
    [
        output_dir / "asp7.xtc",
        output_dir / "glu7.xtc",
    ],
    [
        output_dir / "asp7.pdb",
        output_dir / "glu7.pdb",
    ],
    common_str=["asp7", "glu7"]
)

/home/kevin/git/encoder_map_private/encodermap/kondata.py:477: UserWarning:

EncoderMap's repository at https://sawade.io/encodermap_data/pASP_pGLU is unreachable.

[44]:

trajs_from_files

[44]:

<encodermap.TrajEnsemble object. Current backend is no_load. Containing 2 trajectories. Common str is ['asp7', 'glu7']. Not containing any CVs. Object at 0x7055c7b2fa60>

The common_str argument provides a means to group similar trajectories. It searches for matching substrings in the filename of the provided trajectories and topologies. In the trajs_from_files, there are only one TrajEnsemble with a single trajectory per common_str.

[45]:

for common_str, trajs_subset in trajs_from_files.trajs_by_common_str.items():
    print(f"The common_str '{common_str}' contains a subset of {trajs_subset.n_trajs} trajectories.")

The common_str 'asp7' contains a subset of 1 trajectories.
The common_str 'glu7' contains a subset of 1 trajectories.

Something that EncoderMap does different than packages like MDTraj or MDAnalysis is that trajectories can be grouped together even when their topologies are different.

[46]:

for top, trajs_subset in trajs_from_files.trajs_by_top.items():
    print(f"The topology {top} has {trajs_subset.n_trajs} in the ensemble.")

The topology <mdtraj.Topology with 1 chains, 7 residues, 73 atoms, 71 bonds> has 1 in the ensemble.
The topology <mdtraj.Topology with 1 chains, 7 residues, 80 atoms, 78 bonds> has 1 in the ensemble.

The trajs_from_files instance has trajectories with 2 different topologies in it. This feature is meant to represent the expression of mutations in biological systems.

https://evolution.berkeley.edu/dna-and-mutations/types-of-mutations/

Indexing TrajEnsembles can yield different output types. Indexing with a sinlge int will yield the corresponding trajectory.

[47]:

trajs_from_files[1]

[47]:

<encodermap.SingleTraj object. Currently not in memory. Basename is 'glu7'. Not containing any CVs. Common string is 'glu7'. Object at 0x7055c895a920>

Python lists, NumPy arrays and slices can also be used for indexing. They all return TrajEnsemble object and raise Errors, when they are indexing trajectories that do not exist in the ensemble.

[48]:

trajs_from_files[[0, 1]]

[48]:

<encodermap.TrajEnsemble object. Current backend is no_load. Containing 2 trajectories. Common str is ['glu7', 'asp7']. Not containing any CVs. Object at 0x7055ce6b1a20>

[49]:

trajs_from_files[[0, 1, 2]]

---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
Cell In[49], line 1
----> 1 trajs_from_files[[0, 1, 2]]

File ~/git/encoder_map_private/encodermap/trajinfo/info_all.py:3051, in TrajEnsemble.__getitem__(self, key)
   3049     return self.trajs[key]
   3050 elif isinstance(key, list) and not isinstance(key[0], list):
-> 3051     new_class = self._return_trajs_by_index(key)
   3052     return new_class
   3053 elif isinstance(key, np.ndarray):

File ~/git/encoder_map_private/encodermap/trajinfo/info_all.py:2566, in TrajEnsemble._return_trajs_by_index(self, index)
   2564         trajs_subset = self.trajs[ind]._gen_ensemble()
   2565     else:
-> 2566         new_traj = self.trajs[ind]._gen_ensemble()
   2567         trajs_subset += new_traj
   2568 trajs_subset.common_str = new_common_str

IndexError: list index out of range

Slices can contain start, stop, step or any combination, thereof:

[50]:

trajs_from_files[0:1]

[50]:

<encodermap.TrajEnsemble object. Current backend is no_load. Containing 1 trajectories. Common str is ['asp7']. Not containing any CVs. Object at 0x7055cdfd58a0>

[51]:

trajs_from_files[::3]

[51]:

<encodermap.TrajEnsemble object. Current backend is no_load. Containing 1 trajectories. Common str is ['asp7']. Not containing any CVs. Object at 0x7055c93a5db0>

They will also raise Errors, if trajectories are not accessible.

[52]:

trajs_from_files[5:]

---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
Cell In[52], line 1
----> 1 trajs_from_files[5:]

File ~/git/encoder_map_private/encodermap/trajinfo/info_all.py:3069, in TrajEnsemble.__getitem__(self, key)
   3067     start, stop, step = key.indices(self.n_trajs)
   3068     list_ = list(range(start, stop, step))
-> 3069     new_class = self[list_]
   3070     return new_class
   3071 elif isinstance(key, list) and all(isinstance(k, list) for k in key):

File ~/git/encoder_map_private/encodermap/trajinfo/info_all.py:3050, in TrajEnsemble.__getitem__(self, key)
   3048 if isinstance(key, (int, np.int32, np.int64)):
   3049     return self.trajs[key]
-> 3050 elif isinstance(key, list) and not isinstance(key[0], list):
   3051     new_class = self._return_trajs_by_index(key)
   3052     return new_class

IndexError: list index out of range

Similar to trajectories, which can contain frames with various frame-numbers, TrajEnsemble can contain trajectories with various trajectory-numbers. The regular indexing [] does not use the trajectory number.

[53]:

trajs = trajs_from_files.copy()
trajs[0].traj_num = 10
trajs[1].traj_num = 20
trajs[0].traj_num

[53]:

But the tsel[] selector uses the trajectory-number:

[54]:

trajs.tsel[20].traj_num

[54]:

And raises Errors, if a trajectory with that index is not present.

[55]:

trajs.tsel[0]

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[55], line 1
----> 1 trajs.tsel[0]

File ~/git/encoder_map_private/encodermap/trajinfo/info_all.py:765, in TrajEnsembleTsel.__getitem__(self, item)
    763 if isinstance(item, (int, np.int64)):
    764     if item not in items:
--> 765         raise ValueError(
    766             f"No trajectories with traj_num {item} in TrajEnsemble {self.other} "
    767             f"with trajectories: {items}"
    768         )
    769     return self.other.trajs_by_traj_num[item]
    770 elif isinstance(item, (list, np.ndarray)):

ValueError: No trajectories with traj_num 0 in TrajEnsemble encodermap.TrajEnsemble object. Current backend is no_load. Containing 2 trajectories. Common str is ['asp7', 'glu7']. Not containing any CVs. with trajectories: [10 20]

2-dimensional numpy arrays can be used for indexing TrajEnsembles, by providing trajectory-numbers and frame-numbers arranged in the rows of the array. This type of selection does not respect the actual trajectory-numbers. That’s what - you guessed it - tsel[] is for.

[56]:

index = np.array([
    [0, 0],
    [0, 1],
    [0, 2],
    [0, 3],
    [1, 0],
    [1, 2],
    [1, 4],
])
trajs[index].id

[56]:

array([[10,  0],
       [10,  1],
       [10,  2],
       [10,  3],
       [20,  0],
       [20,  2],
       [20,  4]])

The tsel[] selector does the same, but with trajectory-numbers and frame-numbers.

[57]:

index = np.array([
    [10, 0],
    [10, 10],
    [10, 20],
    [10, 30],
    [20, 0],
    [20, 2],
    [20, 4],
])
trajs.tsel[index].id

[57]:

array([[10,  0],
       [10, 10],
       [10, 20],
       [10, 30],
       [20,  0],
       [20,  2],
       [20,  4]])

Think about the TrajectoryEnsemble as a high-dimensional object, that contains an extra axis. The trajectory axis (The SingleTraj object only has an atom and a time axis).

[58]:

trajs[0] + trajs[1]

[58]:

<encodermap.TrajEnsemble object. Current backend is no_load. Containing 2 trajectories. Common str is ['asp7', 'glu7']. Not containing any CVs. Object at 0x7055c98a69e0>

As we have previously seen, there are some out-of-the box TrajectoryEnsembles that come with EncoderMap. You can get a list of all available projects via:

[59]:

em.ALL_PROJECT_NAMES

[59]:

typing.Union[typing.Literal['linear_dimers'], typing.Literal['pASP_pGLU'], typing.Literal['Ub_K11_mutants'], typing.Literal['cube'], typing.Literal['1am7'], typing.Literal['H1Ub']]

Projects are automatically loaded from either:

The University of Konstanz’ data repository KonDATA
A web-based repository with sample data accessible under: https://134.34.112.158/encodermap_data

[60]:

trajs = em.load_project("pASP_pGLU")
trajs

[60]:

<encodermap.TrajEnsemble object. Current backend is no_load. Containing 7 trajectories. Common str is ['asp10', 'glu8', 'glu6', 'glu7', 'asp8', 'asp6', 'asp7']. CV central_dihedrals with shape (35007, 27) loaded. CV central_angles with shape (35007, 28) loaded. CV side_dihedrals with shape (35007, 28) loaded. CV central_distances with shape (35007, 29) loaded. CV central_cartesians with shape (35007, 30, 3) loaded. Object at 0x7055cdfd5870>

Loading CVs#

After learning about the basics of the SingleTraj and TrajEnsemble class we will come back to collective variables. EncoderMap offers a way of “adding” CV data to your trajectories and treat them as a single instance. A trajectory with CV data. The easiest way to load CV data is to provide an already existing NumPy array. However, you will be asked to also provide the attribute name (attr_name) of the array. With this you could load multiple CV datasets, that differ in ther attribute names. Here’s an example:

Slicing with CVs.#

Slicing keeps your values where they should be.

[67]:

index = np.where(np.array(raman_labels) == 'psi 50')[0]
print(traj[0].raman[0, index])
print(traj[[0, 5, 10, 20]].raman[:,index])

[0.3531243]
[[0.3531243 ]
 [4.37541582]
 [4.62936691]
 [5.31620961]]

Loading from files#

CVs can be loaded by providing a string to files. First, let us save some files.

[68]:

# save numpy
np.save('raman_file.npy', traj.CVs["raman"])

# save text
np.savetxt('cluster_membership_file.txt', traj.cluster_membership)

# save full CV dataset as NetCDF
traj._CVs.to_netcdf('full_CV_dataset.nc', engine="h5netcdf", )

If not providing an attr_name, while loading files, the filename will be used:

[69]:

traj.load_CV('raman_file.npy')
traj.load_CV('cluster_membership_file.txt')

[70]:

print(traj.CVs.keys())

dict_keys(['raman', 'cluster_membership', 'raman_file', 'cluster_membership_file'])

Multiple CVs can be reconstructed from xarray NetCDF files (most end with .nc). If there are conflicts the new data from disk will overwrite the old.

[71]:

traj = em.SingleTraj(traj.traj_file, traj.top_file)
print(traj.CVs.keys())
traj.load_CV('full_CV_dataset.nc')
print(traj.CVs.keys())

dict_keys([])
dict_keys(['raman', 'cluster_membership'])

Wrtiting custom features No 1#

Writing your custom features can be done by subclassing pyemma’s features. Required methods and attributes to make your feature work are:

The methods __init__, describe, and call.
The instance attribute dimension, which defines the shape of the returned array.

If you want to change the name of the feature, as it appears in the xarray.Dataset you can set the attribute name.

In the next cell we will define a Feature that provides a random float for the atoms in a trajectory.

[79]:

class RandomFloatForAtomFeature(em.features.CustomFeature):
    # write an __init__
    def __init__(self, traj, selstr='all'):
        """Init of RandomFloatForAtomFeature.

        Args:
            top (mdtraj.Topology): The topology to select atoms from.

        Keyword Args:
            selstr (str, optional): The string to provide to top.select().
            Defaults to 'all'.

        """
        self.traj = traj
        self.top = traj.top
        self.indexes = self.traj.top.select(selstr)

        # Here, we set a required attribute
        # The dimension is needed and should
        # be an integer value of the length of the feature
        # in our case, that's how many atoms are in `self.indexes`
        self.dimension = len(self.indexes)

    def describe(self):
        """This method is not allowed to take any arguments.

        This method provides labels.

        Returns:
            list: A lsit of str, each str describing one feature.

        """
        # In this method we will build a list of str
        # Each str should describe one of our features
        # We assign floats to atoms, so the labels should tell something about the atoms
        getlbl = (
            lambda at: f"atom {at.name:>4}:{at.index:5} {at.residue.name}:{at.residue.resSeq:>4}"
        )
        labels = []
        for i in self.indexes:
            i = self.traj.top.atom(i)
            labels.append(f"Random float for {getlbl(i)}")
        return labels

    def call(self, traj):
        """This method provides values.

        Args:
            traj (mdtraj.Trajectory): An mdtraj.Trajectory.

        Returns:
            np.ndarray: The values of the features defined in describe.

        """
        # Make sure that the returned array has correct shape
        # It is a good idea, that this array has the same length as
        # the trajectory has frames
        values = traj.xyz[..., 0]
        for i in self.indexes:
            values[:, i] = float(str(hash(str(self.traj.top.atom(i))))[-5:])
        return values

    @property
    def name(self):
        # define the name of the feature to appear in `SingleTraj._CVs`
        return 'MyAwesomeFeature'

[80]:

traj = em.SingleTraj(traj.traj_file, traj.top_file)
print(traj)
featurizer = em.Featurizer(traj)
feat = RandomFloatForAtomFeature(traj)
for i in feat.describe()[:200:25]:
    print(i)

encodermap.SingleTraj object. Currently not in memory. Basename is '01'. Not containing any CVs.
Random float for atom    N:    0 MET:   1
Random float for atom   CA:   25 ILE:   3
Random float for atom    H:   50 VAL:   5
Random float for atom  HG1:   75 THR:   7
Random float for atom    C:  100 GLY:  10
Random float for atom    H:  125 ILE:  13
Random float for atom    O:  150 LEU:  15
Random float for atom  OE1:  175 GLU:  18

[81]:

featurizer.add_custom_feature(feat)

[82]:

traj.load_CV(featurizer)

[83]:

traj._CVs.coords['MYAWESOMEFEATURE'].values

[83]:

array(['Random float for atom    N:    0 MET:   1',
       'Random float for atom    H:    1 MET:   1',
       'Random float for atom   H2:    2 MET:   1', ...,
       'Random float for atom    C: 1518 GLY: 152',
       'Random float for atom    O: 1519 GLY: 152',
       'Random float for atom  OXT: 1520 GLY: 152'], dtype='<U41')

Working with trajectory ensembles#

For Google colab only:#

Primer#

Collective Variables#

Example CVs#

The new `SingleTraj` class#

The `TrajEnsemble` class contains multiple `SingleTraj`s#

Loading CVs#

From numpy#

Slicing with CVs.#

Loading from files#

Loading with EncoderMap’s featurzier#

Wrtiting custom features No 1#

Writing custom features No 2#

Saving trajectory and CVs into one file#

Saving a complete trajectory Ensemble#

Conclusion#

This Page

Working with trajectory ensembles#

For Google colab only:#

Primer#

Collective Variables#

Example CVs#

Sharing MD data#

Classes for working with MD data#

The new SingleTraj class#

The TrajEnsemble class contains multiple SingleTrajs#

Loading CVs#

From numpy#

Slicing with CVs.#

Loading from files#

Loading with EncoderMap’s featurzier#

Wrtiting custom features No 1#

Writing custom features No 2#

Saving trajectory and CVs into one file#

Saving a complete trajectory Ensemble#

Conclusion#

This Page

The new `SingleTraj` class#

The `TrajEnsemble` class contains multiple `SingleTraj`s#