Package gromacs :: Package analysis :: Package plugins :: Module helixbundle
[hide private]
[frames] | no frames]

Source Code for Module gromacs.analysis.plugins.helixbundle

  1  # plugin: helixbundle.py 
  2  # Copyright (c) 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  Helix bundle analysis 
  8  ===================== 
  9   
 10  Analysis of helix bundles with :func:`gromacs.g_bundle`. 
 11   
 12  .. SeeAlso:: HELANAL does some things better than g_bundle, in 
 13                particular it can find kinks automatically. 
 14   
 15  Plugin class 
 16  ------------ 
 17   
 18  .. autoclass:: HelixBundle 
 19     :members: worker_class 
 20     :undoc-members: 
 21   
 22  Example 
 23  ~~~~~~~ 
 24   
 25  .. Warning:: Note that the helices should be sorted in ascending order 
 26               because all indices are sorted that way. If they are in 
 27               mixed order then the correspondence between the name 
 28               column and the helices is lost. This is a shortcoming of 
 29               :program:`make_ndx` that is not easy to rectify. 
 30   
 31  Define helices in a reST table (use column names as shown!), set up 
 32  the plugin, and add it to a simulation instance:: 
 33   
 34     helixtable = ''' 
 35          Table[helices]: helices top = N-term, bot = C-term; kink 
 36          =====  ================ ================ =============== 
 37          name   top              bottom           kink 
 38          =====  ================ ================ =============== 
 39          TM1    G28  P29  F30    T53  S54  S55    I41  Q42  V43 
 40          TM2    Q57  V58  W59    I84  R85  W86    F78  T79  Q80 
 41          TM3    R101 G102 S103   L135 L136 T137   F116 W117 F118 
 42          TM4    L142 P143 L144   T157 F158 Y159   F149 G150 A151 
 43          TM6    F209 S210 T211   D229 I230 V231   G219 W220 I221 
 44          TM7    E243 G244 Q245   L277 V278 G279   V261 P262 A263 
 45          TM8    P297 M298 A299   C327 S328 T329   N314 P315 A316 
 46          TM9    F336 K337 T338   L348 L349 M350   V342 S343 A344 
 47          =====  ================ ================ =============== 
 48          ''' 
 49     hb = HelixBundle(helixtable=helixtable, offset=-9) 
 50     S = Simulation(....) 
 51     S.add_plugin(hb) 
 52   
 53   
 54  Worker class 
 55  ------------ 
 56   
 57  The worker class performs the analysis. 
 58   
 59  .. autoclass:: _HelixBundle 
 60     :members: 
 61   
 62   
 63  """ 
 64  from __future__ import with_statement 
 65   
 66  __docformat__ = "restructuredtext en" 
 67   
 68  import os.path 
 69  import warnings 
 70   
 71  import recsql 
 72   
 73  import gromacs 
 74  from gromacs.utilities import AttributeDict 
 75  from gromacs.analysis.core import Worker, Plugin 
 76   
 77  import logging 
 78  logger = logging.getLogger('gromacs.analysis.plugins.helixbundle') 
 79   
 80  # Worker classes that are registered via Plugins (see below) 
 81  # ---------------------------------------------------------- 
 82  # These must be defined before the plugins. 
 83   
84 -class _HelixBundle(Worker):
85 """HelixBundle worker class.""" 86
87 - def __init__(self,**kwargs):
88 """Set up HelixBundle analysis. 89 90 :Keywords: 91 *helixtable* 92 reST table with columns "name", "top", "bottom", 93 "kink"; see :mod:`gromacs.analysis.plugins.helixbundle` 94 for details 95 *offset* 96 add the *offset* to the residue numbers in *helixtable* [0] 97 *helixndx* 98 provide a index file with the appropriate groups 99 instead of the table; also requires *na* 100 *na* 101 number of helices 102 *with_kinks* 103 take kinks into account [True] 104 *name* 105 plugin name [HelixBundle] 106 *simulation* 107 The :class:`gromacs.analysis.Simulation` instance that 108 owns the plugin [None] 109 """ 110 helixndx = kwargs.pop('helixndx', None) 111 helixtable = kwargs.pop('helixtable', None) 112 na = kwargs.pop('na', None) 113 offset = kwargs.pop('offset', 0) 114 with_kinks = kwargs.pop('with_kinks', True) 115 116 if helixndx is None: 117 if helixtable is None: 118 raise ValueError("HelixBundle: requires a table *helixtable* with residues " 119 "that indicate the helix boundaries and kink or a index " 120 "file *helixndx* appropriate for g_bundle. See the " 121 "docs for details.") 122 self.helices = recsql.rest_table.Table2array(helixtable, autoconvert=True).recarray() 123 na = len(self.helices) 124 elif na is None: 125 raise ValueError("When using *helixndx* one must also provide the total number of " 126 "helices *na*. See g_bundle docs for details.") 127 128 super(_HelixBundle, self).__init__(**kwargs) 129 130 self.parameters.na = na 131 self.parameters.offset = offset 132 self.parameters.with_kinks = with_kinks 133 134 if not self.simulation is None: 135 self._register_hook()
136
137 - def _register_hook(self, **kwargs):
138 """Run when registering; requires simulation.""" 139 140 super(_HelixBundle, self)._register_hook(**kwargs) 141 assert not self.simulation is None 142 143 self.helixndx = self.plugindir('helices.ndx') # special index for g_bundle 144 145 self.parameters.filenames = { # result xvg files 146 'length': self.plugindir('len.xvg'), 147 'distance': self.plugindir('dist.xvg'), 148 'z': self.plugindir('z.xvg'), 149 'tilt': self.plugindir('tilt.xvg'), 150 'tilt_radial': self.plugindir('tilt_radial.xvg'), 151 'tilt_lateral': self.plugindir('tilt_lateral.xvg'), 152 } 153 if self.parameters.with_kinks: 154 self.parameters.filenames.update( 155 {'kink': self.plugindir('kink.xvg'), 156 'kink_radial': self.plugindir('kink_radial.xvg'), 157 'kink_lateral': self.plugindir('kink_lateral.xvg'), 158 }) 159 160 # default filename for the plots -- not used 161 self.parameters.fignames = { 162 'z': self.figdir('z'), 163 'tilt': self.figdir('tilt'), 164 'kink': self.figdir('kink'), 165 }
166
167 - def make_index(self, force=None):
168 """Build g_bundle index file from a record array. 169 170 :Keywords: 171 *force* 172 - ``None`` does nothing if the index file already exists 173 - ``False`` raises exception if it exists 174 - ``True`` always rebuild index file 175 176 Format of the index table: 177 - must contain columns name, top, bottom, kink 178 - a row corresponds to one helis 179 - name contains the name of the helix e.g. 'TM5' 180 - a entry is a *string* of residues which is split on white space; 181 entries in the same column must have the same number of residues 182 (limitation of g_bundle) 183 184 """ 185 186 if self.check_file_exists(self.helixndx, resolve="indicate", force=force): 187 logger.warning("make_index(): The index file %r already exists; " 188 "use force = True to rebuild", self.helixndx) 189 return self.helixndx 190 191 logger.info("Making index file for %d helices %r", len(self.helices), self.helixndx) 192 logger.debug("make_index: tpr = %r", self.simulation.tpr) 193 logger.debug("make_index: offset = %r", self.parameters.offset) 194 195 # quick'n'crappy... 196 tops = [] 197 bottoms = [] 198 kinks = [] 199 200 for h in self.helices: 201 # split entries on white space (Table2array() does not handle it properly yet) 202 tops.extend(h.top.split()) 203 bottoms.extend(h.bottom.split()) 204 kinks.extend(h.kink.split()) 205 206 tpr = self.simulation.tpr 207 offset = self.parameters.offset 208 209 # TODO: use tmp files for indices 210 211 Itop = gromacs.cbook.IndexBuilder(tpr, tops, offset=offset, out_ndx='top.ndx') 212 name, topndx = Itop.combine(name_all='top', defaultgroups=False) 213 214 Ibot = gromacs.cbook.IndexBuilder(tpr, bottoms, offset=offset, out_ndx='bottom.ndx') 215 name, botndx = Ibot.combine(name_all='bottom', defaultgroups=False) 216 217 Ikink = gromacs.cbook.IndexBuilder(tpr, kinks, offset=offset, out_ndx='kinks.ndx') 218 name, kinkndx = Ikink.combine(name_all='kink', defaultgroups=False) 219 220 Iall = gromacs.cbook.IndexBuilder(tpr, ndx=[topndx, botndx, kinkndx]) 221 Iall.cat(self.helixndx) 222 223 return self.helixndx
224 225 226
227 - def run(self, force=None, **gmxargs):
228 """Analyze trajectory and write HelixBundle files. 229 230 :Arguments: 231 - *force*: ``True`` does analysis and overwrites existing files 232 - *gmxargs*: additional keyword arguments for :func:`gromacs.g_bundle` 233 234 .. Note:: The plugin default is *z* = ``True``, i.e. the tilt is computed 235 relative to the box z-axis. 236 """ 237 gmxargs.setdefault('z', True) 238 gmxargs['na'] = self.parameters.na 239 240 if self.check_file_exists(self.parameters.filenames['tilt'], resolve='warning', force=force): 241 return 242 243 self.make_index(force=force) 244 245 logger.info("Analyzing HelixBundle (with_kinks=%r)...", self.parameters.with_kinks) 246 f = self.parameters.filenames 247 if self.parameters.with_kinks: 248 gromacs.g_bundle(s=self.simulation.tpr, f=self.simulation.xtc, n=self.helixndx, 249 ol=f['length'], od=f['distance'], oz=f['z'], 250 ot=f['tilt'], otr=f['tilt_radial'], otl=f['tilt_lateral'], 251 ok=f['kink'], okr=f['kink_radial'], okl=f['kink_lateral'], 252 input=['top','bottom', 'kink'], 253 **gmxargs) 254 else: 255 gromacs.g_bundle(s=self.simulation.tpr, f=self.simulation.xtc, n=self.helixndx, 256 ol=f['length'], od=f['distance'], oz=f['z'], 257 ot=f['tilt'], otr=f['tilt_radial'], otl=f['tilt_lateral'], 258 input=['top','bottom'], 259 **gmxargs)
260 261
262 - def analyze(self,**kwargs):
263 """Collect output xvg files as :class:`gromacs.formats.XVG` objects. 264 265 :Returns: a dictionary of the results and also sets ``self.results``. 266 """ 267 from gromacs.formats import XVG 268 269 logger.info("Preparing HelixBundle graphs as XVG objects.") 270 results = AttributeDict( (k, XVG(fn)) for k,fn in self.parameters.filenames.items() ) 271 self.results = results 272 return results
273
274 - def plot(self, **kwargs):
275 """Plot all results in one graph, labelled by the result keys. 276 277 :Keywords: 278 figure 279 - ``True``: save figures in the given formats 280 - "name.ext": save figure under this filename (``ext`` -> format) 281 - ``False``: only show on screen 282 formats : sequence 283 sequence of all formats that should be saved [('png', 'pdf')] 284 plotargs 285 keyword arguments for pylab.plot() 286 """ 287 288 import pylab 289 figure = kwargs.pop('figure', False) 290 extensions = kwargs.pop('formats', ('pdf','png')) 291 for name,result in self.results.items(): 292 kwargs['label'] = name 293 try: 294 result.plot(**kwargs) # This requires result classes with a plot() method!! 295 except AttributeError: 296 warnings.warn("Sorry, plotting of result %(name)r is not implemented" % vars(), 297 category=UserWarning) 298 pylab.legend(loc='best') 299 if figure is True: 300 for ext in extensions: 301 self.savefig(ext=ext) 302 elif figure: 303 self.savefig(filename=figure)
304 305 306 307 308 # Public classes that register the worker classes 309 #------------------------------------------------ 310
311 -class HelixBundle(Plugin):
312 """*HelixBundle* plugin. 313 314 :func:`gromacs.g_bundle` helix analysis 315 316 .. class:: HelixBundle([helixtable, offset, with_kinks, [name[, simulation]]]) 317 318 :Arguments: 319 *helixtable* 320 reST table with columns "name", "top", "bottom", 321 "kink"; see :mod:`gromacs.analysis.plugins.helixbundle` 322 for details 323 *offset* 324 add the *offset* to the residue numbers in *helixtable* [0] 325 *helixndx* 326 provide a index file with the appropriate groups 327 instead of the table; also requires *na* 328 *na* 329 number of helices 330 *with_kinks* 331 take kinks into account [True] 332 *name* 333 plugin name [HelixBundle] 334 *simulation* 335 The :class:`gromacs.analysis.Simulation` instance that owns 336 the plugin. [None] 337 """ 338 worker_class = _HelixBundle
339