1 """Collections of classes for the determination of dynamics-related properties.
2
3 Classes:
4
5 * MeanSquareDisplacement : sets up a Mean-Square-Displacement analysis.
6 * RootMeanSquareDeviation : sets up a Root Mean-Square-Deviation analysis.
7 * GyrationRadius : sets up a Gyration Radius analysis.
8 * AngularCorrelation : sets up an Angular Correlation analysis.
9 * CartesianVelocityAutoCorrelationFunction : sets up a Cartesian Velocity AutoCorrelation analysis.
10 * DensityOfStates : sets up a Density Of States analysis.
11 * AutoRegressiveAnalysis : sets up an Auto-Regressive analysis.
12 * QuasiHarmonicAnalysis : sets up a Quasi-Harmonic analysis.
13 * PassBandTrajectoryFilter : sets up a Pass-Band Trajectory Filter.
14 * GlobalMotionTrajectoryFilter : sets up a Global Motion Trajectory Filter.
15 * CenterOfMassTrajectory : sets up a Center Of Mass Trajectory.
16 * RigidBodyTrajectory : sets up a Rigid-Body Trajectory.
17 * AngularVelocityAutoCorrelationFunction : sets up an Angular Velocity AutoCorrelation Function.
18 * AngularDensityOfStates : sets up an Angular Density Of States.
19 """
20
21
22 import copy
23 import operator
24 import os
25 import sys
26 from time import asctime
27 from timeit import default_timer
28
29
30 from Scientific import N
31 from Scientific.Functions.Interpolation import InterpolatingFunction
32 from Scientific.Geometry import Vector
33 from Scientific.Geometry.Quaternion import Quaternion
34 from Scientific.Geometry.Transformation import Translation
35 from Scientific.IO.NetCDF import NetCDFFile
36 from Scientific.LA import Heigenvectors
37 from Scientific.Signals.Models import AutoRegressiveModel, AveragedAutoRegressiveModel
38
39
40 from MMTK import Atom, AtomCluster
41 from MMTK import Units
42 from MMTK.Collections import Collection
43 from MMTK.ParticleProperties import ParticleVector
44 from MMTK.Trajectory import SnapshotGenerator, Trajectory, TrajectoryOutput
45
46
47 from nMOLDYN.Analysis.Analysis import Analysis, setUniverseContents
48 from nMOLDYN.Core.Error import Error
49 from nMOLDYN.Core.IOFiles import convertNetCDFToASCII
50 from nMOLDYN.Core.Logger import LogMessage
51 from nMOLDYN.Core.Mathematics import correlation, differentiate, factorial, FFT, gaussianWindow, invFFT, preparePP
52
53
54
55
57 """Sets up a Mean Square Displacement analysis.
58
59 A Subclass of nMOLDYN.Analysis.Analysis.
60
61 Constructor: MeanSquareDisplacement(|parameters| = None)
62
63 Arguments:
64
65 - |parameters| -- a dictionnary of the input parameters, or 'None' to set up the analysis without parameters.
66 * trajectory -- a trajectory file name or an instance of MMTK.Trajectory.Trajectory class.
67 * timeinfo -- a string of the form 'first:last:step' where 'first' is an integer specifying the first frame
68 number to consider, 'last' is an integer specifying the last frame number to consider and
69 'step' is an integer specifying the step number between two frames.
70 * projection -- a string of the form 'vx,vy,vz' specifying the vector along which the analysis
71 will be computed. 'vx', 'vy', and 'vz' are floats specifying respectively the x, y and z value
72 of that vector.
73 * subset -- a selection string specifying the atoms to consider for the analysis.
74 * deuteration -- a selection string specifying the hydrogen atoms whose atomic parameters will be those of the deuterium.
75 * weights -- a string equal to 'equal', 'mass', 'coherent' , 'incoherent' or 'atomicNumber' that specifies the weighting
76 scheme to use.
77 * msd -- the output NetCDF file name. A CDL version of this file will also be generated with the '.cdl' extension
78 instead of the '.nc' extension.
79 * pyroserver -- a string specifying if Pyro will be used and how to run the analysis.
80
81 Running modes:
82
83 - To run the analysis do: a.runAnalysis() where a is the analysis object.
84 - To estimate the analysis do: a.estimateAnalysis() where a is the analysis object.
85 - To save the analysis to 'file' file name do: a.saveAnalysis(file) where a is the analysis object.
86
87 Comments:
88
89 - The algorithm is based on the Fast Correlation Algorithm (FCA) algorithm
90 """
91
92 inputParametersNames = ('trajectory',\
93 'timeinfo',\
94 'projection',\
95 'subset',\
96 'deuteration',\
97 'weights',\
98 'msd',\
99 'pyroserver',\
100 )
101
102 shortName = 'MSD'
103
104
105 canBeEstimated = True
106
108 """The constructor. Insures that the class can not be instanciated directly from here.
109 """
110 raise Error('This class can not be instanciated directly.')
111
113 """Initializes the analysis (e.g. parses and checks input parameters, set some variables ...).
114 """
115
116 self.parseInputParameters()
117
118 if self.pyroServer != 'monoprocessor':
119 if self.statusBar is not None:
120 LogMessage('warning', 'The job status bar will be desactivated for compatibility with Pyro module.')
121 self.statusBar = None
122
123 self.buildTimeInfo()
124
125 self.subset = self.subsetSelection(self.universe, self.subset)
126 self.nAtoms = self.subset.numberOfAtoms()
127
128 self.deuteration = self.deuterationSelection(self.universe, self.deuteration)
129 self.deuteration = Collection(list(set(self.subset) & set(self.deuteration)))
130
131 self.weights = self.weightingScheme(self.universe, self.subset, self.deuteration, self.weights)
132
133 self.preLoadTrajectory('atom')
134
135
136 self.MSD = N.zeros((self.nFrames), N.Float)
137
138 if self.pyroServer != 'monoprocessor':
139
140 delattr(self, 'trajectory')
141
142 - def calc(self, atom, trajname):
143 """Calculates the atomic term.
144
145 @param atom: the atom on which the atomic term has been calculated.
146 @type atom: an instance of MMTK.Atom class.
147
148 @param trajname: the name of the trajectory file name.
149 @type trajname: string
150 """
151
152 if self.preLoadedTraj is None:
153 if self.pyroServer == 'monoprocessor':
154 t = self.trajectory
155 else:
156
157 t = Trajectory(None, trajname, 'r')
158
159
160
161 series = t.readParticleTrajectory(atom, first = self.first, last = self.last, skip = self.skip).array
162
163 else:
164 series = self.preLoadedTraj[atom]
165
166 if self.projection is None:
167
168
169 dsq = N.add.reduce(series * series,1)
170 else:
171
172 series = N.dot(series,self.projection)
173 dsq = series * series
174
175
176 sum_dsq1 = N.add.accumulate(dsq)
177
178
179 sum_dsq2 = N.add.accumulate(dsq[::-1])
180
181
182 sumsq = 2.*sum_dsq1[-1]
183
184
185
186
187 Saabb = sumsq - N.concatenate(([0.], sum_dsq1[:-1])) - N.concatenate(([0.], sum_dsq2[:-1]))
188
189
190
191 Saabb = Saabb / (len(dsq) - N.arange(len(dsq)))
192 Sab = 2.*correlation(series)
193
194 atomicMSD = Saabb - Sab
195
196 if self.preLoadedTraj is not None:
197 del self.preLoadedTraj[atom]
198
199 return atomicMSD
200
202 """
203 """
204 N.add(self.MSD, self.weights[atom]*x, self.MSD)
205
207 """Finalizes the calculations (e.g. averaging the total term, output files creations ...).
208 """
209
210 print self.times
211 print self.MSD
212
213 outputFile = NetCDFFile(self.outputMSD, 'w')
214 outputFile.title = self.__class__.__name__
215 outputFile.jobinfo = self.information + '\nOutput file written on: %s\n\n' % asctime()
216
217
218 outputFile.createDimension('NFRAMES', self.nFrames)
219
220
221
222 TIMES = outputFile.createVariable('time', N.Float, ('NFRAMES',))
223 TIMES[:] = self.times[:]
224 TIMES.units = 'ps'
225
226
227 MSD = outputFile.createVariable('msd', N.Float, ('NFRAMES',))
228 MSD[:] = self.MSD[:]
229 MSD.units = 'nm^2'
230
231 outputFile.close()
232
233 self.toPlot = {'netcdf' : self.outputMSD, 'xVar' : 'time', 'yVar' : 'msd'}
234
235
236 convertNetCDFToASCII(inputFile = self.outputMSD,\
237 outputFile = os.path.splitext(self.outputMSD)[0] + '.cdl',\
238 variables = ['time', 'msd'])
239
241 """Returns the atomic Mean-Square-Displacement.
242
243 @param atom: the atom on which the atomic MSD has been calculated.
244 @type atom: an instance of MMTK.Atom class.
245
246 @param series: a array of dimension (self.nFrames,3) specifying the coordinates of atom |atom|
247 for the selected frames.
248 @type series: NumPy array
249
250 @return: the MSD computed for atom |atom| with trajectory |series|.
251 @rtype: Numpy array
252 """
253
254
255
256
258 """Sets up a Root Mean Square Deviation analysis.
259
260 A Subclass of nMOLDYN.Analysis.Analysis.
261
262 Constructor: RootMeanSquareDeviation(|parameters| = None)
263
264 Arguments:
265
266 - |parameters| -- a dictionnary of the input parameters, or 'None' to set up the analysis without parameters.
267 * trajectory -- a trajectory file name or an instance of MMTK.Trajectory.Trajectory class.
268 * timeinfo -- a string of the form 'first:last:step' where 'first' is an integer specifying the first frame
269 number to consider, 'last' is an integer specifying the last frame number to consider and
270 'step' is an integer specifying the step number between two frames.
271 * referenceframe -- an integer in [1,len(trajectory)] specifying which frame should be the reference.
272 * subset -- a selection string specifying the atoms to consider for the analysis.
273 * deuteration -- a selection string specifying the hydrogen atoms whose atomic parameters will be those of the deuterium.
274 * weights -- a string equal to 'equal', 'mass', 'coherent' , 'incoherent' or 'atomicNumber' that specifies the weighting
275 scheme to use.
276 * rmsd -- the output NetCDF file name. A CDL version of this file will also be generated with the '.cdl' extension
277 instead of the '.nc' extension.
278 * pyroserver -- a string specifying if Pyro will be used and how to run the analysis.
279
280 Running modes:
281
282 - To run the analysis do: a.runAnalysis() where a is the analysis object.
283 - To estimate the analysis do: a.estimateAnalysis() where a is the analysis object.
284 - To save the analysis to 'file' file name do: a.saveAnalysis(file) where a is the analysis object.
285 """
286
287
288 inputParametersNames = ('trajectory',\
289 'timeinfo',\
290 'referenceframe',\
291 'subset',\
292 'deuteration',\
293 'weights',\
294 'rmsd',\
295 'pyroserver',\
296 )
297
298 shortName = 'RMSD'
299
300
301 canBeEstimated = True
302
304 """The constructor. Insures that the class can not be instanciated directly from here.
305 """
306 raise Error('This class can not be instanciated directly.')
307
309 """Initializes the analysis (e.g. parses and checks input parameters, set some variables ...).
310 """
311
312
313 self.parseInputParameters()
314
315 if self.pyroServer != 'monoprocessor':
316 if self.statusBar is not None:
317 LogMessage('warning', 'The job status bar will be desactivated for compatibility with Pyro module.')
318 self.statusBar = None
319
320 self.buildTimeInfo()
321
322 if (self.referenceFrame < 0) or (self.referenceFrame > len(self.trajectory)):
323 self.referenceFrame = 0
324 LogMessage('warning', 'The reference frame must be an integer in [1,%s].\n\
325 It will be set to 1 for the running analysis.' % len(self.trajectory), ['console'])
326
327 self.subset = self.subsetSelection(self.universe, self.subset)
328 self.nAtoms = self.subset.numberOfAtoms()
329
330 self.deuteration = self.deuterationSelection(self.universe, self.deuteration)
331 self.deuteration = Collection(list(set(self.subset) & set(self.deuteration)))
332
333 self.weights = self.weightingScheme(self.universe, self.subset, self.deuteration, self.weights)
334
335 self.preLoadTrajectory('atom')
336
337
338 self.RMSD = N.zeros((self.nFrames), N.Float)
339
340 if self.pyroServer != 'monoprocessor':
341
342 delattr(self, 'trajectory')
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367 - def calc(self, atom, trajname):
368 """Calculates the atomic term.
369
370 @param atom: the atom on which the atomic term has been calculated.
371 @type atom: an instance of MMTK.Atom class.
372
373 @param trajname: the name of the trajectory file name.
374 @type trajname: string
375
376 @note: an atom-by-atom implementation was prefered than a frame-by-frame
377 implementation of the type:
378 msd = t.configuration[frame] - t.configuration[self.referenceFrame]
379 msd = self.weights * msd * msd
380 self.RMSD[frameIndex] = N.sqrt(N.add.reduce(msd))
381 """
382
383 if self.preLoadedTraj is None:
384 if self.pyroServer == 'monoprocessor':
385 t = self.trajectory
386
387 else:
388
389 t = Trajectory(None, trajname, 'r')
390
391
392
393 series = t.readParticleTrajectory(atom, first = self.first, last = self.last, skip = self.skip).array
394
395 else:
396 series = self.preLoadedTraj[atom]
397
398 N.add(series, -series[self.referenceFrame,:], series)
399
400 atomicRMSD = N.add.reduce(series**2,1)
401
402 if self.preLoadedTraj is not None:
403 del self.preLoadedTraj[atom]
404
405 return atomicRMSD
406
408 """
409 """
410
411 N.add(self.RMSD, self.weights[atom]*x, self.RMSD)
412
414 """Finalizes the calculations (e.g. averaging the total term, output files creations ...)
415 """
416
417 self.RMSD = N.sqrt(self.RMSD)
418
419
420 outputFile = NetCDFFile(self.outputRMSD, 'w')
421 outputFile.title = self.__class__.__name__
422 outputFile.jobinfo = self.information + '\nOutput file written on: %s\n\n' % asctime()
423
424
425 outputFile.createDimension('NFRAMES', self.nFrames)
426
427
428
429 TIMES = outputFile.createVariable('time', N.Float, ('NFRAMES',))
430 TIMES[:] = self.times[:]
431 TIMES.units = 'ps'
432
433
434 RMSD = outputFile.createVariable('rmsd', N.Float, ('NFRAMES',))
435 RMSD[:] = self.RMSD[:]
436 RMSD.units = 'nm'
437
438 outputFile.close()
439
440 self.toPlot = {'netcdf' : self.outputRMSD, 'xVar' : 'time', 'yVar' : 'rmsd'}
441
442
443 convertNetCDFToASCII(inputFile = self.outputRMSD,\
444 outputFile = os.path.splitext(self.outputRMSD)[0] + '.cdl',\
445 variables = ['time', 'rmsd'])
446
447
448
449
451 """Sets up a Cartesian Velocity AutoCorrelation analysis.
452
453 A Subclass of nMOLDYN.Analysis.Analysis.
454
455 Constructor: CartesianVelocityAutoCorrelationFunction(|parameters| = None)
456
457 Arguments:
458
459 - |parameters| -- a dictionnary of the input parameters, or 'None' to set up the analysis without parameters.
460 * trajectory -- a trajectory file name or an instance of MMTK.Trajectory.Trajectory class.
461 * timeinfo -- a string of the form 'first:last:step' where 'first' is an integer specifying the first frame
462 number to consider, 'last' is an integer specifying the last frame number to consider and
463 'step' is an integer specifying the step number between two frames.
464 * differentiation -- an integer in [0,5] specifying the order of the differentiation used to get the velocities
465 out of the coordinates. 0 means that the velocities are already present in the trajectory loaded
466 for analysis.
467 * projection -- a string of the form 'vx,vy,vz' specifying the vector along which the analysis
468 will be computed. 'vx', 'vy', and 'vz' are floats specifying respectively the x, y and z value
469 of that vector.
470 * normalize -- a string being one of 'Yes' or 'No' specifying whether the analysis should be normalized to 1
471 at t = 0 ('Yes') or not ('No').
472 * subset -- a selection string specifying the atoms to consider for the analysis.
473 * deuteration -- a selection string specifying the hydrogen atoms whose atomic parameters will be those of the deuterium.
474 * weights -- a string equal to 'equal', 'mass', 'coherent' , 'incoherent' or 'atomicNumber' that specifies the weighting
475 scheme to use.
476 * vacf -- the output NetCDF file name. A CDL version of this file will also be generated with the '.cdl' extension
477 instead of the '.nc' extension.
478 * pyroserver -- a string specifying if Pyro will be used and how to run the analysis.
479
480 Running modes:
481
482 - To run the analysis do: a.runAnalysis() where a is the analysis object.
483 - To estimate the analysis do: a.estimateAnalysis() where a is the analysis object.
484 - To save the analysis to 'file' file name do: a.saveAnalysis(file) where a is the analysis object.
485 """
486
487
488 inputParametersNames = ('trajectory',\
489 'timeinfo',\
490 'differentiation',\
491 'projection',\
492 'normalize',\
493 'subset',\
494 'deuteration',\
495 'weights',\
496 'vacf',\
497 'pyroserver',\
498 )
499
500 shortName = 'VACF'
501
502
503 canBeEstimated = True
504
506 """The constructor. Insures that the class can not be instanciated directly from here.
507 """
508 raise Error('This class can not be instanciated directly.')
509
511 """Initializes the analysis (e.g. parses and checks input parameters, set some variables ...).
512 """
513
514
515 self.parseInputParameters()
516
517 if self.pyroServer != 'monoprocessor':
518 if self.statusBar is not None:
519 LogMessage('warning', 'The job status bar will be desactivated for compatibility with Pyro module.')
520 self.statusBar = None
521
522 self.buildTimeInfo()
523
524 self.subset = self.subsetSelection(self.universe, self.subset)
525 self.nAtoms = self.subset.numberOfAtoms()
526
527 self.deuteration = self.deuterationSelection(self.universe, self.deuteration)
528 self.deuteration = Collection(list(set(self.subset) & set(self.deuteration)))
529
530 self.weights = self.weightingScheme(self.universe, self.subset, self.deuteration, self.weights)
531
532 self.preLoadTrajectory('atom', self.differentiation)
533
534 self.VACF = N.zeros((self.nFrames), N.Float)
535
536 if self.pyroServer != 'monoprocessor':
537
538 delattr(self, 'trajectory')
539
540 - def calc(self, atom, trajname):
541 """Calculates the atomic term.
542
543 @param atom: the atom on which the atomic term has been calculated.
544 @type atom: an instance of MMTK.Atom class.
545
546 @param trajname: the name of the trajectory file name.
547 @type trajname: string
548 """
549
550 if self.preLoadedTraj is None:
551 if self.pyroServer == 'monoprocessor':
552 t = self.trajectory
553
554 else:
555
556 t = Trajectory(None, trajname, 'r')
557
558 if self.differentiation != 0:
559
560 if self.preLoadedTraj is None:
561
562
563 series = t.readParticleTrajectory(atom, first = self.first, last = self.last, skip = self.skip).array
564 else:
565 series = self.preLoadedTraj[atom]
566
567
568 for axis in range(3):
569
570
571
572
573
574 series[:,axis] = differentiate(series[:,axis], self.differentiation, self.dt)
575 else:
576 if self.preLoadedTraj is None:
577
578
579 series = t.readParticleTrajectory(atom, first = self.first, last = self.last, skip = self.skip, variable = 'velocities').array
580 else:
581 series = self.preLoadedTraj[atom]
582
583 if self.projection is None:
584
585 atomicVACF = correlation(series)/3.
586
587 else:
588
589 projectedTrajectory = N.dot(series, self.projection)
590
591 atomicVACF = correlation(projectedTrajectory)
592
593 if self.preLoadedTraj is not None:
594 del self.preLoadedTraj[atom]
595
596 return atomicVACF
597
599 """
600 """
601
602 N.add(self.VACF, self.weights[atom]*x, self.VACF)
603
605 """Finalizes the calculations (e.g. averaging the total term, output files creations ...)
606 """
607
608
609 if self.normalize:
610 if self.VACF[0] == 0:
611 raise Error('The normalization factor is equal to zero !!!')
612 else:
613 self.VACF = self.VACF/self.VACF[0]
614
615
616 outputFile = NetCDFFile(self.outputVACF, 'w')
617 outputFile.title = self.__class__.__name__
618 outputFile.jobinfo = self.information + '\nOutput file written on: %s\n\n' % asctime()
619
620
621 outputFile.createDimension('NFRAMES', self.nFrames)
622
623
624
625 TIMES = outputFile.createVariable('time', N.Float, ('NFRAMES',))
626 TIMES[:] = self.times[:]
627 TIMES.units = 'ps'
628
629
630 VACF = outputFile.createVariable('vacf', N.Float, ('NFRAMES',))
631 VACF[:] = self.VACF[:]
632 if self.normalize:
633 VACF.units = 'Unitless'
634 else:
635 VACF.units = 'nm^2*ps^-2'
636
637 outputFile.close()
638
639 self.toPlot = {'netcdf' : self.outputVACF, 'xVar' : 'time', 'yVar' : 'vacf'}
640
641
642 convertNetCDFToASCII(inputFile = self.outputVACF,\
643 outputFile = os.path.splitext(self.outputVACF)[0] + '.cdl',\
644 variables = ['time', 'vacf'])
645
646
647
648
650 """Sets up a Cartesian Density Of States analysis.
651
652 A Subclass of nMOLDYN.Analysis.Analysis.
653
654 Constructor: CartesianDensityOfStates(|parameters| = None)
655
656 Arguments:
657
658 - |parameters| -- a dictionnary of the input parameters, or 'None' to set up the analysis without parameters.
659 * trajectory -- a trajectory file name or an instance of MMTK.Trajectory.Trajectory class.
660 * timeinfo -- a string of the form 'first:last:step' where 'first' is an integer specifying the first frame
661 number to consider, 'last' is an integer specifying the last frame number to consider and
662 'step' is an integer specifying the step number between two frames.
663 * differentiation -- an integer in [0,5] specifying the order of the differentiation used to get the velocities
664 out of the coordinates. 0 means that the velocities are already present in the trajectory loaded
665 for analysis.
666 * projection -- a string of the form 'vx,vy,vz' specifying the vector along which the analysis
667 will be computed. 'vx', 'vy', and 'vz' are floats specifying respectively the x, y and z value
668 of that vector.
669 * fftwindow -- a float in ]0.0,100.0[ specifying the width of the gaussian, in percentage of the trajectory length
670 that will be used in the smoothing procedure.
671 * subset -- a selection string specifying the atoms to consider for the analysis.
672 * deuteration -- a selection string specifying the hydrogen atoms whose atomic parameters will be those of the deuterium.
673 * weights -- a string equal to 'equal', 'mass', 'coherent' , 'incoherent' or 'atomicNumber' that specifies the weighting
674 scheme to use.
675 * dos -- the output NetCDF file name. A CDL version of this file will also be generated with the '.cdl' extension
676 instead of the '.nc' extension.
677 * pyroserver -- a string specifying if Pyro will be used and how to run the analysis.
678
679 Running modes:
680
681 - To run the analysis do: a.runAnalysis() where a is the analysis object.
682 - To estimate the analysis do: a.estimateAnalysis() where a is the analysis object.
683 - To save the analysis to 'file' file name do: a.saveAnalysis(file) where a is the analysis object.
684 """
685
686
687 inputParametersNames = ('trajectory',\
688 'timeinfo',\
689 'differentiation',\
690 'projection',\
691 'fftwindow',\
692 'subset',\
693 'deuteration',\
694 'weights',\
695 'dos',\
696 'pyroserver',\
697 )
698
699 shortName = 'DOS'
700
701
702 canBeEstimated = True
703
705 """The constructor. Insures that the class can not be instanciated directly from here.
706 """
707 raise Error('This class can not be instanciated directly.')
708
710 """Initializes the analysis (e.g. parses and checks input parameters, set some variables ...).
711 """
712
713
714 self.parseInputParameters()
715
716 if self.pyroServer != 'monoprocessor':
717 if self.statusBar is not None:
718 LogMessage('warning', 'The job status bar will be desactivated for compatibility with Pyro module.')
719 self.statusBar = None
720
721 self.buildTimeInfo()
722
723 self.subset = self.subsetSelection(self.universe, self.subset)
724 self.nAtoms = self.subset.numberOfAtoms()
725
726 self.deuteration = self.deuterationSelection(self.universe, self.deuteration)
727 self.deuteration = Collection(list(set(self.subset) & set(self.deuteration)))
728
729 self.weights = self.weightingScheme(self.universe, self.subset, self.deuteration, self.weights)
730
731 self.preLoadTrajectory('atom', self.differentiation)
732
733 self.DOS = N.zeros((self.nFrames), N.Float)
734
735 if self.pyroServer != 'monoprocessor':
736
737 delattr(self, 'trajectory')
738
739 - def calc(self, atom, trajname):
740 """Calculates the atomic term.
741
742 @param atom: the atom on which the atomic term has been calculated.
743 @type atom: an instance of MMTK.Atom class.
744
745 @param trajname: the name of the trajectory file name.
746 @type trajname: string
747 """
748
749 if self.preLoadedTraj is None:
750 if self.pyroServer == 'monoprocessor':
751 t = self.trajectory
752
753 else:
754
755 t = Trajectory(None, trajname, 'r')
756
757 if self.differentiation != 0:
758
759 if self.preLoadedTraj is None:
760
761
762 series = t.readParticleTrajectory(atom, first = self.first, last = self.last, skip = self.skip).array
763 else:
764 series = self.preLoadedTraj[atom]
765
766
767 for axis in range(3):
768
769
770
771
772
773 series[:,axis] = differentiate(series[:,axis], self.differentiation, self.dt)
774 else:
775 if self.preLoadedTraj is None:
776
777
778 series = t.readParticleTrajectory(atom, first = self.first, last = self.last, skip = self.skip, variable = 'velocities').array
779 else:
780 series = self.preLoadedTraj[atom]
781
782 if self.projection is None:
783
784
785 aVACF = correlation(series)/3.
786 else:
787
788 projectedTrajectory = N.dot(series, self.projection)
789
790 aVACF = correlation(projectedTrajectory)
791
792
793
794 atomicDOS = FFT(gaussianWindow(aVACF, self.fftWindow)).real[:len(aVACF)]
795
796 if self.preLoadedTraj is not None:
797 del self.preLoadedTraj[atom]
798
799 return atomicDOS
800
802 """
803 """
804
805 N.add(self.DOS, self.weights[atom]*x, self.DOS)
806
808 """Finalizes the calculations (e.g. averaging the total term, output files creations ...)
809 """
810
811
812 frequencies = N.arange(self.nFrames)/(2.0*self.nFrames*self.dt)
813
814
815 self.DOS = 0.5*self.dt*self.DOS
816
817
818 outputFile = NetCDFFile(self.outputDOS, 'w')
819 outputFile.title = self.__class__.__name__
820 outputFile.jobinfo = self.information + '\nOutput file written on: %s\n\n' % asctime()
821
822
823 outputFile.createDimension('FREQSTEPS', self.nFrames)
824
825
826
827 FREQUENCIES = outputFile.createVariable('frequency', N.Float, ('FREQSTEPS',))
828 FREQUENCIES[:] = frequencies[:]
829 FREQUENCIES.units = 'THz'
830
831
832 DOS = outputFile.createVariable('dos', N.Float, ('FREQSTEPS',))
833 DOS[:] = self.DOS[:]
834 DOS.units = 'nm^2*ps^-1'
835
836 outputFile.close()
837
838 self.toPlot = {'netcdf' : self.outputDOS, 'xVar' : 'frequency', 'yVar' : 'dos'}
839
840
841 convertNetCDFToASCII(inputFile = self.outputDOS,\
842 outputFile = os.path.splitext(self.outputDOS)[0] + '.cdl',\
843 variables = ['frequency', 'dos'])
844
845
846
847
849 """Sets up an AutoRegressive Analysis analysis.
850
851 A Subclass of nMOLDYN.Analysis.Analysis.
852
853 Constructor: AutoRegressiveAnalysis(|parameters| = None)
854
855 Arguments:
856
857 - |parameters| -- a dictionnary of the input parameters, or 'None' to set up the analysis without parameters.
858 * trajectory -- a trajectory file name or an instance of MMTK.Trajectory.Trajectory class.
859 * timeinfo -- a string of the form 'first:last:step' where 'first' is an integer specifying the first frame
860 number to consider, 'last' is an integer specifying the last frame number to consider and
861 'step' is an integer specifying the step number between two frames.
862 * differentiation -- an integer in [0,5] specifying the order of the differentiation used to get the velocities
863 out of the coordinates. 0 means that the velocities are already present in the trajectory loaded
864 for analysis.
865 * projection -- a string of the form 'vx,vy,vz' specifying the vector along which the analysis
866 will be computed. 'vx', 'vy', and 'vz' are floats specifying respectively the x, y and z value
867 of that vector.
868 * armodelorder -- an integer in [1, len(trajectory)[ specifying the order of the model
869 * subset -- a selection string specifying the atoms to consider for the analysis.
870 * deuteration -- a selection string specifying the hydrogen atoms whose atomic parameters will be those of the deuterium.
871 * weights -- a string equal to 'equal', 'mass', 'coherent' , 'incoherent' or 'atomicNumber' that specifies the weighting
872 scheme to use.
873 * ara -- the output NetCDF file name. A CDL version of this file will also be generated with the '.cdl' extension
874 instead of the '.nc' extension.
875 * pyroserver -- a string specifying if Pyro will be used and how to run the analysis.
876
877 Running modes:
878
879 - To run the analysis do: a.runAnalysis() where a is the analysis object.
880 - To estimate the analysis do: a.estimateAnalysis() where a is the analysis object.
881 - To save the analysis to 'file' file name do: a.saveAnalysis(file) where a is the analysis object.
882 """
883
884
885 inputParametersNames = ('trajectory',\
886 'timeinfo',\
887 'differentiation',\
888 'projection',\
889 'armodelorder',\
890 'subset',\
891 'deuteration',\
892 'weights',\
893 'ara',\
894 'pyroserver',\
895 )
896
897 shortName = 'ARA'
898
899
900 canBeEstimated = True
901
903 """The constructor. Insures that the class can not be instanciated directly from here.
904 """
905 raise Error('This class can not be instanciated directly.')
906
908 """Initializes the analysis (e.g. parses and checks input parameters, set some variables ...).
909 """
910
911
912 self.parseInputParameters()
913
914 if self.pyroServer != 'monoprocessor':
915 if self.statusBar is not None:
916 LogMessage('warning', 'The job status bar will be desactivated for compatibility with Pyro module.')
917 self.statusBar = None
918
919 self.buildTimeInfo()
920
921 if (self.arModelOrder <= 0) or (self.arModelOrder >= self.nFrames):
922 raise Error('The AR order must be an integer in [1,%d[.' % self.nFrames)
923
924 self.subset = self.subsetSelection(self.universe, self.subset)
925 self.nAtoms = self.subset.numberOfAtoms()
926 self.deuteration = self.deuterationSelection(self.universe, self.deuteration)
927 self.deuteration = Collection(list(set(self.subset) & set(self.deuteration)))
928 self.weights = self.weightingScheme(self.universe, self.subset, self.deuteration, self.weights)
929
930 self.preLoadTrajectory('atom', self.differentiation)
931
932 self.ARModel = AveragedAutoRegressiveModel(self.arModelOrder, self.dt)
933
934 if self.pyroServer != 'monoprocessor':
935
936 delattr(self, 'trajectory')
937
938 - def calc(self, atom, trajname):
939 """Calculates the atomic term.
940
941 @param atom: the atom on which the atomic term has been calculated.
942 @type atom: an instance of MMTK.Atom class.
943
944 @param trajname: the name of the trajectory file name.
945 @type trajname: string
946 """
947
948 if self.preLoadedTraj is None:
949 if self.pyroServer == 'monoprocessor':
950 t = self.trajectory
951
952 else:
953
954 t = Trajectory(None, trajname, 'r')
955
956 if self.differentiation != 0:
957
958 if self.preLoadedTraj is None:
959
960
961 series = t.readParticleTrajectory(atom, first = self.first, last = self.last, skip = self.skip).array
962 else:
963 series = self.preLoadedTraj[atom]
964
965
966 for axis in range(3):
967
968
969
970
971
972 series[:,axis] = differentiate(series[:,axis], self.differentiation, self.dt)
973 else:
974 if self.preLoadedTraj is None:
975
976
977 series = t.readParticleTrajectory(atom, first = self.first, last = self.last, skip = self.skip, variable = 'velocities').array
978 else:
979 series = self.preLoadedTraj[atom]
980
981 if self.preLoadedTraj is not None:
982 del self.preLoadedTraj[atom]
983
984 return series
985
987 """
988 """
989
990
991 for axis in range(3):
992 self.ARModel.add(AutoRegressiveModel(self.arModelOrder, x[:, axis], self.dt), self.weights[atom])
993
995 """Finalizes the calculations (e.g. averaging the total term, output files creations ...)
996 """
997
998
999
1000
1001
1002
1003 vacf = self.ARModel.correlation(self.nFrames)
1004 c = vacf.values
1005 c_init = N.fabs(c[0].real)
1006 for i in range(len(c)):
1007 if N.fabs(c[i].real)/c_init < 1.e-10:
1008 break
1009 if i < len(c):
1010 vacf = InterpolatingFunction((vacf.axes[0][:i],), c[:i])
1011
1012
1013 omega_max = 1.1*N.pi/self.ARModel.delta_t
1014 omega = omega_max*N.arange(self.nFrames) / float(self.nFrames)
1015 dos = self.ARModel.spectrum(omega)
1016 dos = InterpolatingFunction((omega/(2.0*N.pi),), dos.values)
1017 dos.axes = (dos.axes[0],)
1018
1019
1020 poles = self.ARModel.poles()
1021 cpoles = N.conjugate(poles)
1022 coeff0 = N.conjugate(self.ARModel.coeff[0])
1023 beta = N.zeros((self.ARModel.order), N.Complex)
1024
1025 for i in range(self.ARModel.order):
1026 pole = poles[i]
1027 beta[i] = -(self.ARModel.sigsq*pole**(self.ARModel.order-1)/coeff0)/\
1028 (N.multiply.reduce((pole-poles)[:i])*\
1029 N.multiply.reduce((pole-poles)[i+1:])*\
1030 N.multiply.reduce(pole-1./cpoles)*\
1031 self.ARModel.variance)
1032
1033 beta = beta/sum(beta)
1034 msd = N.zeros((self.nFrames), N.Float)
1035 n = N.arange(self.nFrames)
1036
1037 for i in range(self.ARModel.order):
1038 pole = poles[i]
1039 msd = msd + (beta[i]*((pole**n-1.)*pole/(1-pole)**2 + n/(1-pole))).real
1040 msd = 6.0 * self.ARModel.delta_t**2 * self.ARModel.variance * msd
1041
1042 msd = InterpolatingFunction((self.ARModel.delta_t*n,), msd)
1043
1044
1045 friction = self.ARModel.frictionConstant()
1046 try:
1047 memory = self.ARModel.memoryFunction(self.ARModel.order+self.ARModel.order/2)
1048 except OverflowError:
1049 null = N.zeros((0), N.Float)
1050 memory = InterpolatingFunction((null,), null)
1051
1052 memory.friction = friction
1053
1054
1055 outputFile = NetCDFFile(self.outputARA, 'w')
1056 outputFile.title = self.__class__.__name__
1057 outputFile.comment = 'Friction constant: %20.15f' % memory.friction
1058 outputFile.comment += 'Variance: %.15f\nSigma: %.15f' % (self.ARModel.variance, self.ARModel.sigma)
1059 outputFile.jobinfo = self.information + '\nOutput file written on: %s\n\n' % asctime()
1060
1061
1062 outputFile.createDimension('NFRAMES_VACF', len(vacf.axes[0]))
1063 outputFile.createDimension('NFRAMES_MSD', len(msd.axes[0]))
1064 outputFile.createDimension('NFRAMES_MEMORY', len(memory.axes[0]))
1065 outputFile.createDimension('FREQSTEPS', len(dos.axes[0]))
1066 outputFile.createDimension('ARMODELORDER', self.arModelOrder)
1067
1068
1069
1070
1071 TIMES = outputFile.createVariable('time_vacf', N.Float, ('NFRAMES_VACF',))
1072 TIMES[:] = vacf.axes[0]
1073 TIMES.units = 'ps'
1074
1075 VACF = outputFile.createVariable('vacf', N.Float, ('NFRAMES_VACF',))
1076 VACF[:] = vacf.values.real
1077 VACF.units = 'nm^2*ps^-2'
1078
1079
1080 FREQUENCIES = outputFile.createVariable('frequency', N.Float, ('FREQSTEPS',))
1081 FREQUENCIES[:] = dos.axes[0]
1082 FREQUENCIES.units = 'THz'
1083
1084 DOS = outputFile.createVariable('dos', N.Float, ('FREQSTEPS',))
1085 DOS[:] = dos.values.real
1086 DOS.units = 'nm2*ps-1'
1087
1088
1089 TIMES = outputFile.createVariable('time_msd', N.Float, ('NFRAMES_MSD',))
1090 TIMES[:] = msd.axes[0]
1091 TIMES.units = 'ps'
1092
1093 MSD = outputFile.createVariable('msd', N.Float, ('NFRAMES_MSD',))
1094 MSD[:] = msd.values.real
1095 MSD.units = 'nm^2'
1096
1097
1098 TIMES = outputFile.createVariable('time_memory', N.Float, ('NFRAMES_MEMORY',))
1099 TIMES[:] = memory.axes[0]
1100 TIMES.units = 'ps'
1101
1102 MEMORY = outputFile.createVariable('memory_function', N.Float, ('NFRAMES_MEMORY',))
1103 MEMORY[:] = memory.values.real
1104 MEMORY.units = 'unitless'
1105
1106
1107 COEFFINDEXES = outputFile.createVariable('n', N.Int, ('ARMODELORDER',))
1108 COEFFINDEXES[:] = 1 + N.arange(0, self.arModelOrder)
1109 COEFFINDEXES.units = 'unitless'
1110
1111 COEFFS = outputFile.createVariable('ar_coefficients', N.Float, ('ARMODELORDER',))
1112 COEFFS[:] = self.ARModel.coeff[::-1]
1113 COEFFS.units = 'unitless'
1114
1115 outputFile.close()
1116
1117 self.toPlot = None
1118
1119
1120 convertNetCDFToASCII(inputFile = self.outputARA,\
1121 outputFile = os.path.splitext(self.outputARA)[0] + '.cdl',\
1122 variables = ['time_vacf',\
1123 'vacf',\
1124 'frequency',\
1125 'dos',\
1126 'time_msd',\
1127 'msd',\
1128 'time_memory',\
1129 'memory_function',\
1130 'n',\
1131 'ar_coefficients'])
1132
1133
1134
1135
1137 """Sets up a Pass-Band Trajectory Filter analysis.
1138
1139 A Subclass of nMOLDYN.Analysis.Analysis.
1140
1141 Constructor: PassBandFilteredTrajectory(|parameters| = None)
1142
1143 Arguments:
1144
1145 - |parameters| -- a dictionnary of the input parameters, or 'None' to set up the analysis without parameters.
1146 * trajectory -- a trajectory file name or an instance of MMTK.Trajectory.Trajectory class.
1147 * timeinfo -- a string of the form 'first:last:step' where 'first' is an integer specifying the first frame
1148 number to consider, 'last' is an integer specifying the last frame number to consider and
1149 'step' is an integer specifying the step number between two frames.
1150 * filter -- a string of the form 'low:high' where 'low' and 'high' are floats specifying respectively
1151 the lower and the upper bounds of the pass-band filter.
1152 * subset -- a selection string specifying the atoms to consider for the analysis.
1153 * pbft -- the output NetCDF file name.
1154 * pyroserver -- a string specifying if Pyro will be used and how to run the analysis.
1155
1156 Running modes:
1157
1158 - To run the analysis do: a.runAnalysis() where a is the analysis object.
1159 - To estimate the analysis do: a.estimateAnalysis() where a is the analysis object.
1160 - To save the analysis to 'file' file name do: a.saveAnalysis(file) where a is the analysis object.
1161 """
1162
1163
1164 inputParametersNames = ('trajectory',\
1165 'timeinfo',\
1166 'filter',\
1167 'subset',\
1168 'pbft',\
1169 'pyroserver',\
1170 )
1171
1172 shortName = 'PBFT'
1173
1174
1175 canBeEstimated = True
1176
1178 """The constructor. Insures that the class can not be instanciated directly from here.
1179 """
1180 raise Error('This class can not be instanciated directly.')
1181
1183 """Initializes the analysis (e.g. parses and checks input parameters, set some variables ...).
1184 """
1185
1186
1187 self.parseInputParameters()
1188
1189 if self.pyroServer != 'monoprocessor':
1190 if self.statusBar is not None:
1191 LogMessage('warning', 'The job status bar will be desactivated for compatibility with Pyro module.')
1192 self.statusBar = None
1193
1194 self.buildTimeInfo()
1195
1196 self.subset = self.subsetSelection(self.universe, self.subset)
1197 self.nAtoms = self.subset.numberOfAtoms()
1198
1199 self.preLoadTrajectory('atom')
1200
1201
1202 frequencies = N.arange(self.nFrames)/(2.0*self.nFrames*self.dt)
1203
1204 indexMin = N.searchsorted(frequencies, self.freqMin)
1205
1206 indexMax = N.searchsorted(frequencies, self.freqMax)
1207
1208 halfFilter = N.zeros(self.nFrames, N.Float)
1209 halfFilter[indexMin:indexMax] = 1
1210
1211 self.pbFilter = N.zeros(2*self.nFrames - 2, N.Float)
1212 self.pbFilter[:self.nFrames] = halfFilter
1213 self.pbFilter[self.nFrames:] = halfFilter[-2:0:-1]
1214
1215
1216
1217 self.filteredTrajectory = {}
1218
1219 if self.pyroServer != 'monoprocessor':
1220
1221 delattr(self, 'trajectory')
1222
1223 - def calc(self, atom, trajname):
1224 """Calculates the atomic term.
1225
1226 @param atom: the atom on which the atomic term has been calculated.
1227 @type atom: an instance of MMTK.Atom class.
1228
1229 @param trajname: the name of the trajectory file name.
1230 @type trajname: string
1231 """
1232
1233 if self.preLoadedTraj is None:
1234 if self.pyroServer == 'monoprocessor':
1235 t = self.trajectory
1236
1237 else:
1238
1239 t = Trajectory(None, trajname, 'r')
1240
1241
1242
1243 series = t.readParticleTrajectory(atom, first = self.first, last = self.last, skip = self.skip).array
1244
1245 else:
1246 series = self.preLoadedTraj[atom]
1247
1248 if self.preLoadedTraj is not None:
1249 del self.preLoadedTraj[atom]
1250
1251 return series
1252
1254 """
1255 """
1256
1257 unfiltered = N.zeros(2*self.nFrames - 2, 'd')
1258 atomicPBFT = N.zeros((self.nFrames,3), 'd')
1259
1260
1261 for coord in range(3):
1262 unfiltered[:self.nFrames] = x[:,coord]
1263 unfiltered[self.nFrames:] = x[:,coord][-2:0:-1]
1264 temp = FFT(unfiltered)
1265 atomicPBFT[:,coord] = invFFT(self.pbFilter*temp)[:self.nFrames].real
1266
1267 self.filteredTrajectory[atom.index] = atomicPBFT
1268
1270 """Finalizes the calculations (e.g. averaging the total term, output files creations ...)
1271 """
1272
1273
1274 outputFile = Trajectory(self.subset, self.outputPBFT, "w")
1275 outputFile.jobinfo = self.information + '\nOutput file written on: %s\n\n' % asctime()
1276
1277 if self.pyroServer == 'monoprocessor':
1278 t = self.trajectory
1279
1280 else:
1281
1282 t = Trajectory(None, self.trajectoryFilename, 'r')
1283
1284
1285 snapshot = SnapshotGenerator(self.universe, actions = [TrajectoryOutput(outputFile, ["all"], 0, None, 1)])
1286
1287
1288 for frameIndex in range(self.nFrames):
1289 frame = self.frameIndexes[frameIndex]
1290 self.universe.setConfiguration(t.configuration[frame])
1291 for at in self.subset.atomList():
1292 at.setPosition(self.filteredTrajectory[at.index][frameIndex,:])
1293
1294 snapshot(data = {'time': self.times[frameIndex]})
1295
1296 outputFile.close()
1297
1298 t.close()
1299
1300 self.toPlot = None
1301
1302
1303
1304
1306 """Sets up a Radius Of Gyration analysis.
1307
1308 A Subclass of nMOLDYN.Analysis.Analysis.
1309
1310 Constructor: RadiusOfGyration(|parameters| = None)
1311
1312 Arguments:
1313
1314 - |parameters| -- a dictionnary of the input parameters, or 'None' to set up the analysis without parameters.
1315 * trajectory -- a trajectory file name or an instance of MMTK.Trajectory.Trajectory class.
1316 * timeinfo -- a string of the form 'first:last:step' where 'first' is an integer specifying the first frame
1317 number to consider, 'last' is an integer specifying the last frame number to consider and
1318 'step' is an integer specifying the step number between two frames.
1319 * subset -- a selection string specifying the atoms to consider for the analysis.
1320 * rog -- the output NetCDF file name. A CDL version of this file will also be generated with the '.cdl' extension
1321 instead of the '.nc' extension.
1322 * pyroserver -- a string specifying if Pyro will be used and how to run the analysis.
1323
1324 Running modes:
1325
1326 - To run the analysis do: a.runAnalysis() where a is the analysis object.
1327 - To estimate the analysis do: a.estimateAnalysis() where a is the analysis object.
1328 - To save the analysis to 'file' file name do: a.saveAnalysis(file) where a is the analysis object.
1329 """
1330
1331
1332 inputParametersNames = ('trajectory',\
1333 'timeinfo',\
1334 'subset',\
1335 'rog',\
1336 'pyroserver',\
1337 )
1338
1339 shortName = 'ROG'
1340
1341
1342 canBeEstimated = True
1343
1345 """The constructor. Insures that the class can not be instanciated directly from here.
1346 """
1347 raise Error('This class can not be instanciated directly.')
1348
1350 """Initializes the analysis (e.g. parses and checks input parameters, set some variables ...).
1351 """
1352
1353
1354 self.parseInputParameters()
1355
1356 if self.pyroServer != 'monoprocessor':
1357 if self.statusBar is not None:
1358 LogMessage('warning', 'The job status bar will be desactivated for compatibility with Pyro module.')
1359 self.statusBar = None
1360
1361 self.buildTimeInfo()
1362
1363 self.subset = self.subsetSelection(self.universe, self.subset)
1364 self.nAtoms = self.subset.numberOfAtoms()
1365
1366 self.preLoadTrajectory('atom')
1367
1368
1369 self.cog = N.zeros((self.nFrames,3), N.Float)
1370
1371
1372 self.sumRi = N.zeros((self.nFrames,3), N.Float)
1373
1374
1375 self.ROG = N.zeros((self.nFrames), N.Float)
1376
1377 if self.pyroServer != 'monoprocessor':
1378
1379 delattr(self, 'trajectory')
1380
1381 - def calc(self, atom, trajname):
1382 """Calculates the atomic term.
1383
1384 @param atom: the atom on which the atomic term has been calculated.
1385 @type atom: an instance of MMTK.Atom class.
1386
1387 @param trajname: the name of the trajectory file name.
1388 @type trajname: string
1389 """
1390
1391 if self.preLoadedTraj is None:
1392 if self.pyroServer == 'monoprocessor':
1393 t = self.trajectory
1394
1395 else:
1396
1397 t = Trajectory(None, trajname, 'r')
1398
1399
1400
1401 series = t.readParticleTrajectory(atom, first = self.first, last = self.last, skip = self.skip).array
1402
1403 else:
1404 series = self.preLoadedTraj[atom]
1405
1406 if self.preLoadedTraj is not None:
1407 del self.preLoadedTraj[atom]
1408
1409 return series
1410
1412 """
1413 """
1414
1415 N.add(self.cog, atom.mass()*x, self.cog)
1416 N.add(self.sumRi, x, self.sumRi)
1417 N.add(self.ROG, N.add.reduce(x**2,1), self.ROG)
1418
1420 """Finalizes the calculations (e.g. averaging the total term, output files creations ...).
1421 """
1422
1423 self.cog = self.cog/self.subset.mass()
1424
1425 N.add(self.ROG, -2.0 * N.add.reduce(self.sumRi*self.cog,1), self.ROG)
1426 N.add(self.ROG, float(self.nAtoms)*N.add.reduce(self.cog**2,1), self.ROG)
1427
1428 self.ROG = N.sqrt(self.ROG/float(self.nAtoms))
1429
1430
1431 outputFile = NetCDFFile(self.outputROG, 'w')
1432 outputFile.title = self.__class__.__name__
1433 outputFile.jobinfo = self.information + '\nOutput file written on: %s\n\n' % asctime()
1434
1435
1436 outputFile.createDimension('NFRAMES', self.nFrames)
1437
1438
1439
1440 TIMES = outputFile.createVariable('time', N.Float, ('NFRAMES',))
1441 TIMES[:] = self.times[:]
1442 TIMES.units = 'ps'
1443
1444
1445 ROG = outputFile.createVariable('rog', N.Float, ('NFRAMES',))
1446 ROG[:] = self.ROG[:]
1447 ROG.units = 'nm'
1448
1449 outputFile.close()
1450
1451 self.toPlot = {'netcdf' : self.outputROG, 'xVar' : 'time', 'yVar' : 'rog'}
1452
1453
1454 convertNetCDFToASCII(inputFile = self.outputROG,\
1455 outputFile = os.path.splitext(self.outputROG)[0] + '.cdl',\
1456 variables = ['time', 'rog'])
1457
1458
1459
1460
1462 """Sets up a Global Motion Trajectory Filter analysis.
1463
1464 A Subclass of nMOLDYN.Analysis.Analysis.
1465
1466 Constructor: GlobalMotionFilteredTrajectory(|parameters| = None)
1467
1468 Arguments:
1469
1470 - |parameters| -- a dictionnary of the input parameters, or 'None' to set up the analysis without parameters.
1471 * trajectory -- a trajectory file name or an instance of MMTK.Trajectory.Trajectory class.
1472 * timeinfo -- a string of the form 'first:last:step' where 'first' is an integer specifying the first frame
1473 number to consider, 'last' is an integer specifying the last frame number to consider and
1474 'step' is an integer specifying the step number between two frames.
1475 * subset -- a selection string specifying the atoms to consider for the analysis.
1476 * gmft -- the output NetCDF file name.
1477 * pyroserver -- a string specifying if Pyro will be used and how to run the analysis.
1478
1479 Running modes:
1480
1481 - To run the analysis do: a.runAnalysis() where a is the analysis object.
1482 - To estimate the analysis do: a.estimateAnalysis() where a is the analysis object.
1483 - To save the analysis to 'file' file name do: a.saveAnalysis(file) where a is the analysis object.
1484 """
1485
1486
1487 inputParametersNames = ('trajectory',\
1488 'timeinfo',\
1489 'subset',\
1490 'gmft',\
1491 'pyroserver',\
1492 )
1493
1494 shortName = 'GMFT'
1495
1496
1497 canBeEstimated = True
1498
1500 """The constructor. Insures that the class can not be instanciated directly from here.
1501 """
1502 raise Error('This class can not be instanciated directly.')
1503
1505 """Initializes the analysis (e.g. parses and checks input parameters, set some variables ...).
1506 """
1507
1508
1509 self.parseInputParameters()
1510
1511 if self.pyroServer != 'monoprocessor':
1512 if self.statusBar is not None:
1513 LogMessage('warning', 'The job status bar will be desactivated for compatibility with Pyro module.')
1514 self.statusBar = None
1515
1516 self.buildTimeInfo()
1517
1518 self.subset = self.subsetSelection(self.universe, self.subset)
1519 self.nAtoms = self.subset.numberOfAtoms()
1520
1521 self.preLoadTrajectory('frame')
1522
1523
1524
1525 self.filteredTrajectory = {}
1526
1527 for at in self.subset:
1528 self.filteredTrajectory[at.index] = N.zeros((self.nFrames, 3), N.Float)
1529
1530 self.initialConf = None
1531
1532 if self.pyroServer != 'monoprocessor':
1533
1534 delattr(self, 'trajectory')
1535
1536 - def calc(self, frameIndex, trajname):
1537 """Calculates the contribution for one frame.
1538
1539 @param frameIndex: the index of the frame in |self.frameIndexes| array.
1540 @type frameIndex: integer.
1541
1542 @param trajname: the name of the trajectory file name.
1543 @type trajname: string
1544 """
1545
1546 frame = self.frameIndexes[frameIndex]
1547
1548 if self.preLoadedTraj is None:
1549 if self.pyroServer == 'monoprocessor':
1550 t = self.trajectory
1551
1552 else:
1553
1554 t = Trajectory(None, trajname, 'r')
1555
1556 conf = t.configuration[frame]
1557 self.universe.setConfiguration(conf)
1558
1559 else:
1560
1561 self.universe.setConfiguration(self.preLoadedTraj[frameIndex])
1562
1563
1564 if frameIndex == 0:
1565 tr = self.subset.normalizingTransformation()
1566
1567
1568
1569 else:
1570 tr, rms = self.subset.findTransformation(self.initialConf)
1571
1572 if self.preLoadedTraj is not None:
1573 del self.preLoadedTraj[frameIndex]
1574
1575 return tr
1576
1577 - def combine(self, frameIndex, x):
1578 """
1579 """
1580
1581
1582 if frameIndex == 0:
1583 self.initialConf = self.universe.copyConfiguration()
1584
1585 self.subset.applyTransformation(x)
1586
1587 for atom in self.subset:
1588 self.filteredTrajectory[atom.index][frameIndex,:] = atom.position()
1589
1590
1592 """Finalizes the calculations (e.g. averaging the total term, output files creations ...).
1593 """
1594
1595
1596 outputFile = Trajectory(self.subset, self.outputGMFT, "w")
1597 outputFile.jobinfo = self.information + '\nOutput file written on: %s\n\n' % asctime()
1598
1599 if self.pyroServer != 'monoprocessor':
1600 t = self.trajectory
1601 else:
1602
1603 t = Trajectory(None, self.trajectoryFilename, 'r')
1604
1605
1606 snapshot = SnapshotGenerator(self.universe, actions = [TrajectoryOutput(outputFile, ["all"], 0, None, 1)])
1607
1608
1609 for frameIndex in range(self.nFrames):
1610 frame = self.frameIndexes[frameIndex]
1611 self.universe.setConfiguration(t.configuration[frame])
1612 for atom in self.subset.atomList():
1613 atom.setPosition(self.filteredTrajectory[atom.index][frameIndex,:])
1614
1615 snapshot(data = {'time': self.times[frameIndex]})
1616
1617 outputFile.close()
1618
1619 t.close()
1620
1621 self.toPlot = None
1622
1623
1624
1625
1627 """Sets up a Center Of Mass Trajectory analysis.
1628
1629 A Subclass of nMOLDYN.Analysis.Analysis.
1630
1631 Constructor: CenterOfMassTrajectory(|parameters| = None)
1632
1633 Arguments:
1634
1635 - |parameters| -- a dictionnary of the input parameters, or 'None' to set up the analysis without parameters.
1636 * trajectory -- a trajectory file name or an instance of MMTK.Trajectory.Trajectory class.
1637 * timeinfo -- a string of the form 'first:last:step' where 'first' is an integer specifying the first frame
1638 number to consider, 'last' is an integer specifying the last frame number to consider and
1639 'step' is an integer specifying the step number between two frames.
1640 * group -- a selection string specifying the groups of atoms on which the center of mass will be defined
1641 (one center of mass per group).
1642 * comt -- the output NetCDF file name.
1643 * pyroserver -- a string specifying if Pyro will be used and how to run the analysis.
1644
1645 Running modes:
1646
1647 - To run the analysis do: a.runAnalysis() where a is the analysis object.
1648 - To estimate the analysis do: a.estimateAnalysis() where a is the analysis object.
1649 - To save the analysis to 'file' file name do: a.saveAnalysis(file) where a is the analysis object.
1650 """
1651
1652
1653 inputParametersNames = ('trajectory',\
1654 'timeinfo',\
1655 'group',\
1656 'comt',\
1657 'pyroserver',\
1658 )
1659
1660 shortName = 'COMT'
1661
1662
1663 canBeEstimated = True
1664
1666 """The constructor. Insures that the class can not be instanciated directly from here.
1667 """
1668 raise Error('This class can not be instanciated directly.')
1669
1671 """Initializes the analysis (e.g. parses and checks input parameters, set some variables ...).
1672 """
1673
1674
1675 self.parseInputParameters()
1676
1677 if self.pyroServer != 'monoprocessor':
1678 if self.statusBar is not None:
1679 LogMessage('warning', 'The job status bar will be desactivated for compatibility with Pyro module.')
1680 self.statusBar = None
1681
1682 self.buildTimeInfo()
1683
1684 self.group = self.groupSelection(self.universe, self.group)
1685 self.nGroups = len(self.group)
1686
1687 self.preLoadTrajectory('frame')
1688
1689
1690 self.superAtomUniverse = self.universe.__copy__()
1691
1692 self.superAtomUniverse.removeObject(self.superAtomUniverse.objectList()[:])
1693
1694
1695
1696 self.comTrajectory = {}
1697
1698
1699 self.comToGroup = {}
1700
1701
1702 for g in range(self.nGroups):
1703 superAtom = Atom('H', name = 'SUPERATOM'+str(g+1))
1704 superAtom._mass = self.group[g].mass()
1705 self.superAtomUniverse.addObject(superAtom)
1706 self.comToGroup[superAtom] = self.group[g]
1707 self.comTrajectory[superAtom] = N.zeros((self.nFrames, 3), N.Float)
1708
1709 if self.pyroServer != 'monoprocessor':
1710
1711 delattr(self, 'trajectory')
1712
1713 - def calc(self, frameIndex, trajname):
1714 """Calculates the contribution for one frame.
1715
1716 @param frameIndex: the index of the frame in |self.frameIndexes| array.
1717 @type frameIndex: integer.
1718
1719 @param trajname: the name of the trajectory file name.
1720 @type trajname: string
1721 """
1722
1723 frame = self.frameIndexes[frameIndex]
1724
1725 if self.preLoadedTraj is None:
1726 if self.pyroServer == 'monoprocessor':
1727 t = self.trajectory
1728
1729 else:
1730
1731 t = Trajectory(None, trajname, 'r')
1732
1733 conf = t.configuration[frame]
1734 self.universe.setConfiguration(conf)
1735
1736 else:
1737
1738 self.universe.setConfiguration(self.preLoadedTraj[frameIndex])
1739
1740 for superAtom, group in self.comToGroup.items():
1741 self.comTrajectory[superAtom][frameIndex,:] = group.centerOfMass()
1742
1743 if self.preLoadedTraj is not None:
1744 del self.preLoadedTraj[frameIndex]
1745
1746 return None
1747
1749 """
1750 """
1751
1752 pass
1753
1755 """Finalizes the calculations (e.g. averaging the total term, output files creations ...).
1756 """
1757
1758
1759 outputFile = Trajectory(self.superAtomUniverse, self.outputCOMT, "w")
1760 outputFile.jobinfo = self.information + '\nOutput file written on: %s\n\n' % asctime()
1761
1762
1763 snapshot = SnapshotGenerator(self.superAtomUniverse,\
1764 actions = [TrajectoryOutput(outputFile, ["all"], 0, None, 1)])
1765
1766
1767 for frameIndex in range(self.nFrames):
1768 self.superAtomUniverse.configuration()
1769 for superAtom, superAtomTraj in self.comTrajectory.items():
1770 superAtom.setPosition(superAtomTraj[frameIndex,:])
1771 snapshot(data = {'time': self.times[frameIndex]})
1772
1773
1774
1775 outputFile.close()
1776
1777 self.toPlot = None
1778
1779
1780
1781
1783 """Sets up a Quasi Harmonic Analysis analysis.
1784
1785 A Subclass of nMOLDYN.Analysis.Analysis.
1786
1787 Constructor: QuasiHarmonicAnalysis(|parameters| = None)
1788
1789 Arguments:
1790
1791 - |parameters| -- a dictionnary of the input parameters, or 'None' to set up the analysis without parameters.
1792 * trajectory -- a trajectory file name or an instance of MMTK.Trajectory.Trajectory class.
1793 * timeinfo -- a string of the form 'first:last:step' where 'first' is an integer specifying the first frame
1794 number to consider, 'last' is an integer specifying the last frame number to consider and
1795 'step' is an integer specifying the step number between two frames.
1796 * temperature -- the temperature at which the MD was performed.
1797 * subset -- a selection string specifying the atoms to consider for the analysis.
1798 * qha -- the output NetCDF file name.
1799
1800 Running modes:
1801
1802 - To run the analysis do: a.runAnalysis() where a is the analysis object.
1803 - To estimate the analysis do: a.estimateAnalysis() where a is the analysis object.
1804 - To save the analysis to 'file' file name do: a.saveAnalysis(file) where a is the analysis object.
1805
1806 Comments:
1807
1808 - This analysis is used to get effective modes of vibration from fluctuations calculated by an MD simulation.
1809 The results of such an analysis can be seen by generating pseudo-trajectories reproducing the vibrations along
1810 a vibration mode.
1811 - For more details: Brooks et al., J. Comp. Chem. 1995, 16, 1522-1542.
1812 """
1813
1814
1815 inputParametersNames = ('trajectory',\
1816 'timeinfo',\
1817 'temperature',\
1818 'subset',\
1819 'qha')
1820
1821 shortName = 'QHA'
1822
1823
1824 canBeEstimated = False
1825
1827 """The constructor. Insures that the class can not be instanciated directly from here.
1828 """
1829 raise Error('This class can not be instanciated directly.')
1830
1832 """Initializes the analysis (e.g. parses and checks input parameters, set some variables ...).
1833 """
1834
1835 self.parseInputParameters()
1836
1837 self.buildTimeInfo()
1838
1839 self.subset = self.subsetSelection(self.universe, self.subset)
1840 self.nAtoms = self.subset.numberOfAtoms()
1841
1843 """Runs the analysis."""
1844
1845
1846 mask = self.subset.booleanMask()
1847
1848 M1_2 = N.zeros((3*self.nAtoms,), N.Float)
1849 weightList = [N.sqrt(el[1]) for el in sorted([(at.index, at.mass()) for at in self.subset])]
1850 for i in range(len(weightList)):
1851 M1_2[3*i:3*(i+1)] = weightList[i]
1852 invM1_2 = (1.0/M1_2)*N.identity(3*self.nAtoms, typecode = N.Float)
1853
1854
1855 initialStructure = self.trajectory.configuration[self.first]
1856
1857 averageStructure = ParticleVector(self.universe)
1858
1859 for conf in self.trajectory.configuration[self.first:self.last:self.skip]:
1860 averageStructure += self.universe.configurationDifference(initialStructure, conf)
1861
1862 averageStructure = averageStructure/self.nFrames
1863 averageStructure = initialStructure + averageStructure
1864
1865 mdr = N.zeros((self.nFrames, 3*self.nAtoms), N.Float)
1866
1867
1868 sigmaPrim = N.zeros((3*self.nAtoms, 3*self.nAtoms), N.Float)
1869 comp = 0
1870 for conf in self.trajectory.configuration[self.first:self.last:self.skip]:
1871 mdr[comp,:] = M1_2*N.ravel(N.compress(mask.array,(conf-averageStructure).array,0))
1872 sigmaPrim += mdr[comp,:, N.NewAxis] * mdr[N.NewAxis, comp, :]
1873 comp += 1
1874
1875 sigmaPrim = sigmaPrim/float(self.nFrames)
1876
1877 try:
1878
1879 omega, dx = Heigenvectors(sigmaPrim)
1880 except MemoryError:
1881 raise Error('Not enough memory to diagonalize the %sx%s fluctuation matrix.' % sigmaPrim.shape)
1882
1883
1884
1885
1886 omega = omega.real/(Units.kg*Units.m**2)
1887 omega = (Units.Hz/Units.invcm)*N.sqrt((Units.k_B*Units.K*self.temperature/Units.J)*(1.0/omega))
1888
1889 dx = dx.real
1890 dx = N.dot(invM1_2, dx)
1891
1892
1893 indices = N.argsort(omega)[::-1]
1894
1895 omega = N.take(omega, indices)
1896 dx = N.take(dx, indices)
1897
1898
1899 mdr = N.take(mdr, indices, axis = 1)
1900 at = N.dot(mdr,N.transpose(dx))
1901
1902
1903 outputFile = NetCDFFile(self.outputQHA, 'w')
1904 outputFile.title = self.__class__.__name__
1905 outputFile.jobinfo = self.information + '\nOutput file written on: %s\n\n' % asctime()
1906
1907
1908 self.universe.removeObject(self.universe.objectList()[:])
1909
1910 atoms = copy.deepcopy(self.subset.atomList())
1911
1912 for a in atoms:
1913 a.parent = None
1914 ac = AtomCluster(atoms,name='QHACluster')
1915 self.universe.addObject(ac)
1916
1917
1918
1919 outputFile.createDimension('NEIGENVALS', len(omega))
1920
1921
1922 outputFile.createDimension('UDESCR', len(self.universe.description()))
1923
1924
1925 outputFile.createDimension('NATOMS', self.nAtoms)
1926
1927
1928 outputFile.createDimension('NFRAMES', self.nFrames)
1929
1930
1931 outputFile.createDimension('NCOORDS', 3)
1932
1933
1934 outputFile.createDimension('3N', 3*self.nAtoms)
1935
1936 if self.universe.cellParameters() is not None:
1937 outputFile.createDimension('BOXDIM', len(self.universe.cellParameters()))
1938
1939
1940
1941 OMEGA = outputFile.createVariable('omega', N.Float, ('NEIGENVALS',))
1942 OMEGA[:] = omega
1943
1944
1945 DX = outputFile.createVariable('dx', N.Float, ('NEIGENVALS','NEIGENVALS'))
1946 DX[:,:] = dx
1947
1948
1949 TIMES = outputFile.createVariable('time', N.Float, ('NFRAMES',))
1950 TIMES[:] = self.times[:]
1951 TIMES.units = 'ps'
1952
1953
1954 MODE = outputFile.createVariable('mode', N.Float, ('3N',))
1955 MODE[:] = 1 + N.arange(3*self.nAtoms)
1956
1957
1958 LCI = outputFile.createVariable('local_character_indicator', N.Float, ('3N',))
1959 LCI[:] = (dx**4).sum(0)
1960
1961
1962 GCI = outputFile.createVariable('global_character_indicator', N.Float, ('3N',))
1963 GCI[:] = (N.sqrt(3.0*self.nAtoms)/(N.absolute(dx)).sum(0))**4
1964
1965
1966 AT = outputFile.createVariable('at', N.Float, ('NFRAMES','3N'))
1967 AT[:,:] = at[:,:]
1968
1969
1970 DESCRIPTION = outputFile.createVariable('description', N.Character, ('UDESCR',))
1971 DESCRIPTION[:] = self.universe.description()
1972
1973
1974 AVGSTRUCT = outputFile.createVariable('avgstruct', N.Float, ('NATOMS','NCOORDS'))
1975 AVGSTRUCT[:,:] = N.compress(mask.array,averageStructure.array,0)
1976
1977
1978 if self.universe.cellParameters() is not None:
1979 BOXSIZE = outputFile.createVariable('box_size', N.Float, ('BOXDIM',))
1980 BOXSIZE[:] = self.universe.cellParameters()
1981
1982 outputFile.close()
1983
1984 self.toPlot = None
1985
1986 return None
1987
1988
1989
1990
1992 """Sets up an Angular Correlation analysis.
1993
1994 A Subclass of nMOLDYN.Analysis.Analysis.
1995
1996 Constructor: AngularCorrelation(|parameters| = None)
1997
1998 Arguments:
1999
2000 - |parameters| -- a dictionnary of the input parameters, or 'None' to set up the analysis without parameters.
2001 * trajectory -- a trajectory file name or an instance of MMTK.Trajectory.Trajectory class.
2002 * timeinfo -- a string of the form 'first:last:step' where 'first' is an integer specifying the first frame
2003 number to consider, 'last' is an integer specifying the last frame number to consider and
2004 'step' is an integer specifying the step number between two frames.
2005 * triplet -- a selection string specifying the groups of three atoms that will define the plane on
2006 which the angular correlation will be computed.
2007 * atomorder -- a string of the form 'atom1,atom2,atom3' where 'atom1', 'atom2' and 'atom3' are
2008 respectively the MMTK atom names of the atoms in the way they should be ordered.
2009 * ac -- the output NetCDF file name. A CDL version of this file will also be generated with the '.cdl' extension
2010 instead of the '.nc' extension.
2011 * pyroserver -- a string specifying if Pyro will be used and how to run the analysis.
2012
2013 Running modes:
2014
2015 - To run the analysis do: a.runAnalysis() where a is the analysis object.
2016 - To estimate the analysis do: a.estimateAnalysis() where a is the analysis object.
2017 - To save the analysis to 'file' file name do: a.saveAnalysis(file) where a is the analysis object.
2018 """
2019
2020
2021 inputParametersNames = ('trajectory',\
2022 'timeinfo',\
2023 'triplet',\
2024 'atomorder',\
2025 'ac',\
2026 'pyroserver',\
2027 )
2028
2029 shortName = 'AC'
2030
2031
2032 canBeEstimated = True
2033
2035 """The constructor. Insures that the class can not be instanciated directly from here.
2036 """
2037 raise Error('This class can not be instanciated directly.')
2038
2040 """Initializes the analysis (e.g. parses and checks input parameters, set some variables ...).
2041 """
2042
2043
2044 self.parseInputParameters()
2045
2046 self.buildTimeInfo()
2047
2048 self.triplet = self.groupSelection(self.universe, self.triplet)
2049 self.triplet = [g for g in self.triplet if g.numberOfAtoms() == 3]
2050
2051 if self.atomOrder is None:
2052 self.triplet = [sorted(g, key = operator.attrgetter('name')) for g in self.triplet]
2053
2054 else:
2055 orderedTriplets = []
2056 for triplet in self.triplet:
2057 gr = []
2058 for atName in self.atomOrder:
2059 found = False
2060 for at in triplet:
2061 if at.name == atName:
2062 found = True
2063 gr.append(at)
2064 if len(gr) == 3:
2065 orderedTriplets.append(gr)
2066 gr = []
2067 else:
2068 if not found:
2069 raise Error('No atom corresponding %s could be found.' % atName)
2070
2071 self.triplet = orderedTriplets
2072
2073 if len(self.triplet) == 0:
2074 raise Error('Your group selection does not contain any triplet.')
2075
2076 self.group = self.triplet
2077
2078
2079 self.nTriplets = self.nGroups = len(self.triplet)
2080
2081 self.preLoadTrajectory('atom')
2082
2083 self.angCorr1 = N.zeros((self.nTriplets,self.nFrames),N.Float)
2084 self.angCorr2 = N.zeros((self.nTriplets,self.nFrames),N.Float)
2085 self.angCorr3 = N.zeros((self.nTriplets,self.nFrames),N.Float)
2086
2087 if self.pyroServer != 'monoprocessor':
2088
2089 delattr(self, 'trajectory')
2090
2091 - def calc(self, tripletIndex, trajname):
2092 """Calculates the contribution for one group.
2093
2094 @param tripletIndex: the index of the triplet in |self.triplet| list.
2095 @type tripletIndex: integer.
2096
2097 @param trajname: the name of the trajectory file name.
2098 @type trajname: string
2099 """
2100
2101 at1, at2, at3 = self.triplet[tripletIndex]
2102
2103 if self.preLoadedTraj is None:
2104 if self.pyroServer == 'monoprocessor':
2105 t = self.trajectory
2106
2107 else:
2108
2109 t = Trajectory(None, trajname, 'r')
2110
2111 at1Traj = t.readParticleTrajectory(at1, first = self.first, last = self.last, skip = self.skip)
2112 at2Traj = t.readParticleTrajectory(at2, first = self.first, last = self.last, skip = self.skip)
2113 at3Traj = t.readParticleTrajectory(at3, first = self.first, last = self.last, skip = self.skip)
2114
2115 else:
2116 at1Traj = self.preLoadedTraj[at1]
2117 at2Traj = self.preLoadedTraj[at2]
2118 at3Traj = self.preLoadedTraj[at3]
2119
2120 v1 = N.zeros((self.nFrames,3),N.Float)
2121 v2 = N.zeros((self.nFrames,3),N.Float)
2122 v3 = N.zeros((self.nFrames,3),N.Float)
2123
2124 for frameIndex in range(self.nFrames):
2125
2126 tempv1 = self.universe.distanceVector(at1Traj[frameIndex], at2Traj[frameIndex]).normal()
2127 tempv2 = self.universe.distanceVector(at1Traj[frameIndex], at3Traj[frameIndex]).normal()
2128
2129 tv1 = tempv1 + tempv2
2130 tv1 = tv1.normal()
2131 tv3 = tv1.cross(tempv2).normal()
2132
2133 if tv3*at1Traj[frameIndex] < 0:
2134 tv3 = -tv3
2135
2136 tv2 = tv3.cross(tv1).normal()
2137
2138 v1[frameIndex,:] = N.array(tv1)
2139 v2[frameIndex,:] = N.array(tv2)
2140 v3[frameIndex,:] = N.array(tv3)
2141
2142 return v1, v2, v3
2143
2144 - def combine(self, tripletIndex, x):
2145 """
2146 """
2147
2148
2149 self.angCorr1[tripletIndex,:] = correlation(x[0])
2150 self.angCorr2[tripletIndex,:] = correlation(x[1])
2151 self.angCorr3[tripletIndex,:] = correlation(x[2])
2152
2153
2155 """Finalizes the calculations (e.g. averaging the total term, output files creations ...).
2156 """
2157
2158 outputFile = NetCDFFile(self.outputAC, 'w')
2159 outputFile.title = self.__class__.__name__
2160 outputFile.jobinfo = self.information + '\nOutput file written on: %s\n\n' % asctime()
2161
2162 outputFile.createDimension('TRIPLET', self.nTriplets)
2163 outputFile.createDimension('TIME', self.nFrames)
2164
2165 TIMES = outputFile.createVariable('time', N.Float, ('TIME',))
2166 TIMES[:] = self.times
2167 TIMES.units = 'ps'
2168
2169 GROUPID = outputFile.createVariable('triplet', N.Int, ('TRIPLET',))
2170 GROUPID = 1 + N.arange(self.nTriplets)
2171
2172 AC1 = outputFile.createVariable('ac1_by_triplet', N.Float, ('TRIPLET','TIME'))
2173 AC1[:,:] = self.angCorr1[:,:]
2174
2175 AC2 = outputFile.createVariable('ac2_by_triplet', N.Float, ('TRIPLET','TIME'))
2176 AC2[:,:] = self.angCorr2[:,:]
2177
2178 AC3 = outputFile.createVariable('ac3_by_triplet', N.Float, ('TRIPLET','TIME'))
2179 AC3[:,:] = self.angCorr3[:,:]
2180
2181 AC1AVG = outputFile.createVariable('ac1', N.Float, ('TIME',))
2182 AC1AVG[:] = self.angCorr1.sum(0)/float(self.nTriplets)
2183
2184 AC2AVG = outputFile.createVariable('ac2', N.Float, ('TIME',))
2185 AC2AVG[:] = self.angCorr2.sum(0)/float(self.nTriplets)
2186
2187 AC3AVG = outputFile.createVariable('ac3', N.Float, ('TIME',))
2188 AC3AVG[:] = self.angCorr3.sum(0)/float(self.nTriplets)
2189
2190 asciiVar = sorted(outputFile.variables.keys())
2191
2192 outputFile.close()
2193
2194 self.toPlot = {'netcdf' : self.outputAC, 'xVar' : 'time', 'yVar' : 'ac3'}
2195
2196
2197 convertNetCDFToASCII(inputFile = self.outputAC,\
2198 outputFile = os.path.splitext(self.outputAC)[0] + '.cdl',\
2199 variables = asciiVar)
2200
2201
2202
2203
2204 -class RigidBodyTrajectory(Analysis):
2205 """Sets up a Rigid Body Trajectory analysis.
2206
2207 A Subclass of nMOLDYN.Analysis.Analysis.
2208
2209 Constructor: RigidBodyTrajectory(|parameters| = None)
2210
2211 Arguments:
2212
2213 - |parameters| -- a dictionnary of the input parameters, or 'None' to set up the analysis without parameters.
2214 * trajectory -- a trajectory file name or an instance of MMTK.Trajectory.Trajectory class.
2215 * timeinfo -- a string of the form 'first:last:step' where 'first' is an integer specifying the first frame
2216 number to consider, 'last' is an integer specifying the last frame number to consider and
2217 'step' is an integer specifying the step number between two frames.
2218 * referenceframe -- an integer in [1,len(trajectory)] specifying which frame should be the reference.
2219 * stepwiserbt -- a string being one of 'Yes' or 'No' specifying whether the reference frame for frame i should be
2220 the frame i - 1 ('Yes') or should be a fixed frame defined with |referenceframe| ('No').
2221 * group -- a selection string specifying the groups of atoms on which the rigid body trajectory will be defined.
2222 (each group being a rigid body).
2223 * rbt -- the output NetCDF file name.
2224 * pyroserver -- a string specifying if Pyro will be used and how to run the analysis.
2225
2226 Running modes:
2227
2228 - To run the analysis do: a.runAnalysis() where a is the analysis object.
2229 - To estimate the analysis do: a.estimateAnalysis() where a is the analysis object.
2230 - To save the analysis to 'file' file name do: a.saveAnalysis(file) where a is the analysis object.
2231 """
2232
2233
2234 inputParametersNames = ('trajectory',\
2235 'timeinfo',\
2236 'referenceframe',\
2237 'removetranslation',\
2238 'stepwiserbt',\
2239 'group',\
2240 'rbt',\
2241 'pyroserver',\
2242 )
2243
2244 shortName = 'RBT'
2245
2246
2247 canBeEstimated = True
2248
2249 - def __init__(self):
2250 """The constructor. Insures that the class can not be instanciated directly from here.
2251 """
2252 raise Error('This class can not be instanciated directly.')
2253
2254 - def initialize(self):
2255 """Initializes the analysis (e.g. parses and checks input parameters, set some variables ...).
2256 """
2257
2258
2259 self.parseInputParameters()
2260
2261 if self.pyroServer != 'monoprocessor':
2262 if self.statusBar is not None:
2263 LogMessage('warning', 'The job status bar will be desactivated for compatibility with Pyro module.')
2264 self.statusBar = None
2265
2266 self.buildTimeInfo()
2267
2268 if (self.referenceFrame < 0) or (self.referenceFrame > len(self.trajectory)):
2269 self.referenceFrame = 0
2270 LogMessage('warning', 'The reference frame must be an integer in [1,%s].\n\
2271 It will be set to 1 for the running analysis.' % len(self.trajectory), ['console'])
2272
2273 self.group = self.groupSelection(self.universe, self.group)
2274 self.nGroups = len(self.group)
2275 self.groupsSize = [g.numberOfAtoms() for g in self.group]
2276
2277
2278 overlap = set()
2279 for g in self.group:
2280 overlap = overlap.union(set(g.atomList()))
2281
2282
2283 if len(overlap) != sum(self.groupsSize):
2284 raise Error('There must not be any common atom between the groups.')
2285
2286
2287 self.refConfig = self.trajectory.configuration[self.referenceFrame]
2288
2289
2290 self.rbtQuaternions = N.zeros((self.nFrames, self.nGroups, 4), typecode = N.Float)
2291 self.rbtCMS = N.zeros((self.nFrames, self.nGroups, 3), typecode = N.Float)
2292 self.rbtFit = N.zeros((self.nFrames, self.nGroups), typecode = N.Float)
2293
2294
2295
2296 self.rbtTraj = {}
2297
2298 if self.pyroServer != 'monoprocessor':
2299
2300 delattr(self, 'trajectory')
2301
2302 - def calc(self, groupIndex, trajname):
2303 """Calculates the contribution for one group.
2304
2305 @param groupIndex: the index of the group in |self.group| list.
2306 @type groupIndex: integer.
2307
2308 @param trajname: the name of the trajectory file name.
2309 @type trajname: string
2310 """
2311
2312 if self.pyroServer == 'monoprocessor':
2313 t = self.trajectory
2314
2315 else:
2316
2317 t = Trajectory(None, trajname, 'r')
2318
2319 group = self.group[groupIndex]
2320
2321
2322 quaternions = N.zeros((self.nFrames, 4), typecode = N.Float)
2323 CMS = N.zeros((self.nFrames, 3), typecode = N.Float)
2324 fit = N.zeros((self.nFrames,), typecode = N.Float)
2325
2326
2327 if self.stepwiseRBT:
2328
2329
2330
2331 for frameIndex in range(self.nFrames):
2332 if frameIndex == 0:
2333 movingRefConfig = t.configuration[self.frameIndexes[frameIndex]]
2334 else:
2335 movingRefConfig = t.configuration[self.frameIndexes[frameIndex-1]]
2336
2337
2338 rbt = t.readRigidBodyTrajectory(group,\
2339 first = self.frameIndexes[frameIndex],\
2340 last = self.frameIndexes[frameIndex] + 1,\
2341 skip = 1,\
2342 reference = movingRefConfig)
2343
2344
2345 quaternions[frameIndex,:] = copy.copy(rbt.quaternions)
2346 CMS[frameIndex,:] = copy.copy(rbt.cms)
2347 fit[frameIndex] = copy.copy(rbt.fit)
2348
2349
2350
2351 else:
2352
2353 rbt = t.readRigidBodyTrajectory(group,\
2354 first = self.first,\
2355 last = self.last,\
2356 skip = self.skip,\
2357 reference = self.refConfig)
2358
2359
2360 quaternions[:,:] = copy.copy(rbt.quaternions)
2361 CMS[:,:] = copy.copy(rbt.cms)
2362 fit[:] = copy.copy(rbt.fit)
2363
2364 return quaternions, CMS, fit
2365
2366 - def combine(self, groupIndex, x):
2367 """
2368 """
2369
2370 self.rbtQuaternions[:,groupIndex,:] = copy.copy(x[0])
2371 self.rbtCMS[:,groupIndex,:] = copy.copy(x[1])
2372 self.rbtFit[:,groupIndex] = copy.copy(x[2])
2373
2374 group = self.group[groupIndex]
2375
2376
2377
2378 centerOfMass = group.centerOfMass(self.refConfig)
2379
2380 atomicRBT = N.zeros((self.nFrames,3), typecode = N.Float)
2381
2382
2383 for atom in group:
2384
2385 xyz = self.refConfig[atom] - centerOfMass
2386
2387
2388 for frameIndex in range(self.nFrames):
2389
2390 transfo = Quaternion(self.rbtQuaternions[frameIndex,groupIndex,:]).asRotation()
2391
2392 if self.removeTranslation:
2393
2394 transfo = Translation(centerOfMass)*transfo
2395
2396
2397 else:
2398
2399 transfo = Translation(Vector(self.rbtCMS[frameIndex,groupIndex,:]))*transfo
2400
2401
2402 atomicRBT[frameIndex,:] = transfo(Vector(xyz))
2403
2404 self.rbtTraj[atom.index] = copy.copy(atomicRBT)
2405
2406 - def finalize(self):
2407 """Finalizes the calculations (e.g. averaging the total term, output files creations ...).
2408 """
2409
2410 selection = Collection(self.group)
2411
2412
2413 outputFile = Trajectory(selection, self.outputRBT, 'w')
2414
2415 if self.pyroServer == 'monoprocessor':
2416 t = self.trajectory
2417
2418 else:
2419
2420 t = Trajectory(None, self.trajectoryFilename, 'r')
2421
2422
2423 snapshot = SnapshotGenerator(self.universe, actions = [TrajectoryOutput(outputFile, ["all"], 0, None, 1)])
2424
2425
2426 for frameIndex in range(self.nFrames):
2427 for atom in selection:
2428 atom.setPosition(self.rbtTraj[atom.index][frameIndex,:])
2429 snapshot(data = {'time' : self.times[frameIndex]})
2430
2431 outputFile.close()
2432
2433 outputFile = NetCDFFile(self.outputRBT, 'a')
2434
2435 outputFile.title = self.__class__.__name__
2436 outputFile.jobinfo = self.information + '\nOutput file written on: %s\n\n' % asctime()
2437
2438
2439
2440 for groupIndex in range(self.nGroups):
2441 outputFile.jobinfo += 'Group %s: %s\n' % (groupIndex+1, [atom.index for atom in self.group[groupIndex]])
2442
2443 outputFile.createDimension('number_of_groups', self.nGroups)
2444 outputFile.createDimension('quaternion_length',4)
2445
2446
2447 quaternions = outputFile.createVariable('quaternion', N.Float32, ('step_number','number_of_groups','quaternion_length'))
2448
2449
2450 cms = outputFile.createVariable('cms', N.Float32, ('step_number','number_of_groups','xyz'))
2451
2452
2453 fit = outputFile.createVariable('fit', N.Float32, ('step_number','number_of_groups'))
2454
2455
2456 for groupIndex in range(self.nGroups):
2457 quaternions[:,groupIndex,:] = copy.copy(self.rbtQuaternions[:,groupIndex,:].astype(N.Float32))
2458 cms[:,groupIndex,:] = copy.copy(self.rbtCMS[:,groupIndex,:].astype(N.Float32))
2459 fit[:,groupIndex] = copy.copy(self.rbtFit[:,groupIndex].astype(N.Float32))
2460
2461 outputFile.close()
2462
2463 self.toPlot = None
2464
2465
2466
2467
2469 """Sets up a Reorientational Correlation Function analysis.
2470
2471 A Subclass of nMOLDYN.Analysis.Analysis.
2472
2473 Constructor: ReorientationalCorrelationFunction(|parameters| = None)
2474
2475 Arguments:
2476
2477 - |parameters| -- a dictionnary of the input parameters, or 'None' to set up the analysis without parameters.
2478 * trajectory -- a trajectory file name or an instance of MMTK.Trajectory.Trajectory class.
2479 * timeinfo -- a string of the form 'first:last:step' where 'first' is an integer specifying the first frame
2480 number to consider, 'last' is an integer specifying the last frame number to consider and
2481 'step' is an integer specifying the step number between two frames.
2482 * referenceframe -- an integer in [1,len(trajectory)] specifying which frame should be the reference.
2483 * stepwiserbt -- a string being one of 'Yes' or 'No' specifying whether the reference frame for frame i should be
2484 the frame i - 1 ('Yes') or should be a fixed frame defined with |referenceframe| ('No').
2485 * wignerindexes -- a string of the form 'j,m,n' where 'j', 'm' and 'n' are respectively the j, m and n indexes of the
2486 Wigner function Djmn.
2487 * group -- a selection string specifying the groups of atoms on which the rigid body trajectory will be defined.
2488 (each group being a rigid body).
2489 * rcf -- the output NetCDF file name. A CDL version of this file will also be generated with the '.cdl' extension
2490 instead of the '.nc' extension.
2491 * pyroserver -- a string specifying if Pyro will be used and how to run the analysis.
2492
2493 Running modes:
2494
2495 - To run the analysis do: a.runAnalysis() where a is the analysis object.
2496 - To estimate the analysis do: a.estimateAnalysis() where a is the analysis object.
2497 - To save the analysis to 'file' file name do: a.saveAnalysis(file) where a is the analysis object.
2498 """
2499
2500
2501 inputParametersNames = ('trajectory',\
2502 'timeinfo',\
2503 'referenceframe',\
2504 'stepwiserbt',\
2505 'wignerindexes',\
2506 'group',\
2507 'rcf',\
2508 'pyroserver',\
2509 )
2510
2511 shortName = 'RCF'
2512
2513
2514 canBeEstimated = True
2515
2517 """The constructor. Insures that the class can not be instanciated directly from here.
2518 """
2519 raise Error('This class can not be instanciated directly.')
2520
2522 """Initializes the analysis (e.g. parses and checks input parameters, set some variables ...).
2523 """
2524
2525
2526 self.parseInputParameters()
2527
2528 if self.pyroServer != 'monoprocessor':
2529 if self.statusBar is not None:
2530 LogMessage('warning', 'The job status bar will be desactivated for compatibility with Pyro module.')
2531 self.statusBar = None
2532
2533 self.buildTimeInfo()
2534
2535 if (self.referenceFrame < 0) or (self.referenceFrame > len(self.trajectory)):
2536 self.referenceFrame = 0
2537 LogMessage('warning', 'The reference frame must be an integer in [1,%s].\n\
2538 It will be set to 1 for the running analysis.' % len(self.trajectory), ['console'])
2539
2540 self.group = self.groupSelection(self.universe, self.group)
2541 self.nGroups = len(self.group)
2542 self.groupsSize = [len(g) for g in self.group]
2543 self.nObjects = sum(self.groupsSize)
2544
2545
2546 self.refConfig = self.trajectory.configuration[self.referenceFrame]
2547
2548
2549 self.RCF = N.zeros((self.nFrames), typecode = N.Float)
2550
2551 if self.pyroServer != 'monoprocessor':
2552
2553 delattr(self, 'trajectory')
2554
2555 - def calc(self, groupIndex, trajname):
2556 """Calculates the contribution for one group.
2557
2558 @param groupIndex: the index of the group in |self.group| list.
2559 @type groupIndex: integer.
2560
2561 @param trajname: the name of the trajectory file name.
2562 @type trajname: string
2563 """
2564
2565 if self.pyroServer == 'monoprocessor':
2566 t = self.trajectory
2567
2568 else:
2569
2570 t = Trajectory(None, trajname, 'r')
2571
2572 group = self.group[groupIndex]
2573
2574 j, m, n = self.wignerIndexes
2575
2576
2577 rbtQuaternions = N.zeros((self.nFrames,4), typecode = N.Float)
2578
2579
2580 if self.stepwiseRBT:
2581
2582
2583 for frameIndex in range(self.nFrames):
2584 if frameIndex == 0:
2585 movingRefConfig = t.configuration[self.frameIndexes[frameIndex]]
2586 else:
2587 movingRefConfig = t.configuration[self.frameIndexes[frameIndex-1]]
2588
2589
2590 rbt = t.readRigidBodyTrajectory(group,\
2591 first = self.frameIndexes[frameIndex],\
2592 last = self.frameIndexes[frameIndex] + 1,\
2593 skip = 1,\
2594 reference = movingRefConfig)
2595
2596 rbtQuaternions[frameIndex,:] = copy.copy(rbt.quaternions)
2597 rbtCMS[frameIndex,:] = copy.copy(rbt.cms)
2598
2599
2600
2601 else:
2602
2603 rbt = t.readRigidBodyTrajectory(group,\
2604 first = self.first,\
2605 last = self.last,\
2606 skip = self.skip,\
2607 reference = self.refConfig)
2608
2609 rbtQuaternions = rbt.quaternions
2610
2611
2612 c1 = N.sqrt(((2.0*j+1)/(4.0*N.pi)))
2613
2614 quat2 = N.zeros(rbtQuaternions.shape, typecode = N.Complex)
2615
2616 quat2[:,0] = rbtQuaternions[:,0] + 1j*rbtQuaternions[:,3]
2617
2618 quat2[:,2] = rbtQuaternions[:,2] + 1j*rbtQuaternions[:, 1]
2619
2620 quat2[:,1] = N.conjugate(quat2[:,0])
2621
2622 quat2[:,3] = N.conjugate(quat2[:,2])
2623
2624 pp = preparePP(j, m, n)
2625 Djmn = N.add.reduce(N.multiply.reduce(quat2[:,N.NewAxis,:]**pp[0][N.NewAxis,:,:],-1)*pp[1],1)
2626
2627 if m == n:
2628 Djnm = Djmn
2629 else:
2630 pp = self.preparePP(j, n, m)
2631 Djnm = N.add.reduce(N.multiply.reduce(quat2[:,N.NewAxis,:]**pp[0][N.NewAxis,:,:],-1)*pp[1],1)
2632
2633 Djmn = Djmn*c1
2634 Djnm = Djnm*c1
2635
2636 return Djmn, Djnm
2637
2638 - def combine(self, groupIndex, x):
2639 """
2640 """
2641
2642 N.add(self.RCF,correlation(x[0],x[1]),self.RCF)
2643
2645 """Finalizes the calculations (e.g. averaging the total term, output files creations ...).
2646 """
2647
2648 self.RCF = 4.0*N.pi*self.RCF/self.nObjects
2649
2650
2651 outputFile = NetCDFFile(self.outputRCF, 'w')
2652 outputFile.title = self.__class__.__name__
2653 outputFile.jobinfo = self.information + '\nOutput file written on: %s\n\n' % asctime()
2654
2655
2656
2657 for groupIndex in range(self.nGroups):
2658 outputFile.jobinfo += 'Group %s: %s\n' % (groupIndex+1,[atom.index for atom in self.group[groupIndex]])
2659
2660
2661 outputFile.createDimension('TIMES', self.nFrames)
2662
2663
2664
2665 TIMES = outputFile.createVariable('time', N.Float, ('TIMES',))
2666 TIMES[:] = self.times
2667 TIMES.units = 'ps'
2668
2669
2670 RCF = outputFile.createVariable('rcf', N.Float, ('TIMES',))
2671 RCF[:] = self.RCF[:]
2672 RCF.units = 'unitless'
2673
2674 outputFile.close()
2675
2676 self.toPlot = {'netcdf' : self.outputRCF, 'xVar' : 'time', 'yVar' : 'rcf'}
2677
2678
2679 convertNetCDFToASCII(inputFile = self.outputRCF,\
2680 outputFile = os.path.splitext(self.outputRCF)[0] + '.cdl',\
2681 variables = ['time', 'rcf'])
2682
2684 """An intermediate class used by |AngularVelocityAutoCorrelationFunction| and |AngularDensityOfStates| classes."""
2685
2689
2691
2692 res = N.zeros((len(data),4,4), N.Float)
2693
2694 res[:,0,0]= data[:,0]
2695 res[:,1,1]= data[:,0]
2696 res[:,2,2]= data[:,0]
2697 res[:,3,3]= data[:,0]
2698
2699 res[:,0,1]= data[:,1]
2700 res[:,1,0]= -data[:,1]
2701 res[:,2,3]= data[:,1]
2702 res[:,3,2]= -data[:,1]
2703
2704 res[:,0,2]= data[:,2]
2705 res[:,1,3]= -data[:,2]
2706 res[:,2,0]= -data[:,2]
2707 res[:,3,1]= data[:,2]
2708
2709 res[:,0,3]= data[:,3]
2710 res[:,1,2]= data[:,3]
2711 res[:,2,1]= -data[:,3]
2712 res[:,3,0]= -data[:,3]
2713
2714 return res
2715
2717 """Computes the Angular Velocity Function for a group |g| (a MMTK Collection)."""
2718
2719 rbtQuaternions = N.zeros((self.nFrames, 4), typecode = N.Float)
2720
2721
2722 refConfig = t.configuration[self.referenceFrame]
2723
2724
2725 if self.stepwiseRBT:
2726
2727
2728 for i in range(self.nFrames):
2729 if i == 0:
2730 movingRefConfig = t.configuration[self.frameIndexes[i]]
2731 else:
2732 movingRefConfig = t.configuration[self.frameIndexes[i-1]]
2733
2734
2735 rbt = t.readRigidBodyTrajectory(g,\
2736 first = self.frameIndexes[i],\
2737 last = self.frameIndexes[i] + 1,\
2738 skip = 1,\
2739 reference = movingRefConfig)
2740
2741 rbtQuaternions[i,:] = copy.copy(rbt.quaternions)
2742
2743
2744
2745 else:
2746
2747 rbt = t.readRigidBodyTrajectory(g,\
2748 first = self.first,\
2749 last = self.last,\
2750 skip = self.skip,\
2751 reference = refConfig)
2752 rbtQuaternions[:,:] = copy.copy(rbt.quaternions)
2753
2754
2755 quaternions_dot = N.zeros((self.nFrames,4), N.Float)
2756 for i in range(4):
2757 quaternions_dot[:,i] = differentiate(rbtQuaternions[:,i], self.differentiation, self.dt)
2758
2759 q = self.qMatrix(rbtQuaternions)
2760
2761 angvel = 2.0*N.add.reduce(q*quaternions_dot[:,N.NewAxis,:], -1)
2762
2763
2764 return angvel[:,1:]
2765
2766
2767
2768
2770 """Sets up an Angular Velocity AutoCorrelation Function analysis.
2771
2772 A Subclass of nMOLDYN.Analysis.Analysis.
2773
2774 Constructor: AngularVelocityAutoCorrelationFunction(|parameters| = None)
2775
2776 Arguments:
2777
2778 - |parameters| -- a dictionnary of the input parameters, or 'None' to set up the analysis without parameters.
2779 * trajectory -- a trajectory file name or an instance of MMTK.Trajectory.Trajectory class.
2780 * timeinfo -- a string of the form 'first:last:step' where 'first' is an integer specifying the first frame
2781 number to consider, 'last' is an integer specifying the last frame number to consider and
2782 'step' is an integer specifying the step number between two frames.
2783 * differentiation -- an integer in [0,5] specifying the order of the differentiation used to get the velocities
2784 out of the coordinates. 0 means that the velocities are already present in the trajectory loaded
2785 for analysis.
2786 * projection -- a string of the form 'vx,vy,vz' specifying the vector along which the analysis
2787 will be computed. 'vx', 'vy', and 'vz' are floats specifying respectively the x, y and z value
2788 of that vector.
2789 * referenceframe -- an integer in [1,len(trajectory)] specifying which frame should be the reference.
2790 * stepwiserbt -- a string being one of 'Yes' or 'No' specifying whether the reference frame for frame i should be
2791 the frame i - 1 ('Yes') or should be a fixed frame defined with |referenceframe| ('No').
2792 * group -- a selection string specifying the groups of atoms on which the rigid body trajectory will be defined.
2793 (each group being a rigid body).
2794 * avacf -- the output NetCDF file name. A CDL version of this file will also be generated with the '.cdl' extension
2795 instead of the '.nc' extension.
2796 * pyroserver -- a string specifying if Pyro will be used and how to run the analysis.
2797
2798 Running modes:
2799
2800 - To run the analysis do: a.runAnalysis() where a is the analysis object.
2801 - To estimate the analysis do: a.estimateAnalysis() where a is the analysis object.
2802 - To save the analysis to 'file' file name do: a.saveAnalysis(file) where a is the analysis object.
2803 """
2804
2805
2806 inputParametersNames = ('trajectory',\
2807 'timeinfo',\
2808 'differentiation',\
2809 'projection',\
2810 'referenceframe',\
2811 'stepwiserbt',\
2812 'group',\
2813 'avacf',\
2814 'pyroserver',\
2815 )
2816
2817 shortName = 'AVACF'
2818
2819
2820 canBeEstimated = True
2821
2823 """The constructor. Insures that the class can not be instanciated directly from here.
2824 """
2825 raise Error('This class can not be instanciated directly.')
2826
2828 """Initializes the analysis (e.g. parses and checks input parameters, set some variables ...).
2829 """
2830
2831
2832 self.parseInputParameters()
2833
2834 if self.pyroServer != 'monoprocessor':
2835 if self.statusBar is not None:
2836 LogMessage('warning', 'The job status bar will be desactivated for compatibility with Pyro module.')
2837 self.statusBar = None
2838
2839 self.buildTimeInfo()
2840
2841 if (self.referenceFrame < 0) or (self.referenceFrame > len(self.trajectory)):
2842 self.referenceFrame = 0
2843 LogMessage('warning', 'The reference frame must be an integer in [1,%s].\n\
2844 It will be set to 1 for the running analysis.' % len(self.trajectory), ['console'])
2845
2846 self.group = self.groupSelection(self.universe, self.group)
2847 self.nGroups = len(self.group)
2848 self.groupsSize = [len(g) for g in self.group]
2849
2850
2851 overlap = set()
2852 for g in self.group:
2853 overlap = overlap.union(set(g.atomList()))
2854
2855
2856 if len(overlap) != sum(self.groupsSize):
2857 raise Error('There must not be any common atom between the groups.')
2858
2859 self.AVACF = {'total' : N.zeros((self.nFrames), typecode = N.Float)}
2860
2861 if self.pyroServer != 'monoprocessor':
2862
2863 delattr(self, 'trajectory')
2864
2865 - def calc(self, groupIndex, trajname):
2866 """Calculates the contribution for one group.
2867
2868 @param groupIndex: the index of the group in |self.group| list.
2869 @type groupIndex: integer.
2870
2871 @param trajname: the name of the trajectory file name.
2872 @type trajname: string
2873 """
2874
2875 if self.pyroServer == 'monoprocessor':
2876 t = self.trajectory
2877
2878 else:
2879
2880 t = Trajectory(None, trajname, 'r')
2881
2882 group = self.group[groupIndex]
2883
2884 AngularVelocity.__init__(self)
2885
2886
2887 angvel = self.getAngularVelocity(t, group)
2888
2889 return angvel
2890
2891 - def combine(self, groupIndex, x):
2892 """
2893 """
2894
2895
2896 if self.projection is None:
2897 self.AVACF['G%d' % (groupIndex+1,)] = correlation(x)/3.0
2898 else:
2899 projectedAngVel = N.dot(x, self.projection)
2900 self.AVACF['G%d' % (groupIndex+1,)] = correlation(projectedAngVel)
2901
2902 self.AVACF['G%d' % (groupIndex+1,)] /= self.AVACF['G%d' % (groupIndex+1,)][0]
2903
2905 """Finalizes the calculations (e.g. averaging the total term, output files creations ...).
2906 """
2907
2908
2909 outputFile = NetCDFFile(self.outputAVACF, 'w')
2910 outputFile.title = self.__class__.__name__
2911 outputFile.jobinfo = self.information + '\nOutput file written on: %s\n\n' % asctime()
2912
2913
2914
2915 for groupIndex in range(self.nGroups):
2916 outputFile.jobinfo += 'Group %d: %s\n' % (groupIndex+1,[atom.index for atom in self.group[groupIndex]])
2917
2918
2919 outputFile.createDimension('NFRAMES', self.nFrames)
2920
2921
2922
2923 TIMES = outputFile.createVariable('time', N.Float, ('NFRAMES',))
2924 TIMES[:] = self.times[:]
2925 TIMES.units = 'ps'
2926
2927 for k in self.AVACF.keys():
2928 if k == 'total':
2929 continue
2930
2931 AVACF = outputFile.createVariable('avacf-' + k, N.Float, ('NFRAMES',))
2932 AVACF[:] = self.AVACF[k][:]
2933 AVACF.units = 'rad^2*ps^-2'
2934
2935 N.add(self.AVACF['total'], self.AVACF[k], self.AVACF['total'])
2936
2937 self.AVACF['total'] /= self.nGroups
2938
2939 AVACF = outputFile.createVariable('avacf-total', N.Float, ('NFRAMES',))
2940 AVACF[:] = self.AVACF['total'][:]
2941 AVACF.units = 'rad^2*ps^-2'
2942
2943 outputFile.close()
2944
2945 self.toPlot = {'netcdf' : self.outputAVACF, 'xVar' : 'time', 'yVar' : 'avacf-total'}
2946
2947
2948 convertNetCDFToASCII(inputFile = self.outputAVACF,\
2949 outputFile = os.path.splitext(self.outputAVACF)[0] + '.cdl',\
2950 variables = ['time', 'avacf-total'])
2951
2952
2953
2954
2956 """Sets up an Angular Density Of States analysis.
2957
2958 A Subclass of nMOLDYN.Analysis.Analysis.
2959
2960 Constructor: AngularDensityOfStates(|parameters| = None)
2961
2962 Arguments:
2963
2964 - |parameters| -- a dictionnary of the input parameters, or 'None' to set up the analysis without parameters.
2965 * trajectory -- a trajectory file name or an instance of MMTK.Trajectory.Trajectory class.
2966 * timeinfo -- a string of the form 'first:last:step' where 'first' is an integer specifying the first frame
2967 number to consider, 'last' is an integer specifying the last frame number to consider and
2968 'step' is an integer specifying the step number between two frames.
2969 * differentiation -- an integer in [0,5] specifying the order of the differentiation used to get the velocities
2970 out of the coordinates. 0 means that the velocities are already present in the trajectory loaded
2971 for analysis.
2972 * projection -- a string of the form 'vx,vy,vz' specifying the vector along which the analysis
2973 will be computed. 'vx', 'vy', and 'vz' are floats specifying respectively the x, y and z value
2974 of that vector.
2975 * referenceframe -- an integer in [1,len(trajectory)] specifying which frame should be the reference.
2976 * stepwiserbt -- a string being one of 'Yes' or 'No' specifying whether the reference frame for frame i should be
2977 the frame i - 1 ('Yes') or should be a fixed frame defined with |referenceframe| ('No').
2978 * fftwindow -- a float in ]0.0,100.0[ specifying the width of the gaussian, in percentage of the trajectory length
2979 that will be used in the smoothing procedure.
2980 * group -- a selection string specifying the groups of atoms on which the rigid body trajectory will be defined.
2981 (each group being a rigid body).
2982 * ados -- the output NetCDF file name. A CDL version of this file will also be generated with the '.cdl' extension
2983 instead of the '.nc' extension.
2984 * pyroserver -- a string specifying if Pyro will be used and how to run the analysis.
2985
2986 Running modes:
2987
2988 - To run the analysis do: a.runAnalysis() where a is the analysis object.
2989 - To estimate the analysis do: a.estimateAnalysis() where a is the analysis object.
2990 - To save the analysis to 'file' file name do: a.saveAnalysis(file) where a is the analysis object.
2991 """
2992
2993
2994 inputParametersNames = ('trajectory',\
2995 'timeInfo',\
2996 'differentiation',\
2997 'projection',\
2998 'referenceframe',\
2999 'stepwiserbt',\
3000 'fftwindow',\
3001 'group',\
3002 'ados',\
3003 'pyroserver',\
3004 )
3005
3006 shortName = 'ADOS'
3007
3008
3009 canBeEstimated = True
3010
3012 """The constructor. Insures that the class can not be instanciated directly from here.
3013 """
3014 raise Error('This class can not be instanciated directly.')
3015
3017 """Initializes the analysis (e.g. parses and checks input parameters, set some variables ...).
3018 """
3019
3020
3021 self.parseInputParameters()
3022
3023 if self.pyroServer != 'monoprocessor':
3024 if self.statusBar is not None:
3025 LogMessage('warning', 'The job status bar will be desactivated for compatibility with Pyro module.')
3026 self.statusBar = None
3027
3028 self.buildTimeInfo()
3029
3030 if (self.referenceFrame < 0) or (self.referenceFrame > len(self.trajectory)):
3031 self.referenceFrame = 0
3032 LogMessage('warning', 'The reference frame must be an integer in [1,%s].\n\
3033 It will be set to 1 for the running analysis.' % len(self.trajectory), ['console'])
3034
3035 self.group = self.groupSelection(self.universe, self.group)
3036 self.nGroups = len(self.group)
3037 self.groupsSize = [len(g) for g in self.group]
3038
3039
3040 overlap = set()
3041 for g in self.group:
3042 overlap = overlap.union(set(g.atomList()))
3043
3044
3045 if len(overlap) != sum(self.groupsSize):
3046 raise Error('There must not be any common atom between the groups.')
3047
3048 self.AVACF = {'total' : N.zeros((self.nFrames), typecode = N.Float)}
3049 self.ADOS = {'total' : N.zeros((self.nFrames), typecode = N.Float)}
3050
3051 if self.pyroServer != 'monoprocessor':
3052
3053 delattr(self, 'trajectory')
3054
3055 - def calc(self, groupIndex, trajname):
3056 """Calculates the contribution for one group.
3057
3058 @param groupIndex: the index of the group in |self.group| list.
3059 @type groupIndex: integer.
3060
3061 @param trajname: the name of the trajectory file name.
3062 @type trajname: string
3063 """
3064
3065 if self.pyroServer == 'monoprocessor':
3066 t = self.trajectory
3067
3068 else:
3069
3070 t = Trajectory(None, trajname, 'r')
3071
3072 group = self.group[groupIndex]
3073
3074 AngularVelocity.__init__(self)
3075
3076
3077 angvel = self.getAngularVelocity(t, group)
3078
3079 return angvel
3080
3081 - def combine(self, groupIndex, x):
3082 """
3083 """
3084
3085 if self.projection is None:
3086 self.AVACF['G%d' % (groupIndex+1,)] = correlation(x)/3.0
3087 else:
3088 projectedAngVel = N.dot(x, self.projection)
3089 self.AVACF['G%d' % (groupIndex+1,)] = correlation(projectedAngVel)
3090
3091 self.ADOS['G%d' % (groupIndex+1,)] = FFT(gaussianWindow(self.AVACF['G%d' % groupIndex+1], self.fftWindow)).real[:self.nFrames]
3092
3094 """Finalizes the calculations (e.g. averaging the total term, output files creations ...).
3095 """
3096
3097
3098 frequencies = N.arange(self.nFrames)/(2.0*self.nFrames*self.dt)
3099
3100
3101 outputFile = NetCDFFile(self.outputADOS, 'w')
3102 outputFile.title = self.__class__.__name__
3103 outputFile.jobinfo = self.information + '\nOutput file written on: %s\n\n' % asctime()
3104
3105
3106
3107 for groupIndex in range(self.nGroups):
3108 outputFile.jobinfo += 'Group %d: %s\n' % (groupIndex+1,[atom.index for atom in self.group[groupIndex]])
3109
3110
3111 outputFile.createDimension('NFRAMES', self.nFrames)
3112
3113
3114
3115 TIMES = outputFile.createVariable('time', N.Float, ('NFRAMES',))
3116 TIMES[:] = self.times[:]
3117 TIMES.units = 'ps'
3118
3119
3120
3121 FREQUENCIES = outputFile.createVariable('frequency', N.Float, ('NFRAMES',))
3122 FREQUENCIES[:] = frequencies[:]
3123 FREQUENCIES.units = 'THz'
3124
3125 for k in self.ADOS.keys():
3126 if k == 'total':
3127 continue
3128
3129 AVACF = outputFile.createVariable('avacf-' + k, N.Float, ('NFRAMES',))
3130 AVACF[:] = self.AVACF[k][:]
3131 AVACF.units = 'rad^2*ps^-2'
3132
3133 N.add(self.AVACF['total'], self.AVACF[k], self.AVACF['total'])
3134
3135 ADOS = outputFile.createVariable('ados-' + k, N.Float, ('NFRAMES',))
3136 ADOS[:] = self.ADOS[k][:]
3137 ADOS.units = 'rad^2*ps^-1'
3138
3139 N.add(self.ADOS['total'], self.ADOS[k], self.ADOS['total'])
3140
3141 AVACF = outputFile.createVariable('avacf-total', N.Float, ('NFRAMES',))
3142 AVACF[:] = self.AVACF['total'][:]
3143 AVACF.units = 'rad^2*ps^-2'
3144
3145 ADOS = outputFile.createVariable('ados-total', N.Float, ('NFRAMES',))
3146 ADOS[:] = self.ADOS['total'][:]
3147 ADOS.units = 'rad^2*ps^-1'
3148
3149 outputFile.close()
3150
3151 self.toPlot = {'netcdf' : self.outputADOS, 'xVar' : 'frequency', 'yVar' : 'ados-total'}
3152
3153
3154 convertNetCDFToASCII(inputFile = self.outputADOS,\
3155 outputFile = os.path.splitext(self.outputADOS)[0] + '.cdl',\
3156 variables = ['frequency', 'ados-total'])
3157