# Source code for MDAnalysis.analysis.rdf

    # -*- 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
    #

    r"""Radial Distribution Functions --- :mod:`MDAnalysis.analysis.rdf`
    ====================================================================

    This module contains two classes to calculate radial
    `pair distribution functions`_ (`radial distribution functions`_ or "RDF").
    The RDF :math:`g_{ab}(r)` between types of particles :math:`a` and :math:`b` is

    .. _equation-gab:

    .. math::

       g_{ab}(r) = \frac{1}{N_{a}} \frac{1}{N_{b}/V} \sum_{i=1}^{N_a} \sum_{j=1}^{N_b}
                   \langle \delta(|\mathbf{r}_i - \mathbf{r}_j| - r) \rangle

    which is normalized so that the RDF becomes 1 for large separations in a
    homogenous system. The RDF effectively counts the average number of :math:`b`
    neighbours in a shell at distance :math:`r` around a :math:`a` particle and
    represents it as a density.

    The radial cumulative distribution function is

    .. math::

       G_{ab}(r) = \int_0^r \!\!dr' 4\pi r'^2 g_{ab}(r')

    and the average number of :math:`b` particles within radius :math:`r`

    .. _equation-countab:

    .. math::

       N_{ab}(r) = \rho G_{ab}(r)

    (with the appropriate density :math:`\rho`). The latter function can be used to
    compute, for instance, coordination numbers such as the number of neighbors in
    the first solvation shell :math:`N(r_1)` where :math:`r_1` is the position of
    the first minimum in :math:`g(r)`.

    We provide options for calculating the density of particle :math:`b`
    in a shell at distance :math:`r` around a :math:`a` particle, which is

    .. _equation-nab:

    .. math::
       n_{ab}(r) = \rho g_{ab}(r)

    .. _`pair distribution functions`:
       https://en.wikipedia.org/wiki/Pair_distribution_function
    .. _`radial distribution functions`:
       https://en.wikipedia.org/wiki/Radial_distribution_function

    .. Not Implemented yet:
    .. - Structure factor?
    .. - Coordination number
    """
    import warnings
    import numpy as np

    from ..lib import distances
    from .base import AnalysisBase, ResultsGroup

    def nested_array_sum(arrs):
        r"""Custom aggregator for nested arrays

        This function takes a nested list or tuple of NumPy arrays, flattens it
        into a single list, and aggregates the elements at alternating indices
        into two separate arrays. The first array accumulates elements at even
        indices, while the second accumulates elements at odd indices.

        Parameters
        ----------
        arrs : list
            List of arrays or nested lists of arrays

        Returns
        -------
        list of ndarray
            A list containing two NumPy arrays:
            - The first array is the sum of all elements at even indices
              in the sum of flattened arrays.
            - The second array is the sum of all elements at odd indices
              in the sum of flattened arrays.
        """

        def flatten(arr):
            if isinstance(arr, (list, tuple)):
                return [item for sublist in arr for item in flatten(sublist)]
            return [arr]

        flat = flatten(arrs)
        aggregated_arr = [np.zeros_like(flat[0]), np.zeros_like(flat[1])]
        for i in range(len(flat) // 2):
            aggregated_arr[0] += flat[2 * i]  # 0, 2, 4, ...
            aggregated_arr[1] += flat[2 * i + 1]  # 1, 3, 5, ...
        return aggregated_arr

    class InterRDF(AnalysisBase):
        r"""Radial distribution function

        :class:`InterRDF` is a tool to calculate average radial distribution
        functions between two groups of atoms. Suppose we have two AtomGroups ``A``
        and ``B``. ``A`` contains atom ``A1``, ``A2``, and ``B`` contains ``B1``,
        ``B2``. Given ``A`` and ``B`` to :class:`InterRDF`, the output will be the
        average of RDFs between ``A1`` and ``B1``, ``A1`` and ``B2``, ``A2`` and
        ``B1``, ``A2`` and ``B2``. A typical application is to calculate the RDF of
        solvent with itself or with another solute.

        The :ref:`radial distribution function<equation-gab>` is calculated by
        histogramming distances between all particles in `g1` and `g2` while taking
        periodic boundary conditions into account via the minimum image
        convention.

        The `exclusion_block` keyword may be used to exclude a set of distances
        from the calculations.

        Results are available in the attributes :attr:`results.rdf`
        and :attr:`results.count`.

        Parameters
        ----------
        g1 : AtomGroup
            First AtomGroup
        g2 : AtomGroup
            Second AtomGroup
        nbins : int
            Number of bins in the histogram
        range : tuple or list
            The size of the RDF
        norm : str, {'rdf', 'density', 'none'}
              For 'rdf' calculate :math:`g_{ab}(r)`. For
              'density' the :ref:`single particle density<equation-nab>`
              :math:`n_{ab}(r)` is computed. 'none' computes the number of
              particles occurences in each spherical shell.

              .. versionadded:: 2.3.0

        exclusion_block : tuple
            A tuple representing the tile to exclude from the distance array.
        exclude_same : str
            Will exclude pairs of atoms that share the same "residue", "segment", or "chain".
            Those are the only valid values. This is intended to remove atoms that are
            spatially correlated due to direct bonded connections.
        verbose : bool
            Show detailed progress of the calculation if set to `True`
        backend : {'serial', 'OpenMP', 'distopia'}, optional
            Keyword selecting the type of acceleration of the distance calculations.
            Default is 'serial'.

            .. versionadded:: 2.10.0

        Attributes
        ----------
        results.bins : numpy.ndarray
           :class:`numpy.ndarray` of the centers of the `nbins` histogram
           bins.

           .. versionadded:: 2.0.0

        bins : numpy.ndarray
           Alias to the :attr:`results.bins` attribute.

           .. deprecated:: 2.0.0
               This attribute will be removed in 3.0.0.
               Use :attr:`results.bins` instead.

        results.edges : numpy.ndarray

          :class:`numpy.ndarray` of the `nbins + 1` edges of the histogram
          bins.

           .. versionadded:: 2.0.0

        edges : numpy.ndarray

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

           .. deprecated:: 2.0.0
               This attribute will be removed in 3.0.0.
               Use :attr:`results.edges` instead.

        results.rdf : numpy.ndarray
          :class:`numpy.ndarray` of the :ref:`radial distribution
          function<equation-gab>` values for the :attr:`results.bins`.

           .. versionadded:: 2.0.0

        rdf : numpy.ndarray
           Alias to the :attr:`results.rdf` attribute.

           .. deprecated:: 2.0.0
               This attribute will be removed in 3.0.0.
               Use :attr:`results.rdf` instead.

        results.count : numpy.ndarray
          :class:`numpy.ndarray` representing the radial histogram, i.e.,
          the raw counts, for all :attr:`results.bins`.

           .. versionadded:: 2.0.0

        count : numpy.ndarray
           Alias to the :attr:`results.count` attribute.

           .. deprecated:: 2.0.0
               This attribute will be removed in 3.0.0.
               Use :attr:`results.count` instead.

        Example
        -------
        First create the :class:`InterRDF` object, by supplying two
        AtomGroups then use the :meth:`run` method ::

          rdf = InterRDF(ag1, ag2)
          rdf.run()

        Results are available through the :attr:`results.bins` and
        :attr:`results.rdf` attributes::

          plt.plot(rdf.results.bins, rdf.results.rdf)

        The `exclusion_block` keyword allows the masking of pairs from
        within the same molecule. For example, if there are 7 of each
        atom in each molecule, the exclusion mask ``(7, 7)`` can be used.

        .. versionadded:: 0.13.0

        .. versionchanged:: 1.0.0
           Support for the `start`, `stop`, and `step` keywords has been
           removed. These should instead be passed to :meth:`InterRDF.run`.

        .. versionchanged:: 2.0.0
           Store results as attributes `bins`, `edges`, `rdf` and `count`
           of the `results` attribute of
           :class:`~MDAnalysis.analysis.AnalysisBase`.

        .. versionchanged:: 2.10.0
           Enabled **parallel execution** with the ``multiprocessing`` and ``dask``
           backends; use the new method :meth:`get_supported_backends` to see all
           supported backends.
        """

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

        _analysis_algorithm_is_parallelizable = True

        def __init__(
            self,
            g1,
            g2,
            nbins=75,
            range=(0.0, 15.0),
            norm="rdf",
            exclusion_block=None,
            exclude_same=None,
            backend="serial",
            **kwargs,
        ):
            super(InterRDF, self).__init__(g1.universe.trajectory, **kwargs)
            self.g1 = g1
            self.g2 = g2
            self.norm = str(norm).lower()

            self.rdf_settings = {"bins": nbins, "range": range}
            self._exclusion_block = exclusion_block
            if exclude_same is not None and exclude_same not in [
                "residue",
                "segment",
                "chain",
            ]:
                raise ValueError(
                    "The exclude_same argument to InterRDF must be None, 'residue', 'segment' "
                    "or 'chain'."
                )
            if exclude_same is not None and exclusion_block is not None:
                raise ValueError(
                    "The exclude_same argument to InterRDF cannot be used with exclusion_block."
                )
            name_to_attr = {
                "residue": "resindices",
                "segment": "segindices",
                "chain": "chainIDs",
            }
            self.exclude_same = name_to_attr.get(exclude_same)

            if self.norm not in ["rdf", "density", "none"]:
                raise ValueError(
                    f"'{self.norm}' is an invalid norm. "
                    "Use 'rdf', 'density' or 'none'."
                )

            self.backend = backend

        def _prepare(self):
            # Empty histogram to store the RDF
            count, edges = np.histogram([-1], **self.rdf_settings)
            count = count.astype(np.float64)
            count *= 0.0
            self.results.count = count
            self.results.edges = edges
            self.results.bins = 0.5 * (edges[:-1] + edges[1:])

            if self.norm == "rdf":
                # Cumulative volume for rdf normalization
                self.results.volume_cum = 0
            # Set the max range to filter the search radius
            self._maxrange = self.rdf_settings["range"][1]

        def _single_frame(self):
            pairs, dist = distances.capped_distance(
                self.g1.positions,
                self.g2.positions,
                self._maxrange,
                box=self._ts.dimensions,
                backend=self.backend,
            )
            # Maybe exclude same molecule distances
            if self._exclusion_block is not None:
                idxA = pairs[:, 0] // self._exclusion_block[0]
                idxB = pairs[:, 1] // self._exclusion_block[1]
                mask = np.where(idxA != idxB)[0]
                dist = dist[mask]

            if self.exclude_same is not None:
                # Ignore distances between atoms in the same attribute
                attr_ix_a = getattr(self.g1, self.exclude_same)[pairs[:, 0]]
                attr_ix_b = getattr(self.g2, self.exclude_same)[pairs[:, 1]]
                mask = np.where(attr_ix_a != attr_ix_b)[0]
                dist = dist[mask]

            count, _ = np.histogram(dist, **self.rdf_settings)
            self.results.count += count

            if self.norm == "rdf":
                self.results.volume_cum += self._ts.volume

        def _get_aggregator(self):
            return ResultsGroup(
                lookup={
                    "count": ResultsGroup.ndarray_sum,
                    "volume_cum": ResultsGroup.ndarray_sum,
                    "bins": ResultsGroup.ndarray_sum,
                    "edges": ResultsGroup.ndarray_mean,
                }
            )

        def _conclude(self):
            norm = self.n_frames
            if self.norm in ["rdf", "density"]:
                # Volume in each radial shell
                vols = np.power(self.results.edges, 3)
                norm *= 4 / 3 * np.pi * np.diff(vols)

            if self.norm == "rdf":
                # Number of each selection
                nA = len(self.g1)
                nB = len(self.g2)
                N = nA * nB

                # If we had exclusions, take these into account
                if self._exclusion_block:
                    xA, xB = self._exclusion_block
                    nblocks = nA / xA
                    N -= xA * xB * nblocks

                # Average number density
                self.volume_cum = self.results.volume_cum
                box_vol = self.volume_cum / self.n_frames
                norm *= N / box_vol

            self.results.rdf = self.results.count / norm

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

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

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

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

    class InterRDF_s(AnalysisBase):
        r"""Site-specific radial distribution function

        Calculates site-specific radial distribution
        functions. Instead of two groups of atoms it takes as input a list of
        pairs of AtomGroup, ``[[A, B], [C, D], ...]``. Given the same ``A`` and
        ``B`` to
        :class:`InterRDF_s`, the output will be a list of individual RDFs between
        ``A1`` and ``B1``, ``A1`` and ``B2``, ``A2`` and ``B1``, ``A2`` and ``B2``
        (and
        similarly for ``C`` and ``D``). These site-specific radial distribution
        functions are typically calculated if one is interested in the solvation
        shells of a ligand in a binding site or the solvation of specific residues
        in a protein.

        Parameters
        ----------
        u : Universe
            a Universe that contains atoms in `ags`

           .. deprecated:: 2.3.0
               This parameter is superflous and will be removed in
               MDAnalysis 3.0.0.

        ags : list
            a list of pairs of :class:`~MDAnalysis.core.groups.AtomGroup`
            instances
        nbins : int
            Number of bins in the histogram
        range : tuple or list
            The size of the RDF
        norm : str, {'rdf', 'density', 'none'}
            For 'rdf' calculate :math:`g_{ab}(r)`. For
            'density' the :ref:`single particle density<equation-nab>`
            :math:`n_{ab}(r)` is computed. 'none' computes the number of
            particles occurences in each spherical shell.

            .. versionadded:: 2.3.0

        density : bool
            `False`: calculate :math:`g_{ab}(r)`; `True`: calculate
            the true :ref:`single particle density<equation-nab>`
            :math:`n_{ab}(r)`. `density` overwrites the `norm` parameter.

            .. versionadded:: 1.0.1

                This keyword was available since 0.19.0 but was not
                documented. Furthermore, it had the opposite
                meaning. Since 1.0.1 it is officially supported as
                documented.

            .. deprecated:: 2.3.0
                Instead of `density=True` use `norm='density'`

        backend : {'serial', 'OpenMP', 'distopia'}, optional
            Keyword selecting the type of acceleration of the distance calculations.
            Default is 'serial'.

            .. versionadded:: 2.10.0

        Attributes
        ----------
        results.bins : numpy.ndarray
            :class:`numpy.ndarray` of the centers of the `nbins` histogram
            bins; all individual site-specific RDFs have the same bins.

            .. versionadded:: 2.0.0

        bins : numpy.ndarray
            Alias to the :attr:`results.bins` attribute.

            .. deprecated:: 2.0.0
                This attribute will be removed in 3.0.0.
                Use :attr:`results.bins` instead.

        results.edges : numpy.ndarray
            array of the ``nbins + 1`` edges of the histogram
            bins; all individual site-specific RDFs have the same bins.

            .. versionadded:: 2.0.0

        edges : numpy.ndarray
            Alias to the :attr:`results.edges` attribute.

            .. deprecated:: 2.0.0
                This attribute will be removed in 3.0.0.
                Use :attr:`results.edges` instead.

        results.rdf : list
            :class:`list` of the site-specific :ref:`radial distribution
            functions<equation-gab>` if `norm='rdf'` or :ref:`density
            functions<equation-nab>` for the :attr:`bins`
            if `norm='density'`. The list contains
            ``len(ags)`` entries. Each entry for the ``i``-th pair `[A, B]
            = ags[i]` in `ags` is a :class:`numpy.ndarray` with shape
            ``(len(A), len(B))``, i.e., a stack of RDFs. For example,
            ``results.rdf[i][0, 2]`` is the RDF between atoms ``A[0]``
            and ``B[2]``.

            .. versionadded:: 2.0.0

        rdf : list
            Alias to the :attr:`results.rdf` attribute.

            .. deprecated:: 2.0.0
                This attribute will be removed in 3.0.0.
                Use :attr:`results.rdf` instead.

        results.count : list
            :class:`list` of the site-specific radial histograms, i.e., the
            raw counts, for all :attr:`results.bins`. The data have the same
            structure as :attr:`results.rdf` except that the arrays contain
            the raw counts.

            .. versionadded:: 2.0.0

        count : list
            Alias to the :attr:`results.count` attribute.

            .. deprecated:: 2.0.0
                This attribute will be removed in 3.0.0.
                Use :attr:`results.count` instead.

        results.cdf : list
            :class:`list` of the site-specific :ref:`cumulative
            counts<equation-countab>`, for all :attr:`results.bins`. The data
            have the same structure as :attr:`results.rdf` except that the arrays
            contain the cumulative counts.

            This attribute only exists after :meth:`get_cdf` has been run.

            .. versionadded:: 2.0.0

        cdf : list
            Alias to the :attr:`results.cdf` attribute.

            .. deprecated:: 2.0.0
                This attribute will be removed in 3.0.0.
                Use :attr:`results.cdf` instead.

        Example
        -------
        First create the :class:`InterRDF_s` object, by supplying one Universe and
        one list of pairs of AtomGroups, then use the :meth:`~InterRDF_s.run`
        method::

          from MDAnalysisTests.datafiles import GRO_MEMPROT, XTC_MEMPROT
          u = mda.Universe(GRO_MEMPROT, XTC_MEMPROT)

          s1 = u.select_atoms('name ZND and resid 289')
          s2 = u.select_atoms('(name OD1 or name OD2) and resid 51 and sphzone 5.0 (resid 289)')
          s3 = u.select_atoms('name ZND and (resid 291 or resid 292)')
          s4 = u.select_atoms('(name OD1 or name OD2) and sphzone 5.0 (resid 291)')
          ags = [[s1, s2], [s3, s4]]

          rdf = InterRDF_s(u, ags)
          rdf.run()

        Results are available through the :attr:`results.bins`
        and :attr:`results.rdf` attributes::

          plt.plot(rdf.results.bins, rdf.results.rdf[0][0, 0])

        (Which plots the rdf between the first atom in ``s1`` and the first atom in
        ``s2``)

        To generate the *cumulative distribution function* (cdf) in the sense of
        "particles within radius :math:`r`", i.e., :math:`N_{ab}(r)`, use the
        :meth:`~InterRDF_s.get_cdf` method ::

          cdf = rdf.get_cdf()

        Results are available through the :attr:`results.cdf` attribute::

          plt.plot(rdf.results.bins, rdf.results.cdf[0][0, 0])

        (Which plots the cdf between the first atom in ``s1`` and the first atom in
        ``s2``)

        .. versionadded:: 0.19.0

        .. versionchanged:: 1.0.0
           Support for the `start`, `stop`, and `step` keywords has been
           removed. These should instead be passed to :meth:`InterRDF_s.run`.

        .. versionchanged:: 2.0.0
           Store results as attributes `bins`, `edges`, `rdf`, `count`
           and `cdf` of the `results` attribute
           of :class:`~MDAnalysis.analysis.AnalysisBase`.

        .. versionchanged:: 2.3.0
           Introduce `norm` and `exclusion_blocks` attributes.
        .. deprecated:: 2.3.0
           Instead of `density=True` use `norm='density'`
        .. deprecated:: 2.3.0
           The `universe` parameter is superflous.
        .. versionchanged:: 2.10.0
           Enabled **parallel execution** with the ``multiprocessing`` and ``dask``
           backends; use the new method :meth:`get_supported_backends` to see all
           supported backends.
        """

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

        _analysis_algorithm_is_parallelizable = True

        def _get_aggregator(self):
            return ResultsGroup(
                lookup={
                    "count": nested_array_sum,
                    "volume_cum": ResultsGroup.ndarray_sum,
                    "bins": ResultsGroup.ndarray_mean,
                    "edges": ResultsGroup.ndarray_mean,
                }
            )

        def __init__(
            self,
            u,
            ags,
            nbins=75,
            range=(0.0, 15.0),
            norm="rdf",
            density=False,
            backend="serial",
            **kwargs,
        ):
            super(InterRDF_s, self).__init__(
                ags[0][0].universe.trajectory, **kwargs
            )

            warnings.warn(
                "The `u` attribute is superflous and will be removed "
                "in MDAnalysis 3.0.0.",
                DeprecationWarning,
            )

            self.ags = ags
            self.norm = str(norm).lower()
            self.rdf_settings = {"bins": nbins, "range": range}

            if self.norm not in ["rdf", "density", "none"]:
                raise ValueError(
                    f"'{self.norm}' is an invalid norm. "
                    "Use 'rdf', 'density' or 'none'."
                )

            if density:
                warnings.warn(
                    "The `density` attribute was deprecated in "
                    "MDAnalysis 2.3.0 and will be removed in "
                    "MDAnalysis 3.0.0. Please use `norm=density` "
                    "instead.",
                    DeprecationWarning,
                )
                self.norm = "density"

            self.backend = backend

        def _prepare(self):
            count, edges = np.histogram([-1], **self.rdf_settings)
            self.results.count = [
                np.zeros((ag1.n_atoms, ag2.n_atoms, len(count)), dtype=np.float64)
                for ag1, ag2 in self.ags
            ]
            self.results.edges = edges
            self.results.bins = 0.5 * (edges[:-1] + edges[1:])

            if self.norm == "rdf":
                # Cumulative volume for rdf normalization
                self.results.volume_cum = 0
            self._maxrange = self.rdf_settings["range"][1]

        def _single_frame(self):
            for i, (ag1, ag2) in enumerate(self.ags):
                pairs, dist = distances.capped_distance(
                    ag1.positions,
                    ag2.positions,
                    self._maxrange,
                    box=self._ts.dimensions,
                    backend=self.backend,
                )

                # Fast manual implementation for distance histogram (equidistant bins)
                nbins = self.rdf_settings["bins"]
                rmin, rmax = self.rdf_settings["range"]
                bin_indices = (dist - rmin) * nbins / (rmax - rmin)
                bin_indices = bin_indices.astype(np.int64)
                counts = (bin_indices >= 0) & (bin_indices < nbins)
                counts = counts.astype(np.int64)
                idx1s = pairs[:, 0]
                idx2s = pairs[:, 1]
                self.results.count[i][idx1s, idx2s, bin_indices] += counts

            if self.norm == "rdf":
                self.results.volume_cum += self._ts.volume

        def _conclude(self):
            norm = self.n_frames
            if self.norm in ["rdf", "density"]:
                # Volume in each radial shell
                vols = np.power(self.results.edges, 3)
                norm *= 4 / 3 * np.pi * np.diff(vols)

            if self.norm == "rdf":
                # Average number density
                self.volume_cum = self.results.volume_cum
                norm *= 1 / (self.volume_cum / self.n_frames)

            # Empty lists to restore indices, RDF
            self.results.indices = []
            self.results.rdf = []

            for i, (ag1, ag2) in enumerate(self.ags):
                # Number of each selection
                self.results.indices.append([ag1.indices, ag2.indices])
                self.results.rdf.append(self.results.count[i] / norm)

        def get_cdf(self):
            r"""Calculate the cumulative counts for all sites.

            This is the :ref:`cumulative count<equation-countab>` within a given
            radius, i.e., :math:`N_{ab}(r)`.

            The result is returned and also stored in the attribute
            :attr:`results.cdf`.

            Returns
            -------
            cdf : list
                  list of arrays with the same structure as :attr:`results.rdf`
            """
            self.results.cdf = []

            for count in self.results.count:
                self.results.cdf.append(np.cumsum(count, axis=2) / self.n_frames)

            return self.results.cdf

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

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

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

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

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