gromacs.analysis.plugins.gridmatmd

1 # $Id$ 2 # Copyright (c) 2009 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.analysis.plugins.gridmatmd` --- Lipid bilayer analysis helper 8 =========================================================================== 9 10 This helper module contains code to drive ``GridMAT-MD.pl``, available 11 from the `GridMAT-MD`_ home page and written by WJ Allen et al 12 [Allen2009]_ . The GromacsWrapper distribution comes with version 13 1.0.2 of ``GridMAT-MD.pl`` and includes a small patch so that it can 14 accept filenames on the command line. 15 16 17 18 References 19 ---------- 20 21 .. [Allen2009] W. J. Allen, J. A. Lemkul, and D. R. Bevan. (2009) "GridMAT-MD: A 22 Grid-based Membrane Analysis Tool for Use With Molecular Dynamics." 23 J. Comput. Chem. 30 (12): 1952-1958. 24 25 .. _GridMAT-MD: http://www.bevanlab.biochem.vt.edu/GridMAT-MD/index.html 26 27 Module contents 28 --------------- 29 30 .. autoclass:: GridMatMD 31 :members: 32 33 .. autoclass:: GridMatData 34 :members: 35 :inherited-members: 36 37 .. autoclass:: Grid2D 38 :members: imshow, _midpoints, __add__, __sub__, __mul__, __div__ 39 40 """ 41 from __future__ import with_statement 42 43 import sys 44 import subprocess 45 import re 46 import glob 47 import cPickle 48 import numpy 49 50 import gromacs.tools 51

52 -class GridMatMD(object):

53 """Analysis of lipid bilayers with GridMAT-MD. 54 55 It requires a configuration file and a list of structure files 56 (gro or pdb) as input. See the documentation (pdf_) for the format of the 57 config file. Note that the *bilayer* keyword will be ignored in 58 the config file. 59 60 .. _pdf: http://www.bevanlab.biochem.vt.edu/GridMAT-MD/doc/GridMAT-MD_ug_v1.0.2.pdf 61 """ 62 63 # regular expression to pull useful information from output 64 pNX = re.compile("There are (?P<NX>\d+) grid points in the X direction, " 65 "spaced every (?P<DX>[0-9.]+) nanometers") 66 pNY = re.compile("There are (?P<NY>\d+) grid points in the Y direction, " 67 "spaced every (?P<DY>[0-9.]+) nanometers") 68 pTOP = re.compile('The top leaflet "thickness" will be printed to (?P<FILENAME>.*_top_.*\.dat)') 69 pAVG = re.compile('The average bilayer "thickness" will be printed to (?P<FILENAME>.*_average_.*\.dat)') 70 pBOT = re.compile('The bottom leaflet "thickness" will be printed to (?P<FILENAME>.*_bottom_.*\.dat)') 71 72

73 - def __init__(self, config, filenames):

74 """Set up GridMAT-MD analysis. 75 76 :Arguments: 77 *config* : filename 78 input file for GridMAT-MD (see docs) 79 *filenames* : list or glob-pattern 80 list of gro or pdb files, or a glob pattern that creates 81 such a list 82 """ 83 self.config = config 84 if type(filenames) is str: 85 # use it as a glob pattern 86 _filenames = glob.glob(filenames) 87 else: 88 _filenames = filenames 89 #: list of filenames (gro or pdb files) 90 self.filenames = _filenames 91 #: Number of bins in the *x* direction. 92 self.nx = None 93 #: Numberof bins in the *y* direction. 94 self.ny = None 95 #: Bin width in *x*. 96 self.dx = None 97 #: Bin width in *y*. 98 self.dy = None 99 #: Thickness results (dict of :class:`GridMatData` instances); 100 #: used to store the results from the last processed frame. 101 self.thickness = {} 102 #: Final result: averaged thickness profiles, indexed by "top", "bottom", "average" 103 self.averages = {} 104 105 #: Instance of :class:`gromacs.tools.GridMAT_MD` which always 106 #: uses *config* for input. 107 self.GridMatMD = gromacs.tools.GridMAT_MD(self.config)

