# -*- Mode: python; tab-width: 4; indent-tabs-mode:nil; coding:utf-8 -*-
# vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4
#
# MDAnalysis --- https://www.mdanalysis.org
# Copyright (c) 2006-2017 The MDAnalysis Development Team and contributors
# (see the file AUTHORS for the full list of names)
#
# Released under the Lesser GNU Public Licence, v2.1 or any higher version
#
# Please cite your use of MDAnalysis in published work:
#
# R. J. Gowers, M. Linke, J. Barnoud, T. J. E. Reddy, M. N. Melo, S. L. Seyler,
# D. L. Dotson, J. Domanski, S. Buchoux, I. M. Kenney, and O. Beckstein.
# MDAnalysis: A Python package for the rapid analysis of molecular dynamics
# simulations. In S. Benthall and S. Rostrup editors, Proceedings of the 15th
# Python in Science Conference, pages 102-109, Austin, TX, 2016. SciPy.
# doi: 10.25080/majora-629e541a-00e
#
# N. Michaud-Agrawal, E. J. Denning, T. B. Woolf, and O. Beckstein.
# MDAnalysis: A Toolkit for the Analysis of Molecular Dynamics Simulations.
# J. Comput. Chem. 32 (2011), 2319--2327, doi:10.1002/jcc.21787
#
"""
Calculating root mean square quantities --- :mod:`MDAnalysis.analysis.rms`
==========================================================================

:Author: Oliver Beckstein, David L. Dotson, John Detlefs
:Year: 2016
:Copyright: Lesser GNU Public License v2.1+

.. versionadded:: 0.7.7
.. versionchanged:: 0.11.0
   Added :class:`RMSF` analysis.
.. versionchanged:: 0.16.0
   Refactored RMSD to fit AnalysisBase API

The module contains code to analyze root mean square quantities such
as the coordinat root mean square distance (:class:`RMSD`) or the
per-residue root mean square fluctuations (:class:`RMSF`).

This module uses the fast QCP algorithm [Theobald2005]_ to calculate
the root mean square distance (RMSD) between two coordinate sets (as
implemented in
:func:`MDAnalysis.lib.qcprot.CalcRMSDRotationalMatrix`).

When using this module in published work please cite [Theobald2005]_.


See Also
--------
:mod:`MDAnalysis.analysis.align`
   aligning structures based on RMSD
:mod:`MDAnalysis.lib.qcprot`
   implements the fast RMSD algorithm.


Example applications
--------------------

Calculating RMSD for multiple domains
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

In this example we will globally fit a protein to a reference
structure and investigate the relative movements of domains by
computing the RMSD of the domains to the reference. The example is a
DIMS trajectory of adenylate kinase, which samples a large
closed-to-open transition. The protein consists of the CORE, LID, and
NMP domain.

* superimpose on the closed structure (frame 0 of the trajectory),
  using backbone atoms

* calculate the backbone RMSD and RMSD for CORE, LID, NMP (backbone atoms)

The trajectory is included with the test data files. The data in
:attr:`RMSD.results.rmsd` is plotted with :func:`matplotlib.pyplot.plot` (see Figure :ref:`RMSD plot figure <figure-RMSD>`)::

   import MDAnalysis
   from MDAnalysis.tests.datafiles import PSF,DCD,CRD
   u = MDAnalysis.Universe(PSF,DCD)
   ref = MDAnalysis.Universe(PSF,DCD)     # reference closed AdK (1AKE) (with the default ref_frame=0)
   #ref = MDAnalysis.Universe(PSF,CRD)    # reference open AdK (4AKE)

   import MDAnalysis.analysis.rms

   R = MDAnalysis.analysis.rms.RMSD(u, ref,
              select="backbone",             # superimpose on whole backbone of the whole protein
              groupselections=["backbone and (resid 1-29 or resid 60-121 or resid 160-214)",   # CORE
                               "backbone and resid 122-159",                                   # LID
                               "backbone and resid 30-59"])                                    # NMP
   R.run()

   import matplotlib.pyplot as plt
   rmsd = R.results.rmsd.T   # transpose makes it easier for plotting
   time = rmsd[1]
   fig = plt.figure(figsize=(4,4))
   ax = fig.add_subplot(111)
   ax.plot(time, rmsd[2], 'k-',  label="all")
   ax.plot(time, rmsd[3], 'k--', label="CORE")
   ax.plot(time, rmsd[4], 'r--', label="LID")
   ax.plot(time, rmsd[5], 'b--', label="NMP")
   ax.legend(loc="best")
   ax.set_xlabel("time (ps)")
   ax.set_ylabel(r"RMSD ($\\AA$)")
   fig.savefig("rmsd_all_CORE_LID_NMP_ref1AKE.pdf")

.. _figure-RMSD:

.. figure:: /images/RSMD_plot.png
      :scale: 50 %
      :alt: RMSD plot

      RMSD plot for backbone and CORE, LID, NMP domain of the protein.


Functions
---------

.. autofunction:: rmsd

Analysis classes
----------------

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

   .. attribute:: results.rmsd

       Contains the time series of the RMSD as an N×3 :class:`numpy.ndarray`
       array with content ``[[frame, time (ps), RMSD (A)], [...], ...]``.

       .. versionadded:: 2.0.0

   .. attribute:: rmsd

       Alias to the :attr:`results.rmsd` attribute.

       .. deprecated:: 2.0.0
          Will be removed in MDAnalysis 3.0.0. Please use :attr:`results.rmsd`
          instead.


.. autoclass:: RMSF
   :members:
   :inherited-members:

   .. attribute:: results.rmsf

      Results are stored in this N-length :class:`numpy.ndarray` array,
      giving RMSFs for each of the given atoms.

      .. versionadded:: 2.0.0

   .. attribute:: rmsf

      Alias to the :attr:`results.rmsf` attribute.

      .. deprecated:: 2.0.0
         Will be removed in MDAnalysis 3.0.0. Please use :attr:`results.rmsf`
         instead.

"""
import numpy as np

