Source code for MDAnalysis.coordinates.GMS

# -*- Mode: python; tab-width: 4; indent-tabs-mode:nil; -*-
# 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
#
#

"""GAMESS trajectory reader --- :mod:`MDAnalysis.coordinates.GMS`
=================================================================

Resources: the GMS output format is a common output format for different
GAMESS distributions: US-GAMESS, Firefly (PC-GAMESS) and GAMESS-UK.

Current version was approbated with US-GAMESS & Firefly only.

There appears to be no rigid format definition so it is likely users
will need to tweak the :class:`GMSReader`.

.. autoclass:: GMSReader
   :members:

"""
import os
import errno
import re

from . import base
import MDAnalysis.lib.util as util


[docs]class GMSReader(base.ReaderBase): """Reads from an GAMESS output file :Data: ts Timestep object containing coordinates of current frame :Methods: ``len(out)`` return number of frames in out ``for ts in out:`` iterate through trajectory Note ---- :class:`GMSReader` can read both uncompressed (``foo.out``) and compressed (``foo.out.bz2`` or ``foo.out.gz``) files; decompression is handled on the fly .. versionchanged:: 0.11.0 Frames now 0-based instead of 1-based. Added `dt` and `time_offset` keywords (passed to :class:`Timestep`). """ format = "GMS" # these are assumed! units = {'time': 'ps', 'length': 'Angstrom'} def __init__(self, outfilename, **kwargs): super(GMSReader, self).__init__(outfilename, **kwargs) # the filename has been parsed to be either b(g)zipped or not self.outfile = util.anyopen(self.filename) # note that, like for xtc and trr files, _n_atoms and _n_frames are used quasi-private variables # to prevent the properties being recalculated # this is because there is no indexing so the way it measures the number of frames is to read the whole file! self._n_atoms = None self._n_frames = None self._runtyp = None self.ts = self._Timestep(0) # need for properties initial calculations # update runtyp property self.runtyp if not self.runtyp in ['optimize', 'surface']: raise AttributeError('Wrong RUNTYP= '+self.runtyp) self.ts = self._Timestep(self.n_atoms, **self._ts_kwargs) # update n_frames property self.n_frames # Read in the first timestep self._read_next_timestep() @property def runtyp(self): """RUNTYP property of the GAMESS run""" if self._runtyp is not None: # return cached value return self._runtyp try: self._runtyp = self._determine_runtyp() except IOError: return 0 else: return self._runtyp def _determine_runtyp(self): with util.openany(self.filename) as out: for line in out: m = re.match(r'^.*RUNTYP=([A-Z]+)\s+.*', line) if m is not None: res = m.group(1).lower() break return res @property def n_atoms(self): """number of atoms in a frame""" if self._n_atoms is not None: # return cached value return self._n_atoms try: self._n_atoms = self._read_out_natoms() except IOError: return 0 else: return self._n_atoms def _read_out_natoms(self): with util.openany(self.filename) as out: for line in out: m = re.match(r'\s*TOTAL NUMBER OF ATOMS\s*=\s*([0-9]+)\s*',line) if m is not None: res = int(m.group(1)) break return res @property def n_frames(self): if self._n_frames is not None: # return cached value return self._n_frames try: self._n_frames = self._read_out_n_frames() except IOError: return 0 else: return self._n_frames def _read_out_n_frames(self): if self.runtyp == 'optimize': trigger = re.compile(b'^.NSERCH=.*') elif self.runtyp == 'surface': trigger = re.compile(b'^.COORD 1=.*') self._offsets = offsets = [] with util.openany(self.filename, 'rb') as out: line = True while not line == b'': # while not EOF line = out.readline() if re.match(trigger, line): offsets.append(out.tell() - len(line)) return len(offsets) def _read_frame(self, frame): self.outfile.seek(self._offsets[frame]) self.ts.frame = frame - 1 # gets +1'd in _read_next return self._read_next_timestep() def _read_next_timestep(self, ts=None): # check that the timestep object exists if ts is None: ts = self.ts # check that the outfile object exists; if not reopen the trajectory if self.outfile is None: self.open_trajectory() x = [] y = [] z = [] flag = 0 counter = 0 for line in self.outfile: if self.runtyp == 'optimize': if (flag == 0) and (re.match(r'^.NSERCH=.*', line) is not None): flag = 1 continue if (flag == 1) and (re.match(r'^ COORDINATES OF ALL ATOMS ARE ',\ line) is not None): flag = 2 continue if (flag == 2) and (re.match(r'^\s*-+\s*', line) is not None): flag = 3 continue if flag == 3 and counter < self.n_atoms: words = line.split() x.append(float(words[2])) y.append(float(words[3])) z.append(float(words[4])) counter += 1 elif self.runtyp == 'surface': if (flag == 0) and (re.match(\ r'^.COORD 1=\s*([-]?[0-9]+\.[0-9]+).*', line) is not None): flag = 1 continue if (flag == 1) and (re.match(\ r'^\s*HAS ENERGY VALUE\s*([-]?[0-9]+\.[0-9]+)\s*', line) is not None): flag = 3 continue if flag == 3 and counter < self.n_atoms: words = line.split() x.append(float(words[1])) y.append(float(words[2])) z.append(float(words[3])) counter += 1 # stop when the cursor has reached the end of that block if counter == self._n_atoms: ts._x[:] = x # more efficient to do it this way to avoid re-creating the numpy arrays ts._y[:] = y ts._z[:] = z #print ts.frame ts.frame += 1 return ts raise EOFError def _reopen(self): self.close() self.open_trajectory() def open_trajectory(self): if self.outfile is not None: raise IOError(errno.EALREADY, 'GMS file already opened', self.filename) if not os.path.exists(self.filename): # must check; otherwise might segmentation fault raise IOError(errno.ENOENT, 'GMS file not found', self.filename) self.outfile = util.anyopen(self.filename) # reset ts ts = self.ts ts.frame = -1 return self.outfile
[docs] def close(self): """Close out trajectory file if it was open.""" if self.outfile is None: return self.outfile.close() self.outfile = None