File: XDR.py

package info (click to toggle)
mdanalysis 2.10.0-1
  • links: PTS, VCS
  • area: main
  • in suites: forky, sid
  • size: 116,696 kB
  • sloc: python: 92,135; ansic: 8,156; makefile: 215; sh: 138
file content (358 lines) | stat: -rw-r--r-- 12,339 bytes parent folder | download
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
# -*- 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 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
#
"""\
XDR based trajectory files --- :mod:`MDAnalysis.coordinates.XDR`
================================================================

This module contains helper function and classes to read the XTC and TRR file
formats.

See Also
--------
MDAnalysis.coordinates.XTC: Read and write GROMACS XTC trajectory files.
MDAnalysis.coordinates.TRR: Read and write GROMACS TRR trajectory files.
MDAnalysis.lib.formats.libmdaxdr: Low level xdr format reader
"""

import errno
import numpy as np
from os.path import getctime, getsize, isfile, split, join
import warnings
from filelock import FileLock

from . import base
from ..lib.mdamath import triclinic_box
from ..lib.util import store_init_arguments


def offsets_filename(filename, ending="npz"):
    """Return offset or its lock filename for XDR files.
    For this the filename is appended
    with `_offsets.{ending}`.

    Parameters
    ----------
    filename : str
        filename of trajectory
    ending : str (optional)
        fileending of offsets file

    Returns
    -------
    offset_filename : str

    """
    head, tail = split(filename)
    return join(head, f".{tail}_offsets.{ending}")


def read_numpy_offsets(filename):
    """read offsets into dictionary.

    This assume offsets have been saved using numpy

    Parameters
    ----------
    filename : str
        filename of offsets

    Returns
    -------
    offsets : dict
        dictionary of offsets information

    """
    try:
        return {k: v for k, v in np.load(filename).items()}

    #  `ValueError` is encountered when the offset file is corrupted.
    except (ValueError, IOError):
        warnings.warn(f"Failed to load offsets file {filename}\n")
        return False


