gromacs.formats

1 # GromacsWrapper: formats.py 2 # Copyright (c) 2009-2010 Oliver Beckstein <orbeckst@gmail.com> 3 # Released under the GNU Public License 3 (or higher, your choice) 4 # See the file COPYING for details. 5 6 """ 7 :mod:`gromacs.formats` -- Accessing various files 8 ================================================= 9 10 This module contains classes that represent data files on 11 disk. Typically one creates an instance and 12 13 - reads from a file using a :meth:`read` method, or 14 15 - populates the instance (in the simplest case with a :meth:`set` 16 method) and the uses the :meth:`write` method to write the data to 17 disk in the appropriate format. 18 19 For function data there typically also exists a :meth:`plot` method 20 which produces a graph (using matplotlib). 21 22 The module defines some classes that are used in other modules; they 23 do *not* make use of :mod:`gromacs.tools` or :mod:`gromacs.cbook` and 24 can be safely imported at any time. 25 26 27 Classes 28 ------- 29 30 .. autoclass:: XVG 31 :members: 32 .. autoclass:: NDX 33 :members: 34 .. autoclass:: uniqueNDX 35 :members: 36 .. autoclass:: GRO 37 :members: 38 39 (Not implemented yet) 40 """ 41 from __future__ import with_statement 42 43 __docformat__ = "restructuredtext en" 44 45 import os 46 import re 47 import warnings 48 import errno 49 import operator 50 51 import numpy 52 53 from odict import odict 54 55 import utilities 56 from gromacs import ParseError, AutoCorrectionWarning 57 58 import logging

59 60 -class XVG(utilities.FileUtils):

61 """Class that represents the numerical data in a grace xvg file. 62 63 All data must be numerical. :const:`NAN` and :const:`INF` values are 64 supported via python's :func:`float` builtin function. 65 66 The :attr:`~XVG.array` attribute can be used to access the the 67 array once it has been read and parsed. The :attr:`~XVG.ma` 68 attribute is a numpy masked array (good for plotting). 69 70 Conceptually, the file on disk and the XVG instance are considered the same 71 data. Whenever the filename for I/O (:meth:`XVG.read` and :meth:`XVG.write`) is 72 changed then the filename associated with the instance is also changed to reflect 73 the association between file and instance. 74 75 With the *permissive* = ``True`` flag one can instruct the file reader to skip 76 unparseable lines. In this case the line numbers of the skipped lines are stored 77 in :attr:`XVG.corrupted_lineno`. 78 79 A number of attributes are defined to give quick access to simple statistics such as 80 81 - :attr:`~XVG.mean`: mean of all data columns 82 - :attr:`~XVG.std`: standard deviation 83 - :attr:`~XVG.min`: minimum of data 84 - :attr:`~XVG.max`: maximum of data 85 - :attr:`~XVG.error`: error on the mean, taking correlation times into 86 account (see also :meth:`XVG.set_correlparameters`) 87 - :attr:`~XVG.tc`: correlation time of the data (assuming a simple 88 exponential decay of the fluctuations around the mean) 89 90 These attributes are numpy arrays that correspond to the data columns, 91 i.e. :attr:`XVG.array`[1:]. 92 93 .. Note:: - Only simple XY or NXY files are currently supported, *not* 94 Grace files that contain multiple data sets separated by '&'. 95 - Any kind of formatting (i.e. :program:`xmgrace` commands) is discarded. 96 """ 97 98 #: Default extension of XVG files. 99 default_extension = "xvg" 100 logger = logging.getLogger('gromacs.formats.XVG') # for pickling: must be class-level 101 102 #: If :attr:`XVG.savedata` is ``False`` then any attributes in 103 #: :attr:`XVG.__pickle_excluded` are *not* pickled as they are but simply 104 #: pickled with the default value. 105 __pickle_excluded = {'__array': None} # note class name un-mangling in __getstate__()! 106

107 - def __init__(self, filename=None, names=None, permissive=False, **kwargs):