import logging
import warnings

from ..lib import qcprot as qcp
from ..analysis.base import AnalysisBase, ResultsGroup
from ..exceptions import SelectionError
from ..lib.util import asiterable, iterable, get_weights


logger = logging.getLogger("MDAnalysis.analysis.rmsd")


def rmsd(a, b, weights=None, center=False, superposition=False):
    r"""Returns RMSD between two coordinate sets `a` and `b`.

    `a` and `b` are arrays of the coordinates of N atoms of shape
    :math:`N times 3` as generated by, e.g.,
    :meth:`MDAnalysis.core.groups.AtomGroup.positions`.

    Note
    ----
    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.


    Parameters
    ----------
    a : array_like
        coordinates to align to `b`
    b : array_like
        coordinates to align to (same shape as `a`)
    weights : array_like (optional)
        1D array with weights, use to compute weighted average
    center : bool (optional)
        subtract center of geometry before calculation. With weights given
        compute weighted average as center.
    superposition : bool (optional)
        perform a rotational and translational superposition with the fast QCP
        algorithm [Theobald2005]_ before calculating the RMSD; implies
        ``center=True``.

    Returns
    -------
    rmsd : float
        RMSD between `a` and `b`

    Notes
    -----
    The RMSD :math:`\rho(t)` as a function of time is calculated as

    .. math::

       \rho(t) = \sqrt{\frac{1}{N} \sum_{i=1}^N w_i \left(\mathbf{x}_i(t)
                                - \mathbf{x}_i^{\text{ref}}\right)^2}

    It is the Euclidean distance in configuration space of the current
    configuration (possibly after optimal translation and rotation) from a
    reference configuration divided by :math:`1/\sqrt{N}` where :math:`N` is
    the number of coordinates.

    The weights :math:`w_i` are calculated from the input weights
    `weights` :math:`w'_i` as relative to the mean:

    .. math::

       w_i = \frac{w'_i}{\langle w' \rangle}


    Example
    -------
    >>> import MDAnalysis as mda
    >>> from MDAnalysis.analysis.rms import rmsd
    >>> from MDAnalysis.tests.datafiles import PSF, DCD
    >>> u = mda.Universe(PSF, DCD)
    >>> bb = u.select_atoms('backbone')
    >>> A = bb.positions.copy()  # coordinates of first frame
    >>> _ = u.trajectory[-1]  # forward to last frame
    >>> B = bb.positions.copy()  # coordinates of last frame
    >>> rmsd(A, B, center=True)
    6.838544558398293


    .. versionchanged:: 0.8.1
       *center* keyword added
    .. versionchanged:: 0.14.0
       *superposition* keyword added

    """

    a = np.asarray(a, dtype=np.float64)
    b = np.asarray(b, dtype=np.float64)
    N = b.shape[0]
    if a.shape != b.shape:
        raise ValueError("a and b must have same shape")

    # superposition only works if structures are centered
    if center or superposition:
        # make copies (do not change the user data!)
        # weights=None is equivalent to all weights 1
        a = a - np.average(a, axis=0, weights=weights)
        b = b - np.average(b, axis=0, weights=weights)

    if weights is not None:
        if len(weights) != len(a):
            raise ValueError("weights must have same length as a and b")
        # weights are constructed as relative to the mean
        weights = np.asarray(weights, dtype=np.float64) / np.mean(weights)

    if superposition:
        return qcp.CalcRMSDRotationalMatrix(a, b, N, None, weights)
    else:
        if weights is not None:
            return np.sqrt(np.sum(weights[:, np.newaxis] * ((a - b) ** 2)) / N)
        else:
            return np.sqrt(np.sum((a - b) ** 2) / N)