108

109 - def run_frame(self, frame):

110 """Run GridMAT-MD on a single *frame* and return results. 111 112 :Arguments: *frame* is a filename (gro or pdb) 113 :Returns: a dict of :class:`GridMapData` objects; the keys 114 are "top", "bottom", "average" 115 """ 116 117 thickness = {} 118 print >> sys.stdout, "Analyzing frame %r (takes a few seconds)..." % frame 119 sys.stdout.flush() 120 rc,output,stderr = self.GridMatMD(frame, stdout=False) 121 for line in output.split('\n'): 122 print >> sys.stdout, '>>> '+line 123 m = self.pNX.match(line) 124 if m: 125 self.nx, self.dx = int(m.group('NX')), float(m.group('DX')) 126 continue 127 m = self.pNY.match(line) 128 if m: 129 self.ny, self.dy = int(m.group('NY')), float(m.group('DY')) 130 continue 131 m = self.pTOP.match(line) 132 if m: 133 thickness['top'] = m.group('FILENAME') 134 continue 135 m = self.pBOT.match(line) 136 if m: 137 thickness['bottom'] = m.group('FILENAME') 138 continue 139 m = self.pAVG.match(line) 140 if m: 141 thickness['average'] = m.group('FILENAME') 142 continue 143 144 d = {} # thickness results 145 for name, filename in thickness.items(): 146 d[name] = GridMatData(filename, shape=(self.nx, self.ny), 147 delta=(self.dx, self.dy)) 148 self.thickness = d # use thickness as common block... useful for interactive work 149 return d

150

151 - def run(self):

152 """Run analysis on all files and average results.""" 153 154 results = {} 155 for f in self.filenames: 156 d = self.run_frame(f) # <-- run GridMAT-MD 157 for name,grid2d in d.items(): 158 try: 159 results[name].append(grid2d) 160 except KeyError: 161 results[name] = [grid2d] 162 163 # averages arrays and bins 164 self.averages = dict([(name,numpy.mean(v)) for name,v in results.items()])

165

166 - def imshow(self, name, **kwargs):

167 """Display array *name* with ``pylab.imshow``.""" 168 from pylab import xlabel, ylabel 169 try: 170 self.averages[name].imshow(**kwargs) 171 except KeyError: 172 raise KeyError("name must be one of %r" % self.arrays.keys()) 173 xlabel(r'$x/$nm') 174 ylabel(r'$y/$nm')

175

176 - def save(self, name):

177 """Save object as pickle.""" 178 with open(name, 'wb') as f: 179 cPickle.dump(self, f, protocol=cPickle.HIGHEST_PROTOCOL)

180

181 -class Grid2D(object):

182 """Represents a 2D array with bin sizes. 183 184 Addition and subtraction of grids is defined for the arrays and 185 the bins. Multiplication and division with scalars is also 186 defined. Each operation returns a new :class:`Grid2D` object. 187 188 (Actually, it should work for arrays of any dimension, not just 2D.) 189 """ 190

191 - def __init__(self, data, bins):

192 """Initialize the Grid2D instance. 193 194 :Arguments: 195 *data* 196 array data, e.g. a list of array 197 *bins* 198 tuple of lists of bin **edges**, one for each dimension 199 """ 200 self.array = numpy.asarray(data) # numpy array 201 self.shape = self.array.shape 202 self.ndim = len(self.shape) 203 self.bins = bins 204 if len(self.bins) != self.ndim: 205 raise ValueError("bins must be a list of bin edges, one for each dimension, D=%(ndim)d" % 206 vars(self)) 207 self.midpoints = [self._midpoints(x) for x in self.bins]

208

209 - def _midpoints(self, x):

210 _x = numpy.asarray(x) 211 return 0.5*(_x[1:] + _x[:-1])

212

213 - def imshow(self, **kwargs):

