Merge pull request #4105 from xhgchen/simple-atomic-distance-analysis

Add simple atomic distance analysis (#3654)

package/CHANGELOG

@@ -52,6 +52,8 @@ Fixes
     (Issue #3336)
 
 Enhancements
+  * Add simple atomic distance analysis to `analysis.atomicdistances` with
+    new class `AtomicDistances` (Issue #3654)
   * Add kwarg `n_frames` to class method `empty()` in
   `MDAnalysis.core.universe`, enabling creation of a `Universe` with
   multiple frames from scratch (PR #4140)

package/MDAnalysis/analysis/atomicdistances.py

@@ -0,0 +1,178 @@
+# -*- 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")
+
+
+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)

package/doc/sphinx/source/documentation_pages/analysis/atomicdistances.rst

@@ -0,0 +1,2 @@
+.. automodule:: MDAnalysis.analysis.atomicdistances
+   :members:

package/doc/sphinx/source/documentation_pages/analysis_modules.rst

@@ -74,6 +74,7 @@ Distances and contacts
    analysis/align
    analysis/contacts
    analysis/distances
+   analysis/atomicdistances
    analysis/rms
    analysis/psa
    analysis/encore

testsuite/MDAnalysisTests/analysis/test_atomicdistances.py

@@ -0,0 +1,136 @@
+# -*- Mode: python; tab-width: 4; indent-tabs-mode:nil; coding:utf-8 -*-
+# vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4 fileencoding=utf-8
+#
+# 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
+#
+import pytest
+
+import MDAnalysis as mda
+
+import MDAnalysis.analysis.atomicdistances as ad
+from MDAnalysis.lib.distances import calc_bonds
+import MDAnalysis.transformations.boxdimensions as bd
+
+from numpy.testing import assert_allclose
+import numpy as np
+
+
+class TestAtomicDistances(object):
+    @staticmethod
+    @pytest.fixture()
+    def ad_u():
+        # Universe from scratch with 100 atoms
+        natoms = 100
+        u = mda.Universe.empty(natoms, trajectory=True)
+
+        # randomly generate 10 frames of x, y, z coordinates of 0 to 100
+        rng = np.random.default_rng()
+        coords = rng.integers(101, size=(10, natoms, 3))
+
+        # load frames and coordinates into Universe
+        u.load_new(coords)
+
+        # set box dimensions to dim of [lx, ly, lz, alpha, beta, gamma]
+        dim = np.array([80, 80, 80, 60, 60, 90])
+        transform = bd.set_dimensions(dim)
+        u.trajectory.add_transformations(transform)
+        return u
+
+    @staticmethod
+    @pytest.fixture()
+    def ad_ag1(ad_u):
+        return ad_u.atoms[10:20]
+
+    @staticmethod
+    @pytest.fixture()
+    def ad_ag2(ad_u):
+        return ad_u.atoms[70:80]
+
+    @staticmethod
+    @pytest.fixture()
+    def ad_ag3(ad_u):
+        # more atoms than expected (exception)
+        return ad_u.atoms[:7]
+
+    @staticmethod
+    @pytest.fixture()
+    def ad_ag4():
+        u_other = mda.Universe.empty(20, trajectory=True)
+        # different trajectory (exception)
+        return u_other.atoms[-10:]
+
+    @staticmethod
+    @pytest.fixture()
+    def expected_dist(ad_ag1, ad_ag2):
+        expected = np.zeros((len(ad_ag1.universe.trajectory),
+                             ad_ag1.atoms.n_atoms))
+
+        # calculate distances without PBCs using dist()
+        for i, ts in enumerate(ad_ag1.universe.trajectory):
+            expected[i] = calc_bonds(ad_ag1, ad_ag2)
+        return expected
+
+    @staticmethod
+    @pytest.fixture()
+    def expected_pbc_dist(ad_ag1, ad_ag2):
+        expected = np.zeros((len(ad_ag1.universe.trajectory),
+                             ad_ag1.atoms.n_atoms))
+
+        # calculate distances with PBCs using dist()
+        for i, ts in enumerate(ad_ag1.universe.trajectory):
+            expected[i] = calc_bonds(ad_ag1, ad_ag2, box=ad_ag1.dimensions)
+        return expected
+
+    def test_ad_exceptions(self, ad_ag1, ad_ag3, ad_ag4):
+        '''Test exceptions raised when number of atoms do not
+        match or AtomGroups come from different trajectories.'''
+
+        # number of atoms do not match
+        match_exp_atoms = "AtomGroups do not"
+        with pytest.raises(ValueError, match=match_exp_atoms):
+            ad_atoms = ad.AtomicDistances(ad_ag1, ad_ag3)
+
+        # AtomGroups come from different trajectories
+        match_exp_traj = "AtomGroups are not"
+        with pytest.raises(ValueError, match=match_exp_traj):
+            ad_traj = ad.AtomicDistances(ad_ag1, ad_ag4)
+
+    # only need to test that this class correctly applies distance calcs
+    # calc_bonds() is tested elsewhere
+    def test_ad_pairwise_dist(self, ad_ag1, ad_ag2,
+                              expected_dist):
+        '''Ensure that pairwise distances between atoms are
+        correctly calculated without PBCs.'''
+        pairwise_no_pbc = (ad.AtomicDistances(ad_ag1, ad_ag2,
+                                              pbc=False).run())
+        actual = pairwise_no_pbc.results
+
+        # compare with expected values from dist()
+        assert_allclose(actual, expected_dist)
+
+    def test_ad_pairwise_dist_pbc(self, ad_ag1, ad_ag2,
+                                  expected_pbc_dist):
+        '''Ensure that pairwise distances between atoms are
+        correctly calculated with PBCs.'''
+        pairwise_pbc = (ad.AtomicDistances(ad_ag1, ad_ag2).run())
+        actual = pairwise_pbc.results
+
+        # compare with expected values from dist()
+        assert_allclose(actual, expected_pbc_dist)