1# -*- Mode: python; tab-width: 4; indent-tabs-mode:nil; coding:utf-8 -*- 2# vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4 fileencoding=utf-8 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========================================================== 25Core objects: Containers --- :mod:`MDAnalysis.core.groups` 26========================================================== 27 28The :class:`~MDAnalysis.core.universe.Universe` instance contains all the 29particles in the system (which MDAnalysis calls :class:`Atom`). Groups of 30:class:`atoms<Atom>` are handled as :class:`AtomGroup` instances. The 31:class:`AtomGroup` is probably the most important object in MDAnalysis because 32virtually everything can be accessed through it. :class:`AtomGroup` instances 33can be easily created (e.g., from an :meth:`AtomGroup.select_atoms` selection or 34simply by slicing). 35 36For convenience, chemically meaningful groups of :class:`Atoms<Atom>` such as a 37:class:`Residue` or a :class:`Segment` (typically a whole molecule or all of the 38solvent) also exist as containers, as well as groups of these units 39(:class:`ResidueGroup`, :class:`SegmentGroup`). 40 41 42Classes 43======= 44 45Collections 46----------- 47 48.. autoclass:: AtomGroup 49 :members: 50 :inherited-members: 51.. autoclass:: ResidueGroup 52 :members: 53 :inherited-members: 54.. autoclass:: SegmentGroup 55 :members: 56 :inherited-members: 57.. autoclass:: UpdatingAtomGroup 58 :members: 59 60Chemical units 61-------------- 62 63.. autoclass:: Atom 64 :members: 65 :inherited-members: 66.. autoclass:: Residue 67 :members: 68 :inherited-members: 69.. autoclass:: Segment 70 :members: 71 :inherited-members: 72 73Levels 74------ 75 76Each of the above classes has a *level* attribute. This can be used to verify 77that two objects are of the same level, or to access a particular class:: 78 79 u = mda.Universe() 80 81 ag = u.atoms[:10] 82 at = u.atoms[11] 83 84 ag.level == at.level # Returns True 85 86 ag.level.singular # Returns Atom class 87 at.level.plural # Returns AtomGroup class 88 89""" 90from __future__ import absolute_import 91from six.moves import zip 92from six import string_types 93 94from collections import namedtuple 95import numpy as np 96import functools 97import itertools 98import numbers 99import os 100import warnings 101 102from numpy.lib.utils import deprecate 103 104from .. import _ANCHOR_UNIVERSES 105from ..lib import util 106from ..lib.util import cached, warn_if_not_unique, unique_int_1d 107from ..lib import distances 108from ..lib import transformations 109from ..selections import get_writer as get_selection_writer_for 110from . import selection 111from . import flags 112from ..exceptions import NoDataError 113from . import topologyobjects 114from ._get_readers import get_writer_for 115 116 117def _unpickle(uhash, ix): 118 try: 119 u = _ANCHOR_UNIVERSES[uhash] 120 except KeyError: 121 # doesn't provide as nice an error message as before as only hash of universe is stored 122 # maybe if we pickled the filename too we could do better... 123 raise RuntimeError( 124 "Couldn't find a suitable Universe to unpickle AtomGroup onto " 125 "with Universe hash '{}'. Available hashes: {}" 126 "".format(uhash, ', '.join([str(k) 127 for k in _ANCHOR_UNIVERSES.keys()]))) 128 return u.atoms[ix] 129 130def _unpickle_uag(basepickle, selections, selstrs): 131 bfunc, bargs = basepickle[0], basepickle[1:][0] 132 basegroup = bfunc(*bargs) 133 return UpdatingAtomGroup(basegroup, selections, selstrs) 134 135 136def make_classes(): 137 """Make a fresh copy of all classes 138 139 Returns 140 ------- 141 Two dictionaries. One with a set of :class:`_TopologyAttrContainer` classes 142 to serve as bases for :class:`~MDAnalysis.core.universe.Universe`\ -specific 143 MDA container classes. Another with the final merged versions of those 144 classes. The classes themselves are used as hashing keys. 145 146 """ 147 bases = {} 148 classes = {} 149 groups = (AtomGroup, ResidueGroup, SegmentGroup) 150 components = (Atom, Residue, Segment) 151 152 # The 'GBase' middle man is needed so that a single topologyattr 153 # patching applies automatically to all groups. 154 GBase = bases[GroupBase] = _TopologyAttrContainer._subclass(is_group=True) 155 for cls in groups: 156 bases[cls] = GBase._subclass(is_group=True) 157 # CBase for patching all components 158 CBase = bases[ComponentBase] = _TopologyAttrContainer._subclass( 159 is_group=False) 160 for cls in components: 161 bases[cls] = CBase._subclass(is_group=False) 162 163 # Initializes the class cache. 164 for cls in groups + components: 165 classes[cls] = bases[cls]._mix(cls) 166 167 return bases, classes 168 169 170class _TopologyAttrContainer(object): 171 """Class factory for receiving sets of :class:`TopologyAttr` objects. 172 173 :class:`_TopologyAttrContainer` is a convenience class to encapsulate the 174 functions that deal with: 175 * the import and namespace transplant of 176 :class:`~MDAnalysis.core.topologyattrs.TopologyAttr` objects; 177 * the copying (subclassing) of itself to create distinct bases for the 178 different container classes (:class:`AtomGroup`, :class:`ResidueGroup`, 179 :class:`SegmentGroup`, :class:`Atom`, :class:`Residue`, :class:`Segment`, 180 and subclasses thereof); 181 * the mixing (subclassing and co-inheritance) with the container classes. 182 The mixed subclasses become the final container classes specific to each 183 :class:`~MDAnalysis.core.universe.Universe`. 184 """ 185 @classmethod 186 def _subclass(cls, is_group): 187 """Factory method returning :class:`_TopologyAttrContainer` subclasses. 188 189 Parameters 190 ---------- 191 is_group : bool 192 The :attr:`_is_group` of the returned class will be set to 193 `is_group`. This is used to distinguish between Groups 194 (:class:`AtomGroup` etc.) and Components (:class:`Atom` etc.) in 195 internal methods when considering actions such as addition of 196 objects or adding 197 :class:`TopologyAttributes<MDAnalysis.core.topologyattrs.TopologyAttr>` 198 to them. 199 200 Returns 201 ------- 202 type 203 A subclass of :class:`_TopologyAttrContainer`, with the same name. 204 """ 205 newcls = type(cls.__name__, (cls,), {'_is_group': bool(is_group)}) 206 if is_group: 207 newcls._SETATTR_WHITELIST = { 208 'positions', 'velocities', 'forces', 'dimensions', 209 'atoms', 'residue', 'residues', 'segment', 'segments', 210 } 211 else: 212 newcls._SETATTR_WHITELIST = { 213 'position', 'velocity', 'force', 'dimensions', 214 'atoms', 'residue', 'residues', 'segment', 215 } 216 217 return newcls 218 219 @classmethod 220 def _mix(cls, other): 221 """Creates a subclass with ourselves and another class as parents. 222 223 Classes mixed at this point override :meth:`__new__`, causing further 224 instantiations to shortcut to :meth:`~object.__new__` (skipping the 225 cache-fetch process for :class:`_MutableBase` subclasses). 226 227 The new class will have an attribute `_derived_class` added, pointing 228 to itself. This pointer instructs which class to use when 229 slicing/adding instances of the new class. At initialization time, the 230 new class may choose to point `_derived_class` to another class (as is 231 done in the initialization of :class:`UpdatingAtomGroup`). 232 233 Parameters 234 ---------- 235 other : type 236 The class to mix with ourselves. 237 238 Returns 239 ------- 240 type 241 A class of parents :class:`_ImmutableBase`, *other* and this class. 242 Its name is the same as *other*'s. 243 """ 244 newcls = type(other.__name__, (_ImmutableBase, other, cls), {}) 245 newcls._derived_class = newcls 246 return newcls 247 248 @classmethod 249 def _add_prop(cls, attr): 250 """Add `attr` into the namespace for this class 251 252 Parameters 253 ---------- 254 attr : A :class:`TopologyAttr` object 255 """ 256 def getter(self): 257 return attr.__getitem__(self) 258 259 def setter(self, values): 260 return attr.__setitem__(self, values) 261 262 if cls._is_group: 263 setattr(cls, attr.attrname, 264 property(getter, setter, None, attr.groupdoc)) 265 cls._SETATTR_WHITELIST.add(attr.attrname) 266 else: 267 setattr(cls, attr.singular, 268 property(getter, setter, None, attr.singledoc)) 269 cls._SETATTR_WHITELIST.add(attr.singular) 270 271 def __setattr__(self, attr, value): 272 # `ag.this = 42` calls setattr(ag, 'this', 42) 273 if not (attr.startswith('_') or # 'private' allowed 274 attr in self._SETATTR_WHITELIST or # known attributes allowed 275 hasattr(self, attr)): # preexisting (eg properties) allowed 276 raise AttributeError( 277 "Cannot set arbitrary attributes to a {}".format( 278 'Group' if self._is_group else 'Component')) 279 # if it is, we allow the setattr to proceed by deferring to the super 280 # behaviour (ie do it) 281 super(_TopologyAttrContainer, self).__setattr__(attr, value) 282 283 284class _MutableBase(object): 285 """ 286 Base class that merges appropriate :class:`_TopologyAttrContainer` classes. 287 288 Implements :meth:`__new__`. In it the instantiating class is fetched from 289 :attr:`~MDAnalysis.core.universe.Universe._classes`. If there is a cache 290 miss, a merged class is made 291 with a base from :attr:`~MDAnalysis.core.universe.Universe._class_bases` 292 and cached. 293 294 The classes themselves are used as the cache dictionary keys for simplcity 295 in cache retrieval. 296 297 """ 298 def __new__(cls, *args, **kwargs): 299 # This pre-initialization wrapper must be pretty generic to 300 # allow for different initialization schemes of the possible classes. 301 # All we really need here is to fish a universe out of the arg list. 302 # The AtomGroup cases get priority and are fished out first. 303 try: 304 u = args[-1].universe 305 except (IndexError, AttributeError): 306 try: 307 # older AtomGroup init method.. 308 u = args[0][0].universe 309 except (TypeError, IndexError, AttributeError): 310 from .universe import Universe 311 # Let's be generic and get the first argument that's either a 312 # Universe, a Group, or a Component, and go from there. 313 # This is where the UpdatingAtomGroup args get matched. 314 for arg in args+tuple(kwargs.values()): 315 if isinstance(arg, (Universe, GroupBase, 316 ComponentBase)): 317 u = arg.universe 318 break 319 else: 320 raise TypeError("No universe, or universe-containing " 321 "object passed to the initialization of " 322 "{}".format(cls.__name__)) 323 try: 324 return object.__new__(u._classes[cls]) 325 except KeyError: 326 # Cache miss. Let's find which kind of class this is and merge. 327 try: 328 parent_cls = next(u._class_bases[parent] 329 for parent in cls.mro() 330 if parent in u._class_bases) 331 except StopIteration: 332 raise TypeError("Attempted to instantiate class '{}' but " 333 "none of its parents are known to the " 334 "universe. Currently possible parent " 335 "classes are: {}".format(cls.__name__, 336 str(sorted(u._class_bases.keys())))) 337 newcls = u._classes[cls] = parent_cls._mix(cls) 338 return object.__new__(newcls) 339 340 341class _ImmutableBase(object): 342 """Class used to shortcut :meth:`__new__` to :meth:`object.__new__`. 343 344 """ 345 # When mixed via _TopologyAttrContainer._mix this class has MRO priority. 346 # Setting __new__ like this will avoid having to go through the 347 # cache lookup if the class is reused (as in ag._derived_class(...)). 348 __new__ = object.__new__ 349 350 351 352def _only_same_level(function): 353 @functools.wraps(function) 354 def wrapped(self, other): 355 if not isinstance(other, (ComponentBase, GroupBase)): # sanity check 356 raise TypeError("Can't perform '{}' between objects:" 357 " '{}' and '{}'".format( 358 function.__name__, 359 type(self).__name__, 360 type(other).__name__)) 361 if self.level != other.level: 362 raise TypeError("Can't perform '{}' on different level objects" 363 "".format(function.__name__)) 364 if self.universe is not other.universe: 365 raise ValueError( 366 "Can't operate on objects from different Universes") 367 return function(self, other) 368 return wrapped 369 370 371class GroupBase(_MutableBase): 372 """Base class from which a :class:`<~MDAnalysis.core.universe.Universe`\ 's 373 Group class is built. 374 375 Instances of :class:`GroupBase` provide the following operations that 376 conserve element repetitions and order: 377 378 +-------------------------------+------------+----------------------------+ 379 | Operation | Equivalent | Result | 380 +===============================+============+============================+ 381 | ``len(s)`` | | number of elements (atoms, | 382 | | | residues or segment) in | 383 | | | the group | 384 +-------------------------------+------------+----------------------------+ 385 | ``s == t`` | | test if ``s`` and ``t`` | 386 | | | contain the same elements | 387 | | | in the same order | 388 +-------------------------------+------------+----------------------------+ 389 | ``x in s`` | | test if component ``x`` is | 390 | | | part of group ``s`` | 391 +-------------------------------+------------+----------------------------+ 392 | ``s.concatenate(t)`` | ``s + t`` | new Group with elements | 393 | | | from ``s`` and from ``t`` | 394 +-------------------------------+------------+----------------------------+ 395 | ``s.subtract(t)`` | | new Group with elements | 396 | | | from ``s`` that are not | 397 | | | in ``t`` | 398 +-------------------------------+------------+----------------------------+ 399 400 The following operations treat the Group as set. Any result will have any 401 duplicate entries removed and the Group will be sorted. 402 403 +-------------------------------+------------+----------------------------+ 404 | Operation | Equivalent | Result | 405 +===============================+============+============================+ 406 | ``s.isdisjoint(t)`` | | ``True`` if ``s`` and | 407 | | | ``t`` do not share | 408 | | | elements | 409 +-------------------------------+------------+----------------------------+ 410 | ``s.issubset(t)`` | | test if all elements of | 411 | | | ``s`` are part of ``t`` | 412 +-------------------------------+------------+----------------------------+ 413 | ``s.is_strict_subset(t)`` | | test if all elements of | 414 | | | ``s`` are part of ``t``, | 415 | | | and ``s != t`` | 416 +-------------------------------+------------+----------------------------+ 417 | ``s.issuperset(t)`` | | test if all elements of | 418 | | | ``t`` are part of ``s`` | 419 +-------------------------------+------------+----------------------------+ 420 | ``s.is_strict_superset(t)`` | | test if all elements of | 421 | | | ``t`` are part of ``s``, | 422 | | | and ``s != t`` | 423 +-------------------------------+------------+----------------------------+ 424 | ``s.union(t)`` | ``s | t`` | new Group with elements | 425 | | | from both ``s`` and ``t`` | 426 +-------------------------------+------------+----------------------------+ 427 | ``s.intersection(t)`` | ``s & t`` | new Group with elements | 428 | | | common to ``s`` and ``t`` | 429 +-------------------------------+------------+----------------------------+ 430 | ``s.difference(t)`` | ``s - t`` | new Group with elements of | 431 | | | ``s`` that are not in ``t``| 432 +-------------------------------+------------+----------------------------+ 433 | ``s.symmetric_difference(t)`` | ``s ^ t`` | new Group with elements | 434 | | | that are part of ``s`` or | 435 | | | ``t`` but not both | 436 +-------------------------------+------------+----------------------------+ 437 """ 438 def __init__(self, *args): 439 try: 440 if len(args) == 1: 441 # list of atoms/res/segs, old init method 442 ix = [at.ix for at in args[0]] 443 u = args[0][0].universe 444 else: 445 # current/new init method 446 ix, u = args 447 except (AttributeError, # couldn't find ix/universe 448 TypeError): # couldn't iterate the object we got 449 raise TypeError( 450 "Can only initialise a Group from an iterable of Atom/Residue/" 451 "Segment objects eg: AtomGroup([Atom1, Atom2, Atom3]) " 452 "or an iterable of indices and a Universe reference " 453 "eg: AtomGroup([0, 5, 7, 8], u).") 454 455 # indices for the objects I hold 456 self._ix = np.asarray(ix, dtype=np.intp) 457 self._u = u 458 self._cache = dict() 459 460 def __hash__(self): 461 return hash((self._u, self.__class__, tuple(self.ix.tolist()))) 462 463 def __len__(self): 464 return len(self.ix) 465 466 def __getitem__(self, item): 467 # supports 468 # - integer access 469 # - boolean slicing 470 # - fancy indexing 471 # because our _ix attribute is a numpy array 472 # it can be sliced by all of these already, 473 # so just return ourselves sliced by the item 474 if isinstance(item, numbers.Integral): 475 return self.level.singular(self.ix[item], self.universe) 476 else: 477 if isinstance(item, list) and item: # check for empty list 478 # hack to make lists into numpy arrays 479 # important for boolean slicing 480 item = np.array(item) 481 # We specify _derived_class instead of self.__class__ to allow 482 # subclasses, such as UpdatingAtomGroup, to control the class 483 # resulting from slicing. 484 return self._derived_class(self.ix[item], self.universe) 485 486 def __repr__(self): 487 name = self.level.name 488 return ("<{}Group with {} {}{}>" 489 "".format(name.capitalize(), len(self), name, 490 "s"[len(self)==1:])) # Shorthand for a conditional plural 's'. 491 492 def __str__(self): 493 name = self.level.name 494 if len(self) <= 10: 495 return '<{}Group {}>'.format(name.capitalize(), repr(list(self))) 496 else: 497 return '<{}Group {}, ..., {}>'.format(name.capitalize(), 498 repr(list(self)[:3])[:-1], 499 repr(list(self)[-3:])[1:]) 500 501 def __add__(self, other): 502 """Concatenate the Group with another Group or Component of the same 503 level. 504 505 Parameters 506 ---------- 507 other : Group or Component 508 Group or Component with `other.level` same as `self.level` 509 510 Returns 511 ------- 512 Group 513 Group with elements of `self` and `other` concatenated 514 515 """ 516 return self.concatenate(other) 517 518 def __radd__(self, other): 519 """Using built-in sum requires supporting 0 + self. If other is 520 anything other 0, an exception will be raised. 521 522 Parameters 523 ---------- 524 other : int 525 Other should be 0, or else an exception will be raised. 526 527 Returns 528 ------- 529 self 530 Group with elements of `self` reproduced 531 532 """ 533 if other == 0: 534 return self._derived_class(self.ix, self.universe) 535 else: 536 raise TypeError("unsupported operand type(s) for +:" 537 " '{}' and '{}'".format(type(self).__name__, 538 type(other).__name__)) 539 def __sub__(self, other): 540 return self.difference(other) 541 542 @_only_same_level 543 def __eq__(self, other): 544 """Test group equality. 545 546 Two groups are equal if they contain the same indices in 547 the same order. Groups that are not at the same level or that belong 548 to different :class:`Universes<MDAnalysis.core.universe.Universe>` 549 cannot be compared. 550 """ 551 o_ix = other.ix 552 return np.array_equal(self.ix, o_ix) 553 554 def __contains__(self, other): 555 if not other.level == self.level: 556 # maybe raise TypeError instead? 557 # eq method raises Error for wrong comparisons 558 return False 559 return other.ix in self.ix 560 561 def __or__(self, other): 562 return self.union(other) 563 564 def __and__(self, other): 565 return self.intersection(other) 566 567 def __xor__(self, other): 568 return self.symmetric_difference(other) 569 570 @property 571 def universe(self): 572 """The underlying :class:`~MDAnalysis.core.universe.Universe` the group 573 belongs to. 574 """ 575 return self._u 576 577 @property 578 def ix(self): 579 """Unique indices of the components in the Group. 580 581 - If this Group is an :class:`AtomGroup`, these are the 582 indices of the :class:`Atom` instances. 583 - If it is a :class:`ResidueGroup`, these are the indices of 584 the :class:`Residue` instances. 585 - If it is a :class:`SegmentGroup`, these are the indices of 586 the :class:`Segment` instances. 587 588 """ 589 return self._ix 590 591 @property 592 def ix_array(self): 593 """Unique indices of the components in the Group. 594 595 For a Group, :attr:`ix_array` is the same as :attr:`ix`. This method 596 gives a consistent API between components and groups. 597 598 See Also 599 -------- 600 :attr:`ix` 601 """ 602 return self._ix 603 604 @property 605 def dimensions(self): 606 """Obtain a copy of the dimensions of the currently loaded Timestep""" 607 return self.universe.trajectory.ts.dimensions.copy() 608 609 @dimensions.setter 610 def dimensions(self, dimensions): 611 self.universe.trajectory.ts.dimensions = dimensions 612 613 @property 614 @cached('isunique') 615 def isunique(self): 616 """Boolean indicating whether all components of the group are unique, 617 i.e., the group contains no duplicates. 618 619 Examples 620 -------- 621 622 >>> ag = u.atoms[[2, 1, 2, 2, 1, 0]] 623 >>> ag 624 <AtomGroup with 6 atoms> 625 >>> ag.isunique 626 False 627 >>> ag2 = ag.unique 628 >>> ag2 629 <AtomGroup with 3 atoms> 630 >>> ag2.isunique 631 True 632 633 634 .. versionadded:: 0.19.0 635 """ 636 if len(self) <= 1: 637 return True 638 # Fast check for uniqueness 639 # 1. get sorted array of component indices: 640 s_ix = np.sort(self._ix) 641 # 2. If the group's components are unique, no pair of adjacent values in 642 # the sorted indices array are equal. We therefore compute a boolean 643 # mask indicating equality of adjacent sorted indices: 644 mask = s_ix[1:] == s_ix[:-1] 645 # 3. The group is unique if all elements in the mask are False. We could 646 # return ``not np.any(mask)`` here but using the following is faster: 647 return not np.count_nonzero(mask) 648 649 @warn_if_not_unique 650 def center(self, weights, pbc=None, compound='group'): 651 """Weighted center of (compounds of) the group 652 653 Computes the weighted center of :class:`Atoms<Atom>` in the group. 654 Weighted centers per :class:`Residue` or per :class:`Segment` can be 655 obtained by setting the `compound` parameter accordingly. 656 657 Parameters 658 ---------- 659 weights : array_like or None 660 Weights to be used. Setting `weights=None` is equivalent to passing 661 identical weights for all atoms of the group. 662 pbc : bool or None, optional 663 If ``True`` and `compound` is ``'group'``, move all atoms to the 664 primary unit cell before calculation. If ``True`` and `compound` is 665 ``'segments'`` or ``'residues'``, the center of each compound will 666 be calculated without moving any :class:`Atoms<Atom>` to keep the 667 compounds intact. Instead, the resulting position vectors will be 668 moved to the primary unit cell after calculation. 669 compound : {'group', 'segments', 'residues'}, optional 670 If ``'group'``, the weighted center of all atoms in the group will 671 be returned as a single position vector. Else, the weighted centers 672 of each :class:`Segment` or :class:`Residue` will be returned as an 673 array of position vectors, i.e. a 2d array. Note that, in any case, 674 *only* the positions of :class:`Atoms<Atom>` *belonging to the 675 group* will be taken into account. 676 677 Returns 678 ------- 679 center : numpy.ndarray 680 Position vector(s) of the weighted center(s) of the group. 681 If `compound` was set to ``'group'``, the output will be a single 682 position vector. 683 If `compound` was set to ``'segments'`` or ``'residues'``, the 684 output will be a 2d array of shape ``(n, 3)`` where ``n`` is the 685 number of compounds. 686 687 Examples 688 -------- 689 690 To find the center of charge of a given :class:`AtomGroup`:: 691 692 >>> sel = u.select_atoms('prop mass > 4.0') 693 >>> sel.center(sel.charges) 694 695 To find the centers of mass per residue of all CA :class:`Atoms<Atom>`:: 696 697 >>> sel = u.select_atoms('name CA') 698 >>> sel.center(sel.masses, compound='residues') 699 700 Notes 701 ----- 702 If the :class:`MDAnalysis.core.flags` flag *use_pbc* is set to 703 ``True`` then the `pbc` keyword is used by default. 704 705 706 .. versionchanged:: 0.19.0 Added `compound` parameter 707 """ 708 709 if pbc is None: 710 pbc = flags['use_pbc'] 711 712 atoms = self.atoms 713 714 # enforce calculations in double precision: 715 dtype = np.float64 716 717 if compound.lower() == 'group': 718 if pbc: 719 coords = atoms.pack_into_box(inplace=False) 720 else: 721 coords = atoms.positions 722 # promote coords or weights to dtype if required: 723 if weights is None: 724 coords = coords.astype(dtype, copy=False) 725 else: 726 weights = weights.astype(dtype, copy=False) 727 return np.average(coords, weights=weights, axis=0) 728 elif compound.lower() == 'residues': 729 compound_indices = atoms.resindices 730 n_compounds = atoms.n_residues 731 elif compound.lower() == 'segments': 732 compound_indices = atoms.segindices 733 n_compounds = atoms.n_segments 734 else: 735 raise ValueError("Unrecognized compound definition: {}\nPlease use" 736 " one of 'group', 'residues', or 'segments'." 737 "".format(compound)) 738 739 # Sort positions and weights by compound index and promote to dtype if 740 # required: 741 sort_indices = np.argsort(compound_indices) 742 compound_indices = compound_indices[sort_indices] 743 coords = atoms.positions[sort_indices] 744 if weights is None: 745 coords = coords.astype(dtype, copy=False) 746 else: 747 weights = weights.astype(dtype, copy=False) 748 weights = weights[sort_indices] 749 # Allocate output array: 750 centers = np.zeros((n_compounds, 3), dtype=dtype) 751 # Get sizes of compounds: 752 unique_compound_indices, compound_sizes = np.unique(compound_indices, 753 return_counts=True) 754 unique_compound_sizes = unique_int_1d(compound_sizes) 755 # Compute centers per compound for each compound size: 756 for compound_size in unique_compound_sizes: 757 compound_mask = compound_sizes == compound_size 758 _compound_indices = unique_compound_indices[compound_mask] 759 atoms_mask = np.in1d(compound_indices, _compound_indices) 760 _coords = coords[atoms_mask].reshape((-1, compound_size, 3)) 761 if weights is None: 762 _centers = _coords.mean(axis=1) 763 else: 764 _weights = weights[atoms_mask].reshape((-1, compound_size)) 765 _centers = (_coords * _weights[:, :, None]).sum(axis=1) 766 _centers /= _weights.sum(axis=1)[:, None] 767 centers[compound_mask] = _centers 768 if pbc: 769 centers = distances.apply_PBC(centers, atoms.dimensions) 770 return centers 771 772 @warn_if_not_unique 773 def center_of_geometry(self, pbc=None, compound='group'): 774 """Center of geometry (also known as centroid) of the group. 775 776 Computes the center of geometry (a.k.a. centroid) of 777 :class:`Atoms<Atom>` in the group. Centers of geometry per 778 :class:`Residue` or per :class:`Segment` can be obtained by setting the 779 `compound` parameter accordingly. 780 781 Parameters 782 ---------- 783 pbc : bool or None, optional 784 If ``True`` and `compound` is ``'group'``, move all atoms to the 785 primary unit cell before calculation. If ``True`` and `compound` is 786 ``'segments'`` or ``'residues'``, the center of each compound will 787 be calculated without moving any :class:`Atoms<Atom>` to keep the 788 compounds intact. Instead, the resulting position vectors will be 789 moved to the primary unit cell after calculation. 790 compound : {'group', 'segments', 'residues'}, optional 791 If ``'group'``, the center of geometry of all :class:`Atoms<Atom>` 792 in the group will be returned as a single position vector. Else, the 793 centers of geometry of each :class:`Segment` or :class:`Residue` 794 will be returned as an array of position vectors, i.e. a 2d array. 795 Note that, in any case, *only* the positions of :class:`Atoms<Atom>` 796 *belonging to the group* will be taken into account. 797 798 Returns 799 ------- 800 center : numpy.ndarray 801 Position vector(s) of the geometric center(s) of the group. 802 If `compound` was set to ``'group'``, the output will be a single 803 position vector. 804 If `compound` was set to ``'segments'`` or ``'residues'``, the 805 output will be a 2d array of shape ``(n, 3)`` where ``n`` is the 806 number of compounds. 807 808 Notes 809 ----- 810 If the :class:`MDAnalysis.core.flags` flag *use_pbc* is set to 811 ``True`` then the `pbc` keyword is used by default. 812 813 814 .. versionchanged:: 0.8 Added `pbc` keyword 815 .. versionchanged:: 0.19.0 Added `compound` parameter 816 """ 817 return self.center(None, pbc=pbc, compound=compound) 818 819 centroid = center_of_geometry 820 821 def bbox(self, **kwargs): 822 """Return the bounding box of the selection. 823 824 The lengths A,B,C of the orthorhombic enclosing box are :: 825 826 L = AtomGroup.bbox() 827 A,B,C = L[1] - L[0] 828 829 Parameters 830 ---------- 831 pbc : bool, optional 832 If ``True``, move all :class:`Atoms<Atom>` to the primary unit cell 833 before calculation. [``False``] 834 835 Returns 836 ------- 837 corners : numpy.ndarray 838 2x3 array giving corners of bounding box as 839 ``[[xmin, ymin, zmin], [xmax, ymax, zmax]]``. 840 841 Note 842 ---- 843 The :class:`MDAnalysis.core.flags` flag *use_pbc* when set to 844 ``True`` allows the *pbc* flag to be used by default. 845 846 847 .. versionadded:: 0.7.2 848 .. versionchanged:: 0.8 Added *pbc* keyword 849 """ 850 atomgroup = self.atoms 851 pbc = kwargs.pop('pbc', flags['use_pbc']) 852 853 if pbc: 854 x = atomgroup.pack_into_box(inplace=False) 855 else: 856 x = atomgroup.positions 857 858 return np.array([x.min(axis=0), x.max(axis=0)]) 859 860 def bsphere(self, **kwargs): 861 """Return the bounding sphere of the selection. 862 863 The sphere is calculated relative to the 864 :meth:`center of geometry<center_of_geometry>`. 865 866 Parameters 867 ---------- 868 pbc : bool, optional 869 If ``True``, move all atoms to the primary unit cell before 870 calculation. [``False``] 871 872 Returns 873 ------- 874 R : float 875 Radius of the bounding sphere. 876 center : numpy.ndarray 877 Coordinates of the sphere center as ``[xcen, ycen, zcen]``. 878 879 Note 880 ---- 881 The :class:`MDAnalysis.core.flags` flag *use_pbc* when set to 882 ``True`` allows the *pbc* flag to be used by default. 883 884 885 .. versionadded:: 0.7.3 886 .. versionchanged:: 0.8 Added *pbc* keyword 887 """ 888 atomgroup = self.atoms.unique 889 pbc = kwargs.pop('pbc', flags['use_pbc']) 890 891 if pbc: 892 x = atomgroup.pack_into_box(inplace=False) 893 centroid = atomgroup.center_of_geometry(pbc=True) 894 else: 895 x = atomgroup.positions 896 centroid = atomgroup.center_of_geometry(pbc=False) 897 898 R = np.sqrt(np.max(np.sum(np.square(x - centroid), axis=1))) 899 900 return R, centroid 901 902 def transform(self, M): 903 r"""Apply homogenous transformation matrix `M` to the coordinates. 904 905 :class:`Atom` coordinates are rotated and translated in-place. 906 907 Parameters 908 ---------- 909 M : array_like 910 4x4 matrix with the rotation in ``R = M[:3, :3]`` and the 911 translation in ``t = M[:3, 3]``. 912 913 Returns 914 ------- 915 self 916 917 See Also 918 -------- 919 MDAnalysis.lib.transformations : module of all coordinate transforms 920 921 Notes 922 ----- 923 The rotation :math:`\mathsf{R}` is about the origin and is applied 924 before the translation :math:`\mathbf{t}`: 925 926 .. math:: 927 928 \mathbf{x}' = \mathsf{R}\mathbf{x} + \mathbf{t} 929 930 """ 931 M = np.asarray(M) 932 R = M[:3, :3] 933 t = M[:3, 3] 934 return self.rotate(R, [0, 0, 0]).translate(t) 935 936 def translate(self, t): 937 r"""Apply translation vector `t` to the selection's coordinates. 938 939 :class:`Atom` coordinates are translated in-place. 940 941 Parameters 942 ---------- 943 t : array_like 944 vector to translate coordinates with 945 946 Returns 947 ------- 948 self 949 950 See Also 951 -------- 952 MDAnalysis.lib.transformations : module of all coordinate transforms 953 954 Notes 955 ----- 956 The method applies a translation to the :class:`AtomGroup` 957 from current coordinates :math:`\mathbf{x}` to new coordinates 958 :math:`\mathbf{x}'`: 959 960 .. math:: 961 962 \mathbf{x}' = \mathbf{x} + \mathbf{t} 963 964 """ 965 atomgroup = self.atoms.unique 966 vector = np.asarray(t) 967 # changes the coordinates in place 968 atomgroup.universe.trajectory.ts.positions[atomgroup.indices] += vector 969 return self 970 971 def rotate(self, R, point=(0, 0, 0)): 972 r"""Apply a rotation matrix `R` to the selection's coordinates. 973 :math:`\mathsf{R}` is a 3x3 orthogonal matrix that transforms a vector 974 :math:`\mathbf{x} \rightarrow \mathbf{x}'`: 975 976 .. math:: 977 978 \mathbf{x}' = \mathsf{R}\mathbf{x} 979 980 :class:`Atom` coordinates are rotated in-place. 981 982 Parameters 983 ---------- 984 R : array_like 985 3x3 rotation matrix 986 point : array_like, optional 987 Center of rotation 988 989 Returns 990 ------- 991 self 992 993 Notes 994 ----- 995 By default, rotates about the origin ``point=(0, 0, 0)``. To rotate 996 a group ``g`` around its center of geometry, use 997 ``g.rotate(R, point=g.center_of_geometry())``. 998 999 See Also 1000 -------- 1001 rotateby : rotate around given axis and angle 1002 MDAnalysis.lib.transformations : module of all coordinate transforms 1003 1004 """ 1005 R = np.asarray(R) 1006 point = np.asarray(point) 1007 1008 # changes the coordinates (in place) 1009 atomgroup = self.atoms.unique 1010 require_translation = bool(np.count_nonzero(point)) 1011 if require_translation: 1012 atomgroup.translate(-point) 1013 x = atomgroup.universe.trajectory.ts.positions 1014 idx = atomgroup.indices 1015 x[idx] = np.dot(x[idx], R.T) 1016 if require_translation: 1017 atomgroup.translate(point) 1018 1019 return self 1020 1021 def rotateby(self, angle, axis, point=None): 1022 r"""Apply a rotation to the selection's coordinates. 1023 1024 Parameters 1025 ---------- 1026 angle : float 1027 Rotation angle in degrees. 1028 axis : array_like 1029 Rotation axis vector. 1030 point : array_like, optional 1031 Center of rotation. If ``None`` then the center of geometry of this 1032 group is used. 1033 1034 Returns 1035 ------- 1036 self 1037 1038 Notes 1039 ----- 1040 The transformation from current coordinates :math:`\mathbf{x}` 1041 to new coordinates :math:`\mathbf{x}'` is 1042 1043 .. math:: 1044 1045 \mathbf{x}' = \mathsf{R}\,(\mathbf{x}-\mathbf{p}) + \mathbf{p} 1046 1047 where :math:`\mathsf{R}` is the rotation by `angle` around the 1048 `axis` going through `point` :math:`\mathbf{p}`. 1049 1050 See Also 1051 -------- 1052 MDAnalysis.lib.transformations.rotation_matrix : 1053 calculate :math:`\mathsf{R}` 1054 1055 """ 1056 alpha = np.radians(angle) 1057 axis = np.asarray(axis) 1058 if point is None: 1059 point = self.center_of_geometry() 1060 point = np.asarray(point) 1061 M = transformations.rotation_matrix(alpha, axis, point=point) 1062 return self.transform(M) 1063 1064 def pack_into_box(self, box=None, inplace=True): 1065 r"""Shift all :class:`Atoms<Atom>` in this group to the primary unit 1066 cell. 1067 1068 Parameters 1069 ---------- 1070 box : array_like 1071 Box dimensions, can be either orthogonal or triclinic information. 1072 Cell dimensions must be in an identical to format to those returned 1073 by :attr:`MDAnalysis.coordinates.base.Timestep.dimensions`, 1074 ``[lx, ly, lz, alpha, beta, gamma]``. If ``None``, uses these 1075 timestep dimensions. 1076 inplace : bool 1077 ``True`` to change coordinates in place. 1078 1079 Returns 1080 ------- 1081 coords : numpy.ndarray 1082 Shifted atom coordinates. 1083 1084 Notes 1085 ----- 1086 All atoms will be moved so that they lie between 0 and boxlength 1087 :math:`L_i` in all dimensions, i.e. the lower left corner of the 1088 simulation box is taken to be at (0,0,0): 1089 1090 .. math:: 1091 1092 x_i' = x_i - \left\lfloor\frac{x_i}{L_i}\right\rfloor 1093 1094 The default is to take unit cell information from the underlying 1095 :class:`~MDAnalysis.coordinates.base.Timestep` instance. The optional 1096 argument `box` can be used to provide alternative unit cell information 1097 (in the MDAnalysis standard format 1098 ``[Lx, Ly, Lz, alpha, beta, gamma]``). 1099 1100 Works with either orthogonal or triclinic box types. 1101 1102 .. note:: 1103 :meth:`AtomGroup.pack_into_box` is currently faster than 1104 :meth:`ResidueGroup.pack_into_box` or 1105 :meth:`SegmentGroup.pack_into_box`. 1106 1107 1108 .. versionadded:: 0.8 1109 """ 1110 # Try and auto detect box dimensions: 1111 if box is None: 1112 box = self.dimensions 1113 1114 # For a vector representation, the diagonal must not be zero, for a 1115 # [x, y, z, alpha, beta, gamma] representation, no element must be zero: 1116 if (box.shape == (3, 3) and (box.diagonal() == 0.0).any()) \ 1117 or (box == 0).any(): 1118 raise ValueError("One or more box dimensions is zero." 1119 " You can specify a box with 'box ='") 1120 # no matter what kind of group we have, we need to work on its atoms: 1121 ag = self.atoms 1122 # If self is unique, so are its atoms. Thus, if the group is unique, we 1123 # can use its atoms as is. 1124 if self.isunique: 1125 packed_coords = distances.apply_PBC(ag.positions, box) 1126 if inplace: 1127 ag.universe.trajectory.ts.positions[ag.ix] = packed_coords 1128 return ag.universe.trajectory.ts.positions[ag.ix] 1129 return packed_coords 1130 else: 1131 unique_ag = ag.unique 1132 unique_ix = unique_ag.ix 1133 restore_mask = ag._unique_restore_mask 1134 coords = ag.universe.trajectory.ts.positions[unique_ix] 1135 packed_coords = distances.apply_PBC(coords, box) 1136 if inplace: 1137 ag.universe.trajectory.ts.positions[unique_ix] = packed_coords 1138 return ag.universe.trajectory.ts.positions[unique_ix[restore_mask]] 1139 return packed_coords[restore_mask] 1140 1141 def wrap(self, compound="atoms", center="com", box=None): 1142 """Shift the contents of this group back into the unit cell. 1143 1144 This is a more powerful version of :meth:`pack_into_box`, allowing 1145 groups of :class:`Atoms<Atom>` to be kept together through the process. 1146 1147 Parameters 1148 ---------- 1149 compound : {'atoms', 'group', 'residues', 'segments', 'fragments'} 1150 The group which will be kept together through the shifting process. 1151 center : {'com', 'cog'} 1152 How to define the center of a given group of atoms. 1153 box : array_like 1154 Box dimensions, can be either orthogonal or triclinic information. 1155 Cell dimensions must be in an identical to format to those returned 1156 by :attr:`MDAnalysis.coordinates.base.Timestep.dimensions`, 1157 ``[lx, ly, lz, alpha, beta, gamma]``. If ``None``, uses these 1158 timestep dimensions. 1159 1160 Notes 1161 ----- 1162 When specifying a `compound`, the translation is calculated based on 1163 each compound. The same translation is applied to all atoms 1164 within this compound, meaning it will not be broken by the shift. 1165 This might however mean that not all atoms from the compound are 1166 inside the unit cell, but rather the center of the compound is. 1167 1168 `center` allows the definition of the center of each group to be 1169 specified. This can be either ``'com'`` for center of mass, or ``'cog'`` 1170 for center of geometry. 1171 1172 `box` allows a unit cell to be given for the transformation. If not 1173 specified, the :attr:`~MDAnalysis.coordinates.base.Timestep.dimensions` 1174 information from the current 1175 :class:`~MDAnalysis.coordinates.base.Timestep` will be used. 1176 1177 .. note:: 1178 :meth:`wrap` with all default keywords is identical to 1179 :meth:`pack_into_box` 1180 1181 1182 .. versionadded:: 0.9.2 1183 """ 1184 atomgroup = self.atoms.unique 1185 if compound.lower() == "atoms": 1186 return atomgroup.pack_into_box(box=box) 1187 1188 if compound.lower() == 'group': 1189 objects = [atomgroup.atoms] 1190 elif compound.lower() == 'residues': 1191 objects = atomgroup.residues 1192 elif compound.lower() == 'segments': 1193 objects = atomgroup.segments 1194 elif compound.lower() == 'fragments': 1195 objects = atomgroup.fragments 1196 else: 1197 raise ValueError("Unrecognized compound definition: {0}" 1198 "Please use one of 'group' 'residues' 'segments'" 1199 "or 'fragments'".format(compound)) 1200 1201 # TODO: ADD TRY-EXCEPT FOR MASSES PRESENCE 1202 if center.lower() in ('com', 'centerofmass'): 1203 centers = np.vstack([o.atoms.center_of_mass() for o in objects]) 1204 elif center.lower() in ('cog', 'centroid', 'centerofgeometry'): 1205 centers = np.vstack([o.atoms.center_of_geometry() for o in objects]) 1206 else: 1207 raise ValueError("Unrecognized center definition: {0}" 1208 "Please use one of 'com' or 'cog'".format(center)) 1209 centers = centers.astype(np.float32) 1210 1211 if box is None: 1212 box = atomgroup.dimensions 1213 1214 # calculate shift per object center 1215 dests = distances.apply_PBC(centers, box=box) 1216 shifts = dests - centers 1217 1218 for o, s in zip(objects, shifts): 1219 # Save some needless shifts 1220 if not all(s == 0.0): 1221 o.atoms.translate(s) 1222 1223 def copy(self): 1224 """Get another group identical to this one. 1225 1226 1227 .. versionadded:: 0.19.0 1228 """ 1229 group = self[:] 1230 # Try to fill the copied group's uniqueness caches: 1231 try: 1232 group._cache['isunique'] = self._cache['isunique'] 1233 if group._cache['isunique']: 1234 group._cache['unique'] = group 1235 except KeyError: 1236 pass 1237 return group 1238 1239 def groupby(self, topattrs): 1240 """Group together items in this group according to values of *topattr* 1241 1242 Parameters 1243 ---------- 1244 topattrs: str or list 1245 One or more topology attributes to group components by. 1246 Single arguments are passed as a string. Multiple arguments 1247 are passed as a list. 1248 1249 Returns 1250 ------- 1251 dict 1252 Unique values of the multiple combinations of topology attributes 1253 as keys, Groups as values. 1254 1255 Example 1256 ------- 1257 To group atoms with the same mass together: 1258 1259 >>> ag.groupby('masses') 1260 {12.010999999999999: <AtomGroup with 462 atoms>, 1261 14.007: <AtomGroup with 116 atoms>, 1262 15.999000000000001: <AtomGroup with 134 atoms>} 1263 1264 To group atoms with the same residue name and mass together: 1265 1266 >>> ag.groupby(['resnames', 'masses']) 1267 {('ALA', 1.008): <AtomGroup with 95 atoms>, 1268 ('ALA', 12.011): <AtomGroup with 57 atoms>, 1269 ('ALA', 14.007): <AtomGroup with 19 atoms>, 1270 ('ALA', 15.999): <AtomGroup with 19 atoms>}, 1271 ('ARG', 1.008): <AtomGroup with 169 atoms>, 1272 ...} 1273 1274 >>> ag.groupby(['resnames', 'masses'])('ALA', 15.999) 1275 <AtomGroup with 19 atoms> 1276 1277 1278 .. versionadded:: 0.16.0 1279 .. versionchanged:: 0.18.0 The function accepts multiple attributes 1280 """ 1281 1282 res = dict() 1283 1284 if isinstance(topattrs, (string_types, bytes)): 1285 attr=topattrs 1286 if isinstance(topattrs, bytes): 1287 attr = topattrs.decode('utf-8') 1288 ta = getattr(self, attr) 1289 1290 return {i: self[ta == i] for i in set(ta)} 1291 1292 else: 1293 attr = topattrs[0] 1294 ta = getattr(self, attr) 1295 for i in set(ta): 1296 if len(topattrs) == 1: 1297 res[i] = self[ta == i] 1298 else: 1299 res[i] = self[ta == i].groupby(topattrs[1:]) 1300 1301 return util.flatten_dict(res) 1302 1303 @_only_same_level 1304 def concatenate(self, other): 1305 """Concatenate with another Group or Component of the same level. 1306 1307 Duplicate entries and original order is preserved. It is synomymous to 1308 the `+` operator. 1309 1310 Parameters 1311 ---------- 1312 other : Group or Component 1313 Group or Component with `other.level` same as `self.level` 1314 1315 Returns 1316 ------- 1317 Group 1318 Group with elements of `self` and `other` concatenated 1319 1320 Example 1321 ------- 1322 The order of the original contents (including duplicates) 1323 are preserved when performing a concatenation. 1324 1325 >>> ag1 = u.select_atoms('name O') 1326 >>> ag2 = u.select_atoms('name N') 1327 >>> ag3 = ag1 + ag2 # or ag1.concatenate(ag2) 1328 >>> ag3[:3].names 1329 array(['O', 'O', 'O'], dtype=object) 1330 >>> ag3[-3:].names 1331 array(['N', 'N', 'N'], dtype=object) 1332 1333 1334 .. versionadded:: 0.16.0 1335 """ 1336 o_ix = other.ix_array 1337 return self._derived_class(np.concatenate([self.ix, o_ix]), 1338 self.universe) 1339 1340 @_only_same_level 1341 def union(self, other): 1342 """Group of elements either in this Group or another 1343 1344 On the contrary to concatenation, this method sort the elements and 1345 removes duplicate ones. It is synomymous to the `|` operator. 1346 1347 Parameters 1348 ---------- 1349 other : Group or Component 1350 Group or Component with `other.level` same as `self.level` 1351 1352 Returns 1353 ------- 1354 Group 1355 Group with the combined elements of `self` and `other`, without 1356 duplicate elements 1357 1358 Example 1359 ------- 1360 In contrast to :meth:`concatenate`, any duplicates are dropped 1361 and the result is sorted. 1362 1363 >>> ag1 = u.select_atoms('name O') 1364 >>> ag2 = u.select_atoms('name N') 1365 >>> ag3 = ag1 | ag2 # or ag1.union(ag2) 1366 >>> ag3[:3].names 1367 array(['N', 'O', 'N'], dtype=object) 1368 1369 See Also 1370 -------- 1371 concatenate, intersection 1372 1373 1374 .. versionadded:: 0.16 1375 """ 1376 o_ix = other.ix_array 1377 return self._derived_class(np.union1d(self.ix, o_ix), self.universe) 1378 1379 @_only_same_level 1380 def intersection(self, other): 1381 """Group of elements which are in both this Group and another 1382 1383 This method removes duplicate elements and sorts the result. It is 1384 synomymous to the `&` operator. 1385 1386 Parameters 1387 ---------- 1388 other : Group or Component 1389 Group or Component with `other.level` same as `self.level` 1390 1391 Returns 1392 ------- 1393 Group 1394 Group with the common elements of `self` and `other`, without 1395 duplicate elements 1396 1397 Example 1398 ------- 1399 Intersections can be used when the select atoms string would 1400 become too complicated. For example to find the water atoms 1401 which are within 4.0A of two segments: 1402 1403 >>> water = u.select_atoms('resname SOL') 1404 >>> shell1 = water.select_atoms('around 4.0 segid 1') 1405 >>> shell2 = water.select_atoms('around 4.0 segid 2') 1406 >>> common = shell1 & shell2 # or shell1.intersection(shell2) 1407 1408 See Also 1409 -------- 1410 union 1411 1412 1413 .. versionadded:: 0.16 1414 """ 1415 o_ix = other.ix_array 1416 return self._derived_class(np.intersect1d(self.ix, o_ix), self.universe) 1417 1418 @_only_same_level 1419 def subtract(self, other): 1420 """Group with elements from this Group that don't appear in other 1421 1422 The original order of this group is kept, as well as any duplicate 1423 elements. If an element of this Group is duplicated and appears in 1424 the other Group or Component, then all the occurences of that element 1425 are removed from the returned Group. 1426 1427 Parameters 1428 ---------- 1429 other : Group or Component 1430 Group or Component with `other.level` same as `self.level` 1431 1432 Returns 1433 ------- 1434 Group 1435 Group with the elements of `self` that are not in `other`, 1436 conserves order and duplicates. 1437 1438 Example 1439 ------- 1440 Unlike :meth:`difference` this method will not sort or remove 1441 duplicates. 1442 1443 >>> ag1 = u.atoms[[3, 3, 2, 2, 1, 1]] 1444 >>> ag2 = u.atoms[2] 1445 >>> ag3 = ag1 - ag2 # or ag1.subtract(ag2) 1446 >>> ag1.indices 1447 array([3, 3, 1, 1]) 1448 1449 See Also 1450 -------- 1451 concatenate, difference 1452 1453 1454 .. versionadded:: 0.16 1455 """ 1456 o_ix = other.ix_array 1457 in_other = np.in1d(self.ix, o_ix) # mask of in self.ix AND other 1458 return self[~in_other] # ie inverse of previous mask 1459 1460 @_only_same_level 1461 def difference(self, other): 1462 """Elements from this Group that do not appear in another 1463 1464 This method removes duplicate elements and sorts the result. As such, 1465 it is different from :meth:`subtract`. :meth:`difference` is synomymous 1466 to the `-` operator. 1467 1468 Parameters 1469 ---------- 1470 other : Group or Component 1471 Group or Component with `other.level` same as `self.level` 1472 1473 Returns 1474 ------- 1475 Group 1476 Group with the elements of `self` that are not in `other`, without 1477 duplicate elements 1478 1479 See Also 1480 -------- 1481 subtract, symmetric_difference 1482 1483 1484 .. versionadded:: 0.16 1485 """ 1486 o_ix = other.ix_array 1487 return self._derived_class(np.setdiff1d(self._ix, o_ix), self._u) 1488 1489 @_only_same_level 1490 def symmetric_difference(self, other): 1491 """Group of elements which are only in one of this Group or another 1492 1493 This method removes duplicate elements and the result is sorted. It is 1494 synomym to the `^` operator. 1495 1496 Parameters 1497 ---------- 1498 other : Group or Component 1499 Group or Component with `other.level` same as `self.level` 1500 1501 Returns 1502 ------- 1503 Group 1504 Group with the elements that are in `self` or in `other` but not in 1505 both, without duplicate elements 1506 1507 Example 1508 ------- 1509 1510 >>> ag1 = u.atoms[[0, 1, 5, 3, 3, 2]] 1511 >>> ag2 = u.atoms[[4, 4, 6, 2, 3, 5]] 1512 >>> ag3 = ag1 ^ ag2 # or ag1.symmetric_difference(ag2) 1513 >>> ag3.indices # 0 and 1 are only in ag1, 4 and 6 are only in ag2 1514 [0, 1, 4, 6] 1515 1516 See Also 1517 -------- 1518 difference 1519 1520 1521 .. versionadded:: 0.16 1522 """ 1523 o_ix = other.ix_array 1524 return self._derived_class(np.setxor1d(self._ix, o_ix), self._u) 1525 1526 def isdisjoint(self, other): 1527 """If the Group has no elements in common with the other Group 1528 1529 Parameters 1530 ---------- 1531 other : Group or Component 1532 Group or Component with `other.level` same as `self.level` 1533 1534 Returns 1535 ------- 1536 bool 1537 ``True`` if the two Groups do not have common elements 1538 1539 1540 .. versionadded:: 0.16 1541 """ 1542 return len(self.intersection(other)) == 0 1543 1544 @_only_same_level 1545 def issubset(self, other): 1546 """If all elements of this Group are part of another Group 1547 1548 Note that an empty group is a subset of any group of the same level. 1549 1550 Parameters 1551 ---------- 1552 other : Group or Component 1553 Group or Component with `other.level` same as `self.level` 1554 1555 Returns 1556 ------- 1557 bool 1558 ``True`` if this Group is a subset of the other one 1559 1560 1561 .. versionadded:: 0.16 1562 """ 1563 o_ix = set(other.ix_array) 1564 s_ix = set(self.ix) 1565 return s_ix.issubset(o_ix) 1566 1567 def is_strict_subset(self, other): 1568 """If this Group is a subset of another Group but not identical 1569 1570 Parameters 1571 ---------- 1572 other : Group or Component 1573 Group or Component with `other.level` same as `self.level` 1574 1575 Returns 1576 ------- 1577 bool 1578 ``True`` if this Group is a strict subset of the other one 1579 1580 1581 .. versionadded:: 0.16 1582 """ 1583 return self.issubset(other) and not self == other 1584 1585 @_only_same_level 1586 def issuperset(self, other): 1587 """If all elements of another Group are part of this Group 1588 1589 Parameters 1590 ---------- 1591 other : Group or Component 1592 Group or Component with `other.level` same as `self.level` 1593 1594 Returns 1595 ------- 1596 bool 1597 ``True`` if this Group is a subset of the other one 1598 1599 1600 .. versionadded:: 0.16 1601 """ 1602 o_ix = set(other.ix_array) 1603 s_ix = set(self.ix) 1604 return s_ix.issuperset(o_ix) 1605 1606 def is_strict_superset(self, other): 1607 """If this Group is a superset of another Group but not identical 1608 1609 Parameters 1610 ---------- 1611 other : Group or Component 1612 Group or Component with `other.level` same as `self.level` 1613 1614 Returns 1615 ------- 1616 bool 1617 ``True`` if this Group is a strict superset of the other one 1618 1619 1620 .. versionadded:: 0.16 1621 """ 1622 return self.issuperset(other) and not self == other 1623 1624 1625class AtomGroup(GroupBase): 1626 """An ordered array of atoms. 1627 1628 Can be initiated from an iterable of :class:`Atoms<Atom>`:: 1629 1630 ag = AtomGroup([Atom1, Atom2, Atom3]) 1631 1632 Or from providing a list of indices and the 1633 :class:`~MDAnalysis.core.universe.Universe` it should belong to:: 1634 1635 ag = AtomGroup([72, 14, 25], u) 1636 1637 Alternatively, an :class:`AtomGroup` is generated by indexing/slicing 1638 another :class:`AtomGroup`, such as the group of all :class:`Atoms<Atom>` in 1639 the :class:`~MDAnalysis.core.universe.Universe` at 1640 :attr:`MDAnalysis.core.universe.Universe.atoms`. 1641 1642 An :class:`AtomGroup` can be indexed and sliced like a list:: 1643 1644 ag[0], ag[-1] 1645 1646 will return the first and the last :class:`Atom` in the group whereas the 1647 slice:: 1648 1649 ag[0:6:2] 1650 1651 returns an :class:`AtomGroup` of every second element, corresponding to 1652 indices 0, 2, and 4. 1653 1654 It also supports "advanced slicing" when the argument is a 1655 :class:`numpy.ndarray` or a :class:`list`:: 1656 1657 aslice = [0, 3, -1, 10, 3] 1658 ag[aslice] 1659 1660 will return a new :class:`AtomGroup` of :class:`Atoms<Atom>` with those 1661 indices in the old :class:`AtomGroup`. 1662 1663 Finally, :class:`AtomGroups<AtomGroup>` can be created from a selection. 1664 See :meth:`select_atoms`. 1665 1666 .. note:: 1667 1668 :class:`AtomGroups<AtomGroup>` originating from a selection are sorted 1669 and duplicate elements are removed. This is not true for 1670 :class:`AtomGroups<AtomGroup>` produced by slicing. Thus, slicing can be 1671 used when the order of atoms is crucial (for instance, in order to 1672 define angles or dihedrals). 1673 1674 :class:`AtomGroups<AtomGroup>` can be compared and combined using group 1675 operators. For instance, :class:`AtomGroups<AtomGroup>` can be concatenated 1676 using `+` or :meth:`concatenate`:: 1677 1678 ag_concat = ag1 + ag2 # or ag_concat = ag1.concatenate(ag2) 1679 1680 When groups are concatenated, the order of the :class:`Atoms<Atom>` is 1681 conserved. If :class:`Atoms<Atom>` appear several times in one of the 1682 groups, the duplicates are kept in the resulting group. On the contrary to 1683 :meth:`concatenate`, :meth:`union` treats the :class:`AtomGroups<AtomGroup>` 1684 as sets so that duplicates are removed from the resulting group, and 1685 :class:`Atoms<Atom>` are ordered. The `|` operator is synomymous to 1686 :meth:`union`:: 1687 1688 ag_union = ag1 | ag2 # or ag_union = ag1.union(ag2) 1689 1690 The opposite operation to :meth:`concatenate` is :meth:`subtract`. This 1691 method creates a new group with all the :class:`Atoms<Atom>` of the group 1692 that are not in a given other group; the order of the :class:`Atoms<Atom>` 1693 is kept, and so are duplicates. :meth:`difference` is the set version of 1694 :meth:`subtract`. The resulting group is sorted and deduplicated. 1695 1696 All set methods are listed in the table below. These methods treat the 1697 groups as sorted and deduplicated sets of :class:`Atoms<Atom>`. 1698 1699 +-------------------------------+------------+----------------------------+ 1700 | Operation | Equivalent | Result | 1701 +===============================+============+============================+ 1702 | ``s.isdisjoint(t)`` | | ``True`` if ``s`` and | 1703 | | | ``t`` do not share | 1704 | | | elements | 1705 +-------------------------------+------------+----------------------------+ 1706 | ``s.issubset(t)`` | | test if all elements of | 1707 | | | ``s`` are part of ``t`` | 1708 +-------------------------------+------------+----------------------------+ 1709 | ``s.is_strict_subset(t)`` | | test if all elements of | 1710 | | | ``s`` are part of ``t``, | 1711 | | | and ``s != t`` | 1712 +-------------------------------+------------+----------------------------+ 1713 | ``s.issuperset(t)`` | | test if all elements of | 1714 | | | ``t`` are part of ``s`` | 1715 +-------------------------------+------------+----------------------------+ 1716 | ``s.is_strict_superset(t)`` | | test if all elements of | 1717 | | | ``t`` are part of ``s``, | 1718 | | | and ``s != t`` | 1719 +-------------------------------+------------+----------------------------+ 1720 | ``s.union(t)`` | ``s | t`` | new Group with elements | 1721 | | | from both ``s`` and ``t`` | 1722 +-------------------------------+------------+----------------------------+ 1723 | ``s.intersection(t)`` | ``s & t`` | new Group with elements | 1724 | | | common to ``s`` and ``t`` | 1725 +-------------------------------+------------+----------------------------+ 1726 | ``s.difference(t)`` | ``s - t`` | new Group with elements of | 1727 | | | ``s`` that are not in ``t``| 1728 +-------------------------------+------------+----------------------------+ 1729 | ``s.symmetric_difference(t)`` | ``s ^ t`` | new Group with elements | 1730 | | | that are part of ``s`` or | 1731 | | | ``t`` but not both | 1732 +-------------------------------+------------+----------------------------+ 1733 1734 The following methods keep the order of the atoms as well as duplicates. 1735 1736 +-------------------------------+------------+----------------------------+ 1737 | Operation | Equivalent | Result | 1738 +===============================+============+============================+ 1739 | ``len(s)`` | | number of elements (atoms, | 1740 | | | residues or segment) in | 1741 | | | the group | 1742 +-------------------------------+------------+----------------------------+ 1743 | ``s == t`` | | test if ``s`` and ``t`` | 1744 | | | contain the same elements | 1745 | | | in the same order | 1746 +-------------------------------+------------+----------------------------+ 1747 | ``s.concatenate(t)`` | ``s + t`` | new Group with elements | 1748 | | | from ``s`` and from ``t`` | 1749 +-------------------------------+------------+----------------------------+ 1750 | ``s.subtract(t)`` | | new Group with elements | 1751 | | | from ``s`` that are not | 1752 | | | in ``t`` | 1753 +-------------------------------+------------+----------------------------+ 1754 1755 The `in` operator allows to test if an :class:`Atom` is in the 1756 :class:`AtomGroup`. 1757 1758 :class:`AtomGroup` instances are always bound to a 1759 :class:`MDAnalysis.core.universe.Universe`. They cannot exist in isolation. 1760 1761 .. rubric:: Deprecated functionality 1762 1763 *Instant selectors* will be removed in the 1.0 release. See issue `#1377 1764 <https://github.com/MDAnalysis/mdanalysis/issues/1377>`_ for more details. 1765 1766 :class:`Atoms<Atom>` can also be accessed in a Pythonic fashion by using the 1767 :class:`Atom` name as an attribute. For instance, :: 1768 1769 ag.CA 1770 1771 will provide an :class:`AtomGroup` of all CA :class:`Atoms<Atom>` in the 1772 group. These *instant selector* attributes are auto-generated for 1773 each atom name encountered in the group. 1774 1775 Notes 1776 ----- 1777 The name-attribute instant selector access to :class:`Atoms<Atom>` is mainly 1778 meant for quick interactive work. Thus it either returns a single 1779 :class:`Atom` if there is only one matching :class:`Atom`, *or* a 1780 new :class:`AtomGroup` for multiple matches. This makes it difficult to use 1781 the feature consistently in scripts. 1782 1783 1784 See Also 1785 -------- 1786 :class:`MDAnalysis.core.universe.Universe` 1787 1788 1789 .. deprecated:: 0.16.2 1790 *Instant selectors* of :class:`AtomGroup` will be removed in the 1.0 1791 release. See :ref:`Instant selectors <instance-selectors>` for details 1792 and alternatives. 1793 """ 1794 def __getitem__(self, item): 1795 # DEPRECATED in 0.16.2 1796 # REMOVE in 1.0 1797 # 1798 # u.atoms['HT1'] access, otherwise default 1799 if isinstance(item, string_types): 1800 try: 1801 return self._get_named_atom(item) 1802 except (AttributeError, selection.SelectionError): 1803 pass 1804 return super(AtomGroup, self).__getitem__(item) 1805 1806 def __getattr__(self, attr): 1807 # DEPRECATED in 0.16.2 1808 # REMOVE in 1.0 1809 # 1810 # is this a known attribute failure? 1811 if attr in ('fragments',): # TODO: Generalise this to cover many attributes 1812 # eg: 1813 # if attr in _ATTR_ERRORS: 1814 # raise NDE(_ATTR_ERRORS[attr]) 1815 raise NoDataError("AtomGroup has no fragments; this requires Bonds") 1816 elif hasattr(self.universe._topology, 'names'): 1817 # Ugly hack to make multiple __getattr__s work 1818 try: 1819 return self._get_named_atom(attr) 1820 except selection.SelectionError: 1821 pass 1822 raise AttributeError("{cls} has no attribute {attr}".format( 1823 cls=self.__class__.__name__, attr=attr)) 1824 1825 def __reduce__(self): 1826 return (_unpickle, (self.universe.anchor_name, self.ix)) 1827 1828 @property 1829 def atoms(self): 1830 """The :class:`AtomGroup` itself. 1831 1832 See Also 1833 -------- 1834 copy : return a true copy of the :class:`AtomGroup` 1835 1836 1837 .. versionchanged:: 0.19.0 1838 In previous versions, this returned a copy, but now 1839 the :class:`AtomGroup` itself is returned. This should 1840 not affect any code but only speed up calculations. 1841 1842 """ 1843 return self 1844 1845 @property 1846 def n_atoms(self): 1847 """Number of atoms in the :class:`AtomGroup`. 1848 1849 Equivalent to ``len(self)``. 1850 """ 1851 return len(self) 1852 1853 @property 1854 def residues(self): 1855 """A sorted :class:`ResidueGroup` of the unique 1856 :class:`Residues<Residue>` present in the :class:`AtomGroup`. 1857 """ 1858 rg = self.universe.residues[unique_int_1d(self.resindices)] 1859 rg._cache['isunique'] = True 1860 rg._cache['unique'] = rg 1861 return rg 1862 1863 @residues.setter 1864 def residues(self, new): 1865 # Can set with Res, ResGroup or list/tuple of Res 1866 if isinstance(new, Residue): 1867 r_ix = itertools.cycle((new.resindex,)) 1868 elif isinstance(new, ResidueGroup): 1869 r_ix = new.resindices 1870 else: 1871 try: 1872 r_ix = [r.resindex for r in new] 1873 except AttributeError: 1874 raise TypeError("Can only set AtomGroup residues to Residue " 1875 "or ResidueGroup not {}".format( 1876 ', '.join(type(r) for r in new 1877 if not isinstance(r, Residue)) 1878 )) 1879 if not isinstance(r_ix, itertools.cycle) and len(r_ix) != len(self): 1880 raise ValueError("Incorrect size: {} for AtomGroup of size: {}" 1881 "".format(len(new), len(self))) 1882 # Optimisation TODO: 1883 # This currently rebuilds the tt len(self) times 1884 # Ideally all changes would happen and *afterwards* tables are built 1885 # Alternatively, if the changes didn't rebuild table, this list 1886 # comprehension isn't terrible. 1887 for at, r in zip(self, r_ix): 1888 self.universe._topology.tt.move_atom(at.ix, r) 1889 1890 @property 1891 def n_residues(self): 1892 """Number of unique :class:`Residues<Residue>` present in the 1893 :class:`AtomGroup`. 1894 1895 Equivalent to ``len(self.residues)``. 1896 1897 """ 1898 return len(self.residues) 1899 1900 @property 1901 def segments(self): 1902 """A sorted :class:`SegmentGroup` of the unique segments present in the 1903 :class:`AtomGroup`. 1904 """ 1905 sg = self.universe.segments[unique_int_1d(self.segindices)] 1906 sg._cache['isunique'] = True 1907 sg._cache['unique'] = sg 1908 return sg 1909 1910 @segments.setter 1911 def segments(self, new): 1912 raise NotImplementedError("Cannot assign Segments to AtomGroup. " 1913 "Segments are assigned to Residues") 1914 1915 @property 1916 def n_segments(self): 1917 """Number of unique segments present in the :class:`AtomGroup`. 1918 1919 Equivalent to ``len(self.segments)``. 1920 """ 1921 return len(self.segments) 1922 1923 @property 1924 @cached('unique_restore_mask') 1925 def _unique_restore_mask(self): 1926 # The _unique_restore_mask property's cache is populated whenever the 1927 # AtomGroup.unique property of a *non-unique* AtomGroup is accessed. 1928 # If _unique_restore_mask is not cached, it is *definitely* used in the 1929 # wrong place, so we raise an exception here. In principle, the 1930 # exception should be an AttributeError, but the error message would 1931 # then be replaced by the __getattr__() error message. To prevent the 1932 # message from being overridden, we raise a RuntimeError instead. 1933 if self.isunique: 1934 msg = ("{0}._unique_restore_mask is not available if the {0} is " 1935 "unique. ".format(self.__class__.__name__)) 1936 else: 1937 msg = ("{0}._unique_restore_mask is only available after " 1938 "accessing {0}.unique. ".format(self.__class__.__name__)) 1939 msg += ("If you see this error message in an unmodified release " 1940 "version of MDAnalysis, this is almost certainly a bug!") 1941 raise RuntimeError(msg) 1942 1943 @_unique_restore_mask.setter 1944 def _unique_restore_mask(self, mask): 1945 self._cache['unique_restore_mask'] = mask 1946 1947 @property 1948 @cached('unique') 1949 def unique(self): 1950 """An :class:`AtomGroup` containing sorted and unique 1951 :class:`Atoms<Atom>` only. 1952 1953 If the :class:`AtomGroup` is unique, this is the group itself. 1954 1955 Examples 1956 -------- 1957 1958 >>> ag = u.atoms[[2, 1, 2, 2, 1, 0]] 1959 >>> ag 1960 <AtomGroup with 6 atoms> 1961 >>> ag.ix 1962 array([2, 1, 2, 2, 1, 0]) 1963 >>> ag2 = ag.unique 1964 >>> ag2 1965 <AtomGroup with 3 atoms> 1966 >>> ag2.ix 1967 array([0, 1, 2]) 1968 >>> ag2.unique is ag2 1969 True 1970 1971 1972 .. versionadded:: 0.16.0 1973 .. versionchanged:: 0.19.0 If the :class:`AtomGroup` is already unique, 1974 :attr:`AtomGroup.unique` now returns the group itself instead of a 1975 copy. 1976 """ 1977 if self.isunique: 1978 return self 1979 unique_ix, restore_mask = np.unique(self.ix, return_inverse=True) 1980 _unique = self.universe.atoms[unique_ix] 1981 self._unique_restore_mask = restore_mask 1982 # Since we know that _unique is a unique AtomGroup, we set its 1983 # uniqueness caches from here: 1984 _unique._cache['isunique'] = True 1985 _unique._cache['unique'] = _unique 1986 return _unique 1987 1988 @property 1989 def positions(self): 1990 """Coordinates of the :class:`Atoms<Atom>` in the :class:`AtomGroup`. 1991 1992 A :class:`numpy.ndarray` with 1993 :attr:`~numpy.ndarray.shape`\ ``=(``\ :attr:`~AtomGroup.n_atoms`\ ``, 3)`` 1994 and :attr:`~numpy.ndarray.dtype`\ ``=numpy.float32``. 1995 1996 The positions can be changed by assigning an array of the appropriate 1997 shape, i.e., either ``(``\ :attr:`~AtomGroup.n_atoms`\ ``, 3)`` to 1998 assign individual coordinates, or ``(3,)`` to assign the *same* 1999 coordinate to all :class:`Atoms<Atom>` (e.g., 2000 ``ag.positions = array([0,0,0])`` will move all :class:`Atoms<Atom>` 2001 to the origin). 2002 2003 .. note:: Changing positions is not reflected in any files; reading any 2004 frame from the 2005 :attr:`~MDAnalysis.core.universe.Universe.trajectory` will 2006 replace the change with that from the file *except* if the 2007 :attr:`~MDAnalysis.core.universe.Universe.trajectory` is held 2008 in memory, e.g., when the 2009 :meth:`~MDAnalysis.core.universe.Universe.transfer_to_memory` 2010 method was used. 2011 2012 Raises 2013 ------ 2014 ~MDAnalysis.exceptions.NoDataError 2015 If the underlying :class:`~MDAnalysis.coordinates.base.Timestep` 2016 does not contain 2017 :attr:`~MDAnalysis.coordinates.base.Timestep.positions`. 2018 """ 2019 return self.universe.trajectory.ts.positions[self.ix] 2020 2021 @positions.setter 2022 def positions(self, values): 2023 ts = self.universe.trajectory.ts 2024 ts.positions[self.ix, :] = values 2025 2026 @property 2027 def velocities(self): 2028 """Velocities of the :class:`Atoms<Atom>` in the :class:`AtomGroup`. 2029 2030 A :class:`numpy.ndarray` with 2031 :attr:`~numpy.ndarray.shape`\ ``=(``\ :attr:`~AtomGroup.n_atoms`\ ``, 3)`` 2032 and :attr:`~numpy.ndarray.dtype`\ ``=numpy.float32``. 2033 2034 The velocities can be changed by assigning an array of the appropriate 2035 shape, i.e. either ``(``\ :attr:`~AtomGroup.n_atoms`\ ``, 3)`` to assign 2036 individual velocities or ``(3,)`` to assign the *same* velocity to all 2037 :class:`Atoms<Atom>` (e.g. ``ag.velocities = array([0,0,0])`` will give 2038 all :class:`Atoms<Atom>` zero :attr:`~Atom.velocity`). 2039 2040 Raises 2041 ------ 2042 ~MDAnalysis.exceptions.NoDataError 2043 If the underlying :class:`~MDAnalysis.coordinates.base.Timestep` 2044 does not contain 2045 :attr:`~MDAnalysis.coordinates.base.Timestep.velocities`. 2046 """ 2047 ts = self.universe.trajectory.ts 2048 try: 2049 return np.array(ts.velocities[self.ix]) 2050 except (AttributeError, NoDataError): 2051 raise NoDataError("Timestep does not contain velocities") 2052 2053 @velocities.setter 2054 def velocities(self, values): 2055 ts = self.universe.trajectory.ts 2056 try: 2057 ts.velocities[self.ix, :] = values 2058 except (AttributeError, NoDataError): 2059 raise NoDataError("Timestep does not contain velocities") 2060 2061 @property 2062 def forces(self): 2063 """Forces on each :class:`Atom` in the :class:`AtomGroup`. 2064 2065 A :class:`numpy.ndarray` with 2066 :attr:`~numpy.ndarray.shape`\ ``=(``\ :attr:`~AtomGroup.n_atoms`\ ``, 3)`` 2067 and :attr:`~numpy.ndarray.dtype`\ ``=numpy.float32``. 2068 2069 The forces can be changed by assigning an array of the appropriate 2070 shape, i.e. either ``(``\ :attr:`~AtomGroup.n_atoms`\ ``, 3)`` to assign 2071 individual forces or ``(3,)`` to assign the *same* force to all 2072 :class:`Atoms<Atom>` (e.g. ``ag.forces = array([0,0,0])`` will give all 2073 :class:`Atoms<Atom>` a zero :attr:`~Atom.force`). 2074 2075 Raises 2076 ------ 2077 ~MDAnalysis.exceptions.NoDataError 2078 If the :class:`~MDAnalysis.coordinates.base.Timestep` does not 2079 contain :attr:`~MDAnalysis.coordinates.base.Timestep.forces`. 2080 """ 2081 ts = self.universe.trajectory.ts 2082 try: 2083 return ts.forces[self.ix] 2084 except (AttributeError, NoDataError): 2085 raise NoDataError("Timestep does not contain forces") 2086 2087 @forces.setter 2088 def forces(self, values): 2089 ts = self.universe.trajectory.ts 2090 try: 2091 ts.forces[self.ix, :] = values 2092 except (AttributeError, NoDataError): 2093 raise NoDataError("Timestep does not contain forces") 2094 2095 @property 2096 def ts(self): 2097 """Temporary Timestep that contains the selection coordinates. 2098 2099 A :class:`~MDAnalysis.coordinates.base.Timestep` instance, 2100 which can be passed to a trajectory writer. 2101 2102 If :attr:`~AtomGroup.ts` is modified then these modifications 2103 will be present until the frame number changes (which 2104 typically happens when the underlying 2105 :attr:`~MDAnalysis.core.universe.Universe.trajectory` frame changes). 2106 2107 It is not possible to assign a new 2108 :class:`~MDAnalysis.coordinates.base.Timestep` to the 2109 :attr:`AtomGroup.ts` attribute; change attributes of the object. 2110 """ 2111 trj_ts = self.universe.trajectory.ts # original time step 2112 2113 return trj_ts.copy_slice(self.indices) 2114 2115 # As with universe.select_atoms, needing to fish out specific kwargs 2116 # (namely, 'updating') doesn't allow a very clean signature. 2117 def select_atoms(self, sel, *othersel, **selgroups): 2118 """Select :class:`Atoms<Atom>` using a selection string. 2119 2120 Returns an :class:`AtomGroup` with :class:`Atoms<Atom>` sorted according 2121 to their index in the topology (this is to ensure that there are no 2122 duplicates, which can happen with complicated selections). 2123 2124 Raises 2125 ------ 2126 TypeError 2127 If the arbitrary groups passed are not of type 2128 :class:`MDAnalysis.core.groups.AtomGroup` 2129 2130 Examples 2131 -------- 2132 All simple selection listed below support multiple arguments which are 2133 implicitly combined with an or operator. For example 2134 2135 >>> sel = universe.select_atoms('resname MET GLY') 2136 2137 is equivalent to 2138 2139 >>> sel = universe.select_atoms('resname MET or resname GLY') 2140 2141 Will select all atoms with a residue name of either MET or GLY. 2142 2143 Subselections can be grouped with parentheses. 2144 2145 >>> sel = universe.select_atoms("segid DMPC and not ( name H* O* )") 2146 >>> sel 2147 <AtomGroup with 3420 atoms> 2148 2149 2150 Existing :class:`AtomGroup` objects can be passed as named arguments, 2151 which will then be available to the selection parser. 2152 2153 >>> universe.select_atoms("around 10 group notHO", notHO=sel) 2154 <AtomGroup with 1250 atoms> 2155 2156 Selections can be set to update automatically on frame change, by 2157 setting the `updating` keyword argument to `True`. This will return 2158 a :class:`UpdatingAtomGroup` which can represent the solvation shell 2159 around another object. 2160 2161 >>> universe.select_atoms("resname SOL and around 2.0 protein", updating=True) 2162 <Updating AtomGroup with 100 atoms> 2163 2164 Notes 2165 ----- 2166 2167 If exact ordering of atoms is required (for instance, for 2168 :meth:`~AtomGroup.angle` or :meth:`~AtomGroup.dihedral` calculations) 2169 then one supplies selections *separately* in the required order. Also, 2170 when multiple :class:`AtomGroup` instances are concatenated with the 2171 ``+`` operator, then the order of :class:`Atom` instances is preserved 2172 and duplicates are *not* removed. 2173 2174 2175 See Also 2176 -------- 2177 :ref:`selection-commands-label` for further details and examples. 2178 2179 2180 .. rubric:: Selection syntax 2181 2182 2183 The selection parser understands the following CASE SENSITIVE 2184 *keywords*: 2185 2186 **Simple selections** 2187 2188 protein, backbone, nucleic, nucleicbackbone 2189 selects all atoms that belong to a standard set of residues; 2190 a protein is identfied by a hard-coded set of residue names so 2191 it may not work for esoteric residues. 2192 segid *seg-name* 2193 select by segid (as given in the topology), e.g. ``segid 4AKE`` 2194 or ``segid DMPC`` 2195 resid *residue-number-range* 2196 resid can take a single residue number or a range of numbers. A 2197 range consists of two numbers separated by a colon (inclusive) 2198 such as ``resid 1:5``. A residue number ("resid") is taken 2199 directly from the topology. 2200 If icodes are present in the topology, then these will be 2201 taken into account. Ie 'resid 163B' will only select resid 2202 163 with icode B while 'resid 163' will select only residue 163. 2203 Range selections will also respect icodes, so 'resid 162-163B' 2204 will select all residues in 162 and those in 163 up to icode B. 2205 resnum *resnum-number-range* 2206 resnum is the canonical residue number; typically it is set to 2207 the residue id in the original PDB structure. 2208 resname *residue-name* 2209 select by residue name, e.g. ``resname LYS`` 2210 name *atom-name* 2211 select by atom name (as given in the topology). Often, this is 2212 force field dependent. Example: ``name CA`` (for Cα atoms) 2213 or ``name OW`` (for SPC water oxygen) 2214 type *atom-type* 2215 select by atom type; this is either a string or a number and 2216 depends on the force field; it is read from the topology file 2217 (e.g. the CHARMM PSF file contains numeric atom types). It has 2218 non-sensical values when a PDB or GRO file is used as a topology 2219 atom *seg-name* *residue-number* *atom-name* 2220 a selector for a single atom consisting of segid resid atomname, 2221 e.g. ``DMPC 1 C2`` selects the C2 carbon of the first residue of 2222 the DMPC segment 2223 altloc *alternative-location* 2224 a selection for atoms where alternative locations are available, 2225 which is often the case with high-resolution crystal structures 2226 e.g. `resid 4 and resname ALA and altloc B` selects only the 2227 atoms of ALA-4 that have an altloc B record. 2228 moltype *molecule-type* 2229 select by molecule type, e.g. ``moltype Protein_A``. At the 2230 moment, only the TPR format defines the molecule type. 2231 2232 **Boolean** 2233 2234 not 2235 all atoms not in the selection, e.g. ``not protein`` selects 2236 all atoms that aren't part of a protein 2237 2238 and, or 2239 combine two selections according to the rules of boolean 2240 algebra, e.g. ``protein and not resname ALA LYS`` 2241 selects all atoms that belong to a protein, but are not in a 2242 lysine or alanine residue 2243 2244 **Geometric** 2245 2246 around *distance* *selection* 2247 selects all atoms a certain cutoff away from another selection, 2248 e.g. ``around 3.5 protein`` selects all atoms not belonging to 2249 protein that are within 3.5 Angstroms from the protein 2250 point *x* *y* *z* *distance* 2251 selects all atoms within a cutoff of a point in space, make sure 2252 coordinate is separated by spaces, 2253 e.g. ``point 5.0 5.0 5.0 3.5`` selects all atoms within 3.5 2254 Angstroms of the coordinate (5.0, 5.0, 5.0) 2255 prop [abs] *property* *operator* *value* 2256 selects atoms based on position, using *property* **x**, **y**, 2257 or **z** coordinate. Supports the **abs** keyword (for absolute 2258 value) and the following *operators*: **<, >, <=, >=, ==, !=**. 2259 For example, ``prop z >= 5.0`` selects all atoms with z 2260 coordinate greater than 5.0; ``prop abs z <= 5.0`` selects all 2261 atoms within -5.0 <= z <= 5.0. 2262 sphzone *radius* *selection* 2263 Selects all atoms that are within *radius* of the center of 2264 geometry of *selection* 2265 sphlayer *inner radius* *outer radius* *selection* 2266 Similar to sphzone, but also excludes atoms that are within 2267 *inner radius* of the selection COG 2268 cyzone *externalRadius* *zMax* *zMin* *selection* 2269 selects all atoms within a cylindric zone centered in the 2270 center of geometry (COG) of a given selection, 2271 e.g. ``cyzone 15 4 -8 protein and resid 42`` selects the 2272 center of geometry of protein and resid 42, and creates a 2273 cylinder of external radius 15 centered on the COG. In z, the 2274 cylinder extends from 4 above the COG to 8 below. Positive 2275 values for *zMin*, or negative ones for *zMax*, are allowed. 2276 cylayer *innerRadius* *externalRadius* *zMax* *zMin* *selection* 2277 selects all atoms within a cylindric layer centered in the 2278 center of geometry (COG) of a given selection, 2279 e.g. ``cylayer 5 10 10 -8 protein`` selects the center of 2280 geometry of protein, and creates a cylindrical layer of inner 2281 radius 5, external radius 10 centered on the COG. In z, the 2282 cylinder extends from 10 above the COG to 8 below. Positive 2283 values for *zMin*, or negative ones for *zMax*, are allowed. 2284 2285 **Connectivity** 2286 2287 byres *selection* 2288 selects all atoms that are in the same segment and residue as 2289 selection, e.g. specify the subselection after the byres keyword 2290 bonded *selection* 2291 selects all atoms that are bonded to selection 2292 eg: ``select name H and bonded name O`` selects only hydrogens 2293 bonded to oxygens 2294 2295 **Index** 2296 2297 bynum *index-range* 2298 selects all atoms within a range of (1-based) inclusive indices, 2299 e.g. ``bynum 1`` selects the first atom in the universe; 2300 ``bynum 5:10`` selects atoms 5 through 10 inclusive. All atoms 2301 in the :class:`~MDAnalysis.core.universe.Universe` are 2302 consecutively numbered, and the index runs from 1 up to the 2303 total number of atoms. 2304 2305 **Preexisting selections** 2306 2307 group `group-name` 2308 selects the atoms in the :class:`AtomGroup` passed to the 2309 function as an argument named `group-name`. Only the atoms 2310 common to `group-name` and the instance 2311 :meth:`~MDAnalysis.core.groups.AtomGroup.select_atoms` 2312 was called from will be considered, unless ``group`` is 2313 preceded by the ``global`` keyword. `group-name` will be 2314 included in the parsing just by comparison of atom indices. 2315 This means that it is up to the user to make sure the 2316 `group-name` group was defined in an appropriate 2317 :class:`~MDAnalysis.core.universe.Universe`. 2318 2319 global *selection* 2320 by default, when issuing 2321 :meth:`~MDAnalysis.core.groups.AtomGroup.select_atoms` from an 2322 :class:`~MDAnalysis.core.groups.AtomGroup`, selections and 2323 subselections are returned intersected with the atoms of that 2324 instance. Prefixing a selection term with ``global`` causes its 2325 selection to be returned in its entirety. As an example, the 2326 ``global`` keyword allows for 2327 ``lipids.select_atoms("around 10 global protein")`` --- where 2328 ``lipids`` is a group that does not contain any proteins. Were 2329 ``global`` absent, the result would be an empty selection since 2330 the ``protein`` subselection would itself be empty. When issuing 2331 :meth:`~MDAnalysis.core.groups.AtomGroup.select_atoms` from a 2332 :class:`~MDAnalysis.core.universe.Universe`, ``global`` is 2333 ignored. 2334 2335 **Dynamic selections** 2336 If :meth:`~MDAnalysis.core.groups.AtomGroup.select_atoms` is 2337 invoked with named argument `updating` set to `True`, an 2338 :class:`~MDAnalysis.core.groups.UpdatingAtomGroup` instance will be 2339 returned, instead of a regular 2340 :class:`~MDAnalysis.core.groups.AtomGroup`. It behaves just like 2341 the latter, with the difference that the selection expressions are 2342 re-evaluated every time the trajectory frame changes (this happens 2343 lazily, only when the 2344 :class:`~MDAnalysis.core.groups.UpdatingAtomGroup` is accessed so 2345 that there is no redundant updating going on). 2346 Issuing an updating selection from an already updating group will 2347 cause later updates to also reflect the updating of the base group. 2348 A non-updating selection or a slicing operation made on an 2349 :class:`~MDAnalysis.core.groups.UpdatingAtomGroup` will return a 2350 static :class:`~MDAnalysis.core.groups.AtomGroup`, which will no 2351 longer update across frames. 2352 2353 2354 .. versionchanged:: 0.7.4 Added *resnum* selection. 2355 .. versionchanged:: 0.8.1 Added *group* and *fullgroup* selections. 2356 .. deprecated:: 0.11 The use of *fullgroup* has been deprecated in favor 2357 of the equivalent *global group* selections. 2358 .. versionchanged:: 0.13.0 Added *bonded* selection. 2359 .. versionchanged:: 0.16.0 Resid selection now takes icodes into account 2360 where present. 2361 .. versionchanged:: 0.16.0 Updating selections now possible by setting 2362 the `updating` argument. 2363 .. versionchanged:: 0.17.0 Added *moltype* and *molnum* selections. 2364 .. versionchanged:: 0.19.0 2365 Added strict type checking for passed groups. 2366 Added periodic kwarg (default True) 2367 """ 2368 # once flags removed, replace with default=True 2369 periodic = selgroups.pop('periodic', flags['use_periodic_selections']) 2370 2371 updating = selgroups.pop('updating', False) 2372 sel_strs = (sel,) + othersel 2373 2374 for group, thing in selgroups.items(): 2375 if not isinstance(thing, AtomGroup): 2376 raise TypeError("Passed groups must be AtomGroups. " 2377 "You provided {} for group '{}'".format( 2378 thing.__class__.__name__, group)) 2379 2380 selections = tuple((selection.Parser.parse(s, selgroups, periodic=periodic) 2381 for s in sel_strs)) 2382 if updating: 2383 atomgrp = UpdatingAtomGroup(self, selections, sel_strs) 2384 else: 2385 # Apply the first selection and sum to it 2386 atomgrp = sum([sel.apply(self) for sel in selections[1:]], 2387 selections[0].apply(self)) 2388 return atomgrp 2389 2390 def split(self, level): 2391 """Split :class:`AtomGroup` into a :class:`list` of 2392 :class:`AtomGroups<AtomGroup>` by `level`. 2393 2394 Parameters 2395 ---------- 2396 level : {'atom', 'residue', 'molecule', 'segment'} 2397 2398 2399 .. versionadded:: 0.9.0 2400 .. versionchanged:: 0.17.0 Added the 'molecule' level. 2401 """ 2402 accessors = {'segment': 'segindices', 2403 'residue': 'resindices', 2404 'molecule': 'molnums'} 2405 2406 if level == "atom": 2407 return [self.universe.atoms[[a.ix]] for a in self] 2408 2409 # higher level groupings 2410 try: 2411 levelindices = getattr(self, accessors[level]) 2412 except AttributeError: 2413 raise AttributeError('This universe does not have {} ' 2414 'information. Maybe it is not provided in the ' 2415 'topology format in use.'.format(level)) 2416 except KeyError: 2417 raise ValueError("level = '{0}' not supported, " 2418 "must be one of {1}".format(level, 2419 accessors.keys())) 2420 2421 return [self[levelindices == index] for index in 2422 unique_int_1d(levelindices)] 2423 2424 def guess_bonds(self, vdwradii=None): 2425 """Guess bonds that exist within this :class:`AtomGroup` and add them to 2426 the underlying :attr:`~AtomGroup.universe`. 2427 2428 Parameters 2429 ---------- 2430 vdwradii : dict, optional 2431 Dict relating atom types: vdw radii 2432 2433 2434 See Also 2435 -------- 2436 :func:`MDAnalysis.topology.guessers.guess_bonds` 2437 2438 2439 .. versionadded:: 0.10.0 2440 """ 2441 from ..topology.core import guess_bonds, guess_angles, guess_dihedrals 2442 from .topologyattrs import Bonds, Angles, Dihedrals 2443 2444 def get_TopAttr(u, name, cls): 2445 """either get *name* or create one from *cls*""" 2446 try: 2447 return getattr(u._topology, name) 2448 except AttributeError: 2449 attr = cls([]) 2450 u.add_TopologyAttr(attr) 2451 return attr 2452 2453 # indices of bonds 2454 b = guess_bonds(self.atoms, self.atoms.positions, vdwradii=vdwradii) 2455 bondattr = get_TopAttr(self.universe, 'bonds', Bonds) 2456 bondattr.add_bonds(b, guessed=True) 2457 2458 a = guess_angles(self.bonds) 2459 angleattr = get_TopAttr(self.universe, 'angles', Angles) 2460 angleattr.add_bonds(a, guessed=True) 2461 2462 d = guess_dihedrals(self.angles) 2463 diheattr = get_TopAttr(self.universe, 'dihedrals', Dihedrals) 2464 diheattr.add_bonds(d) 2465 2466 @property 2467 def bond(self): 2468 """This :class:`AtomGroup` represented as a 2469 :class:`MDAnalysis.core.topologyobjects.Bond` object 2470 2471 Raises 2472 ------ 2473 ValueError 2474 If the :class:`AtomGroup` is not length 2 2475 2476 2477 .. versionadded:: 0.11.0 2478 """ 2479 if len(self) != 2: 2480 raise ValueError( 2481 "bond only makes sense for a group with exactly 2 atoms") 2482 return topologyobjects.Bond(self.ix, self.universe) 2483 2484 @property 2485 def angle(self): 2486 """This :class:`AtomGroup` represented as an 2487 :class:`MDAnalysis.core.topologyobjects.Angle` object 2488 2489 Raises 2490 ------ 2491 ValueError 2492 If the :class:`AtomGroup` is not length 3 2493 2494 2495 .. versionadded:: 0.11.0 2496 """ 2497 if len(self) != 3: 2498 raise ValueError( 2499 "angle only makes sense for a group with exactly 3 atoms") 2500 return topologyobjects.Angle(self.ix, self.universe) 2501 2502 @property 2503 def dihedral(self): 2504 """This :class:`AtomGroup` represented as a 2505 :class:`~MDAnalysis.core.topologyobjects.Dihedral` object 2506 2507 Raises 2508 ------ 2509 ValueError 2510 If the :class:`AtomGroup` is not length 4 2511 2512 2513 .. versionadded:: 0.11.0 2514 """ 2515 if len(self) != 4: 2516 raise ValueError( 2517 "dihedral only makes sense for a group with exactly 4 atoms") 2518 return topologyobjects.Dihedral(self.ix, self.universe) 2519 2520 @property 2521 def improper(self): 2522 """This :class:`AtomGroup` represented as an 2523 :class:`MDAnalysis.core.topologyobjects.ImproperDihedral` object 2524 2525 Raises 2526 ------ 2527 ValueError 2528 If the :class:`AtomGroup` is not length 4 2529 2530 2531 .. versionadded:: 0.11.0 2532 """ 2533 if len(self) != 4: 2534 raise ValueError( 2535 "improper only makes sense for a group with exactly 4 atoms") 2536 return topologyobjects.ImproperDihedral(self.ix, self.universe) 2537 2538 def write(self, filename=None, file_format="PDB", 2539 filenamefmt="{trjname}_{frame}", frames=None, **kwargs): 2540 """Write `AtomGroup` to a file. 2541 2542 The output can either be a coordinate file or a selection, depending on 2543 the format. 2544 2545 Examples 2546 -------- 2547 2548 >>> ag = u.atoms 2549 >>> ag.write('selection.ndx') # Write a gromacs index file 2550 >>> ag.write('coordinates.pdb') # Write the current frame as PDB 2551 >>> # Write the trajectory in XTC format 2552 >>> ag.write('trajectory.xtc', frames='all') 2553 >>> # Write every other frame of the trajectory in PBD format 2554 >>> ag.write('trajectory.pdb', frames=u.trajectory[::2]) 2555 2556 Parameters 2557 ---------- 2558 filename : str, optional 2559 ``None``: create TRJNAME_FRAME.FORMAT from filenamefmt [``None``] 2560 file_format : str, optional 2561 The name or extension of a coordinate, trajectory, or selection 2562 file format such as PDB, CRD, GRO, VMD (tcl), PyMol (pml), Gromacs 2563 (ndx) CHARMM (str) or Jmol (spt); case-insensitive [PDB] 2564 filenamefmt : str, optional 2565 format string for default filename; use substitution tokens 2566 'trjname' and 'frame' ["%(trjname)s_%(frame)d"] 2567 bonds : str, optional 2568 how to handle bond information, especially relevant for PDBs. 2569 ``"conect"``: write only the CONECT records defined in the original 2570 file. ``"all"``: write out all bonds, both the original defined and 2571 those guessed by MDAnalysis. ``None``: do not write out bonds. 2572 Default is ``"conect"``. 2573 frames: array-like or slice or FrameIteratorBase or str, optional 2574 An ensemble of frames to write. The ensemble can be an list or 2575 array of frame indices, a mask of booleans, an instance of 2576 :class:`slice`, or the value returned when a trajectory is indexed. 2577 By default, `frames` is set to ``None`` and only the current frame 2578 is written. If `frames` is set to "all", then all the frame from 2579 trajectory are written. 2580 2581 2582 .. versionchanged:: 0.9.0 Merged with write_selection. This method can 2583 now write both selections out. 2584 .. versionchanged:: 0.19.0 2585 Can write multiframe trajectories with the 'frames' argument. 2586 """ 2587 # TODO: Add a 'verbose' option alongside 'frames'. 2588 2589 # check that AtomGroup actually has any atoms (Issue #434) 2590 if len(self.atoms) == 0: 2591 raise IndexError("Cannot write an AtomGroup with 0 atoms") 2592 2593 trj = self.universe.trajectory # unified trajectory API 2594 if frames is None or frames == 'all': 2595 trj_frames = trj[::] 2596 elif isinstance(frames, numbers.Integral): 2597 # We accept everything that indexes a trajectory and returns a 2598 # subset of it. Though, numbers return a Timestep instead. 2599 raise TypeError('The "frames" argument cannot be a number.') 2600 else: 2601 try: 2602 test_trajectory = frames.trajectory 2603 except AttributeError: 2604 trj_frames = trj[frames] 2605 else: 2606 if test_trajectory is not trj: 2607 raise ValueError( 2608 'The trajectory of {} provided to the frames keyword ' 2609 'attribute is different from the trajectory of the ' 2610 'AtomGroup.'.format(frames) 2611 ) 2612 trj_frames = frames 2613 2614 if filename is None: 2615 trjname, ext = os.path.splitext(os.path.basename(trj.filename)) 2616 filename = filenamefmt.format(trjname=trjname, frame=trj.frame) 2617 filename = util.filename(filename, ext=file_format.lower(), keep=True) 2618 2619 # Some writer behave differently when they are given a "multiframe" 2620 # argument. It is the case of the PDB writer tht writes models when 2621 # "multiframe" is True. 2622 # We want to honor what the user provided with the argument if 2623 # provided explicitly. If not, then we need to figure out if we write 2624 # multiple frames or not. 2625 multiframe = kwargs.pop('multiframe', None) 2626 if len(trj_frames) > 1 and multiframe == False: 2627 raise ValueError( 2628 'Cannot explicitely set "multiframe" to False and request ' 2629 'more than 1 frame with the "frames" keyword argument.' 2630 ) 2631 elif multiframe is None: 2632 if frames is None: 2633 # By default we only write the current frame. 2634 multiframe = False 2635 else: 2636 multiframe = len(trj_frames) > 1 2637 2638 # From the following blocks, one must pass. 2639 # Both can't pass as the extensions don't overlap. 2640 # Try and select a Class using get_ methods (becomes `writer`) 2641 # Once (and if!) class is selected, use it in with block 2642 try: 2643 # format keyword works differently in get_writer and get_selection_writer 2644 # here it overrides everything, in get_sel it is just a default 2645 # apply sparingly here! 2646 format = os.path.splitext(filename)[1][1:] # strip initial dot! 2647 format = format or file_format 2648 format = format.strip().upper() 2649 2650 writer = get_writer_for(filename, format=format, multiframe=multiframe) 2651 except (ValueError, TypeError): 2652 pass 2653 else: 2654 with writer(filename, n_atoms=self.n_atoms, **kwargs) as w: 2655 if frames is None: 2656 w.write(self.atoms) 2657 else: 2658 current_frame = trj.ts.frame 2659 try: 2660 for _ in trj_frames: 2661 w.write(self.atoms) 2662 finally: 2663 trj[current_frame] 2664 return 2665 2666 try: 2667 # here `file_format` is only used as default, 2668 # anything pulled off `filename` will be used preferentially 2669 writer = get_selection_writer_for(filename, file_format) 2670 except (TypeError, NotImplementedError): 2671 pass 2672 else: 2673 with writer(filename, n_atoms=self.n_atoms, **kwargs) as w: 2674 w.write(self.atoms) 2675 return 2676 2677 raise ValueError("No writer found for format: {}".format(filename)) 2678 2679 2680class ResidueGroup(GroupBase): 2681 """ResidueGroup base class. 2682 2683 This class is used by a :class:`~MDAnalysis.core.universe.Universe` for 2684 generating its Topology-specific :class:`ResidueGroup` class. All the 2685 :class:`~MDAnalysis.core.topologyattrs.TopologyAttr` components are obtained 2686 from :class:`GroupBase`, so this class only includes ad-hoc methods 2687 specific to :class:`ResidueGroups<ResidueGroup>`. 2688 2689 ResidueGroups can be compared and combined using group operators. See the 2690 list of these operators on :class:`GroupBase`. 2691 2692 .. deprecated:: 0.16.2 2693 *Instant selectors* of Segments will be removed in the 1.0 release. 2694 See :ref:`Instant selectors <instance-selectors>` for details and 2695 alternatives. 2696 """ 2697 2698 @property 2699 def atoms(self): 2700 """An :class:`AtomGroup` of :class:`Atoms<Atom>` present in this 2701 :class:`ResidueGroup`. 2702 2703 The :class:`Atoms<Atom>` are ordered locally by :class:`Residue` in the 2704 :class:`ResidueGroup`. Duplicates are *not* removed. 2705 """ 2706 # If indices is an empty list np.concatenate will fail (Issue #1999). 2707 try: 2708 ag = self.universe.atoms[np.concatenate(self.indices)] 2709 except ValueError: 2710 ag = self.universe.atoms[self.indices] 2711 # If the ResidueGroup is known to be unique, this also holds for the 2712 # atoms therein, since atoms can only belong to one residue at a time. 2713 # On the contrary, if the ResidueGroup is not unique, this does not 2714 # imply non-unique atoms, since residues might be empty. 2715 try: 2716 if self._cache['isunique']: 2717 ag._cache['isunique'] = True 2718 ag._cache['unique'] = ag 2719 except KeyError: 2720 pass 2721 return ag 2722 2723 @property 2724 def n_atoms(self): 2725 """Number of :class:`Atoms<Atom>` present in this :class:`ResidueGroup`, 2726 including duplicate residues (and thus, duplicate atoms). 2727 2728 Equivalent to ``len(self.atoms)``. 2729 """ 2730 return len(self.atoms) 2731 2732 @property 2733 def residues(self): 2734 """The :class:`ResidueGroup` itself. 2735 2736 See Also 2737 -------- 2738 copy : return a true copy of the :class:`ResidueGroup` 2739 2740 2741 .. versionchanged:: 0.19.0 2742 In previous versions, this returned a copy, but now 2743 the :class:`ResidueGroup` itself is returned. This should 2744 not affect any code but only speed up calculations. 2745 2746 """ 2747 return self 2748 2749 @property 2750 def n_residues(self): 2751 """Number of residues in the :class:`ResidueGroup`. 2752 2753 Equivalent to ``len(self)``. 2754 """ 2755 return len(self) 2756 2757 @property 2758 def segments(self): 2759 """Get sorted :class:`SegmentGroup` of the unique segments present in 2760 the :class:`ResidueGroup`. 2761 """ 2762 sg = self.universe.segments[unique_int_1d(self.segindices)] 2763 sg._cache['isunique'] = True 2764 sg._cache['unique'] = sg 2765 return sg 2766 2767 @segments.setter 2768 def segments(self, new): 2769 # Can set with Seg, SegGroup or list/tuple of Seg 2770 if isinstance(new, Segment): 2771 s_ix = itertools.cycle((new.segindex,)) 2772 elif isinstance(new, SegmentGroup): 2773 s_ix = new.segindices 2774 else: 2775 try: 2776 s_ix = [s.segindex for s in new] 2777 except AttributeError: 2778 raise TypeError("Can only set ResidueGroup segments to Segment " 2779 "or SegmentGroup, not {}".format( 2780 ', '.join(type(r) for r in new 2781 if not isinstance(r, Segment)) 2782 )) 2783 if not isinstance(s_ix, itertools.cycle) and len(s_ix) != len(self): 2784 raise ValueError("Incorrect size: {} for ResidueGroup of size: {}" 2785 "".format(len(new), len(self))) 2786 # Optimisation TODO: 2787 # This currently rebuilds the tt len(self) times 2788 # Ideally all changes would happen and *afterwards* tables are built 2789 # Alternatively, if the changes didn't rebuild table, this list 2790 # comprehension isn't terrible. 2791 for r, s in zip(self, s_ix): 2792 self.universe._topology.tt.move_residue(r.ix, s) 2793 2794 @property 2795 def n_segments(self): 2796 """Number of unique segments present in the ResidueGroup. 2797 2798 Equivalent to ``len(self.segments)``. 2799 """ 2800 return len(self.segments) 2801 2802 @property 2803 @cached('unique') 2804 def unique(self): 2805 """Return a :class:`ResidueGroup` containing sorted and unique 2806 :class:`Residues<Residue>` only. 2807 2808 If the :class:`ResidueGroup` is unique, this is the group itself. 2809 2810 Examples 2811 -------- 2812 2813 >>> rg = u.residues[[2, 1, 2, 2, 1, 0]] 2814 >>> rg 2815 <ResidueGroup with 6 residues> 2816 >>> rg.ix 2817 array([2, 1, 2, 2, 1, 0]) 2818 >>> rg2 = rg.unique 2819 >>> rg2 2820 <ResidueGroup with 3 residues> 2821 >>> rg2.ix 2822 array([0, 1, 2]) 2823 >>> rg2.unique is rg2 2824 True 2825 2826 2827 .. versionadded:: 0.16.0 2828 .. versionchanged:: 0.19.0 If the :class:`ResidueGroup` is already 2829 unique, :attr:`ResidueGroup.unique` now returns the group itself 2830 instead of a copy. 2831 """ 2832 if self.isunique: 2833 return self 2834 _unique = self.universe.residues[unique_int_1d(self.ix)] 2835 # Since we know that _unique is a unique ResidueGroup, we set its 2836 # uniqueness caches from here: 2837 _unique._cache['isunique'] = True 2838 _unique._cache['unique'] = _unique 2839 return _unique 2840 2841 2842class SegmentGroup(GroupBase): 2843 """:class:`SegmentGroup` base class. 2844 2845 This class is used by a :class:`~MDAnalysis.core.universe.Universe` for 2846 generating its Topology-specific :class:`SegmentGroup` class. All the 2847 :class:`~MDAnalysis.core.topologyattrs.TopologyAttr` components are obtained 2848 from :class:`GroupBase`, so this class only includes ad-hoc methods specific 2849 to :class:`SegmentGroups<SegmentGroup>`. 2850 2851 :class:`SegmentGroups<SegmentGroup>` can be compared and combined using 2852 group operators. See the list of these operators on :class:`GroupBase`. 2853 2854 .. deprecated:: 0.16.2 2855 *Instant selectors* of Segments will be removed in the 1.0 release. 2856 See :ref:`Instant selectors <instance-selectors>` for details and 2857 alternatives. 2858 """ 2859 2860 @property 2861 def atoms(self): 2862 """An :class:`AtomGroup` of :class:`Atoms<Atom>` present in this 2863 :class:`SegmentGroup`. 2864 2865 The :class:`Atoms<Atom>` are ordered locally by :class:`Residue`, which 2866 are further ordered by :class:`Segment` in the :class:`SegmentGroup`. 2867 Duplicates are *not* removed. 2868 """ 2869 # If indices is an empty list np.concatenate will fail (Issue #1999). 2870 try: 2871 ag = self.universe.atoms[np.concatenate(self.indices)] 2872 except ValueError: 2873 ag = self.universe.atoms[self.indices] 2874 # If the SegmentGroup is known to be unique, this also holds for the 2875 # residues therein, and thus, also for the atoms in those residues. 2876 # On the contrary, if the SegmentGroup is not unique, this does not 2877 # imply non-unique atoms, since segments or residues might be empty. 2878 try: 2879 if self._cache['isunique']: 2880 ag._cache['isunique'] = True 2881 ag._cache['unique'] = ag 2882 except KeyError: 2883 pass 2884 return ag 2885 2886 @property 2887 def n_atoms(self): 2888 """Number of atoms present in the :class:`SegmentGroup`, including 2889 duplicate segments (and thus, duplicate atoms). 2890 2891 Equivalent to ``len(self.atoms)``. 2892 """ 2893 return len(self.atoms) 2894 2895 @property 2896 def residues(self): 2897 """A :class:`ResidueGroup` of :class:`Residues<Residue>` present in this 2898 :class:`SegmentGroup`. 2899 2900 The :class:`Residues<Residue>` are ordered locally by 2901 :class:`Segment` in the :class:`SegmentGroup`. Duplicates are *not* 2902 removed. 2903 """ 2904 rg = self.universe.residues[np.concatenate(self.resindices)] 2905 # If the SegmentGroup is known to be unique, this also holds for the 2906 # residues therein. On the contrary, if the SegmentGroup is not unique, 2907 # this does not imply non-unique residues, since segments might be 2908 # empty. 2909 try: 2910 if self._cache['isunique']: 2911 rg._cache['isunique'] = True 2912 rg._cache['unique'] = rg 2913 except KeyError: 2914 pass 2915 return rg 2916 2917 @property 2918 def n_residues(self): 2919 """Number of residues present in this :class:`SegmentGroup`, including 2920 duplicate segments (and thus, residues). 2921 2922 Equivalent to ``len(self.residues)``. 2923 """ 2924 return len(self.residues) 2925 2926 @property 2927 def segments(self): 2928 """The :class:`SegmentGroup` itself. 2929 2930 See Also 2931 -------- 2932 copy : return a true copy of the :class:`SegmentGroup` 2933 2934 2935 .. versionchanged:: 0.19.0 2936 In previous versions, this returned a copy, but now 2937 the :class:`SegmentGroup` itself is returned. This should 2938 not affect any code but only speed up calculations. 2939 2940 """ 2941 return self 2942 2943 @property 2944 def n_segments(self): 2945 """Number of segments in the :class:`SegmentGroup`. 2946 2947 Equivalent to ``len(self)``. 2948 """ 2949 return len(self) 2950 2951 @property 2952 @cached('unique') 2953 def unique(self): 2954 """Return a :class:`SegmentGroup` containing sorted and unique 2955 :class:`Segments<Segment>` only. 2956 2957 If the :class:`SegmentGroup` is unique, this is the group itself. 2958 2959 Examples 2960 -------- 2961 2962 >>> sg = u.segments[[2, 1, 2, 2, 1, 0]] 2963 >>> sg 2964 <SegmentGroup with 6 segments> 2965 >>> sg.ix 2966 array([2, 1, 2, 2, 1, 0]) 2967 >>> sg2 = sg.unique 2968 >>> sg2 2969 <SegmentGroup with 3 segments> 2970 >>> sg2.ix 2971 array([0, 1, 2]) 2972 >>> sg2.unique is sg2 2973 True 2974 2975 2976 .. versionadded:: 0.16.0 2977 .. versionchanged:: 0.19.0 If the :class:`SegmentGroup` is already 2978 unique, :attr:`SegmentGroup.unique` now returns the group itself 2979 instead of a copy. 2980 """ 2981 if self.isunique: 2982 return self 2983 _unique = self.universe.segments[unique_int_1d(self.ix)] 2984 # Since we know that _unique is a unique SegmentGroup, we set its 2985 # uniqueness caches from here: 2986 _unique._cache['isunique'] = True 2987 _unique._cache['unique'] = _unique 2988 return _unique 2989 2990 2991@functools.total_ordering 2992class ComponentBase(_MutableBase): 2993 """Base class from which a :class:`~MDAnalysis.core.universe.Universe`\ 's 2994 Component class is built. 2995 2996 Components are the individual objects that are found in Groups. 2997 """ 2998 def __init__(self, ix, u): 2999 # index of component 3000 self._ix = ix 3001 self._u = u 3002 3003 def __lt__(self, other): 3004 if self.level != other.level: 3005 raise TypeError("Can't compare different level objects") 3006 return self.ix < other.ix 3007 3008 def __eq__(self, other): 3009 if self.level != other.level: 3010 raise TypeError("Can't compare different level objects") 3011 return self.ix == other.ix 3012 3013 def __ne__(self, other): 3014 return not self == other 3015 3016 def __hash__(self): 3017 return hash(self.ix) 3018 3019 @_only_same_level 3020 def __add__(self, other): 3021 """Concatenate the Component with another Component or Group of the 3022 same level. 3023 3024 Parameters 3025 ---------- 3026 other : Component or Group 3027 Component or Group with `other.level` same as `self.level` 3028 3029 Returns 3030 ------- 3031 Group 3032 Group with elements of `self` and `other` concatenated 3033 """ 3034 o_ix = other.ix_array 3035 3036 return self.level.plural( 3037 np.concatenate([self.ix_array, o_ix]), self.universe) 3038 3039 def __radd__(self, other): 3040 """Using built-in sum requires supporting 0 + self. If other is 3041 anything other 0, an exception will be raised. 3042 3043 Parameters 3044 ---------- 3045 other : int 3046 Other should be 0, or else an exception will be raised. 3047 3048 Returns 3049 ------- 3050 self 3051 Group with elements of `self` reproduced 3052 """ 3053 if other == 0: 3054 return self.level.plural(self.ix_array, self.universe) 3055 else: 3056 raise TypeError("unsupported operand type(s) for +:" 3057 " '{}' and '{}'".format(type(self).__name__, 3058 type(other).__name__)) 3059 3060 @property 3061 def universe(self): 3062 return self._u 3063 3064 @property 3065 def ix(self): 3066 """Unique index of this component. 3067 3068 If this component is an :class:`Atom`, this is the index of the 3069 :class:`Atom`. 3070 If it is a :class:`Residue`, this is the index of the :class:`Residue`. 3071 If it is a :class:`Segment`, this is the index of the :class:`Segment`. 3072 """ 3073 return self._ix 3074 3075 @property 3076 def ix_array(self): 3077 """Unique index of this component as an array. 3078 3079 This method gives a consistent API between components and groups. 3080 3081 See Also 3082 -------- 3083 ix 3084 """ 3085 return np.array([self.ix], dtype=np.intp) 3086 3087 3088class Atom(ComponentBase): 3089 """:class:`Atom` base class. 3090 3091 This class is used by a :class:`~MDAnalysis.core.universe.Universe` for 3092 generating its Topology-specific :class:`Atom` class. All the 3093 :class:`~MDAnalysis.core.topologyattrs.TopologyAttr` components are obtained 3094 from :class:`ComponentBase`, so this class only includes ad-hoc methods 3095 specific to :class:`Atoms<Atom>`. 3096 """ 3097 def __getattr__(self, attr): 3098 """Try and catch known attributes and give better error message""" 3099 if attr in ('fragment',): 3100 raise NoDataError("Atom has no fragment data, this requires Bonds") 3101 else: 3102 raise AttributeError("{cls} has no attribute {attr}".format( 3103 cls=self.__class__.__name__, attr=attr)) 3104 3105 def __repr__(self): 3106 me = '<Atom {}:'.format(self.ix + 1) 3107 if hasattr(self, 'name'): 3108 me += ' {}'.format(self.name) 3109 if hasattr(self, 'type'): 3110 me += ' of type {}'.format(self.type) 3111 if hasattr(self, 'resname'): 3112 me += ' of resname {},'.format(self.resname) 3113 if hasattr(self, 'resid'): 3114 me += ' resid {}'.format(self.resid) 3115 if hasattr(self, 'segid'): 3116 me += ' and segid {}'.format(self.segid) 3117 if hasattr(self, 'altLoc'): 3118 me += ' and altLoc {}'.format(self.altLoc) 3119 return me + '>' 3120 3121 @property 3122 def residue(self): 3123 return self.universe.residues[self.universe._topology.resindices[self]] 3124 3125 @residue.setter 3126 def residue(self, new): 3127 if not isinstance(new, Residue): 3128 raise TypeError("Can only set Atom residue to Residue, not {}" 3129 "".format(type(new))) 3130 self.universe._topology.tt.move_atom(self.ix, new.resindex) 3131 3132 @property 3133 def segment(self): 3134 return self.universe.segments[self.universe._topology.segindices[self]] 3135 3136 @segment.setter 3137 def segment(self, new): 3138 raise NotImplementedError("Cannot set atom segment. " 3139 "Segments are assigned to Residues") 3140 3141 @property 3142 def position(self): 3143 """Coordinates of the atom. 3144 3145 The position can be changed by assigning an array of length (3,). 3146 3147 .. note:: changing the position is not reflected in any files; reading 3148 any frame from the trajectory will replace the change with 3149 that from the file 3150 3151 Raises 3152 ------ 3153 ~MDAnalysis.exceptions.NoDataError 3154 If the underlying :class:`~MDAnalysis.coordinates.base.Timestep` 3155 does not contain 3156 :attr:`~MDAnalysis.coordinates.base.Timestep.positions`. 3157 """ 3158 return self.universe.trajectory.ts.positions[self.ix].copy() 3159 3160 @position.setter 3161 def position(self, values): 3162 self.universe.trajectory.ts.positions[self.ix, :] = values 3163 3164 @property 3165 def velocity(self): 3166 """Velocity of the atom. 3167 3168 The velocity can be changed by assigning an array of shape ``(3,)``. 3169 3170 .. note:: changing the velocity is not reflected in any files; reading 3171 any frame from the trajectory will replace the change with 3172 that from the file 3173 3174 Raises 3175 ------ 3176 ~MDAnalysis.exceptions.NoDataError 3177 If the underlying :class:`~MDAnalysis.coordinates.base.Timestep` 3178 does not contain 3179 :attr:`~MDAnalysis.coordinates.base.Timestep.velocities`. 3180 """ 3181 ts = self.universe.trajectory.ts 3182 try: 3183 return ts.velocities[self.ix].copy() 3184 except (AttributeError, NoDataError): 3185 raise NoDataError("Timestep does not contain velocities") 3186 3187 @velocity.setter 3188 def velocity(self, values): 3189 ts = self.universe.trajectory.ts 3190 try: 3191 ts.velocities[self.ix, :] = values 3192 except (AttributeError, NoDataError): 3193 raise NoDataError("Timestep does not contain velocities") 3194 3195 @property 3196 def force(self): 3197 """Force on the atom. 3198 3199 The force can be changed by assigning an array of shape ``(3,)``. 3200 3201 .. note:: changing the force is not reflected in any files; reading any 3202 frame from the trajectory will replace the change with that 3203 from the file 3204 3205 Raises 3206 ------ 3207 ~MDAnalysis.exceptions.NoDataError 3208 If the underlying :class:`~MDAnalysis.coordinates.base.Timestep` 3209 does not contain 3210 :attr:`~MDAnalysis.coordinates.base.Timestep.forces`. 3211 """ 3212 ts = self.universe.trajectory.ts 3213 try: 3214 return ts.forces[self.ix].copy() 3215 except (AttributeError, NoDataError): 3216 raise NoDataError("Timestep does not contain forces") 3217 3218 @force.setter 3219 def force(self, values): 3220 ts = self.universe.trajectory.ts 3221 try: 3222 ts.forces[self.ix, :] = values 3223 except (AttributeError, NoDataError): 3224 raise NoDataError("Timestep does not contain forces") 3225 3226 3227class Residue(ComponentBase): 3228 """:class:`Residue` base class. 3229 3230 This class is used by a :class:`~MDAnalysis.core.universe.Universe` for 3231 generating its Topology-specific :class:`Residue` class. All the 3232 :class:`~MDAnalysis.core.topologyattrs.TopologyAttr` components are obtained 3233 from :class:`ComponentBase`, so this class only includes ad-hoc methods 3234 specific to :class:`Residues<Residue>`. 3235 """ 3236 def __repr__(self): 3237 me = '<Residue' 3238 if hasattr(self, 'resname'): 3239 me += ' {},'.format(self.resname) 3240 if hasattr(self, 'resid'): 3241 me += ' {}'.format(self.resid) 3242 3243 return me + '>' 3244 3245 @property 3246 def atoms(self): 3247 """An :class:`AtomGroup` of :class:`Atoms<Atom>` present in this 3248 :class:`Residue`. 3249 """ 3250 ag = self.universe.atoms[self.universe._topology.indices[self][0]] 3251 ag._cache['isunique'] = True 3252 ag._cache['unique'] = ag 3253 return ag 3254 3255 @property 3256 def segment(self): 3257 """The :class:`Segment` this :class:`Residue` belongs to. 3258 """ 3259 return self.universe.segments[self.universe._topology.segindices[self]] 3260 3261 @segment.setter 3262 def segment(self, new): 3263 if not isinstance(new, Segment): 3264 raise TypeError("Can only set Residue segment to Segment, not {}" 3265 "".format(type(new))) 3266 self.universe._topology.tt.move_residue(self.ix, new.segindex) 3267 3268 3269class Segment(ComponentBase): 3270 """:class:`Segment` base class. 3271 3272 This class is used by a :class:`~MDAnalysis.core.universe.Universe` for 3273 generating its Topology-specific :class:`Segment` class. All the 3274 :class:`~MDAnalysis.core.topologyattrs.TopologyAttr` components are obtained 3275 from :class:`ComponentBase`, so this class only includes ad-hoc methods 3276 specific to :class:`Segments<Segment>`. 3277 3278 .. deprecated:: 0.16.2 3279 *Instant selectors* of :class:`Segments<Segment>` will be removed in the 3280 1.0 release. See :ref:`Instant selectors <instance-selectors>` for 3281 details and alternatives. 3282 """ 3283 def __repr__(self): 3284 me = '<Segment' 3285 if hasattr(self, 'segid'): 3286 me += ' {}'.format(self.segid) 3287 return me + '>' 3288 3289 @property 3290 def atoms(self): 3291 """An :class:`AtomGroup` of :class:`Atoms<Atom>` present in this 3292 :class:`Segment`. 3293 """ 3294 ag = self.universe.atoms[self.universe._topology.indices[self][0]] 3295 ag._cache['isunique'] = True 3296 ag._cache['unique'] = ag 3297 return ag 3298 3299 @property 3300 def residues(self): 3301 """A :class:`ResidueGroup` of :class:`Residues<Residue>` present in this 3302 :class:`Segment`. 3303 """ 3304 rg = self.universe.residues[self.universe._topology.resindices[self][0]] 3305 rg._cache['isunique'] = True 3306 rg._cache['unique'] = rg 3307 return rg 3308 3309 def __getattr__(self, attr): 3310 # DEPRECATED in 0.16.2 3311 # REMOVE in 1.0 3312 # 3313 # Segment.r1 access 3314 if attr.startswith('r') and attr[1:].isdigit(): 3315 resnum = int(attr[1:]) 3316 rg = self.residues[resnum - 1] # convert to 0 based 3317 warnings.warn("Instant selectors Segment.r<N> will be removed in " 3318 "1.0. Use Segment.residues[N-1] instead.", 3319 DeprecationWarning) 3320 return rg 3321 # Resname accesss 3322 if hasattr(self.residues, 'resnames'): 3323 try: 3324 return self.residues._get_named_residue(attr) 3325 except selection.SelectionError: 3326 pass 3327 raise AttributeError("{cls} has no attribute {attr}" 3328 "".format(cls=self.__class__.__name__, attr=attr)) 3329 3330# Accessing these attrs doesn't trigger an update. The class and instance 3331# methods of UpdatingAtomGroup that are used during __init__ must all be 3332# here, otherwise we get __getattribute__ infinite loops. 3333_UAG_SHORTCUT_ATTRS = { 3334 # Class information of the UAG 3335 "__class__", "_derived_class", 3336 # Metadata of the UAG 3337 "_base_group", "_selections", "_lastupdate", 3338 "level", "_u", "universe", 3339 # Methods of the UAG 3340 "_ensure_updated", 3341 "is_uptodate", 3342 "update_selection", 3343} 3344 3345class UpdatingAtomGroup(AtomGroup): 3346 """:class:`AtomGroup` subclass that dynamically updates its selected atoms. 3347 3348 Accessing any attribute/method of an :class:`UpdatingAtomGroup` instance 3349 triggers a check for the last frame the group was updated. If the last 3350 frame matches the current trajectory frame, the attribute is returned 3351 normally; otherwise the group is updated (the stored selections are 3352 re-applied), and only then is the attribute returned. 3353 3354 3355 .. versionadded:: 0.16.0 3356 """ 3357 # WARNING: This class has __getattribute__ and __getattr__ methods (the 3358 # latter inherited from AtomGroup). Because of this bugs introduced in the 3359 # class that cause an AttributeError may be very hard to diagnose and 3360 # debug: the most obvious symptom is an infinite loop going through both 3361 # __getattribute__ and __getattr__, and a solution might be to add said 3362 # attribute to _UAG_SHORTCUT_ATTRS. 3363 3364 def __init__(self, base_group, selections, strings): 3365 """ 3366 3367 Parameters 3368 ---------- 3369 base_group : :class:`AtomGroup` 3370 group on which *selections* are to be applied. 3371 selections : a tuple of :class:`~MDAnalysis.core.selection.Selection` 3372 instances selections ready to be applied to *base_group*. 3373 """ 3374 # Because we're implementing __getattribute__, which needs _u for 3375 # its check, no self.attribute access can be made before this line 3376 self._u = base_group.universe 3377 self._selections = selections 3378 self._selection_strings = strings 3379 self._base_group = base_group 3380 self._lastupdate = None 3381 self._derived_class = base_group._derived_class 3382 if self._selections: 3383 # Allows the creation of a cheap placeholder UpdatingAtomGroup 3384 # by passing an empty selection tuple. 3385 self._ensure_updated() 3386 3387 def update_selection(self): 3388 """ 3389 Forces the reevaluation and application of the group's selection(s). 3390 3391 This method is triggered automatically when accessing attributes, if 3392 the last update occurred under a different trajectory frame. 3393 """ 3394 bg = self._base_group 3395 sels = self._selections 3396 if sels: 3397 # As with select_atoms, we select the first sel and then sum to it. 3398 ix = sum([sel.apply(bg) for sel in sels[1:]], 3399 sels[0].apply(bg)).ix 3400 else: 3401 ix = np.array([], dtype=np.intp) 3402 # Run back through AtomGroup init with this information to remake 3403 # ourselves 3404 super(UpdatingAtomGroup, self).__init__(ix, self.universe) 3405 self.is_uptodate = True 3406 3407 @property 3408 def is_uptodate(self): 3409 """ 3410 Checks whether the selection needs updating based on frame number only. 3411 3412 Modifications to the coordinate data that render selections stale are 3413 not caught, and in those cases :attr:`is_uptodate` may return an 3414 erroneous value. 3415 3416 Returns 3417 ------- 3418 bool 3419 ``True`` if the group's selection is up-to-date, ``False`` 3420 otherwise. 3421 """ 3422 try: 3423 return self.universe.trajectory.frame == self._lastupdate 3424 except AttributeError: # self.universe has no trajectory 3425 return self._lastupdate == -1 3426 3427 @is_uptodate.setter 3428 def is_uptodate(self, value): 3429 if value: 3430 try: 3431 self._lastupdate = self.universe.trajectory.frame 3432 except AttributeError: # self.universe has no trajectory 3433 self._lastupdate = -1 3434 else: 3435 # This always marks the selection as outdated 3436 self._lastupdate = None 3437 3438 def _ensure_updated(self): 3439 """ 3440 Checks whether the selection needs updating and updates it if needed. 3441 3442 Returns 3443 ------- 3444 bool 3445 ``True`` if the group was already up-to-date, ``False`` otherwise. 3446 """ 3447 status = self.is_uptodate 3448 if not status: 3449 self.update_selection() 3450 return status 3451 3452 def __getattribute__(self, name): 3453 # ALL attribute access goes through here 3454 # If the requested attribute is public (not starting with '_') and 3455 # isn't in the shortcut list, update ourselves 3456 if not (name.startswith('_') or name in _UAG_SHORTCUT_ATTRS): 3457 self._ensure_updated() 3458 # Going via object.__getattribute__ then bypasses this check stage 3459 return object.__getattribute__(self, name) 3460 3461 def __reduce__(self): 3462 # strategy for unpickling is: 3463 # - unpickle base group 3464 # - recreate UAG as created through select_atoms (basegroup and selstrs) 3465 # even if base_group is a UAG this will work through recursion 3466 return (_unpickle_uag, 3467 (self._base_group.__reduce__(), self._selections, 3468 self._selection_strings)) 3469 3470 def __repr__(self): 3471 basestr = super(UpdatingAtomGroup, self).__repr__() 3472 if not self._selection_strings: 3473 return basestr 3474 sels = "'{}'".format("' + '".join(self._selection_strings)) 3475 # Cheap comparison. Might fail for corner cases but this is 3476 # mostly cosmetic. 3477 if self._base_group is self.universe.atoms: 3478 basegrp = "the entire Universe." 3479 else: 3480 basegrp = "another AtomGroup." 3481 # With a shorthand to conditionally append the 's' in 'selections'. 3482 return "{}, with selection{} {} on {}>".format(basestr[:-1], 3483 "s"[len(self._selection_strings)==1:], sels, basegrp) 3484 3485 @property 3486 def atoms(self): 3487 """Get a *static* :class:`AtomGroup` identical to the group of currently 3488 selected :class:`Atoms<Atom>` in the :class:`UpdatingAtomGroup`. 3489 3490 3491 By returning a *static* :class:`AtomGroup` it becomes possible to 3492 compare the contents of the group *between* trajectory frames. See the 3493 Example below. 3494 3495 3496 Note 3497 ---- 3498 The :attr:`atoms` attribute of an :class:`UpdatingAtomGroup` behaves 3499 differently from :attr:`AtomGroup.atoms`: the latter returns the 3500 :class:`AtomGroup` itself whereas the former returns a 3501 :class:`AtomGroup` and not an :class:`UpdatingAtomGroup` (for this, use 3502 :meth:`UpdatingAtomGroup.copy`). 3503 3504 3505 Example 3506 ------- 3507 The static :attr:`atoms` allows comparison of groups of atoms between 3508 frames. For example, track water molecules that move in and out of a 3509 solvation shell of a protein:: 3510 3511 u = mda.Universe(TPR, XTC) 3512 water_shell = u.select_atoms("name OW and around 3.5 protein", updating=True) 3513 water_shell_prev = water_shell.atoms 3514 3515 for ts in u.trajectory: 3516 exchanged = water_shell - water_shell_prev 3517 3518 print(ts.time, "waters in shell =", water_shell.n_residues) 3519 print(ts.time, "waters that exchanged = ", exchanged.n_residues) 3520 print(ts.time, "waters that remained bound = ", 3521 water_shell.n_residues - exchanged.n_residues) 3522 3523 water_shell_prev = water_shell.atoms 3524 3525 By remembering the atoms of the current time step in 3526 `water_shell_prev`, it becomes possible to use the :meth:`subtraction 3527 of AtomGroups<AtomGroup.subtract>` to find the water molecules that 3528 changed. 3529 3530 3531 See Also 3532 -------- 3533 copy : return a true copy of the :class:`UpdatingAtomGroup` 3534 3535 """ 3536 return self[:] 3537 3538 def copy(self): 3539 """Get another :class:`UpdatingAtomGroup` identical to this one. 3540 3541 3542 .. versionadded:: 0.19.0 3543 """ 3544 return UpdatingAtomGroup(self._base_group, self._selections, 3545 self._selection_strings) 3546 3547 3548# Define relationships between these classes 3549# with Level objects 3550_Level = namedtuple('Level', ['name', 'singular', 'plural']) 3551ATOMLEVEL = _Level('atom', Atom, AtomGroup) 3552RESIDUELEVEL = _Level('residue', Residue, ResidueGroup) 3553SEGMENTLEVEL = _Level('segment', Segment, SegmentGroup) 3554 3555Atom.level = ATOMLEVEL 3556AtomGroup.level = ATOMLEVEL 3557Residue.level = RESIDUELEVEL 3558ResidueGroup.level = RESIDUELEVEL 3559Segment.level = SEGMENTLEVEL 3560SegmentGroup.level = SEGMENTLEVEL 3561 3562def requires(*attrs): 3563 """Decorator to check if all :class:`AtomGroup` arguments have certain 3564 attributes 3565 3566 Example 3567 ------- 3568 When used to wrap a function, will check all :class:`AtomGroup` arguments 3569 for the listed requirements 3570 3571 @requires('masses', 'charges') 3572 def mass_times_charge(atomgroup): 3573 return atomgroup.masses * atomgroup.charges 3574 """ 3575 def require_dec(func): 3576 @functools.wraps(func) 3577 def check_args(*args, **kwargs): 3578 for a in args: # for each argument 3579 if isinstance(a, AtomGroup): 3580 # Make list of missing attributes 3581 missing = [attr for attr in attrs 3582 if not hasattr(a, attr)] 3583 if missing: 3584 raise NoDataError( 3585 "{funcname} failed. " 3586 "AtomGroup is missing the following required " 3587 "attributes: {attrs}".format( 3588 funcname=func.__name__, 3589 attrs=', '.join(missing))) 3590 return func(*args, **kwargs) 3591 return check_args 3592 return require_dec 3593