1# -*- Mode: python; tab-width: 4; indent-tabs-mode:nil; -*-
2# vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4
3#
4# MDAnalysis --- https://www.mdanalysis.org
5# Copyright (c) 2006-2017 The MDAnalysis Development Team and contributors
6# (see the file AUTHORS for the full list of names)
7#
8# Released under the GNU Public Licence, v2 or any higher version
9#
10# Please cite your use of MDAnalysis in published work:
11#
12# R. J. Gowers, M. Linke, J. Barnoud, T. J. E. Reddy, M. N. Melo, S. L. Seyler,
13# D. L. Dotson, J. Domanski, S. Buchoux, I. M. Kenney, and O. Beckstein.
14# MDAnalysis: A Python package for the rapid analysis of molecular dynamics
15# simulations. In S. Benthall and S. Rostrup editors, Proceedings of the 15th
16# Python in Science Conference, pages 102-109, Austin, TX, 2016. SciPy.
17#
18# N. Michaud-Agrawal, E. J. Denning, T. B. Woolf, and O. Beckstein.
19# MDAnalysis: A Toolkit for the Analysis of Molecular Dynamics Simulations.
20# J. Comput. Chem. 32 (2011), 2319--2327, doi:10.1002/jcc.21787
21#
22#
23
24
25"""
26Distance analysis --- :mod:`MDAnalysis.analysis.distances`
27==========================================================
28
29This module provides functions to rapidly compute distances between
30atoms or groups of atoms.
31
32:func:`dist` and :func:`between` can take atom groups that do not even
33have to be from the same :class:`~MDAnalysis.core.universe.Universe`.
34
35See Also
36--------
37:mod:`MDAnalysis.lib.distances`
38"""
39
40from __future__ import absolute_import
41
42__all__ = ['distance_array', 'self_distance_array',
43           'contact_matrix', 'dist', 'between']
44
45import numpy as np
46import scipy.sparse
47
48from MDAnalysis.lib.distances import distance_array, self_distance_array
49from MDAnalysis.lib.c_distances import contact_matrix_no_pbc, contact_matrix_pbc
50from MDAnalysis.lib.NeighborSearch import AtomNeighborSearch
51from MDAnalysis.lib.distances import calc_bonds
52
53
54import warnings
55import logging
56logger = logging.getLogger("MDAnalysis.analysis.distances")
57
58
59def contact_matrix(coord, cutoff=15.0, returntype="numpy", box=None):
60    '''Calculates a matrix of contacts.
61
62    There is a fast, high-memory-usage version for small systems
63    (*returntype* = 'numpy'), and a slower, low-memory-usage version for
64    larger systems (*returntype* = 'sparse').
65
66    If *box* dimensions are passed then periodic boundary conditions
67    are applied.
68
69    Parameters
70    ---------
71    coord : array
72       Array of coordinates of shape ``(N, 3)`` and dtype float32.
73    cutoff : float, optional, default 15
74       Particles within `cutoff` are considered to form a contact.
75    returntype : string, optional, default "numpy"
76       Select how the contact matrix is returned.
77       * ``"numpy"``: return as an ``(N. N)`` :class:`numpy.ndarray`
78       * ``"sparse"``: return as a :class:`scipy.sparse.lil_matrix`
79    box : array-like or ``None``, optional, default ``None``
80       Simulation cell dimensions in the form of
81       :attr:`MDAnalysis.trajectory.base.Timestep.dimensions` when
82       periodic boundary conditions should be taken into account for
83       the calculation of contacts.
84
85    Returns
86    -------
87    array or sparse matrix
88       The contact matrix is returned in a format determined by the `returntype`
89       keyword.
90
91    See Also
92    --------
93    :mod:`MDAnalysis.analysis.contacts` for native contact analysis
94
95
96    .. versionchanged:: 0.11.0
97       Keyword *suppress_progmet* and *progress_meter_freq* were removed.
98    '''
99
100    if returntype == "numpy":
101        adj = (distance_array(coord, coord, box=box) < cutoff)
102        return adj
103    elif returntype == "sparse":
104        # Initialize square List of Lists matrix of dimensions equal to number
105        # of coordinates passed
106        sparse_contacts = scipy.sparse.lil_matrix((len(coord), len(coord)), dtype='bool')
107        if box is not None:
108            # with PBC
109            contact_matrix_pbc(coord, sparse_contacts, box, cutoff)
110        else:
111            # without PBC
112            contact_matrix_no_pbc(coord, sparse_contacts, cutoff)
113        return sparse_contacts
114
115
116def dist(A, B, offset=0, box=None):
117    """Return distance between atoms in two atom groups.
118
119    The distance is calculated atom-wise. The residue ids are also
120    returned because a typical use case is to look at CA distances
121    before and after an alignment. Using the `offset` keyword one can
122    also add a constant offset to the resids which facilitates
123    comparison with PDB numbering.
124
125    Arguments
126    ---------
127    A, B : AtomGroup
128       :class:`~MDAnalysis.core.groups.AtomGroup` with the
129       same number of atoms
130    offset : integer or tuple, optional, default 0
131       An integer `offset` is added to *resids_A* and *resids_B* (see
132       below) in order to produce PDB numbers.
133
134       If `offset` is :class:`tuple` then ``offset[0]`` is added to
135       *resids_A* and ``offset[1]`` to *resids_B*. Note that one can
136       actually supply numpy arrays of the same length as the atom
137       group so that an individual offset is added to each resid.
138
139    Returns
140    -------
141    resids_A : array
142        residue ids of the `A` group (possibly changed with `offset`)
143    resids_B : array
144       residue ids of the `B` group (possibly changed with `offset`)
145    distances : array
146       distances between the atoms
147    """
148
149    if A.atoms.n_atoms != B.atoms.n_atoms:
150        raise ValueError("AtomGroups A and B do not have the same number of atoms")
151    try:
152        off_A, off_B = offset
153    except (TypeError, ValueError):
154        off_A = off_B = int(offset)
155    residues_A = np.array(A.resids) + off_A
156    residues_B = np.array(B.resids) + off_B
157
158    d = calc_bonds(A.positions, B.positions, box)
159    return np.array([residues_A, residues_B, d])
160
161
162def between(group, A, B, distance):
163    """Return sub group of `group` that is within `distance` of both `A` and `B`
164
165    This function is not aware of periodic boundary conditions.
166
167    Can be used to find bridging waters or molecules in an interface.
168
169    Similar to "*group* and (AROUND *A* *distance* and AROUND *B* *distance*)".
170
171    Parameters
172    ----------
173    group : AtomGroup
174        Find members of `group` that are between `A` and `B`
175    A : AtomGroup
176    B : AtomGroup
177        `A` and `B` are the groups of atoms between which atoms in
178        `group` are searched for.  The function works is more
179        efficient if `group` is bigger than either `A` or `B`.
180    distance : float
181        maximum distance for an atom to be counted as in the vicinity of
182        `A` or `B`
183
184    Returns
185    -------
186    AtomGroup
187        :class:`~MDAnalysis.core.groups.AtomGroup` of atoms that
188        fulfill the criterion
189
190
191    .. versionadded: 0.7.5
192
193    """
194    ns_group = AtomNeighborSearch(group)
195    resA = set(ns_group.search(A, distance))
196    resB = set(ns_group.search(B, distance))
197    return sum(sorted(resB.intersection(resA)))
198