Source code for MDAnalysis.analysis.atomicdistances

# -*- Mode: python; tab-width: 4; indent-tabs-mode:nil; -*-
# 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 GNU Public Licence, v2 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"""
Simple atomic distance analysis --- :mod:`MDAnalysis.analysis.atomicdistances`
==============================================================================

:Author: Xu Hong Chen
:Year: 2023
:Copyright: GNU Public License v3

This module provides a class to efficiently compute distances between
two groups of atoms with an equal number of atoms over a trajectory.
Specifically, for two atom groups ``ag1`` and ``ag2``, it will return
the distances

.. math::
    |ag1[i] - ag2[i]|

for all :math:`i` from :math:`0` to `n_atoms` :math:`- 1`, where
`n_atoms` is the number of atoms in each atom group. By default,
this computation is done with periodic boundary conditions, but this
can be easily turned off. These distances are grouped by time step in
a NumPy array.

For more general functions on computing distances between atoms or
groups of atoms, please see :class:`MDAnalysis.analysis.distances`.

See Also
--------
:mod:`MDAnalysis.analysis.distances`
:mod:`MDAnalysis.lib.distances`


Basic usage
-----------

This example uses files from the MDAnalysis test suite
(:data:`~MDAnalysis.tests.datafiles.GRO` and
:data:`~MDAnalysis.tests.datafiles.XTC`). To get started, execute  ::

   >>> import MDAnalysis as mda
   >>> from MDAnalysis.tests.datafiles import GRO, XTC
   >>> import MDAnalysis.analysis.atomicdistances as ad

We will calculate the distances between an atom group of atoms 101-105
and an atom group of atoms 4001-4005 with periodic boundary conditions.
To select these atoms:

   >>> u = mda.Universe(GRO, XTC)
   >>> ag1 = u.atoms[100:105]
   >>> ag2 = u.atoms[4000:4005]

We can run the calculations using any variable of choice such as
``my_dists`` and access our results using ``my_dists.results``:

   >>> my_dists = ad.AtomicDistances(ag1, ag2).run()
   >>> my_dists.results
   array([[37.80813681, 33.2594864 , 34.93676414, 34.51183299, 34.96340209],
       [27.11746625, 31.19878079, 31.69439435, 32.63446126, 33.10451345],
       [23.27210749, 30.38714688, 32.48269361, 31.91444505, 31.84583838],
       [18.40607922, 39.21993135, 39.33468192, 41.0133789 , 39.46885946],
       [26.26006981, 37.9966713 , 39.14991106, 38.13423586, 38.95451427],
       [26.83845081, 34.66255735, 35.59335027, 34.8926705 , 34.27175056],
       [37.51994763, 38.12161091, 37.56481743, 36.8488121 , 35.75278065],
       [37.27275501, 37.7831456 , 35.74359073, 34.54893794, 34.76495816],
       [38.76272761, 41.31816555, 38.81588421, 39.82491432, 38.890219  ],
       [39.20012515, 40.00563374, 40.83857688, 38.77886735, 41.45775864]])

To do the computation without periodic boundary conditions, we can enter
the keyword argument ``pbc=False`` after ``ag2``. The result is different
in this case:

   >>> my_dists_nopbc = ad.AtomicDistances(ag1, ag2, pbc=False).run()
   >>> my_dists_nopbc.results
   array([[37.80813681, 33.2594864 , 34.93676414, 34.51183299, 34.96340209],
       [27.11746625, 31.19878079, 31.69439435, 32.63446126, 33.10451345],
       [23.27210749, 30.38714688, 32.482695  , 31.91444505, 31.84583838],
       [18.40607922, 39.21992825, 39.33468192, 41.0133757 , 39.46885946],
       [26.26006981, 37.99666906, 39.14990985, 38.13423708, 38.95451311],
       [26.83845081, 34.66255625, 35.59335027, 34.8926705 , 34.27174827],
       [51.86981409, 48.10347964, 48.39570072, 49.14423513, 50.44804292],
       [37.27275501, 37.7831456 , 35.74359073, 34.54893794, 34.76495816],
       [56.39657447, 41.31816555, 38.81588421, 39.82491432, 38.890219  ],
       [39.20012515, 40.00563374, 40.83857688, 38.77886735, 41.45775864]])

"""

import numpy as np

from MDAnalysis.lib.distances import calc_bonds


import warnings
import logging
from .base import AnalysisBase

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


[docs]class AtomicDistances(AnalysisBase): r"""Class to calculate atomic distances between two AtomGroups over a trajectory. Parameters ---------- ag1, ag2 : AtomGroup :class:`~MDAnalysis.core.groups.AtomGroup` with the same number of atoms pbc : bool, optional If ``True``, calculates atomic distances with periodic boundary conditions (PBCs). Setting `pbc` to ``False``, calculates atomic distances without considering PBCs. Defaults to ``True``. Attributes ---------- results : :class:`numpy.ndarray` The distances :math:`|ag1[i] - ag2[i]|` for all :math:`i` from :math:`0` to `n_atoms` :math:`- 1` for each frame over the trajectory. n_frames : int Number of frames included in the analysis. n_atoms : int Number of atoms in each atom group. .. versionadded:: 2.5.0 """ def __init__(self, ag1, ag2, pbc=True, **kwargs): # check ag1 and ag2 have the same number of atoms if ag1.atoms.n_atoms != ag2.atoms.n_atoms: raise ValueError("AtomGroups do not " "have the same number of atoms") # check ag1 and ag2 are from the same trajectory elif ag1.universe.trajectory != ag2.universe.trajectory: raise ValueError("AtomGroups are not " "from the same trajectory") super(AtomicDistances, self).__init__(ag1.universe.trajectory, **kwargs) self._ag1 = ag1 self._ag2 = ag2 self._pbc = pbc def _prepare(self): # initialize NumPy array of frames x distances for results self.results = np.zeros((self.n_frames, self._ag1.atoms.n_atoms)) def _single_frame(self): # if PBCs considered, get box size box = self._ag1.dimensions if self._pbc else None self.results[self._frame_index] = calc_bonds(self._ag1.positions, self._ag2.positions, box)