def process_selection(select):
    """Return a canonical selection dictionary.

    Parameters
    ----------
    select : str or tuple or dict

        - `str` -> Any valid string selection
        - `dict` -> ``{'mobile':sel1, 'reference':sel2}``
        - `tuple` -> ``(sel1, sel2)``

    Returns
    -------
    dict
        selections for 'reference' and 'mobile'. Values are guarenteed to be
        iterable (so that one can provide selections to retain order)

    Notes
    -----
    The dictionary input for `select` can be generated by
    :func:`fasta2select` based on a ClustalW_ or STAMP_ sequence alignment.
    """

    if isinstance(select, str):
        select = {"reference": str(select), "mobile": str(select)}
    elif type(select) is tuple:
        try:
            select = {"mobile": select[0], "reference": select[1]}
        except IndexError:
            raise IndexError(
                "select must contain two selection strings "
                "(reference, mobile)"
            ) from None
    elif type(select) is dict:
        # compatability hack to use new nomenclature
        try:
            select["mobile"]
            select["reference"]
        except KeyError:
            raise KeyError(
                "select dictionary must contain entries for keys "
                "'mobile' and 'reference'."
            ) from None
    else:
        raise TypeError("'select' must be either a string, 2-tuple, or dict")
    select["mobile"] = asiterable(select["mobile"])
    select["reference"] = asiterable(select["reference"])
    return select