108 """Initialize the class from a xvg file. 109 110 :Arguments: 111 *filename* 112 is the xvg file; it can only be of type XY or 113 NXY. If it is supplied then it is read and parsed 114 when :attr:`XVG.array` is accessed. 115 *names* 116 optional labels for the columns (currently only 117 written as comments to file); string with columns 118 separated by commas or a list of strings 119 *permissive* 120 ``False`` raises a :exc:`ValueError` and logs and errior 121 when encountering data lines that it cannot parse. 122 ``True`` ignores those lines and logs a warning---this is 123 a risk because it might read a corrupted input file [``False``] 124 *savedata* 125 ``True`` includes the data (:attr:`XVG.array`` and 126 associated caches) when the instance is pickled (see 127 :mod:`pickle`); this is oftens not desirable because the 128 data are already on disk (the xvg file *filename*) and the 129 resulting pickle file can become very big. ``False`` omits 130 those data from a pickle. [``False``] 131 132 """ 133 self.__array = None # cache for array (BIG) (used by XVG.array) 134 self.__cache = {} # cache for computed results 135 self.savedata = kwargs.pop('savedata', False) 136 if not filename is None: 137 self._init_filename(filename) # note: reading data from file is delayed until required 138 if names is None: 139 self.names = [] 140 else: 141 try: 142 self.names = names.split(',') 143 except AttributeError: 144 self.names = names 145 self.permissive = permissive 146 self.corrupted_lineno = None # must parse() first before this makes sense 147 # default number of data points for calculating correlation times via FFT 148 self.ncorrel = kwargs.pop('ncorrel', 25000) 149 self.__correlkwargs = {} # see set_correlparameters()

150

151 - def read(self, filename=None):

152 """Read and parse xvg file *filename*.""" 153 self._init_filename(filename) 154 self.parse()

155

156 - def write(self, filename=None):

157 """Write array to xvg file *filename* in NXY format. 158 159 .. Note:: Only plain files working at the moment, not compressed. 160 """ 161 self._init_filename(filename) 162 with utilities.openany(self.real_filename, 'w') as xvg: 163 xvg.write("# xmgrace compatible NXY data file\n" 164 "# Written by gromacs.formats.XVG()\n") 165 xvg.write("# :columns: %r\n" % self.names) 166 for xyy in self.array.T: 167 xyy.tofile(xvg, sep=" ", format="%-8s") # quick and dirty ascii output...--no compression! 168 xvg.write('\n')

169 170 @property

171 - def array(self):

172 """Represent xvg data as a (cached) numpy array. 173 174 The array is returned with column-first indexing, i.e. for a data file with 175 columns X Y1 Y2 Y3 ... the array a will be a[0] = X, a[1] = Y1, ... . 176 """ 177 if self.__array is None: 178 self.parse() 179 return self.__array

180 181 @property

182 - def ma(self):

183 """Represent data as a masked array. 184 185 The array is returned with column-first indexing, i.e. for a data file with 186 columns X Y1 Y2 Y3 ... the array a will be a[0] = X, a[1] = Y1, ... . 187 188 inf and nan are filtered via :func:`numpy.isfinite`. 189 """ 190 a = self.array 191 return numpy.ma.MaskedArray(a, mask=numpy.logical_not(numpy.isfinite(a)))

192 193 @property

194 - def mean(self):

195 """Mean value of all data columns.""" 196 return self.array[1:].mean(axis=1)

197 198 @property

199 - def std(self):

200 """Standard deviation from the mean of all data columns.""" 201 return self.array[1:].std(axis=1)

202 203 @property

204 - def min(self):

205 """Minimum of the data columns.""" 206 return self.array[1:].min(axis=1)

207 208 @property

209 - def max(self):

210 """Maximum of the data columns.""" 211 return self.array[1:].max(axis=1)

212

213 - def _tcorrel(self, nstep=100, **kwargs):