class XDRBaseReader(base.ReaderBase):
    """Base class for libmdaxdr file formats xtc and trr

    This class handles integration of XDR based formats into MDAnalysis. The
    XTC and TRR classes only implement `_write_next_frame` and
    `_frame_to_ts`.

    .. _offsets-label:

    Notes
    -----
    XDR based readers store persistent offsets on disk. The offsets are used to
    enable access to random frames efficiently. These offsets will be generated
    automatically the  first time the  trajectory is opened.  Generally offsets
    are stored  in hidden  `*_offsets.npz` files.  Afterwards opening  the same
    file again is fast. It sometimes can happen that the stored offsets get out
    off sync with the trajectory they refer to. For this the offsets also store
    the number of atoms, size of the file and last modification time. If any of
    them change  the offsets are recalculated.  Writing of the offset  file can
    fail when the  directory where the trajectory file resides  is not writable
    or if the  disk is full. In this  case a warning message will  be shown but
    the offsets will nevertheless be used during the lifetime of the trajectory
    Reader. However, the  next time the trajectory is opened,  the offsets will
    have to be rebuilt again.

    .. versionchanged:: 1.0.0
       XDR offsets read from trajectory if offsets file read-in fails
    .. versionchanged:: 2.0.0
       Add a InterProcessLock when generating offsets
    .. versionchanged:: 2.4.0
       Use a direct read into ts attributes
    .. versionchanged:: 2.9.0
       Changed fasteners.InterProcessLock() to filelock.FileLock
    """

    @store_init_arguments
    def __init__(
        self,
        filename,
        convert_units=True,
        sub=None,
        refresh_offsets=False,
        dt=None,
        **kwargs,
    ):
        """
        Parameters
        ----------
        filename : str
            trajectory filename
        convert_units : bool (optional)
            convert units to MDAnalysis units
        sub : array_like (optional)
            `sub` is an array of indices to pick out the corresponding
            coordinates and load only them; this requires that the topology
            itself is that of the sub system.
        refresh_offsets : bool (optional)
            force refresh of offsets
        dt : float (optional)
            timestep in MDAnalysis units to load trajectory with;
            if `dt` is ``None``, the time is taken from the xdr file;
            else, the time is set to `dt` * frame
        **kwargs : dict
            General reader arguments.

        """
        super(XDRBaseReader, self).__init__(
            filename, convert_units=convert_units, **kwargs
        )
        self._xdr = self._file(self.filename)

        self._sub = sub
        if self._sub is not None:
            self.n_atoms = len(self._sub)
        else:
            self.n_atoms = self._xdr.n_atoms

        if not refresh_offsets:
            self._load_offsets()
        else:
            self._read_offsets(store=True)
        frame = self._xdr.read()
        try:
            xdr_frame = self._xdr.read()
            if dt is None:
                dt = xdr_frame.time - frame.time
            else:
                self._ts_kwargs["dt"] = dt
            self._xdr.seek(1)
        except StopIteration:
            dt = 0

        self.ts = self._Timestep(self.n_atoms, **self._ts_kwargs)
        self._frame = 0
        self._frame_to_ts(frame, self.ts)
        # these should only be initialized once
        self.ts.dt = dt
        self.ts.dimensions = triclinic_box(*frame.box)
        if self.convert_units:
            if self.ts.dimensions is not None:
                self.convert_pos_from_native(self.ts.dimensions[:3])

    @classmethod
    def parse_n_atoms(cls, filename, **kwargs):
        with cls._file(filename) as f:
            n_atoms = f.n_atoms
        return n_atoms

    def close(self):
        """close reader"""
        self._xdr.close()

    def _load_offsets(self):
        """load frame offsets from file, reread them from the trajectory if that
        fails. To prevent the competition of generating the same offset file
        from multiple processes, an `InterProcessLock` is used."""
        fname = offsets_filename(self.filename)
        lock_name = offsets_filename(self.filename, ending="lock")

        #  check if the location of the lock is writable.
        try:
            with FileLock(lock_name) as filelock:
                pass
        except OSError as e:
            if isinstance(e, PermissionError) or e.errno == errno.EROFS:
                warnings.warn(
                    f"Cannot write lock/offset file in same location as "
                    f"{self.filename}. Using slow offset calculation."
                )
                self._read_offsets(store=False)
                return
            else:
                raise

        with FileLock(lock_name) as filelock:
            if not isfile(fname):
                self._read_offsets(store=True)
                return

            # if offsets file read correctly, data will be a dictionary of offsets
            # if not, data will be False
            # if False, offsets should be read from the trajectory
            # this warning can be avoided by loading Universe like:
            # u = mda.Universe(data.TPR, data.XTC, refresh_offsets=True)
            # refer to Issue #1893
            data = read_numpy_offsets(fname)
            if not data:
                warnings.warn(
                    f"Reading offsets from {fname} failed, "
                    "reading offsets from trajectory instead.\n"
                    "Consider setting 'refresh_offsets=True' "
                    "when loading your Universe."
                )
                self._read_offsets(store=True)
                return

            ctime_ok = size_ok = n_atoms_ok = False

            try:
                ctime_ok = getctime(self.filename) == data["ctime"]
                size_ok = getsize(self.filename) == data["size"]
                n_atoms_ok = self._xdr.n_atoms == data["n_atoms"]
            except KeyError:
                # we tripped over some old offset formated file
                pass

            if not (ctime_ok and size_ok and n_atoms_ok):
                warnings.warn(
                    "Reload offsets from trajectory\n "
                    "ctime or size or n_atoms did not match"
                )
                self._read_offsets(store=True)
            else:
                self._xdr.set_offsets(data["offsets"])

    def _read_offsets(self, store=False):
        """read frame offsets from trajectory"""
        fname = offsets_filename(self.filename)
        offsets = self._xdr.offsets
        if store:
            ctime = getctime(self.filename)
            size = getsize(self.filename)
            try:
                np.savez(
                    fname,
                    offsets=offsets,
                    size=size,
                    ctime=ctime,
                    n_atoms=self._xdr.n_atoms,
                )
            except Exception as e:
                warnings.warn(f"Couldn't save offsets because: {e}")

    @property
    def n_frames(self):
        """number of frames in trajectory"""
        return len(self._xdr)

    def _reopen(self):
        """reopen trajectory"""
        self.ts.frame = 0
        self._frame = -1
        offsets = self._xdr.offsets.copy()
        self._xdr.close()
        self._xdr.open(self.filename.encode("utf-8"), "r")
        # only restore in case we actually had offsets
        if len(offsets) != 0:
            self._xdr.set_offsets(offsets)

    def _read_frame(self, i):
        """read frame i"""
        self._frame = i - 1
        try:
            self._xdr.seek(i)
            timestep = self._read_next_timestep()
        except IOError:
            warnings.warn("seek failed, recalculating offsets and retrying")
            offsets = self._xdr.calc_offsets()
            self._xdr.set_offsets(offsets)
            self._read_offsets(store=True)
            self._xdr.seek(i)
            timestep = self._read_next_timestep()
        return timestep

    def Writer(self, filename, n_atoms=None, **kwargs):
        """Return writer for trajectory format"""
        if n_atoms is None:
            n_atoms = self.n_atoms
        return self._writer(filename, n_atoms=n_atoms, **kwargs)


class XDRBaseWriter(base.WriterBase):
    """Base class for libmdaxdr file formats xtc and trr"""

    def __init__(
        self, filename, n_atoms, convert_units=True, dt=None, **kwargs
    ):
        """
        Parameters
        ----------
        filename : str
            filename of trajectory
        n_atoms : int
            number of atoms to be written
        convert_units : bool (optional)
            convert from MDAnalysis units to format specific units
        dt : float (optional)
            timestep in MDAnalysis units to write trajectory with;
            if `dt` is ``None``, time for a frame is set from the timestep;
            else, the time for a frame is `dt` * frame
        **kwargs : dict
            General writer arguments
        """
        self.filename = filename
        self._convert_units = convert_units
        self.n_atoms = n_atoms
        self._xdr = self._file(self.filename, "w")
        self._dt = dt

    def close(self):
        """close trajectory"""
        self._xdr.close()

    def __del__(self):
        self.close()