Source code for MDAnalysis.auxiliary.XVG

# -*- 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
#

"""
XVG auxiliary reader --- :mod:`MDAnalysis.auxiliary.XVG`
========================================================

xvg files are produced by Gromacs during simulation or analysis, formatted
for plotting data with Grace.

Data is column-formatted; time/data selection is enabled by providing column
indices.

Note
----
By default, the time of each step is assumed to be stored in the first column,
in units of ps.


.. autoclass:: XVGStep
   :members:


XVG Readers
-----------
The default :class:`XVGReader` reads and stores the full contents of the .xvg
file on initialisation, while a second reader (:class:`XVGFileReader`) that
reads steps one at a time as required is also provided for when a lower memory
footprint is desired.

Note
----
Data is assumed to be time-ordered.

Multiple datasets, separated in the .xvg file by '&', are currently not
supported (the readers will stop at the first line starting '&').


.. autoclass:: XVGReader
   :members:

.. autoclass:: XVGFileReader
   :members:


.. autofunction:: uncomment

"""
import numbers
import os
import numpy as np
from . import base
from ..lib.util import anyopen


[docs] def uncomment(lines): """Remove comments from lines in an .xvg file Parameters ---------- lines : list of str Lines as directly read from .xvg file Yields ------ str The next non-comment line, with any trailing comments removed """ for line in lines: stripped_line = line.strip() # ignore blank lines if not stripped_line: continue # '@' must be at the beginning of a line to be a grace instruction if stripped_line[0] == "@": continue # '#' can be anywhere in the line, everything after is a comment comment_position = stripped_line.find("#") if comment_position > 0 and stripped_line[:comment_position]: yield stripped_line[:comment_position] elif comment_position < 0 and stripped_line: yield stripped_line
# if comment_position == 0, then the line is empty
[docs] class XVGStep(base.AuxStep): """AuxStep class for .xvg file format. Extends the base AuxStep class to allow selection of time and data-of-interest fields (by column index) from the full set of data read each step. Parameters ---------- time_selector : int | None, optional Index of column in .xvg file storing time, assumed to be in ps. Default value is 0 (i.e. first column). data_selector : list of int | None, optional List of indices of columns in .xvg file containing data of interest to be stored in ``data``. Default value is ``None``. **kwargs Other AuxStep options. See Also -------- :class:`~MDAnalysis.auxiliary.base.AuxStep` """ def __init__(self, time_selector=0, data_selector=None, **kwargs): super(XVGStep, self).__init__( time_selector=time_selector, data_selector=data_selector, **kwargs ) def _select_time(self, key): if key is None: # here so that None is a valid value; just return return if isinstance(key, numbers.Integral): return self._select_data(key) else: raise ValueError("Time selector must be single index") def _select_data(self, key): if key is None: # here so that None is a valid value; just return return if isinstance(key, numbers.Integral): try: return self._data[key] except IndexError: errmsg = ( f"{key} not a valid index for data with " f"{len(self._data)} columns" ) raise ValueError(errmsg) from None else: return np.array([self._select_data(i) for i in key])
[docs] class XVGReader(base.AuxReader): """Auxiliary reader to read data from an .xvg file. Default reader for .xvg files. All data from the file will be read and stored on initialisation. Parameters ---------- filename : str Location of the file containing the auxiliary data. **kwargs Other AuxReader options. See Also -------- :class:`~MDAnalysis.auxiliary.base.AuxReader` Note ---- The file is assumed to be of a size such that reading and storing the full contents is practical. """ format = "XVG" _Auxstep = XVGStep def __init__(self, filename, **kwargs): self._auxdata = os.path.abspath(filename) with anyopen(filename) as xvg_file: lines = xvg_file.readlines() auxdata_values = [] # remove comments before storing for i, line in enumerate(uncomment(lines)): if line.lstrip()[0] == "&": # multiple data sets not supported; stop at end of first set break auxdata_values.append([float(val) for val in line.split()]) # check the number of columns is consistent if len(auxdata_values[i]) != len(auxdata_values[0]): raise ValueError( "Step {0} has {1} columns instead of " "{2}".format(i, auxdata_values[i], auxdata_values[0]) ) self._auxdata_values = np.array(auxdata_values) self._n_steps = len(self._auxdata_values) super(XVGReader, self).__init__(**kwargs) def _memory_usage(self): return self._auxdata_values.nbytes def _read_next_step(self): """Read next auxiliary step and update ``auxstep``. Returns ------- AuxStep object Updated with the data for the new step. Raises ------ StopIteration When end of auxiliary data set is reached. """ auxstep = self.auxstep new_step = self.step + 1 if new_step < self.n_steps: auxstep._data = self._auxdata_values[new_step] auxstep.step = new_step return auxstep else: self.rewind() raise StopIteration def _go_to_step(self, i): """Move to and read i-th auxiliary step. Parameters ---------- i: int Step number (0-indexed) to move to Returns ------- :class:`XVGStep` Raises ------ ValueError If step index not in valid range. """ if i >= self.n_steps or i < 0: raise ValueError( "Step index {0} is not valid for auxiliary " "(num. steps {1})".format(i, self.n_steps) ) self.auxstep.step = i - 1 self.next() return self.auxstep
[docs] def read_all_times(self): """Get NumPy array of time at each step. Returns ------- NumPy array of float Time at each step. """ return self._auxdata_values[:, self.time_selector]
[docs] class XVGFileReader(base.AuxFileReader): """Auxiliary reader to read (one step at a time) from an .xvg file. An alternative XVG reader which reads each step from the .xvg file as needed (rather than reading and storing all from the start), for a lower memory footprint. Parameters ---------- filename : str Location of the file containing the auxiliary data. **kwargs Other AuxReader options. See Also -------- :class:`~MDAnalysis.auxiliary.base.AuxFileReader` Note ---- The default reader for .xvg files is :class:`XVGReader`. """ format = "XVG-F" _Auxstep = XVGStep def __init__(self, filename, **kwargs): super(XVGFileReader, self).__init__(filename, **kwargs) def _read_next_step(self): """Read next recorded step in xvg file and update ``auxstep``. Returns ------- AuxStep object Updated with the data for the new step. Raises ------ StopIteration When end of file or end of first data set is reached. """ line = next(self.auxfile) while True: if not line or (line.strip() and line.strip()[0] == "&"): # at end of file or end of first set of data (multiple sets # currently not supported) self.rewind() raise StopIteration # uncomment the line for uncommented in uncomment([line]): # line has data in it; add to auxstep + return auxstep = self.auxstep auxstep.step = self.step + 1 auxstep._data = [float(i) for i in uncommented.split()] # see if we've set n_cols yet... try: auxstep._n_cols except AttributeError: # haven't set n_cols yet; set now auxstep._n_cols = len(auxstep._data) if len(auxstep._data) != auxstep._n_cols: raise ValueError( f"Step {self.step} has " f"{len(auxstep._data)} columns instead " f"of {auxstep._n_cols}" ) return auxstep # line is comment only - move to next line = next(self.auxfile) def _count_n_steps(self): """Iterate through all steps to count total number. Returns ------- int Total number of steps """ if not self.constant_dt: # check if we've already iterated through to build _times list try: return len(self._times) except AttributeError: # might as well build _times now, since we'll need to iterate # through anyway self._times = self.read_all_times() return len(self.read_all_times()) else: # don't need _times; iterate here instead self._restart() count = 0 for step in self: count = count + 1 return count
[docs] def read_all_times(self): """Iterate through all steps to build times list. Returns ------- NumPy array of float Time of each step """ self._restart() times = [] for step in self: times.append(self.time) return np.array(times)