gromacs.analysis.plugins.dihedrals

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 Dihedrals plugin 8 ================ 9 10 Analyze one or more dihedral angles. 11 12 Plugin class 13 ------------ 14 15 .. autoclass:: Dihedrals 16 :members: worker_class 17 :undoc-members: 18 19 Worker class 20 ------------ 21 22 The worker class performs the analysis. 23 24 .. autoclass:: _Dihedrals 25 :members: 26 27 28 """ 29 from __future__ import with_statement 30 31 __docformat__ = "restructuredtext en" 32 33 import os.path 34 import warnings 35 36 import numpy 37 38 import gromacs 39 from gromacs import MissingDataError 40 from gromacs.utilities import AttributeDict 41 from gromacs.analysis.core import Worker, Plugin 42 from gromacs.formats import XVG, NDX 43 44 45 # Worker classes that are registered via Plugins (see below) 46 # ---------------------------------------------------------- 47 # These must be defined before the plugins. 48

49 -class _Dihedrals(Worker):

50 """Dihedrals worker class.""" 51

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

53 """Set up dihedral analysis. 54 55 :Arguments: 56 *dihedrals* 57 list of tuples; each tuple contains :class:`gromacs.cbook.IndexBuilder` 58 atom selection commands. 59 *labels* 60 optional list of labels for the dihedrals. Must have as many entries as 61 *dihedrals*. 62 63 """ 64 # specific arguments: take them before calling the super class that 65 # does not know what to do with them 66 dihedralgroups = kwargs.pop('dihedrals', []) 67 for dih in dihedralgroups: 68 if len(dih) != 4 or type(dih) is str: 69 raise ValueError("Each dihedral index group must contain exactly four atomnumbers, " 70 "but this is not the case for %r." % dih) 71 labels = kwargs.pop('labels', None) 72 if not labels is None: 73 if len(labels) != len(dihedralgroups): 74 raise ValueError("Provide one label in labels for each dihedral in dihedrals.") 75 76 # super class init: do this before doing anything else 77 # (also sets up self.parameters and self.results) 78 super(_Dihedrals, self).__init__(**kwargs) 79 80 # process specific parameters now and set instance variables 81 self.parameters.dihedrals = dihedralgroups 82 self.parameters.labels = labels 83 84 # self.simulation might have been set by the super class 85 # already; just leave this snippet at the end. Do all 86 # initialization that requires the simulation class in the 87 # _register_hook() method. 88 if not self.simulation is None: 89 self._register_hook()

90

91 - def _register_hook(self, **kwargs):

92 """Run when registering; requires simulation.""" 93 94 super(_Dihedrals, self)._register_hook(**kwargs) 95 assert not self.simulation is None 96 97 def P(*args): 98 return self.plugindir(*args)

99 100 self.parameters.ndx = P('dih.ndx') 101 self.parameters.filenames = { 102 'comb_distribution': P('dih_combdist.xvg'), # combined distribution from ALL dih! 103 'timeseries': P('dih_timeseries.xvg'), # with -all, ie (phi, avg, dh1, dh2, ...) 104 'transitions': P('dih_trans.xvg'), # only works for multiplicity 3 105 'transition_histogram': P('dih_transhisto.xvg'), 106 'acf': P('dih_acf.xvg'), 107 'distributions' : P('dih_distributions.xvg'), # from analyze() 108 'PMF': P('dih_pmf.xvg'), # from analyze() 109 } 110 self.parameters.figname = self.figdir('dihedrals_PMF') 111 112 # just make index now... we have everything and it's really simple 113 self.make_index()

114

115 - def make_index(self):

116 """Make one index group *dihedrals* from the selections. 117 118 Right now this requires **raw atom numbers** from the topology. 119 """ 120 ndx = NDX() 121 ndx['dihedrals'] = numpy.concatenate(self.parameters.dihedrals) 122 ndx.write(filename=self.parameters.ndx, ncol=4)

123 124 125 # override 'API' methods of base class 126

127 - def run(self, force=False, **gmxargs):

128 """Collect dihedral data from trajectory with ``g_angle`` and save to data files. 129 130 Note that the ``-all`` flag is always enabled and hence all time series 131 can be found under the *timeseries* results key. 132 """ 133 if not force and \ 134 self.check_file_exists(self.parameters.filenames['timeseries'], resolve='warn'): 135 return 136 137 angle_type = gmxargs.setdefault('type','dihedral') 138 allowed_types = ("dihedral", "improper", "ryckaert-bellemans") 139 if angle_type not in allowed_types: 140 raise ValueError("The *type* can not be %r, only be one of\n\t%r\n" % 141 (angle_type, allowed_types)) 142 gmxargs['all'] = True # required because we analyze the time series file 143 gmxargs.setdefault('periodic', True) 144 F = self.parameters.filenames 145 gromacs.g_angle(f=self.simulation.xtc, n=self.parameters.ndx, 146 od=F['comb_distribution'], ov=F['timeseries'], ot=F['transitions'], 147 oc=F['acf'], oh=F['transition_histogram'], **gmxargs)

148 149

150 - def analyze(self, **kwargs):

151 """Load results from disk into :attr:`_Dihedrals.results` and compute PMF. 152 153 The PMF W(phi) in kT is computed from each dihedral 154 probability distribution P(phi) as 155 156 W(phi) = -kT ln P(phi) 157 158 It is stored in :attr:`_Dihedrals.results` with the key *PMF*. 159 160 :Keywords: 161 *bins* 162 bins for histograms (passed to numpy.histogram(new=True)) 163 164 :Returns: a dictionary of the results and also sets 165 :attr:`_Dihedrals.results`. 166 """ 167 168 bins = kwargs.pop('bins', 361) 169 170 results = AttributeDict() 171 172 # get graphs that were produced by g_angle 173 for name, f in self.parameters.filenames.items(): 174 try: 175 results[name] = XVG(f) 176 except IOError: 177 pass # either not computed (yet) or some failure 178 179 # compute individual distributions 180 ts = results['timeseries'].array # ts[0] = time, ts[1] = avg 181 dih = ts[2:] 182 183 phi_range = (-180., 180.) 184 185 Ndih = len(dih) 186 p = Ndih * [None] # histograms (prob. distributions), one for each dihedral i 187 for i in xrange(Ndih): 188 phis = dih[i] 189 p[i],e = numpy.histogram(phis, bins=bins, range=phi_range, normed=True, new=True) 190 191 P = numpy.array(p) 192 phi = 0.5*(e[:-1]+e[1:]) # midpoints of bin edges 193 distributions = numpy.concatenate((phi[numpy.newaxis, :], P)) # phi, P[0], P[1], ... 194 195 xvg = XVG() 196 xvg.set(distributions) 197 xvg.write(self.parameters.filenames['distributions']) 198 results['distributions'] = xvg 199 del xvg 200 201 # compute PMF (from individual distributions) 202 W = -numpy.log(P) # W(phi)/kT = -ln P 203 W -= W.min(axis=1)[:, numpy.newaxis] # minimum at 0 kT 204 pmf = numpy.concatenate((phi[numpy.newaxis, :], W), axis=0) 205 xvg = XVG() 206 xvg.set(pmf) 207 xvg.write(self.parameters.filenames['PMF']) 208 results['PMF'] = xvg 209 210 self.results = results 211 return results

212

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

214 """Plot all results in one graph, labelled by the result keys. 215 216 :Keywords: 217 figure 218 - ``True``: save figures in the given formats 219 - "name.ext": save figure under this filename (``ext`` -> format) 220 - ``False``: only show on screen 221 formats : sequence 222 sequence of all formats that should be saved [('png', 'pdf')] 223 with_legend : bool 224 add legend from *labels* to graphs 225 plotargs 226 keyword arguments for pylab.plot() 227 """ 228 229 from pylab import subplot, xlabel, ylabel, legend 230 figure = kwargs.pop('figure', False) 231 extensions = kwargs.pop('formats', ('pdf','png')) 232 233 kwargs.setdefault('linestyle', '-') 234 kwargs.setdefault('linewidth', 3) 235 with_legend = kwargs.pop('with_legend', True) 236 237 # data 238 try: 239 P = self.results['distributions'] 240 W = self.results['PMF'] 241 except KeyError: 242 raise MissingDataError("No post-processed data available: " 243 "Run the analysis() method first.") 244 245 # plot probability and PMF 246 subplot(211) 247 P.plot(**kwargs) 248 subplot(212) 249 W.plot(**kwargs) 250 xlabel('dihedral angle $\phi/\degree$') 251 ylabel(r'PMF $\mathcal{W}/kT$') 252 253 # use legend labels from init 254 if with_legend and len(self.parameters.labels) > 0: 255 if 'columns' in kwargs: 256 # sneaky option for XVG.plot: selects columns to plot, eg [0,1, 3] so we 257 # need to get the correct labels, namely [0, 2] 258 c = numpy.asarray(kwargs['columns']) 259 not_zero = (c != 0) 260 c = c[not_zero] # filter 0th column (phi) which does not have a label 261 c -= 1 # it's now an index into labels 262 if numpy.all(not_zero): 263 # special case when 0th column was not included (eg correlation plots): 264 # we have to discard the first label (and *should* change the xlabel...) 265 xlabel(self.parameters.labels[c[0]]) 266 c = c[1:] 267 _labels = tuple([self.parameters.labels[i] for i in c]) 268 else: 269 _labels = tuple(self.parameters.labels) 270 subplot(211) 271 legend(_labels, loc='best') 272 273 # write graphics file(s) 274 if figure is True: 275 for ext in extensions: 276 self.savefig(ext=ext) 277 elif figure: 278 self.savefig(filename=figure)

279 280 281 282 # Public classes that register the worker classes 283 #------------------------------------------------ 284

285 -class Dihedrals(Plugin):

286 """*Dihedrals* plugin. 287 288 .. class:: Dihedrals(dihedrals[, labels[,name[, simulation]]]) 289 290 :Keywords: 291 *dihedrals* 292 list of tuples; each tuple contains atom indices that define the dihedral. 293 *labels* 294 optional list of labels for the dihedrals. Must have as many entries as 295 *dihedrals*. 296 *name* : string 297 plugin name (used to access it) 298 *simulation* : instance 299 The :class:`gromacs.analysis.Simulation` instance that owns the plugin. 300 301 """ 302 worker_class = _Dihedrals

303

Source Code for Module gromacs.analysis.plugins.dihedrals