1# -*- Mode: python; tab-width: 4; indent-tabs-mode:nil; coding: utf-8 -*- 2# vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4 3# 4# MDAnalysis --- https://www.mdanalysis.org 5# Copyright (c) 2006-2017 The MDAnalysis Development Team and contributors 6# (see the file AUTHORS for the full list of names) 7# 8# Released under the GNU Public Licence, v2 or any higher version 9# 10# Please cite your use of MDAnalysis in published work: 11# 12# R. J. Gowers, M. Linke, J. Barnoud, T. J. E. Reddy, M. N. Melo, S. L. Seyler, 13# D. L. Dotson, J. Domanski, S. Buchoux, I. M. Kenney, and O. Beckstein. 14# MDAnalysis: A Python package for the rapid analysis of molecular dynamics 15# simulations. In S. Benthall and S. Rostrup editors, Proceedings of the 15th 16# Python in Science Conference, pages 102-109, Austin, TX, 2016. SciPy. 17# 18# N. Michaud-Agrawal, E. J. Denning, T. B. Woolf, and O. Beckstein. 19# MDAnalysis: A Toolkit for the Analysis of Molecular Dynamics Simulations. 20# J. Comput. Chem. 32 (2011), 2319--2327, doi:10.1002/jcc.21787 21# 22 23"""\ 24========================================================= 25Core object: Universe --- :mod:`MDAnalysis.core.universe` 26========================================================= 27 28The :class:`~MDAnalysis.core.universe.Universe` class ties a topology 29and a trajectory together. Almost all code in MDAnalysis starts with a 30``Universe``. 31 32Normally, a ``Universe`` is created from files:: 33 34 import MDAnalysis as mda 35 u = mda.Universe("topology.psf", "trajectory.dcd") 36 37In order to construct new simulation system it is also convenient to 38construct a ``Universe`` from existing 39:class:`~MDAnalysis.core.group.AtomGroup` instances with the 40:func:`Merge` function. 41 42 43Working with Universes 44====================== 45 46 47Quick segid selection 48--------------------- 49 50.. deprecated:: 0.16.2 51 Instant selectors will be removed in the 1.0 release. See issue `#1377 52 <https://github.com/MDAnalysis/mdanalysis/issues/1377>`_ for more details. 53 54 55If the loaded topology provided segids, then these are made accessible 56as attributes of the Universe. If the segid starts with a number such 57as '4AKE', the letter 's' will be prepended to the segid. 58For example:: 59 60 import MDAnalysis as mda 61 from MDAnalysisTests.datafiles import PSF, DCD 62 63 u = mda.Universe(PSF, DCD) 64 u.select_atoms('segid 4AKE') # selects all segments with segid 4AKE 65 66If only a single segment has that segid then a Segment object will 67be returned, otherwise a SegmentGroup will be returned. 68 69 70Classes 71======= 72 73.. autoclass:: Universe 74 :members: 75 76Functions 77========= 78 79.. autofunction:: Merge 80 81""" 82from __future__ import absolute_import 83from six.moves import range 84import six 85 86import errno 87import numpy as np 88import logging 89import copy 90import warnings 91 92import MDAnalysis 93import sys 94 95# When used in an MPI environment with Infiniband, importing MDAnalysis may 96# trigger an MPI warning because importing the uuid module triggers a call to 97# os.fork(). This happens if MPI_Init() has been called prior to importing 98# MDAnalysis. The problem is actually caused by the uuid module and not by 99# MDAnalysis itself. Python 3.7 fixes the problem. However, for Python < 3.7, 100# the uuid module works perfectly fine with os.fork() disabled during import. 101# A clean solution is therefore simply to disable os.fork() prior to importing 102# the uuid module and to re-enable it afterwards. 103import os 104 105# Windows doesn't have os.fork so can ignore 106# the issue for that platform 107if sys.version_info >= (3, 7) or os.name == 'nt': 108 import uuid 109else: 110 _os_dot_fork, os.fork = os.fork, None 111 import uuid 112 os.fork = _os_dot_fork 113 del _os_dot_fork 114 115from .. import _ANCHOR_UNIVERSES, _TOPOLOGY_ATTRS, _PARSERS 116from ..exceptions import NoDataError 117from ..lib import util 118from ..lib.log import ProgressMeter 119from ..lib.util import cached, NamedStream, isstream 120from ..lib.mdamath import find_fragments 121from . import groups 122from ._get_readers import get_reader_for, get_parser_for 123from .groups import (ComponentBase, GroupBase, 124 Atom, Residue, Segment, 125 AtomGroup, ResidueGroup, SegmentGroup) 126from .topology import Topology 127from .topologyattrs import AtomAttr, ResidueAttr, SegmentAttr 128 129logger = logging.getLogger("MDAnalysis.core.universe") 130 131 132class Universe(object): 133 """The MDAnalysis Universe contains all the information describing the system. 134 135 The system always requires a *topology* file --- in the simplest case just 136 a list of atoms. This can be a CHARMM/NAMD PSF file or a simple coordinate 137 file with atom informations such as XYZ, PDB, Gromacs GRO, or CHARMM 138 CRD. See :ref:`Supported topology formats` for what kind of topologies can 139 be read. 140 141 A trajectory provides coordinates; the coordinates have to be ordered in 142 the same way as the list of atoms in the topology. A trajectory can be a 143 single frame such as a PDB, CRD, or GRO file, or it can be a MD trajectory 144 (in CHARMM/NAMD/LAMMPS DCD, Gromacs XTC/TRR, or generic XYZ format). See 145 :ref:`Supported coordinate formats` for what can be read as a 146 "trajectory". 147 148 As a special case, when the topology is a file that contains atom 149 information *and* coordinates (such as XYZ, PDB, GRO or CRD, see 150 :ref:`Supported coordinate formats`) then the coordinates are immediately 151 loaded from the "topology" file unless a trajectory is supplied. 152 153 Examples for setting up a universe:: 154 155 u = Universe(topology, trajectory) # read system from file(s) 156 u = Universe(pdbfile) # read atoms and coordinates from PDB or GRO 157 u = Universe(topology, [traj1, traj2, ...]) # read from a list of trajectories 158 u = Universe(topology, traj1, traj2, ...) # read from multiple trajectories 159 160 Load new data into a universe (replaces old trajectory and does *not* append):: 161 162 u.load_new(trajectory) # read from a new trajectory file 163 164 Select atoms, with syntax similar to CHARMM (see 165 :class:`~Universe.select_atoms` for details):: 166 167 u.select_atoms(...) 168 169 Parameters 170 ---------- 171 topology : str, Topology object or stream 172 A CHARMM/XPLOR PSF topology file, PDB file or Gromacs GRO file; used to 173 define the list of atoms. If the file includes bond information, 174 partial charges, atom masses, ... then these data will be available to 175 MDAnalysis. A "structure" file (PSF, PDB or GRO, in the sense of a 176 topology) is always required. Alternatively, an existing 177 :class:`MDAnalysis.core.topology.Topology` instance may also be given. 178 topology_format 179 Provide the file format of the topology file; ``None`` guesses it from 180 the file extension [``None``] Can also pass a subclass of 181 :class:`MDAnalysis.topology.base.TopologyReaderBase` to define a custom 182 reader to be used on the topology file. 183 format 184 Provide the file format of the coordinate or trajectory file; ``None`` 185 guesses it from the file extension. Note that this keyword has no 186 effect if a list of file names is supplied because the "chained" reader 187 has to guess the file format for each individual list member. 188 [``None``] Can also pass a subclass of 189 :class:`MDAnalysis.coordinates.base.ProtoReader` to define a custom 190 reader to be used on the trajectory file. 191 all_coordinates : bool 192 If set to ``True`` specifies that if more than one filename is passed 193 they are all to be used, if possible, as coordinate files (employing a 194 :class:`MDAnalysis.coordinates.chain.ChainReader`). [``False``] The 195 default behavior is to take the first file as a topology and the 196 remaining as coordinates. The first argument will always always be used 197 to infer a topology regardless of *all_coordinates*. This parameter is 198 ignored if only one argument is passed. 199 guess_bonds : bool, optional 200 Once Universe has been loaded, attempt to guess the connectivity 201 between atoms. This will populate the .bonds .angles and .dihedrals 202 attributes of the Universe. 203 vdwradii : dict, optional 204 For use with *guess_bonds*. Supply a dict giving a vdwradii for each 205 atom type which are used in guessing bonds. 206 is_anchor : bool, optional 207 When unpickling instances of 208 :class:`MDAnalysis.core.groups.AtomGroup` existing Universes are 209 searched for one where to anchor those atoms. Set to ``False`` to 210 prevent this Universe from being considered. [``True``] 211 anchor_name : str, optional 212 Setting to other than ``None`` will cause 213 :class:`MDAnalysis.core.groups.AtomGroup` instances pickled from the 214 Universe to only unpickle if a compatible Universe with matching 215 *anchor_name* is found. Even if *anchor_name* is set *is_anchor* will 216 still be honored when unpickling. 217 transformations: function or list, optional 218 Provide a list of transformations that you wish to apply to the 219 trajectory upon reading. Transformations can be found in 220 :mod:`MDAnalysis.transformations`, or can be user-created. 221 in_memory 222 After reading in the trajectory, transfer it to an in-memory 223 representations, which allow for manipulation of coordinates. 224 in_memory_step 225 Only read every nth frame into in-memory representation. 226 continuous : bool, optional 227 The `continuous` option is used by the 228 :mod:`ChainReader<MDAnalysis.coordinates.chain>`, which contains the 229 functionality to treat independent trajectory files as a single virtual 230 trajectory. 231 232 Attributes 233 ---------- 234 trajectory 235 currently loaded trajectory reader; 236 dimensions 237 current system dimensions (simulation unit cell, if set in the 238 trajectory) 239 atoms, residues, segments 240 master Groups for each topology level 241 bonds, angles, dihedrals 242 master ConnectivityGroups for each connectivity type 243 244 """ 245 246 def __init__(self, *args, **kwargs): 247 # Store the segments for the deprecated instant selector feature. 248 # This attribute has to be defined early to avoid recursion in 249 # __getattr__. 250 self._instant_selectors = {} 251 # hold on to copy of kwargs; used by external libraries that 252 # reinitialize universes 253 self._kwargs = copy.deepcopy(kwargs) 254 255 # managed attribute holding Reader 256 self._trajectory = None 257 self._cache = {} 258 259 if not args: 260 # create an empty universe 261 self._topology = None 262 self.atoms = None 263 else: 264 topology_format = kwargs.pop('topology_format', None) 265 if len(args) == 1: 266 # special hacks to treat a coordinate file as a coordinate AND 267 # topology file 268 if kwargs.get('format', None) is None: 269 kwargs['format'] = topology_format 270 elif topology_format is None: 271 topology_format = kwargs.get('format', None) 272 273 # if we're given a Topology object, we don't need to parse anything 274 if isinstance(args[0], Topology): 275 self._topology = args[0] 276 self.filename = None 277 else: 278 if isinstance(args[0], NamedStream): 279 self.filename = args[0] 280 elif isstream(args[0]): 281 filename = None 282 if hasattr(args[0], 'name'): 283 filename = args[0].name 284 self.filename = NamedStream(args[0], filename) 285 else: 286 self.filename = args[0] 287 parser = get_parser_for(self.filename, format=topology_format) 288 try: 289 with parser(self.filename) as p: 290 self._topology = p.parse(**kwargs) 291 except (IOError, OSError) as err: 292 # There are 2 kinds of errors that might be raised here: 293 # one because the file isn't present 294 # or the permissions are bad, second when the parser fails 295 if (err.errno is not None and 296 errno.errorcode[err.errno] in ['ENOENT', 'EACCES']): 297 # Runs if the error is propagated due to no permission / file not found 298 six.reraise(*sys.exc_info()) 299 else: 300 # Runs when the parser fails 301 raise IOError( 302 "Failed to load from the topology file {0}" 303 " with parser {1}.\n" 304 "Error: {2}".format(self.filename, parser, err)) 305 except (ValueError, NotImplementedError) as err: 306 raise ValueError( 307 "Failed to construct topology from file {0}" 308 " with parser {1}.\n" 309 "Error: {2}".format(self.filename, parser, err)) 310 311 # generate and populate Universe version of each class 312 self._generate_from_topology() 313 314 # Load coordinates 315 if len(args) == 1 or kwargs.get('all_coordinates', False): 316 if self.filename is None: 317 # If we got the topology as a Topology object, then we 318 # cannot read coordinates from it. 319 coordinatefile = args[1:] 320 else: 321 # Can the topology file also act as coordinate file? 322 try: 323 _ = get_reader_for(self.filename, 324 format=kwargs.get('format', None)) 325 except ValueError: 326 coordinatefile = args[1:] 327 else: 328 coordinatefile = (self.filename,) + args[1:] 329 else: 330 coordinatefile = args[1:] 331 332 if not coordinatefile: 333 coordinatefile = None 334 335 self.load_new(coordinatefile, **kwargs) 336 # parse transformations 337 trans_arg = kwargs.pop('transformations', None) 338 if trans_arg: 339 transforms =[trans_arg] if callable(trans_arg) else trans_arg 340 self.trajectory.add_transformations(*transforms) 341 342 # Check for guess_bonds 343 if kwargs.pop('guess_bonds', False): 344 self.atoms.guess_bonds(vdwradii=kwargs.pop('vdwradii', None)) 345 346 # None causes generic hash to get used. 347 # We store the name ieven if is_anchor is False in case the user later 348 # wants to make the universe an anchor. 349 self._anchor_name = kwargs.get('anchor_name', None) 350 # Universes are anchors by default 351 self.is_anchor = kwargs.get('is_anchor', True) 352 353 354 def copy(self): 355 """Return an independent copy of this Universe""" 356 new = self.__class__(self._topology.copy()) 357 new.trajectory = self.trajectory.copy() 358 return new 359 360 def _generate_from_topology(self): 361 # generate Universe version of each class 362 # AG, RG, SG, A, R, S 363 self._class_bases, self._classes = groups.make_classes() 364 365 # Put Group level stuff from topology into class 366 for attr in self._topology.attrs: 367 self._process_attr(attr) 368 369 # Generate atoms, residues and segments. 370 # These are the first such groups generated for this universe, so 371 # there are no cached merged classes yet. Otherwise those could be 372 # used directly to get a (very) small speedup. (Only really pays off 373 # the readability loss if instantiating millions of AtomGroups at 374 # once.) 375 self.atoms = AtomGroup(np.arange(self._topology.n_atoms), self) 376 377 self.residues = ResidueGroup( 378 np.arange(self._topology.n_residues), self) 379 380 self.segments = SegmentGroup( 381 np.arange(self._topology.n_segments), self) 382 383 # Update Universe namespace with segids 384 # Many segments can have same segid, so group together first 385 # 386 # DEPRECATED in 0.16.2 387 # REMOVE in 1.0 388 # See https://github.com/MDAnalysis/mdanalysis/issues/1377 389 try: 390 # returns dict of segid:segment 391 segids = self.segments.groupby('segids') 392 except AttributeError: 393 # no segids, don't do this step 394 pass 395 else: 396 for segid, segment in segids.items(): 397 if not segid: # ignore blank segids 398 continue 399 400 # cannot start attribute with number 401 if segid[0].isdigit(): 402 # prefix 's' if starts with number 403 name = 's' + segid 404 else: 405 name = segid 406 # if len 1 SegmentGroup, convert to Segment 407 if len(segment) == 1: 408 segment = segment[0] 409 self._instant_selectors[name] = segment 410 411 @classmethod 412 def empty(cls, n_atoms, n_residues=None, n_segments=None, 413 atom_resindex=None, residue_segindex=None, 414 trajectory=False, velocities=False, forces=False): 415 """Create a blank Universe 416 417 Useful for building a Universe without requiring existing files, 418 for example for system building. 419 420 If `trajectory` is set to True, a 421 :class:`MDAnalysis.coordinates.memory.MemoryReader` will be 422 attached to the Universe. 423 424 Parameters 425 ---------- 426 n_atoms : int 427 number of Atoms in the Universe 428 n_residues : int, optional 429 number of Residues in the Universe, defaults to 1 430 n_segments : int, optional 431 number of Segments in the Universe, defaults to 1 432 atom_resindex : array like, optional 433 mapping of atoms to residues, e.g. with 6 atoms, 434 `atom_resindex=[0, 0, 1, 1, 2, 2]` would put 2 atoms 435 into each of 3 residues. 436 residue_segindex : array like, optional 437 mapping of residues to segments 438 trajectory : bool, optional 439 if True, attaches a :class:`MDAnalysis.coordinates.memory.MemoryReader` 440 allowing coordinates to be set and written. Default is False 441 velocities : bool, optional 442 include velocities in the :class:`MDAnalysis.coordinates.memory.MemoryReader` 443 forces : bool, optional 444 include forces in the :class:`MDAnalysis.coordinates.memory.MemoryReader` 445 446 Returns 447 ------- 448 MDAnalysis.Universe object 449 450 Examples 451 -------- 452 For example to create a new Universe with 6 atoms in 2 residues, with 453 positions for the atoms and a mass attribute: 454 455 >>> u = mda.Universe.empty(6, 2, 456 atom_resindex=np.array([0, 0, 0, 1, 1, 1]), 457 trajectory=True, 458 ) 459 >>> u.add_TopologyAttr('masses') 460 461 .. versionadded:: 0.17.0 462 .. versionchanged:: 0.19.0 463 The attached Reader when trajectory=True is now a MemoryReader 464 """ 465 if n_residues is None: 466 n_residues = 1 467 elif atom_resindex is None: 468 warnings.warn( 469 'Multiple residues specified but no atom_resindex given. ' 470 'All atoms will be placed in first Residue.', 471 UserWarning) 472 if n_segments is None: 473 n_segments = 1 474 elif residue_segindex is None: 475 warnings.warn( 476 'Multiple segments specified but no segment_resindex given. ' 477 'All residues will be placed in first Segment', 478 UserWarning) 479 480 top = Topology(n_atoms, n_residues, n_segments, 481 atom_resindex=atom_resindex, 482 residue_segindex=residue_segindex, 483 ) 484 485 u = cls(top) 486 487 if trajectory: 488 coords = np.zeros((1, n_atoms, 3), dtype=np.float32) 489 dims = np.zeros(6, dtype=np.float64) 490 vels = np.zeros_like(coords) if velocities else None 491 forces = np.zeros_like(coords) if forces else None 492 493 # grab and attach a MemoryReader 494 u.trajectory = get_reader_for(coords)( 495 coords, order='fac', n_atoms=n_atoms, 496 dimensions=dims, velocities=vels, forces=forces) 497 498 return u 499 500 def __getattr__(self, key): 501 # This implements the instant selector of segments from a Universe. 502 # It is implemented as __getattr__ so a deprecation warning can be 503 # issued when the feature is used. Instant selectors are deprecated 504 # since version 0.16.2 and are tareted to be deleted in version 1.0. 505 # self._instant_selectors is populated in self._process_attr and 506 # created at the beginning of __init__. 507 try: 508 segment = self._instant_selectors[key] 509 except KeyError: 510 raise AttributeError('No attribute "{}".'.format(key)) 511 else: 512 warnings.warn("Instant selector Universe.<segid> " 513 "is deprecated and will be removed in 1.0. " 514 "Use SegmentGroup[SegmentGroup.segids == '<segid>'] " 515 "instead.", 516 DeprecationWarning) 517 return segment 518 519 @property 520 def universe(self): 521 # for Writer.write(universe), see Issue 49 522 # Encapsulation in an accessor prevents the Universe from 523 # having to keep a reference to itself, 524 # which might be undesirable if it has a __del__ method. 525 # It is also cleaner than a weakref. 526 return self 527 528 def load_new(self, filename, format=None, in_memory=False, **kwargs): 529 """Load coordinates from `filename`. 530 531 The file format of `filename` is autodetected from the file name suffix 532 or can be explicitly set with the `format` keyword. A sequence of files 533 can be read as a single virtual trajectory by providing a list of 534 filenames. 535 536 537 Parameters 538 ---------- 539 filename : str or list 540 the coordinate file (single frame or trajectory) *or* a list of 541 filenames, which are read one after another. 542 format : str or list or object (optional) 543 provide the file format of the coordinate or trajectory file; 544 ``None`` guesses it from the file extension. Note that this 545 keyword has no effect if a list of file names is supplied because 546 the "chained" reader has to guess the file format for each 547 individual list member [``None``]. Can also pass a subclass of 548 :class:`MDAnalysis.coordinates.base.ProtoReader` to define a custom 549 reader to be used on the trajectory file. 550 in_memory : bool (optional) 551 Directly load trajectory into memory with the 552 :class:`~MDAnalysis.coordinates.memory.MemoryReader` 553 554 .. versionadded:: 0.16.0 555 556 **kwargs : dict 557 Other kwargs are passed to the trajectory reader (only for 558 advanced use) 559 560 Returns 561 ------- 562 universe : Universe 563 564 Raises 565 ------ 566 TypeError if trajectory format can not be 567 determined or no appropriate trajectory reader found 568 569 570 .. versionchanged:: 0.8 571 If a list or sequence that is provided for `filename` only contains 572 a single entry then it is treated as single coordinate file. This 573 has the consequence that it is not read by the 574 :class:`~MDAnalysis.coordinates.chain.ChainReader` but directly by 575 its specialized file format reader, which typically has more 576 features than the 577 :class:`~MDAnalysis.coordinates.chain.ChainReader`. 578 579 .. versionchanged:: 0.17.0 580 Now returns a :class:`Universe` instead of the tuple of file/array 581 and detected file type. 582 """ 583 # filename==None happens when only a topology is provided 584 if filename is None: 585 return self 586 587 if len(util.asiterable(filename)) == 1: 588 # make sure a single filename is not handed to the ChainReader 589 filename = util.asiterable(filename)[0] 590 logger.debug("Universe.load_new(): loading {0}...".format(filename)) 591 592 try: 593 reader = get_reader_for(filename, format=format) 594 except ValueError as err: 595 raise TypeError( 596 "Cannot find an appropriate coordinate reader for file '{0}'.\n" 597 " {1}".format(filename, err)) 598 599 # supply number of atoms for readers that cannot do it for themselves 600 kwargs['n_atoms'] = self.atoms.n_atoms 601 602 self.trajectory = reader(filename, **kwargs) 603 if self.trajectory.n_atoms != len(self.atoms): 604 raise ValueError("The topology and {form} trajectory files don't" 605 " have the same number of atoms!\n" 606 "Topology number of atoms {top_n_atoms}\n" 607 "Trajectory: {fname} Number of atoms {trj_n_atoms}".format( 608 form=self.trajectory.format, 609 top_n_atoms=len(self.atoms), 610 fname=filename, 611 trj_n_atoms=self.trajectory.n_atoms)) 612 613 if in_memory: 614 self.transfer_to_memory(step=kwargs.get("in_memory_step", 1)) 615 616 return self 617 618 def transfer_to_memory(self, start=None, stop=None, step=None, 619 verbose=False): 620 """Transfer the trajectory to in memory representation. 621 622 Replaces the current trajectory reader object with one of type 623 :class:`MDAnalysis.coordinates.memory.MemoryReader` to support in-place 624 editing of coordinates. 625 626 Parameters 627 ---------- 628 start: int, optional 629 start reading from the nth frame. 630 stop: int, optional 631 read upto and excluding the nth frame. 632 step : int, optional 633 Read in every nth frame. [1] 634 verbose : bool, optional 635 Will print the progress of loading trajectory to memory, if 636 set to True. Default value is False. 637 638 639 .. versionadded:: 0.16.0 640 """ 641 from ..coordinates.memory import MemoryReader 642 643 if not isinstance(self.trajectory, MemoryReader): 644 n_frames = len(range( 645 *self.trajectory.check_slice_indices(start, stop, step) 646 )) 647 n_atoms = len(self.atoms) 648 pm_format = '{step}/{numsteps} frames copied to memory (frame {frame})' 649 pm = ProgressMeter(n_frames, interval=1, 650 verbose=verbose, format=pm_format) 651 coordinates = np.zeros((n_frames, n_atoms, 3), dtype=np.float32) 652 ts = self.trajectory.ts 653 has_vels = ts.has_velocities 654 has_fors = ts.has_forces 655 has_dims = ts.dimensions is not None 656 657 velocities = np.zeros_like(coordinates) if has_vels else None 658 forces = np.zeros_like(coordinates) if has_fors else None 659 dimensions = (np.zeros((n_frames, 6), dtype=np.float64) 660 if has_dims else None) 661 662 for i, ts in enumerate(self.trajectory[start:stop:step]): 663 np.copyto(coordinates[i], ts.positions) 664 if has_vels: 665 np.copyto(velocities[i], ts.velocities) 666 if has_fors: 667 np.copyto(forces[i], ts.forces) 668 if has_dims: 669 np.copyto(dimensions[i], ts.dimensions) 670 pm.echo(i, frame=ts.frame) 671 672 # Overwrite trajectory in universe with an MemoryReader 673 # object, to provide fast access and allow coordinates 674 # to be manipulated 675 if step is None: 676 step = 1 677 self.trajectory = MemoryReader( 678 coordinates, 679 dimensions=dimensions, 680 dt=self.trajectory.ts.dt * step, 681 filename=self.trajectory.filename, 682 velocities=velocities, 683 forces=forces, 684 ) 685 686 # python 2 doesn't allow an efficient splitting of kwargs in function 687 # argument signatures. 688 # In python3-only we'd be able to explicitly define this function with 689 # something like (sel, *othersels, updating=False, **selgroups) 690 def select_atoms(self, *args, **kwargs): 691 """Select atoms. 692 693 See Also 694 -------- 695 :meth:`MDAnalysis.core.groups.AtomGroup.select_atoms` 696 """ 697 return self.atoms.select_atoms(*args, **kwargs) 698 699 @property 700 def bonds(self): 701 """Bonds between atoms""" 702 return self.atoms.bonds 703 704 @property 705 def angles(self): 706 """Angles between atoms""" 707 return self.atoms.angles 708 709 @property 710 def dihedrals(self): 711 """Dihedral angles between atoms""" 712 return self.atoms.dihedrals 713 714 @property 715 def impropers(self): 716 """Improper dihedral angles between atoms""" 717 return self.atoms.impropers 718 719 @property 720 def anchor_name(self): 721 return self._gen_anchor_hash() 722 723 @anchor_name.setter 724 def anchor_name(self, name): 725 self.remove_anchor() # clear any old anchor 726 self._anchor_name = str(name) if not name is None else name 727 self.make_anchor() # add anchor again 728 729 def _gen_anchor_hash(self): 730 # hash used for anchoring. 731 # Try and use anchor_name, else use (and store) uuid 732 if self._anchor_name is not None: 733 return self._anchor_name 734 else: 735 try: 736 return self._anchor_uuid 737 except AttributeError: 738 # store this so we can later recall it if needed 739 self._anchor_uuid = uuid.uuid4() 740 return self._anchor_uuid 741 742 @property 743 def is_anchor(self): 744 """Is this Universe an anchoring for unpickling AtomGroups""" 745 return self._gen_anchor_hash() in _ANCHOR_UNIVERSES 746 747 @is_anchor.setter 748 def is_anchor(self, new): 749 if new: 750 self.make_anchor() 751 else: 752 self.remove_anchor() 753 754 def remove_anchor(self): 755 """Remove this Universe from the possible anchor list for unpickling""" 756 _ANCHOR_UNIVERSES.pop(self._gen_anchor_hash(), None) 757 758 def make_anchor(self): 759 _ANCHOR_UNIVERSES[self._gen_anchor_hash()] = self 760 761 def __repr__(self): 762 # return "<Universe with {n_atoms} atoms{bonds}>".format( 763 # n_atoms=len(self.atoms), 764 # bonds=" and {0} bonds".format(len(self.bonds)) if self.bonds else "") 765 766 return "<Universe with {n_atoms} atoms>".format( 767 n_atoms=len(self.atoms)) 768 769 def __getstate__(self): 770 raise NotImplementedError 771 772 def __setstate__(self, state): 773 raise NotImplementedError 774 775 # Properties 776 @property 777 def dimensions(self): 778 """Current dimensions of the unitcell""" 779 return self.coord.dimensions 780 781 @dimensions.setter 782 def dimensions(self, box): 783 """Set dimensions if the Timestep allows this 784 785 .. versionadded:: 0.9.0 786 """ 787 # Add fancy error handling here or use Timestep? 788 self.coord.dimensions = box 789 790 @property 791 def coord(self): 792 """Reference to current timestep and coordinates of universe. 793 794 The raw trajectory coordinates are :attr:`Universe.coord.positions`, 795 represented as a :class:`numpy.float32` array. 796 797 Because :attr:`coord` is a reference to a 798 :class:`~MDAnalysis.coordinates.base.Timestep`, it changes its contents 799 while one is stepping through the trajectory. 800 801 .. Note:: 802 803 In order to access the coordinates it is better to use the 804 :meth:`AtomGroup.positions` method; for instance, all coordinates of 805 the Universe as a numpy array: :meth:`Universe.atoms.positions`. 806 807 """ 808 return self.trajectory.ts 809 810 @property 811 def kwargs(self): 812 """keyword arguments used to initialize this universe""" 813 return copy.deepcopy(self._kwargs) 814 815 @property 816 def trajectory(self): 817 """Reference to trajectory reader object containing trajectory data.""" 818 if self._trajectory is not None: 819 return self._trajectory 820 else: 821 raise AttributeError("No trajectory loaded into Universe") 822 823 @trajectory.setter 824 def trajectory(self, value): 825 del self._trajectory # guarantees that files are closed (?) 826 self._trajectory = value 827 828 def add_TopologyAttr(self, topologyattr, values=None): 829 """Add a new topology attribute to the Universe 830 831 Adding a TopologyAttribute to the Universe makes it available to 832 all AtomGroups etc throughout the Universe. 833 834 Parameters 835 ---------- 836 topologyattr : TopologyAttr or string 837 Either a MDAnalysis TopologyAttr object or the name of a possible 838 topology attribute. 839 values : np.ndarray, optional 840 If initiating an attribute from a string, the initial values to 841 use. If not supplied, the new TopologyAttribute will have empty 842 or zero values. 843 844 Example 845 ------- 846 For example to add bfactors to a Universe: 847 848 >>> u.add_TopologyAttr('bfactors') 849 >>> u.atoms.bfactors 850 array([ 0., 0., 0., ..., 0., 0., 0.]) 851 852 .. versionchanged:: 0.17.0 853 Can now also add TopologyAttrs with a string of the name of the 854 attribute to add (eg 'charges'), can also supply initial values 855 using values keyword. 856 """ 857 if isinstance(topologyattr, six.string_types): 858 try: 859 tcls = _TOPOLOGY_ATTRS[topologyattr] 860 except KeyError: 861 raise ValueError( 862 "Unrecognised topology attribute name: '{}'." 863 " Possible values: '{}'\n" 864 "To raise an issue go to: http://issues.mdanalysis.org" 865 "".format( 866 topologyattr, ', '.join(sorted(_TOPOLOGY_ATTRS.keys()))) 867 ) 868 else: 869 topologyattr = tcls.from_blank( 870 n_atoms=self._topology.n_atoms, 871 n_residues=self._topology.n_residues, 872 n_segments=self._topology.n_segments, 873 values=values) 874 self._topology.add_TopologyAttr(topologyattr) 875 self._process_attr(topologyattr) 876 877 def _process_attr(self, attr): 878 """Squeeze a topologyattr for its information 879 880 Grabs: 881 - Group properties (attribute access) 882 - Component properties 883 - Transplant methods 884 """ 885 n_dict = {'atom': self._topology.n_atoms, 886 'residue': self._topology.n_residues, 887 'segment': self._topology.n_segments} 888 logger.debug("_process_attr: Adding {0} to topology".format(attr)) 889 if (attr.per_object is not None and len(attr) != n_dict[attr.per_object]): 890 raise ValueError('Length of {attr} does not' 891 ' match number of {obj}s.\n' 892 'Expect: {n:d} Have: {m:d}'.format( 893 attr=attr.attrname, 894 obj=attr.per_object, 895 n=n_dict[attr.per_object], 896 m=len(attr))) 897 898 for cls in attr.target_classes: 899 self._class_bases[cls]._add_prop(attr) 900 901 # TODO: Try and shove this into cls._add_prop 902 # Group transplants 903 for cls in (Atom, Residue, Segment, GroupBase, 904 AtomGroup, ResidueGroup, SegmentGroup): 905 for funcname, meth in attr.transplants[cls]: 906 setattr(self._class_bases[cls], funcname, meth) 907 # Universe transplants 908 for funcname, meth in attr.transplants['Universe']: 909 setattr(self.__class__, funcname, meth) 910 911 def add_Residue(self, segment=None, **attrs): 912 """Add a new Residue to this Universe 913 914 New Residues will not contain any Atoms, but can be assigned to Atoms 915 as per usual. If the Universe contains multiple segments, this must 916 be specified as a keyword. 917 918 Parameters 919 ---------- 920 segment : MDAnalysis.Segment 921 If there are multiple segments, then the Segment that the new 922 Residue will belong in must be specified. 923 attrs : dict 924 For each Residue attribute, the value for the new Residue must be 925 specified 926 927 Returns 928 ------- 929 A reference to the new Residue 930 931 Raises 932 ------ 933 NoDataError 934 If any information was missing. This happens before any changes have 935 been made, ie the change is rolled back. 936 937 938 Example 939 ------- 940 941 Adding a new GLY residue, then placing atoms within it: 942 943 >>> newres = u.add_Residue(segment=u.segments[0], resid=42, resname='GLY') 944 >>> u.atoms[[1, 2, 3]].residues = newres 945 >>> u.select_atoms('resname GLY and resid 42') 946 <AtomGroup with 3 atoms> 947 948 """ 949 if len(self.segments) == 1: # if only one segment, use this 950 segment = self.segments[0] 951 if segment is None: 952 raise NoDataError("") 953 # pass this information to the topology 954 residx = self._topology.add_Residue(segment, **attrs) 955 # resize my residues 956 self.residues = ResidueGroup(np.arange(self._topology.n_residues), self) 957 958 # return the new residue 959 return self.residues[residx] 960 961 def add_Segment(self, **attrs): 962 """Add a new Segment to this Universe 963 964 Parameters 965 ---------- 966 attrs : dict 967 For each Segment attribute as a key, give the value in the new 968 Segment 969 970 Returns 971 ------- 972 A reference to the new Segment 973 974 Raises 975 ------ 976 NoDataError 977 If any attributes were not specified as a keyword. 978 979 """ 980 # pass this information to the topology 981 segidx = self._topology.add_Segment(**attrs) 982 # resize my segments 983 self.segments = SegmentGroup(np.arange(self._topology.n_segments), self) 984 # return the new segment 985 return self.segments[segidx] 986 987 # TODO: Maybe put this as a Bond attribute transplant 988 # Problems: Can we transplant onto Universe? 989 # Probably a smarter way to do this too, could generate 990 # these on demand *per atom*. 991 # Wouldn't then need the Universe linkage here 992 # 993 # Alternate idea: Bonds Attribute generates a Fragments 994 # Attribute (ie, 2 for the price of 1) 995 # Fragments then gets its own Class/namespace/jazz. 996 @property 997 @cached('fragments') 998 def _fragdict(self): 999 """ 1000 .. versionadded:: 0.9.0 1001 .. versionchanged:: 0.16.0 1002 Fragment atoms are sorted by their index, and framgents are sorted 1003 by their first atom index so their order is predictable. 1004 .. versionchanged:: 0.19.0 1005 Uses faster C++ implementation 1006 """ 1007 atoms = self.atoms.ix 1008 bonds = self.atoms.bonds.to_indices() 1009 1010 frag_indices = find_fragments(atoms, bonds) 1011 frags = tuple([AtomGroup(np.sort(ix), self) for ix in frag_indices]) 1012 1013 fragdict = {} 1014 for f in frags: 1015 for a in f: 1016 fragdict[a] = f 1017 1018 return fragdict 1019 1020 1021# TODO: what is the point of this function??? 1022def as_Universe(*args, **kwargs): 1023 """Return a universe from the input arguments. 1024 1025 1. If the first argument is a universe, just return it:: 1026 1027 as_Universe(universe) --> universe 1028 1029 2. Otherwise try to build a universe from the first or the first 1030 and second argument:: 1031 1032 as_Universe(PDB, **kwargs) --> Universe(PDB, **kwargs) 1033 as_Universe(PSF, DCD, **kwargs) --> Universe(PSF, DCD, **kwargs) 1034 as_Universe(*args, **kwargs) --> Universe(*args, **kwargs) 1035 1036 Returns 1037 ------- 1038 :class:`~MDAnalysis.core.groups.Universe` 1039 """ 1040 if len(args) == 0: 1041 raise TypeError("as_Universe() takes at least one argument (%d given)" % len(args)) 1042 elif len(args) == 1 and isinstance(args[0], Universe): 1043 return args[0] 1044 return Universe(*args, **kwargs) 1045 1046 1047def Merge(*args): 1048 """Create a new new :class:`Universe` from one or more 1049 :class:`~MDAnalysis.core.groups.AtomGroup` instances. 1050 1051 Parameters 1052 ---------- 1053 *args : :class:`~MDAnalysis.core.groups.AtomGroup` 1054 One or more AtomGroups. 1055 1056 Returns 1057 ------- 1058 universe : :class:`Universe` 1059 1060 Raises 1061 ------ 1062 ValueError 1063 Too few arguments or an AtomGroup is empty and 1064 TypeError 1065 Arguments are not :class:`AtomGroup` instances. 1066 1067 Notes 1068 ----- 1069 The resulting :class:`Universe` will only inherit the common topology 1070 attributes that all merged universes share. 1071 1072 :class:`AtomGroup` instances can come from different Universes, or can come 1073 directly from a :meth:`~Universe.select_atoms` call. 1074 1075 :class:`Merge` can also be used with a single :class:`AtomGroup` if the 1076 user wants to, for example, re-order the atoms in the :class:`Universe`. 1077 1078 If multiple :class:`AtomGroup` instances from the same :class:`Universe` 1079 are given, the merge will first simply "add" together the 1080 :class:`AtomGroup` instances. 1081 1082 Merging does not create a full trajectory but only a single structure even 1083 if the input consists of one or more trajectories. However, one can use 1084 the :class:`~MDAnalysis.coordinates.memory.MemoryReader` to construct a 1085 trajectory for the new Universe as described under 1086 :ref:`creating-in-memory-trajectory-label`. 1087 1088 Example 1089 ------- 1090 In this example, protein, ligand, and solvent were externally prepared in 1091 three different PDB files. They are loaded into separate :class:`Universe` 1092 objects (where they could be further manipulated, e.g. renumbered, 1093 relabeled, rotated, ...) The :func:`Merge` command is used to combine all 1094 of them together:: 1095 1096 u1 = Universe("protein.pdb") 1097 u2 = Universe("ligand.pdb") 1098 u3 = Universe("solvent.pdb") 1099 u = Merge(u1.select_atoms("protein"), u2.atoms, u3.atoms) 1100 u.atoms.write("system.pdb") 1101 1102 The complete system is then written out to a new PDB file. 1103 1104 1105 .. versionchanged:: 0.9.0 1106 Raises exceptions instead of assertion errors. 1107 1108 .. versionchanged:: 0.16.0 1109 The trajectory is now a 1110 :class:`~MDAnalysis.coordinates.memory.MemoryReader`. 1111 1112 """ 1113 from ..topology.base import squash_by 1114 1115 if len(args) == 0: 1116 raise ValueError("Need at least one AtomGroup for merging") 1117 1118 for a in args: 1119 if not isinstance(a, groups.AtomGroup): 1120 raise TypeError(repr(a) + " is not an AtomGroup") 1121 for a in args: 1122 if len(a) == 0: 1123 raise ValueError("cannot merge empty AtomGroup") 1124 1125 # Create a new topology using the intersection of topology attributes 1126 blank_topology_attrs = set(dir(Topology(attrs=[]))) 1127 common_attrs = set.intersection(*[set(dir(ag.universe._topology)) 1128 for ag in args]) 1129 tops = set(['bonds', 'angles', 'dihedrals', 'impropers']) 1130 1131 attrs = [] 1132 1133 # Create set of attributes which are array-valued and can be simply 1134 # concatenated together 1135 common_array_attrs = common_attrs - blank_topology_attrs - tops 1136 # Build up array-valued topology attributes including only attributes 1137 # that all arguments' universes have 1138 for attrname in common_array_attrs: 1139 for ag in args: 1140 attr = getattr(ag, attrname) 1141 attr_class = type(getattr(ag.universe._topology, attrname)) 1142 if issubclass(attr_class, AtomAttr): 1143 pass 1144 elif issubclass(attr_class, ResidueAttr): 1145 attr = getattr(ag.residues, attrname) 1146 elif issubclass(attr_class, SegmentAttr): 1147 attr = getattr(ag.segments, attrname) 1148 else: 1149 raise NotImplementedError("Don't know how to handle" 1150 " TopologyAttr not subclassed" 1151 " from AtomAttr, ResidueAttr," 1152 " or SegmentAttr.") 1153 if type(attr) != np.ndarray: 1154 raise TypeError('Encountered unexpected topology ' 1155 'attribute of type {}'.format(type(attr))) 1156 try: 1157 attr_array.extend(attr) 1158 except NameError: 1159 attr_array = list(attr) 1160 attrs.append(attr_class(np.array(attr_array, dtype=attr.dtype))) 1161 del attr_array 1162 1163 # Build up topology groups including only those that all arguments' 1164 # universes have 1165 for t in (tops & common_attrs): 1166 offset = 0 1167 bondidx = [] 1168 types = [] 1169 for ag in args: 1170 # create a mapping scheme for this atomgroup 1171 mapping = {a.index: i for i, a in enumerate(ag, start=offset)} 1172 offset += len(ag) 1173 1174 tg = getattr(ag, t) 1175 bonds_class = type(getattr(ag.universe._topology, t)) 1176 # Create a topology group of only bonds that are within this ag 1177 # ie we don't want bonds that extend out of the atomgroup 1178 tg = tg.atomgroup_intersection(ag, strict=True) 1179 1180 # Map them so they refer to our new indices 1181 new_idx = [tuple([mapping[x] for x in entry]) for entry in tg.indices] 1182 bondidx.extend(new_idx) 1183 if hasattr(tg, '_bondtypes'): 1184 types.extend(tg._bondtypes) 1185 else: 1186 types.extend([None]*len(tg)) 1187 if any(t is None for t in types): 1188 attrs.append(bonds_class(bondidx)) 1189 else: 1190 types = np.array(types, dtype='|S8') 1191 attrs.append(bonds_class(bondidx, types)) 1192 1193 # Renumber residue and segment indices 1194 n_atoms = sum([len(ag) for ag in args]) 1195 residx = [] 1196 segidx = [] 1197 res_offset = 0 1198 seg_offset = 0 1199 for ag in args: 1200 # create a mapping scheme for this atomgroup's parents 1201 res_mapping = {r.resindex: i for i, r in enumerate(ag.residues, 1202 start=res_offset)} 1203 seg_mapping = {r.segindex: i for i, r in enumerate(ag.segments, 1204 start=seg_offset)} 1205 res_offset += len(ag.residues) 1206 seg_offset += len(ag.segments) 1207 1208 # Map them so they refer to our new indices 1209 residx.extend([res_mapping[x] for x in ag.resindices]) 1210 segidx.extend([seg_mapping[x] for x in ag.segindices]) 1211 1212 residx = np.array(residx, dtype=np.int32) 1213 segidx = np.array(segidx, dtype=np.int32) 1214 1215 _, _, [segidx] = squash_by(residx, segidx) 1216 1217 n_residues = len(set(residx)) 1218 n_segments = len(set(segidx)) 1219 1220 top = Topology(n_atoms, n_residues, n_segments, 1221 attrs=attrs, 1222 atom_resindex=residx, 1223 residue_segindex=segidx) 1224 1225 # Create and populate a universe 1226 coords = np.vstack([a.positions for a in args]) 1227 u = Universe(top, coords[None, :, :], 1228 format=MDAnalysis.coordinates.memory.MemoryReader) 1229 1230 return u 1231