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
17 import os
18 import re
19 from time import asctime
20 from timeit import default_timer
21
22
23 from Scientific import N
24 from Scientific.IO.NetCDF import _NetCDFFile, NetCDFFile
25 from Scientific.Signals.Models import AutoRegressiveModel, AveragedAutoRegressiveModel
26
27
28 from MMTK.Collections import Collection
29 from MMTK_forcefield import NonbondedList
30 from MMTK.Trajectory import Trajectory
31
32
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
41
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
61 tim = SF.variables['time']
62 dt = tim[1] - tim[0]
63
64 nQValues = len(SF.variables['q'][:])
65
66
67 frequencies = N.arange(len(tim))/(2.0*len(tim)*dt)
68
69
70 SF.createDimension('NFREQUENCIES', len(frequencies[:]))
71
72
73
74 FREQ = SF.createVariable('frequency', N.Float, ('NFREQUENCIES',))
75 FREQ.units = 'THz'
76
77
78 FREQ[:] = frequencies[:]
79
80
81 partialTerms = [re.findall('Fqt-(.*)',v)[0] for v in SF.variables if 'Fqt-' in v]
82 for pTerm in partialTerms:
83
84 DSF = SF.createVariable('Sqw-'+pTerm, N.Float, ('NQVALUES','NFREQUENCIES'))
85
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
94 convertNetCDFToASCII(inputFile = netcdf,\
95 outputFile = os.path.splitext(netcdf)[0] + '.cdl',\
96 variables = asciiVar)
97
98
99
100
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
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
163 canBeEstimated = False
164
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
171 """Initializes the analysis (e.g. parses and checks input parameters, set some variables ...).
172 """
173
174
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
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
228 for symbol1 in self.atomicSubsets.keys():
229
230 for symbol2 in self.atomicSubsets.keys():
231
232 pairName = tuple(sorted((symbol1,symbol2)))
233
234 self.FQT[pairName] = N.zeros((self.nQValues,self.nFrames), typecode = N.Float)
235
236
237 for qIndex in range(self.nQValues):
238
239 qRadius = self.qRadii[qIndex]
240
241
242 qarray = N.array(self.qVectors[qRadius])
243 qarray = self.universe._boxToRealPointArray(qarray)
244
245 self.qVectors[qRadius] = N.transpose(qarray)
246
247
248
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
256 for frameIndex in range(self.nFrames):
257
258 frame = self.frameIndexes[frameIndex]
259
260
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
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
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
301 for frameIndex in range(self.nFrames):
302
303 frame = self.frameIndexes[frameIndex]
304
305 if self.allConf is None:
306
307
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
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
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
342
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
348 outputFile = NetCDFFile(self.outputDCSF, 'w')
349
350 outputFile.title = self.__class__.__name__
351 outputFile.jobinfo = self.information + '\nOutput file written on: %s\n\n' % asctime()
352
353
354 outputFile.createDimension('NQVALUES', self.nQValues)
355 outputFile.createDimension('NTIMES', self.nFrames)
356 outputFile.createDimension('NOCTANS', 8)
357 outputFile.createDimension('OCTANNAME', 8)
358
359
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
392
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
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
454 canBeEstimated = False
455
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
462 """Initializes the analysis (e.g. parses and checks input parameters, set some variables ...).
463 """
464
465
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
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
519 for symbol1 in self.atomicSubsets.keys():
520
521 for symbol2 in self.atomicSubsets.keys():
522
523 pairName = tuple(sorted((symbol1,symbol2)))
524
525 self.SQ[pairName] = N.zeros((self.nQValues,), typecode = N.Float)
526
527
528 for qIndex in range(self.nQValues):
529
530 qRadius = self.qRadii[qIndex]
531
532
533 qarray = N.array(self.qVectors[qRadius])
534 qarray = self.universe._boxToRealPointArray(qarray)
535
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
544 for frameIndex in range(self.nFrames):
545
546 frame = self.frameIndexes[frameIndex]
547
548
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
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
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
589 for frameIndex in range(self.nFrames):
590
591 frame = self.frameIndexes[frameIndex]
592
593 if self.allConf is None:
594
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
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
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
639 outputFile = NetCDFFile(self.outputSCSF, 'w')
640
641 outputFile.title = self.__class__.__name__
642 outputFile.jobinfo = self.information + '\nOutput file written on: %s\n\n' % asctime()
643
644
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
675
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
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
735 canBeEstimated = False
736
737 default = {'weights' : 'coherent'}
738
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
745 """Initializes the analysis (e.g. parses and checks input parameters, set some variables ...).
746 """
747
748
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
783 for qIndex in range(self.nQValues):
784
785 qRadius = self.qRadii[qIndex]
786
787
788 qarray = N.array(self.qVectors[qRadius])
789 qarray = self.universe._boxToRealPointArray(qarray)
790
791 self.qVectors[qRadius] = N.transpose(qarray)
792
793
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
808 for frameIndex in range(self.nFrames):
809
810 frame = self.frameIndexes[frameIndex]
811
812
813 conf = self.trajectory.configuration[frame]
814 conf.convertToBoxCoordinates()
815
816
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
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
842 t = Trajectory(None, trajname, 'r')
843
844 qRadius = self.qRadii[qIndex]
845
846 model = AveragedAutoRegressiveModel(self.arModelOrder, self.dt)
847
848
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
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
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
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
960
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
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
1023 canBeEstimated = True
1024
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
1032 """Initializes the analysis (e.g. parses and checks input parameters, set some variables ...).
1033 """
1034
1035
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
1069 self.symToInd = {}
1070
1071
1072 self.symToWeight = {}
1073
1074
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
1096 for symbol in self.symToInd.keys():
1097
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
1103 qArray = N.array(self.qVectors[qRadius])
1104 self.qVectors[qRadius] = N.transpose(qArray)
1105
1106 if self.pyroServer != 'monoprocessor':
1107
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
1126 t = Trajectory(None, trajname, 'r')
1127
1128
1129
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
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
1155 """
1156 """
1157 symbol = self.indToSym[atom.index]
1158 N.add(self.ISF[symbol], x, self.ISF[symbol])
1159
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
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
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
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
1213
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
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
1273 canBeEstimated = True
1274
1275 default = {'weights' : 'incoherent'}
1276
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
1283 """Initializes the analysis (e.g. parses and checks input parameters, set some variables ...).
1284 """
1285
1286
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
1321 for qIndex in range(self.nQValues):
1322
1323 qRadius = self.qRadii[qIndex]
1324
1325
1326 qarray = N.array(self.qVectors[qRadius])
1327
1328 self.qVectors[qRadius] = N.transpose(qarray)
1329
1330
1331 self.ARModel = {}
1332
1333 for qIndex in range(self.nQValues):
1334 qRadius = self.qRadii[qIndex]
1335 self.ARModel[qRadius] = AveragedAutoRegressiveModel(self.arModelOrder, self.dt)
1336
1337
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
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
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
1401 """
1402 """
1403
1404 N.add(self.FQT, x[0], self.FQT)
1405 N.add(self.SQW, x[1], self.SQW)
1406
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
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
1490
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
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
1542 canBeEstimated = True
1543
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
1551 """Initializes the analysis (e.g. parses and checks input parameters, set some variables ...).
1552 """
1553
1554
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
1574
1575 self.qSquare = N.array(self.qShellValues)*N.array(self.qShellValues)
1576
1577 self.preLoadTrajectory('atom')
1578
1579
1580 self.symToInd = {}
1581
1582
1583 self.symToWeight = {}
1584
1585
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
1608 for symbol in self.symToInd.keys():
1609
1610 self.ISFG[symbol] = N.zeros((self.nQValues,self.nFrames), typecode = N.Float)
1611
1612 if self.pyroServer != 'monoprocessor':
1613
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
1632 t = Trajectory(None, trajname, 'r')
1633
1634
1635
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
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
1659 """
1660 """
1661 symbol = self.indToSym[atom.index]
1662 N.add(self.ISFG[symbol], x, self.ISFG[symbol])
1663
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
1676 outputFile = NetCDFFile(self.outputDISFG, 'w')
1677
1678 outputFile.title = self.__class__.__name__
1679 outputFile.jobinfo = self.information + '\nOutput file written on: %s\n\n' % asctime()
1680
1681
1682 outputFile.createDimension('NQVALUES', self.nQValues)
1683 outputFile.createDimension('NTIMES', self.nFrames)
1684
1685
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
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
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
1717
1718 dsq = N.add.reduce(series * series,1)
1719
1720
1721 sum_dsq1 = N.add.accumulate(dsq)
1722
1723
1724 sum_dsq2 = N.add.accumulate(dsq[::-1])
1725
1726
1727 sumsq = 2.*sum_dsq1[-1]
1728
1729
1730
1731
1732 Saabb = sumsq - N.concatenate(([0.], sum_dsq1[:-1])) - N.concatenate(([0.], sum_dsq2[:-1]))
1733
1734
1735
1736 Saabb = Saabb / (len(dsq) - N.arange(len(dsq)))
1737 Sab = 2.*correlation(series)
1738
1739
1740 return (Saabb - Sab)
1741
1742
1743
1744
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
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
1805 canBeEstimated = True
1806
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
1814 """Initializes the analysis (e.g. parses and checks input parameters, set some variables ...).
1815 """
1816
1817
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
1851 self.symToInd = {}
1852
1853
1854 self.symToWeight = {}
1855
1856
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
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
1884 qArray = N.array(self.qVectors[qRadius])
1885
1886
1887
1888 self.qVectors[qRadius] = N.transpose(qArray)
1889
1890 if self.pyroServer != 'monoprocessor':
1891
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
1910 t = Trajectory(None, trajname, 'r')
1911
1912
1913
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
1936 """
1937 """
1938 symbol = self.indToSym[atom.index]
1939 N.add(self.EISF[symbol], x, self.EISF[symbol])
1940
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
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
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
1971
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
1988 convertNetCDFToASCII(inputFile = self.outputEISF,\
1989 outputFile = os.path.splitext(self.outputEISF)[0] + '.cdl',\
1990 variables = asciiVar)
1991
1992
1993
1994
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
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
2050 canBeEstimated = True
2051
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
2058 """Initializes the analysis (e.g. parses and checks input parameters, set some variables ...).
2059 """
2060
2061
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
2087 self.population = {}
2088
2089
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
2113 self.speciesList = sorted(self.population.keys())
2114
2115
2116 self.nSpecies = len(self.speciesList)
2117
2118
2119 self.indexes = N.zeros((self.nAtoms,), typecode = N.Int)
2120
2121
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
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
2164 t = Trajectory(None, trajname, 'r')
2165
2166 conf = t.configuration[frame]
2167 self.universe.setConfiguration(conf)
2168
2169 else:
2170
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
2199 """Finalizes the calculations (e.g. averaging the total term, output files creations ...).
2200 """
2201
2202
2203 SSF = {'total' : N.zeros((self.nQValues), typecode = N.Float)}
2204 for i in range(self.nSpecies):
2205 a = self.speciesList[i]
2206
2207 nA = self.population[a]
2208
2209 wA = self.symToWeight[a]
2210
2211 for j in range(i, self.nSpecies):
2212 b = self.speciesList[j]
2213
2214 nB = self.population[b]
2215
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
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
2240 outputFile.createDimension('NQVALUES', self.nQValues)
2241
2242
2243
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
2260 convertNetCDFToASCII(inputFile = self.outputSCSF,\
2261 outputFile = os.path.splitext(self.outputSCSF)[0] + '.cdl',\
2262 variables = asciiVar)
2263