Source code for MDAnalysis.core.topologyobjects

# -*- Mode: python; tab-width: 4; indent-tabs-mode:nil; coding:utf-8 -*-
# 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
#

"""
Core Topology Objects --- :mod:`MDAnalysis.core.topologyobjects`
================================================================

The building blocks for MDAnalysis' description of topology

"""
from __future__ import print_function, absolute_import, division

from six.moves import zip
from six import raise_from
import numbers
import numpy as np
import functools

from ..lib import mdamath
from ..lib.util import cached
from ..lib import util
from ..lib import distances


[docs]@functools.total_ordering class TopologyObject(object): """Base class for all Topology items. Defines the behaviour by which Bonds/Angles/etc in MDAnalysis should behave. .. versionadded:: 0.9.0 .. versionchanged:: 0.10.0 All TopologyObject now keep track of if they were guessed or not via the ``is_guessed`` managed property. .. versionadded:: 0.11.0 Added the `value` method to return the size of the object """ __slots__ = ("_ix", "_u", "btype", "_bondtype", "_guessed", "order") def __init__(self, ix, universe, type=None, guessed=False, order=None): """Create a topology object Parameters ---------- ix : numpy array indices of the Atoms universe : MDAnalysis.Universe type : optional Type of the bond guessed : optional If the Bond is guessed """ self._ix = ix self._u = universe self._bondtype = type self._guessed = guessed self.order = order @property def atoms(self): """Atoms within this Bond""" return self._u.atoms[self._ix] @property def indices(self): """Tuple of indices describing this object .. versionadded:: 0.10.0 """ return self._ix @property def universe(self): return self._u @property def type(self): """Type of the bond as a tuple Note ---- When comparing types, it is important to consider the reverse of the type too, i.e.:: a.type == b.type or a.type == b.type[::-1] """ if self._bondtype is not None: return self._bondtype else: return tuple(self.atoms.types) @property def is_guessed(self): return bool(self._guessed) def __hash__(self): return hash((self._u, tuple(self.indices))) def __repr__(self): indices = (self.indices if self.indices[0] < self.indices[-1] else self.indices[::-1]) return "<{cname} between: {conts}>".format( cname=self.__class__.__name__, conts=", ".join([ "Atom {0}".format(i) for i in indices])) def __contains__(self, other): """Check whether an atom is in this :class:`TopologyObject`""" return other in self.atoms def __eq__(self, other): """Check whether two bonds have identical contents""" if not self.universe == other.universe: return False return (np.array_equal(self.indices, other.indices) or np.array_equal(self.indices[::-1], other.indices)) def __ne__(self, other): return not self == other def __lt__(self, other): return tuple(self.indices) < tuple(other.indices) def __getitem__(self, item): """Can retrieve a given Atom from within""" return self.atoms[item] def __iter__(self): return iter(self.atoms) def __len__(self): return len(self._ix)
[docs]class Bond(TopologyObject): """A bond between two :class:`~MDAnalysis.core.groups.Atom` instances. Two :class:`Bond` instances can be compared with the ``==`` and ``!=`` operators. A bond is equal to another if the same atom numbers are connected and they have the same bond order. The ordering of the two atom numbers is ignored as is the fact that a bond was guessed. The presence of a particular atom can also be queried:: >>> Atom in Bond will return either ``True`` or ``False``. .. versionchanged:: 0.9.0 Now a subclass of :class:`TopologyObject`. Changed class to use :attr:`__slots__` and stores atoms in :attr:`atoms` attribute. """ btype = 'bond'
[docs] def partner(self, atom): """Bond.partner(Atom) Returns ------- the other :class:`~MDAnalysis.core.groups.Atom` in this bond """ if atom == self.atoms[0]: return self.atoms[1] elif atom == self.atoms[1]: return self.atoms[0] else: raise ValueError("Unrecognised Atom")
[docs] def length(self, pbc=True): """Length of the bond. .. versionchanged:: 0.11.0 Added pbc keyword .. versionchanged:: 0.19.0 Changed default of pbc to True """ box = self.universe.dimensions if pbc else None return distances.calc_bonds(self[0].position, self[1].position, box)
value = length
[docs]class Angle(TopologyObject): """An angle between three :class:`~MDAnalysis.core.groups.Atom` instances. Atom 2 is the apex of the angle .. versionadded:: 0.8 .. versionchanged:: 0.9.0 Now a subclass of :class:`TopologyObject`; now uses :attr:`__slots__` and stores atoms in :attr:`atoms` attribute """ btype = 'angle'
[docs] def angle(self, pbc=True): """Returns the angle in degrees of this Angle. Angle between atoms 0 and 2 with apex at 1:: 2 / / 1------0 Note ---- The numerical precision is typically not better than 4 decimals (and is only tested to 3 decimals). .. versionadded:: 0.9.0 .. versionchanged:: 0.17.0 Fixed angles close to 180 giving NaN .. versionchanged:: 0.19.0 Added pbc keyword, default True """ box = self.universe.dimensions if pbc else None return np.rad2deg(distances.calc_angles( self[0].position, self[1].position, self[2].position, box))
value = angle
[docs]class Dihedral(TopologyObject): """Dihedral (dihedral angle) between four :class:`~MDAnalysis.core.groups.Atom` instances. The dihedral is defined as the angle between the planes formed by Atoms (1, 2, 3) and (2, 3, 4). .. versionadded:: 0.8 .. versionchanged:: 0.9.0 Now a subclass of :class:`TopologyObject`; now uses :attr:`__slots__` and stores atoms in :attr:`atoms` attribute. .. versionchanged:: 0.11.0 Renamed to Dihedral (was Torsion) """ # http://cbio.bmt.tue.nl/pumma/uploads/Theory/dihedral.png btype = 'dihedral'
[docs] def dihedral(self, pbc=True): """Calculate the dihedral angle in degrees. Dihedral angle around axis connecting atoms 1 and 2 (i.e. the angle between the planes spanned by atoms (0,1,2) and (1,2,3)):: 3 | 1-----2 / 0 Note ---- The numerical precision is typically not better than 4 decimals (and is only tested to 3 decimals). .. versionadded:: 0.9.0 .. versionchanged:: 0.19.0 Added pbc keyword, default True """ box = self.universe.dimensions if pbc else None A, B, C, D = self.atoms return np.rad2deg(distances.calc_dihedrals( A.position, B.position, C.position, D.position, box))
value = dihedral
# subclass Dihedral to inherit dihedral method
[docs]class ImproperDihedral(Dihedral): """ Improper Dihedral (improper dihedral angle) between four :class:`~MDAnalysis.core.groups.Atom` instances. MDAnalysis treats the improper dihedral angle as the angle between the planes formed by Atoms (1, 2, 3) and (2, 3, 4). .. warning:: Definitions of Atom ordering in improper dihedrals can change. Check the definitions here against your software. .. versionadded:: 0.9.0 .. versionchanged:: 0.11.0 Renamed to ImproperDihedral (was Improper_Torsion) """ # http://cbio.bmt.tue.nl/pumma/uploads/Theory/improper.png btype = 'improper'
[docs] def improper(self): """Improper dihedral angle in degrees. Note ---- The numerical precision is typically not better than 4 decimals (and is only tested to 3 decimals). """ return self.dihedral()
[docs]class UreyBradley(TopologyObject): """A Urey-Bradley angle between two :class:`~MDAnalysis.core.groups.Atom` instances. Two :class:`UreyBradley` instances can be compared with the ``==`` and ``!=`` operators. A UreyBradley angle is equal to another if the same atom numbers are involved. .. versionadded:: 1.0.0 """ btype = 'ureybradley'
[docs] def partner(self, atom): """UreyBradley.partner(Atom) Returns ------- the other :class:`~MDAnalysis.core.groups.Atom` in this interaction """ if atom == self.atoms[0]: return self.atoms[1] elif atom == self.atoms[1]: return self.atoms[0] else: raise ValueError("Unrecognised Atom")
[docs] def distance(self, pbc=True): """Distance between the atoms. """ box = self.universe.dimensions if pbc else None return distances.calc_bonds(self[0].position, self[1].position, box)
value = distance
[docs]class CMap(TopologyObject): """ Coupled-torsion correction map term between five :class:`~MDAnalysis.core.groups.Atom` instances. .. versionadded:: 1.0.0 """ btype = 'cmap'
[docs]class TopologyDict(object): """A customised dictionary designed for sorting the bonds, angles and dihedrals present in a group of atoms. Usage:: topologydict = TopologyDict(members) TopologyDicts are also built lazily from a :class:`TopologyGroup.topDict` attribute. The :class:`TopologyDict` collects all the selected topology type from the atoms and categorises them according to the types of the atoms within. A :class:`TopologyGroup` containing all of a given bond type can be made by querying with the appropriate key. The keys to the :class:`TopologyDict` are a tuple of the atom types that the bond represents and can be viewed using the :meth:`keys` method. For example, from a system containing pure ethanol :: >>> td = u.bonds.topDict >>> td.keys() [('C', 'C'), ('C', 'H'), ('O', 'H'), ('C', 'O')] >>> td['C', 'O'] < TopologyGroup containing 912 bonds > .. Note:: The key for a bond is taken from the type attribute of the atoms. Getting and setting types of bonds is done smartly, so a C-C-H angle is considered identical to a H-C-C angle. Duplicate entries are automatically removed upon creation and combination of different Dicts. This means a bond between atoms 1 and 2 will only ever appear once in a dict despite both atoms 1 and 2 having the bond in their :attr:`bond` attribute. Two :class:`TopologyDict` instances can be combined using addition and it will not create any duplicate bonds in the process. Arguments --------- members : A list of :class:`TopologyObject` instances .. versionadded:: 0.8 .. versionchanged:: 0.9.0 Changed initialisation to use a list of :class:`TopologyObject` instances instead of list of atoms; now used from within :class:`TopologyGroup` instead of accessed from :class:`AtomGroup`. """ def __init__(self, topologygroup): if not isinstance(topologygroup, TopologyGroup): raise TypeError("Can only construct from TopologyGroup") self.dict = dict() self._u = topologygroup.universe self.toptype = topologygroup.btype for b in topologygroup: btype = b.type try: self.dict[btype] += [b] except KeyError: self.dict[btype] = [b] self._removeDupes() def _removeDupes(self): """Sorts through contents and makes sure that there are no duplicate keys (through type reversal) """ newdict = dict() # Go through all keys, if the reverse of the key exists add this to # that entry else make a new entry for k in self.dict: if not k[::-1] in newdict: newdict[k] = self.dict[k] else: newdict[k[::-1]] += self.dict[k] self.dict = newdict @property def universe(self): return self._u def __len__(self): """Returns the number of types of bond in the topology dictionary""" return len(self.dict.keys())
[docs] def keys(self): """Returns a list of the different types of available bonds""" return self.dict.keys()
def __iter__(self): """Iterator over keys in this dictionary""" return iter(self.dict) def __repr__(self): return "<TopologyDict with {num} unique {type}s>".format( num=len(self), type=self.toptype) def __getitem__(self, key): """Returns a TopologyGroup matching the criteria if possible, otherwise returns ``None`` """ if key in self: if key in self.dict: selection = self.dict[key] else: selection = self.dict[key[::-1]] bix = np.vstack([s.indices for s in selection]) return TopologyGroup(bix, self._u, btype=self.toptype) else: raise KeyError(key) def __contains__(self, other): """ Returns boolean on whether a given type exists within this dictionary For topology groups the key (1,2,3) is considered the same as (3,2,1) """ return other in self.dict or other[::-1] in self.dict
_BTYPE_TO_SHAPE = {'bond': 2, 'ureybradley': 2, 'angle': 3, 'dihedral': 4, 'improper': 4, 'cmap': 5}
[docs]class TopologyGroup(object): """A container for a groups of bonds. All bonds of a certain types can be retrieved from within the :class:`TopologyGroup` by querying with a tuple of types:: tg2 = tg.select_bonds([key]) Where *key* describes the desired bond as a tuple of the involved :class:`~MDAnalysis.core.groups.Atom` types, as defined by the .type Atom attribute). A list of available keys can be displayed using the :meth:`types` method. Alternatively, all the bonds which are in a given :class:`~MDAnalysis.core.groups.AtomGroup` can be extracted using :meth:`atomgroup_intersection`:: tg2 = tg.atomgroup_intersection(ag) This allows the keyword *strict* to be given, which forces all members of all bonds to be inside the AtomGroup passed to it. Finally, a TopologyGroup can be sliced similarly to AtomGroups:: tg2 = tg[5:10] The :meth:`bonds`, :meth:`angles` and :meth:`dihedrals` methods offer a "shortcut" to the Cython distance calculation functions in :class:`MDAnalysis.lib.distances`. TopologyGroups can be combined with TopologyGroups of the same bond type (ie can combine two angle containing TopologyGroups). .. versionadded:: 0.8 .. versionchanged:: 0.9.0 Overhauled completely: (1) Added internal :class:`TopologyDict` accessible by the :attr:`topDict` attribute. (2) :meth:`selectBonds` allows the :attr:`topDict` to be queried with tuple of types. (3) Added :meth:`atomgroup_intersection` to allow bonds which are in a given :class:`AtomGroup` to be retrieved. .. versionchanged:: 0.10.0 Added :func:`from_indices` constructor, allowing class to be created from indices. Can now create empty Group. Renamed :meth:`dump_contents` to :meth:`to_indices` .. versionchanged:: 0.11.0 Added `values` method to return the size of each object in this group Deprecated selectBonds method in favour of select_bonds .. versionchanged:: 0.19.0 Empty TopologyGroup now returns correctly shaped empty array via indices property and to_indices() .. versionchanged::1.0.0 ``type``, ``guessed``, and ``order`` are no longer reshaped to arrays with an extra dimension """ def __init__(self, bondidx, universe, btype=None, type=None, guessed=None, order=None): if btype is None: # guess what I am # difference between dihedral and improper # not really important self.btype = {2: 'bond', 3: 'angle', 4: 'dihedral'}[len(bondidx[0])] elif btype in _BTYPE_TO_SHAPE: self.btype = btype else: raise ValueError("Unsupported btype, use one of '{}'" "".format(', '.join(_BTYPE_TO_SHAPE))) bondidx = np.asarray(bondidx) nbonds = len(bondidx) # remove duplicate bonds if type is None: type = np.repeat(None, nbonds) if guessed is None: guessed = np.repeat(True, nbonds) elif guessed is True or guessed is False: guessed = np.repeat(guessed, nbonds) else: guessed = np.asarray(guessed, dtype=np.bool) if order is None: order = np.repeat(None, nbonds) if nbonds > 0: uniq, uniq_idx = util.unique_rows(bondidx, return_index=True) self._bix = uniq self._bondtypes = type[uniq_idx] self._guessed = guessed[uniq_idx] self._order = order[uniq_idx] # Create vertical AtomGroups self._ags = [universe.atoms[self._bix[:, i]] for i in range(self._bix.shape[1])] else: # Empty TopologyGroup self._bix = np.array([]) self._bondtypes = np.array([]) self._guessed = np.array([]) self._order = np.array([]) self._ags = [] self._u = universe self._cache = dict() # used for topdict saving @property def universe(self): """The Universe that we belong to""" return self._u
[docs] def select_bonds(self, selection): """Retrieves a selection from this topology group based on types. .. seeAlso :meth:`types` .. versionadded 0.9.0 """ return self.topDict[selection]
selectBonds = select_bonds
[docs] def types(self): """Return a list of the bond types in this TopologyGroup .. versionadded 0.9.0 """ return list(self.topDict.keys())
@property @cached('dict') def topDict(self): """ Returns the TopologyDict for this topology group. This is used for the select_bonds method when fetching a certain type of bond. This is a cached property so will be generated the first time it is accessed. .. versionadded 0.9.0 """ return TopologyDict(self)
[docs] def atomgroup_intersection(self, ag, **kwargs): """Retrieve all bonds from within this TopologyGroup that are within the AtomGroup which is passed. Parameters ---------- ag : AtomGroup The `:class:~MDAnalysis.core.groups.AtomGroup` to intersect with. strict : bool Only retrieve bonds which are completely contained within the AtomGroup. [``False``] .. versionadded:: 0.9.0 """ # Issue #780 - if self is empty, return self to avoid invalid mask if not self: return self # Strict requires all items in a row to be seen, # otherwise any item in a row func = np.all if kwargs.get('strict', False) else np.any atom_idx = ag.indices # Create a list of boolean arrays, # each representing a column of bond indices. seen = [np.in1d(col, atom_idx) for col in self._bix.T] # Create final boolean mask by summing across rows mask = func(seen, axis=0) return self[mask]
@property def indices(self): """all bond indices See Also -------- to_indices : function that just returns `indices` """ if not self: # empty TG shape = _BTYPE_TO_SHAPE[self.btype] return np.zeros((0, shape), dtype=np.int32) else: return self._bix
[docs] def to_indices(self): """Return a data structure with atom indices describing the bonds. This format should be identical to the original contents of the entries in universe._topology. Note that because bonds are sorted as they are initialised, the order that atoms are defined in each entry might be reversed. Returns ------- indices : tuple A tuple of tuples which define the contents of this TopologyGroup in terms of the atom numbers. (0 based index within u.atoms) .. versionadded:: 0.9.0 .. versionchanged:: 0.10.0 Renamed from "dump_contents" to "to_indices" """ return self.indices
dump_contents = to_indices def __len__(self): """Number of bonds in the topology group""" return self._bix.shape[0] def __add__(self, other): """Combine two TopologyGroups together. Can combined two TopologyGroup of the same type, or add a single TopologyObject to a TopologyGroup. """ # check addition is sane if not isinstance(other, (TopologyObject, TopologyGroup)): raise TypeError("Can only combine TopologyObject or " "TopologyGroup to TopologyGroup, not {0}" "".format(type(other))) # cases where either other or self is empty TG if not other: # adding empty TG to me return self if not self: if isinstance(other, TopologyObject): # Reshape indices to be 2d array return TopologyGroup(other.indices[None, :], other.universe, btype=other.btype, type=np.array([other._bondtype]), guessed=np.array([other.is_guessed]), order=np.array([other.order]), ) else: return TopologyGroup(other.indices, other.universe, btype=other.btype, type=other._bondtypes, guessed=other._guessed, order=other._order, ) else: if not other.btype == self.btype: raise TypeError("Cannot add different types of " "TopologyObjects together") if isinstance(other, TopologyObject): # add TO to me return TopologyGroup( np.concatenate([self.indices, other.indices[None, :]]), self.universe, btype=self.btype, type=np.concatenate([self._bondtypes, np.array([other._bondtype])]), guessed=np.concatenate([self._guessed, np.array([other.is_guessed])]), order=np.concatenate([self._order, np.array([other.order])]), ) else: # add TG to me return TopologyGroup( np.concatenate([self.indices, other.indices]), self.universe, btype=self.btype, type=np.concatenate([self._bondtypes, other._bondtypes]), guessed=np.concatenate([self._guessed, other._guessed]), order=np.concatenate([self._order, other._order]), ) def __getitem__(self, item): """Returns a particular bond as single object or a subset of this TopologyGroup as another TopologyGroup .. versionchanged:: 0.10.0 Allows indexing via boolean numpy array """ # Grab a single Item, similar to Atom/AtomGroup relationship if isinstance(item, numbers.Integral): outclass = {'bond': Bond, 'angle': Angle, 'dihedral': Dihedral, 'improper': ImproperDihedral, 'ureybradley': UreyBradley, 'cmap': CMap}[self.btype] return outclass(self._bix[item], self._u, type=self._bondtypes[item], guessed=self._guessed[item], order=self._order[item]) else: # Slice my index array with the item return self.__class__(self._bix[item], self._u, btype=self.btype, type=self._bondtypes[item], guessed=self._guessed[item], order=self._order[item],) def __contains__(self, item): """Tests if this TopologyGroup contains a bond""" return item.indices in self._bix def __repr__(self): return "<TopologyGroup containing {num} {type}s>".format( num=len(self), type=self.btype) def __eq__(self, other): """Test if contents of TopologyGroups are equal""" return np.array_equal(self.indices, other.indices) def __ne__(self, other): return not self == other def __nonzero__(self): return not len(self) == 0 @property def atom1(self): """The first atom in each TopologyObject in this Group""" return self._ags[0] @property def atom2(self): """The second atom in each TopologyObject in this Group""" return self._ags[1] @property def atom3(self): """The third atom in each TopologyObject in this Group""" try: return self._ags[2] except IndexError: nvert = _BTYPE_TO_SHAPE[self.btype] raise_from( IndexError( "TopologyGroup of {}s only has {} vertical AtomGroups".format( self.btype, nvert)), None) @property def atom4(self): """The fourth atom in each TopologyObject in this Group""" try: return self._ags[3] except IndexError: nvert = _BTYPE_TO_SHAPE[self.btype] raise_from( IndexError( "TopologyGroup of {}s only has {} vertical AtomGroups".format( self.btype, nvert)), None) # Distance calculation methods below # "Slow" versions exist as a way of testing the Cython implementations
[docs] def values(self, **kwargs): """Return the size of each object in this Group :Keywords: *pbc* apply periodic boundary conditions when calculating distance [``False``] *result* allows a predefined results array to be used, note that this will be overwritten .. versionadded:: 0.11.0 """ if self.btype == 'bond': return self.bonds(**kwargs) elif self.btype == 'angle': return self.angles(**kwargs) elif self.btype == 'dihedral': return self.dihedrals(**kwargs) elif self.btype == 'improper': return self.dihedrals(**kwargs)
def _bondsSlow(self, pbc=False): # pragma: no cover """Slow version of bond (numpy implementation)""" if not self.btype == 'bond': return TypeError("TopologyGroup is not of type 'bond'") else: bond_dist = self._ags[0].positions - self._ags[1].positions if pbc: box = self._ags[0].dimensions # orthogonal and divide by zero check if (box[6:9] == 90.).all() and not (box[0:3] == 0).any(): bond_dist -= np.rint(bond_dist / box[0:3]) * box[0:3] else: raise ValueError("Only orthogonal boxes supported") return np.array([mdamath.norm(a) for a in bond_dist])
[docs] def bonds(self, pbc=False, result=None): """Calculates the distance between all bonds in this TopologyGroup :Keywords: *pbc* apply periodic boundary conditions when calculating distance [False] *result* allows a predefined results array to be used, note that this will be overwritten Uses cython implementation """ if not self.btype == 'bond': raise TypeError("TopologyGroup is not of type 'bond'") if not result: result = np.zeros(len(self), np.float64) if pbc: return distances.calc_bonds(self._ags[0].positions, self._ags[1].positions, box=self._ags[0].dimensions, result=result) else: return distances.calc_bonds(self._ags[0].positions, self._ags[1].positions, result=result)
def _anglesSlow(self): # pragma: no cover """Slow version of angle (numpy implementation)""" if not self.btype == 'angle': raise TypeError("TopologyGroup is not of type 'angle'") vec1 = self._ags[0].positions - self._ags[1].positions vec2 = self._ags[2].positions - self._ags[1].positions angles = np.array([mdamath.angle(a, b) for a, b in zip(vec1, vec2)]) return angles
[docs] def angles(self, result=None, pbc=False): """Calculates the angle in radians formed between a bond between atoms 1 and 2 and a bond between atoms 2 & 3 Parameters ---------- result : array_like allows a predefined results array to be used, note that this will be overwritten pbc : bool apply periodic boundary conditions when calculating angles [``False``] this is important when connecting vectors between atoms might require minimum image convention Returns ------- angles : ndarray .. versionchanged :: 0.9.0 Added *pbc* option (default ``False``) """ if not self.btype == 'angle': raise TypeError("TopologyGroup is not of type 'angle'") if not result: result = np.zeros(len(self), np.float64) if pbc: return distances.calc_angles(self._ags[0].positions, self._ags[1].positions, self._ags[2].positions, box=self._ags[0].dimensions, result=result) else: return distances.calc_angles(self._ags[0].positions, self._ags[1].positions, self._ags[2].positions, result=result)
def _dihedralsSlow(self): # pragma: no cover """Slow version of dihedral (numpy implementation)""" if self.btype not in ['dihedral', 'improper']: raise TypeError("TopologyGroup is not of type 'dihedral' or " "'improper'") ab = self._ags[0].positions - self._ags[1].positions bc = self._ags[1].positions - self._ags[2].positions cd = self._ags[2].positions - self._ags[3].positions return np.array([mdamath.dihedral(a, b, c) for a, b, c in zip(ab, bc, cd)])
[docs] def dihedrals(self, result=None, pbc=False): """Calculate the dihedral angle in radians for this topology group. Defined as the angle between a plane formed by atoms 1, 2 and 3 and a plane formed by atoms 2, 3 and 4. Parameters ---------- result : array_like allows a predefined results array to be used, note that this will be overwritten pbc : bool apply periodic boundary conditions when calculating angles [``False``] this is important when connecting vectors between atoms might require minimum image convention Returns ------- angles : ndarray .. versionchanged:: 0.9.0 Added *pbc* option (default ``False``) """ if self.btype not in ['dihedral', 'improper']: raise TypeError("TopologyGroup is not of type 'dihedral' or " "'improper'") if not result: result = np.zeros(len(self), np.float64) if pbc: return distances.calc_dihedrals(self._ags[0].positions, self._ags[1].positions, self._ags[2].positions, self._ags[3].positions, box=self._ags[0].dimensions, result=result) else: return distances.calc_dihedrals(self._ags[0].positions, self._ags[1].positions, self._ags[2].positions, self._ags[3].positions, result=result)