214 """Correlation "time" of data. 215 216 The 0-th column of the data is interpreted as a time and the 217 decay of the data is computed from the autocorrelation 218 function (using FFT). 219 """ 220 from numkit.timeseries import tcorrel 221 from gromacs.analysis.collections import Collection 222 t = self.array[0,::nstep] 223 r = Collection([tcorrel(t, Y, nstep=1, **kwargs) for Y in self.array[1:,::nstep]]) 224 return r

225

226 - def set_correlparameters(self, **kwargs):

227 """Set and change the parameters for calculations involving correlation functions. 228 229 :Keywords: 230 *nstep* 231 only process every *nstep* data point to speed up the FFT; if 232 left empty a default is chosen that produces roughly 25,000 data 233 points (or whatever is set in *ncorrel*) 234 *ncorrel* 235 If no *nstep* is supplied, aim at using *ncorrel* data points for 236 the FFT; sets :attr:`XVG.ncorrel`. 237 *force* 238 force recalculating correlation data even if cached values are 239 available 240 *kwargs* 241 see :func:`numkit.timeseries.tcorrel` for other options 242 243 .. SeeAlso: :attr:`XVG.error` for details and references. 244 """ 245 self.ncorrel = kwargs.pop('ncorrel', self.ncorrel) or 25000 246 nstep = kwargs.pop('nstep', None) 247 if nstep is None: 248 # good step size leads to ~25,000 data points 249 nstep = len(self.array[0])/self.ncorrel # needs the array so can take a while when loading 250 kwargs['nstep'] = nstep 251 self.__correlkwargs.update(kwargs) # only contains legal kw for numkit.timeseries.tcorrel or force 252 return self.__correlkwargs

253

254 - def _correlprop(self, key, **kwargs):

255 kwargs = self.set_correlparameters(**kwargs) 256 if not self.__cache.get('tcorrel',None) or kwargs.pop('force', False): 257 self.__cache['tcorrel'] = self._tcorrel(**kwargs) 258 return numpy.array(self.__cache['tcorrel'].get(key).tolist())

259 260 @property

261 - def error(self):

262 """Error on the mean of the data, taking the correlation time into account. 263 264 See Frenkel and Smit, Academic Press, San Diego 2002, p526: 265 266 error = sqrt(2*tc*acf[0]/T) 267 268 where acf() is the autocorrelation function of the fluctuations around 269 the mean, y-<y>, tc is the correlation time, and T the total length of 270 the simulation. 271 """ 272 return self._correlprop('sigma')

273 274 @property

275 - def tc(self):

276 """Correlation time of the data. 277 278 See :meth:`XVG.error` for details. 279 """ 280 return self._correlprop('tc')

281

282 - def parse(self):

283 """Read and cache the file as a numpy array. 284 285 The array is returned with column-first indexing, i.e. for a data file with 286 columns X Y1 Y2 Y3 ... the array a will be a[0] = X, a[1] = Y1, ... . 287 """ 288 self.corrupted_lineno = [] 289 # cannot use numpy.loadtxt() because xvg can have two types of 'comment' lines 290 with utilities.openany(self.real_filename) as xvg: 291 rows = [] 292 ncol = None 293 for lineno,line in enumerate(xvg): 294 line = line.strip() 295 if line.startswith(('#', '@')) or len(line) == 0: 296 continue 297 if line.startswith('&'): 298 raise NotImplementedError('%s: Multi-data not supported, only simple NXY format.' 299 % self.real_filename) 300 # parse line as floats 301 try: 302 row = map(float, line.split()) 303 except: 304 if self.permissive: 305 self.logger.warn("%s: SKIPPING unparsable line %d: %r", 306 self.real_filename, lineno+1, line) 307 self.corrupted_lineno.append(lineno+1) 308 continue 309 self.logger.error("%s: Cannot parse line %d: %r", 310 self.real_filename, lineno+1, line) 311 raise 312 # check for same number of columns as in previous step 313 if not ncol is None and len(row) != ncol: 314 if self.permissive: 315 self.logger.warn("%s: SKIPPING line %d with wrong number of columns: %r", 316 self.real_filename, lineno+1, line) 317 self.corrupted_lineno.append(lineno+1) 318 continue 319 errmsg = "%s: Wrong number of columns in line %d: %r" % (self.real_filename, lineno+1, line) 320 self.logger.error(errmsg) 321 raise IOError(errno.ENODATA, errmsg, self.real_filename) 322 # finally: a good line 323 ncol = len(row) 324 rows.append(row) 325 try: 326 self.__array = numpy.array(rows).transpose() # cache result 327 except: 328 self.logger.error("%s: Failed reading XVG file, possibly data corrupted. " 329 "Check the last line of the file...", self.real_filename) 330 raise