class RMSD(AnalysisBase):
    r"""Class to perform RMSD analysis on a trajectory.

    The RMSD will be computed for two groups of atoms and all frames in the
    trajectory belonging to `atomgroup`. The groups of atoms are obtained by
    applying the selection selection `select` to the changing `atomgroup` and
    the fixed `reference`.

    Note
    ----
    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 selected
    and reference structure are properly superimposed.


    Run the analysis with :meth:`RMSD.run`, which stores the results
    in the array :attr:`RMSD.results.rmsd`.


    .. versionchanged:: 1.0.0
       ``save()`` method was removed, use ``np.savetxt()`` on
       :attr:`RMSD.results.rmsd` instead.
    .. versionchanged:: 2.0.0
       :attr:`rmsd` results are now stored in a
       :class:`MDAnalysis.analysis.base.Results` instance.
    .. versionchanged:: 2.8.0
       introduced :meth:`get_supported_backends` allowing for parallel
       execution on ``multiprocessing`` and ``dask`` backends.
    """

    _analysis_algorithm_is_parallelizable = True

    @classmethod
    def get_supported_backends(cls):
        return (
            "serial",
            "multiprocessing",
            "dask",
        )

    def __init__(
        self,
        atomgroup,
        reference=None,
        select="all",
        groupselections=None,
        weights=None,
        weights_groupselections=False,
        tol_mass=0.1,
        ref_frame=0,
        **kwargs,
    ):
        r"""Parameters
        ----------
        atomgroup : AtomGroup or Universe
            Group of atoms for which the RMSD is calculated. If a trajectory is
            associated with the atoms then the computation iterates over the
            trajectory.
        reference : AtomGroup or Universe (optional)
            Group of reference atoms; if ``None`` then the current frame of
            `atomgroup` is used.
        select : str or dict or tuple (optional)
            The selection to operate on; can be one of:

            1. any valid selection string for
               :meth:`~MDAnalysis.core.groups.AtomGroup.select_atoms` that
               produces identical selections in `atomgroup` and `reference`; or

            2. a dictionary ``{'mobile': sel1, 'reference': sel2}`` where *sel1*
               and *sel2* are valid selection strings that are applied to
               `atomgroup` and `reference` respectively (the
               :func:`MDAnalysis.analysis.align.fasta2select` function returns such
               a dictionary based on a ClustalW_ or STAMP_ sequence alignment); or

            3. a tuple ``(sel1, sel2)``

            When using 2. or 3. with *sel1* and *sel2* then these selection strings
            are applied to `atomgroup` and `reference` respectively and should
            generate *groups of equivalent atoms*.  *sel1* and *sel2* can each also
            be a *list of selection strings* to generate a
            :class:`~MDAnalysis.core.groups.AtomGroup` with defined atom order as
            described under :ref:`ordered-selections-label`).

        groupselections : list (optional)
            A list of selections as described for `select`, with the difference
            that these selections are *always applied to the full universes*,
            i.e., ``atomgroup.universe.select_atoms(sel1)`` and
            ``reference.universe.select_atoms(sel2)``. Each selection describes
            additional RMSDs to be computed *after the structures have been
            superimposed* according to `select`. No additional fitting is
            performed.The output contains one additional column for each
            selection.

            .. Note:: Experimental feature. Only limited error checking
                      implemented.

        weights : {"mass", ``None``} or array_like (optional)
             1. "mass" will use masses as weights for both `select` and `groupselections`.

             2. ``None`` will weigh each atom equally for both `select` and `groupselections`.

             3. If 1D float array of the same length as `atomgroup` is provided,
             use each element of the `array_like` as a weight for the
             corresponding atom in `select`, and assumes ``None`` for `groupselections`.

        weights_groupselections : False or list of {"mass", ``None`` or array_like} (optional)
             1. ``False`` will apply imposed weights to `groupselections` from
             ``weights`` option if ``weights`` is either ``"mass"`` or ``None``.
             Otherwise will assume a list of length equal to length of
             `groupselections` filled with ``None`` values.

             2. A list of {"mass", ``None`` or array_like} with the length of `groupselections`
             will apply the weights to `groupselections` correspondingly.

        tol_mass : float (optional)
             Reject match if the atomic masses for matched atoms differ by more
             than `tol_mass`.
        ref_frame : int (optional)
             frame index to select frame from `reference`
        verbose : bool (optional)
             Show detailed progress of the calculation if set to ``True``; the
             default is ``False``.

        Raises
        ------
        SelectionError
             If the selections from `atomgroup` and `reference` do not match.
        TypeError
             If `weights` or `weights_groupselections` is not of the appropriate type;
             see also :func:`MDAnalysis.lib.util.get_weights`
        ValueError
             If `weights` are not compatible with `atomgroup` (not the same
             length) or if it is not a 1D array (see
             :func:`MDAnalysis.lib.util.get_weights`).

             A :exc:`ValueError` is also raised if the length of `weights_groupselections`
             are not compatible with `groupselections`.

        Notes
        -----
        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 w_i \left(\mathbf{x}_i(t)
                                    - \mathbf{x}_i^{\text{ref}}\right)^2}

        The weights :math:`w_i` are calculated from the input weights `weights`
        :math:`w'_i` as relative to the mean of the input weights:

        .. math::

           w_i = \frac{w'_i}{\langle w' \rangle}

        The selected coordinates from `atomgroup` are optimally superimposed
        (translation and rotation) on the `reference` coordinates at each time step
        as to minimize the RMSD. Douglas Theobald's fast QCP algorithm
        [Theobald2005]_ is used for the rotational superposition and to calculate
        the RMSD (see :mod:`MDAnalysis.lib.qcprot` for implementation details).

        The class runs various checks on the input to ensure that the two atom
        groups can be compared. This includes a comparison of atom masses (i.e.,
        only the positions of atoms of the same mass will be considered to be
        correct for comparison). If masses should not be checked, just set
        `tol_mass` to a large value such as 1000.

        .. _ClustalW: http://www.clustal.org/
        .. _STAMP: http://www.compbio.dundee.ac.uk/manuals/stamp.4.2/


        See Also
        --------
        rmsd


        .. versionadded:: 0.7.7
        .. versionchanged:: 0.8
           `groupselections` added
        .. versionchanged:: 0.16.0
           Flexible weighting scheme with new `weights` keyword.
        .. deprecated:: 0.16.0
           Instead of ``mass_weighted=True`` (removal in 0.17.0) use new
           ``weights='mass'``; refactored to fit with AnalysisBase API
        .. versionchanged:: 0.17.0
           removed deprecated `mass_weighted` keyword; `groupselections`
           are *not* rotationally superimposed any more.
        .. versionchanged:: 1.0.0
           `filename` keyword was removed.

        """
        super(RMSD, self).__init__(atomgroup.universe.trajectory, **kwargs)
        self.atomgroup = atomgroup
        self.reference = reference if reference is not None else self.atomgroup

        select = process_selection(select)
        self.groupselections = (
            [process_selection(s) for s in groupselections]
            if groupselections is not None
            else []
        )
        self.weights = weights
        self.tol_mass = tol_mass
        self.ref_frame = ref_frame
        self.weights_groupselections = weights_groupselections
        self.ref_atoms = self.reference.select_atoms(*select["reference"])
        self.mobile_atoms = self.atomgroup.select_atoms(*select["mobile"])

        if len(self.ref_atoms) != len(self.mobile_atoms):
            err = (
                "Reference and trajectory atom selections do "
                "not contain the same number of atoms: "
                "N_ref={0:d}, N_traj={1:d}".format(
                    self.ref_atoms.n_atoms, self.mobile_atoms.n_atoms
                )
            )
            logger.exception(err)
            raise SelectionError(err)
        logger.info(
            "RMS calculation " "for {0:d} atoms.".format(len(self.ref_atoms))
        )
        mass_mismatches = (
            np.absolute((self.ref_atoms.masses - self.mobile_atoms.masses))
            > self.tol_mass
        )

        if np.any(mass_mismatches):
            # diagnostic output:
            logger.error("Atoms: reference | mobile")
            for ar, at in zip(self.ref_atoms, self.mobile_atoms):
                if ar.name != at.name:
                    logger.error(
                        "{0!s:>4} {1:3d} {2!s:>3} {3!s:>3} {4:6.3f}"
                        "|  {5!s:>4} {6:3d} {7!s:>3} {8!s:>3}"
                        "{9:6.3f}".format(
                            ar.segid,
                            ar.resid,
                            ar.resname,
                            ar.name,
                            ar.mass,
                            at.segid,
                            at.resid,
                            at.resname,
                            at.name,
                            at.mass,
                        )
                    )
            errmsg = (
                "Inconsistent selections, masses differ by more than"
                "{0:f}; mis-matching atoms"
                "are shown above.".format(self.tol_mass)
            )
            logger.error(errmsg)
            raise SelectionError(errmsg)
        del mass_mismatches

        # TODO:
        # - make a group comparison a class that contains the checks above
        # - use this class for the *select* group and the additional
        #   *groupselections* groups each a dict with reference/mobile
        self._groupselections_atoms = [
            {
                "reference": self.reference.universe.select_atoms(
                    *s["reference"]
                ),
                "mobile": self.atomgroup.universe.select_atoms(*s["mobile"]),
            }
            for s in self.groupselections
        ]
        # sanity check
        for igroup, (sel, atoms) in enumerate(
            zip(self.groupselections, self._groupselections_atoms)
        ):
            if len(atoms["mobile"]) != len(atoms["reference"]):
                logger.exception("SelectionError: Group Selection")
                raise SelectionError(
                    "Group selection {0}: {1} | {2}: Reference and trajectory "
                    "atom selections do not contain the same number of atoms: "
                    "N_ref={3}, N_traj={4}".format(
                        igroup,
                        sel["reference"],
                        sel["mobile"],
                        len(atoms["reference"]),
                        len(atoms["mobile"]),
                    )
                )

        # check weights type
        acceptable_dtypes = (np.dtype("float64"), np.dtype("int64"))
        msg = (
            "weights should only be 'mass', None or 1D float array."
            "For weights on groupselections, "
            "use **weight_groupselections**"
        )

        if iterable(self.weights):
            element_lens = []
            for element in self.weights:
                if iterable(element):
                    element_lens.append(len(element))
                else:
                    element_lens.append(1)
            if np.unique(element_lens).size > 1:
                # jagged data structure
                raise TypeError(msg)
            if np.array(element).dtype not in acceptable_dtypes:
                raise TypeError(msg)

        if iterable(self.weights) or self.weights != "mass":
            get_weights(self.mobile_atoms, self.weights)

        if self.weights_groupselections:
            if len(self.weights_groupselections) != len(self.groupselections):
                raise ValueError(
                    "Length of weights_groupselections is not equal to "
                    "length of groupselections "
                )
            for weights, atoms, selection in zip(
                self.weights_groupselections,
                self._groupselections_atoms,
                self.groupselections,
            ):
                try:
                    if iterable(weights) or weights != "mass":
                        get_weights(atoms["mobile"], weights)
                except Exception as e:
                    raise type(e)(
                        str(e)
                        + " happens in selection %s" % selection["mobile"]
                    )

    def _prepare(self):
        self._n_atoms = self.mobile_atoms.n_atoms
        if not self.weights_groupselections:
            if not iterable(
                self.weights
            ):  # apply 'mass' or 'None' to groupselections
                self.weights_groupselections = [self.weights] * len(
                    self.groupselections
                )
            else:
                self.weights_groupselections = [None] * len(
                    self.groupselections
                )

        for igroup, (weights, atoms) in enumerate(
            zip(self.weights_groupselections, self._groupselections_atoms)
        ):
            if str(weights) == "mass":
                self.weights_groupselections[igroup] = atoms["mobile"].masses
            if weights is not None:
                self.weights_groupselections[igroup] = np.asarray(
                    self.weights_groupselections[igroup], dtype=np.float64
                ) / np.mean(self.weights_groupselections[igroup])
        # add the array of weights to weights_select
        self.weights_select = get_weights(self.mobile_atoms, self.weights)
        self.weights_ref = get_weights(self.ref_atoms, self.weights)
        if self.weights_select is not None:
            self.weights_select = np.asarray(
                self.weights_select, dtype=np.float64
            ) / np.mean(self.weights_select)
            self.weights_ref = np.asarray(
                self.weights_ref, dtype=np.float64
            ) / np.mean(self.weights_ref)

        current_frame = self.reference.universe.trajectory.ts.frame

        try:
            # Move to the ref_frame
            # (coordinates MUST be stored in case the ref traj is advanced
            # elsewhere or if ref == mobile universe)
            self.reference.universe.trajectory[self.ref_frame]
            self._ref_com = self.ref_atoms.center(self.weights_ref)
            # makes a copy
            self._ref_coordinates = self.ref_atoms.positions - self._ref_com
            if self._groupselections_atoms:
                self._groupselections_ref_coords64 = [
                    (
                        self.reference.select_atoms(
                            *s["reference"]
                        ).positions.astype(np.float64)
                    )
                    for s in self.groupselections
                ]
        finally:
            # Move back to the original frame
            self.reference.universe.trajectory[current_frame]

        self._ref_coordinates64 = self._ref_coordinates.astype(np.float64)

        if self._groupselections_atoms:
            # Only carry out a rotation if we want to calculate secondary
            # RMSDs.
            # R: rotation matrix that aligns r-r_com, x~-x~com
            #    (x~: selected coordinates, x: all coordinates)
            # Final transformed traj coordinates: x' = (x-x~_com)*R + ref_com
            self._rot = np.zeros(9, dtype=np.float64)  # allocate space
            self._R = self._rot.reshape(3, 3)
        else:
            self._rot = None

        self.results.rmsd = np.zeros(
            (self.n_frames, 3 + len(self._groupselections_atoms))
        )

        self._mobile_coordinates64 = self.mobile_atoms.positions.copy().astype(
            np.float64
        )

    def _get_aggregator(self):
        return ResultsGroup(lookup={"rmsd": ResultsGroup.ndarray_vstack})

    def _single_frame(self):
        mobile_com = self.mobile_atoms.center(self.weights_select).astype(
            np.float64
        )
        self._mobile_coordinates64[:] = self.mobile_atoms.positions
        self._mobile_coordinates64 -= mobile_com

        self.results.rmsd[self._frame_index, :2] = (
            self._ts.frame,
            self._trajectory.time,
        )

        if self._groupselections_atoms:
            # superimpose structures: MDAnalysis qcprot needs Nx3 coordinate
            # array with float64 datatype (float32 leads to errors up to 1e-3 in
            # RMSD). Note that R is defined in such a way that it acts **to the
            # left** so that we can easily use broadcasting and save one
            # expensive numpy transposition.

            self.results.rmsd[self._frame_index, 2] = (
                qcp.CalcRMSDRotationalMatrix(
                    self._ref_coordinates64,
                    self._mobile_coordinates64,
                    self._n_atoms,
                    self._rot,
                    self.weights_select,
                )
            )

            self._R[:, :] = self._rot.reshape(3, 3)
            # Transform each atom in the trajectory (use inplace ops to
            # avoid copying arrays) (Marginally (~3%) faster than
            # "ts.positions[:] = (ts.positions - x_com) * R + ref_com".)
            self._ts.positions[:] -= mobile_com

            # R acts to the left & is broadcasted N times.
            self._ts.positions[:] = np.dot(self._ts.positions, self._R)
            self._ts.positions += self._ref_com

            # 2) calculate secondary RMSDs (without any further
            #    superposition)
            for igroup, (refpos, atoms) in enumerate(
                zip(
                    self._groupselections_ref_coords64,
                    self._groupselections_atoms,
                ),
                3,
            ):
                self.results.rmsd[self._frame_index, igroup] = rmsd(
                    refpos,
                    atoms["mobile"].positions,
                    weights=self.weights_groupselections[igroup - 3],
                    center=False,
                    superposition=False,
                )
        else:
            # only calculate RMSD by setting the Rmatrix to None (no need
            # to carry out the rotation as we already get the optimum RMSD)
            self.results.rmsd[self._frame_index, 2] = (
                qcp.CalcRMSDRotationalMatrix(
                    self._ref_coordinates64,
                    self._mobile_coordinates64,
                    self._n_atoms,
                    None,
                    self.weights_select,
                )
            )

    @property
    def rmsd(self):
        wmsg = (
            "The `rmsd` attribute was deprecated in MDAnalysis 2.0.0 and "
            "will be removed in MDAnalysis 3.0.0. Please use "
            "`results.rmsd` instead."
        )
        warnings.warn(wmsg, DeprecationWarning)
        return self.results.rmsd


