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"""Fast distance array computation --- :mod:`MDAnalysis.lib.distances` 25=================================================================== 26 27Fast C-routines to calculate arrays of distances or angles from coordinate 28arrays. Many of the functions also exist in parallel versions, which typically 29provide higher performance than the serial code. 30The boolean attribute `MDAnalysis.lib.distances.USED_OPENMP` can be checked to 31see if OpenMP was used in the compilation of MDAnalysis. 32 33Selection of acceleration ("backend") 34------------------------------------- 35 36All functions take the optional keyword `backend`, which determines the type of 37acceleration. Currently, the following choices are implemented (`backend` is 38case-insensitive): 39 40.. Table:: Available *backends* for accelerated distance functions. 41 42 ========== ========================= ====================================== 43 *backend* module description 44 ========== ========================= ====================================== 45 "serial" :mod:`c_distances` serial implementation in C/Cython 46 47 "OpenMP" :mod:`c_distances_openmp` parallel implementation in C/Cython 48 with OpenMP 49 ========== ========================= ====================================== 50 51.. versionadded:: 0.13.0 52 53Functions 54--------- 55.. autofunction:: distance_array 56.. autofunction:: self_distance_array 57.. autofunction:: capped_distance 58.. autofunction:: self_capped_distance 59.. autofunction:: calc_bonds 60.. autofunction:: calc_angles 61.. autofunction:: calc_dihedrals 62.. autofunction:: apply_PBC 63.. autofunction:: transform_RtoS 64.. autofunction:: transform_StoR 65.. autofunction:: augment_coordinates(coordinates, box, r) 66.. autofunction:: undo_augment(results, translation, nreal) 67""" 68from __future__ import division, absolute_import 69from six.moves import range 70 71import numpy as np 72from numpy.lib.utils import deprecate 73 74from .util import check_coords, check_box 75from .mdamath import triclinic_vectors, triclinic_box 76from ._augment import augment_coordinates, undo_augment 77from .nsgrid import FastNS 78 79# hack to select backend with backend=<backend> kwarg. Note that 80# the cython parallel code (prange) in parallel.distances is 81# independent from the OpenMP code 82import importlib 83_distances = {} 84_distances['serial'] = importlib.import_module(".c_distances", 85 package="MDAnalysis.lib") 86try: 87 _distances['openmp'] = importlib.import_module(".c_distances_openmp", 88 package="MDAnalysis.lib") 89except ImportError: 90 pass 91del importlib 92 93def _run(funcname, args=None, kwargs=None, backend="serial"): 94 """Helper function to select a backend function `funcname`.""" 95 args = args if args is not None else tuple() 96 kwargs = kwargs if kwargs is not None else dict() 97 backend = backend.lower() 98 try: 99 func = getattr(_distances[backend], funcname) 100 except KeyError: 101 raise ValueError("Function {0} not available with backend {1}; try one " 102 "of: {2}".format(funcname, backend, _distances.keys())) 103 return func(*args, **kwargs) 104 105# serial versions are always available (and are typically used within 106# the core and topology modules) 107from .c_distances import (calc_distance_array, 108 calc_distance_array_ortho, 109 calc_distance_array_triclinic, 110 calc_self_distance_array, 111 calc_self_distance_array_ortho, 112 calc_self_distance_array_triclinic, 113 coord_transform, 114 calc_bond_distance, 115 calc_bond_distance_ortho, 116 calc_bond_distance_triclinic, 117 calc_angle, 118 calc_angle_ortho, 119 calc_angle_triclinic, 120 calc_dihedral, 121 calc_dihedral_ortho, 122 calc_dihedral_triclinic, 123 ortho_pbc, 124 triclinic_pbc) 125 126from .c_distances_openmp import OPENMP_ENABLED as USED_OPENMP 127 128 129def _check_result_array(result, shape): 130 """Check if the result array is ok to use. 131 132 The `result` array must meet the following requirements: 133 * Must have a shape equal to `shape`. 134 * Its dtype must be ``numpy.float64``. 135 136 Paramaters 137 ---------- 138 result : numpy.ndarray or None 139 The result array to check. If `result` is `None``, a newly created 140 array of correct shape and dtype ``numpy.float64`` will be returned. 141 shape : tuple 142 The shape expected for the `result` array. 143 144 Returns 145 ------- 146 result : numpy.ndarray (``dtype=numpy.float64``, ``shape=shape``) 147 The input array or a newly created array if the input was ``None``. 148 149 Raises 150 ------ 151 ValueError 152 If `result` is of incorrect shape. 153 TypeError 154 If the dtype of `result` is not ``numpy.float64``. 155 """ 156 if result is None: 157 return np.zeros(shape, dtype=np.float64) 158 if result.shape != shape: 159 raise ValueError("Result array has incorrect shape, should be {0}, got " 160 "{1}.".format(shape, result.shape)) 161 if result.dtype != np.float64: 162 raise TypeError("Result array must be of type numpy.float64, got {}." 163 "".format(result.dtype)) 164# The following two lines would break a lot of tests. WHY?! 165# if not coords.flags['C_CONTIGUOUS']: 166# raise ValueError("{0} is not C-contiguous.".format(desc)) 167 return result 168 169 170@check_coords('reference', 'configuration', reduce_result_if_single=False, 171 check_lengths_match=False) 172def distance_array(reference, configuration, box=None, result=None, 173 backend="serial"): 174 """Calculate all possible distances between a reference set and another 175 configuration. 176 177 If there are ``n`` positions in `reference` and ``m`` positions in 178 `configuration`, a distance array of shape ``(n, m)`` will be computed. 179 180 If the optional argument `box` is supplied, the minimum image convention is 181 applied when calculating distances. Either orthogonal or triclinic boxes are 182 supported. 183 184 If a 2D numpy array of dtype ``numpy.float64`` with the shape ``(n, m)`` 185 is provided in `result`, then this preallocated array is filled. This can 186 speed up calculations. 187 188 Parameters 189 ---------- 190 reference : numpy.ndarray 191 Reference coordinate array of shape ``(3,)`` or ``(n, 3)`` (dtype is 192 arbitrary, will be converted to ``numpy.float32`` internally). 193 configuration : numpy.ndarray 194 Configuration coordinate array of shape ``(3,)`` or ``(m, 3)`` (dtype is 195 arbitrary, will be converted to ``numpy.float32`` internally). 196 box : array_like, optional 197 The unitcell dimensions of the system, which can be orthogonal or 198 triclinic and must be provided in the same format as returned by 199 :attr:`MDAnalysis.coordinates.base.Timestep.dimensions`:\n 200 ``[lx, ly, lz, alpha, beta, gamma]``. 201 result : numpy.ndarray, optional 202 Preallocated result array which must have the shape ``(n, m)`` and dtype 203 ``numpy.float64``. 204 Avoids creating the array which saves time when the function 205 is called repeatedly. 206 backend : {'serial', 'OpenMP'}, optional 207 Keyword selecting the type of acceleration. 208 209 Returns 210 ------- 211 d : numpy.ndarray (``dtype=numpy.float64``, ``shape=(n, m)``) 212 Array containing the distances ``d[i,j]`` between reference coordinates 213 ``i`` and configuration coordinates ``j``. 214 215 216 .. versionchanged:: 0.13.0 217 Added *backend* keyword. 218 .. versionchanged:: 0.19.0 219 Internal dtype conversion of input coordinates to ``numpy.float32``. 220 Now also accepts single coordinates as input. 221 """ 222 confnum = configuration.shape[0] 223 refnum = reference.shape[0] 224 225 distances = _check_result_array(result, (refnum, confnum)) 226 if len(distances) == 0: 227 return distances 228 229 if box is not None: 230 boxtype, box = check_box(box) 231 if boxtype == 'ortho': 232 _run("calc_distance_array_ortho", 233 args=(reference, configuration, box, distances), 234 backend=backend) 235 else: 236 _run("calc_distance_array_triclinic", 237 args=(reference, configuration, box, distances), 238 backend=backend) 239 else: 240 _run("calc_distance_array", 241 args=(reference, configuration, distances), 242 backend=backend) 243 244 return distances 245 246 247@check_coords('reference', reduce_result_if_single=False) 248def self_distance_array(reference, box=None, result=None, backend="serial"): 249 """Calculate all possible distances within a configuration `reference`. 250 251 If the optional argument `box` is supplied, the minimum image convention is 252 applied when calculating distances. Either orthogonal or triclinic boxes are 253 supported. 254 255 If a 1D numpy array of dtype ``numpy.float64`` with the shape 256 ``(n*(n-1)/2,)`` is provided in `result`, then this preallocated array is 257 filled. This can speed up calculations. 258 259 Parameters 260 ---------- 261 reference : numpy.ndarray 262 Reference coordinate array of shape ``(3,)`` or ``(n, 3)`` (dtype is 263 arbitrary, will be converted to ``numpy.float32`` internally). 264 box : array_like, optional 265 The unitcell dimensions of the system, which can be orthogonal or 266 triclinic and must be provided in the same format as returned by 267 :attr:`MDAnalysis.coordinates.base.Timestep.dimensions`:\n 268 ``[lx, ly, lz, alpha, beta, gamma]``. 269 result : numpy.ndarray, optional 270 Preallocated result array which must have the shape ``(n*(n-1)/2,)`` and 271 dtype ``numpy.float64``. Avoids creating the array which saves time when 272 the function is called repeatedly. 273 backend : {'serial', 'OpenMP'}, optional 274 Keyword selecting the type of acceleration. 275 276 Returns 277 ------- 278 d : numpy.ndarray (``dtype=numpy.float64``, ``shape=(n*(n-1)/2,)``) 279 Array containing the distances ``dist[i,j]`` between reference 280 coordinates ``i`` and ``j`` at position ``d[k]``. Loop through ``d``: 281 282 .. code-block:: python 283 284 for i in range(n): 285 for j in range(i + 1, n): 286 k += 1 287 dist[i, j] = d[k] 288 289 290 .. versionchanged:: 0.13.0 291 Added *backend* keyword. 292 .. versionchanged:: 0.19.0 293 Internal dtype conversion of input coordinates to ``numpy.float32``. 294 """ 295 refnum = reference.shape[0] 296 distnum = refnum * (refnum - 1) // 2 297 298 distances = _check_result_array(result, (distnum,)) 299 if len(distances) == 0: 300 return distances 301 302 if box is not None: 303 boxtype, box = check_box(box) 304 if boxtype == 'ortho': 305 _run("calc_self_distance_array_ortho", 306 args=(reference, box, distances), 307 backend=backend) 308 else: 309 _run("calc_self_distance_array_triclinic", 310 args=(reference, box, distances), 311 backend=backend) 312 else: 313 _run("calc_self_distance_array", 314 args=(reference, distances), 315 backend=backend) 316 317 return distances 318 319 320def capped_distance(reference, configuration, max_cutoff, min_cutoff=None, 321 box=None, method=None, return_distances=True): 322 """Calculates pairs of indices corresponding to entries in the `reference` 323 and `configuration` arrays which are separated by a distance lying within 324 the specified cutoff(s). Optionally, these distances can be returned as 325 well. 326 327 If the optional argument `box` is supplied, the minimum image convention is 328 applied when calculating distances. Either orthogonal or triclinic boxes are 329 supported. 330 331 An automatic guessing of the optimal method to calculate the distances is 332 included in the function. An optional keyword for the method is also 333 provided. Users can enforce a particular method with this functionality. 334 Currently brute force, grid search, and periodic KDtree methods are 335 implemented. 336 337 Parameters 338 ----------- 339 reference : numpy.ndarray 340 Reference coordinate array with shape ``(3,)`` or ``(n, 3)``. 341 configuration : numpy.ndarray 342 Configuration coordinate array with shape ``(3,)`` or ``(m, 3)``. 343 max_cutoff : float 344 Maximum cutoff distance between the reference and configuration. 345 min_cutoff : float, optional 346 Minimum cutoff distance between reference and configuration. 347 box : array_like, optional 348 The unitcell dimensions of the system, which can be orthogonal or 349 triclinic and must be provided in the same format as returned by 350 :attr:`MDAnalysis.coordinates.base.Timestep.dimensions`:\n 351 ``[lx, ly, lz, alpha, beta, gamma]``. 352 method : {'bruteforce', 'nsgrid', 'pkdtree'}, optional 353 Keyword to override the automatic guessing of the employed search 354 method. 355 return_distances : bool, optional 356 If set to ``True``, distances will also be returned. 357 358 Returns 359 ------- 360 pairs : numpy.ndarray (``dtype=numpy.int64``, ``shape=(n_pairs, 2)``) 361 Pairs of indices, corresponding to coordinates in the `reference` and 362 `configuration` arrays such that the distance between them lies within 363 the interval (`min_cutoff`, `max_cutoff`]. 364 Each row in `pairs` is an index pair ``[i, j]`` corresponding to the 365 ``i``-th coordinate in `reference` and the ``j``-th coordinate in 366 `configuration`. 367 distances : numpy.ndarray (``dtype=numpy.float64``, ``shape=(n_pairs,)``), optional 368 Distances corresponding to each pair of indices. Only returned if 369 `return_distances` is ``True``. ``distances[k]`` corresponds to the 370 ``k``-th pair returned in `pairs` and gives the distance between the 371 coordinates ``reference[pairs[k, 0]]`` and 372 ``configuration[pairs[k, 1]]``. 373 374 .. code-block:: python 375 376 pairs, distances = capped_distances(reference, configuration, 377 max_cutoff, return_distances=True) 378 for k, [i, j] in enumerate(pairs): 379 coord1 = reference[i] 380 coord2 = configuration[j] 381 distance = distances[k] 382 383 Note 384 ----- 385 Currently supports brute force, grid-based, and periodic KDtree search 386 methods. 387 388 See Also 389 -------- 390 distance_array 391 MDAnalysis.lib.pkdtree.PeriodicKDTree.search 392 MDAnalysis.lib.nsgrid.FastNS.search 393 """ 394 if box is not None: 395 box = np.asarray(box, dtype=np.float32) 396 if box.shape[0] != 6: 397 raise ValueError("Box Argument is of incompatible type. The " 398 "dimension should be either None or of the form " 399 "[lx, ly, lz, alpha, beta, gamma]") 400 method = _determine_method(reference, configuration, max_cutoff, 401 min_cutoff=min_cutoff, box=box, method=method) 402 return method(reference, configuration, max_cutoff, min_cutoff=min_cutoff, 403 box=box, return_distances=return_distances) 404 405 406def _determine_method(reference, configuration, max_cutoff, min_cutoff=None, 407 box=None, method=None): 408 """Guesses the fastest method for capped distance calculations based on the 409 size of the coordinate sets and the relative size of the target volume. 410 411 Parameters 412 ---------- 413 reference : numpy.ndarray 414 Reference coordinate array with shape ``(3,)`` or ``(n, 3)``. 415 configuration : numpy.ndarray 416 Configuration coordinate array with shape ``(3,)`` or ``(m, 3)``. 417 max_cutoff : float 418 Maximum cutoff distance between `reference` and `configuration` 419 coordinates. 420 min_cutoff : float, optional 421 Minimum cutoff distance between `reference` and `configuration` 422 coordinates. 423 box : numpy.ndarray 424 The unitcell dimensions of the system, which can be orthogonal or 425 triclinic and must be provided in the same format as returned by 426 :attr:`MDAnalysis.coordinates.base.Timestep.dimensions`:\n 427 ``[lx, ly, lz, alpha, beta, gamma]``. 428 method : {'bruteforce', 'nsgrid', 'pkdtree'}, optional 429 Keyword to override the automatic guessing of the employed search 430 method. 431 432 Returns 433 ------- 434 function : callable 435 The function implementing the guessed (or deliberatly chosen) method. 436 """ 437 methods = {'bruteforce': _bruteforce_capped, 438 'pkdtree': _pkdtree_capped, 439 'nsgrid': _nsgrid_capped} 440 441 if method is not None: 442 return methods[method.lower()] 443 444 if len(reference) < 10 or len(configuration) < 10: 445 return methods['bruteforce'] 446 elif len(reference) * len(configuration) >= 1e8: 447 # CAUTION : for large datasets, shouldnt go into 'bruteforce' 448 # in any case. Arbitrary number, but can be characterized 449 return methods['nsgrid'] 450 else: 451 if box is None: 452 min_dim = np.array([reference.min(axis=0), 453 configuration.min(axis=0)]) 454 max_dim = np.array([reference.max(axis=0), 455 configuration.max(axis=0)]) 456 size = max_dim.max(axis=0) - min_dim.min(axis=0) 457 elif np.all(box[3:] == 90.0): 458 size = box[:3] 459 else: 460 tribox = triclinic_vectors(box) 461 size = tribox.max(axis=0) - tribox.min(axis=0) 462 if np.any(max_cutoff > 0.3*size): 463 return methods['bruteforce'] 464 else: 465 return methods['nsgrid'] 466 467 468@check_coords('reference', 'configuration', enforce_copy=False, 469 reduce_result_if_single=False, check_lengths_match=False) 470def _bruteforce_capped(reference, configuration, max_cutoff, min_cutoff=None, 471 box=None, return_distances=True): 472 """Capped distance evaluations using a brute force method. 473 474 Computes and returns an array containing pairs of indices corresponding to 475 entries in the `reference` and `configuration` arrays which are separated by 476 a distance lying within the specified cutoff(s). Employs naive distance 477 computations (brute force) to find relevant distances. 478 479 Optionally, these distances can be returned as well. 480 481 If the optional argument `box` is supplied, the minimum image convention is 482 applied when calculating distances. Either orthogonal or triclinic boxes are 483 supported. 484 485 Parameters 486 ---------- 487 reference : numpy.ndarray 488 Reference coordinate array with shape ``(3,)`` or ``(n, 3)`` (dtype will 489 be converted to ``numpy.float32`` internally). 490 configuration : array 491 Configuration coordinate array with shape ``(3,)`` or ``(m, 3)`` (dtype 492 will be converted to ``numpy.float32`` internally). 493 max_cutoff : float 494 Maximum cutoff distance between `reference` and `configuration` 495 coordinates. 496 min_cutoff : float, optional 497 Minimum cutoff distance between `reference` and `configuration` 498 coordinates. 499 box : numpy.ndarray, optional 500 The unitcell dimensions of the system, which can be orthogonal or 501 triclinic and must be provided in the same format as returned by 502 :attr:`MDAnalysis.coordinates.base.Timestep.dimensions`:\n 503 ``[lx, ly, lz, alpha, beta, gamma]``. 504 return_distances : bool, optional 505 If set to ``True``, distances will also be returned. 506 507 Returns 508 ------- 509 pairs : numpy.ndarray (``dtype=numpy.int64``, ``shape=(n_pairs, 2)``) 510 Pairs of indices, corresponding to coordinates in the `reference` and 511 `configuration` arrays such that the distance between them lies within 512 the interval (`min_cutoff`, `max_cutoff`]. 513 Each row in `pairs` is an index pair ``[i, j]`` corresponding to the 514 ``i``-th coordinate in `reference` and the ``j``-th coordinate in 515 `configuration`. 516 distances : numpy.ndarray (``dtype=numpy.float64``, ``shape=(n_pairs,)``), optional 517 Distances corresponding to each pair of indices. Only returned if 518 `return_distances` is ``True``. ``distances[k]`` corresponds to the 519 ``k``-th pair returned in `pairs` and gives the distance between the 520 coordinates ``reference[pairs[k, 0]]`` and 521 ``configuration[pairs[k, 1]]``. 522 """ 523 # Default return values (will be overwritten only if pairs are found): 524 pairs = np.empty((0, 2), dtype=np.int64) 525 distances = np.empty((0,), dtype=np.float64) 526 527 if len(reference) > 0 and len(configuration) > 0: 528 _distances = distance_array(reference, configuration, box=box) 529 if min_cutoff is not None: 530 mask = np.where((_distances <= max_cutoff) & \ 531 (_distances > min_cutoff)) 532 else: 533 mask = np.where((_distances <= max_cutoff)) 534 if mask[0].size > 0: 535 pairs = np.c_[mask[0], mask[1]] 536 if return_distances: 537 distances = _distances[mask] 538 539 if return_distances: 540 return pairs, distances 541 else: 542 return pairs 543 544 545@check_coords('reference', 'configuration', enforce_copy=False, 546 reduce_result_if_single=False, check_lengths_match=False) 547def _pkdtree_capped(reference, configuration, max_cutoff, min_cutoff=None, 548 box=None, return_distances=True): 549 """Capped distance evaluations using a KDtree method. 550 551 Computes and returns an array containing pairs of indices corresponding to 552 entries in the `reference` and `configuration` arrays which are separated by 553 a distance lying within the specified cutoff(s). Employs a (periodic) KDtree 554 algorithm to find relevant distances. 555 556 Optionally, these distances can be returned as well. 557 558 If the optional argument `box` is supplied, the minimum image convention is 559 applied when calculating distances. Either orthogonal or triclinic boxes are 560 supported. 561 562 Parameters 563 ---------- 564 reference : numpy.ndarray 565 Reference coordinate array with shape ``(3,)`` or ``(n, 3)`` (dtype will 566 be converted to ``numpy.float32`` internally). 567 configuration : array 568 Configuration coordinate array with shape ``(3,)`` or ``(m, 3)`` (dtype 569 will be converted to ``numpy.float32`` internally). 570 max_cutoff : float 571 Maximum cutoff distance between `reference` and `configuration` 572 coordinates. 573 min_cutoff : float, optional 574 Minimum cutoff distance between `reference` and `configuration` 575 coordinates. 576 box : numpy.ndarray, optional 577 The unitcell dimensions of the system, which can be orthogonal or 578 triclinic and must be provided in the same format as returned by 579 :attr:`MDAnalysis.coordinates.base.Timestep.dimensions`:\n 580 ``[lx, ly, lz, alpha, beta, gamma]``. 581 return_distances : bool, optional 582 If set to ``True``, distances will also be returned. 583 584 Returns 585 ------- 586 pairs : numpy.ndarray (``dtype=numpy.int64``, ``shape=(n_pairs, 2)``) 587 Pairs of indices, corresponding to coordinates in the `reference` and 588 `configuration` arrays such that the distance between them lies within 589 the interval (`min_cutoff`, `max_cutoff`]. 590 Each row in `pairs` is an index pair ``[i, j]`` corresponding to the 591 ``i``-th coordinate in `reference` and the ``j``-th coordinate in 592 `configuration`. 593 distances : numpy.ndarray (``dtype=numpy.float64``, ``shape=(n_pairs,)``), optional 594 Distances corresponding to each pair of indices. Only returned if 595 `return_distances` is ``True``. ``distances[k]`` corresponds to the 596 ``k``-th pair returned in `pairs` and gives the distance between the 597 coordinates ``reference[pairs[k, 0]]`` and 598 ``configuration[pairs[k, 1]]``. 599 """ 600 from .pkdtree import PeriodicKDTree # must be here to avoid circular import 601 602 # Default return values (will be overwritten only if pairs are found): 603 pairs = np.empty((0, 2), dtype=np.int64) 604 distances = np.empty((0,), dtype=np.float64) 605 606 if len(reference) > 0 and len(configuration) > 0: 607 kdtree = PeriodicKDTree(box=box) 608 cut = max_cutoff if box is not None else None 609 kdtree.set_coords(configuration, cutoff=cut) 610 _pairs = kdtree.search_tree(reference, max_cutoff) 611 if _pairs.size > 0: 612 pairs = _pairs 613 if (return_distances or (min_cutoff is not None)): 614 refA, refB = pairs[:, 0], pairs[:, 1] 615 distances = calc_bonds(reference[refA], configuration[refB], 616 box=box) 617 if min_cutoff is not None: 618 mask = np.where(distances > min_cutoff) 619 pairs, distances = pairs[mask], distances[mask] 620 621 if return_distances: 622 return pairs, distances 623 else: 624 return pairs 625 626 627@check_coords('reference', 'configuration', enforce_copy=False, 628 reduce_result_if_single=False, check_lengths_match=False) 629def _nsgrid_capped(reference, configuration, max_cutoff, min_cutoff=None, 630 box=None, return_distances=True): 631 """Capped distance evaluations using a grid-based search method. 632 633 Computes and returns an array containing pairs of indices corresponding to 634 entries in the `reference` and `configuration` arrays which are separated by 635 a distance lying within the specified cutoff(s). Employs a grid-based search 636 algorithm to find relevant distances. 637 638 Optionally, these distances can be returned as well. 639 640 If the optional argument `box` is supplied, the minimum image convention is 641 applied when calculating distances. Either orthogonal or triclinic boxes are 642 supported. 643 644 Parameters 645 ---------- 646 reference : numpy.ndarray 647 Reference coordinate array with shape ``(3,)`` or ``(n, 3)`` (dtype will 648 be converted to ``numpy.float32`` internally). 649 configuration : array 650 Configuration coordinate array with shape ``(3,)`` or ``(m, 3)`` (dtype 651 will be converted to ``numpy.float32`` internally). 652 max_cutoff : float 653 Maximum cutoff distance between `reference` and `configuration` 654 coordinates. 655 min_cutoff : float, optional 656 Minimum cutoff distance between `reference` and `configuration` 657 coordinates. 658 box : numpy.ndarray (``dtype=numpy.float64``, ``shape=(n_pairs,)``), optional 659 The unitcell dimensions of the system, which can be orthogonal or 660 triclinic and must be provided in the same format as returned by 661 :attr:`MDAnalysis.coordinates.base.Timestep.dimensions`:\n 662 ``[lx, ly, lz, alpha, beta, gamma]``. 663 return_distances : bool, optional 664 If set to ``True``, distances will also be returned. 665 666 Returns 667 ------- 668 pairs : numpy.ndarray (``dtype=numpy.int64``, ``shape=(n_pairs, 2)``) 669 Pairs of indices, corresponding to coordinates in the `reference` and 670 `configuration` arrays such that the distance between them lies within 671 the interval (`min_cutoff`, `max_cutoff`]. 672 Each row in `pairs` is an index pair ``[i, j]`` corresponding to the 673 ``i``-th coordinate in `reference` and the ``j``-th coordinate in 674 `configuration`. 675 distances : numpy.ndarray, optional 676 Distances corresponding to each pair of indices. Only returned if 677 `return_distances` is ``True``. ``distances[k]`` corresponds to the 678 ``k``-th pair returned in `pairs` and gives the distance between the 679 coordinates ``reference[pairs[k, 0]]`` and 680 ``configuration[pairs[k, 1]]``. 681 """ 682 # Default return values (will be overwritten only if pairs are found): 683 pairs = np.empty((0, 2), dtype=np.int64) 684 distances = np.empty((0,), dtype=np.float64) 685 686 if len(reference) > 0 and len(configuration) > 0: 687 if box is None: 688 # create a pseudobox 689 # define the max range 690 # and supply the pseudobox 691 # along with only one set of coordinates 692 pseudobox = np.zeros(6, dtype=np.float32) 693 all_coords = np.concatenate([reference, configuration]) 694 lmax = all_coords.max(axis=0) 695 lmin = all_coords.min(axis=0) 696 # Using maximum dimension as the box size 697 boxsize = (lmax-lmin).max() 698 # to avoid failures for very close particles but with 699 # larger cutoff 700 boxsize = np.maximum(boxsize, 2 * max_cutoff) 701 pseudobox[:3] = 1.2 * boxsize 702 pseudobox[3:] = 90. 703 shiftref, shiftconf = reference.copy(), configuration.copy() 704 # Extra padding near the origin 705 shiftref -= lmin - 0.1*boxsize 706 shiftconf -= lmin - 0.1*boxsize 707 gridsearch = FastNS(max_cutoff, shiftconf, box=pseudobox, pbc=False) 708 results = gridsearch.search(shiftref) 709 else: 710 gridsearch = FastNS(max_cutoff, configuration, box=box) 711 results = gridsearch.search(reference) 712 713 pairs = results.get_pairs() 714 if return_distances or (min_cutoff is not None): 715 distances = results.get_pair_distances() 716 if min_cutoff is not None: 717 idx = distances > min_cutoff 718 pairs, distances = pairs[idx], distances[idx] 719 720 if return_distances: 721 return pairs, distances 722 else: 723 return pairs 724 725 726def self_capped_distance(reference, max_cutoff, min_cutoff=None, box=None, 727 method=None): 728 """Calculates pairs of indices corresponding to entries in the `reference` 729 array which are separated by a distance lying within the specified 730 cutoff(s). The respective distances are returned as well. 731 732 If the optional argument `box` is supplied, the minimum image convention is 733 applied when calculating distances. Either orthogonal or triclinic boxes are 734 supported. 735 736 An automatic guessing of the optimal method to calculate the distances is 737 included in the function. An optional keyword for the method is also 738 provided. Users can enforce a particular method with this functionality. 739 Currently brute force, grid search, and periodic KDtree methods are 740 implemented. 741 742 Parameters 743 ----------- 744 reference : numpy.ndarray 745 Reference coordinate array with shape ``(3,)`` or ``(n, 3)``. 746 max_cutoff : float 747 Maximum cutoff distance between `reference` coordinates. 748 min_cutoff : float, optional 749 Minimum cutoff distance between `reference` coordinates. 750 box : array_like, optional 751 The unitcell dimensions of the system, which can be orthogonal or 752 triclinic and must be provided in the same format as returned by 753 :attr:`MDAnalysis.coordinates.base.Timestep.dimensions`:\n 754 ``[lx, ly, lz, alpha, beta, gamma]``. 755 method : {'bruteforce', 'nsgrid', 'pkdtree'}, optional 756 Keyword to override the automatic guessing of the employed search 757 method. 758 759 Returns 760 ------- 761 pairs : numpy.ndarray (``dtype=numpy.int64``, ``shape=(n_pairs, 2)``) 762 Pairs of indices, corresponding to coordinates in the `reference` array 763 such that the distance between them lies within the interval 764 (`min_cutoff`, `max_cutoff`]. 765 Each row in `pairs` is an index pair ``[i, j]`` corresponding to the 766 ``i``-th and the ``j``-th coordinate in `reference`. 767 distances : numpy.ndarray (``dtype=numpy.float64``, ``shape=(n_pairs,)``) 768 Distances corresponding to each pair of indices. ``distances[k]`` 769 corresponds to the ``k``-th pair returned in `pairs` and gives the 770 distance between the coordinates ``reference[pairs[k, 0]]`` and 771 ``reference[pairs[k, 1]]``. 772 773 .. code-block:: python 774 775 pairs, distances = self_capped_distances(reference, max_cutoff) 776 for k, [i, j] in enumerate(pairs): 777 coord1 = reference[i] 778 coord2 = reference[j] 779 distance = distances[k] 780 781 782 Note 783 ----- 784 Currently supports brute force, grid-based, and periodic KDtree search 785 methods. 786 787 See Also 788 -------- 789 self_distance_array 790 MDAnalysis.lib.pkdtree.PeriodicKDTree.search 791 MDAnalysis.lib.nsgrid.FastNS.self_search 792 """ 793 if box is not None: 794 box = np.asarray(box, dtype=np.float32) 795 if box.shape[0] != 6: 796 raise ValueError("Box Argument is of incompatible type. The " 797 "dimension should be either None or of the form " 798 "[lx, ly, lz, alpha, beta, gamma]") 799 method = _determine_method_self(reference, max_cutoff, 800 min_cutoff=min_cutoff, 801 box=box, method=method) 802 return method(reference, max_cutoff, min_cutoff=min_cutoff, box=box) 803 804 805def _determine_method_self(reference, max_cutoff, min_cutoff=None, box=None, 806 method=None): 807 """Guesses the fastest method for capped distance calculations based on the 808 size of the `reference` coordinate set and the relative size of the target 809 volume. 810 811 Parameters 812 ---------- 813 reference : numpy.ndarray 814 Reference coordinate array with shape ``(3,)`` or ``(n, 3)``. 815 max_cutoff : float 816 Maximum cutoff distance between `reference` coordinates. 817 min_cutoff : float, optional 818 Minimum cutoff distance between `reference` coordinates. 819 box : numpy.ndarray 820 The unitcell dimensions of the system, which can be orthogonal or 821 triclinic and must be provided in the same format as returned by 822 :attr:`MDAnalysis.coordinates.base.Timestep.dimensions`:\n 823 ``[lx, ly, lz, alpha, beta, gamma]``. 824 method : {'bruteforce', 'nsgrid', 'pkdtree'}, optional 825 Keyword to override the automatic guessing of the employed search 826 method. 827 828 Returns 829 ------- 830 function : callable 831 The function implementing the guessed (or deliberatly chosen) method. 832 """ 833 methods = {'bruteforce': _bruteforce_capped_self, 834 'pkdtree': _pkdtree_capped_self, 835 'nsgrid': _nsgrid_capped_self} 836 837 if method is not None: 838 return methods[method.lower()] 839 840 if len(reference) < 100: 841 return methods['bruteforce'] 842 843 if box is None: 844 min_dim = np.array([reference.min(axis=0)]) 845 max_dim = np.array([reference.max(axis=0)]) 846 size = max_dim.max(axis=0) - min_dim.min(axis=0) 847 elif np.all(box[3:] == 90.0): 848 size = box[:3] 849 else: 850 tribox = triclinic_vectors(box) 851 size = tribox.max(axis=0) - tribox.min(axis=0) 852 853 if max_cutoff < 0.03*size.min(): 854 return methods['pkdtree'] 855 else: 856 return methods['nsgrid'] 857 858 859@check_coords('reference', enforce_copy=False, reduce_result_if_single=False) 860def _bruteforce_capped_self(reference, max_cutoff, min_cutoff=None, box=None): 861 """Capped distance evaluations using a brute force method. 862 863 Computes and returns an array containing pairs of indices corresponding to 864 entries in the `reference` array which are separated by a distance lying 865 within the specified cutoff(s). Employs naive distance computations (brute 866 force) to find relevant distances. These distances are returned as well. 867 868 If the optional argument `box` is supplied, the minimum image convention is 869 applied when calculating distances. Either orthogonal or triclinic boxes are 870 supported. 871 872 Parameters 873 ---------- 874 reference : numpy.ndarray 875 Reference coordinate array with shape ``(3,)`` or ``(n, 3)`` (dtype will 876 be converted to ``numpy.float32`` internally). 877 max_cutoff : float 878 Maximum cutoff distance between `reference` coordinates. 879 min_cutoff : float, optional 880 Minimum cutoff distance between `reference` coordinates. 881 box : numpy.ndarray, optional 882 The unitcell dimensions of the system, which can be orthogonal or 883 triclinic and must be provided in the same format as returned by 884 :attr:`MDAnalysis.coordinates.base.Timestep.dimensions`:\n 885 ``[lx, ly, lz, alpha, beta, gamma]``. 886 887 Returns 888 ------- 889 pairs : numpy.ndarray (``dtype=numpy.int64``, ``shape=(n_pairs, 2)``) 890 Pairs of indices, corresponding to coordinates in the `reference` array 891 such that the distance between them lies within the interval 892 (`min_cutoff`, `max_cutoff`]. 893 Each row in `pairs` is an index pair ``[i, j]`` corresponding to the 894 ``i``-th and the ``j``-th coordinate in `reference`. 895 distances : numpy.ndarray (``dtype=numpy.float64``, ``shape=(n_pairs,)``) 896 Distances corresponding to each pair of indices. ``distances[k]`` 897 corresponds to the ``k``-th pair returned in `pairs` and gives the 898 distance between the coordinates ``reference[pairs[k, 0]]`` and 899 ``reference[pairs[k, 1]]``. 900 """ 901 # Default return values (will be overwritten only if pairs are found): 902 pairs = np.empty((0, 2), dtype=np.int64) 903 distances = np.empty((0,), dtype=np.float64) 904 905 N = len(reference) 906 # We're searching within a single coordinate set, so we need at least two 907 # coordinates to find distances between them. 908 if N > 1: 909 distvec = self_distance_array(reference, box=box) 910 dist = np.full((N, N), np.finfo(np.float64).max, dtype=np.float64) 911 dist[np.triu_indices(N, 1)] = distvec 912 913 if min_cutoff is not None: 914 mask = np.where((dist <= max_cutoff) & (dist > min_cutoff)) 915 else: 916 mask = np.where((dist <= max_cutoff)) 917 918 if mask[0].size > 0: 919 pairs = np.c_[mask[0], mask[1]] 920 distances = dist[mask] 921 return pairs, distances 922 923 924@check_coords('reference', enforce_copy=False, reduce_result_if_single=False) 925def _pkdtree_capped_self(reference, max_cutoff, min_cutoff=None, box=None): 926 """Capped distance evaluations using a KDtree method. 927 928 Computes and returns an array containing pairs of indices corresponding to 929 entries in the `reference` array which are separated by a distance lying 930 within the specified cutoff(s). Employs a (periodic) KDtree algorithm to 931 find relevant distances. These distances are returned as well. 932 933 If the optional argument `box` is supplied, the minimum image convention is 934 applied when calculating distances. Either orthogonal or triclinic boxes are 935 supported. 936 937 Parameters 938 ---------- 939 reference : numpy.ndarray 940 Reference coordinate array with shape ``(3,)`` or ``(n, 3)`` (dtype will 941 be converted to ``numpy.float32`` internally). 942 max_cutoff : float 943 Maximum cutoff distance between `reference` coordinates. 944 min_cutoff : float, optional 945 Minimum cutoff distance between `reference` coordinates. 946 box : numpy.ndarray, optional 947 The unitcell dimensions of the system, which can be orthogonal or 948 triclinic and must be provided in the same format as returned by 949 :attr:`MDAnalysis.coordinates.base.Timestep.dimensions`:\n 950 ``[lx, ly, lz, alpha, beta, gamma]``. 951 952 Returns 953 ------- 954 pairs : numpy.ndarray (``dtype=numpy.int64``, ``shape=(n_pairs, 2)``) 955 Pairs of indices, corresponding to coordinates in the `reference` array 956 such that the distance between them lies within the interval 957 (`min_cutoff`, `max_cutoff`]. 958 Each row in `pairs` is an index pair ``[i, j]`` corresponding to the 959 ``i``-th and the ``j``-th coordinate in `reference`. 960 distances : numpy.ndarray (``dtype=numpy.float64``, ``shape=(n_pairs,)``) 961 Distances corresponding to each pair of indices. ``distances[k]`` 962 corresponds to the ``k``-th pair returned in `pairs` and gives the 963 distance between the coordinates ``reference[pairs[k, 0]]`` and 964 ``reference[pairs[k, 1]]``. 965 """ 966 from .pkdtree import PeriodicKDTree # must be here to avoid circular import 967 968 # Default return values (will be overwritten only if pairs are found): 969 pairs = np.empty((0, 2), dtype=np.int64) 970 distances = np.empty((0,), dtype=np.float64) 971 972 # We're searching within a single coordinate set, so we need at least two 973 # coordinates to find distances between them. 974 if len(reference) > 1: 975 kdtree = PeriodicKDTree(box=box) 976 cut = max_cutoff if box is not None else None 977 kdtree.set_coords(reference, cutoff=cut) 978 _pairs = kdtree.search_pairs(max_cutoff) 979 if _pairs.size > 0: 980 pairs = _pairs 981 refA, refB = pairs[:, 0], pairs[:, 1] 982 distances = calc_bonds(reference[refA], reference[refB], box=box) 983 if min_cutoff is not None: 984 idx = distances > min_cutoff 985 pairs, distances = pairs[idx], distances[idx] 986 return pairs, distances 987 988@check_coords('reference', enforce_copy=False, reduce_result_if_single=False) 989def _nsgrid_capped_self(reference, max_cutoff, min_cutoff=None, box=None): 990 """Capped distance evaluations using a grid-based search method. 991 992 Computes and returns an array containing pairs of indices corresponding to 993 entries in the `reference` array which are separated by a distance lying 994 within the specified cutoff(s). Employs a grid-based search algorithm to 995 find relevant distances. These distances are returned as well. 996 997 If the optional argument `box` is supplied, the minimum image convention is 998 applied when calculating distances. Either orthogonal or triclinic boxes are 999 supported. 1000 1001 Parameters 1002 ---------- 1003 reference : numpy.ndarray 1004 Reference coordinate array with shape ``(3,)`` or ``(n, 3)`` (dtype will 1005 be converted to ``numpy.float32`` internally). 1006 max_cutoff : float 1007 Maximum cutoff distance between `reference` coordinates. 1008 min_cutoff : float, optional 1009 Minimum cutoff distance between `reference` coordinates. 1010 box : numpy.ndarray, optional 1011 The unitcell dimensions of the system, which can be orthogonal or 1012 triclinic and must be provided in the same format as returned by 1013 :attr:`MDAnalysis.coordinates.base.Timestep.dimensions`:\n 1014 ``[lx, ly, lz, alpha, beta, gamma]``. 1015 1016 Returns 1017 ------- 1018 pairs : numpy.ndarray (``dtype=numpy.int64``, ``shape=(n_pairs, 2)``) 1019 Pairs of indices, corresponding to coordinates in the `reference` array 1020 such that the distance between them lies within the interval 1021 (`min_cutoff`, `max_cutoff`]. 1022 Each row in `pairs` is an index pair ``[i, j]`` corresponding to the 1023 ``i``-th and the ``j``-th coordinate in `reference`. 1024 distances : numpy.ndarray (``dtype=numpy.float64``, ``shape=(n_pairs,)``) 1025 Distances corresponding to each pair of indices. ``distances[k]`` 1026 corresponds to the ``k``-th pair returned in `pairs` and gives the 1027 distance between the coordinates ``reference[pairs[k, 0]]`` and 1028 ``reference[pairs[k, 1]]``. 1029 """ 1030 # Default return values (will be overwritten only if pairs are found): 1031 pairs = np.empty((0, 2), dtype=np.int64) 1032 distances = np.empty((0,), dtype=np.float64) 1033 1034 # We're searching within a single coordinate set, so we need at least two 1035 # coordinates to find distances between them. 1036 if len(reference) > 1: 1037 if box is None: 1038 # create a pseudobox 1039 # define the max range 1040 # and supply the pseudobox 1041 # along with only one set of coordinates 1042 pseudobox = np.zeros(6, dtype=np.float32) 1043 lmax = reference.max(axis=0) 1044 lmin = reference.min(axis=0) 1045 # Using maximum dimension as the box size 1046 boxsize = (lmax-lmin).max() 1047 # to avoid failures of very close particles 1048 # but with larger cutoff 1049 if boxsize < 2*max_cutoff: 1050 # just enough box size so that NSGrid doesnot fails 1051 sizefactor = 2.2*max_cutoff/boxsize 1052 else: 1053 sizefactor = 1.2 1054 pseudobox[:3] = sizefactor*boxsize 1055 pseudobox[3:] = 90. 1056 shiftref = reference.copy() 1057 # Extra padding near the origin 1058 shiftref -= lmin - 0.1*boxsize 1059 gridsearch = FastNS(max_cutoff, shiftref, box=pseudobox, pbc=False) 1060 results = gridsearch.self_search() 1061 else: 1062 gridsearch = FastNS(max_cutoff, reference, box=box) 1063 results = gridsearch.self_search() 1064 1065 pairs = results.get_pairs()[::2, :] 1066 distances = results.get_pair_distances()[::2] 1067 if min_cutoff is not None: 1068 idx = distances > min_cutoff 1069 pairs, distances = pairs[idx], distances[idx] 1070 1071 return pairs, distances 1072 1073 1074@check_coords('coords') 1075def transform_RtoS(coords, box, backend="serial"): 1076 """Transform an array of coordinates from real space to S space (a.k.a. 1077 lambda space) 1078 1079 S space represents fractional space within the unit cell for this system. 1080 1081 Reciprocal operation to :meth:`transform_StoR`. 1082 1083 Parameters 1084 ---------- 1085 coords : numpy.ndarray 1086 A ``(3,)`` or ``(n, 3)`` array of coordinates (dtype is arbitrary, will 1087 be converted to ``numpy.float32`` internally). 1088 box : numpy.ndarray 1089 The unitcell dimensions of the system, which can be orthogonal or 1090 triclinic and must be provided in the same format as returned by 1091 :attr:`MDAnalysis.coordinates.base.Timestep.dimensions`:\n 1092 ``[lx, ly, lz, alpha, beta, gamma]``. 1093 backend : {'serial', 'OpenMP'}, optional 1094 Keyword selecting the type of acceleration. 1095 1096 Returns 1097 ------- 1098 newcoords : numpy.ndarray (``dtype=numpy.float32``, ``shape=coords.shape``) 1099 An array containing fractional coordiantes. 1100 1101 1102 .. versionchanged:: 0.13.0 1103 Added *backend* keyword. 1104 .. versionchanged:: 0.19.0 1105 Internal dtype conversion of input coordinates to ``numpy.float32``. 1106 Now also accepts (and, likewise, returns) a single coordinate. 1107 """ 1108 if len(coords) == 0: 1109 return coords 1110 boxtype, box = check_box(box) 1111 if boxtype == 'ortho': 1112 box = np.diag(box) 1113 1114 # Create inverse matrix of box 1115 # need order C here 1116 inv = np.array(np.linalg.inv(box), dtype=np.float32, order='C') 1117 1118 _run("coord_transform", args=(coords, inv), backend=backend) 1119 1120 return coords 1121 1122 1123@check_coords('coords') 1124def transform_StoR(coords, box, backend="serial"): 1125 """Transform an array of coordinates from S space into real space. 1126 1127 S space represents fractional space within the unit cell for this system. 1128 1129 Reciprocal operation to :meth:`transform_RtoS` 1130 1131 Parameters 1132 ---------- 1133 coords : numpy.ndarray 1134 A ``(3,)`` or ``(n, 3)`` array of coordinates (dtype is arbitrary, will 1135 be converted to ``numpy.float32`` internally). 1136 box : numpy.ndarray 1137 The unitcell dimensions of the system, which can be orthogonal or 1138 triclinic and must be provided in the same format as returned by 1139 :attr:`MDAnalysis.coordinates.base.Timestep.dimensions`:\n 1140 ``[lx, ly, lz, alpha, beta, gamma]``. 1141 backend : {'serial', 'OpenMP'}, optional 1142 Keyword selecting the type of acceleration. 1143 1144 Returns 1145 ------- 1146 newcoords : numpy.ndarray (``dtype=numpy.float32``, ``shape=coords.shape``) 1147 An array containing real space coordiantes. 1148 1149 1150 .. versionchanged:: 0.13.0 1151 Added *backend* keyword. 1152 .. versionchanged:: 0.19.0 1153 Internal dtype conversion of input coordinates to ``numpy.float32``. 1154 Now also accepts (and, likewise, returns) a single coordinate. 1155 """ 1156 if len(coords) == 0: 1157 return coords 1158 boxtype, box = check_box(box) 1159 if boxtype == 'ortho': 1160 box = np.diag(box) 1161 1162 _run("coord_transform", args=(coords, box), backend=backend) 1163 return coords 1164 1165 1166@check_coords('coords1', 'coords2') 1167def calc_bonds(coords1, coords2, box=None, result=None, backend="serial"): 1168 """Calculates the bond lengths between pairs of atom positions from the two 1169 coordinate arrays `coords1` and `coords2`, which must contain the same 1170 number of coordinates. ``coords1[i]`` and ``coords2[i]`` represent the 1171 positions of atoms connected by the ``i``-th bond. If single coordinates are 1172 supplied, a single distance will be returned. 1173 1174 In comparison to :meth:`distance_array` and :meth:`self_distance_array`, 1175 which calculate distances between all possible combinations of coordinates, 1176 :meth:`calc_bonds` only calculates distances between pairs of coordinates, 1177 similar to:: 1178 1179 numpy.linalg.norm(a - b) for a, b in zip(coords1, coords2) 1180 1181 If the optional argument `box` is supplied, the minimum image convention is 1182 applied when calculating distances. Either orthogonal or triclinic boxes are 1183 supported. 1184 1185 If a numpy array of dtype ``numpy.float64`` with shape ``(n,)`` (for ``n`` 1186 coordinate pairs) is provided in `result`, then this preallocated array is 1187 filled. This can speed up calculations. 1188 1189 Parameters 1190 ---------- 1191 coords1 : numpy.ndarray 1192 Coordinate array of shape ``(3,)`` or ``(n, 3)`` for one half of a 1193 single or ``n`` bonds, respectively (dtype is arbitrary, will be 1194 converted to ``numpy.float32`` internally). 1195 coords2 : numpy.ndarray 1196 Coordinate array of shape ``(3,)`` or ``(n, 3)`` for the other half of 1197 a single or ``n`` bonds, respectively (dtype is arbitrary, will be 1198 converted to ``numpy.float32`` internally). 1199 box : numpy.ndarray, optional 1200 The unitcell dimensions of the system, which can be orthogonal or 1201 triclinic and must be provided in the same format as returned by 1202 :attr:`MDAnalysis.coordinates.base.Timestep.dimensions`:\n 1203 ``[lx, ly, lz, alpha, beta, gamma]``. 1204 result : numpy.ndarray, optional 1205 Preallocated result array of dtype ``numpy.float64`` and shape ``(n,)`` 1206 (for ``n`` coordinate pairs). Avoids recreating the array in repeated 1207 function calls. 1208 backend : {'serial', 'OpenMP'}, optional 1209 Keyword selecting the type of acceleration. 1210 1211 Returns 1212 ------- 1213 bondlengths : numpy.ndarray (``dtype=numpy.float64``, ``shape=(n,)``) or numpy.float64 1214 Array containing the bond lengths between each pair of coordinates. If 1215 two single coordinates were supplied, their distance is returned as a 1216 single number instead of an array. 1217 1218 1219 .. versionadded:: 0.8 1220 .. versionchanged:: 0.13.0 1221 Added *backend* keyword. 1222 .. versionchanged:: 0.19.0 1223 Internal dtype conversion of input coordinates to ``numpy.float32``. 1224 Now also accepts single coordinates as input. 1225 """ 1226 numatom = coords1.shape[0] 1227 bondlengths = _check_result_array(result, (numatom,)) 1228 1229 if numatom > 0: 1230 if box is not None: 1231 boxtype, box = check_box(box) 1232 if boxtype == 'ortho': 1233 _run("calc_bond_distance_ortho", 1234 args=(coords1, coords2, box, bondlengths), 1235 backend=backend) 1236 else: 1237 _run("calc_bond_distance_triclinic", 1238 args=(coords1, coords2, box, bondlengths), 1239 backend=backend) 1240 else: 1241 _run("calc_bond_distance", 1242 args=(coords1, coords2, bondlengths), 1243 backend=backend) 1244 1245 return bondlengths 1246 1247 1248@check_coords('coords1', 'coords2', 'coords3') 1249def calc_angles(coords1, coords2, coords3, box=None, result=None, 1250 backend="serial"): 1251 """Calculates the angles formed between triplets of atom positions from the 1252 three coordinate arrays `coords1`, `coords2`, and `coords3`. All coordinate 1253 arrays must contain the same number of coordinates. 1254 1255 The coordinates in `coords2` represent the apices of the angles:: 1256 1257 2---3 1258 / 1259 1 1260 1261 Configurations where the angle is undefined (e.g., when coordinates 1 or 3 1262 of a triplet coincide with coordinate 2) result in a value of **zero** for 1263 that angle. 1264 1265 If the optional argument `box` is supplied, periodic boundaries are taken 1266 into account when constructing the connecting vectors between coordinates, 1267 i.e., the minimum image convention is applied for the vectors forming the 1268 angles. Either orthogonal or triclinic boxes are supported. 1269 1270 If a numpy array of dtype ``numpy.float64`` with shape ``(n,)`` (for ``n`` 1271 coordinate triplets) is provided in `result`, then this preallocated array 1272 is filled. This can speed up calculations. 1273 1274 Parameters 1275 ---------- 1276 coords1 : numpy.ndarray 1277 Array of shape ``(3,)`` or ``(n, 3)`` containing the coordinates of one 1278 side of a single or ``n`` angles, respectively (dtype is arbitrary, will 1279 be converted to ``numpy.float32`` internally) 1280 coords2 : numpy.ndarray 1281 Array of shape ``(3,)`` or ``(n, 3)`` containing the coordinates of the 1282 apices of a single or ``n`` angles, respectively (dtype is arbitrary, 1283 will be converted to ``numpy.float32`` internally) 1284 coords3 : numpy.ndarray 1285 Array of shape ``(3,)`` or ``(n, 3)`` containing the coordinates of the 1286 other side of a single or ``n`` angles, respectively (dtype is 1287 arbitrary, will be converted to ``numpy.float32`` internally) 1288 box : numpy.ndarray, optional 1289 The unitcell dimensions of the system, which can be orthogonal or 1290 triclinic and must be provided in the same format as returned by 1291 :attr:`MDAnalysis.coordinates.base.Timestep.dimensions`:\n 1292 ``[lx, ly, lz, alpha, beta, gamma]``. 1293 result : numpy.ndarray, optional 1294 Preallocated result array of dtype ``numpy.float64`` and shape ``(n,)`` 1295 (for ``n`` coordinate triplets). Avoids recreating the array in repeated 1296 function calls. 1297 backend : {'serial', 'OpenMP'}, optional 1298 Keyword selecting the type of acceleration. 1299 1300 Returns 1301 ------- 1302 angles : numpy.ndarray (``dtype=numpy.float64``, ``shape=(n,)``) or numpy.float64 1303 Array containing the angles between each triplet of coordinates. Values 1304 are returned in radians (rad). If three single coordinates were 1305 supplied, the angle is returned as a single number instead of an array. 1306 1307 1308 .. versionadded:: 0.8 1309 .. versionchanged:: 0.9.0 1310 Added optional box argument to account for periodic boundaries in 1311 calculation 1312 .. versionchanged:: 0.13.0 1313 Added *backend* keyword. 1314 .. versionchanged:: 0.19.0 1315 Internal dtype conversion of input coordinates to ``numpy.float32``. 1316 Now also accepts single coordinates as input. 1317 """ 1318 numatom = coords1.shape[0] 1319 angles = _check_result_array(result, (numatom,)) 1320 1321 if numatom > 0: 1322 if box is not None: 1323 boxtype, box = check_box(box) 1324 if boxtype == 'ortho': 1325 _run("calc_angle_ortho", 1326 args=(coords1, coords2, coords3, box, angles), 1327 backend=backend) 1328 else: 1329 _run("calc_angle_triclinic", 1330 args=(coords1, coords2, coords3, box, angles), 1331 backend=backend) 1332 else: 1333 _run("calc_angle", 1334 args=(coords1, coords2, coords3, angles), 1335 backend=backend) 1336 1337 return angles 1338 1339 1340@check_coords('coords1', 'coords2', 'coords3', 'coords4') 1341def calc_dihedrals(coords1, coords2, coords3, coords4, box=None, result=None, 1342 backend="serial"): 1343 """Calculates the dihedral angles formed between quadruplets of positions 1344 from the four coordinate arrays `coords1`, `coords2`, `coords3`, and 1345 `coords4`, which must contain the same number of coordinates. 1346 1347 The dihedral angle formed by a quadruplet of positions (1,2,3,4) is 1348 calculated around the axis connecting positions 2 and 3 (i.e., the angle 1349 between the planes spanned by positions (1,2,3) and (2,3,4)):: 1350 1351 4 1352 | 1353 2-----3 1354 / 1355 1 1356 1357 If all coordinates lie in the same plane, the cis configuration corresponds 1358 to a dihedral angle of zero, and the trans configuration to :math:`\pi` 1359 radians (180 degrees). Configurations where the dihedral angle is undefined 1360 (e.g., when all coordinates lie on the same straight line) result in a value 1361 of ``nan`` (not a number) for that dihedral. 1362 1363 If the optional argument `box` is supplied, periodic boundaries are taken 1364 into account when constructing the connecting vectors between coordinates, 1365 i.e., the minimum image convention is applied for the vectors forming the 1366 dihedral angles. Either orthogonal or triclinic boxes are supported. 1367 1368 If a numpy array of dtype ``numpy.float64`` with shape ``(n,)`` (for ``n`` 1369 coordinate quadruplets) is provided in `result` then this preallocated array 1370 is filled. This can speed up calculations. 1371 1372 Parameters 1373 ---------- 1374 coords1 : numpy.ndarray 1375 Coordinate array of shape ``(3,)`` or ``(n, 3)`` containing the 1st 1376 positions in dihedrals (dtype is arbitrary, will be converted to 1377 ``numpy.float32`` internally) 1378 coords2 : numpy.ndarray 1379 Coordinate array of shape ``(3,)`` or ``(n, 3)`` containing the 2nd 1380 positions in dihedrals (dtype is arbitrary, will be converted to 1381 ``numpy.float32`` internally) 1382 coords3 : numpy.ndarray 1383 Coordinate array of shape ``(3,)`` or ``(n, 3)`` containing the 3rd 1384 positions in dihedrals (dtype is arbitrary, will be converted to 1385 ``numpy.float32`` internally) 1386 coords4 : numpy.ndarray 1387 Coordinate array of shape ``(3,)`` or ``(n, 3)`` containing the 4th 1388 positions in dihedrals (dtype is arbitrary, will be converted to 1389 ``numpy.float32`` internally) 1390 box : numpy.ndarray, optional 1391 The unitcell dimensions of the system, which can be orthogonal or 1392 triclinic and must be provided in the same format as returned by 1393 :attr:`MDAnalysis.coordinates.base.Timestep.dimensions`:\n 1394 ``[lx, ly, lz, alpha, beta, gamma]``. 1395 result : numpy.ndarray, optional 1396 Preallocated result array of dtype ``numpy.float64`` and shape ``(n,)`` 1397 (for ``n`` coordinate quadruplets). Avoids recreating the array in 1398 repeated function calls. 1399 backend : {'serial', 'OpenMP'}, optional 1400 Keyword selecting the type of acceleration. 1401 1402 Returns 1403 ------- 1404 dihedrals : numpy.ndarray (``dtype=numpy.float64``, ``shape=(n,)``) or numpy.float64 1405 Array containing the dihedral angles formed by each quadruplet of 1406 coordinates. Values are returned in radians (rad). If four single 1407 coordinates were supplied, the dihedral angle is returned as a single 1408 number instead of an array. 1409 1410 1411 .. versionadded:: 0.8 1412 .. versionchanged:: 0.9.0 1413 Added optional box argument to account for periodic boundaries in 1414 calculation 1415 .. versionchanged:: 0.11.0 1416 Renamed from calc_torsions to calc_dihedrals 1417 .. versionchanged:: 0.13.0 1418 Added *backend* keyword. 1419 .. versionchanged:: 0.19.0 1420 Internal dtype conversion of input coordinates to ``numpy.float32``. 1421 Now also accepts single coordinates as input. 1422 """ 1423 numatom = coords1.shape[0] 1424 dihedrals = _check_result_array(result, (numatom,)) 1425 1426 if numatom > 0: 1427 if box is not None: 1428 boxtype, box = check_box(box) 1429 if boxtype == 'ortho': 1430 _run("calc_dihedral_ortho", 1431 args=(coords1, coords2, coords3, coords4, box, dihedrals), 1432 backend=backend) 1433 else: 1434 _run("calc_dihedral_triclinic", 1435 args=(coords1, coords2, coords3, coords4, box, dihedrals), 1436 backend=backend) 1437 else: 1438 _run("calc_dihedral", 1439 args=(coords1, coords2, coords3, coords4, dihedrals), 1440 backend=backend) 1441 1442 return dihedrals 1443 1444 1445@check_coords('coords') 1446def apply_PBC(coords, box, backend="serial"): 1447 """Moves coordinates into the primary unit cell. 1448 1449 Parameters 1450 ---------- 1451 coords : numpy.ndarray 1452 Coordinate array of shape ``(3,)`` or ``(n, 3)`` (dtype is arbitrary, 1453 will be converted to ``numpy.float32`` internally). 1454 box : numpy.ndarray 1455 The unitcell dimensions of the system, which can be orthogonal or 1456 triclinic and must be provided in the same format as returned by 1457 :attr:`MDAnalysis.coordinates.base.Timestep.dimensions`:\n 1458 ``[lx, ly, lz, alpha, beta, gamma]``. 1459 backend : {'serial', 'OpenMP'}, optional 1460 Keyword selecting the type of acceleration. 1461 1462 Returns 1463 ------- 1464 newcoords : numpy.ndarray (``dtype=numpy.float32``, ``shape=coords.shape``) 1465 Array containing coordinates that all lie within the primary unit cell 1466 as defined by `box`. 1467 1468 1469 .. versionadded:: 0.8 1470 .. versionchanged:: 0.13.0 1471 Added *backend* keyword. 1472 .. versionchanged:: 0.19.0 1473 Internal dtype conversion of input coordinates to ``numpy.float32``. 1474 Now also accepts (and, likewise, returns) single coordinates. 1475 """ 1476 if len(coords) == 0: 1477 return coords 1478 boxtype, box = check_box(box) 1479 if boxtype == 'ortho': 1480 box_inv = box ** (-1) 1481 _run("ortho_pbc", args=(coords, box, box_inv), backend=backend) 1482 else: 1483 box_inv = np.diagonal(box) ** (-1) 1484 _run("triclinic_pbc", args=(coords, box, box_inv), backend=backend) 1485 1486 return coords 1487