331

332 - def set(self, a):

333 """Set the *array* data from *a* (i.e. completely replace). 334 335 No sanity checks at the moment... 336 """ 337 self.__array = numpy.asarray(a)

338

339 - def plot(self, **kwargs):

340 """Plot xvg file data. 341 342 The first column of the data is always taken as the abscissa 343 X. Additional columns are plotted as ordinates Y1, Y2, ... 344 345 In the special case that there is only a single column then this column 346 is plotted against the index, i.e. (N, Y). 347 348 :Keywords: 349 *columns* : list 350 Select the columns of the data to be plotted; the list 351 is used as a numpy.array extended slice. The default is 352 to use all columns. Columns are selected *after* a transform. 353 *transform* : function 354 function ``transform(array) -> array`` which transforms 355 the original array; must return a 2D numpy array of 356 shape [X, Y1, Y2, ...] where X, Y1, ... are column 357 vectors. By default the transformation is the 358 identity [``lambda x: x``]. 359 *maxpoints* : int 360 limit the total number of data points; matplotlib has issues processing 361 png files with >100,000 points and pdfs take forever to display. Set to 362 ``None`` if really all data should be displayed. At the moment we simply 363 subsample the data at regular intervals. [10000] 364 *kwargs* 365 All other keyword arguments are passed on to :func:`pylab.plot`. 366 """ 367 import pylab 368 369 maxpoints_default = 10000 370 columns = kwargs.pop('columns', Ellipsis) # slice for everything 371 maxpoints = kwargs.pop('maxpoints', maxpoints_default) 372 transform = kwargs.pop('transform', lambda x: x) # default is identity transformation 373 a = numpy.asarray(transform(self.array))[columns] # (slice o transform)(array) 374 375 ny = a.shape[-1] # assume 1D or 2D array with last dimension varying fastest 376 if not maxpoints is None and ny > maxpoints: 377 # reduce size by subsampling (primitive --- can leave out 378 # bits at the end or end up with almost twice of maxpoints) 379 stepsize = int(ny / maxpoints) 380 a = a[..., ::stepsize] 381 if maxpoints == maxpoints_default: # only warn if user did not set maxpoints 382 warnings.warn("Plot had %d datapoints > maxpoints = %d; subsampled to %d regularly spaced points." 383 % (ny, maxpoints, a.shape[-1]), category=AutoCorrectionWarning) 384 385 if len(a.shape) == 1: 386 # special case: plot against index; plot would do this automatically but 387 # we'll just produce our own xdata and pretend that this was X all along 388 X = numpy.arange(len(a)) 389 a = numpy.concatenate([[X], [a]]) # does NOT overwrite original a but make a new one 390 391 # now deal with infs, nans etc AFTER all transformations (needed for plotting across inf/nan) 392 ma = numpy.ma.MaskedArray(a, mask=numpy.logical_not(numpy.isfinite(a))) 393 394 # finally plot 395 kwargs['xdata'] = ma[0] # abscissa set separately 396 pylab.plot(ma[1:].T, **kwargs) # plot all other columns in parallel

397

398 - def errorbar(self, **kwargs):

