/
universe.py
1658 lines (1375 loc) · 63 KB
/
universe.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
# -*- Mode: python; tab-width: 4; indent-tabs-mode:nil; coding: utf-8 -*-
# vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4
#
# MDAnalysis --- https://www.mdanalysis.org
# Copyright (c) 2006-2017 The MDAnalysis Development Team and contributors
# (see the file AUTHORS for the full list of names)
#
# Released under the GNU Public Licence, v2 or any higher version
#
# Please cite your use of MDAnalysis in published work:
#
# R. J. Gowers, M. Linke, J. Barnoud, T. J. E. Reddy, M. N. Melo, S. L. Seyler,
# D. L. Dotson, J. Domanski, S. Buchoux, I. M. Kenney, and O. Beckstein.
# MDAnalysis: A Python package for the rapid analysis of molecular dynamics
# simulations. In S. Benthall and S. Rostrup editors, Proceedings of the 15th
# Python in Science Conference, pages 102-109, Austin, TX, 2016. SciPy.
# doi: 10.25080/majora-629e541a-00e
#
# N. Michaud-Agrawal, E. J. Denning, T. B. Woolf, and O. Beckstein.
# MDAnalysis: A Toolkit for the Analysis of Molecular Dynamics Simulations.
# J. Comput. Chem. 32 (2011), 2319--2327, doi:10.1002/jcc.21787
#
"""\
=========================================================
Core object: Universe --- :mod:`MDAnalysis.core.universe`
=========================================================
The :class:`~MDAnalysis.core.universe.Universe` class ties a topology
and a trajectory together. Almost all code in MDAnalysis starts with a
``Universe``.
Normally, a ``Universe`` is created from files::
import MDAnalysis as mda
u = mda.Universe("topology.psf", "trajectory.dcd")
In order to construct new simulation system it is also convenient to
construct a ``Universe`` from existing
:class:`~MDAnalysis.core.group.AtomGroup` instances with the
:func:`Merge` function.
Classes
=======
.. autoclass:: Universe
:members:
Functions
=========
.. autofunction:: Merge
"""
import errno
import numpy as np
import logging
import copy
import warnings
import contextlib
import collections
import MDAnalysis
import sys
# When used in an MPI environment with Infiniband, importing MDAnalysis may
# trigger an MPI warning because importing the uuid module triggers a call to
# os.fork(). This happens if MPI_Init() has been called prior to importing
# MDAnalysis. The problem is actually caused by the uuid module and not by
# MDAnalysis itself. Python 3.7 fixes the problem.
import os
import uuid
from .. import _TOPOLOGY_ATTRS, _PARSERS
from ..exceptions import NoDataError
from ..lib import util
from ..lib.log import ProgressBar
from ..lib.util import cached, NamedStream, isstream
from ..lib.mdamath import find_fragments
from . import groups
from ._get_readers import get_reader_for, get_parser_for
from .groups import (ComponentBase, GroupBase,
Atom, Residue, Segment,
AtomGroup, ResidueGroup, SegmentGroup)
from .topology import Topology
from .topologyattrs import AtomAttr, ResidueAttr, SegmentAttr, BFACTOR_WARNING
from .topologyobjects import TopologyObject
logger = logging.getLogger("MDAnalysis.core.universe")
def _check_file_like(topology):
if isstream(topology):
if hasattr(topology, 'name'):
_name = topology.name
else:
_name = None
return NamedStream(topology, _name)
return topology
def _topology_from_file_like(topology_file, topology_format=None,
**kwargs):
parser = get_parser_for(topology_file, format=topology_format)
try:
with parser(topology_file) as p:
topology = p.parse(**kwargs)
except (IOError, OSError) as err:
# There are 2 kinds of errors that might be raised here:
# one because the file isn't present
# or the permissions are bad, second when the parser fails
if (err.errno is not None and
errno.errorcode[err.errno] in ['ENOENT', 'EACCES']):
# Runs if the error is propagated due to no permission / file not found
raise sys.exc_info()[1] from err
else:
# Runs when the parser fails
raise IOError("Failed to load from the topology file {0}"
" with parser {1}.\n"
"Error: {2}".format(topology_file, parser, err))
except (ValueError, NotImplementedError) as err:
raise ValueError(
"Failed to construct topology from file {0}"
" with parser {1}.\n"
"Error: {2}".format(topology_file, parser, err))
return topology
def _resolve_formats(*coordinates, format=None, topology_format=None):
if not coordinates:
if format is None:
format = topology_format
elif topology_format is None:
topology_format = format
return format, topology_format
def _resolve_coordinates(filename, *coordinates, format=None,
all_coordinates=False):
if all_coordinates or not coordinates and filename is not None:
try:
get_reader_for(filename, format=format)
except (ValueError, TypeError):
warnings.warn('No coordinate reader found for {}. Skipping '
'this file.'.format(filename))
else:
coordinates = (filename,) + coordinates
return coordinates
def _generate_from_topology(universe):
# generate Universe version of each class
# AG, RG, SG, A, R, S
universe._class_bases, universe._classes = groups.make_classes()
# Put Group level stuff from topology into class
for attr in universe._topology.attrs:
universe._process_attr(attr)
# Generate atoms, residues and segments.
# These are the first such groups generated for this universe, so
# there are no cached merged classes yet. Otherwise those could be
# used directly to get a (very) small speedup. (Only really pays off
# the readability loss if instantiating millions of AtomGroups at
# once.)
universe.atoms = AtomGroup(np.arange(universe._topology.n_atoms), universe)
universe.residues = ResidueGroup(
np.arange(universe._topology.n_residues), universe)
universe.segments = SegmentGroup(
np.arange(universe._topology.n_segments), universe)
class Universe(object):
"""The MDAnalysis Universe contains all the information describing the system.
The system always requires a *topology file* --- in the simplest case just
a list of atoms. This can be a CHARMM/NAMD PSF file or a simple coordinate
file with atom informations such as XYZ, PDB, GROMACS GRO or TPR, or CHARMM
CRD. See :ref:`Supported topology formats` for what kind of topologies can
be read.
A *trajectory file* provides coordinates; the coordinates have to be
ordered in the same way as the list of atoms in the topology. A trajectory
can be a single frame such as a PDB, CRD, or GRO file, or it can be a MD
trajectory (in CHARMM/NAMD/LAMMPS DCD, GROMACS XTC/TRR, AMBER nc, generic
XYZ format, ...). See :ref:`Supported coordinate formats` for what can be
read as a "trajectory".
As a special case, when the topology is a file that contains atom
information *and* coordinates (such as XYZ, PDB, GRO or CRD, see
:ref:`Supported coordinate formats`) then the coordinates are immediately
loaded from the "topology" file unless a trajectory is supplied.
Parameters
----------
topology: str, stream, Topology, numpy.ndarray, None
A CHARMM/XPLOR PSF topology file, PDB file or Gromacs GRO file; used to
define the list of atoms. If the file includes bond information,
partial charges, atom masses, ... then these data will be available to
MDAnalysis. Alternatively, an existing
:class:`MDAnalysis.core.topology.Topology` instance may be given,
numpy coordinates, or ``None`` for an empty universe.
coordinates: str, stream, list of str, list of stream (optional)
Coordinates can be provided as files of
a single frame (eg a PDB, CRD, or GRO file); a list of single
frames; or a trajectory file (in CHARMM/NAMD/LAMMPS DCD, Gromacs
XTC/TRR, or generic XYZ format). The coordinates must be
ordered in the same way as the list of atoms in the topology.
See :ref:`Supported coordinate formats` for what can be read
as coordinates. Alternatively, streams can be given.
topology_format: str, ``None``, default ``None``
Provide the file format of the topology file; ``None`` guesses it from
the file extension. Can also pass a subclass of
:class:`MDAnalysis.topology.base.TopologyReaderBase` to define a custom
reader to be used on the topology file.
format: str, ``None``, default ``None``
Provide the file format of the coordinate or trajectory file; ``None``
guesses it from the file extension. Note that this keyword has no
effect if a list of file names is supplied because the "chained" reader
has to guess the file format for each individual list member.
Can also pass a subclass of :class:`MDAnalysis.coordinates.base.ProtoReader`
to define a custom reader to be used on the trajectory file.
all_coordinates: bool, default ``False``
If set to ``True`` specifies that if more than one filename is passed
they are all to be used, if possible, as coordinate files (employing a
:class:`MDAnalysis.coordinates.chain.ChainReader`). The
default behavior is to take the first file as a topology and the
remaining as coordinates. The first argument will always always be used
to infer a topology regardless of *all_coordinates*.
guess_bonds: bool, default ``False``
Once Universe has been loaded, attempt to guess the connectivity
between atoms. This will populate the .bonds, .angles, and .dihedrals
attributes of the Universe.
vdwradii: dict, ``None``, default ``None``
For use with *guess_bonds*. Supply a dict giving a vdwradii for each
atom type which are used in guessing bonds.
fudge_factor: float, default [0.55]
For use with *guess_bonds*. Supply the factor by which atoms must
overlap each other to be considered a bond.
lower_bound: float, default [0.1]
For use with *guess_bonds*. Supply the minimum bond length.
transformations: function or list, ``None``, default ``None``
Provide a list of transformations that you wish to apply to the
trajectory upon reading. Transformations can be found in
:mod:`MDAnalysis.transformations`, or can be user-created.
in_memory: bool, default ``False``
After reading in the trajectory, transfer it to an in-memory
representations, which allow for manipulation of coordinates.
in_memory_step: int, default 1
Only read every nth frame into in-memory representation.
continuous: bool, default ``False``
The `continuous` option is used by the
:mod:`ChainReader<MDAnalysis.coordinates.chain>`, which contains the
functionality to treat independent trajectory files as a single virtual
trajectory.
**kwargs: extra arguments are passed to the topology parser.
Attributes
----------
trajectory : base.ReaderBase or base.SingleFrameReaderBase
currently loaded trajectory reader; readers are described in
:ref:`Coordinates`
dimensions : numpy.ndarray
system dimensions (simulation unit cell, if set in the
trajectory) at the *current time step*
(see :attr:`MDAnalysis.coordinates.timestep.Timestep.dimensions`).
The unit cell can be set for the current time step (but the change is
not permanent unless written to a file).
atoms : AtomGroup
all particles (:class:`~MDAnalysis.core.groups.Atom`) in the system,
as read from the `topology` file
residues : ResidueGroup
all residues (:class:`~MDAnalysis.core.groups.Residue`) in the system
segments : SegmentGroup
all segments (:class:`~MDAnalysis.core.groups.Segment`) in the system
bonds : topologyattrs.Bonds
all bonds (if defined in the `topology`) as provided by
:attr:`Universe.atoms.bonds`
angles : topologyattrs.Angles
all angles (if defined in the `topology`), same as
:attr:`Universe.atoms.angles`
dihedrals : topologyattrs.Dihedrals
all dihedral angles (if defined in the `topology`), same as
:attr:`Universe.atoms.dihedrals`
impropers : topologyattrs.Impropers
all improper dihedral angles (if defined in the `topology`), same as
:attr:`Universe.atoms.impropers`
Examples
--------
Examples for setting up a :class:`Universe`::
u = Universe(topology, trajectory) # read system from file(s)
u = Universe(pdbfile) # read atoms and coordinates from PDB or GRO
u = Universe(topology, [traj1, traj2, ...]) # read from a list of trajectories
u = Universe(topology, traj1, traj2, ...) # read from multiple trajectories
Load new data into a universe (replaces old trajectory and does *not* append)::
u.load_new(trajectory) # read from a new trajectory file
Selecting atoms with :meth:`~Universe.select_atoms` ::
ag = u.select_atoms(...)
returns an :class:`~MDAnalysis.core.groups.AtomGroup`.
.. versionchanged:: 1.0.0
Universe() now raises an error. Use Universe(None) or :func:`Universe.empty()` instead.
Removed instant selectors.
.. versionchanged:: 2.0.0
Universe now can be (un)pickled.
``topology`` and ``trajectory`` are reserved upon unpickle.
.. versionchanged:: 2.5.0
Added fudge_factor and lower_bound parameters for use with
*guess_bonds*.
"""
def __init__(self, topology=None, *coordinates, all_coordinates=False,
format=None, topology_format=None, transformations=None,
guess_bonds=False, vdwradii=None, fudge_factor=0.55,
lower_bound=0.1, in_memory=False,
in_memory_step=1, **kwargs):
self._trajectory = None # managed attribute holding Reader
self._cache = {'_valid': {}}
self.atoms = None
self.residues = None
self.segments = None
self.filename = None
self._kwargs = {
'transformations': transformations,
'guess_bonds': guess_bonds,
'vdwradii': vdwradii,
'fudge_factor': fudge_factor,
'lower_bound': lower_bound,
'in_memory': in_memory,
'in_memory_step': in_memory_step,
'format': format,
'topology_format': topology_format,
'all_coordinates': all_coordinates
}
self._kwargs.update(kwargs)
format, topology_format = _resolve_formats(*coordinates, format=format,
topology_format=topology_format)
if not isinstance(topology, Topology) and not topology is None:
self.filename = _check_file_like(topology)
topology = _topology_from_file_like(self.filename,
topology_format=topology_format,
**kwargs)
if topology is not None:
self._topology = topology
else:
# point to Universe.empty instead of making empty universe
raise TypeError('Topology argument required to make Universe. '
'Try Universe.empty(n_atoms, ...) to construct '
'your own Universe.')
_generate_from_topology(self) # make real atoms, res, segments
coordinates = _resolve_coordinates(self.filename, *coordinates,
format=format,
all_coordinates=all_coordinates)
if coordinates:
self.load_new(coordinates, format=format, in_memory=in_memory,
in_memory_step=in_memory_step, **kwargs)
if transformations:
if callable(transformations):
transformations = [transformations]
self._trajectory.add_transformations(*transformations)
if guess_bonds:
self.atoms.guess_bonds(vdwradii=vdwradii, fudge_factor=fudge_factor,
lower_bound=lower_bound)
def copy(self):
"""Return an independent copy of this Universe"""
new = self.__class__(self._topology.copy())
new.trajectory = self.trajectory.copy()
return new
@classmethod
def empty(cls, n_atoms, n_residues=1, n_segments=1, n_frames=1,
atom_resindex=None, residue_segindex=None,
trajectory=False, velocities=False, forces=False):
"""Create a blank Universe
Useful for building a Universe without requiring existing files,
for example for system building.
If `trajectory` is set to ``True``, a
:class:`MDAnalysis.coordinates.memory.MemoryReader` will be
attached to the Universe.
Parameters
----------
n_atoms: int
number of Atoms in the Universe
n_residues: int, default 1
number of Residues in the Universe, defaults to 1
n_segments: int, default 1
number of Segments in the Universe, defaults to 1
n_frames: int, default 1
number of Frames in the Universe, defaults to 1
atom_resindex: array like, optional
mapping of atoms to residues, e.g. with 6 atoms,
`atom_resindex=[0, 0, 1, 1, 2, 2]` would put 2 atoms
into each of 3 residues.
residue_segindex: array like, optional
mapping of residues to segments
trajectory: bool, optional
if ``True``, attaches a
:class:`MDAnalysis.coordinates.memory.MemoryReader` allowing
coordinates to be set and written.
velocities: bool, optional
include velocities in the
:class:`MDAnalysis.coordinates.memory.MemoryReader`
forces: bool, optional
include forces in the
:class:`MDAnalysis.coordinates.memory.MemoryReader`
Returns
-------
Universe
:class:`~MDAnalysis.core.universe.Universe` instance with dummy
values for atoms and undefined coordinates/velocities/forces
Examples
--------
For example to create a new Universe with 6 atoms in 2 residues, with
positions for the atoms and a mass attribute::
u = mda.Universe.empty(6, 2,
atom_resindex=np.array([0, 0, 0, 1, 1, 1]),
trajectory=True,
)
u.add_TopologyAttr('masses')
.. versionadded:: 0.17.0
.. versionchanged:: 0.19.0
The attached Reader when trajectory=True is now a MemoryReader
.. versionchanged:: 1.0.0
Universes can now be created with 0 atoms
"""
if not n_atoms:
n_residues = 0
n_segments = 0
if atom_resindex is None and n_residues > 1:
warnings.warn(
'Residues specified but no atom_resindex given. '
'All atoms will be placed in first Residue.',
UserWarning)
if residue_segindex is None and n_segments > 1:
warnings.warn(
'Segments specified but no segment_resindex given. '
'All residues will be placed in first Segment',
UserWarning)
top = Topology(n_atoms, n_residues, n_segments,
atom_resindex=atom_resindex,
residue_segindex=residue_segindex,
)
u = cls(top)
if n_frames > 1 or trajectory:
coords = np.zeros((n_frames, n_atoms, 3), dtype=np.float32)
vels = np.zeros_like(coords) if velocities else None
forces = np.zeros_like(coords) if forces else None
# grab and attach a MemoryReader
u.trajectory = get_reader_for(coords)(
coords, order='fac', n_atoms=n_atoms,
velocities=vels, forces=forces)
return u
@property
def universe(self):
# for Writer.write(universe), see Issue 49
# Encapsulation in an accessor prevents the Universe from
# having to keep a reference to itself,
# which might be undesirable if it has a __del__ method.
# It is also cleaner than a weakref.
return self
def load_new(self, filename, format=None, in_memory=False,
in_memory_step=1, **kwargs):
"""Load coordinates from `filename`.
The file format of `filename` is autodetected from the file name suffix
or can be explicitly set with the `format` keyword. A sequence of files
can be read as a single virtual trajectory by providing a list of
filenames.
Parameters
----------
filename: str or list
the coordinate file (single frame or trajectory) *or* a list of
filenames, which are read one after another.
format: str or list or object (optional)
provide the file format of the coordinate or trajectory file;
``None`` guesses it from the file extension. Note that this
keyword has no effect if a list of file names is supplied because
the "chained" reader has to guess the file format for each
individual list member [``None``]. Can also pass a subclass of
:class:`MDAnalysis.coordinates.base.ProtoReader` to define a custom
reader to be used on the trajectory file.
in_memory: bool (optional)
Directly load trajectory into memory with the
:class:`~MDAnalysis.coordinates.memory.MemoryReader`
.. versionadded:: 0.16.0
**kwargs: dict
Other kwargs are passed to the trajectory reader (only for
advanced use)
Returns
-------
universe: Universe
Raises
------
TypeError
if trajectory format can not be determined or no appropriate
trajectory reader found
.. versionchanged:: 0.8
If a list or sequence that is provided for `filename` only contains
a single entry then it is treated as single coordinate file. This
has the consequence that it is not read by the
:class:`~MDAnalysis.coordinates.chain.ChainReader` but directly by
its specialized file format reader, which typically has more
features than the
:class:`~MDAnalysis.coordinates.chain.ChainReader`.
.. versionchanged:: 0.17.0
Now returns a :class:`Universe` instead of the tuple of file/array
and detected file type.
.. versionchanged:: 2.4.0
Passes through kwargs if `in_memory=True`.
"""
# filename==None happens when only a topology is provided
if filename is None:
return self
if len(util.asiterable(filename)) == 1:
# make sure a single filename is not handed to the ChainReader
filename = util.asiterable(filename)[0]
logger.debug("Universe.load_new(): loading {0}...".format(filename))
try:
reader = get_reader_for(filename, format=format)
except ValueError as err:
raise TypeError(
"Cannot find an appropriate coordinate reader for file '{0}'.\n"
" {1}".format(filename, err))
# supply number of atoms for readers that cannot do it for themselves
kwargs['n_atoms'] = self.atoms.n_atoms
self.trajectory = reader(filename, format=format, **kwargs)
if self.trajectory.n_atoms != len(self.atoms):
raise ValueError("The topology and {form} trajectory files don't"
" have the same number of atoms!\n"
"Topology number of atoms {top_n_atoms}\n"
"Trajectory: {fname} Number of atoms {trj_n_atoms}".format(
form=self.trajectory.format,
top_n_atoms=len(self.atoms),
fname=filename,
trj_n_atoms=self.trajectory.n_atoms))
if in_memory:
self.transfer_to_memory(step=in_memory_step, **kwargs)
return self
def transfer_to_memory(self, start=None, stop=None, step=None,
verbose=False, **kwargs):
"""Transfer the trajectory to in memory representation.
Replaces the current trajectory reader object with one of type
:class:`MDAnalysis.coordinates.memory.MemoryReader` to support in-place
editing of coordinates.
Parameters
----------
start: int, optional
start reading from the nth frame.
stop: int, optional
read upto and excluding the nth frame.
step: int, optional
Read in every nth frame. [1]
verbose: bool, optional
Will print the progress of loading trajectory to memory, if
set to True. Default value is False.
.. versionadded:: 0.16.0
.. versionchanged:: 2.4.0
Passes through kwargs to MemoryReader
"""
from ..coordinates.memory import MemoryReader
if not isinstance(self.trajectory, MemoryReader):
n_frames = len(range(
*self.trajectory.check_slice_indices(start, stop, step)
))
n_atoms = len(self.atoms)
coordinates = np.zeros((n_frames, n_atoms, 3), dtype=np.float32)
ts = self.trajectory.ts
has_vels = ts.has_velocities
has_fors = ts.has_forces
has_dims = ts.dimensions is not None
velocities = np.zeros_like(coordinates) if has_vels else None
forces = np.zeros_like(coordinates) if has_fors else None
dimensions = (np.zeros((n_frames, 6), dtype=np.float32)
if has_dims else None)
for i, ts in enumerate(ProgressBar(self.trajectory[start:stop:step],
verbose=verbose,
desc="Loading frames")):
np.copyto(coordinates[i], ts.positions)
if has_vels:
np.copyto(velocities[i], ts.velocities)
if has_fors:
np.copyto(forces[i], ts.forces)
if has_dims:
np.copyto(dimensions[i], ts.dimensions)
# Overwrite trajectory in universe with an MemoryReader
# object, to provide fast access and allow coordinates
# to be manipulated
if step is None:
step = 1
self.trajectory = MemoryReader(
coordinates,
dimensions=dimensions,
dt=self.trajectory.ts.dt * step,
filename=self.trajectory.filename,
velocities=velocities,
forces=forces, **kwargs)
# python 2 doesn't allow an efficient splitting of kwargs in function
# argument signatures.
# In python3-only we'd be able to explicitly define this function with
# something like (sel, *othersels, updating=False, **selgroups)
def select_atoms(self, *args, **kwargs):
"""Select atoms.
See Also
--------
:meth:`MDAnalysis.core.groups.AtomGroup.select_atoms`
"""
return self.atoms.select_atoms(*args, **kwargs)
@property
def bonds(self):
"""Bonds between atoms.
:meta private:
"""
return self.atoms.bonds
@property
def angles(self):
"""Angles between atoms.
:meta private:
"""
return self.atoms.angles
@property
def dihedrals(self):
"""Dihedral angles between atoms.
:meta private:
"""
return self.atoms.dihedrals
@property
def impropers(self):
"""Improper dihedral angles between atoms.
:meta private:
"""
return self.atoms.impropers
def __repr__(self):
# return "<Universe with {n_atoms} atoms{bonds}>".format(
# n_atoms=len(self.atoms),
# bonds=" and {0} bonds".format(len(self.bonds)) if self.bonds else "")
return "<Universe with {n_atoms} atoms>".format(
n_atoms=len(self.atoms))
@classmethod
def _unpickle_U(cls, top, traj):
"""Special method used by __reduce__ to deserialise a Universe"""
# top is a Topology obj at this point, but Universe can handle that.
u = cls(top)
u.trajectory = traj
return u
def __reduce__(self):
# __setstate__/__getstate__ will raise an error when Universe has a
# transformation (that has AtomGroup inside). Use __reduce__ instead.
# Universe's two "legs" of top and traj both serialise themselves.
return (self._unpickle_U, (self._topology,
self._trajectory))
# Properties
@property
def dimensions(self):
"""Current dimensions of the unitcell.
:meta private:
"""
return self.coord.dimensions
@dimensions.setter
def dimensions(self, box):
"""Set dimensions if the Timestep allows this
.. versionadded:: 0.9.0
"""
# Add fancy error handling here or use Timestep?
self.coord.dimensions = box
@property
def coord(self):
"""Reference to current timestep and coordinates of universe.
The raw trajectory coordinates are :attr:`Universe.coord.positions`,
represented as a :class:`numpy.float32` array.
Because :attr:`coord` is a reference to a
:class:`~MDAnalysis.coordinates.timestep.Timestep`, it changes its contents
while one is stepping through the trajectory.
.. Note::
In order to access the coordinates it is better to use the
:meth:`AtomGroup.positions` method; for instance, all coordinates of
the Universe as a numpy array: :meth:`Universe.atoms.positions`.
"""
return self.trajectory.ts
@property
def kwargs(self):
"""keyword arguments used to initialize this universe"""
return copy.deepcopy(self._kwargs)
@property
def trajectory(self):
"""Reference to trajectory reader object containing trajectory data.
:meta private:
"""
if self._trajectory is not None:
return self._trajectory
else:
raise AttributeError("No trajectory loaded into Universe")
@trajectory.setter
def trajectory(self, value):
del self._trajectory # guarantees that files are closed (?)
self._trajectory = value
def add_TopologyAttr(self, topologyattr, values=None):
"""Add a new topology attribute to the Universe
Adding a TopologyAttribute to the Universe makes it available to
all AtomGroups etc throughout the Universe.
Parameters
----------
topologyattr: TopologyAttr or string
Either a MDAnalysis TopologyAttr object or the name of a possible
topology attribute.
values: np.ndarray, optional
If initiating an attribute from a string, the initial values to
use. If not supplied, the new TopologyAttribute will have empty
or zero values.
Example
-------
For example to add bfactors to a Universe:
>>> import MDAnalysis as mda
>>> from MDAnalysis.tests.datafiles import PSF, DCD
>>> u = mda.Universe(PSF, DCD)
>>> u.add_TopologyAttr('tempfactors')
>>> u.atoms.tempfactors
array([0., 0., 0., ..., 0., 0., 0.])
.. versionchanged:: 0.17.0
Can now also add TopologyAttrs with a string of the name of the
attribute to add (eg 'charges'), can also supply initial values
using values keyword.
.. versionchanged:: 1.1.0
Now warns when adding bfactors to a Universe with
existing tempfactors, or adding tempfactors to a
Universe with existing bfactors.
In version 2.0, MDAnalysis will stop treating
tempfactors and bfactors as separate attributes. Instead,
they will be aliases of the same attribute.
"""
if isinstance(topologyattr, str):
if topologyattr in ("bfactor", "bfactors"):
warnings.warn(BFACTOR_WARNING, DeprecationWarning)
try:
tcls = _TOPOLOGY_ATTRS[topologyattr]
except KeyError:
errmsg = (
"Unrecognised topology attribute name: '{}'."
" Possible values: '{}'\n"
"To raise an issue go to: "
"https://github.com/MDAnalysis/mdanalysis/issues"
"".format(
topologyattr, ', '.join(
sorted(_TOPOLOGY_ATTRS.keys()))))
raise ValueError(errmsg) from None
else:
topologyattr = tcls.from_blank(
n_atoms=self._topology.n_atoms,
n_residues=self._topology.n_residues,
n_segments=self._topology.n_segments,
values=values)
self._topology.add_TopologyAttr(topologyattr)
self._process_attr(topologyattr)
def del_TopologyAttr(self, topologyattr):
"""Remove a topology attribute from the Universe
Removing a TopologyAttribute from the Universe makes it unavailable to
all AtomGroups etc throughout the Universe.
Parameters
----------
topologyattr: TopologyAttr or string
Either a MDAnalysis TopologyAttr object or the name of a possible
topology attribute.
Example
-------
For example to remove bfactors to a Universe:
>>> import MDAnalysis as mda
>>> from MDAnalysis.tests.datafiles import PSF, DCD
>>> u = mda.Universe(PSF, DCD)
>>> u.add_TopologyAttr('tempfactors')
>>> hasattr(u.atoms[:3], 'tempfactors')
True
>>>
>>> u.del_TopologyAttr('tempfactors')
>>> hasattr(u.atoms[:3], 'tempfactors')
False
.. versionadded:: 2.0.0
"""
if not isinstance(topologyattr, str):
try:
topologyattr = topologyattr.attrname
except AttributeError:
# either TopologyGroup or not valid
try:
# this may not end well
# e.g. matrix -> matrices
topologyattr = topologyattr.btype + "s"
except AttributeError:
raise ValueError("Topology attribute must be str or "
"TopologyAttr object or class. "
f"Given: {type(topologyattr)}") from None
try:
topologyattr = _TOPOLOGY_ATTRS[topologyattr].attrname
except KeyError:
attrs = ', '.join(sorted(_TOPOLOGY_ATTRS))
errmsg = (f"Unrecognised topology attribute: '{topologyattr}'."
f" Possible values: '{attrs}'\n"
"To raise an issue go to: "
"https://github.com/MDAnalysis/mdanalysis/issues")
raise ValueError(errmsg) from None
try:
topattr = getattr(self._topology, topologyattr)
except AttributeError:
raise ValueError(f"Topology attribute {topologyattr} "
"not in Universe.") from None
self._topology.del_TopologyAttr(topattr)
self._unprocess_attr(topattr)
def _process_attr(self, attr):
"""Squeeze a topologyattr for its information
Grabs:
- Group properties (attribute access)
- Component properties
- Transplant methods
"""
n_dict = {'atom': self._topology.n_atoms,
'residue': self._topology.n_residues,
'segment': self._topology.n_segments}
logger.debug("_process_attr: Adding {0} to topology".format(attr))
if (attr.per_object is not None and len(attr) != n_dict[attr.per_object]):
raise ValueError('Length of {attr} does not'
' match number of {obj}s.\n'
'Expect: {n:d} Have: {m:d}'.format(
attr=attr.attrname,
obj=attr.per_object,
n=n_dict[attr.per_object],
m=len(attr)))
for cls in attr.target_classes:
self._class_bases[cls]._add_prop(attr)
# TODO: Try and shove this into cls._add_prop
# Group transplants
for cls in (Atom, Residue, Segment, GroupBase,
AtomGroup, ResidueGroup, SegmentGroup):
for funcname, meth in attr.transplants[cls]:
setattr(self._class_bases[cls], funcname, meth)
# Universe transplants
for funcname, meth in attr.transplants['Universe']:
setattr(self.__class__, funcname, meth)
def _unprocess_attr(self, attr):
"""
Undo all the stuff in _process_attr.
If the topology attribute is not present, nothing happens
(silent fail).
"""
for cls in attr.target_classes:
self._class_bases[cls]._del_prop(attr)
# Universe transplants
for funcname, _ in attr.transplants.pop("Universe", []):
delattr(self.__class__, funcname)
# Group transplants
for cls, transplants in attr.transplants.items():
for funcname, _ in transplants:
delattr(self._class_bases[cls], funcname)
def add_Residue(self, segment=None, **attrs):
"""Add a new Residue to this Universe
New Residues will not contain any Atoms, but can be assigned to Atoms
as per usual. If the Universe contains multiple segments, this must
be specified as a keyword.
Parameters
----------
segment: MDAnalysis.Segment
If there are multiple segments, then the Segment that the new
Residue will belong in must be specified.
attrs: dict
For each Residue attribute, the value for the new Residue must be
specified
Returns
-------
A reference to the new Residue
Raises
------
NoDataError
If any information was missing. This happens before any changes have
been made, ie the change is rolled back.
Example
-------
Adding a new GLY residue, then placing atoms within it: