Page on trajectories and on-the-fly transformations

doc/source/conf.py

@@ -98,10 +98,8 @@ def sort_authors(filename):
 # List of patterns, relative to source directory, that match files and
 # directories to ignore when looking for source files.
 # This pattern also affects html_static_path and html_extra_path.
-
 exclude_patterns = ['_build', 'Thumbs.db', '.DS_Store',
-                    '.ipynb_checkpoints', '**/.ipynb_checkpoints', 'scripts', '.*.ipynb']
-
+                    '.ipynb_checkpoints', '**/.ipynb_checkpoints', 'scripts', '.*.ipynb', '.*']
 
 # -- Options for HTML output -------------------------------------------------
 
@@ -191,11 +189,10 @@ def sort_authors(filename):
     <script src='http://cdnjs.cloudflare.com/ajax/libs/require.js/2.1.10/require.min.js'></script>
     <script>require=requirejs;</script>
 
-
 """
 
 # substitutions
-MDAnalysis_version = '0.20.1'
+MDAnalysis_version = '1.0.0'
 
 # rst-epilog implements substitutions
 rst_epilog = """

doc/source/groups_of_atoms.rst

@@ -21,10 +21,11 @@ A :class:`~MDAnalysis.core.groups.Residue` is composed of :class:`~MDAnalysis.co
 The corresponding container groups are :class:`~MDAnalysis.core.groups.ResidueGroup` and :class:`~MDAnalysis.core.groups.SegmentGroup`. These have similar properties and available methods as :class:`~MDAnalysis.core.groups.AtomGroup`.
 
 .. ipython:: python
+    :okwarning:
 
     import MDAnalysis as mda
-    from MDAnalysis.tests.datafiles import TPR2019B3
-    u = mda.Universe(TPR2019B3)
+    from MDAnalysis.tests.datafiles import TPR, XTC
+    u = mda.Universe(TPR, XTC)
     ag = u.atoms.select_atoms('resname ARG and name CA')
     ag
 

doc/source/index.rst

@@ -84,7 +84,15 @@ Wherever possible, do not take these conversations to private channels, includin
    groups_of_atoms
    selections
    topology_system
-   trajectories
+
+.. toctree::
+   :maxdepth: 1
+   :caption: Trajectories
+   :hidden:
+
+   trajectories/trajectories
+   trajectories/slicing_trajectories
+   trajectories/transformations
    units
 
 .. toctree::

doc/source/topology_system.rst

@@ -88,6 +88,7 @@ Each of the attributes above can be added to a Universe if it was not available
     * :code:`values` (optional) : if :code:`topologyattr` is a string, the values for that attribute. This can be :code:`None` if the attribute has default values defined, e.g. :code:`bfactors`. 
 
 .. ipython:: python
+    :okwarning:
 
     import MDAnalysis as mda
     from MDAnalysis.tests.datafiles import PSF

doc/source/trajectories/slicing_trajectories.rst

@@ -0,0 +1,69 @@
+.. -*- coding: utf-8 -*-
+.. _slicing-trajectories:
+
+====================
+Slicing trajectories
+====================
+
+MDAnalysis trajectories can be indexed to return a :class:`~MDAnalysis.coordinates.base.Timestep`, or sliced to give a :class:`~MDAnalysis.coordinates.base.FrameIterator`. 
+
+.. ipython:: python
+
+    import MDAnalysis as mda
+    from MDAnalysis.tests.datafiles import PSF, DCD
+
+    u = mda.Universe(PSF, DCD)
+    u.trajectory[4]
+
+
+Indexing a trajectory shifts the :class:`~MDAnalysis.core.universe.Universe` to point towards that particular frame, updating dynamic data such as ``Universe.atoms.positions``. 
+
+.. note::
+
+    The trajectory frame is not read from the MD data. It is the internal index assigned by MDAnalysis.
+
+.. ipython:: python
+
+    u.trajectory.frame
+
+*Creating* a :class:`~MDAnalysis.coordinates.base.FrameIterator` by slicing a trajectory does not shift the :class:`~MDAnalysis.core.universe.Universe` to a new frame, but *iterating* over the sliced trajectory will rewind the trajectory back to the first frame.
+
+.. ipython::: python
+
+    fiter = u.trajectory[10::2]
+    u.trajectory.frame
+
+.. ipython:: python
+
+    fiter = u.trajectory[10::10]
+    frames = [ts.frame for ts in fiter]
+    print(frames, u.trajectory.frame)
+
+You can also create a sliced trajectory with boolean indexing and fancy indexing. Boolean indexing allows you to select only frames that meet a certain condition, by passing a :class:`~numpy.ndarray` with the same length as the original trajectory. Only frames that have a boolean value of ``True`` will be in the resulting :class:`~MDAnalysis.coordinates.base.FrameIterator`. For example, to select only the frames of the trajectory with an RMSD under 2 angstrom:
+
+.. ipython:: python
+
+    from MDAnalysis.analysis import rms
+
+    protein = u.select_atoms('protein')
+    rmsd = rms.RMSD(protein, protein).run()
+    bools = rmsd.rmsd.T[-1] < 2
+    print(bools)
+
+.. ipython:: python
+
+    fiter = u.trajectory[bools]
+    print([ts.frame for ts in fiter])
+
+You can also use fancy indexing to control the order of specific frames.
+
+.. ipython:: python
+
+    indices = [10, 2, 3, 9, 4, 55, 2]
+    print([ts.frame for ts in u.trajectory[indices]])
+
+You can even slice a :class:`~MDAnalysis.coordinates.base.FrameIterator` to create a new :class:`~MDAnalysis.coordinates.base.FrameIterator`.
+
+.. ipython:: python
+
+    print([ts.frame for ts in fiter[::3]])

doc/source/trajectories/trajectories.rst

@@ -0,0 +1,73 @@
+.. -*- coding: utf-8 -*-
+.. _trajectories:
+
+============
+Trajectories
+============
+
+In MDAnalysis, static data is contained in your universe Topology, while dynamic data is drawn from its trajectory at ``Universe.trajectory``. This is typically loaded from a trajectory file and includes information such as:
+
+    * atom coordinates (``Universe.atoms.positions``)
+    * box size (``Universe.dimensions``)
+    * velocities and forces (if your file format contains the data) (``Universe.atoms.velocities``)
+
+Although these properties look static, they are actually dynamic, and the data contained within can change.
+In order to remain memory-efficient, MDAnalysis does not load every frame of your trajectory into memory at once. Instead, a Universe has a state: the particular timestep that it is currently associated with in the trajectory. When the timestep changes, the data in the properties above shifts accordingly.
+
+The typical way to change a timestep is to index it. ``Universe.trajectory`` can be thought of as a list of :class:`~MDAnalysis.coordinates.base.Timestep`\ s, a data structure that holds information for the current time frame. For example, you can query its length.
+
+.. ipython:: python
+
+    import MDAnalysis as mda
+    from MDAnalysis.tests.datafiles import PSF, DCD
+
+    u = mda.Universe(PSF, DCD)
+    len(u.trajectory)
+
+When a trajectory is first loaded from a file, it is set to the first frame (with index 0), by default.
+
+.. ipython:: python
+
+    print(u.trajectory.ts, u.trajectory.time)
+
+Indexing the trajectory returns the timestep for that frame, and sets the Universe to point to that frame until the timestep next changes.
+
+.. ipython:: python
+
+    u.trajectory[3]
+
+.. ipython:: python
+
+    print('Time of fourth frame', u.trajectory.time)
+
+Many tasks involve applying a function to each frame of a trajectory. For these, you need to iterate through the frames, *even if you don't directly use the timestep*. This is because the act of iterating moves the Universe onto the next frame, changing the dynamic atom coordinates. 
+
+Trajectories can also be :ref:`sliced <slicing-trajectories>` if you only want to work on a subset of frames.
+
+.. ipython:: python
+
+    protein = u.select_atoms('protein')
+    for ts in u.trajectory[:20:4]:
+        # the radius of gyration depends on the changing positions
+        rad = protein.radius_of_gyration()
+        print('frame={}: radgyr={}'.format(ts.frame, rad))
+
+Note that after iterating over the trajectory, the frame is always set back to the first frame, even if your loop stopped before the trajectory end.
+
+.. ipython:: python
+
+    u.trajectory.frame
+
+Because MDAnalysis will pull trajectory data directly from the file it is reading from, changes to atom coordinates and box dimensions will not persist once the frame is changed. The only way to make these changes permanent is to load the trajectory into memory, or to write a new trajectory to file for every frame. For example, to set a cubic box size for every frame and write it out to a file::
+
+    with mda.Writer('with_box.trr', 'w', n_atoms=u.atoms.n_atoms) as w:
+        for ts in u.trajectory:
+            ts.dimensions = [10, 10, 10, 90, 90, 90]
+            w.write(u.atoms)
+
+    u_with_box = mda.Universe(PSF, 'with_box.trr')
+
+
+Sometimes you may wish to only transform part of the trajectory, or to not write a file out. In these cases, MDAnalysis supports :ref:`"on-the-fly" transformations <transformations>` that are performed on a frame when it is read. 
+
+

doc/source/trajectories/transformations.rst

@@ -0,0 +1,105 @@
+.. -*- coding: utf-8 -*-
+.. _transformations:
+
+On-the-fly transformations
+==========================
+
+An on-the-fly transformation is a function that silently modifies the dynamic data contained in a trajectory :class:`~MDAnalysis.coordinates.base.Timestep` (typically coordinates) as it is loaded into memory. It is called for each current time step to transform data into your desired representation. A transformation function must also return the current :class:`~MDAnalysis.coordinates.base.Timestep`, as transformations are often chained together.
+
+The :mod:`MDAnalysis.transformations` module contains a collection of transformations. For example, :func:`~MDAnalysis.transformations.fit.fit_rot_trans` can perform a mass-weighted alignment on an :class:`~MDAnalysis.core.groups.AtomGroup` to a reference.
+
+.. ipython:: python
+
+    import MDAnalysis as mda
+    from MDAnalysis.tests.datafiles import TPR, XTC
+    from MDAnalysis import transformations as trans
+
+    u = mda.Universe(TPR, XTC)
+    protein = u.select_atoms('protein')
+    align_transform = trans.fit_rot_trans(protein, protein, weights='mass')
+    u.trajectory.add_transformations(align_transform)
+
+Other implemented transformations include functions to :mod:`~MDAnalysis.transformations.translate`, :mod:`~MDAnalysis.transformations.rotate`, :mod:`~MDAnalysis.transformations.fit` an :class:`~MDAnalysis.core.groups.AtomGroup` to a reference, and :mod:`~MDAnalysis.transformations.wrap` or unwrap groups in the unit cell. 
+
+Although you can only call :meth:`~MDAnalysis.coordinates.base.ProtoReader.add_transformations` *once*, you can pass in multiple transformations in a list, which will be executed in order. For example, the below workflow:
+
+* makes all molecules whole (unwraps them over periodic boundary conditions)
+* centers the protein in the center of the box
+* wraps water back into the box
+
+.. ipython:: python
+
+    # create new Universe for new transformations
+    u = mda.Universe(TPR, XTC)
+    protein = u.select_atoms('protein')
+    water = u.select_atoms('resname SOL')
+    workflow = [trans.unwrap(u.atoms),
+                trans.center_in_box(protein, center='geometry'),
+                trans.wrap(water, compound='residues')]
+    u.trajectory.add_transformations(*workflow)
+
+If your transformation does not depend on something within the :class:`~MDAnalysis.core.universe.Universe` (e.g. a chosen :class:`~MDAnalysis.core.groups.AtomGroup`), you can also create a :class:`~MDAnalysis.core.universe.Universe` directly with transformations. The code below translates coordinates 1 angstrom up on the z-axis:
+
+.. ipython:: python
+
+    u = mda.Universe(TPR, XTC, transformations=[trans.translate([0, 0, 1])])
+
+If you need a different transformation, it is easy to implement your own.
+
+----------------------
+Custom transformations
+----------------------
+
+At its core, a transformation function must only take a :class:`~MDAnalysis.coordinates.base.Timestep` as its input and return the :class:`~MDAnalysis.coordinates.base.Timestep` as the output.
+
+.. ipython:: python
+
+    def up_by_2(ts):
+        """Translates atoms up by 2 angstrom"""
+        ts.positions += np.array([0.0, 0.0, 0.2])
+        return ts
+
+    u = mda.Universe(TPR, XTC, transformations=[up_by_2])
+
+
+If your transformation needs other arguments, you will need to wrap your core transformation with a wrapper function that can accept the other arguments.
+
+.. ipython:: python
+
+    def up_by_x(x):
+        """Translates atoms up by x angstrom"""
+        def wrapped(ts):
+            """Handles the actual Timestep"""
+            ts.positions += np.array([0.0, 0.0, float(x)])
+            return ts
+        return wrapped
+
+    # load Universe with transformations that move it up by 7 angstrom
+    u = mda.Universe(TPR, XTC, transformations=[up_by_x(5), up_by_x(2)])
+
+
+Alternatively, you can use :func:`functools.partial` to substitute the other arguments.
+
+.. ipython:: python
+
+    import functools
+
+    def up_by_x(ts, x):
+        ts.positions += np.array([0.0, 0.0, float(x)])
+        return ts
+
+    up_by_5 = functools.partial(up_by_x, x=5)
+    u = mda.Universe(TPR, XTC, transformations=[up_by_5])
+
+On-the-fly transformation functions can be applied to any property of a Timestep, not just the atom positions. For example, to give each frame of a trajectory a box:
+
+.. ipython:: python
+
+    def set_box(ts):
+        # creates box of length 10 on x-axis, 20 on y-axis, 30 on z-axis
+        # angles are all 90 degrees
+        ts.dimensions = [10, 20, 30, 90, 90, 90]
+        return ts
+
+    u = mda.Universe(TPR, XTC, transformations=[set_box])
+