399 """Quick hack: errorbar plot. 400 401 Set columns to select [x, y, dy]. 402 """ 403 import pylab 404 405 kwargs.setdefault('capsize', 0) 406 kwargs.setdefault('elinewidth', 1) 407 kwargs.setdefault('alpha', 0.3) 408 kwargs.setdefault('fmt', None) 409 410 maxpoints_default = 10000 411 columns = kwargs.pop('columns', Ellipsis) # slice for everything 412 maxpoints = kwargs.pop('maxpoints', maxpoints_default) 413 transform = kwargs.pop('transform', lambda x: x) # default is identity transformation 414 a = numpy.asarray(transform(self.array))[columns] # (slice o transform)(array) 415 416 ny = a.shape[-1] # assume 1D or 2D array with last dimension varying fastest 417 if not maxpoints is None and ny > maxpoints: 418 # reduce size by subsampling (primitive --- can leave out 419 # bits at the end or end up with almost twice of maxpoints) 420 stepsize = int(ny / maxpoints) 421 a = a[..., ::stepsize] 422 if maxpoints == maxpoints_default: # only warn if user did not set maxpoints 423 warnings.warn("Plot had %d datapoints > maxpoints = %d; subsampled to %d regularly spaced points." 424 % (ny, maxpoints, a.shape[-1]), category=AutoCorrectionWarning) 425 426 if len(a.shape) == 1: 427 # special case: plot against index; plot would do this automatically but 428 # we'll just produce our own xdata and pretend that this was X all along 429 X = numpy.arange(len(a)) 430 a = numpy.concatenate([[X], [a]]) # does NOT overwrite original a but make a new one 431 432 # now deal with infs, nans etc AFTER all transformations (needed for plotting across inf/nan) 433 ma = numpy.ma.MaskedArray(a, mask=numpy.logical_not(numpy.isfinite(a))) 434 435 # finally plot 436 X = ma[0] # abscissa set separately 437 Y = ma[1] 438 try: 439 kwargs['yerr'] = ma[3] 440 kwargs['xerr'] = ma[2] 441 except IndexError: 442 kwargs['yerr'] = ma[2] 443 444 pylab.errorbar(X, Y, **kwargs)

445

446 - def __getstate__(self):

447 """custom pickling protocol: http://docs.python.org/library/pickle.html 448 449 If :attr:`XVG.savedata` is ``False`` then any attributes in 450 :attr:`XVG.__pickle_excluded` are *not* pickled as they are but simply 451 pickled with the default value. 452 """ 453 if self.savedata: 454 d = self.__dict__ 455 else: 456 # do not pickle the big array cache 457 mangleprefix = '_'+self.__class__.__name__ 458 def demangle(k): 459 """_XVG__array --> __array""" 460 if k.startswith(mangleprefix): 461 k = k.replace(mangleprefix,'') 462 return k

463 d = {} 464 for k in self.__dict__: 465 d[k] = self.__pickle_excluded.get(demangle(k), self.__dict__[k]) 466 return d

467

468 - def __setstate__(self, d):

469 # compatibility with older (pre 0.1.13) pickled instances 470 if not 'savedata' in d: 471 wmsg = "Reading pre 0.1.13 pickle file: setting savedata=False" 472 warnings.warn(wmsg, category=DeprecationWarning) 473 self.logger.warn(wmsg) 474 d['savedata'] = False # new default 475 self.__dict__.update(d)

476

477 478 -class NDX(odict, utilities.FileUtils):

