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

# Source Code for Module nMOLDYN.Analysis.Scattering

       1  """Collections of classes for the determination of scattering-related properties.
       2
       3  Classes:
       4      * DynamicCoherentStructureFactor           : sets up a Dynamic Coherent Structure Factor analysis.
       5      * DynamicCoherentStructureFactorARModel    : sets up a Dynamic Coherent Structure Factor analysis using an Auto Regressive model.
       6      * DynamicIncoherentStructureFactor         : sets up an Dynamic Incoherent Structure Factor analysis.
       7      * DynamicIncoherentStructureFactorGaussian : sets up an Dynamic Incoherent Structure Factor analysis using a Gaussian approximation.
       8      * IncoherentStructureFactorARModel         : sets up an Dynamic Incoherent Structure Factor analysis using an Auto Regressive model.
       9      * ElasticIncoherentStructureFactor         : sets up an Elastic Incoherent Structure Factor analysis.
      10      * StaticCoherentStructureFactor            : sets up a Static Coherent Structure Factor analysis.
      11
      12  Procedures:
      13      * DynamicStructureFactor : returns the Dynamic Structure Factor.
      14  """
      15
      16  # The python distribution modules
      17  import os
      18  import re
      19  from time import asctime
      20  from timeit import default_timer
      21
      22  # The ScientificPython modules
      23  from Scientific import N
      24  from Scientific.IO.NetCDF import _NetCDFFile, NetCDFFile
      25  from Scientific.Signals.Models import AutoRegressiveModel, AveragedAutoRegressiveModel
      26
      27  # The MMTK distribution modules
      28  from MMTK.Collections import Collection
      29  from MMTK_forcefield import NonbondedList
      30  from MMTK.Trajectory import Trajectory
      31
      32  # The nMOLDYN modules
      33  from nMOLDYN.Analysis.Analysis import Analysis, QVectors
      34  from nMOLDYN.Core.Error import Error
      35  from nMOLDYN.Core.IOFiles import convertNetCDFToASCII
      36  from nMOLDYN.Core.Mathematics import correlation, FFT, gaussianWindow
      37  from nMOLDYN_SSF import computeSSF
      38
      39  #####################################################################################
      40  # DYNAMIC STRUCTURE FACTOR ANALYSIS
      41  #####################################################################################
      42 -def DynamicStructureFactor(netcdf, alpha):
      43      """Computes the dynamic structure factor from an intermediate scattering function.
      44
      45      @param netcdf: the intermediate scattering function from which the dynamic structure factor will be computed..
      46      @type netcdf: string or instance of _NetCDFFile
      47
      48      @param alpha: the width, in percentage of the trajectory length, of the gaussian used in the smoothing procedure.
      49      @type alpha: float
      50      """
      51
      52      if isinstance(netcdf, str):
      53          try:
      54              SF = NetCDFFile(netcdf, 'a')
      55          except IOError:
      56              raise Error('The file %s could not be loaded.' % netcdf)
      57      else:
      58          return
      59
      60      # The tim NetCDF variable is extracted of the NC file.
      61      tim = SF.variables['time']
      62      dt = tim[1] - tim[0]
      63
      64      nQValues = len(SF.variables['q'][:])
      65
      66      # frequencies = 1D Numeric array. Contains the frequencies at which the dynamic structure is calculated
      67      frequencies = N.arange(len(tim))/(2.0*len(tim)*dt)
      68
      69      # Creation of the NetCDF dimension FREQUENCY. The number of frequency point present in the NC output file.
      70      SF.createDimension('NFREQUENCIES', len(frequencies[:]))
      71
      72      # Creation of the NetCDF variable |frequency| that will store the frequency at which the
      73      # dynamic structure factor is calculated
      74      FREQ  = SF.createVariable('frequency', N.Float, ('NFREQUENCIES',))
      75      FREQ.units = 'THz'
      76
      77      # FREQ is the subsampled version of |frequencies| list.
      78      FREQ[:] = frequencies[:]
      79
      80      # Is there some psf terms in the intermdiate scattering function input file ?
      81      partialTerms = [re.findall('Fqt-(.*)',v)[0] for v in SF.variables if 'Fqt-' in v]
      82      for pTerm in partialTerms:
      83          # Creation of the NetCDF variable |dsf| that will store the dynamic structure factor
      84          DSF = SF.createVariable('Sqw-'+pTerm, N.Float, ('NQVALUES','NFREQUENCIES'))
      85          # A loop is done for each q where the scattering function was calculated.
      86          for qVal in range(nQValues):
      87              DSF[qVal] = 0.5 * dt * FFT(gaussianWindow(SF.variables['Fqt-'+pTerm][qVal], alpha)).real[:len(frequencies)]
      88
      89      asciiVar = sorted(SF.variables.keys())
      90
      91      SF.close()
      92
      93      # Creates an ASCII version of the NetCDF output file.
      94      convertNetCDFToASCII(inputFile = netcdf,\
      95                           outputFile = os.path.splitext(netcdf)[0] + '.cdl',\
      96                           variables = asciiVar)
      97
      98  #####################################################################################
      99  # COHERENT STRUCTURE FACTOR ANALYSIS
     100  #####################################################################################
     101 -class DynamicCoherentStructureFactor(Analysis):
     102      """Sets up a Dynamic Coherent Structure Factor analysis.
     103
     104      A Subclass of nMOLDYN.Analysis.Analysis.
     105
     106      Constructor: DynamicCoherentStructureFactor(|parameters| = None)
     107
     108      Arguments:
     109
     110          - |parameters| -- a dictionnary of the input parameters, or 'None' to set up the analysis without parameters.
     111              * trajectory        -- a trajectory file name or an instance of MMTK.Trajectory.Trajectory class.
     112              * timeinfo          -- a string of the form 'first:last:step' where 'first' is an integer specifying the first frame
     113                                     number to consider, 'last' is an integer specifying the last frame number to consider and
     114                                     'step' is an integer specifying the step number between two frames.
     115              * qshellvalues      -- a string of the form 'qmin1:qmax1:dq1;qmin2:qmax2:dq2...' where 'qmin1', 'qmin2' ... ,
     116                                     'qmax1', 'qmax2' ... and 'dq1', 'dq2' ... are floats that represents respectively
     117                                     the q minimum, the q maximum and the q steps for q interval 1, 2 ...
     118              * qshellwidth       -- a float specifying the width of the q shells.
     119              * qvectorspershell  -- a float specifying the number of q vectors to generate per q shell.
     120              * qvectorsgenerator -- a string being one of 'isotropic', 'anisotropic' or 'explicit' specifying the way the q vectors
     121                                     will be generated.
     122              * qvectorsdirection -- a string of the form 'v1x,v1y,v1z;v2x,v2y,v2z...' where 'v1x', 'v2x' ..., 'v1y', 'v2y' ... and
     123                                     'v1z', 'v2z' ... are floats that represents respectively the x, y and z values of the vectord along
     124                                     which the q vectors should be generated.
     125              * fftwindow         -- a float in ]0.0,100.0[ specifying the width of the gaussian, in percentage of the trajectory length
     126                                     that will be used in the smoothing procedure.
     127              * subset            -- a selection string specifying the atoms to consider for the analysis.
     128              * deuteration       -- a selection string specifying the hydrogen atoms whose atomic parameters will be those of the deuterium.
     129              * weights           -- a string equal to 'equal', 'mass', 'coherent' , 'incoherent' or 'atomicNumber' that specifies the weighting
     130                                     scheme to use.
     131              * dcsf              -- the output NetCDF file name for the intermediate scattering function.
     132              * pyroserver        -- a string specifying if Pyro will be used and how to run the analysis.
     133
     134      Running modes:
     135
     136          - To run the analysis do: a.runAnalysis() where a is the analysis object.
     137          - To estimate the analysis do: a.estimateAnalysis() where a is the analysis object.
     138          - To save the analysis to 'file' file name do: a.saveAnalysis(file) where a is the analysis object.
     139
     140      """
     141
     142      # The minimal set of input parameters names necessary to perform the analysis.
     143      inputParametersNames = ('trajectory',\
     144                              'timeinfo',\
     145                              'qshellvalues',\
     146                              'qshellwidth',\
     147                              'qvectorspershell',\
     148                              'qvectorsgenerator',\
     149                              'qvectorsdirection',\
     150                              'fftwindow',\
     151                              'subset',\
     152                              'deuteration',\
     153                              'weights',\
     154                              'dcsf',\
     155                              'pyroserver',\
     156                              )
     157
     158      default = {'weights' : 'coherent'}
     159
     160      shortName = 'DCSF'
     161
     162      # Indicate whether this analysis can be estimated or not.
     163      canBeEstimated = False
     164
     165 -    def __init__(self):
     166          """The constructor. Insures that the class can not be instanciated directly from here.
     167          """
     168          raise Error('This class can not be instanciated.')
     169
     170 -    def initialize(self):
     171          """Initializes the analysis (e.g. parses and checks input parameters, set some variables ...).
     172          """
     173
     174          # The input parameters are parsed.
     175          self.parseInputParameters()
     176
     177          if self.pyroServer != 'monoprocessor':
     178              if self.statusBar is not None:
     179                  LogMessage('warning', 'The job status bar will be desactivated for compatibility with Pyro module.')
     180                  self.statusBar = None
     181
     182          self.buildTimeInfo()
     183
     184          self.subset = self.subsetSelection(self.universe, self.subset)
     185          self.nAtoms = self.subset.numberOfAtoms()
     186
     187          self.deuteration = self.deuterationSelection(self.universe, self.deuteration)
     188          self.deuteration = Collection(list(set(self.subset) & set(self.deuteration)))
     189
     190          self.weights = self.weightingScheme(self.universe, self.subset, self.deuteration, self.weights)
     191
     192          qv = QVectors(self.universe,\
     193                        self.qVectorsGenerator,\
     194                        self.qShellValues,\
     195                        self.qShellWidth,\
     196                        self.qVectorsPerShell,\
     197                        self.qVectorsDirection)
     198
     199          self.qRadii = qv.qRadii
     200          self.qVectors = qv.qVectors
     201          self.qVectorsStatistics = qv.statistics
     202
     203          self.nQValues = len(self.qRadii)
     204
     205          # Dictionnary between the atom symbols and their corresponding weight.
     206          self.symToWeight = {}
     207
     208          self.atomicSubsets = {}
     209
     210          for at in self.deuteration.atomList():
     211              if self.atomicSubsets.has_key('D'):
     212                  self.atomicSubsets['D'].addObject(at)
     213              else:
     214                  self.atomicSubsets['D'] = Collection(at)
     215                  self.symToWeight['D'] = self.weights[at]
     216
     217          for at in set(self.subset).difference(set(self.deuteration)):
     218
     219              if self.atomicSubsets.has_key(at.symbol):
     220                  self.atomicSubsets[at.symbol].addObject(at)
     221              else:
     222                  self.atomicSubsets[at.symbol] = Collection(at)
     223                  self.symToWeight[at.symbol] = self.weights[at]
     224
     225          self.FQT = {'total' : N.zeros((self.nQValues,self.nFrames), typecode = N.Float)}
     226
     227          # Loop over all the atom symbol.
     228          for symbol1 in self.atomicSubsets.keys():
     229              # Loop over all the atom symbol.
     230              for symbol2 in self.atomicSubsets.keys():
     231                  # The symbol pair tuple that will be used as the key for the histogram dictionnary.
     232                  pairName = tuple(sorted((symbol1,symbol2)))
     233                  # For each |pairName| key, the entry is a subdictionnary that stores the intra and intermolecular distances histograms.
     234                  self.FQT[pairName] = N.zeros((self.nQValues,self.nFrames), typecode = N.Float)
     235
     236          # Loop over the user-defined q values
     237          for qIndex in range(self.nQValues):
     238
     239              qRadius = self.qRadii[qIndex]
     240
     241              # For each q vector length, the corresponding set of q vector is transform into an array
     242              qarray = N.array(self.qVectors[qRadius])
     243              qarray = self.universe._boxToRealPointArray(qarray)
     244              # qVect[1][j] is now of shape (3,Qcount)
     245              self.qVectors[qRadius] = N.transpose(qarray)
     246
     247
     248          # Try to store all the frame if there is enough memory.
     249          try:
     250
     251              self.allConf = {}
     252              for symbol in self.atomicSubsets.keys():
     253                  self.allConf[symbol] = N.zeros((self.nFrames, self.atomicSubsets[symbol].numberOfAtoms(), 3), N.Float)
     254
     255              # loop over the trajectory time steps
     256              for frameIndex in range(self.nFrames):
     257
     258                  frame = self.frameIndexes[frameIndex]
     259
     260                  # conf contains the positions of all the atoms of the system for time i.
     261                  conf = self.trajectory.configuration[frame]
     262                  conf.convertToBoxCoordinates()
     263
     264                  for symbol in self.atomicSubsets.keys():
     265
     266                      mask = self.atomicSubsets[symbol].booleanMask()
     267                      self.allConf[symbol][frameIndex, :, :] = N.compress(mask.array, conf.array, 0)
     268
     269          except:
     270              self.allConf = None
     271
     272          if self.pyroServer != 'monoprocessor':
     273              # The attribute trajectory is removed because it can not be pickled by Pyro.
     274              delattr(self, 'trajectory')
     275
     276 -    def calc(self, qIndex, trajname):
     277          """Calculates the contribution for one Q-shell.
     278
     279          @param qIndex: the index of the Q-shell in |self.qRadii| list.
     280          @type qIndex: integer.
     281
     282          @param trajname: the name of the trajectory file name.
     283          @type trajname: string
     284          """
     285
     286          if self.allConf is None:
     287              if self.pyroServer == 'monoprocessor':
     288                  t = self.trajectory
     289
     290              else:
     291                  # Load the whole trajectory set.
     292                  t = Trajectory(None, trajname, 'r')
     293
     294          qRadius = self.qRadii[qIndex]
     295
     296          rho = {}
     297          for symbol in self.atomicSubsets.keys():
     298              rho[symbol] = N.zeros((self.nFrames, self.qVectors[qRadius].shape[1]), N.Complex)
     299
     300          # loop over the trajectory time steps
     301          for frameIndex in range(self.nFrames):
     302
     303              frame = self.frameIndexes[frameIndex]
     304
     305              if self.allConf is None:
     306
     307                  # conf contains the positions of all the atoms of the system for time i.
     308                  conf = t.configuration[frame]
     309                  conf.convertToBoxCoordinates()
     310
     311                  for symbol in self.atomicSubsets.keys():
     312                      mask = self.atomicSubsets[symbol].booleanMask()
     313                      selAtoms = N.compress(mask.array, conf.array, 0)
     314                      rho[symbol][frameIndex,:] = N.add.reduce(N.exp(1j*N.dot(selAtoms, self.qVectors[qRadius])))
     315
     316              else:
     317                  for symbol in self.atomicSubsets.keys():
     318                      rho[symbol][frameIndex,:] = N.add.reduce(N.exp(1j*N.dot(self.allConf[symbol][frameIndex,:,:], self.qVectors[qRadius])))
     319
     320          return rho
     321
     322 -    def combine(self, qIndex, x):
     323          """
     324          """
     325          qRadius = self.qRadii[qIndex]
     326          for symbol1 in self.atomicSubsets.keys():
     327              for symbol2 in self.atomicSubsets.keys():
     328                  pairName = tuple(sorted((symbol1,symbol2)))
     329                  corr = correlation(x[symbol1],x[symbol2])[:]
     330                  self.FQT[pairName][qIndex,:] += corr/self.qVectors[qRadius].shape[1]
     331
     332 -    def finalize(self):
     333          """Finalizes the calculations (e.g. averaging the total term, output files creations ...).
     334          """
     335
     336          for pairName in self.FQT.keys():
     337              if pairName == 'total':
     338                  continue
     339
     340              nAB = N.sqrt(self.atomicSubsets[pairName[0]].numberOfAtoms()*self.atomicSubsets[pairName[1]].numberOfAtoms())
     341  #            if pairName[0] != pairName[1]:
     342  #                nAB = 2.0*nAB
     343              self.FQT[pairName] = self.FQT[pairName]/nAB
     344              wAB = self.symToWeight[pairName[0]]*self.symToWeight[pairName[1]]
     345              N.add(self.FQT['total'], nAB*wAB*self.FQT[pairName], self.FQT['total'])
     346
     347          # The NetCDF output file intermediate scattering function.
     348          outputFile       = NetCDFFile(self.outputDCSF, 'w')
     349          # The title for the NetCDF file.
     350          outputFile.title = self.__class__.__name__
     351          outputFile.jobinfo = self.information + '\nOutput file written on: %s\n\n' % asctime()
     352
     353          # Some NetCDF dimensions.
     354          outputFile.createDimension('NQVALUES', self.nQValues)
     355          outputFile.createDimension('NTIMES', self.nFrames)
     356          outputFile.createDimension('NOCTANS', 8)
     357          outputFile.createDimension('OCTANNAME', 8)
     358
     359          # Creation of the NetCDF variable |sf| that will store the structure factor
     360          TIMES    = outputFile.createVariable('time', N.Float, ('NTIMES',))
     361          TIMES[:] = self.times[:]
     362          TIMES.units = 'ps'
     363
     364          QLENGTH = outputFile.createVariable('q', N.Float, ('NQVALUES',))
     365          QLENGTH[:] = self.qRadii
     366          QLENGTH.units = 'nm^-1'
     367
     368          OCTANS = outputFile.createVariable('octan', N.Character, ('NOCTANS','OCTANNAME'))
     369          OCTANS[:,:] = N.array([list('X+.Y+.Z+'),list('X+.Y+.Z-'),list('X+.Y-.Z+'),list('X+.Y-.Z-'),\
     370                                  list('X-.Y+.Z+'),list('X-.Y+.Z-'),list('X-.Y-.Z+'),list('X-.Y-.Z-')])
     371          OCTANS.units = 'unitless'
     372
     373          STAT = outputFile.createVariable('qvectors_statistics', N.Int, ('NQVALUES', 'NOCTANS'))
     374          STAT[:,:] = self.qVectorsStatistics
     375
     376          for k in self.FQT.keys():
     377              if isinstance(k,tuple):
     378                  DCSF = outputFile.createVariable('Fqt-%s' % ''.join(k), N.Float, ('NQVALUES','NTIMES'))
     379              else:
     380                  DCSF = outputFile.createVariable('Fqt-%s' % k, N.Float, ('NQVALUES','NTIMES'))
     381              DCSF[:] = self.FQT[k]
     382              DCSF.units = 'unitless'
     383
     384          outputFile.close()
     385
     386          self.toPlot = {'netcdf' : self.outputDCSF, 'xVar' : 'q', 'yVar' : 'time', 'zVar' : 'Fqt-total'}
     387
     388          DynamicStructureFactor(self.outputDCSF, self.fftWindow)
     389
     390  #####################################################################################
     391  # STATIC COHERENT STRUCTURE FACTOR ANALYSIS
     392  #####################################################################################
     393 -class StaticCoherentStructureFactor(Analysis):
     394      """Sets up a Coherent Structure Factor analysis.
     395
     396      A Subclass of nMOLDYN.Analysis.Analysis.
     397
     398      Constructor: StaticCoherentStructureFactor(|parameters| = None)
     399
     400      Arguments:
     401
     402          - |parameters| -- a dictionnary of the input parameters, or 'None' to set up the analysis without parameters.
     403              * trajectory        -- a trajectory file name or an instance of MMTK.Trajectory.Trajectory class.
     404              * timeinfo          -- a string of the form 'first:last:step' where 'first' is an integer specifying the first frame
     405                                     number to consider, 'last' is an integer specifying the last frame number to consider and
     406                                     'step' is an integer specifying the step number between two frames.
     407              * qshellvalues      -- a string of the form 'qmin1:qmax1:dq1;qmin2:qmax2:dq2...' where 'qmin1', 'qmin2' ... ,
     408                                     'qmax1', 'qmax2' ... and 'dq1', 'dq2' ... are floats that represents respectively
     409                                     the q minimum, the q maximum and the q steps for q interval 1, 2 ...
     410              * qshellwidth       -- a float specifying the width of the q shells.
     411              * qvectorspershell  -- a float specifying the number of q vectors to generate per q shell.
     412              * qvectorsgenerator -- a string being one of 'isotropic', 'anisotropic' or 'explicit' specifying the way the q vectors
     413                                     will be generated.
     414              * qvectorsdirection -- a string of the form 'v1x,v1y,v1z;v2x,v2y,v2z...' where 'v1x', 'v2x' ..., 'v1y', 'v2y' ... and
     415                                     'v1z', 'v2z' ... are floats that represents respectively the x, y and z values of the vectord along
     416                                     which the q vectors should be generated.
     417              * fftwindow         -- a float in ]0.0,100.0[ specifying the width of the gaussian, in percentage of the trajectory length
     418                                     that will be used in the smoothing procedure.
     419              * subset            -- a selection string specifying the atoms to consider for the analysis.
     420              * deuteration       -- a selection string specifying the hydrogen atoms whose atomic parameters will be those of the deuterium.
     421              * weights           -- a string equal to 'equal', 'mass', 'coherent' , 'incoherent' or 'atomicNumber' that specifies the weighting
     422                                     scheme to use.
     423              * csf               -- the output NetCDF file name for the intermediate scattering function.
     424              * pyroserver        -- a string specifying if Pyro will be used and how to run the analysis.
     425
     426      Running modes:
     427
     428          - To run the analysis do: a.runAnalysis() where a is the analysis object.
     429          - To estimate the analysis do: a.estimateAnalysis() where a is the analysis object.
     430          - To save the analysis to 'file' file name do: a.saveAnalysis(file) where a is the analysis object.
     431
     432      """
     433
     434      # The minimal set of input parameters names necessary to perform the analysis.
     435      inputParametersNames = ('trajectory',\
     436                              'timeinfo',\
     437                              'qshellvalues',\
     438                              'qshellwidth',\
     439                              'qvectorspershell',\
     440                              'qvectorsgenerator',\
     441                              'qvectorsdirection',\
     442                              'subset',\
     443                              'deuteration',\
     444                              'weights',\
     445                              'scsf',\
     446                              'pyroserver',\
     447                              )
     448
     449      default = {'weights' : 'coherent'}
     450
     451      shortName = 'SCSF'
     452
     453      # Indicate whether this analysis can be estimated or not.
     454      canBeEstimated = False
     455
     456 -    def __init__(self):
     457          """The constructor. Insures that the class can not be instanciated directly from here.
     458          """
     459          raise Error('This class can not be instanciated.')
     460
     461 -    def initialize(self):
     462          """Initializes the analysis (e.g. parses and checks input parameters, set some variables ...).
     463          """
     464
     465          # The input parameters are parsed.
     466          self.parseInputParameters()
     467
     468          if self.pyroServer != 'monoprocessor':
     469              if self.statusBar is not None:
     470                  LogMessage('warning', 'The job status bar will be desactivated for compatibility with Pyro module.')
     471                  self.statusBar = None
     472
     473          self.buildTimeInfo()
     474
     475          self.subset = self.subsetSelection(self.universe, self.subset)
     476          self.nAtoms = self.subset.numberOfAtoms()
     477
     478          self.deuteration = self.deuterationSelection(self.universe, self.deuteration)
     479          self.deuteration = Collection(list(set(self.subset) & set(self.deuteration)))
     480
     481          self.weights = self.weightingScheme(self.universe, self.subset, self.deuteration, self.weights)
     482
     483          qv = QVectors(self.universe,\
     484                        self.qVectorsGenerator,\
     485                        self.qShellValues,\
     486                        self.qShellWidth,\
     487                        self.qVectorsPerShell,\
     488                        self.qVectorsDirection)
     489
     490          self.qRadii = qv.qRadii
     491          self.qVectors = qv.qVectors
     492          self.qVectorsStatistics = qv.statistics
     493
     494          self.nQValues = len(self.qRadii)
     495
     496          # Dictionnary between the atom symbols and their corresponding weight.
     497          self.symToWeight = {}
     498
     499          self.atomicSubsets = {}
     500
     501          for at in self.deuteration.atomList():
     502              if self.atomicSubsets.has_key('D'):
     503                  self.atomicSubsets['D'].addObject(at)
     504              else:
     505                  self.atomicSubsets['D'] = Collection(at)
     506                  self.symToWeight['D'] = self.weights[at]
     507
     508          for at in set(self.subset).difference(set(self.deuteration)):
     509
     510              if self.atomicSubsets.has_key(at.symbol):
     511                  self.atomicSubsets[at.symbol].addObject(at)
     512              else:
     513                  self.atomicSubsets[at.symbol] = Collection(at)
     514                  self.symToWeight[at.symbol] = self.weights[at]
     515
     516          self.SQ = {'total' : N.zeros((self.nQValues,), typecode = N.Float)}
     517
     518          # Loop over all the atom symbol.
     519          for symbol1 in self.atomicSubsets.keys():
     520              # Loop over all the atom symbol.
     521              for symbol2 in self.atomicSubsets.keys():
     522                  # The symbol pair tuple that will be used as the key for the histogram dictionnary.
     523                  pairName = tuple(sorted((symbol1,symbol2)))
     524                  # For each |pairName| key, the entry is a subdictionnary that stores the intra and intermolecular distances histograms.
     525                  self.SQ[pairName] = N.zeros((self.nQValues,), typecode = N.Float)
     526
     527          # Loop over the user-defined q values
     528          for qIndex in range(self.nQValues):
     529
     530              qRadius = self.qRadii[qIndex]
     531
     532              # For each q vector length, the corresponding set of q vector is transform into an array
     533              qarray = N.array(self.qVectors[qRadius])
     534              qarray = self.universe._boxToRealPointArray(qarray)
     535              # qVect[1][j] is now of shape (3,Qcount)
     536              self.qVectors[qRadius] = N.transpose(qarray)
     537
     538          try:
     539              self.allConf = {}
     540              for symbol in self.atomicSubsets.keys():
     541                  self.allConf[symbol] = N.zeros((self.nFrames, self.atomicSubsets[symbol].numberOfAtoms(), 3), N.Float)
     542
     543              # loop over the trajectory time steps
     544              for frameIndex in range(self.nFrames):
     545
     546                  frame = self.frameIndexes[frameIndex]
     547
     548                  # conf contains the positions of all the atoms of the system for time i.
     549                  conf = self.trajectory.configuration[frame]
     550                  conf.convertToBoxCoordinates()
     551
     552                  for symbol in self.atomicSubsets.keys():
     553
     554                      mask = self.atomicSubsets[symbol].booleanMask()
     555                      self.allConf[symbol][frameIndex, :, :] = N.compress(mask.array, conf.array, 0)
     556
     557          except:
     558              self.allConf = None
     559
     560          if self.pyroServer != 'monoprocessor':
     561              # The attribute trajectory is removed because it can not be pickled by Pyro.
     562              delattr(self, 'trajectory')
     563
     564 -    def calc(self, qIndex, trajname):
     565          """Calculates the contribution for one Q-shell.
     566
     567          @param qIndex: the index of the Q-shell in |self.qRadii| list.
     568          @type qIndex: integer.
     569
     570          @param trajname: the name of the trajectory file name.
     571          @type trajname: string
     572          """
     573
     574          if self.allConf is None:
     575              if self.pyroServer == 'monoprocessor':
     576                  t = self.trajectory
     577
     578              else:
     579                  # Load the whole trajectory set.
     580                  t = Trajectory(None, trajname, 'r')
     581
     582          qRadius = self.qRadii[qIndex]
     583
     584          rho = {}
     585          for symbol in self.atomicSubsets.keys():
     586              rho[symbol] = N.zeros((self.nFrames, self.qVectors[qRadius].shape[1]), N.Complex)
     587
     588          # loop over the trajectory time steps
     589          for frameIndex in range(self.nFrames):
     590
     591              frame = self.frameIndexes[frameIndex]
     592
     593              if self.allConf is None:
     594                  # conf contains the positions of all the atoms of the system for time i.
     595                  conf = t.configuration[frame]
     596                  conf.convertToBoxCoordinates()
     597
     598                  for symbol in self.atomicSubsets.keys():
     599
     600                      mask = self.atomicSubsets[symbol].booleanMask()
     601                      selAtoms = N.compress(mask.array, conf.array, 0)
     602                      rho[symbol][frameIndex,:] = N.add.reduce(N.exp(1j*N.dot(selAtoms, self.qVectors[qRadius])))
     603              else:
     604                  for symbol in self.atomicSubsets.keys():
     605                      rho[symbol][frameIndex,:] = N.add.reduce(N.exp(1j*N.dot(self.allConf[symbol][frameIndex,:,:], self.qVectors[qRadius])))
     606
     607          return rho
     608
     609 -    def combine(self, qIndex, x):
     610          """
     611          """
     612
     613          qRadius = self.qRadii[qIndex]
     614          for symbol1 in self.atomicSubsets.keys():
     615              for symbol2 in self.atomicSubsets.keys():
     616                  pairName = tuple(sorted((symbol1,symbol2)))
     617
     618                  rhoirhoj = N.conjugate(x[symbol1]) * x[symbol2]
     619                  corrt0 = N.add.reduce(N.add.reduce(rhoirhoj,1)).real/self.nFrames
     620
     621                  self.SQ[pairName][qIndex] += corrt0/self.qVectors[qRadius].shape[1]
     622
     623 -    def finalize(self):
     624          """Finalizes the calculations (e.g. averaging the total term, output files creations ...).
     625          """
     626
     627          for pairName in self.SQ.keys():
     628              if pairName == 'total':
     629                  continue
     630
     631              nAB = N.sqrt(self.atomicSubsets[pairName[0]].numberOfAtoms()*self.atomicSubsets[pairName[1]].numberOfAtoms())
     632              if pairName[0] != pairName[1]:
     633                  nAB = 2.0*nAB
     634              self.SQ[pairName] = self.SQ[pairName]/nAB
     635              wAB = self.symToWeight[pairName[0]]*self.symToWeight[pairName[1]]
     636              N.add(self.SQ['total'], nAB*wAB*self.SQ[pairName], self.SQ['total'])
     637
     638          # The NetCDF output file intermediate scattering function.
     639          outputFile       = NetCDFFile(self.outputSCSF, 'w')
     640          # The title for the NetCDF file.
     641          outputFile.title = self.__class__.__name__
     642          outputFile.jobinfo = self.information + '\nOutput file written on: %s\n\n' % asctime()
     643
     644          # Some NetCDF dimensions.
     645          outputFile.createDimension('NQVALUES', None)
     646          outputFile.createDimension('NOCTANS', 8)
     647          outputFile.createDimension('OCTANNAME', 8)
     648
     649          QLENGTH = outputFile.createVariable('q', N.Float, ('NQVALUES',))
     650          QLENGTH[:] = self.qRadii
     651          QLENGTH.units = 'nm^-1'
     652
     653          OCTANS = outputFile.createVariable('octan', N.Character, ('NOCTANS','OCTANNAME'))
     654          OCTANS[:,:] = N.array([list('X+.Y+.Z+'),list('X+.Y+.Z-'),list('X+.Y-.Z+'),list('X+.Y-.Z-'),\
     655                                  list('X-.Y+.Z+'),list('X-.Y+.Z-'),list('X-.Y-.Z+'),list('X-.Y-.Z-')])
     656          OCTANS.units = 'unitless'
     657
     658          STAT = outputFile.createVariable('qvectors_statistics', N.Int, ('NQVALUES', 'NOCTANS'))
     659          STAT[:,:] = self.qVectorsStatistics
     660
     661          for k in self.SQ.keys():
     662              if isinstance(k,tuple):
     663                  SCSF = outputFile.createVariable('Sq-%s' % ''.join(k), N.Float, ('NQVALUES',))
     664              else:
     665                  SCSF = outputFile.createVariable('Sq-%s' % k, N.Float, ('NQVALUES',))
     666              SCSF[:] = self.SQ[k]
     667              SCSF.units = 'unitless'
     668
     669          outputFile.close()
     670
     671          self.toPlot = {'netcdf' : self.outputSCSF, 'xVar' : 'q', 'yVar' : 'Sq-total'}
     672
     673  #####################################################################################
     674  # COHERENT STRUCTURE FACTOR ANALYSIS WITHIN AR MODEL FRAMEWORK
     675  #####################################################################################
     676 -class DynamicCoherentStructureFactorAR(Analysis):
     677      """Sets up a Dynamic Coherent Structure Factor analysis using an Auto Regressive model.
     678
     679      A Subclass of nMOLDYN.Analysis.Analysis.
     680
     681      Constructor: DynamicCoherentStructureFactorARModel(|parameters| = None)
     682
     683      Arguments:
     684
     685          - |parameters| -- a dictionnary of the input parameters, or 'None' to set up the analysis without parameters.
     686              * trajectory        -- a trajectory file name or an instance of MMTK.Trajectory.Trajectory class.
     687              * timeinfo          -- a string of the form 'first:last:step' where 'first' is an integer specifying the first frame
     688                                     number to consider, 'last' is an integer specifying the last frame number to consider and
     689                                     'step' is an integer specifying the step number between two frames.
     690              * armodelorder      -- an integer in [1, len(trajectory)[ specifying the order of the model
     691              * qshellvalues      -- a string of the form 'qmin1:qmax1:dq1;qmin2:qmax2:dq2...' where 'qmin1', 'qmin2' ... ,
     692                                     'qmax1', 'qmax2' ... and 'dq1', 'dq2' ... are floats that represents respectively
     693                                     the q minimum, the q maximum and the q steps for q interval 1, 2 ...
     694              * qshellwidth       -- a float specifying the width of the q shells.
     695              * qvectorspershell  -- a float specifying the number of q vectors to generate per q shell.
     696              * qvectorsgenerator -- a string being one of 'isotropic', 'anisotropic' or 'explicit' specifying the way the q vectors
     697                                     will be generated.
     698              * qvectorsdirection -- a string of the form 'v1x,v1y,v1z;v2x,v2y,v2z...' where 'v1x', 'v2x' ..., 'v1y', 'v2y' ... and
     699                                     'v1z', 'v2z' ... are floats that represents respectively the x, y and z values of the vectord along
     700                                     which the q vectors should be generated.
     701              * subset            -- a selection string specifying the atoms to consider for the analysis.
     702              * deuteration       -- a selection string specifying the hydrogen atoms whose atomic parameters will be those of the deuterium.
     703              * weights           -- a string equal to 'equal', 'mass', 'coherent' , 'incoherent' or 'atomicNumber' that specifies the weighting
     704                                     scheme to use.
     705              * dcsfar            -- the output NetCDF file name.
     706              * pyroserver        -- a string specifying if Pyro will be used and how to run the analysis.
     707
     708      Running modes:
     709
     710          - To run the analysis do: a.runAnalysis() where a is the analysis object.
     711          - To estimate the analysis do: a.estimateAnalysis() where a is the analysis object.
     712          - To save the analysis to 'file' file name do: a.saveAnalysis(file) where a is the analysis object.
     713
     714      """
     715
     716      # The minimal set of input parameters names necessary to perform the analysis.
     717      inputParametersNames = ('trajectory',\
     718                              'timeinfo',\
     719                              'armodelorder',\
     720                              'qshellvalues',\
     721                              'qshellwidth',\
     722                              'qvectorspershell',\
     723                              'qvectorsgenerator',\
     724                              'qvectorsdirection',\
     725                              'subset',\
     726                              'deuteration',\
     727                              'weights',\
     728                              'dcsfar',\
     729                              'pyroserver',\
     730                              )
     731
     732      shortName = 'DCSFAR'
     733
     734      # Indicate whether this analysis can be estimated or not.
     735      canBeEstimated = False
     736
     737      default = {'weights' : 'coherent'}
     738
     739 -    def __init__(self):
     740          """The constructor. Insures that the class can not be instanciated directly from here.
     741          """
     742          raise Error('This class can not be instanciated.')
     743
     744 -    def initialize(self):
     745          """Initializes the analysis (e.g. parses and checks input parameters, set some variables ...).
     746          """
     747
     748          # The input parameters are parsed.
     749          self.parseInputParameters()
     750
     751          if self.pyroServer != 'monoprocessor':
     752              if self.statusBar is not None:
     753                  LogMessage('warning', 'The job status bar will be desactivated for compatibility with Pyro module.')
     754                  self.statusBar = None
     755
     756          self.buildTimeInfo()
     757
     758          self.subset = self.subsetSelection(self.universe, self.subset)
     759          self.nAtoms = self.subset.numberOfAtoms()
     760
     761          self.deuteration = self.deuterationSelection(self.universe, self.deuteration)
     762          self.deuteration = Collection(list(set(self.subset) & set(self.deuteration)))
     763
     764          self.weights = self.weightingScheme(self.universe, self.subset, self.deuteration, self.weights)
     765
     766          if (self.arModelOrder <= 0) or (self.arModelOrder >= self.nFrames):
     767              raise Error('The AR order must be an integer in [1,%d[.' % self.nFrames)
     768
     769          qv = QVectors(self.universe,\
     770                        self.qVectorsGenerator,\
     771                        self.qShellValues,\
     772                        self.qShellWidth,\
     773                        self.qVectorsPerShell,\
     774                        self.qVectorsDirection)
     775
     776          self.qRadii = qv.qRadii
     777          self.qVectors = qv.qVectors
     778          self.qVectorsStatistics = qv.statistics
     779
     780          self.nQValues = len(self.qRadii)
     781
     782          # Loop over the user-defined q values
     783          for qIndex in range(self.nQValues):
     784
     785              qRadius = self.qRadii[qIndex]
     786
     787              # For each q vector length, the corresponding set of q vector is transform into an array
     788              qarray = N.array(self.qVectors[qRadius])
     789              qarray = self.universe._boxToRealPointArray(qarray)
     790              # qVect[1][j] is now of shape (3,Qcount)
     791              self.qVectors[qRadius] = N.transpose(qarray)
     792
     793          # The frequency values.
     794          freqMax = 1.0/(2.0*self.dt)
     795          df = freqMax/self.nFrames
     796          self.frequencies = N.arange(self.nFrames + 1) * df
     797
     798          self.modelRealPart = N.zeros((self.nQValues, self.arModelOrder + 2), typecode = N.Float)
     799          self.modelImagPart = N.zeros((self.nQValues, self.arModelOrder + 2), typecode = N.Float)
     800          self.FQT = N.zeros((self.nQValues, self.nFrames), typecode = N.Float)
     801          self.FQTMemory = N.zeros((self.nQValues, self.arModelOrder + self.arModelOrder/2), typecode = N.Float)
     802          self.SQW = N.zeros((self.nQValues, self.nFrames + 1), typecode = N.Float)
     803
     804          try:
     805              self.allConf = N.zeros((self.nFrames, self.nAtoms, 3), N.Float)
     806
     807              # loop over the trajectory time steps
     808              for frameIndex in range(self.nFrames):
     809
     810                  frame = self.frameIndexes[frameIndex]
     811
     812                  # conf contains the positions of all the atoms of the system for time i.
     813                  conf = self.trajectory.configuration[frame]
     814                  conf.convertToBoxCoordinates()
     815
     816                  # mask = 1D Numeric array. 0 for the unselected atoms, 1 for the selected atoms.
     817                  mask = self.subset.booleanMask()
     818                  self.allConf[frameIndex, :, :] = N.compress(mask.array, conf.array, 0)
     819
     820          except:
     821              self.allConf = None
     822
     823          if self.pyroServer != 'monoprocessor':
     824              # The attribute trajectory is removed because it can not be pickled by Pyro.
     825              delattr(self, 'trajectory')
     826
     827 -    def calc(self, qIndex, trajname):
     828          """Calculates the contribution for one Q-shell.
     829
     830          @param qIndex: the index of the Q-shell in |self.qRadii| list.
     831          @type qIndex: integer.
     832
     833          @param trajname: the name of the trajectory file name.
     834          @type trajname: string
     835          """
     836
     837          if self.pyroServer == 'monoprocessor':
     838              t = self.trajectory
     839
     840          else:
     841              # Load the whole trajectory set.
     842              t = Trajectory(None, trajname, 'r')
     843
     844          qRadius = self.qRadii[qIndex]
     845
     846          model = AveragedAutoRegressiveModel(self.arModelOrder, self.dt)
     847
     848          # mask = 1D Numeric array. 0 for the unselected atoms, 1 for the selected atoms.
     849          mask = self.subset.booleanMask()
     850          weights = N.compress(mask.array, self.weights.array)[:, N.NewAxis]
     851
     852          rho = N.zeros((self.nFrames, self.qVectors[qRadius].shape[1]), N.Complex)
     853          for frameIndex in range(self.nFrames):
     854              frame = self.frameIndexes[frameIndex]
     855
     856              if self.allConf is None:
     857                  conf = t.configuration[frame]
     858                  conf.convertToBoxCoordinates()
     859                  selAtoms = N.compress(mask.array, conf.array, 0)
     860                  rho[frameIndex,:] = N.add.reduce(weights*N.exp(1j*N.dot(selAtoms, self.qVectors[qRadius])))
     861              else:
     862                  rho[frameIndex,:] = N.add.reduce(weights*N.exp(1j*N.dot(self.allConf[frameIndex,:,:], self.qVectors[qRadius])))
     863
     864          return rho
     865
     866 -    def combine(self, qIndex, x):
     867          """
     868          """
     869
     870          for i in range(x.shape[1]):
     871              data = x[:, i]
     872              data = data - N.add.reduce(data)/len(data)
     873              m = AutoRegressiveModel(self.arModelOrder, data, self.dt)
     874              model.add(m)
     875
     876          parameters = N.concatenate((model.coeff[::-1], N.array([model.sigma, model.variance])))
     877          self.modelRealPart[qIndex, :] = parameters.real
     878          self.modelImagPart[qIndex, :] = parameters.imag
     879
     880          average = N.add.reduce(N.add.reduce(x))/N.multiply.reduce(x.shape)
     881          self.FQT[qIndex,:] = model.correlation(self.nFrames).values.real + (average*N.conjugate(average)).real
     882
     883          mem = model.memoryFunction(self.arModelOrder + self.arModelOrder/2).values.real
     884          self.FQTMemory[qIndex,:] = mem[:]
     885
     886          spectrum = model.spectrum(2.0*N.pi*self.frequencies)
     887          self.SQW[qIndex,:] = spectrum.values
     888
     889 -    def finalize(self):
     890          """Finalizes the calculations (e.g. averaging the total term, output files creations ...).
     891          """
     892
     893          outputFile = NetCDFFile(self.outputDCSFAR, 'w')
     894          outputFile.title = '%s (order %d)' % (self.__class__.__name__, self.arModelOrder)
     895          outputFile.jobinfo = self.information + '\nOutput file written on: %s\n\n' % asctime()
     896
     897          outputFile.createDimension('NQVALUES', self.nQValues)
     898          outputFile.createDimension('NTIMES', self.nFrames)
     899          outputFile.createDimension('NTIMES_MEMORY', self.arModelOrder + self.arModelOrder/2)
     900          outputFile.createDimension('NPOLES', self.arModelOrder + 2)
     901          outputFile.createDimension('NFREQUENCIES', self.nFrames + 1)
     902          outputFile.createDimension('NOCTANS', 8)
     903          outputFile.createDimension('OCTANNAME', 8)
     904
     905          OCTANS = outputFile.createVariable('octan', N.Character, ('NOCTANS','OCTANNAME'))
     906          OCTANS[:,:] = N.array([list('X+.Y+.Z+'),list('X+.Y+.Z-'),list('X+.Y-.Z+'),list('X+.Y-.Z-'),\
     907                                  list('X-.Y+.Z+'),list('X-.Y+.Z-'),list('X-.Y-.Z+'),list('X-.Y-.Z-')])
     908          OCTANS.units = 'unitless'
     909
     910          STAT = outputFile.createVariable('qvectors_statistics', N.Int, ('NQVALUES', 'NOCTANS'))
     911          STAT[:,:] = self.qVectorsStatistics
     912
     913          # Creation of the NetCDF variable |sf| that will store the structure factor
     914          TIMES = outputFile.createVariable('time', N.Float, ('NTIMES',))
     915          TIMES[:] = self.times[:]
     916          TIMES.units = 'ps'
     917
     918          QLENGTH = outputFile.createVariable('q', N.Float, ('NQVALUES',))
     919          QLENGTH[:] = self.qRadii[:]
     920          QLENGTH.units = 'nm^-1'
     921
     922          TIMEMEMORY = outputFile.createVariable('time_memory', N.Float, ('NTIMES_MEMORY',))
     923          TIMEMEMORY[:] = self.dt * N.arange(self.arModelOrder + self.arModelOrder/2)
     924          TIMEMEMORY.units = 'ps'
     925
     926          FREQUENCY = outputFile.createVariable('frequency', N.Float, ('NFREQUENCIES',))
     927          FREQUENCY[:] = self.frequencies[:]
     928          FREQUENCY.units = 'THz'
     929
     930          MODELORDER = outputFile.createVariable('n', N.Int, ('NPOLES',))
     931          MODELORDER[:] = N.arange(0, self.arModelOrder + 2)
     932          MODELORDER.units = 'unitless'
     933
     934          MODELREAL = outputFile.createVariable('ar_coefficients_real', N.Float, ('NQVALUES', 'NPOLES'))
     935          MODELREAL[:,:] = self.modelRealPart[:,:]
     936          MODELREAL.units = 'unitless'
     937
     938          MODELIMAG = outputFile.createVariable('ar_coefficients_imag', N.Float, ('NQVALUES', 'NPOLES'))
     939          MODELIMAG[:,:] = self.modelImagPart[:,:]
     940          MODELIMAG.units = 'unitless'
     941
     942          FQT     = outputFile.createVariable('Fqt', N.Float, ('NQVALUES','NTIMES'))
     943          FQT[:,:] = self.FQT[:,:]
     944          FQT.units = 'unitless'
     945
     946          FQTMEM  = outputFile.createVariable('Fqt_memory', N.Float, ('NQVALUES','NTIMES_MEMORY'))
     947          FQTMEM[:,:] = self.FQTMemory[:,:]
     948          FQTMEM.units = 'unitless'
     949
     950          SQW = outputFile.createVariable('Sqw', N.Float, ('NQVALUES','NFREQUENCIES'))
     951          SQW[:,:] = self.SQW[:,:]
     952          SQW.units = 'unitless'
     953
     954          outputFile.close()
     955
     956          self.toPlot = None
     957
     958  #####################################################################################
     959  # DYNAMIC INCOHERENT STRUCTURE FACTOR ANALYSIS
     960  #####################################################################################
     961 -class DynamicIncoherentStructureFactor(Analysis):
     962      """Sets up an Dynamic Incoherent Structure Factor analysis.
     963
     964      A Subclass of nMOLDYN.Analysis.Analysis.
     965
     966      Constructor: DynamicIncoherentStructureFactorARModel(|parameters| = None)
     967
     968      Arguments:
     969
     970          - |parameters| -- a dictionnary of the input parameters, or 'None' to set up the analysis without parameters.
     971              * trajectory        -- a trajectory file name or an instance of MMTK.Trajectory.Trajectory class.
     972              * timeinfo          -- a string of the form 'first:last:step' where 'first' is an integer specifying the first frame
     973                                     number to consider, 'last' is an integer specifying the last frame number to consider and
     974                                     'step' is an integer specifying the step number between two frames.
     975              * qshellvalues      -- a string of the form 'qmin1:qmax1:dq1;qmin2:qmax2:dq2...' where 'qmin1', 'qmin2' ... ,
     976                                     'qmax1', 'qmax2' ... and 'dq1', 'dq2' ... are floats that represents respectively
     977                                     the q minimum, the q maximum and the q steps for q interval 1, 2 ...
     978              * qshellwidth       -- a float specifying the width of the q shells.
     979              * qvectorspershell  -- a float specifying the number of q vectors to generate per q shell.
     980              * qvectorsgenerator -- a string being one of 'isotropic', 'anisotropic' or 'explicit' specifying the way the q vectors
     981                                     will be generated.
     982              * qvectorsdirection -- a string of the form 'v1x,v1y,v1z;v2x,v2y,v2z...' where 'v1x', 'v2x' ..., 'v1y', 'v2y' ... and
     983                                     'v1z', 'v2z' ... are floats that represents respectively the x, y and z values of the vectord along
     984                                     which the q vectors should be generated.
     985              * fftwindow         -- a float in ]0.0,100.0[ specifying the width of the gaussian, in percentage of the trajectory length
     986                                     that will be used in the smoothing procedure.
     987              * subset            -- a selection string specifying the atoms to consider for the analysis.
     988              * deuteration       -- a selection string specifying the hydrogen atoms whose atomic parameters will be those of the deuterium.
     989              * weights           -- a string equal to 'equal', 'mass', 'coherent' , 'incoherent' or 'atomicNumber' that specifies the weighting
     990                                     scheme to use.
     991              * disf              -- the output NetCDF file name for the intermediate scattering function.
     992              * pyroserver        -- a string specifying if Pyro will be used and how to run the analysis.
     993
     994      Running modes:
     995
     996          - To run the analysis do: a.runAnalysis() where a is the analysis object.
     997          - To estimate the analysis do: a.estimateAnalysis() where a is the analysis object.
     998          - To save the analysis to 'file' file name do: a.saveAnalysis(file) where a is the analysis object.
     999
    1000      """
    1001
    1002      # The minimal set of input parameters names necessary to perform the analysis.
    1003      inputParametersNames = ('trajectory',\
    1004                              'timeinfo',\
    1005                              'qshellvalues',\
    1006                              'qshellwidth',\
    1007                              'qvectorspershell',\
    1008                              'qvectorsgenerator', \
    1009                              'qvectorsdirection',\
    1010                              'fftwindow',\
    1011                              'subset',\
    1012                              'deuteration',\
    1013                              'weights',\
    1014                              'disf',\
    1015                              'pyroserver',\
    1016                              )
    1017
    1018      default = {'weights' : 'incoherent'}
    1019
    1020      shortName = 'DISF'
    1021
    1022      # Indicate whether this analysis can be estimated or not.
    1023      canBeEstimated = True
    1024
    1025 -    def __init__(self):
    1026          """The constructor. Insures that the class can not be instanciated directly from here.
    1027          """
    1028
    1029          raise Error('This class can not be instanciated.')
    1030
    1031 -    def initialize(self):
    1032          """Initializes the analysis (e.g. parses and checks input parameters, set some variables ...).
    1033          """
    1034
    1035          # The input parameters are parsed.
    1036          self.parseInputParameters()
    1037
    1038          if self.pyroServer != 'monoprocessor':
    1039              if self.statusBar is not None:
    1040                  LogMessage('warning', 'The job status bar will be desactivated for compatibility with Pyro module.')
    1041                  self.statusBar = None
    1042
    1043          self.buildTimeInfo()
    1044
    1045          self.subset = self.subsetSelection(self.universe, self.subset)
    1046          self.nAtoms = self.subset.numberOfAtoms()
    1047
    1048          self.deuteration = self.deuterationSelection(self.universe, self.deuteration)
    1049          self.deuteration = Collection(list(set(self.subset) & set(self.deuteration)))
    1050
    1051          self.weights = self.weightingScheme(self.universe, self.subset, self.deuteration, self.weights)
    1052
    1053          qv = QVectors(self.universe,\
    1054                        self.qVectorsGenerator,\
    1055                        self.qShellValues,\
    1056                        self.qShellWidth,\
    1057                        self.qVectorsPerShell,\
    1058                        self.qVectorsDirection)
    1059
    1060          self.qRadii = qv.qRadii
    1061          self.qVectors = qv.qVectors
    1062          self.qVectorsStatistics = qv.statistics
    1063
    1064          self.nQValues = len(self.qRadii)
    1065
    1066          self.preLoadTrajectory('atom')
    1067
    1068          # Dictionnary between the atom symbols and the index of the corresponding atoms in the selected atoms for the analysis.
    1069          self.symToInd = {}
    1070
    1071          # Dictionnary between the atom symbols and their corresponding weight.
    1072          self.symToWeight = {}
    1073
    1074          # Dictionnary between the atom index and their corresponding atom symbol.
    1075          self.indToSym = {}
    1076
    1077          for at in self.deuteration.atomList():
    1078              self.indToSym[at.index] = 'D'
    1079              if self.symToInd.has_key('D'):
    1080                  self.symToInd['D'].append(at.index)
    1081              else:
    1082                  self.symToInd['D'] = [at.index]
    1083                  self.symToWeight['D'] = self.weights[at]
    1084
    1085          for at in set(self.subset).difference(set(self.deuteration)):
    1086              self.indToSym[at.index] = at.symbol
    1087              if self.symToInd.has_key(at.symbol):
    1088                  self.symToInd[at.symbol].append(at.index)
    1089              else:
    1090                  self.symToInd[at.symbol] = [at.index]
    1091                  self.symToWeight[at.symbol] = self.weights[at]
    1092
    1093          self.ISF = {}
    1094
    1095          # Loop over all the atom symbol.
    1096          for symbol in self.symToInd.keys():
    1097              # For each |pairName| key, the entry is a subdictionnary that stores the intra and intermolecular distances histograms.
    1098              self.ISF[symbol] = N.zeros((self.nQValues,self.nFrames), typecode = N.Float)
    1099
    1100          for qIndex in range(len(self.qRadii)):
    1101              qRadius = self.qRadii[qIndex]
    1102              # For each q vector length, the corresponding set of q vector is transform into an array
    1103              qArray = N.array(self.qVectors[qRadius])
    1104              self.qVectors[qRadius] = N.transpose(qArray)
    1105
    1106          if self.pyroServer != 'monoprocessor':
    1107              # The attribute trajectory is removed because it can not be pickled by Pyro.
    1108              delattr(self, 'trajectory')
    1109
    1110 -    def calc(self, atom, trajname):
    1111          """Calculates the atomic term.
    1112
    1113          @param atom: the atom on which the atomic term has been calculated.
    1114          @type atom: an instance of MMTK.Atom class.
    1115
    1116          @param trajname: the name of the trajectory file name.
    1117          @type trajname: string
    1118          """
    1119
    1120          if self.preLoadedTraj is None:
    1121              if self.pyroServer == 'monoprocessor':
    1122                  t = self.trajectory
    1123
    1124              else:
    1125                  # Load the whole trajectory set.
    1126                  t = Trajectory(None, trajname, 'r')
    1127
    1128              # series = 2D Numeric array. The positions of the selected atom |at| from the first step to the
    1129              # last step with the selected step increment.
    1130              series = t.readParticleTrajectory(atom, first = self.first, last = self.last, skip = self.skip).array
    1131
    1132          else:
    1133              series = self.preLoadedTraj[atom]
    1134
    1135          symbol = self.indToSym[atom.index]
    1136
    1137          atomicISF = N.zeros((self.nQValues,self.nFrames), typecode = N.Float)
    1138
    1139          for qIndex in range(self.nQValues):
    1140              qRadius = self.qRadii[qIndex]
    1141
    1142              expSeries = N.exp(1j*N.dot(series, self.qVectors[qRadius]))
    1143
    1144              # The ensemble average of (exp(-iqR(0))*exp(iqR(t))) is replaced by an autocorrelation.
    1145              res = correlation(expSeries, expSeries)[:self.nFrames]/self.qVectors[qRadius].shape[1]
    1146
    1147              N.add(atomicISF[qIndex,:], res, atomicISF[qIndex,:])
    1148
    1149          if self.preLoadedTraj is not None:
    1150              del self.preLoadedTraj[atom]
    1151
    1152          return atomicISF
    1153
    1154 -    def combine(self, atom, x):
    1155          """
    1156          """
    1157          symbol = self.indToSym[atom.index]
    1158          N.add(self.ISF[symbol], x, self.ISF[symbol])
    1159
    1160 -    def finalize(self):
    1161          """Finalizes the calculations (e.g. averaging the total term, output files creations ...)
    1162          """
    1163
    1164          self.ISF['total'] = N.zeros((self.nQValues,self.nFrames), typecode = N.Float)
    1165
    1166          for symbol in self.symToInd.keys():
    1167              self.ISF[symbol] = self.ISF[symbol]/float(len(self.symToInd[symbol]))
    1168              N.add(self.ISF['total'],\
    1169                    float(len(self.symToInd[symbol]))*self.symToWeight[symbol]*self.ISF[symbol],\
    1170                    self.ISF['total'])
    1171
    1172          # The NetCDF output file is opened for writing.
    1173          outputFile       = NetCDFFile(self.outputDISF, 'w')
    1174          outputFile.title = self.__class__.__name__
    1175          outputFile.jobinfo = self.information + '\nOutput file written on: %s\n\n' % asctime()
    1176
    1177          # Some NetCDF dimensions.
    1178          outputFile.createDimension('NQVALUES', self.nQValues)
    1179          outputFile.createDimension('NTIMES', self.nFrames)
    1180          outputFile.createDimension('NOCTANS', 8)
    1181          outputFile.createDimension('OCTANNAME', 8)
    1182
    1183          QLENGTH = outputFile.createVariable('q', N.Float, ('NQVALUES',))
    1184          QLENGTH[:] = self.qRadii
    1185          QLENGTH.units = 'nm^-1'
    1186
    1187          TIMES = outputFile.createVariable('time', N.Float, ('NTIMES',))
    1188          TIMES[:] = self.times[:]
    1189          TIMES.units = 'ps'
    1190
    1191          OCTANS = outputFile.createVariable('octan', N.Character, ('NOCTANS','OCTANNAME'))
    1192          OCTANS[:,:] = N.array([list('X+.Y+.Z+'),list('X+.Y+.Z-'),list('X+.Y-.Z+'),list('X+.Y-.Z-'),\
    1193                                  list('X-.Y+.Z+'),list('X-.Y+.Z-'),list('X-.Y-.Z+'),list('X-.Y-.Z-')])
    1194          OCTANS.units = 'unitless'
    1195
    1196          STAT = outputFile.createVariable('qvectors_statistics', N.Int, ('NQVALUES', 'NOCTANS'))
    1197          STAT[:,:] = self.qVectorsStatistics
    1198
    1199          for k in self.ISF.keys():
    1200              DISF = outputFile.createVariable('Fqt-%s' % k, N.Float, ('NQVALUES','NTIMES'))
    1201              DISF[:,:] = self.ISF[k]
    1202              DISF.units = 'unitless'
    1203
    1204          # The NetCDF output file is closed.
    1205          outputFile.close()
    1206
    1207          self.toPlot = {'netcdf' : self.outputDISF, 'xVar' : 'q', 'yVar' : 'time', 'zVar' : 'Fqt-total'}
    1208
    1209          DynamicStructureFactor(self.outputDISF, self.fftWindow)
    1210
    1211  #####################################################################################
    1212  # DYNAMIC INCOHERENT STRUCTURE FACTOR WITHIN AR MODEL FRAMEWORK
    1213  #####################################################################################
    1214 -class DynamicIncoherentStructureFactorAR(Analysis):
    1215      """Sets up an Dynamic Incoherent Structure Factor analysis using an Auto Regressive model.
    1216
    1217      A Subclass of nMOLDYN.Analysis.Analysis.
    1218
    1219      Constructor: DynamicIncoherentStructureFactorARModel(|parameters| = None)
    1220
    1221      Arguments:
    1222
    1223          - |parameters| -- a dictionnary of the input parameters, or 'None' to set up the analysis without parameters.
    1224              * trajectory        -- a trajectory file name or an instance of MMTK.Trajectory.Trajectory class.
    1225              * timeinfo          -- a string of the form 'first:last:step' where 'first' is an integer specifying the first frame
    1226                                     number to consider, 'last' is an integer specifying the last frame number to consider and
    1227                                     'step' is an integer specifying the step number between two frames.
    1228              * armodelorder      -- an integer in [1, len(trajectory)[ specifying the order of the model
    1229              * qshellvalues      -- a string of the form 'qmin1:qmax1:dq1;qmin2:qmax2:dq2...' where 'qmin1', 'qmin2' ... ,
    1230                                     'qmax1', 'qmax2' ... and 'dq1', 'dq2' ... are floats that represents respectively
    1231                                     the q minimum, the q maximum and the q steps for q interval 1, 2 ...
    1232              * qshellwidth       -- a float specifying the width of the q shells.
    1233              * qvectorspershell  -- a float specifying the number of q vectors to generate per q shell.
    1234              * qvectorsgenerator -- a string being one of 'isotropic', 'anisotropic' or 'explicit' specifying the way the q vectors
    1235                                     will be generated.
    1236              * qvectorsdirection -- a string of the form 'v1x,v1y,v1z;v2x,v2y,v2z...' where 'v1x', 'v2x' ..., 'v1y', 'v2y' ... and
    1237                                     'v1z', 'v2z' ... are floats that represents respectively the x, y and z values of the vectord along
    1238                                     which the q vectors should be generated.
    1239              * subset            -- a selection string specifying the atoms to consider for the analysis.
    1240              * deuteration       -- a selection string specifying the hydrogen atoms whose atomic parameters will be those of the deuterium.
    1241              * weights           -- a string equal to 'equal', 'mass', 'coherent' , 'incoherent' or 'atomicNumber' that specifies the weighting
    1242                                     scheme to use.
    1243              * disfar            -- the output NetCDF file name for the intermediate scattering function.
    1244              * pyroserver        -- a string specifying if Pyro will be used and how to run the analysis.
    1245
    1246      Running modes:
    1247
    1248          - To run the analysis do: a.runAnalysis() where a is the analysis object.
    1249          - To estimate the analysis do: a.estimateAnalysis() where a is the analysis object.
    1250          - To save the analysis to 'file' file name do: a.saveAnalysis(file) where a is the analysis object.
    1251
    1252      """
    1253
    1254      # The minimal set of input parameters names necessary to perform the analysis.
    1255      inputParametersNames = ('trajectory',\
    1256                              'timeInfo',\
    1257                              'armodelorder',\
    1258                              'qshellvalues',\
    1259                              'qshellwidth',\
    1260                              'qvectorspershell',
    1261                              'qvectorsgenerator',\
    1262                              'qvectorsdirection',\
    1263                              'subset',\
    1264                              'deuteration',\
    1265                              'weights',\
    1266                              'disfar',\
    1267                              'pyroserver',\
    1268                              )
    1269
    1270      shortName = 'DISFAR'
    1271
    1272      # Indicate whether this analysis can be estimated or not.
    1273      canBeEstimated = True
    1274
    1275      default = {'weights' : 'incoherent'}
    1276
    1277 -    def __init__(self):
    1278          """The constructor. Insures that the class can not be instanciated directly from here.
    1279          """
    1280          raise Error('This class can not be instanciated.')
    1281
    1282 -    def initialize(self):
    1283          """Initializes the analysis (e.g. parses and checks input parameters, set some variables ...).
    1284          """
    1285
    1286          # The input parameters are parsed.
    1287          self.parseInputParameters()
    1288
    1289          if self.pyroServer != 'monoprocessor':
    1290              if self.statusBar is not None:
    1291                  LogMessage('warning', 'The job status bar will be desactivated for compatibility with Pyro module.')
    1292                  self.statusBar = None
    1293
    1294          self.buildTimeInfo()
    1295
    1296          self.subset = self.subsetSelection(self.universe, self.subset)
    1297          self.nAtoms = self.subset.numberOfAtoms()
    1298
    1299          self.deuteration = self.deuterationSelection(self.universe, self.deuteration)
    1300          self.deuteration = Collection(list(set(self.subset) & set(self.deuteration)))
    1301
    1302          self.weights = self.weightingScheme(self.universe, self.subset, self.deuteration, self.weights)
    1303
    1304          if (self.arModelOrder <= 0) or (self.arModelOrder >= self.nFrames):
    1305              raise Error('The order of the AR model must be in [1,%s[' % self.nFrames)
    1306
    1307          qv = QVectors(self.universe,\
    1308                        self.qVectorsGenerator,\
    1309                        self.qShellValues,\
    1310                        self.qShellWidth,\
    1311                        self.qVectorsPerShell,\
    1312                        self.qVectorsDirection)
    1313
    1314          self.qRadii = qv.qRadii
    1315          self.qVectors = qv.qVectors
    1316          self.qVectorsStatistics = qv.statistics
    1317
    1318          self.nQValues = len(self.qRadii)
    1319
    1320          # Loop over the user-defined q values
    1321          for qIndex in range(self.nQValues):
    1322
    1323              qRadius = self.qRadii[qIndex]
    1324
    1325              # For each q vector length, the corresponding set of q vector is transform into an array
    1326              qarray = N.array(self.qVectors[qRadius])
    1327              # qVect[1][j] is now of shape (3,Qcount)
    1328              self.qVectors[qRadius] = N.transpose(qarray)
    1329
    1330          # Dictionnary that will store the AR model for each Q-shell.
    1331          self.ARModel = {}
    1332          # Loop over the user-defined q values
    1333          for qIndex in range(self.nQValues):
    1334              qRadius = self.qRadii[qIndex]
    1335              self.ARModel[qRadius] = AveragedAutoRegressiveModel(self.arModelOrder, self.dt)
    1336
    1337          # The frequency values.
    1338          freqMax = 1.0/(2.0*self.dt)
    1339          df = freqMax/self.nFrames
    1340          self.frequencies = N.arange(self.nFrames + 1) * df
    1341
    1342          self.modelRealPart = N.zeros((self.nQValues, self.arModelOrder + 2), typecode = N.Float)
    1343          self.modelImagPart = N.zeros((self.nQValues, self.arModelOrder + 2), typecode = N.Float)
    1344          self.FQT = N.zeros((self.nQValues, self.nFrames), typecode = N.Float)
    1345          self.FQTMemory = N.zeros((self.nQValues, self.arModelOrder + self.arModelOrder/2), typecode = N.Float)
    1346          self.SQW = N.zeros((self.nQValues, self.nFrames + 1), typecode = N.Float)
    1347
    1348          if self.pyroServer != 'monoprocessor':
    1349              # The attribute trajectory is removed because it can not be pickled by Pyro.
    1350              delattr(self, 'trajectory')
    1351
    1352 -    def calc(self, atom, trajname):
    1353          """Calculates the atomic term.
    1354
    1355          @param atom: the atom on which the atomic term has been calculated.
    1356          @type atom: an instance of MMTK.Atom class.
    1357
    1358          @param trajname: the name of the trajectory file name.
    1359          @type trajname: string
    1360          """
    1361
    1362          if self.weights[atom] == 0.0:
    1363              return
    1364
    1365          if self.pyroServer == 'monoprocessor':
    1366              t = self.trajectory
    1367
    1368          else:
    1369              # Load the whole trajectory set.
    1370              t = Trajectory(None, trajname, 'r')
    1371
    1372          series = t.readParticleTrajectory(atom, first = self.first, last = self.last, skip = self.skip).array
    1373
    1374          atomicFQT = N.zeros((self.nQValues,self.nFrames), typecode = N.Float)
    1375          atomicSQW = N.zeros((self.nQValues,self.nFrames + 1), typecode = N.Float)
    1376
    1377          for qIndex in range(self.nQValues):
    1378              qRadius = self.qRadii[qIndex]
    1379              rho = N.exp(1j*N.dot(series, self.qVectors[qRadius]))
    1380              rho_av = N.add.reduce(rho)/rho.shape[0]
    1381              qModel = AveragedAutoRegressiveModel(self.arModelOrder, self.dt)
    1382
    1383              for i in range(rho.shape[1]):
    1384                  data = rho[:, i] - rho_av[i]
    1385                  m = AutoRegressiveModel(self.arModelOrder, data, self.dt)
    1386                  qModel.add(m, self.weights[atom])
    1387                  self.ARModel[qRadius].add(m, self.weights[atom])
    1388
    1389              rho_av_sq = (rho_av*N.conjugate(rho_av)).real
    1390              average = N.add.reduce(rho_av_sq)/rho.shape[1]
    1391
    1392              res = self.weights[atom]*(qModel.correlation(self.nFrames).values.real + average)
    1393              N.add(atomicFQT[qIndex,:], res, atomicFQT[qIndex,:])
    1394
    1395              res = self.weights[atom]*(qModel.spectrum(2.0*N.pi*self.frequencies).values)
    1396              N.add(atomicSQW[qIndex,:], res, atomicSQW[qIndex,:])
    1397
    1398          return atomicFQT, atomicSQW
    1399
    1400 -    def combine(self, atom, x):
    1401          """
    1402          """
    1403
    1404          N.add(self.FQT, x[0], self.FQT)
    1405          N.add(self.SQW, x[1], self.SQW)
    1406
    1407 -    def finalize(self):
    1408          """Finalizes the calculations (e.g. averaging the total term, output files creations ...).
    1409          """
    1410
    1411          for qIndex in range(self.nQValues):
    1412              qRadius = self.qRadii[qIndex]
    1413
    1414              parameters = N.concatenate((self.ARModel[qRadius].coeff[::-1],\
    1415                                          N.array([self.ARModel[qRadius].sigma, self.ARModel[qRadius].variance])))
    1416
    1417              self.modelRealPart[qIndex, :] = parameters.real
    1418              self.modelImagPart[qIndex, :] = parameters.imag
    1419
    1420              mem = self.ARModel[qRadius].memoryFunction(self.arModelOrder + self.arModelOrder/2).values.real
    1421              self.FQTMemory[qIndex,:] = mem
    1422
    1423          outputFile = NetCDFFile(self.outputDISFAR, 'w')
    1424          outputFile.title = '%s (order %d)' % (self.__class__.__name__, self.arModelOrder)
    1425          outputFile.jobinfo = self.information + '\nOutput file written on: %s\n\n' % asctime()
    1426
    1427          outputFile.createDimension('NQVALUES', self.nQValues)
    1428          outputFile.createDimension('NTIMES', self.nFrames)
    1429          outputFile.createDimension('NTIMES_MEMORY', self.arModelOrder + self.arModelOrder/2)
    1430          outputFile.createDimension('NPOLES', self.arModelOrder + 2)
    1431          outputFile.createDimension('NFREQUENCIES', self.nFrames + 1)
    1432          outputFile.createDimension('NOCTANS', 8)
    1433          outputFile.createDimension('OCTANNAME', 8)
    1434
    1435          OCTANS = outputFile.createVariable('octan', N.Character, ('NOCTANS','OCTANNAME'))
    1436          OCTANS[:,:] = N.array([list('X+.Y+.Z+'),list('X+.Y+.Z-'),list('X+.Y-.Z+'),list('X+.Y-.Z-'),\
    1437                                  list('X-.Y+.Z+'),list('X-.Y+.Z-'),list('X-.Y-.Z+'),list('X-.Y-.Z-')])
    1438          OCTANS.units = 'unitless'
    1439
    1440          STAT = outputFile.createVariable('qvectors_statistics', N.Int, ('NQVALUES', 'NOCTANS'))
    1441          STAT[:,:] = self.qVectorsStatistics
    1442
    1443          # Creation of the NetCDF variable |sf| that will store the structure factor
    1444          TIMES = outputFile.createVariable('time', N.Float, ('NTIMES',))
    1445          TIMES[:] = self.times[:]
    1446          TIMES.units = 'ps'
    1447
    1448          QLENGTH = outputFile.createVariable('q', N.Float, ('NQVALUES',))
    1449          QLENGTH[:] = self.qRadii[:]
    1450          QLENGTH.units = 'nm^-1'
    1451
    1452          TIMEMEMORY = outputFile.createVariable('time_memory', N.Float, ('NTIMES_MEMORY',))
    1453          TIMEMEMORY[:] = self.dt * N.arange(self.arModelOrder + self.arModelOrder/2)
    1454          TIMEMEMORY.units = 'ps'
    1455
    1456          FREQUENCY = outputFile.createVariable('frequency', N.Float, ('NFREQUENCIES',))
    1457          FREQUENCY[:] = self.frequencies[:]
    1458          FREQUENCY.units = 'THz'
    1459
    1460          MODELORDER = outputFile.createVariable('n', N.Int, ('NPOLES',))
    1461          MODELORDER[:] = N.arange(0, self.arModelOrder + 2)
    1462          MODELORDER.units = 'unitless'
    1463
    1464          MODELREAL = outputFile.createVariable('ar_model_real', N.Float, ('NQVALUES', 'NPOLES'))
    1465          MODELREAL[:,:] = self.modelRealPart[:,:]
    1466          MODELREAL.units = 'unitless'
    1467
    1468          MODELIMAG = outputFile.createVariable('ar_model_imag', N.Float, ('NQVALUES', 'NPOLES'))
    1469          MODELIMAG[:,:] = self.modelImagPart[:,:]
    1470          MODELIMAG.units = 'unitless'
    1471
    1472          FQT     = outputFile.createVariable('Fqt', N.Float, ('NQVALUES','NTIMES'))
    1473          FQT[:,:] = self.FQT[:,:]
    1474          FQT.units = 'unitless'
    1475
    1476          FQTMEM  = outputFile.createVariable('Fqt_memory', N.Float, ('NQVALUES','NTIMES_MEMORY'))
    1477          FQTMEM[:,:] = self.FQTMemory[:,:]
    1478          FQTMEM.units = 'unitless'
    1479
    1480          SQW = outputFile.createVariable('Sqw', N.Float, ('NQVALUES','NFREQUENCIES'))
    1481          SQW[:,:] = self.SQW[:,:]
    1482          SQW.units = 'unitless'
    1483
    1484          outputFile.close()
    1485
    1486          self.toPlot = None
    1487
    1488  #####################################################################################
    1489  # DYNAMIC INCOHERENT STRUCTURE FACTOR ANALYSIS USING GAUSSIAN APPROXIMATION
    1490  #####################################################################################
    1491 -class DynamicIncoherentStructureFactorGaussian(Analysis):
    1492      """Sets up an Dynamic Incoherent Structure Factor analysis within Gaussian approximation.
    1493
    1494      A Subclass of nMOLDYN.Analysis.Analysis.
    1495
    1496      Constructor: DynamicIncoherentStructureFactorGaussian(|parameters| = None)
    1497
    1498      Arguments:
    1499
    1500          - |parameters| -- a dictionnary of the input parameters, or 'None' to set up the analysis without parameters.
    1501              * trajectory   -- a trajectory file name or an instance of MMTK.Trajectory.Trajectory class.
    1502              * timeinfo     -- a string of the form 'first:last:step' where 'first' is an integer specifying the first frame
    1503                                number to consider, 'last' is an integer specifying the last frame number to consider and
    1504                                'step' is an integer specifying the step number between two frames.
    1505              * qshellvalues -- a string of the form 'qmin1:qmax1:dq1;qmin2:qmax2:dq2...' where 'qmin1', 'qmin2' ... ,
    1506                                'qmax1', 'qmax2' ... and 'dq1', 'dq2' ... are floats that represents respectively
    1507                                the q minimum, the q maximum and the q steps for q interval 1, 2 ...
    1508              * fftwindow    -- a float in ]0.0,100.0[ specifying the width of the gaussian, in percentage of the trajectory length
    1509                                that will be used in the smoothing procedure.
    1510              * subset       -- a selection string specifying the atoms to consider for the analysis.
    1511              * deuteration  -- a selection string specifying the hydrogen atoms whose atomic parameters will be those of the deuterium.
    1512              * weights      -- a string equal to 'equal', 'mass', 'coherent' , 'incoherent' or 'atomicNumber' that specifies the weighting
    1513                                scheme to use.
    1514              * disfg        -- the output NetCDF file name for the intermediate scattering function.
    1515              * pyroserver   -- a string specifying if Pyro will be used and how to run the analysis.
    1516
    1517      Running modes:
    1518
    1519          - To run the analysis do: a.runAnalysis() where a is the analysis object.
    1520          - To estimate the analysis do: a.estimateAnalysis() where a is the analysis object.
    1521          - To save the analysis to 'file' file name do: a.saveAnalysis(file) where a is the analysis object.
    1522
    1523      """
    1524
    1525      # The minimal set of input parameters names necessary to perform the analysis.
    1526      inputParametersNames = ('trajectory',\
    1527                              'timeinfo',\
    1528                              'qshellvalues',\
    1529                              'fftwindow',\
    1530                              'subset',\
    1531                              'deuteration',\
    1532                              'weights',\
    1533                              'disfg',\
    1534                              'pyroserver',\
    1535                              )
    1536
    1537      default = {'weights' : 'incoherent'}
    1538
    1539      shortName = 'DISFG'
    1540
    1541      # Indicate whether this analysis can be estimated or not.
    1542      canBeEstimated = True
    1543
    1544 -    def __init__(self):
    1545          """The constructor. Insures that the class can not be instanciated directly from here.
    1546          """
    1547
    1548          raise Error('This class can not be instanciated.')
    1549
    1550 -    def initialize(self):
    1551          """Initializes the analysis (e.g. parses and checks input parameters, set some variables ...).
    1552          """
    1553
    1554          # The input parameters are parsed.
    1555          self.parseInputParameters()
    1556
    1557          if self.pyroServer != 'monoprocessor':
    1558              if self.statusBar is not None:
    1559                  LogMessage('warning', 'The job status bar will be desactivated for compatibility with Pyro module.')
    1560                  self.statusBar = None
    1561
    1562          self.buildTimeInfo()
    1563
    1564          self.subset = self.subsetSelection(self.universe, self.subset)
    1565          self.nAtoms = self.subset.numberOfAtoms()
    1566
    1567          self.deuteration = self.deuterationSelection(self.universe, self.deuteration)
    1568          self.deuteration = Collection(list(set(self.subset) & set(self.deuteration)))
    1569
    1570          self.weights = self.weightingScheme(self.universe, self.subset, self.deuteration, self.weights)
    1571
    1572          self.nQValues = len(self.qShellValues)
    1573          # |qSquare| is an array whose values are the squared of the q values list. These are the q^2 values used in
    1574          # equation 3.35
    1575          self.qSquare    = N.array(self.qShellValues)*N.array(self.qShellValues)
    1576
    1577          self.preLoadTrajectory('atom')
    1578
    1579          # Dictionnary between the atom symbols and the index of the corresponding atoms in the selected atoms for the analysis.
    1580          self.symToInd = {}
    1581
    1582          # Dictionnary between the atom symbols and their corresponding weight.
    1583          self.symToWeight = {}
    1584
    1585          # Dictionnary between the atom index and their corresponding atom symbol.
    1586          self.indToSym = {}
    1587
    1588          for at in self.deuteration.atomList():
    1589              self.indToSym[at.index] = 'D'
    1590              if self.symToInd.has_key('D'):
    1591                  self.symToInd['D'].append(at.index)
    1592
    1593              else:
    1594                  self.symToInd['D'] = [at.index]
    1595                  self.symToWeight['D'] = self.weights[at]
    1596
    1597          for at in set(self.subset).difference(set(self.deuteration)):
    1598              self.indToSym[at.index] = at.symbol
    1599              if self.symToInd.has_key(at.symbol):
    1600                  self.symToInd[at.symbol].append(at.index)
    1601              else:
    1602                  self.symToInd[at.symbol] = [at.index]
    1603                  self.symToWeight[at.symbol] = self.weights[at]
    1604
    1605          self.ISFG = {}
    1606
    1607          # Loop over all the atom symbol.
    1608          for symbol in self.symToInd.keys():
    1609              # For each |pairName| key, the entry is a subdictionnary that stores the intra and intermolecular distances histograms.
    1610              self.ISFG[symbol] = N.zeros((self.nQValues,self.nFrames), typecode = N.Float)
    1611
    1612          if self.pyroServer != 'monoprocessor':
    1613              # The attribute trajectory is removed because it can not be pickled by Pyro.
    1614              delattr(self, 'trajectory')
    1615
    1616 -    def calc(self, atom, trajname):
    1617          """Calculates the atomic term.
    1618
    1619          @param atom: the atom on which the atomic term has been calculated.
    1620          @type atom: an instance of MMTK.Atom class.
    1621
    1622          @param trajname: the name of the trajectory file name.
    1623          @type trajname: string
    1624          """
    1625
    1626          if self.preLoadedTraj is None:
    1627              if self.pyroServer == 'monoprocessor':
    1628                  t = self.trajectory
    1629
    1630              else:
    1631                  # Load the whole trajectory set.
    1632                  t = Trajectory(None, trajname, 'r')
    1633
    1634              # series = 2D Numeric array. The positions of the selected atom |at| from the first step to the
    1635              # last step with the selected step increment.
    1636              series = t.readParticleTrajectory(atom, first = self.first, last = self.last, skip = self.skip).array
    1637
    1638          else:
    1639              series = self.preLoadedTraj[atom]
    1640
    1641          symbol = self.indToSym[atom.index]
    1642
    1643          atomicISFG = N.zeros((self.nQValues,self.nFrames), typecode = N.Float)
    1644
    1645          res  = self.getMSD(series)[:self.nFrames]
    1646
    1647          # This is the formula 3.36. The factor 6 is in fact 2*3 where 3 account for the isotropic calculation.
    1648          for qVal in range(self.nQValues):
    1649              gaussian = N.exp(-res*self.qSquare[qVal]/6.)
    1650
    1651              N.add(atomicISFG[qVal,:], gaussian, atomicISFG[qVal,:])
    1652
    1653          if self.preLoadedTraj is not None:
    1654              del self.preLoadedTraj[atom]
    1655
    1656          return atomicISFG
    1657
    1658 -    def combine(self, atom, x):
    1659          """
    1660          """
    1661          symbol = self.indToSym[atom.index]
    1662          N.add(self.ISFG[symbol], x, self.ISFG[symbol])
    1663
    1664 -    def finalize(self):
    1665          """Finalizes the calculations (e.g. averaging the total term, output files creations ...)
    1666          """
    1667
    1668          self.ISFG['total'] = N.zeros((self.nQValues,self.nFrames), typecode = N.Float)
    1669          for symbol in self.symToInd.keys():
    1670              self.ISFG[symbol] = self.ISFG[symbol]/float(len(self.symToInd[symbol]))
    1671              N.add(self.ISFG['total'],\
    1672                    float(len(self.symToInd[symbol]))*self.symToWeight[symbol]*self.ISFG[symbol],\
    1673                    self.ISFG['total'])
    1674
    1675          # The NetCDF output file intermediate scattering function.
    1676          outputFile       = NetCDFFile(self.outputDISFG, 'w')
    1677          # The title for the NetCDF file.
    1678          outputFile.title = self.__class__.__name__
    1679          outputFile.jobinfo = self.information + '\nOutput file written on: %s\n\n' % asctime()
    1680
    1681          # Some NetCDF dimensions.
    1682          outputFile.createDimension('NQVALUES', self.nQValues)
    1683          outputFile.createDimension('NTIMES', self.nFrames)
    1684
    1685          # Creation of the NetCDF variable |sf| that will store the structure factor
    1686          Qlength = outputFile.createVariable('q', N.Float, ('NQVALUES',))
    1687          Q2length = outputFile.createVariable('q2', N.Float, ('NQVALUES',))
    1688          Time    = outputFile.createVariable('time', N.Float, ('NTIMES',))
    1689
    1690          Qlength[:] = N.sqrt(self.qSquare)
    1691          Q2length[:] = self.qSquare
    1692          Time[:] = self.times[:]
    1693
    1694          Time.units = 'ps'
    1695          Qlength.units = 'nm^-1'
    1696
    1697          for key, value in self.ISFG.items():
    1698              DISFG = outputFile.createVariable('Fqt-%s' % key, N.Float, ('NQVALUES','NTIMES'))
    1699              DISFG[:,:] = value
    1700              DISFG.units = 'unitless'
    1701
    1702          # The NetCDF output file is closed.
    1703          outputFile.close()
    1704
    1705          self.toPlot = {'netcdf' : self.outputDISFG, 'xVar' : 'q', 'yVar' : 'time', 'zVar' : 'Fqt-total'}
    1706
    1707          DynamicStructureFactor(self.outputDISFG, self.fftWindow)
    1708
    1709 -    def getMSD(self, series):
    1710          """
    1711          Computes the atomic component of the Mean-Square-Displacement.
    1712          This is the exact copy of the version written in nMOLDYN.Simulations.Dynamics but rewritten here
    1713          for to keep the module Scattering independant from module Dynamics.
    1714          """
    1715
    1716          # dsq is the squared norm of the position for each time step
    1717          # dsq refers to DSQ(k) in the published algorithm
    1718          dsq = N.add.reduce(series * series,1)
    1719
    1720          # sum_dsq1 is the cumulative sum of dsq
    1721          sum_dsq1 = N.add.accumulate(dsq)
    1722
    1723          # sum_dsq1 is the reversed cumulative sum of dsq
    1724          sum_dsq2 = N.add.accumulate(dsq[::-1])
    1725
    1726          # sumsq refers to SUMSQ in the published algorithm
    1727          sumsq = 2.*sum_dsq1[-1]
    1728
    1729          # this line refers to the instruction SUMSQ <-- SUMSQ - DSQ(m-1) - DSQ(N - m) of the published algorithm
    1730          # In this case, msd is an array because the instruction is computed for each m ranging from 0 to len(traj) - 1
    1731          # So, this single instruction is performing the loop in the published algorithm
    1732          Saabb  = sumsq - N.concatenate(([0.], sum_dsq1[:-1])) - N.concatenate(([0.], sum_dsq2[:-1]))
    1733
    1734          # Saabb refers to SAA+BB/(N-m) in the published algorithm
    1735          # Sab refers to SAB(m)/(N-m) in the published algorithm
    1736          Saabb = Saabb / (len(dsq) - N.arange(len(dsq)))
    1737          Sab   = 2.*correlation(series)
    1738
    1739          # this line refers to the instruction MSD(m) <-- SUMSQ - 2Sab(m) in the published algorithm
    1740          return (Saabb - Sab)
    1741
    1742  #####################################################################################
    1743  # ELASTIC INCOHERENT STRUCTURE FACTOR
    1744  #####################################################################################
    1745 -class ElasticIncoherentStructureFactor(Analysis):
    1746      """Sets up an Elastic Incoherent Structure Factor.
    1747
    1748      A Subclass of nMOLDYN.Analysis.Analysis.
    1749
    1750      Constructor: ElasticIncoherentStructureFactor(|parameters| = None)
    1751
    1752      Arguments:
    1753
    1754          - |parameters| -- a dictionnary of the input parameters, or 'None' to set up the analysis without parameters.
    1755              * trajectory        -- a trajectory file name or an instance of MMTK.Trajectory.Trajectory class.
    1756              * timeinfo          -- a string of the form 'first:last:step' where 'first' is an integer specifying the first frame
    1757                                     number to consider, 'last' is an integer specifying the last frame number to consider and
    1758                                     'step' is an integer specifying the step number between two frames.
    1759              * qshellvalues      -- a string of the form 'qmin1:qmax1:dq1;qmin2:qmax2:dq2...' where 'qmin1', 'qmin2' ... ,
    1760                                     'qmax1', 'qmax2' ... and 'dq1', 'dq2' ... are floats that represents respectively
    1761                                     the q minimum, the q maximum and the q steps for q interval 1, 2 ...
    1762              * qshellwidth       -- a float specifying the width of the q shells.
    1763              * qvectorspershell  -- a float specifying the number of q vectors to generate per q shell.
    1764              * qvectorsgenerator -- a string being one of 'isotropic', 'anisotropic' or 'explicit' specifying the way the q vectors
    1765                                     will be generated.
    1766              * qvectorsdirection -- a string of the form 'v1x,v1y,v1z;v2x,v2y,v2z...' where 'v1x', 'v2x' ..., 'v1y', 'v2y' ... and
    1767                                     'v1z', 'v2z' ... are floats that represents respectively the x, y and z values of the vectord along
    1768                                     which the q vectors should be generated.
    1769              * subset            -- a selection string specifying the atoms to consider for the analysis.
    1770              * deuteration       -- a selection string specifying the hydrogen atoms whose atomic parameters will be those of the deuterium.
    1771              * weights           -- a string equal to 'equal', 'mass', 'coherent' , 'incoherent' or 'atomicNumber' that specifies the weighting
    1772                                     scheme to use.
    1773              * eisf              -- the output NetCDF file name. A CDL version of this file will also be generated with the '.cdl' extension
    1774                                     instead of the '.nc' extension.
    1775              * pyroserver        -- a string specifying if Pyro will be used and how to run the analysis.
    1776
    1777      Running modes:
    1778
    1779          - To run the analysis do: a.runAnalysis() where a is the analysis object.
    1780          - To estimate the analysis do: a.estimateAnalysis() where a is the analysis object.
    1781          - To save the analysis to 'file' file name do: a.saveAnalysis(file) where a is the analysis object.
    1782
    1783      """
    1784
    1785      # The minimal set of input parameters names necessary to perform the analysis.
    1786      inputParametersNames = ('trajectory',\
    1787                              'timeinfo',\
    1788                              'qshellvalues',\
    1789                              'qshellwidth',\
    1790                              'qvectorspershell',\
    1791                              'qvectorsgenerator',\
    1792                              'qvectorsdirection',\
    1793                              'subset',\
    1794                              'deuteration',\
    1795                              'weights',\
    1796                              'eisf',\
    1797                              'pyroserver',\
    1798                              )
    1799
    1800      default = {'weights' : 'incoherent'}
    1801
    1802      shortName = 'EISF'
    1803
    1804      # Indicate whether this analysis can be estimated or not.
    1805      canBeEstimated = True
    1806
    1807 -    def __init__(self):
    1808          """The constructor. Insures that the class can not be instanciated directly from here.
    1809          """
    1810
    1811          raise Error('This class can not be instanciated.')
    1812
    1813 -    def initialize(self):
    1814          """Initializes the analysis (e.g. parses and checks input parameters, set some variables ...).
    1815          """
    1816
    1817          # The input parameters are parsed.
    1818          self.parseInputParameters()
    1819
    1820          if self.pyroServer != 'monoprocessor':
    1821              if self.statusBar is not None:
    1822                  LogMessage('warning', 'The job status bar will be desactivated for compatibility with Pyro module.')
    1823                  self.statusBar = None
    1824
    1825          self.buildTimeInfo()
    1826
    1827          self.subset = self.subsetSelection(self.universe, self.subset)
    1828          self.nAtoms = self.subset.numberOfAtoms()
    1829
    1830          self.deuteration = self.deuterationSelection(self.universe, self.deuteration)
    1831          self.deuteration = Collection(list(set(self.subset) & set(self.deuteration)))
    1832
    1833          self.weights = self.weightingScheme(self.universe, self.subset, self.deuteration, self.weights)
    1834
    1835          qv = QVectors(self.universe,\
    1836                        self.qVectorsGenerator,\
    1837                        self.qShellValues,\
    1838                        self.qShellWidth,\
    1839                        self.qVectorsPerShell,\
    1840                        self.qVectorsDirection)
    1841
    1842          self.qRadii = qv.qRadii
    1843          self.qVectors = qv.qVectors
    1844          self.qVectorsStatistics = qv.statistics
    1845
    1846          self.nQValues = len(self.qRadii)
    1847
    1848          self.preLoadTrajectory('atom')
    1849
    1850          # Dictionnary between the atom symbols and the index of the corresponding atoms in the selected atoms for the analysis.
    1851          self.symToInd = {}
    1852
    1853          # Dictionnary between the atom symbols and their corresponding weight.
    1854          self.symToWeight = {}
    1855
    1856          # Dictionnary between the atom index and their corresponding atom symbol.
    1857          self.indToSym = {}
    1858
    1859          for at in self.deuteration.atomList():
    1860              self.indToSym[at.index] = 'D'
    1861              if self.symToInd.has_key('D'):
    1862                  self.symToInd['D'].append(at.index)
    1863              else:
    1864                  self.symToInd['D'] = [at.index]
    1865                  self.symToWeight['D'] = self.weights[at]
    1866
    1867          for at in set(self.subset).difference(set(self.deuteration)):
    1868              self.indToSym[at.index] = at.symbol
    1869              if self.symToInd.has_key(at.symbol):
    1870                  self.symToInd[at.symbol].append(at.index)
    1871              else:
    1872                  self.symToInd[at.symbol] = [at.index]
    1873                  self.symToWeight[at.symbol] = self.weights[at]
    1874
    1875          self.EISF = {}
    1876          # Loop over all the atom symbol.
    1877          for symbol in self.symToInd.keys():
    1878              self.EISF[symbol] = N.zeros((self.nQValues,), typecode = N.Float)
    1879
    1880          for qIndex in range(self.nQValues):
    1881              qRadius = self.qRadii[qIndex]
    1882
    1883              # For each q vector length, the corresponding set of q vector is transform into an array
    1884              qArray = N.array(self.qVectors[qRadius])
    1885
    1886              # For each q vector length, the corresponding set of q vector is transform into an array
    1887              # qVect[1][qVal] is now of shape (3,Qcount)
    1888              self.qVectors[qRadius] = N.transpose(qArray)
    1889
    1890          if self.pyroServer != 'monoprocessor':
    1891              # The attribute trajectory is removed because it can not be pickled by Pyro.
    1892              delattr(self, 'trajectory')
    1893
    1894 -    def calc(self, atom, trajname):
    1895          """Calculates the atomic term.
    1896
    1897          @param atom: the atom on which the atomic term has been calculated.
    1898          @type atom: an instance of MMTK.Atom class.
    1899
    1900          @param trajname: the name of the trajectory file name.
    1901          @type trajname: string
    1902          """
    1903
    1904          if self.preLoadedTraj is None:
    1905              if self.pyroServer == 'monoprocessor':
    1906                  t = self.trajectory
    1907
    1908              else:
    1909                  # Load the whole trajectory set.
    1910                  t = Trajectory(None, trajname, 'r')
    1911
    1912              # series = 2D Numeric array. The positions of the selected atom |at| from the first step to the
    1913              # last step with the selected step increment.
    1914              series = t.readParticleTrajectory(atom, first = self.first, last = self.last, skip = self.skip).array
    1915
    1916          else:
    1917              series = self.preLoadedTraj[atom]
    1918
    1919          symbol = self.indToSym[atom.index]
    1920
    1921          atomicEISF = N.zeros((self.nQValues,), typecode = N.Float)
    1922
    1923          for qIndex in range(self.nQValues):
    1924              qRadius = self.qRadii[qIndex]
    1925
    1926              a       = N.add.reduce(N.exp(1j * N.dot(series, self.qVectors[qRadius])))/self.nFrames
    1927              a = (a*N.conjugate(a)).real
    1928              atomicEISF[qIndex] = N.add.reduce(a)/self.qVectors[qRadius].shape[1]
    1929
    1930          if self.preLoadedTraj is not None:
    1931              del self.preLoadedTraj[atom]
    1932
    1933          return atomicEISF
    1934
    1935 -    def combine(self, atom, x):
    1936          """
    1937          """
    1938          symbol = self.indToSym[atom.index]
    1939          N.add(self.EISF[symbol], x, self.EISF[symbol])
    1940
    1941 -    def finalize(self):
    1942          """Finalizes the calculations (e.g. averaging the total term, output files creations ...)
    1943          """
    1944
    1945          self.EISF['total'] = N.zeros((self.nQValues,), typecode = N.Float)
    1946          for symbol in self.symToInd.keys():
    1947              self.EISF[symbol] = self.EISF[symbol]/float(len(self.symToInd[symbol]))
    1948              N.add(self.EISF['total'],\
    1949                    float(len(self.symToInd[symbol]))*self.symToWeight[symbol]*self.EISF[symbol],\
    1950                    self.EISF['total'])
    1951
    1952          # The NetCDF output file is opened for writing.
    1953          outputFile       = NetCDFFile(self.outputEISF, 'w')
    1954          outputFile.title = self.__class__.__name__
    1955          outputFile.jobinfo = self.information + '\nOutput file written on: %s\n\n' % asctime()
    1956
    1957          # Some dimensions are created.
    1958          outputFile.createDimension('NQVALUES', self.nQValues)
    1959          outputFile.createDimension('NOCTANS', 8)
    1960          outputFile.createDimension('OCTANNAME', 8)
    1961
    1962          OCTANS = outputFile.createVariable('octan', N.Character, ('NOCTANS','OCTANNAME'))
    1963          OCTANS[:,:] = N.array([list('X+.Y+.Z+'),list('X+.Y+.Z-'),list('X+.Y-.Z+'),list('X+.Y-.Z-'),\
    1964                                  list('X-.Y+.Z+'),list('X-.Y+.Z-'),list('X-.Y-.Z+'),list('X-.Y-.Z-')])
    1965          OCTANS.units = 'unitless'
    1966
    1967          STAT = outputFile.createVariable('qvectors_statistics', N.Int, ('NQVALUES', 'NOCTANS'))
    1968          STAT[:,:] = self.qVectorsStatistics
    1969
    1970          # Creation of the NetCDF output variables.
    1971          # The Q.
    1972          QRADII = outputFile.createVariable('q', N.Float, ('NQVALUES',))
    1973          QRADII[:] = self.qRadii
    1974          QRADII.units = 'nm^-1'
    1975
    1976          for key, value in self.EISF.items():
    1977              EISF = outputFile.createVariable('eisf-%s' % ''.join(key), N.Float, ('NQVALUES',))
    1978              EISF[:] = value[:]
    1979              EISF.units = 'unitless'
    1980
    1981          asciiVar = sorted(outputFile.variables.keys())
    1982
    1983          outputFile.close()
    1984
    1985          self.toPlot = {'netcdf' : self.outputEISF, 'xVar' : 'q', 'yVar' : 'eisf-total'}
    1986
    1987          # Creates an ASCII version of the NetCDF output file.
    1988          convertNetCDFToASCII(inputFile = self.outputEISF,\
    1989                               outputFile = os.path.splitext(self.outputEISF)[0] + '.cdl',\
    1990                               variables = asciiVar)
    1991
    1992  #####################################################################################
    1993  # SMOOTHED STATIC COHERENT STRUCTURE FACTOR
    1994  #####################################################################################
    1995 -class SmoothedStaticCoherentStructureFactor(Analysis):
    1996      """Sets up an Smoothed Static Coherent Structure Factor.
    1997
    1998      A Subclass of nMOLDYN.Analysis.Analysis.
    1999
    2000      Constructor: SmoothedStaticCoherentStructureFactor(|parameters| = None)
    2001
    2002      Arguments:
    2003
    2004          - |parameters| -- a dictionnary of the input parameters, or 'None' to set up the analysis without parameters.
    2005              * trajectory        -- a trajectory file name or an instance of MMTK.Trajectory.Trajectory class.
    2006              * timeinfo          -- a string of the form 'first:last:step' where 'first' is an integer specifying the first frame
    2007                                     number to consider, 'last' is an integer specifying the last frame number to consider and
    2008                                     'step' is an integer specifying the step number between two frames.
    2009              * qshellvalues      -- a string of the form 'qmin1:qmax1:dq1;qmin2:qmax2:dq2...' where 'qmin1', 'qmin2' ... ,
    2010                                     'qmax1', 'qmax2' ... and 'dq1', 'dq2' ... are floats that represents respectively
    2011                                     the q minimum, the q maximum and the q steps for q interval 1, 2 ...
    2012              * subset            -- a selection string specifying the atoms to consider for the analysis.
    2013              * deuteration       -- a selection string specifying the hydrogen atoms whose atomic parameters will be those of the deuterium.
    2014              * weights           -- a string equal to 'equal', 'mass', 'coherent' , 'incoherent' or 'atomicNumber' that specifies the weighting
    2015                                     scheme to use.
    2016              * scsf              -- the output NetCDF file name. A CDL version of this file will also be generated with the '.cdl' extension
    2017                                     instead of the '.nc' extension.
    2018              * pyroserver        -- a string specifying if Pyro will be used and how to run the analysis.
    2019
    2020      Running modes:
    2021
    2022          - To run the analysis do: a.runAnalysis() where a is the analysis object.
    2023          - To estimate the analysis do: a.estimateAnalysis() where a is the analysis object.
    2024          - To save the analysis to 'file' file name do: a.saveAnalysis(file) where a is the analysis object.
    2025
    2026      Comments:
    2027
    2028          - The analysis is based on the angular averaged coherent static structure factor formula where the
    2029            summation over the q vectors is replaced by an integral over the q space. The formula used is taken
    2030            from equation 2.35 of Fischer et al. Rep. Prog. Phys. 69 (2006) 233-299.
    2031
    2032      """
    2033
    2034      # The minimal set of input parameters names necessary to perform the analysis.
    2035      inputParametersNames = ('trajectory',\
    2036                              'timeinfo',\
    2037                              'qshellvalues',\
    2038                              'subset',\
    2039                              'deuteration',\
    2040                              'weights',\
    2041                              'sscsf',\
    2042                              'pyroserver',\
    2043                              )
    2044
    2045      default = {'weights' : 'coherent'}
    2046
    2047      shortName = 'SSCSF'
    2048
    2049      # Indicate whether this analysis can be estimated or not.
    2050      canBeEstimated = True
    2051
    2052 -    def __init__(self):
    2053          """The constructor. Insures that the class can not be instanciated directly from here.
    2054          """
    2055          raise Error('This class can not be instanciated.')
    2056
    2057 -    def initialize(self):
    2058          """Initializes the analysis (e.g. parses and checks input parameters, set some variables ...).
    2059          """
    2060
    2061          # The input parameters are parsed.
    2062          self.parseInputParameters()
    2063
    2064          if self.universe.basisVectors() is None:
    2065              raise Error('The universe must be periodic for this kind of calculation.')
    2066
    2067          if self.pyroServer != 'monoprocessor':
    2068              if self.statusBar is not None:
    2069                  LogMessage('warning', 'The job status bar will be desactivated for compatibility with Pyro module.')
    2070                  self.statusBar = None
    2071
    2072          self.buildTimeInfo()
    2073
    2074          self.subset = self.subsetSelection(self.universe, self.subset)
    2075          self.nAtoms = self.subset.numberOfAtoms()
    2076
    2077          self.deuteration = self.deuterationSelection(self.universe, self.deuteration)
    2078          self.deuteration = Collection(list(set(self.subset) & set(self.deuteration)))
    2079
    2080          self.weights = self.weightingScheme(self.universe, self.subset, self.deuteration, self.weights)
    2081
    2082          self.nQValues = len(self.qShellValues)
    2083
    2084          self.preLoadTrajectory('frame')
    2085
    2086          # The number of atoms for each specie found in the universe.
    2087          self.population = {}
    2088
    2089          # The weight corresponding to each specie found in the universe..
    2090          self.symToWeight = {}
    2091
    2092          for at in self.deuteration.atomList():
    2093              if not self.symToWeight.has_key('D'):
    2094                  self.symToWeight['D'] = self.weights[at]
    2095
    2096              if self.population.has_key('D'):
    2097                  self.population['D'] += 1.0
    2098
    2099              else:
    2100                  self.population['D'] = 1.0
    2101
    2102          for at in set(self.subset).difference(set(self.deuteration)):
    2103              if not self.symToWeight.has_key(at.symbol):
    2104                  self.symToWeight[at.symbol] = self.weights[at]
    2105
    2106              if self.population.has_key(at.symbol):
    2107                  self.population[at.symbol] += 1.0
    2108
    2109              else:
    2110                  self.population[at.symbol] = 1.0
    2111
    2112          # The list of the species name found in the universe.
    2113          self.speciesList = sorted(self.population.keys())
    2114
    2115          # The number of species found in the universe.
    2116          self.nSpecies = len(self.speciesList)
    2117
    2118          # The indexes of the selected atoms.
    2119          self.indexes = N.zeros((self.nAtoms,), typecode = N.Int)
    2120
    2121          # The id of the specie to which belong each selected atom.
    2122          self.species = N.zeros((self.nAtoms,), typecode = N.Int)
    2123
    2124          self.rIndex = {}
    2125          comp = 0
    2126          for at in self.subset:
    2127              self.rIndex[at.index] = comp
    2128              comp += 1
    2129
    2130          for at in self.deuteration.atomList():
    2131              ind = self.rIndex[at.index]
    2132              self.species[ind] = self.speciesList.index('D')
    2133              self.indexes[ind] = at.index
    2134
    2135          for at in set(self.subset).difference(set(self.deuteration)):
    2136              ind = self.rIndex[at.index]
    2137              self.species[ind] = self.speciesList.index(at.symbol)
    2138              self.indexes[ind] = at.index
    2139
    2140          self.ssfIJ = N.zeros((self.nSpecies,self.nSpecies,self.nQValues), typecode = N.Float)
    2141
    2142          if self.pyroServer != 'monoprocessor':
    2143              # The attribute trajectory is removed because it can not be pickled by Pyro.
    2144              delattr(self, 'trajectory')
    2145
    2146 -    def calc(self, frameIndex, trajname):
    2147          """Calculates the contribution for one frame.
    2148
    2149          @param frameIndex: the index of the frame in |self.frameIndexes| array.
    2150          @type frameIndex: integer.
    2151
    2152          @param trajname: the name of the trajectory file name.
    2153          @type trajname: string
    2154          """
    2155
    2156          frame = self.frameIndexes[frameIndex]
    2157
    2158          if self.preLoadedTraj is None:
    2159              if self.pyroServer == 'monoprocessor':
    2160                  t = self.trajectory
    2161
    2162              else:
    2163                  # Load the whole trajectory set.
    2164                  t = Trajectory(None, trajname, 'r')
    2165
    2166              conf = t.configuration[frame]
    2167              self.universe.setConfiguration(conf)
    2168
    2169          else:
    2170              # The configuration is updated to the current trajectory step.
    2171              self.universe.setConfiguration(self.preLoadedTraj[frameIndex])
    2172
    2173          directCell = N.ravel(N.array([v for v in self.universe.basisVectors()]))
    2174          reverseCell = N.ravel(N.transpose(N.array([v for v in self.universe.reciprocalBasisVectors()])))
    2175
    2176          ssfIJTemp = N.zeros((self.nSpecies,self.nSpecies,self.nQValues), typecode = N.Float)
    2177
    2178          computeSSF(self.universe.contiguousObjectConfiguration().array,\
    2179                     directCell,\
    2180                     reverseCell,\
    2181                     N.array(self.qShellValues),\
    2182                     self.indexes,\
    2183                     self.species,\
    2184                     ssfIJTemp)
    2185
    2186          ssfIJTemp = 2.0*ssfIJTemp
    2187
    2188          if self.preLoadedTraj is not None:
    2189              del self.preLoadedTraj[frameIndex]
    2190
    2191          return ssfIJTemp
    2192
    2193 -    def combine(self, frameIndex, x):
    2194          """
    2195          """
    2196          N.add(self.ssfIJ, x, self.ssfIJ)
    2197
    2198 -    def finalize(self):
    2199          """Finalizes the calculations (e.g. averaging the total term, output files creations ...).
    2200          """
    2201
    2202          # The SSF.
    2203          SSF = {'total' : N.zeros((self.nQValues), typecode = N.Float)}
    2204          for i in range(self.nSpecies):
    2205              a = self.speciesList[i]
    2206              # The number of atoms of specie a.
    2207              nA = self.population[a]
    2208              # The weight corresponding to specie a.
    2209              wA = self.symToWeight[a]
    2210
    2211              for j in range(i, self.nSpecies):
    2212                  b = self.speciesList[j]
    2213                  # The number of atoms of specie b.
    2214                  nB = self.population[b]
    2215                  # The weight corresponding to specie b.
    2216                  wB = self.symToWeight[b]
    2217
    2218                  pair = a + '-' + b
    2219                  nAB = N.sqrt(nA*nB)
    2220                  wAB = wA*wB
    2221                  deltaAB = 1.0
    2222
    2223                  if i == j:
    2224                      SSF[pair] = self.ssfIJ[i,j,:][:]
    2225
    2226                  else:
    2227                      deltaAB = 0.0
    2228                      SSF[pair] = self.ssfIJ[i,j,:][:] + self.ssfIJ[j,i,:][:]
    2229
    2230                  SSF[pair] = deltaAB + SSF[pair]/(float(self.nFrames)*nAB)
    2231
    2232                  N.add(SSF['total'],nAB*wAB*SSF[pair],SSF['total'])
    2233
    2234          # The NetCDF output file is opened for writing.
    2235          outputFile       = NetCDFFile(self.outputSCSF, 'w')
    2236          outputFile.title = self.__class__.__name__
    2237          outputFile.jobinfo = self.information + '\nOutput file written on: %s\n\n' % asctime()
    2238
    2239          # Some dimensions are created.
    2240          outputFile.createDimension('NQVALUES', self.nQValues)
    2241
    2242          # Creation of the NetCDF output variables.
    2243          # The Q.
    2244          QRADII = outputFile.createVariable('q', N.Float, ('NQVALUES',))
    2245          QRADII[:] = self.qShellValues
    2246          QRADII.units = 'nm^-1'
    2247
    2248          for key, value in SSF.items():
    2249              SCSF = outputFile.createVariable('Sq-%s' % ''.join(key), N.Float, ('NQVALUES',))
    2250              SCSF[:] = value[:]
    2251              SCSF.units = 'unitless'
    2252
    2253          asciiVar = sorted(outputFile.variables.keys())
    2254
    2255          outputFile.close()
    2256
    2257          self.toPlot = {'netcdf' : self.outputSCSF, 'xVar' : 'q', 'yVar' : 'ssf-total'}
    2258
    2259          # Creates an ASCII version of the NetCDF output file.
    2260          convertNetCDFToASCII(inputFile = self.outputSCSF,\
    2261                               outputFile = os.path.splitext(self.outputSCSF)[0] + '.cdl',\
    2262                               variables = asciiVar)
    2263