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
|
# -*- Mode: python; tab-width: 4; indent-tabs-mode:nil; coding:utf-8 -*-
# vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4 fileencoding=utf-8
#
# 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 Lesser GNU Public Licence, v2.1 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
#
"""
TPR file format --- :mod:`MDAnalysis.coordinates.TPR`
======================================================
Class for reading positions and velocities from GROMACS TPR files.
Reading TPR files
-----------------
MDAnalysis can read positions and velocities from GROMACS TPR files,
and can also do so using different versions of the topology and
coordinate file for a given system.
For example, both of these are supported modes for reading
in positions and velocities from GROMACS TPR files::
>>> u = mda.Universe(TPR2020, TPR2024_4)
>>> u = mda.Universe(TPR2020)
Classes
-------
.. autoclass:: TPRReader
:members:
:inherited-members:
"""
from . import base
from ..lib import util
from .timestep import Timestep
import MDAnalysis.topology.tpr.utils as tpr_utils
import MDAnalysis.topology.tpr.setting as S
import logging
logger = logging.getLogger("MDAnalysis.coordinates.TPR")
import numpy as np
class TPRReader(base.SingleFrameReaderBase):
# TODO: reduce duplication with `TPRparser`;
# we could share some state for the position
# in the binary file to avoid re-reading topology
# or perhaps combine the topology and coordinate reading
# with some inheritance shenanigans?
"""Class supporting read in of positions and velocities from GROMACS TPR files.
.. versionadded:: 2.10.0
"""
format = "TPR"
units = {"length": "nm", "velocity": "nm/ps"}
_Timestep = Timestep
def _read_first_frame(self):
# Read header/move over topology
# TODO: reduce duplication with TPRparser perhaps...
with util.openany(self.filename, mode="rb") as infile:
tprf = infile.read()
data = tpr_utils.TPXUnpacker(tprf)
try:
th = tpr_utils.read_tpxheader(data) # tpxheader
except (EOFError, ValueError):
msg = f"{self.filename}: Invalid tpr coordinate file or cannot be recognized"
logger.critical(msg)
raise IOError(msg)
self.ts = ts = self._Timestep(th.natoms, **self._ts_kwargs)
self.n_atoms = th.natoms
# Starting with gromacs 2020 (tpx version 119), the body of the file
# is encoded differently. We change the unpacker accordingly.
if th.fver >= S.tpxv_AddSizeField and th.fgen >= 27:
actual_body_size = len(data.get_buffer()) - data.get_position()
if actual_body_size == 4 * th.sizeOfTprBody:
# See issue #2428.
msg = (
"TPR files produced with beta versions of gromacs 2020 "
"are not supported."
)
logger.critical(msg)
raise IOError(msg)
data = tpr_utils.TPXUnpacker2020.from_unpacker(data)
state_ngtc = th.ngtc # done init_state() in src/gmxlib/tpxio.c
if th.bBox:
tpr_utils.extract_box_info(data, th.fver)
if state_ngtc > 0:
if th.fver < 69: # redundancy due to different versions
tpr_utils.ndo_real(data, state_ngtc)
tpr_utils.ndo_real(
data, state_ngtc
) # relevant to Berendsen tcoupl_lambda
tpr_top = tpr_utils.do_mtop(
data, th.fver, tpr_resid_from_one=True, precision=th.precision
)
if th.bX:
self.ts._pos = np.asarray(
tpr_utils.ndo_rvec(data, th.natoms), dtype=np.float32
)
if th.bV:
self.ts.velocities = np.asarray(
tpr_utils.ndo_rvec(data, th.natoms), dtype=np.float32
)
self.ts.has_velocities = True
|