479 """Gromacs index file. 480 481 Represented as a ordered dict where the keys are index group names and 482 values are numpy arrays of atom numbers. 483 484 Use the :meth:`NDX.read` and :meth:`NDX.write` methods for 485 I/O. Access groups by name via the :meth:`NDX.get` and 486 :meth:`NDX.set` methods. 487 488 Alternatively, simply treat the :class:`NDX` instance as a 489 dictionary. Setting a key automatically transforms the new value 490 into a integer 1D numpy array (*not* a set, as would be the 491 :program:`make_ndx` behaviour). 492 493 .. Note:: The index entries themselves are ordered and can contain 494 duplicates so that output from NDX can be easily used for 495 :program:`g_dih` and friends. If you need set-like behaviour 496 you will have do use :class:`gromacs.formats.uniqueNDX` or 497 :class:`gromacs.cbook.IndexBuilder` (which uses 498 :program:`make_ndx` throughout). 499 500 **Example** 501 502 Read index file, make new group and write to disk:: 503 504 ndx = NDX() 505 ndx.read('system.ndx') 506 print ndx['Protein'] 507 ndx['my_group'] = [2, 4, 1, 5] # add new group 508 ndx.write('new.ndx') 509 510 Or quicker (replacing the input file ``system.ndx``):: 511 512 ndx = NDX('system') # suffix .ndx is automatically added 513 ndx['chi1'] = [2, 7, 8, 10] 514 ndx.write() 515 516 """ 517 default_extension = "ndx" 518 519 # match: [ index_groupname ] 520 SECTION = re.compile("""\s*\[\s*(?P<name>\S.*\S)\s*\]\s*""") 521 522 #: standard ndx file format: 15 columns 523 ncol = 15 524 #: standard ndx file format: '%6d' 525 format = '%6d' 526

527 - def __init__(self, filename=None, **kwargs):

528 super(NDX, self).__init__(**kwargs) # can use kwargs to set dict! (but no sanity checks!) 529 530 if not filename is None: 531 self._init_filename(filename) 532 self.read(filename)

533

534 - def read(self, filename=None):

535 """Read and parse index file *filename*.""" 536 self._init_filename(filename) 537 538 data = odict() 539 with open(self.real_filename) as ndx: 540 current_section = None 541 for line in ndx: 542 line = line.strip() 543 if len(line) == 0: 544 continue 545 m = self.SECTION.match(line) 546 if m: 547 current_section = m.group('name') 548 data[current_section] = [] # can fail if name not legal python key 549 continue 550 if not current_section is None: 551 data[current_section].extend(map(int, line.split())) 552 553 super(NDX,self).update(odict([(name, self._transform(atomnumbers)) 554 for name, atomnumbers in data.items()]))

555

556 - def write(self, filename=None, ncol=ncol, format=format):

557 """Write index file to *filename* (or overwrite the file that the index was read from)""" 558 with open(self.filename(filename, ext='ndx'), 'w') as ndx: 559 for name in self: 560 atomnumbers = self._getarray(name) # allows overriding 561 ndx.write('[ %s ]\n' % name) 562 for k in xrange(0, len(atomnumbers), ncol): 563 line = atomnumbers[k:k+ncol].astype(int) # nice formatting in ncol-blocks 564 n = len(line) 565 ndx.write((" ".join(n*[format])+'\n') % tuple(line)) 566 ndx.write('\n')

567

568 - def get(self, name):

569 """Return index array for index group *name*.""" 570 return self[name]

571

572 - def set(self, name, value):

573 """Set or add group *name* as a 1D numpy array.""" 574 self[name] = value

575

576 - def size(self, name):

577 """Return number of entries for group *name*.""" 578 return len(self[name])

579 580 @property

581 - def groups(self):

582 """Return a list of all groups.""" 583 return self.keys()

584 585 @property

586 - def sizes(self):

587 """Return a dict with group names and number of entries,""" 588 return dict([(name, len(atomnumbers)) for name, atomnumbers in self.items()])

589 590 @property

591 - def ndxlist(self):

592 """Return a list of groups in the same format as :func:`gromacs.cbook.get_ndx_groups`. 593 594 Format: 595 [ {'name': group_name, 'natoms': number_atoms, 'nr': # group_number}, ....] 596 """ 597 return [{'name': name, 'natoms': len(atomnumbers), 'nr': nr+1} for 598 nr,(name,atomnumbers) in enumerate(self.items())]

599

600 - def _getarray(self, name):

601 """Helper getter that is used in write(). 602 Override when using a _transform that stores something that 603 cannot be indexed, e.g. when using set()s. 604 """ 605 return self[name]

