Source code for pmda.rms.rmsd

# -*- Mode: python; tab-width: 4; indent-tabs-mode:nil; coding:utf-8 -*-
# vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4
#
# PMDA
# Copyright (c) 2019 The MDAnalysis Development Team and contributors
# (see the file AUTHORS for the full list of names)
#
# Released under the GNU Public Licence, v2 or any higher version

"""
Calculating Root-Mean-Square Deviations (RMSD) --- :mod:`pmda.rms`
=====================================================================

This module contains parallel versions of analysis tasks in
:mod:`MDAnalysis.analysis.rms.RMSD`.

See Also
--------
:mod:`MDAnalysis.analysis.rms`


.. autoclass:: RMSD
    :members:
    :inherited-members:

"""
from __future__ import absolute_import

from MDAnalysis.analysis import rms

import numpy as np

from pmda.parallel import ParallelAnalysisBase


[docs]class RMSD(ParallelAnalysisBase): r"""Parallel RMSD analysis. Optimally superimpose the coordinates in the :class:`~MDAnalysis.core.groups.AtomGroup` `mobile` onto `ref` for each frame in the trajectory of `mobile`, and calculate the time series of the RMSD. The single frame calculation is performed with :func:`MDAnalysis.analysis.rms.rmsd` (with ``superposition=True`` by default). Attributes ---------- rmsd : array Contains the time series of the RMSD as a `Tx3` :class:`numpy.ndarray` array with content ``[[frame, time (ps), RMSD (Å)], [...], ...]``, where `T` is the number of time steps selected in the :meth:`run` method. Parameters ---------- mobile : AtomGroup atoms that are optimally superimposed on `ref` before the RMSD is calculated for all atoms. The coordinates of `mobile` change with each frame in the trajectory. ref : AtomGroup fixed reference coordinates superposition : bool, optional ``True`` perform a RMSD-superposition, ``False`` only calculates the RMSD. The default is ``True``. See Also -------- MDAnalysis.analysis.rms.RMSD Notes ----- If you use trajectory data from simulations performed under *periodic boundary conditions* then you **must make your molecules whole** before performing RMSD calculations so that the centers of mass of the mobile and reference structure are properly superimposed. Run the analysis with :meth:`RMSD.run`, which stores the results in the array :attr:`RMSD.rmsd`. The root mean square deviation :math:`\rho(t)` of a group of :math:`N` atoms relative to a reference structure as a function of time is calculated as: .. math:: \rho(t) = \sqrt{\frac{1}{N} \sum_{i=1}^N \left(\mathbf{x}_i(t) - \mathbf{x}_i^{\text{ref}}\right)^2} The selected coordinates from `atomgroup` are optimally superimposed (translation and rotation) on the `reference` coordinates at each time step as to minimize the RMSD. At the moment, this class has far fewer features than its serial counterpart, :class:`MDAnalysis.analysis.rms.RMSD`. Examples -------- In this example we will globally fit a protein to a reference structure. The example is a DIMS trajectory of adenylate kinase, which samples a large closed-to-open transition. The trajectory is included in the MDAnalysis test data files. The data in :attr:`RMSD.rmsd` is plotted with :func:`matplotlib.pyplot.plot`:: import MDAnalysis from MDAnalysis.tests.datafiles import PSF, DCD, CRD mobile = MDAnalysis.Universe(PSF,DCD).atoms # reference closed AdK (1AKE) (with the default ref_frame=0) ref = MDAnalysis.Universe(PSF,DCD).atoms from pmda.rms import RMSD R = RMSD(mobile, ref) R.run() import matplotlib.pyplot as plt rmsd = R.rmsd.T[2] # transpose makes it easier for plotting time = rmsd[0] fig = plt.figure(figsize=(4,4)) ax = fig.add_subplot(111) ax.plot(time, rmsd, label="all") ax.legend(loc="best") ax.set_xlabel("time (ps)") ax.set_ylabel(r"RMSD ($\\AA$)") fig.savefig("rmsd_all_CORE_LID_NMP_ref1AKE.pdf") """ def __init__(self, mobile, ref, superposition=True): universe = mobile.universe super(RMSD, self).__init__(universe, (mobile, )) self._ref_pos = ref.positions.copy() self.superposition = superposition def _prepare(self): self.rmsd = None def _conclude(self): self.rmsd = np.vstack(self._results) def _single_frame(self, ts, atomgroups): return (ts.frame, ts.time, rms.rmsd(atomgroups[0].positions, self._ref_pos, superposition=self.superposition))