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