606

607 - def _transform(self, v):

608 """Transform input to the stored representation. 609 610 Override eg with ``return set(v)`` for index lists as sets. 611 """ 612 return numpy.ravel(v).astype(int)

613

614 - def __setitem__(self, k, v):

615 super(NDX, self).__setitem__(k, self._transform(v))

616

617 - def setdefault(*args,**kwargs):

618 raise NotImplementedError

619

620 621 -class IndexSet(set):

622 """set which defines '+' as union (OR) and '-' as intersection (AND)."""

623 - def __add__(self, x):

624 return self.union(x)

625 - def __sub__(self, x):

626 return self.intersection(x)

627

628 629 -class uniqueNDX(NDX):

630 """Index that behaves like make_ndx, i.e. entries behaves as sets, 631 not lists. 632 633 The index lists behave like sets: 634 - adding sets with '+' is equivalent to a logical OR: x + y == "x | y" 635 - subtraction '-' is AND: x - y == "x & y" 636 - see :meth:`~gromacs.formats.join` for ORing multiple groups (x+y+z+...) 637 638 **Example** :: 639 I = uniqueNDX('system.ndx') 640 I['SOLVENT'] = I['SOL'] + I['NA+'] + I['CL-'] 641 """ 642

643 - def join(self, *groupnames):

644 """Return an index group that contains atoms from all *groupnames*. 645 646 The method will silently ignore any groups that are not in the 647 index. 648 649 **Example** 650 651 Always make a solvent group from water and ions, even if not 652 all ions are present in all simulations:: 653 654 I['SOLVENT'] = I.join('SOL', 'NA+', 'K+', 'CL-') 655 """ 656 return self._sum([self[k] for k in groupnames if k in self])

657

658 - def _sum(self, sequence):

659 return reduce(operator.add, sequence)

660

661 - def _transform(self, v):

662 return IndexSet(v)

663

664 - def _getarray(self, k):

665 return numpy.sort(numpy.fromiter(self[k],dtype=int,count=len(self[k])))

666

667 668 669 # or use list of these? 670 # class IndexGroup(dict): 671 # def __init__(self, groupnumber=None, name="empty", atomnumbers=None, **kwargs): 672 # atomnumbers = atomnumbers or [] 673 # _atomnumbers = numpy.asarray(atomnumbers).astype(int) 674 # super(IndexGroup, self).__init__(name=str(name), 675 # atomnumbers=_atomnumbers, 676 # nr=groupnumber) 677 678 -class GRO(utilities.FileUtils):

679 """Class that represents a GROMOS (gro) structure file. 680 681 682 File format: 683 """ 684 default_extension = "gro" 685 logger = logging.getLogger('gromacs.formats.GRO') 686

687 - def __init__(self, **kwargs):

688 689 raise NotImplementedError 690 691 filename = kwargs.pop('filename',None) 692 super(GRO, self).__init__(**kwargs) 693 694 if not filename is None: 695 self._init_filename(filename) 696 self.read(filename)

697

698 - def read(self, filename=None):

699 """Read and parse index file *filename*.""" 700 self._init_filename(filename) 701 702 with open(self.real_filename) as gro: 703 pass

704

705 706 707 -class MDP(odict, utilities.FileUtils):

708 """Class that represents a Gromacs mdp run input file. 709 710 The MDP instance is an ordered dictionary. 711 712 - *Parameter names* are keys in the dictionary. 713 - *Comments* are sequentially numbered with keys Comment0001, 714 Comment0002, ... 715 - *Empty lines* are similarly preserved as Blank0001, .... 716 717 When writing, the dictionary is dumped in the recorded order to a 718 file. Inserting keys at a specific position is not possible. 719 720 Currently, comments after a parameter on the same line are 721 discarded. Leading and trailing spaces are always stripped. 722 723 .. SeeAlso:: For editing a mdp file one can also use 724 :func:`gromacs.cbook.edit_mdp` (which works like a 725 poor replacement for sed). 726 """ 727 default_extension = "mdp" 728 logger = logging.getLogger('gromacs.formats.MDP') 729 730 COMMENT = re.compile("""\s*;\s*(?P<value>.*)""") # eat initial ws 731 # see regex in cbook.edit_mdp() 732 PARAMETER = re.compile(""" 733 \s*(?P<parameter>[^=]+?)\s*=\s* # parameter (ws-stripped), before '=' 734 (?P<value>[^;]*) # value (stop before comment=;) 735 (?P<comment>\s*;.*)? # optional comment 736 """, re.VERBOSE) 737

