diff --git a/package/CHANGELOG b/package/CHANGELOG index 5b536f697c7..0da243a8564 100644 --- a/package/CHANGELOG +++ b/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) diff --git a/package/MDAnalysis/analysis/atomicdistances.py b/package/MDAnalysis/analysis/atomicdistances.py new file mode 100644 index 00000000000..e13c9543bd4 --- /dev/null +++ b/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) diff --git a/package/doc/sphinx/source/documentation_pages/analysis/atomicdistances.rst b/package/doc/sphinx/source/documentation_pages/analysis/atomicdistances.rst new file mode 100644 index 00000000000..42c884e7932 --- /dev/null +++ b/package/doc/sphinx/source/documentation_pages/analysis/atomicdistances.rst @@ -0,0 +1,2 @@ +.. automodule:: MDAnalysis.analysis.atomicdistances + :members: diff --git a/package/doc/sphinx/source/documentation_pages/analysis_modules.rst b/package/doc/sphinx/source/documentation_pages/analysis_modules.rst index c69745b5920..96ca0010dc3 100644 --- a/package/doc/sphinx/source/documentation_pages/analysis_modules.rst +++ b/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 diff --git a/testsuite/MDAnalysisTests/analysis/test_atomicdistances.py b/testsuite/MDAnalysisTests/analysis/test_atomicdistances.py new file mode 100644 index 00000000000..fb247a32ea9 --- /dev/null +++ b/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)