class RMSF(AnalysisBase):
    r"""Calculate RMSF of given atoms across a trajectory.

    Note
    ----
    No RMSD-superposition is performed; it is assumed that the user is
    providing a trajectory where the protein of interest has been structurally
    aligned to a reference structure (see the Examples section below). The
    protein also has be whole because periodic boundaries are not taken into
    account.


    Run the analysis with :meth:`RMSF.run`, which stores the results
    in the array :attr:`RMSF.results.rmsf`.

    """

    @classmethod
    def get_supported_backends(cls):
        return ("serial",)

    def __init__(self, atomgroup, **kwargs):
        r"""Parameters
        ----------
        atomgroup : AtomGroup
            Atoms for which RMSF is calculated
        verbose : bool (optional)
             Show detailed progress of the calculation if set to ``True``; the
             default is ``False``.

        Raises
        ------
        ValueError
             raised if negative values are calculated, which indicates that a
             numerical overflow or underflow occured


        Notes
        -----
        The root mean square fluctuation of an atom :math:`i` is computed as the
        time average

        .. math::

          \rho_i = \sqrt{\left\langle (\mathbf{x}_i - \langle\mathbf{x}_i\rangle)^2 \right\rangle}

        No mass weighting is performed.

        This method implements an algorithm for computing sums of squares while
        avoiding overflows and underflows :footcite:p:`Welford1962`.


        Examples
        --------

        In this example we calculate the residue RMSF fluctuations by analyzing
        the :math:`\text{C}_\alpha` atoms. First we need to fit the trajectory
        to the average structure as a reference. That requires calculating the
        average structure first. Because we need to analyze and manipulate the
        same trajectory multiple times, we are going to load it into memory
        using the :mod:`~MDAnalysis.coordinates.MemoryReader`. (If your
        trajectory does not fit into memory, you will need to :ref:`write out
        intermediate trajectories <writing-trajectories>` to disk or
        :ref:`generate an in-memory universe
        <creating-in-memory-trajectory-label>` that only contains, say, the
        protein)::

           import MDAnalysis as mda
           from MDAnalysis.analysis import align

           from MDAnalysis.tests.datafiles import TPR, XTC

           u = mda.Universe(TPR, XTC, in_memory=True)
           protein = u.select_atoms("protein")

           # 1) the current trajectory contains a protein split across
           #    periodic boundaries, so we first make the protein whole and
           #    center it in the box using on-the-fly transformations
           import MDAnalysis.transformations as trans

           not_protein = u.select_atoms('not protein')
           transforms = [trans.unwrap(protein),
                         trans.center_in_box(protein, wrap=True),
                         trans.wrap(not_protein)]
           u.trajectory.add_transformations(*transforms)

           # 2) fit to the initial frame to get a better average structure
           #    (the trajectory is changed in memory)
           prealigner = align.AlignTraj(u, u, select="protein and name CA",
                                        in_memory=True).run()

           # 3) reference = average structure
           ref_coordinates = u.trajectory.timeseries(asel=protein).mean(axis=1)
           # make a reference structure (need to reshape into a 1-frame
           # "trajectory")
           reference = mda.Merge(protein).load_new(ref_coordinates[:, None, :],
                                                   order="afc")

        We created a new universe ``reference`` that contains a single frame
        with the averaged coordinates of the protein.  Now we need to fit the
        whole trajectory to the reference by minimizing the RMSD. We use
        :class:`MDAnalysis.analysis.align.AlignTraj`::

           aligner = align.AlignTraj(u, reference,
                                     select="protein and name CA",
                                     in_memory=True).run()

        The trajectory is now fitted to the reference (the RMSD is stored as
        `aligner.results.rmsd` for further inspection). Now we can calculate
        the RMSF::

           from MDAnalysis.analysis.rms import RMSF

           calphas = protein.select_atoms("name CA")
           rmsfer = RMSF(calphas, verbose=True).run()

        and plot::

           import matplotlib.pyplot as plt

           plt.plot(calphas.resnums, rmsfer.results.rmsf)



        References
        ----------
        .. footbibliography::


        .. versionadded:: 0.11.0
        .. versionchanged:: 0.16.0
           refactored to fit with AnalysisBase API
        .. deprecated:: 0.16.0
           the keyword argument `quiet` is deprecated in favor of `verbose`.
        .. versionchanged:: 0.17.0
           removed unused keyword `weights`
        .. versionchanged:: 1.0.0
           Support for the ``start``, ``stop``, and ``step`` keywords has been
           removed. These should instead be passed to :meth:`RMSF.run`.

        """
        super(RMSF, self).__init__(atomgroup.universe.trajectory, **kwargs)
        self.atomgroup = atomgroup

    def _prepare(self):
        self.sumsquares = np.zeros((self.atomgroup.n_atoms, 3))
        self.mean = self.sumsquares.copy()

    def _single_frame(self):
        k = self._frame_index
        self.sumsquares += (k / (k + 1.0)) * (
            self.atomgroup.positions - self.mean
        ) ** 2
        self.mean = (k * self.mean + self.atomgroup.positions) / (k + 1)

    def _conclude(self):
        k = self._frame_index
        self.results.rmsf = np.sqrt(self.sumsquares.sum(axis=1) / (k + 1))

        if not (self.results.rmsf >= 0).all():
            raise ValueError(
                "Some RMSF values negative; overflow "
                + "or underflow occurred"
            )

    @property
    def rmsf(self):
        wmsg = (
            "The `rmsf` attribute was deprecated in MDAnalysis 2.0.0 and "
            "will be removed in MDAnalysis 3.0.0. Please use "
            "`results.rmsd` instead."
        )
        warnings.warn(wmsg, DeprecationWarning)
        return self.results.rmsf