738 - def __init__(self, filename=None, autoconvert=True, **kwargs):

739 """Initialize mdp structure. 740 741 :Arguments: 742 *filename* 743 read from mdp file 744 *autoconvert* : boolean 745 ``True`` converts numerical values to python numerical types; 746 ``False`` keeps everything as strings [``True``] 747 *kwargs* 748 Populate the MDP with key=value pairs. (NO SANITY CHECKS; and also 749 does not work for keys that are not legal python variable names such 750 as anything that includes a minus '-' sign or starts with a number). 751 """ 752 super(MDP, self).__init__(**kwargs) # can use kwargs to set dict! (but no sanity checks!) 753 754 self.autoconvert = autoconvert 755 756 if not filename is None: 757 self._init_filename(filename) 758 self.read(filename)

759

760 - def _transform(self, value):

761 if self.autoconvert: 762 return autoconvert(value) 763 else: 764 return value

765

766 - def read(self, filename=None):

767 """Read and parse mdp file *filename*.""" 768 self._init_filename(filename) 769 770 def BLANK(i): 771 return "B%04d" % i

772 def COMMENT(i): 773 return "C%04d" % i

774 775 data = odict() 776 iblank = icomment = 0 777 with open(self.real_filename) as mdp: 778 for line in mdp: 779 line = line.strip() 780 if len(line) == 0: 781 iblank += 1 782 data[BLANK(iblank)] = '' 783 continue 784 m = self.COMMENT.match(line) 785 if m: 786 icomment += 1 787 data[COMMENT(icomment)] = m.group('value') 788 continue 789 # parameter 790 m = self.PARAMETER.match(line) 791 if m: 792 # check for comments after parameter?? -- currently discarded 793 parameter = m.group('parameter') 794 value = self._transform(m.group('value')) 795 data[parameter] = value 796 else: 797 errmsg = '%(filename)r: unknown line in mdp file, %(line)r' % vars() 798 self.logger.error(errmsg) 799 raise ParseError(errmsg) 800 801 super(MDP,self).update(data) 802 803

804 - def write(self, filename=None, skipempty=False):

805 """Write mdp file to *filename*. 806 807 :Keywords: 808 *filename* 809 output mdp file; default is the filename the mdp 810 was read from 811 *skipempty* : boolean 812 ``True`` removes any parameter lines from output that 813 contain empty values [``False``] 814 815 .. Note:: Overwrites the file that the mdp was read from if no 816 *filename* supplied. 817 """ 818 819 with open(self.filename(filename, ext='mdp'), 'w') as mdp: 820 for k,v in self.items(): 821 if k[0] == 'B': # blank line 822 mdp.write("\n") 823 elif k[0] == 'C': # comment 824 mdp.write("; %(v)s\n" % vars()) 825 else: # parameter = value 826 if skipempty and (v == '' or v is None): 827 continue 828 mdp.write("%(k)s = %(v)s\n" % vars())

829

830 831 -def autoconvert(s):

832 """Convert input to a numerical type if possible. 833 834 1. A non-string object is returned as it is 835 2. Try conversion to int, float, str. 836 """ 837 if not type(s) is str: 838 return s 839 for converter in int, float, str: # try them in increasing order of lenience 840 try: 841 return converter(s) 842 except ValueError: 843 pass 844 raise ValueError("Failed to autoconvert %r" % s)

845

Source Code for Module gromacs.formats