214 """Display data as a 2D image using :func:`pylab.imshow`.""" 215 import pylab 216 # extent: [ None | scalars (left, right, bottom, top) ] 217 extent = numpy.concatenate([self.bins[0][[0,-1]], self.bins[1][[0,-1]]]) 218 kwargs.setdefault('extent', extent) 219 kwargs['origin'] = 'lower' 220 pylab.imshow(self.array, **kwargs)

221

222 - def __add__(self, other):

223 """Add arrays and bins (really only makes sense when averaging).""" 224 if self.array.shape != other.array.shape: 225 raise TypeError("arrays are of incompatible shape") 226 _bins = [self.bins[dim] + other.bins[dim] for dim in xrange(len(self.bins))] 227 _array = self.array + other.array 228 return Grid2D(data=_array, bins=_bins)

229

230 - def __sub__(self, other):

231 """Subtract other from self (also subtracts bins... which is odd but consistent).""" 232 if self.array.shape != other.array.shape: 233 raise TypeError("arrays are of incompatible shape") 234 _bins = [self.bins[dim] - other.bins[dim] for dim in xrange(len(self.bins))] 235 _array = self.array - other.array 236 return Grid2D(data=_array, bins=_bins)

237

238 - def __mul__(self, x):

239 """Multiply arrays (and bins) by a scalar *x*.""" 240 _bins = [self.bins[dim] * x for dim in xrange(len(self.bins))] 241 _array = self.array * x 242 return Grid2D(data=_array, bins=_bins)

243

244 - def __div__(self, x):

245 """Divide arrays (and bins) by a scalar *x*.""" 246 _bins = [self.bins[dim] / x for dim in xrange(len(self.bins))] 247 _array = self.array / x 248 return Grid2D(data=_array, bins=_bins)

249 250 251

252 -class GridMatData(Grid2D):

253 """Represent GridMatMD data file. 254 255 The loaded array data is accessible as a numpy array in 256 :attr:`GridMatData.array` and bins and midpoints as 257 :attr:`GridMatData.bins` and :attr:`GridMatData.midpoints` 258 respectively. 259 """ 260 261 DATANAME = re.compile(""" 262 (?P<NX>\d+) # number of grid points in X 263 x 264 (?P<NY>\d+) # number of grid points in Y 265 _(top|bottom|average)_.*\.dat # all the rest 266 """, re.VERBOSE) 267 268

269 - def __init__(self, filename, shape=None, delta=None):

270 """Load the data into a numpy array. 271 272 The *filename* is an output file from GridMAT-MD. *shape* and 273 *delta* are optional. The *shape* of the array is parsed from 274 the filename if not provided. The spacing is set to (1,1) if 275 not provided. 276 277 :Arguments: 278 *filename* 279 2D grid as written by GridMAT-MD 280 *shape* 281 Shape tuple (NX, NY) of the array in filename. 282 *delta* 283 Tuple of bin sizes of grid (DX, DY). 284 """ 285 self.filename = filename 286 self.shape = shape 287 if shape is None: 288 self.shape = self.parse_filename(filename) 289 if delta is None: 290 self.delta = tuple([1.0] * len(self.shape)) # arbitrary spacing 291 else: 292 self.delta = tuple(delta) 293 294 #: 2D data from file is made available as a numpy array. 295 _array = numpy.loadtxt(self.filename).reshape(self.shape) 296 297 #: reconstructed bins (always start at 0,0 as the offset is lost) 298 _bins = [numpy.linspace(0,n*dx,n+1) for n,dx in zip(self.shape, self.delta)] 299 300 Grid2D.__init__(self, data=_array, bins=_bins)

301 302 # reconstructed midpoints (always start at 0,0 as the offset 303 # is lost) 304

305 - def parse_filename(self, filename):

306 """Get dimensions from filename""" 307 m = self.DATANAME.match(filename) 308 if m is None: 309 raise ValueError("filename %s does not appear to be a standard GridMat-MD output filename") 310 return map(int, m.group('NX', 'NY'))

311

Source Code for Module gromacs.analysis.plugins.gridmatmd