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