| Trees | Indices | Help |
|
|---|
|
|
1 # plugin: com.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 Centre of mass
8 ==============
9
10 Calculate the centre of mass of index groups.
11
12 Plugin class
13 ------------
14
15 .. autoclass:: COM
16 :members: worker_class
17 :undoc-members:
18
19
20 Worker class
21 ------------
22
23 The worker class performs the analysis.
24
25 .. autoclass:: _COM
26 :members:
27
28
29 """
30 from __future__ import with_statement
31
32 __docformat__ = "restructuredtext en"
33
34 import os.path
35 import warnings
36
37 import numpy
38
39 import gromacs
40 from gromacs.utilities import AttributeDict, asiterable
41 from gromacs.analysis.core import Worker, Plugin
42
43 import logging
44 logger = logging.getLogger('gromacs.analysis.plugins.com')
45
46 # Worker classes that are registered via Plugins (see below)
47 # ----------------------------------------------------------
48 # These must be defined before the plugins.
49
51 """COM worker class."""
52
54 """Set up COM analysis.
55
56 :Keywords:
57 *group_names*
58 list of index group names
59 *ndx*
60 index file if groups are not in the default index
61 *offset*
62 add the *offset* to the residue numbers [0]
63 *name*
64 plugin name [COM]
65 *simulation*
66 The :class:`gromacs.analysis.Simulation` instance that
67 owns the plugin [None]
68 """
69 group_names = asiterable(kwargs.pop('group_names', []))
70 ndx = kwargs.pop('ndx', None)
71 offset = kwargs.pop('offset', 0)
72
73 super(_COM, self).__init__(**kwargs)
74
75 self.parameters.group_names = group_names
76 self.parameters.offset = offset
77 self.ndx = ndx
78
79 if not self.simulation is None:
80 self._register_hook()
81
83 """Run when registering; requires simulation."""
84
85 super(_COM, self)._register_hook(**kwargs)
86 assert not self.simulation is None
87
88 if self.ndx is None:
89 self.ndx = self.simulation.ndx
90
91 self.parameters.filenames = { # result xvg files
92 'com': self.plugindir('com.xvg'),
93 }
94
95 # default filename for the plots -- not used
96 self.parameters.fignames = {
97 'com': self.figdir('com'),
98 }
99
101 """Analyze trajectory and write COM file.
102
103 All three components of the COM coordinate are written.
104
105 :Arguments:
106 - *force*: ``True`` does analysis and overwrites existing files
107 - *gmxargs*: additional keyword arguments for :func:`gromacs.g_bundle`
108 """
109 gmxargs['com'] = True
110 gmxargs['mol'] = False
111 gmxargs['ng'] = len(self.parameters.group_names)
112 gmxargs['x'] = True
113 gmxargs['y'] = True
114 gmxargs['z'] = True
115
116 if gmxargs['ng'] == 0:
117 errmsg = "No index group name(s) provided. Use group_name with the constructor."
118 logger.error(errmsg)
119 raise ValueError(errmsg)
120
121 if self.check_file_exists(self.parameters.filenames['com'], resolve='warning', force=force):
122 return
123
124 logger.info("Analyzing COM ...")
125 f = self.parameters.filenames
126 gromacs.g_traj(s=self.simulation.tpr, f=self.simulation.xtc, n=self.ndx,
127 ox=f['com'], input=self.parameters.group_names, **gmxargs)
128
129
131 """Collect output xvg files as :class:`gromacs.formats.XVG` objects.
132
133 - Make COM as a function of time available as XVG files and
134 objects.
135 - Compute RMSD of the COM of each group (from average
136 position, "rmsd").
137 - Compute distance whic encompasses 50% of observations ("median")
138 - Compute drift of COM, i.e. length of the vector between
139 initial and final position. Initial and final position are
140 computed as averages over *nframesavg* frames ("drift").
141
142 RMSD, median, and drift are columns in an xvg file. The rows correspond
143 to the groups in :attr:`gromacs.analysis.plugins.com.Worker.results.group_names`.
144
145 :Keywords:
146 *nframesavg*
147 number of initial and final frames that are averaged in
148 order to compute the drift of the COM of each group
149 [5000]
150 *refgroup*
151 group name whose com is taken as the reference and subtracted from
152 all other coms for the distance calculations. If supplied,
153 additional result 'com_relative_*refgroup*' is created.
154
155 :Returns: a dictionary of the results and also sets
156 :attr:`gromacs.analysis.plugins.com.Worker.results`.
157 """
158 from gromacs.formats import XVG
159
160 logger.info("Preparing COM graphs as XVG objects.")
161 self.results = AttributeDict( (k, XVG(fn)) for k,fn in self.parameters.filenames.items() )
162
163 # compute RMSD of COM and shift of COM (drift) between avg pos
164 # over first/last 5,000 frames
165 nframesavg = kwargs.pop('nframesavg', 5000)
166 ngroups = len(self.parameters.group_names)
167 xcom = self.results['com'].array
168
169 refgroup = kwargs.pop('refgroup', None)
170 if not refgroup is None:
171 if not refgroup in self.parameters.group_names:
172 errmsg = "refgroup=%s must be one of %r" % (refgroup, self.parameters.group_names)
173 logger.error(errmsg)
174 raise ValueError(errmsg)
175 nreference = 1 + 3 * self.parameters.group_names.index(refgroup) # 1-based !!
176 reference_com = xcom[nreference:nreference+3]
177 xcom[1:] -= numpy.vstack(ngroups * [reference_com]) # can't use broadcast
178 logger.debug("distances computed with refgroup %r", refgroup)
179
180 self.store_xvg('com_relative_%s' % refgroup, xcom,
181 names=['time']+self.parameters.group_names)
182
183
184 def vlength(v):
185 return numpy.sqrt(numpy.sum(v**2, axis=0)) # distances over time step
186
187 logger.debug("drift calculated between %d-frame averages at beginning and end",nframesavg)
188 records = []
189 for i in xrange(1, 3*ngroups+1, 3):
190 x = xcom[i:i+3]
191 r = vlength(x - x.mean(axis=1)[:,numpy.newaxis]) # distances over time step
192 #r0 = vlength(r - r[:,0][:,numpy.newaxis]) # distances over time step from r(t=0)
193 #h,edges = numpy.histogram(r, bins=kwargs.get('bins', 100), normed=True)
194 #m = 0.5*(edges[1:]+edges[:-1])
195 #c = h.cumsum(dtype=float) # integral
196 #c /= c[-1] # normalized (0 to 1)
197 #median = m[c < 0.5][-1]
198 #g = h/(4*numpy.pi*m**2)
199 #import scipy.integrate
200 #radint = lambda y: 4*numpy.pi*scipy.integrate.simps(m**2*y, x=m)
201 #g /= radint(g) # properly normalized radial distribution function
202 rmsd = numpy.sqrt(numpy.mean(r**2)) # radial spread sqrt(radint(m**2 * g))
203 median = numpy.median(r) # radius that contains 50% of the observations
204 dx = x[:,:nframesavg].mean(axis=1) - x[:,-nframesavg:].mean(axis=1)
205 drift = vlength(dx)
206 records.append((rmsd, median, drift))
207 self.store_xvg('distance', numpy.transpose(records), names="rmsd,median,drift")
208
209 return self.results
210
212 """Plot all results in one graph, labelled by the result keys.
213
214 :Keywords:
215 observables
216 select one or more of the stored results. Can be a list
217 or a string (a key into the results dict). ``None``
218 plots everything [``None``]
219 figure
220 - ``True``: save figures in the given formats
221 - "name.ext": save figure under this filename (``ext`` -> format)
222 - ``False``: only show on screen [``False``]
223 formats : sequence
224 sequence of all formats that should be saved [('png', 'pdf')]
225 plotargs
226 keyword arguments for pylab.plot()
227 """
228
229 import pylab
230 figure = kwargs.pop('figure', False)
231 observables = asiterable(kwargs.pop('observables', self.results.keys()))
232 extensions = kwargs.pop('formats', ('pdf','png'))
233
234 for name in observables:
235 result = self.results[name]
236 try:
237 result.plot(**kwargs) # This requires result classes with a plot() method!!
238 except AttributeError:
239 warnings.warn("Sorry, plotting of result %(name)r is not implemented" % vars(),
240 category=UserWarning)
241
242 # quick labels -- relies on the proper ordering
243 labels = [str(n)+" "+dim for n in self.parameters.group_names
244 for dim in 'xyz']
245 if not kwargs.get('columns', None) is None:
246 # select labels according to columns; only makes sense
247 # if plotting against the time (col 0)
248 if kwargs['columns'][0] == 0:
249 labels = numpy.array([None]+labels)[kwargs['columns'][1:]]
250 else:
251 labels = ()
252
253 pylab.legend(labels, loc='best')
254 if figure is True:
255 for ext in extensions:
256 self.savefig(ext=ext)
257 elif figure:
258 self.savefig(filename=figure)
259
260
261
262
263 # Public classes that register the worker classes
264 #------------------------------------------------
265
267 """*COM* plugin.
268
269 Calculate the centre of mass (COM) of various index groups.
270
271 .. class:: COM(group_names, [ndx[, offset [, name[, simulation]]]])
272
273 """
274 worker_class = _COM
275
| Trees | Indices | Help |
|
|---|
| Generated by Epydoc 3.0.1 on Sat Jun 12 15:59:40 2010 | http://epydoc.sourceforge.net |