Package nMOLDYN :: Package Analysis :: Module Dynamics
[hide private]
[frames] | no frames]

Source Code for Module nMOLDYN.Analysis.Dynamics

   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  # The python distribution modules 
  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  # The ScientificPython modules 
  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  # The MMTK distribution modules 
  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  # The nMOLDYN modules 
  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  # MEAN-SQUARE DISPLACEMENT ANALYSIS 
  55  ##################################################################################### 
56 -class MeanSquareDisplacement(Analysis):
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 # Indicate whether this analysis can be estimated or not. 105 canBeEstimated = True 106
107 - def __init__(self):
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
112 - def initialize(self):
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 # MSD = 1D Numeric array. Contains the whole system MSD. 136 self.MSD = N.zeros((self.nFrames), N.Float) 137 138 if self.pyroServer != 'monoprocessor': 139 # The attribute trajectory is removed because it can not be pickled by Pyro. 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 # Load the whole trajectory set. 157 t = Trajectory(None, trajname, 'r') 158 159 # series = 2D Numeric array. The positions of the selected atom |at| from the first step to the 160 # last step with the selected step increment. 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 # dsq is the squared norm of the position for each time step 168 # dsq refers to DSQ(k) in the published algorithm 169 dsq = N.add.reduce(series * series,1) 170 else: 171 # if a projection vector is given, the trajectory is then projected onto the projection vector. 172 series = N.dot(series,self.projection) 173 dsq = series * series 174 175 # sum_dsq1 is the cumulative sum of dsq 176 sum_dsq1 = N.add.accumulate(dsq) 177 178 # sum_dsq1 is the reversed cumulative sum of dsq 179 sum_dsq2 = N.add.accumulate(dsq[::-1]) 180 181 # sumsq refers to SUMSQ in the published algorithm 182 sumsq = 2.*sum_dsq1[-1] 183 184 # this line refers to the instruction SUMSQ <-- SUMSQ - DSQ(m-1) - DSQ(N - m) of the published algorithm 185 # In this case, msd is an array because the instruction is computed for each m ranging from 0 to len(traj) - 1 186 # So, this single instruction is performing the loop in the published algorithm 187 Saabb = sumsq - N.concatenate(([0.], sum_dsq1[:-1])) - N.concatenate(([0.], sum_dsq2[:-1])) 188 189 # Saabb refers to SAA+BB/(N-m) in the published algorithm 190 # Sab refers to SAB(m)/(N-m) in the published algorithm 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
201 - def combine(self, atom, x):
202 """ 203 """ 204 N.add(self.MSD, self.weights[atom]*x, self.MSD)
205
206 - def finalize(self):
207 """Finalizes the calculations (e.g. averaging the total term, output files creations ...). 208 """ 209 210 print self.times 211 print self.MSD 212 # The NetCDF output file is opened for writing. 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 # Some dimensions are created. 218 outputFile.createDimension('NFRAMES', self.nFrames) 219 220 # Creation of the NetCDF output variables. 221 # The time. 222 TIMES = outputFile.createVariable('time', N.Float, ('NFRAMES',)) 223 TIMES[:] = self.times[:] 224 TIMES.units = 'ps' 225 226 # The MSD. 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 # Creates an ASCII version of the NetCDF output file. 236 convertNetCDFToASCII(inputFile = self.outputMSD,\ 237 outputFile = os.path.splitext(self.outputMSD)[0] + '.cdl',\ 238 variables = ['time', 'msd'])
239
240 - def atomicMSD(self, atom, series):
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 # ROOT MEAN-SQUARE DEVIATION ANALYSIS 256 #####################################################################################
257 -class RootMeanSquareDeviation(Analysis):
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 # The minimal set of input parameters names necessary to perform the analysis. 288 inputParametersNames = ('trajectory',\ 289 'timeinfo',\ 290 'referenceframe',\ 291 'subset',\ 292 'deuteration',\ 293 'weights',\ 294 'rmsd',\ 295 'pyroserver',\ 296 ) 297 298 shortName = 'RMSD' 299 300 # Indicate whether this analysis can be estimated or not. 301 canBeEstimated = True 302
303 - def __init__(self):
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
308 - def initialize(self):
309 """Initializes the analysis (e.g. parses and checks input parameters, set some variables ...). 310 """ 311 312 # The input parameters are parsed. 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 # RMSD = 1D Numeric array. 338 self.RMSD = N.zeros((self.nFrames), N.Float) 339 340 if self.pyroServer != 'monoprocessor': 341 # The attribute trajectory is removed because it can not be pickled by Pyro. 342 delattr(self, 'trajectory')
343 344 # def calc(self, frameIndex, trajname): 345 # """Calculates the contribution for one frame. 346 # 347 # @param frameIndex: the index of the frame in |self.frameIndexes| array. 348 # @type frameIndex: integer. 349 # 350 # @param trajname: the name of the trajectory file name. 351 # @type trajname: string 352 # """ 353 354 # if self.pyroServer[0] == 'no': 355 # t = self.trajectory 356 # else: 357 # # Load the whole trajectory set. 358 # t = Trajectory(None, trajname, 'r') 359 360 # frame = self.frameIndexes[frameIndex] 361 # msd = t.configuration[frame] - t.configuration[self.referenceFrame] 362 # msd = self.weights * msd * msd 363 # self.RMSD[frameIndex] = N.add.reduce(msd) 364 365 # return None 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 # Load the whole trajectory set. 389 t = Trajectory(None, trajname, 'r') 390 391 # series = 2D Numeric array. The positions of the selected atom |at| from the first step to the 392 # last step with the selected step increment. 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
407 - def combine(self, atom, x):
408 """ 409 """ 410 411 N.add(self.RMSD, self.weights[atom]*x, self.RMSD)
412
413 - def finalize(self):
414 """Finalizes the calculations (e.g. averaging the total term, output files creations ...) 415 """ 416 417 self.RMSD = N.sqrt(self.RMSD) 418 419 # The NetCDF output file is opened for writing. 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 # Some dimensions are created. 425 outputFile.createDimension('NFRAMES', self.nFrames) 426 427 # Creation of the NetCDF output variables. 428 # The time. 429 TIMES = outputFile.createVariable('time', N.Float, ('NFRAMES',)) 430 TIMES[:] = self.times[:] 431 TIMES.units = 'ps' 432 433 # The RMSD. 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 # Create an ASCII version of the NetCDF output file. 443 convertNetCDFToASCII(inputFile = self.outputRMSD,\ 444 outputFile = os.path.splitext(self.outputRMSD)[0] + '.cdl',\ 445 variables = ['time', 'rmsd'])
446 447 ##################################################################################### 448 # CARTESIAN VELOCITY AUTOCORRELATION FUNCTION ANALYSIS 449 #####################################################################################
450 -class CartesianVelocityAutoCorrelationFunction(Analysis):
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 # The minimal set of input parameters names necessary to perform the analysis. 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 # Indicate whether this analysis can be estimated or not. 503 canBeEstimated = True 504
505 - def __init__(self):
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
510 - def initialize(self):
511 """Initializes the analysis (e.g. parses and checks input parameters, set some variables ...). 512 """ 513 514 # The input parameters are parsed. 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 # The attribute trajectory is removed because it can not be pickled by Pyro. 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 # Load the whole trajectory set. 556 t = Trajectory(None, trajname, 'r') 557 558 if self.differentiation != 0: 559 560 if self.preLoadedTraj is None: 561 # series = 2D Numeric array. The positions of the selected atom |at| from the first step to the 562 # last step with the selected step increment. 563 series = t.readParticleTrajectory(atom, first = self.first, last = self.last, skip = self.skip).array 564 else: 565 series = self.preLoadedTraj[atom] 566 567 # the x, y and z axis 568 for axis in range(3): 569 # diff is a function of CalcFunctions.py module of nMoldyn package 570 # it generates the velocities out of positions 571 # diffScheme is the algorithm used to derive velocities from position and can be one 572 # of the keywords 1 (default), 2, 3, 4, 5 573 # dt is the differentiation time 574 series[:,axis] = differentiate(series[:,axis], self.differentiation, self.dt) 575 else: 576 if self.preLoadedTraj is None: 577 # series contains the velocities of the selected atom at from timeInfo[O] to timeInfo[1] with a time 578 # increment of timeInfo[2] 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 # if the trajectory is not projected, the calculations is done on the actual dependant variable 585 atomicVACF = correlation(series)/3. 586 587 else: 588 # if a projection vector is given, the trajectory is then projected onto the projection vector. 589 projectedTrajectory = N.dot(series, self.projection) 590 # if the trajectory has been projected, the calculations is performed on the projected trajectory 591 atomicVACF = correlation(projectedTrajectory) 592 593 if self.preLoadedTraj is not None: 594 del self.preLoadedTraj[atom] 595 596 return atomicVACF
597
598 - def combine(self, atom, x):
599 """ 600 """ 601 602 N.add(self.VACF, self.weights[atom]*x, self.VACF)
603
604 - def finalize(self):
605 """Finalizes the calculations (e.g. averaging the total term, output files creations ...) 606 """ 607 608 # It is of common use to normalize the VACF by its value at t=0. 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 # The NetCDF output file is opened for writing. 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 # Some dimensions are created. 621 outputFile.createDimension('NFRAMES', self.nFrames) 622 623 # Creation of the NetCDF output variables. 624 # The time. 625 TIMES = outputFile.createVariable('time', N.Float, ('NFRAMES',)) 626 TIMES[:] = self.times[:] 627 TIMES.units = 'ps' 628 629 # The VACF. 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 # Create an ASCII version of the NetCDF output file. 642 convertNetCDFToASCII(inputFile = self.outputVACF,\ 643 outputFile = os.path.splitext(self.outputVACF)[0] + '.cdl',\ 644 variables = ['time', 'vacf'])
645 646 ##################################################################################### 647 # CARTESIAN DENSITY OF STATES ANALYSIS 648 #####################################################################################
649 -class CartesianDensityOfStates(Analysis):
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 # The minimal set of input parameters names necessary to perform the analysis. 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 # Indicate whether this analysis can be estimated or not. 702 canBeEstimated = True 703
704 - def __init__(self):
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
709 - def initialize(self):
710 """Initializes the analysis (e.g. parses and checks input parameters, set some variables ...). 711 """ 712 713 # The input parameters are parsed. 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 # The attribute trajectory is removed because it can not be pickled by Pyro. 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 # Load the whole trajectory set. 755 t = Trajectory(None, trajname, 'r') 756 757 if self.differentiation != 0: 758 759 if self.preLoadedTraj is None: 760 # series = 2D Numeric array. The positions of the selected atom |at| from the first step to the 761 # last step with the selected step increment. 762 series = t.readParticleTrajectory(atom, first = self.first, last = self.last, skip = self.skip).array 763 else: 764 series = self.preLoadedTraj[atom] 765 766 # the x, y and z axis 767 for axis in range(3): 768 # diff is a function of CalcFunctions.py module of nMoldyn package 769 # it generates the velocities out of positions 770 # diffScheme is the algorithm used to derive velocities from position and can be one 771 # of the keywords 1 (default), 2, 3, 4, 5 772 # dt is the differentiation time 773 series[:,axis] = differentiate(series[:,axis], self.differentiation, self.dt) 774 else: 775 if self.preLoadedTraj is None: 776 # series contains the velocities of the selected atom at from timeInfo[O] to timeInfo[1] with a time 777 # increment of timeInfo[2] 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 # if the trajectory is not projected, the calculations is done on the actual dependant variable 784 # The 1/3 factor is used to normalize the VACF 785 aVACF = correlation(series)/3. 786 else: 787 # if a projection vector is given, the trajectory is then projected onto the projection vector. 788 projectedTrajectory = N.dot(series, self.projection) 789 # if the trajectory has been projected, the calculations is performed on the projected trajectory 790 aVACF = correlation(projectedTrajectory) 791 792 # the dos is the fft of the product between the gaussian kernel and Sab 793 # The DOS for atom 'at' is added to the dos array 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
801 - def combine(self, atom, x):
802 """ 803 """ 804 805 N.add(self.DOS, self.weights[atom]*x, self.DOS)
806
807 - def finalize(self):
808 """Finalizes the calculations (e.g. averaging the total term, output files creations ...) 809 """ 810 811 # 'freqencies' = 1D Numeric array. Frequencies at which the DOS was computed 812 frequencies = N.arange(self.nFrames)/(2.0*self.nFrames*self.dt) 813 814 # The final form for the DOS. 815 self.DOS = 0.5*self.dt*self.DOS 816 817 # The NetCDF output file is opened for writing. 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 # Some dimensions are created. 823 outputFile.createDimension('FREQSTEPS', self.nFrames) 824 825 # Creation of the NetCDF output variables. 826 # The frequencies. 827 FREQUENCIES = outputFile.createVariable('frequency', N.Float, ('FREQSTEPS',)) 828 FREQUENCIES[:] = frequencies[:] 829 FREQUENCIES.units = 'THz' 830 831 # The DOS. 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 # Create an ASCII version of the NetCDF output file. 841 convertNetCDFToASCII(inputFile = self.outputDOS,\ 842 outputFile = os.path.splitext(self.outputDOS)[0] + '.cdl',\ 843 variables = ['frequency', 'dos'])
844 845 ##################################################################################### 846 # AUTO-REGRESSIVE ANALYSIS ALYSIS 847 #####################################################################################
848 -class AutoRegressiveAnalysis(Analysis):
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 # The minimal set of input parameters names necessary to perform the analysis. 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 # Indicate whether this analysis can be estimated or not. 900 canBeEstimated = True 901
902 - def __init__(self):
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
907 - def initialize(self):
908 """Initializes the analysis (e.g. parses and checks input parameters, set some variables ...). 909 """ 910 911 # The input parameters are parsed. 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 # The attribute trajectory is removed because it can not be pickled by Pyro. 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 # Load the whole trajectory set. 954 t = Trajectory(None, trajname, 'r') 955 956 if self.differentiation != 0: 957 958 if self.preLoadedTraj is None: 959 # series = 2D Numeric array. The positions of the selected atom |at| from the first step to the 960 # last step with the selected step increment. 961 series = t.readParticleTrajectory(atom, first = self.first, last = self.last, skip = self.skip).array 962 else: 963 series = self.preLoadedTraj[atom] 964 965 # the x, y and z axis 966 for axis in range(3): 967 # diff is a function of CalcFunctions.py module of nMoldyn package 968 # it generates the velocities out of positions 969 # diffScheme is the algorithm used to derive velocities from position and can be one 970 # of the keywords 1 (default), 2, 3, 4, 5 971 # dt is the differentiation time 972 series[:,axis] = differentiate(series[:,axis], self.differentiation, self.dt) 973 else: 974 if self.preLoadedTraj is None: 975 # series contains the velocities of the selected atom at from timeInfo[O] to timeInfo[1] with a time 976 # increment of timeInfo[2] 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
986 - def combine(self, atom, x):
987 """ 988 """ 989 990 # The x, y and z axis 991 for axis in range(3): 992 self.ARModel.add(AutoRegressiveModel(self.arModelOrder, x[:, axis], self.dt), self.weights[atom])
993
994 - def finalize(self):
995 """Finalizes the calculations (e.g. averaging the total term, output files creations ...) 996 """ 997 998 # This is the former fonction "evaluate_model". This fonction was used only twice 999 # in the AutoRegressiveAnalysis and the AutoRegressiveAnalysisXYZ functions. As these 1000 # functions has been fused in the current implementation no need for a function. 1001 1002 # The VACF within the ARA framework. 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 # The DOS within the ARA framework. 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 # The MSD within the ARA framework. 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 # The Memory function 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 # The NetCDF output file is opened for writing. 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 # Some dimensions are created. 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 # Creation of the NetCDF output variables. 1069 1070 # VACF related variables. 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 # DOS related variables. 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 # MSD related variables. 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 # Memory function related variables. 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 # ARA parameters related variables. 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 # Create an ASCII version of the NetCDF output file. 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 # PASSBANDFILTEREDTRAJECTORY ANALYSIS 1135 #####################################################################################
1136 -class PassBandFilteredTrajectory(Analysis):
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 # The minimal set of input parameters names necessary to perform the analysis. 1164 inputParametersNames = ('trajectory',\ 1165 'timeinfo',\ 1166 'filter',\ 1167 'subset',\ 1168 'pbft',\ 1169 'pyroserver',\ 1170 ) 1171 1172 shortName = 'PBFT' 1173 1174 # Indicate whether this analysis can be estimated or not. 1175 canBeEstimated = True 1176
1177 - def __init__(self):
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
1182 - def initialize(self):
1183 """Initializes the analysis (e.g. parses and checks input parameters, set some variables ...). 1184 """ 1185 1186 # The input parameters are parsed. 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 # The frequencies. 1202 frequencies = N.arange(self.nFrames)/(2.0*self.nFrames*self.dt) 1203 # Computation of the index corresponding to the lowest frequency of the pass-band filter. 1204 indexMin = N.searchsorted(frequencies, self.freqMin) 1205 # Computation of the index corresponding to the highest frequency of the pass-band filter. 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 # This dictionnary works a little bit like a Particle Tensor object. Its keys are the index of the atoms of 1216 # the selection. Its value are the filtered atomic trajectories. 1217 self.filteredTrajectory = {} 1218 1219 if self.pyroServer != 'monoprocessor': 1220 # The attribute trajectory is removed because it can not be pickled by Pyro. 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 # Load the whole trajectory set. 1239 t = Trajectory(None, trajname, 'r') 1240 1241 # series = 2D Numeric array. The positions of the selected atom |at| from the first step to the 1242 # last step with the selected step increment. 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
1253 - def combine(self, atom, x):
1254 """ 1255 """ 1256 1257 unfiltered = N.zeros(2*self.nFrames - 2, 'd') 1258 atomicPBFT = N.zeros((self.nFrames,3), 'd') 1259 1260 # Loop over the x, y and z coordinates. 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
1269 - def finalize(self):
1270 """Finalizes the calculations (e.g. averaging the total term, output files creations ...) 1271 """ 1272 1273 # traj_new is the filtered trajectory 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 # Load the whole trajectory set. 1282 t = Trajectory(None, self.trajectoryFilename, 'r') 1283 1284 # Create the snapshot generator 1285 snapshot = SnapshotGenerator(self.universe, actions = [TrajectoryOutput(outputFile, ["all"], 0, None, 1)]) 1286 1287 # Loop over the output frames. 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 # GYRATION RADIUS ANALYSIS 1304 #####################################################################################
1305 -class RadiusOfGyration(Analysis):
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 # The minimal set of input parameters names necessary to perform the analysis. 1332 inputParametersNames = ('trajectory',\ 1333 'timeinfo',\ 1334 'subset',\ 1335 'rog',\ 1336 'pyroserver',\ 1337 ) 1338 1339 shortName = 'ROG' 1340 1341 # Indicate whether this analysis can be estimated or not. 1342 canBeEstimated = True 1343
1344 - def __init__(self):
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
1349 - def initialize(self):
1350 """Initializes the analysis (e.g. parses and checks input parameters, set some variables ...). 1351 """ 1352 1353 # The input parameters are parsed. 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 # The center of gravity. 1369 self.cog = N.zeros((self.nFrames,3), N.Float) 1370 1371 # The sum of each position. 1372 self.sumRi = N.zeros((self.nFrames,3), N.Float) 1373 1374 # ROG = 1D Numeric array. The array that stores the Gyration Radius over time. 1375 self.ROG = N.zeros((self.nFrames), N.Float) 1376 1377 if self.pyroServer != 'monoprocessor': 1378 # The attribute trajectory is removed because it can not be pickled by Pyro. 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 # Load the whole trajectory set. 1397 t = Trajectory(None, trajname, 'r') 1398 1399 # series = 2D Numeric array. The positions of the selected atom |at| from the first step to the 1400 # last step with the selected step increment. 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
1411 - def combine(self, atom, x):
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
1419 - def finalize(self):
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 # The NetCDF output file is opened for writing. 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 # Some dimensions are created. 1436 outputFile.createDimension('NFRAMES', self.nFrames) 1437 1438 # Creation of the NetCDF output variables. 1439 # The time. 1440 TIMES = outputFile.createVariable('time', N.Float, ('NFRAMES',)) 1441 TIMES[:] = self.times[:] 1442 TIMES.units = 'ps' 1443 1444 # The Gyration radius. 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 # Create an ASCII version of the NetCDF output file. 1454 convertNetCDFToASCII(inputFile = self.outputROG,\ 1455 outputFile = os.path.splitext(self.outputROG)[0] + '.cdl',\ 1456 variables = ['time', 'rog'])
1457 1458 ##################################################################################### 1459 # GlobalMotionTrajectoryFilter ANALYSIS 1460 #####################################################################################
1461 -class GlobalMotionFilteredTrajectory(Analysis):
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 # The minimal set of input parameters names necessary to perform the analysis. 1487 inputParametersNames = ('trajectory',\ 1488 'timeinfo',\ 1489 'subset',\ 1490 'gmft',\ 1491 'pyroserver',\ 1492 ) 1493 1494 shortName = 'GMFT' 1495 1496 # Indicate whether this analysis can be estimated or not. 1497 canBeEstimated = True 1498
1499 - def __init__(self):
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
1504 - def initialize(self):
1505 """Initializes the analysis (e.g. parses and checks input parameters, set some variables ...). 1506 """ 1507 1508 # The input parameters are parsed. 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 # This dictionnary works a little bit like a Particle Tensor object. Its keys are the index of the atoms of 1524 # the selection. Its value are the filtered atomic trajectories. 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 # The attribute trajectory is removed because it can not be pickled by Pyro. 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 # Load the whole trajectory set. 1554 t = Trajectory(None, trajname, 'r') 1555 1556 conf = t.configuration[frame] 1557 self.universe.setConfiguration(conf) 1558 1559 else: 1560 # The configuration is updated to the current trajectory step. 1561 self.universe.setConfiguration(self.preLoadedTraj[frameIndex]) 1562 1563 # For the fist frame, performs a principal axes transformation. 1564 if frameIndex == 0: 1565 tr = self.subset.normalizingTransformation() 1566 1567 # For the nest frames, find and apply the transfo that wll minimize the RMS with the configuration 1568 # resulting from the principal axes transformation. 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 # For the fist frame, performs a principal axes transformation. 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
1591 - def finalize(self):
1592 """Finalizes the calculations (e.g. averaging the total term, output files creations ...). 1593 """ 1594 1595 # traj_new is the filtered trajectory 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 # Load the whole trajectory set. 1603 t = Trajectory(None, self.trajectoryFilename, 'r') 1604 1605 # Create the snapshot generator 1606 snapshot = SnapshotGenerator(self.universe, actions = [TrajectoryOutput(outputFile, ["all"], 0, None, 1)]) 1607 1608 # Loop over the output frames. 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 # CENTER OF MASS TRAJECTORY ANALYSIS 1625 #####################################################################################
1626 -class CenterOfMassTrajectory(Analysis):
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 # The minimal set of input parameters names necessary to perform the analysis. 1653 inputParametersNames = ('trajectory',\ 1654 'timeinfo',\ 1655 'group',\ 1656 'comt',\ 1657 'pyroserver',\ 1658 ) 1659 1660 shortName = 'COMT' 1661 1662 # Indicate whether this analysis can be estimated or not. 1663 canBeEstimated = True 1664
1665 - def __init__(self):
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
1670 - def initialize(self):
1671 """Initializes the analysis (e.g. parses and checks input parameters, set some variables ...). 1672 """ 1673 1674 # The input parameters are parsed. 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 # The universe that will contain the com trajectory is build from the input trajectory. 1690 self.superAtomUniverse = self.universe.__copy__() 1691 # All the atoms of the universe are deleted. 1692 self.superAtomUniverse.removeObject(self.superAtomUniverse.objectList()[:]) 1693 1694 # This dictionnary works a little bit like a Particle Tensor object. Its keys are the index of the atoms of 1695 # the selection. Its value are the filtered atomic trajectories. 1696 self.comTrajectory = {} 1697 1698 # This dictionnary is a lookup table between the superAtom representing the com and its associated group. 1699 self.comToGroup = {} 1700 # Each COM trajectory is represented by the trajectory of a 'virtual' hydrogen atom. The name of the 1701 # hydrogen is COMn where n is the COM number. 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 # The attribute trajectory is removed because it can not be pickled by Pyro. 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 # Load the whole trajectory set. 1731 t = Trajectory(None, trajname, 'r') 1732 1733 conf = t.configuration[frame] 1734 self.universe.setConfiguration(conf) 1735 1736 else: 1737 # The configuration is updated to the current trajectory step. 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
1748 - def combine(frameIndex, x):
1749 """ 1750 """ 1751 1752 pass
1753
1754 - def finalize(self):
1755 """Finalizes the calculations (e.g. averaging the total term, output files creations ...). 1756 """ 1757 1758 # traj_new is the filtered trajectory 1759 outputFile = Trajectory(self.superAtomUniverse, self.outputCOMT, "w") 1760 outputFile.jobinfo = self.information + '\nOutput file written on: %s\n\n' % asctime() 1761 1762 # Each time |snapshot| is called, the universe contents i flushed into the output file. 1763 snapshot = SnapshotGenerator(self.superAtomUniverse,\ 1764 actions = [TrajectoryOutput(outputFile, ["all"], 0, None, 1)]) 1765 1766 # Loop over the output frames. 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 # The output COM trajectory is closed. 1775 outputFile.close() 1776 1777 self.toPlot = None
1778 1779 ##################################################################################### 1780 # QUASI-HARMONIC ANALYSIS ANALYSIS 1781 #####################################################################################
1782 -class QuasiHarmonicAnalysis(Analysis):
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 # The minimal set of input parameters names necessary to perform the analysis. 1815 inputParametersNames = ('trajectory',\ 1816 'timeinfo',\ 1817 'temperature',\ 1818 'subset',\ 1819 'qha') 1820 1821 shortName = 'QHA' 1822 1823 # Indicate whether this analysis can be estimated or not. 1824 canBeEstimated = False 1825
1826 - def __init__(self):
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
1831 - def initialize(self):
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
1842 - def internalRun(self):
1843 """Runs the analysis.""" 1844 1845 # mask for the selected atoms. 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 # The initial structure configuration. 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 # Calculate the fluctuation matrix. 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 # Calculate the quasiharmonic modes 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 # Due to numerical imprecisions, the result can have imaginary parts. 1884 # In that case, throw the imaginary parts away. 1885 # Conversion from uma*nm2 to kg*m2 (SI) 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 # Sort eigen vectors by decreasing fluctuation amplitude. 1893 indices = N.argsort(omega)[::-1] 1894 1895 omega = N.take(omega, indices) 1896 dx = N.take(dx, indices) 1897 1898 # Eq 66 of the reference paper. 1899 mdr = N.take(mdr, indices, axis = 1) 1900 at = N.dot(mdr,N.transpose(dx)) 1901 1902 # The NetCDF output file is opened. 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 # The universe is emptied from its objects keeping just its topology. 1908 self.universe.removeObject(self.universe.objectList()[:]) 1909 # The atoms of the subset are copied 1910 atoms = copy.deepcopy(self.subset.atomList()) 1911 # And their parent attribute removed to allow their transfer in the empty universe. 1912 for a in atoms: 1913 a.parent = None 1914 ac = AtomCluster(atoms,name='QHACluster') 1915 self.universe.addObject(ac) 1916 1917 # Some dimensions are created. 1918 # NEIVALS = the number eigen values 1919 outputFile.createDimension('NEIGENVALS', len(omega)) 1920 1921 # UDESCR = the universe description length 1922 outputFile.createDimension('UDESCR', len(self.universe.description())) 1923 1924 # NATOMS = the number of atoms of the universe 1925 outputFile.createDimension('NATOMS', self.nAtoms) 1926 1927 # NFRAMES = the number of frames. 1928 outputFile.createDimension('NFRAMES', self.nFrames) 1929 1930 # NCOORDS = the number of coordinates (always = 3). 1931 outputFile.createDimension('NCOORDS', 3) 1932 1933 # 3N. 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 # Creation of the NetCDF output variables. 1940 # EIVALS = the eigen values. 1941 OMEGA = outputFile.createVariable('omega', N.Float, ('NEIGENVALS',)) 1942 OMEGA[:] = omega 1943 1944 # EIVECS = array of eigen vectors. 1945 DX = outputFile.createVariable('dx', N.Float, ('NEIGENVALS','NEIGENVALS')) 1946 DX[:,:] = dx 1947 1948 # The time. 1949 TIMES = outputFile.createVariable('time', N.Float, ('NFRAMES',)) 1950 TIMES[:] = self.times[:] 1951 TIMES.units = 'ps' 1952 1953 # MODE = the mode number. 1954 MODE = outputFile.createVariable('mode', N.Float, ('3N',)) 1955 MODE[:] = 1 + N.arange(3*self.nAtoms) 1956 1957 # LCI = local character indicator. See eq 56. 1958 LCI = outputFile.createVariable('local_character_indicator', N.Float, ('3N',)) 1959 LCI[:] = (dx**4).sum(0) 1960 1961 # GCI = global character indicator. See eq 57. 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 # Projection of MD traj onto normal modes. See eq 66.. 1966 AT = outputFile.createVariable('at', N.Float, ('NFRAMES','3N')) 1967 AT[:,:] = at[:,:] 1968 1969 # DESCRIPTION = the universe description. 1970 DESCRIPTION = outputFile.createVariable('description', N.Character, ('UDESCR',)) 1971 DESCRIPTION[:] = self.universe.description() 1972 1973 # AVGSTRUCT = the average structure. 1974 AVGSTRUCT = outputFile.createVariable('avgstruct', N.Float, ('NATOMS','NCOORDS')) 1975 AVGSTRUCT[:,:] = N.compress(mask.array,averageStructure.array,0) 1976 1977 # If the universe is periodic, create an extra variable storing the box size. 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 # ANGULAR CORRELATION ANALYSIS 1990 #####################################################################################
1991 -class AngularCorrelation(Analysis):
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 # The minimal set of input parameters names necessary to perform the analysis. 2021 inputParametersNames = ('trajectory',\ 2022 'timeinfo',\ 2023 'triplet',\ 2024 'atomorder',\ 2025 'ac',\ 2026 'pyroserver',\ 2027 ) 2028 2029 shortName = 'AC' 2030 2031 # Indicate whether this analysis can be estimated or not. 2032 canBeEstimated = True 2033
2034 - def __init__(self):
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
2039 - def initialize(self):
2040 """Initializes the analysis (e.g. parses and checks input parameters, set some variables ...). 2041 """ 2042 2043 # The input parameters are parsed. 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 # The self.nGroups is needed when dealing with the template defined for analysis whose mainloop is over groups 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 # The attribute trajectory is removed because it can not be pickled by Pyro. 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 # Load the whole trajectory set. 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 # The autocorrelation for v1, v2, v3 are computed. As they are unit vectors AC1(0) = AC2(0) = AC3(0)= 1 2149 self.angCorr1[tripletIndex,:] = correlation(x[0]) 2150 self.angCorr2[tripletIndex,:] = correlation(x[1]) 2151 self.angCorr3[tripletIndex,:] = correlation(x[2])
2152 2153
2154 - def finalize(self):
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 # Creates an ASCII version of the NetCDF output file. 2197 convertNetCDFToASCII(inputFile = self.outputAC,\ 2198 outputFile = os.path.splitext(self.outputAC)[0] + '.cdl',\ 2199 variables = asciiVar)
2200 2201 ##################################################################################### 2202 # RIGID BODY TRAJECTORY ANALYSIS 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 # The minimal set of input parameters names necessary to perform the analysis. 2234 inputParametersNames = ('trajectory',\ 2235 'timeinfo',\ 2236 'referenceframe',\ 2237 'removetranslation',\ 2238 'stepwiserbt',\ 2239 'group',\ 2240 'rbt',\ 2241 'pyroserver',\ 2242 ) 2243 2244 shortName = 'RBT' 2245 2246 # Indicate whether this analysis can be estimated or not. 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 # The input parameters are parsed. 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 # Check that there is no overlap between the groups. 2278 overlap = set() 2279 for g in self.group: 2280 overlap = overlap.union(set(g.atomList())) 2281 2282 # It would be no sense to have some atoms belonging to several groups. 2283 if len(overlap) != sum(self.groupsSize): 2284 raise Error('There must not be any common atom between the groups.') 2285 2286 # If a fixed reference has been set. We can already set the reference configuration here. 2287 self.refConfig = self.trajectory.configuration[self.referenceFrame] 2288 2289 # Those matrix will store the quaternions and the CMS coming from the RBT trajectory. 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 # This dictionnary works a little bit like a Particle Tensor object. Its keys are the index of the atoms of 2295 # the selection. Its value are the rigid-body atomic trajectories. 2296 self.rbtTraj = {} 2297 2298 if self.pyroServer != 'monoprocessor': 2299 # The attribute trajectory is removed because it can not be pickled by Pyro. 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 # Load the whole trajectory set. 2317 t = Trajectory(None, trajname, 'r') 2318 2319 group = self.group[groupIndex] 2320 2321 # Those matrix will store the quaternions and the CMS coming from the RBT trajectory. 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 # Case of a moving reference. 2327 if self.stepwiseRBT: 2328 2329 # The reference configuration is always the one of the previous frame excepted for the first frame 2330 # where it is set by definition to the first frame (could we think about a cyclic alternative way ?). 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 # The RBT is created just for the current step. 2338 rbt = t.readRigidBodyTrajectory(group,\ 2339 first = self.frameIndexes[frameIndex],\ 2340 last = self.frameIndexes[frameIndex] + 1,\ 2341 skip = 1,\ 2342 reference = movingRefConfig) 2343 2344 # The corresponding quaternions and cms are stored in their corresponding matrix. 2345 quaternions[frameIndex,:] = copy.copy(rbt.quaternions) 2346 CMS[frameIndex,:] = copy.copy(rbt.cms) 2347 fit[frameIndex] = copy.copy(rbt.fit) 2348 2349 # The simplest case, the reference frame is fixed. 2350 # A unique RBT is performed from first to last skipping skip steps and using refConfig as the reference. 2351 else: 2352 # The RBT is created. 2353 rbt = t.readRigidBodyTrajectory(group,\ 2354 first = self.first,\ 2355 last = self.last,\ 2356 skip = self.skip,\ 2357 reference = self.refConfig) 2358 2359 # The corresponding quaternions and cms are stored in their corresponding matrix. 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 # I can not use the centers of mass defined by rbt.cms because the reference frame 2377 # selected can be out of the selected frames for the Rigid Body Trajectory. 2378 centerOfMass = group.centerOfMass(self.refConfig) 2379 2380 atomicRBT = N.zeros((self.nFrames,3), typecode = N.Float) 2381 2382 # Loop over the atoms of the group to set the RBT trajectory. 2383 for atom in group: 2384 # The coordinates of the atoms are centered around the center of mass of the group. 2385 xyz = self.refConfig[atom] - centerOfMass 2386 2387 # Loop over the selected frames. 2388 for frameIndex in range(self.nFrames): 2389 # The rotation matrix corresponding to the selected frame in the RBT. 2390 transfo = Quaternion(self.rbtQuaternions[frameIndex,groupIndex,:]).asRotation() 2391 2392 if self.removeTranslation: 2393 # The transformation matrix corresponding to the selected frame in the RBT. 2394 transfo = Translation(centerOfMass)*transfo 2395 2396 # Compose with the CMS translation if the removeTranslation flag is set off. 2397 else: 2398 # The transformation matrix corresponding to the selected frame in the RBT. 2399 transfo = Translation(Vector(self.rbtCMS[frameIndex,groupIndex,:]))*transfo 2400 2401 # The RBT is performed on the CMS centered coordinates of atom at. 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 # Create trajectory 2413 outputFile = Trajectory(selection, self.outputRBT, 'w') 2414 2415 if self.pyroServer == 'monoprocessor': 2416 t = self.trajectory 2417 2418 else: 2419 # Load the whole trajectory set. 2420 t = Trajectory(None, self.trajectoryFilename, 'r') 2421 2422 # Create the snapshot generator 2423 snapshot = SnapshotGenerator(self.universe, actions = [TrajectoryOutput(outputFile, ["all"], 0, None, 1)]) 2424 2425 # The output is written 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 # Dictionnary whose keys are of the form Gi where i is the group number 2439 # and the entries are the list of the index of the atoms building the group. 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 # The NetCDF variable that stores the quaternions. 2447 quaternions = outputFile.createVariable('quaternion', N.Float32, ('step_number','number_of_groups','quaternion_length')) 2448 2449 # The NetCDF variable that stores the centers of mass. 2450 cms = outputFile.createVariable('cms', N.Float32, ('step_number','number_of_groups','xyz')) 2451 2452 # The NetCDF variable that stores the rigid-body fit. 2453 fit = outputFile.createVariable('fit', N.Float32, ('step_number','number_of_groups')) 2454 2455 # Loop over the groups. 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 # REORIENTATIONAL CORRELATION FUNCTION ANALYSIS 2467 #####################################################################################
2468 -class ReorientationalCorrelationFunction(Analysis):
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 # The minimal set of input parameters names necessary to perform the analysis. 2501 inputParametersNames = ('trajectory',\ 2502 'timeinfo',\ 2503 'referenceframe',\ 2504 'stepwiserbt',\ 2505 'wignerindexes',\ 2506 'group',\ 2507 'rcf',\ 2508 'pyroserver',\ 2509 ) 2510 2511 shortName = 'RCF' 2512 2513 # Indicate whether this analysis can be estimated or not. 2514 canBeEstimated = True 2515
2516 - def __init__(self):
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
2521 - def initialize(self):
2522 """Initializes the analysis (e.g. parses and checks input parameters, set some variables ...). 2523 """ 2524 2525 # The input parameters are parsed. 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 # If a fixed reference has been set. We can already set the reference configuration here. 2546 self.refConfig = self.trajectory.configuration[self.referenceFrame] 2547 2548 # The array that stores the reorientational function. 2549 self.RCF = N.zeros((self.nFrames), typecode = N.Float) 2550 2551 if self.pyroServer != 'monoprocessor': 2552 # The attribute trajectory is removed because it can not be pickled by Pyro. 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 # Load the whole trajectory set. 2570 t = Trajectory(None, trajname, 'r') 2571 2572 group = self.group[groupIndex] 2573 2574 j, m, n = self.wignerIndexes 2575 2576 # This matrix will store the quaternions coming from the RBT trajectory. 2577 rbtQuaternions = N.zeros((self.nFrames,4), typecode = N.Float) 2578 2579 # Case of a moving reference. 2580 if self.stepwiseRBT: 2581 # The reference configuration is always the one of the previous frame excepted for the first frame 2582 # where it is set by definition to the first frame (could we think about a cyclic alternative way ?). 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 # The RBT is created just for the current step. 2590 rbt = t.readRigidBodyTrajectory(group,\ 2591 first = self.frameIndexes[frameIndex],\ 2592 last = self.frameIndexes[frameIndex] + 1,\ 2593 skip = 1,\ 2594 reference = movingRefConfig) 2595 # The corresponding quaternions and cms are stored in their corresponding matrix. 2596 rbtQuaternions[frameIndex,:] = copy.copy(rbt.quaternions) 2597 rbtCMS[frameIndex,:] = copy.copy(rbt.cms) 2598 2599 # The simplest case, the reference frame is fixed. 2600 # A unique RBT is performed from first to last skipping skip steps and using refConfig as the reference. 2601 else: 2602 # The RBT is created. 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 # c1 is the scaling factor converting Wigner function into spherical harmonics. It depends only on j. 2612 c1 = N.sqrt(((2.0*j+1)/(4.0*N.pi))) 2613 2614 quat2 = N.zeros(rbtQuaternions.shape, typecode = N.Complex) 2615 # quat2[:,0] refers to the (q0+iq3) of equation 3.55 2616 quat2[:,0] = rbtQuaternions[:,0] + 1j*rbtQuaternions[:,3] 2617 # quat2[:,2] refers to the (q2+iq1) of equation 3.55 2618 quat2[:,2] = rbtQuaternions[:,2] + 1j*rbtQuaternions[:, 1] 2619 # quat2[:,1] refers to the (q0-iq3) of equation 3.55 2620 quat2[:,1] = N.conjugate(quat2[:,0]) 2621 # quat2[:,3] refers to the (q2-iq1) of equation 3.55 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
2644 - def finalize(self):
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 # The NetCDF output file is opened for writing. 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 # Dictionnary whose keys are of the form Gi where i is the group number 2656 # and the entries are the list of the index of the atoms building the group. 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 # Some dimensions are created. 2661 outputFile.createDimension('TIMES', self.nFrames) 2662 2663 # Creation of the NetCDF output variables. 2664 # The time. 2665 TIMES = outputFile.createVariable('time', N.Float, ('TIMES',)) 2666 TIMES[:] = self.times 2667 TIMES.units = 'ps' 2668 2669 # The RCF. 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 # Create an ASCII version of the NetCDF output file. 2679 convertNetCDFToASCII(inputFile = self.outputRCF,\ 2680 outputFile = os.path.splitext(self.outputRCF)[0] + '.cdl',\ 2681 variables = ['time', 'rcf'])
2682
2683 -class AngularVelocity:
2684 """An intermediate class used by |AngularVelocityAutoCorrelationFunction| and |AngularDensityOfStates| classes.""" 2685
2686 - def __init__(self):
2687 2688 pass
2689
2690 - def qMatrix(self, data):
2691 2692 res = N.zeros((len(data),4,4), N.Float) 2693 # qs 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 # qx 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 # qy 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 # qz 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
2716 - def getAngularVelocity(self, t, g):
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 # If a fixed reference has been set. We can already set the reference configuration here. 2722 refConfig = t.configuration[self.referenceFrame] 2723 2724 # Case of a moving reference. 2725 if self.stepwiseRBT: 2726 # The reference configuration is always the one of the previous frame excepted for the first frame 2727 # where it is set by definition to the first frame (could we think about a cyclic alternative way ?). 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 # The RBT is created just for the current step. 2735 rbt = t.readRigidBodyTrajectory(g,\ 2736 first = self.frameIndexes[i],\ 2737 last = self.frameIndexes[i] + 1,\ 2738 skip = 1,\ 2739 reference = movingRefConfig) 2740 # The corresponding quaternions and cms are stored in their corresponding matrix. 2741 rbtQuaternions[i,:] = copy.copy(rbt.quaternions) 2742 2743 # The simplest case, the reference frame is fixed. 2744 # A unique RBT is performed from first to last skipping skip steps and using refConfig as the reference. 2745 else: 2746 # The RBT is created. 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 # The quaternions derivatives. 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 # 1: refers to the vectorial part of the quaternion matrix. 2764 return angvel[:,1:]
2765 2766 ##################################################################################### 2767 # ANGULAR VELOCITY AUTOCORRELATION FUNCTION ANALYSIS 2768 #####################################################################################
2769 -class AngularVelocityAutoCorrelationFunction(Analysis,AngularVelocity):
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 # The minimal set of input parameters names necessary to perform the analysis. 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 # Indicate whether this analysis can be estimated or not. 2820 canBeEstimated = True 2821
2822 - def __init__(self):
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
2827 - def initialize(self):
2828 """Initializes the analysis (e.g. parses and checks input parameters, set some variables ...). 2829 """ 2830 2831 # The input parameters are parsed. 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 # Check that there is no overlap between the groups. 2851 overlap = set() 2852 for g in self.group: 2853 overlap = overlap.union(set(g.atomList())) 2854 2855 # It would be no sense to have some atoms belonging to several groups. 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 # The attribute trajectory is removed because it can not be pickled by Pyro. 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 # Load the whole trajectory set. 2880 t = Trajectory(None, trajname, 'r') 2881 2882 group = self.group[groupIndex] 2883 2884 AngularVelocity.__init__(self) 2885 2886 # angvel = 1D Numeric array. Angular velocity for the 'group'. 2887 angvel = self.getAngularVelocity(t, group) 2888 2889 return angvel
2890
2891 - def combine(self, groupIndex, x):
2892 """ 2893 """ 2894 2895 # the AVACF for group g 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
2904 - def finalize(self):
2905 """Finalizes the calculations (e.g. averaging the total term, output files creations ...). 2906 """ 2907 2908 # The NetCDF output file is opened for writing. 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 # Dictionnary whose keys are of the form Gi where i is the group number 2914 # and the entries are the list of the index of the atoms building the group. 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 # Some dimensions are created. 2919 outputFile.createDimension('NFRAMES', self.nFrames) 2920 2921 # Creation of the NetCDF output variables. 2922 # The time. 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 # Create an ASCII version of the NetCDF output file. 2948 convertNetCDFToASCII(inputFile = self.outputAVACF,\ 2949 outputFile = os.path.splitext(self.outputAVACF)[0] + '.cdl',\ 2950 variables = ['time', 'avacf-total'])
2951 2952 ##################################################################################### 2953 # ANGULAR DENSITY OF STATES ANALYSIS 2954 #####################################################################################
2955 -class AngularDensityOfStates(Analysis,AngularVelocity):
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 # The minimal set of input parameters names necessary to perform the analysis. 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 # Indicate whether this analysis can be estimated or not. 3009 canBeEstimated = True 3010
3011 - def __init__(self):
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
3016 - def initialize(self):
3017 """Initializes the analysis (e.g. parses and checks input parameters, set some variables ...). 3018 """ 3019 3020 # The input parameters are parsed. 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 # Check that there is no overlap between the groups. 3040 overlap = set() 3041 for g in self.group: 3042 overlap = overlap.union(set(g.atomList())) 3043 3044 # It would be no sense to have some atoms belonging to several groups. 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 # The attribute trajectory is removed because it can not be pickled by Pyro. 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 # Load the whole trajectory set. 3070 t = Trajectory(None, trajname, 'r') 3071 3072 group = self.group[groupIndex] 3073 3074 AngularVelocity.__init__(self) 3075 3076 # angvel = 1D Numeric array. Angular velocity for the 'group'. 3077 angvel = self.getAngularVelocity(t, group) 3078 3079 return angvel
3080
3081 - def combine(self, groupIndex, x):
3082 """ 3083 """ 3084 # the AVACF for group g 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
3093 - def finalize(self):
3094 """Finalizes the calculations (e.g. averaging the total term, output files creations ...). 3095 """ 3096 3097 # 'freqencies' = 1D Numeric array. Frequencies at which the DOS was computed 3098 frequencies = N.arange(self.nFrames)/(2.0*self.nFrames*self.dt) 3099 3100 # The NetCDF output file is opened for writing. 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 # Dictionnary whose keys are of the form Gi where i is the group number 3106 # and the entries are the list of the index of the atoms building the group. 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 # Some dimensions are created. 3111 outputFile.createDimension('NFRAMES', self.nFrames) 3112 3113 # Creation of the NetCDF output variables. 3114 # The time. 3115 TIMES = outputFile.createVariable('time', N.Float, ('NFRAMES',)) 3116 TIMES[:] = self.times[:] 3117 TIMES.units = 'ps' 3118 3119 # Creation of the NetCDF output variables. 3120 # The frequencies. 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 # Create an ASCII version of the NetCDF output file. 3154 convertNetCDFToASCII(inputFile = self.outputADOS,\ 3155 outputFile = os.path.splitext(self.outputADOS)[0] + '.cdl',\ 3156 variables = ['frequency', 'ados-total'])
3157