1 """Collections of classes for the determination of structure-related properties.
2
3 Classes:
4
5 * PairDistributionFunction : sets up a Pair Distribution Function Analysis.
6 * CoordinationNumber : sets up a Coordination Number Analysis.
7 * SpatialDensity : sets up a Spatial Density Analysis.
8 * ScrewFit : sets up a Screw Fit Analysis.
9 """
10
11
12 import copy
13 from numpy import linalg
14 import operator
15 import os
16 import sys
17 from time import asctime
18 from timeit import default_timer
19
20
21 from Scientific import N
22 from Scientific.Geometry import Quaternion, Transformation, Vector
23 from Scientific.Geometry.Objects3D import Line
24 from Scientific.IO.NetCDF import NetCDFFile
25
26
27 from MMTK import Atom
28 from MMTK.Collections import Collection
29 from MMTK import Units
30 from MMTK.NucleicAcids import NucleotideChain
31 from MMTK.Proteins import PeptideChain, Protein
32 from MMTK.Trajectory import Trajectory
33 from MMTK.Universe import InfiniteUniverse
34
35
36 from nMOLDYN.Analysis.Analysis import Analysis
37 from nMOLDYN_CN import coordinationnumber
38 from nMOLDYN.Core.Error import Error
39 from nMOLDYN.Core.IOFiles import convertNetCDFToASCII
40 from nMOLDYN.Core.Logger import LogMessage
41 from nMOLDYN.Core.Mathematics import changeBasis, correlation, sphericalCoordinates
42 from nMOLDYN_RDF import distancehistogram
43
44
45
46
48 """Sets up a Pair Distribution Function analysis.
49
50 A Subclass of nMOLDYN.Analysis.Analysis.
51
52 Constructor: PairDistributionFunction(|parameters| = None)
53
54 Arguments:
55
56 - |parameters| -- a dictionnary of the input parameters, or 'None' to set up the analysis without parameters.
57 * trajectory -- a trajectory file name or an instance of MMTK.Trajectory.Trajectory class.
58 * timeinfo -- a string of the form 'first:last:step' where 'first' is an integer specifying the first frame
59 number to consider, 'last' is an integer specifying the last frame number to consider and
60 'step' is an integer specifying the step number between two frames.
61 * rvalues -- a string of the form 'rmin:rmax:dr' where 'rmin' is a float specifying the minimum distance to
62 consider, 'rmax' is a float specifying the maximum distance value to consider and 'dr' is a float
63 specifying the distance increment.
64 * subset -- a selection string specifying the atoms to consider for the analysis.
65 * deuteration -- a selection string specifying the hydrogen atoms whose atomic parameters will be those of the deuterium.
66 * weights -- a string equal to 'equal', 'mass', 'coherent' , 'incoherent' or 'atomicNumber' that specifies the weighting
67 scheme to use.
68 * pdf -- the output NetCDF file name. A CDL version of this file will also be generated with the '.cdl' extension
69 instead of the '.nc' extension.
70 * pyroserver -- a string specifying if Pyro will be used and how to run the analysis.
71
72 Running modes:
73
74 - To run the analysis do: a.runAnalysis() where a is the analysis object.
75 - To estimate the analysis do: a.estimateAnalysis() where a is the analysis object.
76 - To save the analysis to 'file' file name do: a.saveAnalysis(file) where a is the analysis object.
77
78 Comments:
79
80 - This code contains a pyrex function for the distance histogram calculation that is based on a FORTRAN code
81 written by Miguel Gonzalez, Insitut Laue Langevin, Grenoble, France.
82 """
83
84
85 inputParametersNames = ('trajectory',\
86 'timeinfo',\
87 'rvalues',\
88 'subset',\
89 'deuteration',\
90 'weights',\
91 'pdf',\
92 'pyroserver',\
93 )
94
95 shortName = 'PDF'
96
97
98 canBeEstimated = True
99
101 """The constructor. Insures that the class can not be instanciated directly from here.
102 """
103 raise Error('This class can not be instanciated.')
104
106 """Initializes the analysis (e.g. parses and checks input parameters, set some variables ...).
107 """
108
109
110 self.parseInputParameters()
111
112
113 if self.universe.basisVectors() is None:
114 raise Error('The universe must be periodic for this kind of calculation.')
115
116 if self.pyroServer != 'monoprocessor':
117 if self.statusBar is not None:
118 LogMessage('warning', 'The job status bar will be desactivated for compatibility with Pyro module.')
119 self.statusBar = None
120
121 self.buildTimeInfo()
122
123 self.subset = self.subsetSelection(self.universe, self.subset)
124 self.nAtoms = self.subset.numberOfAtoms()
125
126 self.deuteration = self.deuterationSelection(self.universe, self.deuteration)
127 self.deuteration = Collection(list(set(self.subset) & set(self.deuteration)))
128
129 self.weights = self.weightingScheme(self.universe, self.subset, self.deuteration, self.weights)
130
131 self.preLoadTrajectory('frame')
132
133
134 self.population = {}
135
136
137 self.symToWeight = {}
138
139 for at in self.deuteration.atomList():
140 if not self.symToWeight.has_key('D'):
141 self.symToWeight['D'] = self.weights[at]
142
143 if self.population.has_key('D'):
144 self.population['D'] += 1.0
145
146 else:
147 self.population['D'] = 1.0
148
149 for at in set(self.subset).difference(set(self.deuteration)):
150 if not self.symToWeight.has_key(at.symbol):
151 self.symToWeight[at.symbol] = self.weights[at]
152
153 if self.population.has_key(at.symbol):
154 self.population[at.symbol] += 1.0
155
156 else:
157 self.population[at.symbol] = 1.0
158
159
160 self.speciesList = sorted(self.population.keys())
161
162
163 self.nSpecies = len(self.speciesList)
164
165
166 self.indexes = N.zeros((self.nAtoms,), typecode = N.Int)
167
168
169 self.species = N.zeros((self.nAtoms,), typecode = N.Int)
170
171 self.molIdList = []
172
173 self.molecules = N.zeros((self.nAtoms,), typecode = N.Int)
174
175 self.rIndex = {}
176 comp = 0
177 for at in self.subset:
178 self.rIndex[at.index] = comp
179 comp += 1
180
181 for at in self.deuteration.atomList():
182 ind = self.rIndex[at.index]
183 self.species[ind] = self.speciesList.index('D')
184 self.indexes[ind] = at.index
185 molId = id(at.topLevelChemicalObject())
186 if molId not in self.molIdList:
187 self.molIdList.append(molId)
188 self.molecules[ind] = self.molIdList.index(molId)
189
190 for at in set(self.subset).difference(set(self.deuteration)):
191 ind = self.rIndex[at.index]
192 self.species[ind] = self.speciesList.index(at.symbol)
193 self.indexes[ind] = at.index
194 molId = id(at.topLevelChemicalObject())
195 if molId not in self.molIdList:
196 self.molIdList.append(molId)
197 self.molecules[ind] = self.molIdList.index(molId)
198
199
200 self.hIntra = N.zeros((self.nSpecies,self.nSpecies,self.nBins), typecode = N.Float)
201
202 self.hInter = N.zeros((self.nSpecies,self.nSpecies,self.nBins), typecode = N.Float)
203
204 self.averageDensity = 0.0
205
206 if self.pyroServer != 'monoprocessor':
207
208 delattr(self, 'trajectory')
209
210 - def calc(self, frameIndex, trajname):
211 """Calculates the contribution for one frame.
212
213 @param frameIndex: the index of the frame in |self.frameIndexes| array.
214 @type frameIndex: integer.
215
216 @param trajname: the name of the trajectory file name.
217 @type trajname: string
218 """
219
220 frame = self.frameIndexes[frameIndex]
221
222 if self.preLoadedTraj is None:
223 if self.pyroServer == 'monoprocessor':
224 t = self.trajectory
225
226 else:
227
228 t = Trajectory(None, trajname, 'r')
229
230 conf = t.configuration[frame]
231 self.universe.setConfiguration(conf)
232
233 else:
234
235 self.universe.setConfiguration(self.preLoadedTraj[frameIndex])
236
237 directCell = N.ravel(N.array([v for v in self.universe.basisVectors()]))
238 reverseCell = N.ravel(N.transpose(N.array([v for v in self.universe.reciprocalBasisVectors()])))
239 cellVolume = self.universe.cellVolume()
240
241 hIntraTemp = N.zeros(self.hIntra.shape, typecode = N.Float)
242 hInterTemp = N.zeros(self.hInter.shape, typecode = N.Float)
243
244 distancehistogram(self.universe.contiguousObjectConfiguration().array,\
245 directCell,\
246 reverseCell,\
247 self.indexes,\
248 self.molecules,\
249 self.species,\
250 hIntraTemp,\
251 hInterTemp,\
252 self.rMin,\
253 self.rMax,\
254 self.dr)
255
256 hIntraTemp = cellVolume * hIntraTemp
257 hInterTemp = cellVolume * hInterTemp
258
259 if self.preLoadedTraj is not None:
260 del self.preLoadedTraj[frameIndex]
261
262 return cellVolume, hIntraTemp, hInterTemp
263
265 """
266 """
267
268 self.averageDensity += self.nAtoms/x[0]
269
270
271
272
273 N.add(self.hIntra, x[1], self.hIntra)
274 N.add(self.hInter, x[2], self.hInter)
275
277 """Finalizes the calculations (e.g. averaging the total term, output files creations ...).
278 """
279
280 self.averageDensity /= self.nFrames
281 densityFactor = 4.0*N.pi*self.rValues
282 shellSurfaces = densityFactor*self.rValues
283 shellVolumes = shellSurfaces*self.dr
284
285
286 tcfIntra = {}
287
288 tcfInter = {}
289
290 tcfTotal = {'total' : N.zeros((self.nBins), typecode = N.Float)}
291
292
293 rdfIntra = {}
294
295 rdfInter = {}
296
297 rdfTotal = {'total' : N.zeros((self.nBins), typecode = N.Float)}
298
299
300 pdfIntra = {}
301
302 pdfInter = {}
303
304 pdfTotal = {'total' : N.zeros((self.nBins), typecode = N.Float)}
305 for i in range(self.nSpecies):
306 a = self.speciesList[i]
307
308 nA = self.population[a]
309
310 wA = self.symToWeight[a]
311
312 for j in range(i, self.nSpecies):
313 b = self.speciesList[j]
314
315 nB = self.population[b]
316
317 wB = self.symToWeight[b]
318
319 pair = a + '-' + b
320
321 if i == j:
322 pdfIntra[pair] = self.hIntra[i,j,:][:]
323 pdfInter[pair] = self.hInter[i,j,:][:]
324 wei = wA*wA*nA*nA
325 nAnB = nA*(nA-1)/2.0
326
327 else:
328 pdfIntra[pair] = self.hIntra[i,j,:][:] + self.hIntra[j,i,:][:]
329 pdfInter[pair] = self.hInter[i,j,:][:] + self.hInter[j,i,:][:]
330 wei = 2.0*wA*wB*nA*nB
331 nAnB = nA*nB
332
333 pdfIntra[pair] /= (float(self.nFrames)*shellVolumes*nAnB)
334 pdfInter[pair] /= (float(self.nFrames)*shellVolumes*nAnB)
335 pdfTotal[pair] = pdfIntra[pair] + pdfInter[pair]
336
337 rdfIntra[pair] = shellSurfaces*self.averageDensity*pdfIntra[pair]
338 rdfInter[pair] = shellSurfaces*self.averageDensity*pdfInter[pair]
339 rdfTotal[pair] = rdfIntra[pair] + rdfInter[pair]
340
341 tcfIntra[pair] = densityFactor*self.averageDensity*(pdfIntra[pair] - 1.0)
342 tcfInter[pair] = densityFactor*self.averageDensity*(pdfInter[pair] - 1.0)
343 tcfTotal[pair] = tcfIntra[pair] + tcfInter[pair]
344
345 N.add(rdfTotal['total'],wei*rdfTotal[pair],rdfTotal['total'])
346 N.add(pdfTotal['total'],wei*pdfTotal[pair],pdfTotal['total'])
347 N.add(tcfTotal['total'],wei*tcfTotal[pair],tcfTotal['total'])
348
349
350 outputFile = NetCDFFile(self.outputPDF, 'w')
351 outputFile.title = self.__class__.__name__
352 outputFile.jobinfo = self.information + '\nOutput file written on: %s\n\n' % asctime()
353
354
355 outputFile.createDimension('NRVALUES', int(self.nBins))
356
357
358
359 RVALUES = outputFile.createVariable('r', N.Float, ('NRVALUES',))
360 RVALUES[:] = self.rValues
361 RVALUES.units = 'nm'
362
363 for k in pdfIntra.keys():
364
365 PDFINTRA = outputFile.createVariable('pdf-%s-intra' % k, N.Float, ('NRVALUES',))
366 PDFINTRA[:] = pdfIntra[k]
367 PDFINTRA.units = 'nm^3'
368
369 PDFINTER = outputFile.createVariable('pdf-%s-inter' % k, N.Float, ('NRVALUES',))
370 PDFINTER[:] = pdfInter[k]
371 PDFINTER.units = 'nm^3'
372
373 PDFTOTAL = outputFile.createVariable('pdf-%s' % k, N.Float, ('NRVALUES',))
374 PDFTOTAL[:] = pdfTotal[k]
375 PDFTOTAL.units = 'nm^3'
376
377 RDFINTRA = outputFile.createVariable('rdf-%s-intra' % k, N.Float, ('NRVALUES',))
378 RDFINTRA[:] = rdfIntra[k]
379 RDFINTRA.units = 'unitless'
380
381 RDFINTER = outputFile.createVariable('rdf-%s-inter' % k, N.Float, ('NRVALUES',))
382 RDFINTER[:] = rdfInter[k]
383 RDFINTER.units = 'unitless'
384
385 RDFTOTAL = outputFile.createVariable('rdf-%s' % k, N.Float, ('NRVALUES',))
386 RDFTOTAL[:] = rdfTotal[k]
387 RDFTOTAL.units = 'unitless'
388
389 TCFINTRA = outputFile.createVariable('tcf-%s-intra' % k, N.Float, ('NRVALUES',))
390 TCFINTRA[:] = tcfIntra[k]
391 TCFINTRA.units = 'unitless'
392
393 TCFINTER = outputFile.createVariable('tcf-%s-inter' % k, N.Float, ('NRVALUES',))
394 TCFINTER[:] = tcfInter[k]
395 TCFINTER.units = 'unitless'
396
397 TCFTOTAL = outputFile.createVariable('tcf-%s' % k, N.Float, ('NRVALUES',))
398 TCFTOTAL[:] = tcfTotal[k]
399 TCFTOTAL.units = 'unitless'
400
401 PDF = outputFile.createVariable('pdf-total', N.Float, ('NRVALUES',))
402 PDF[:] = pdfTotal['total']
403 PDF.units = 'nm^3'
404
405 RDF = outputFile.createVariable('rdf-total', N.Float, ('NRVALUES',))
406 RDF[:] = rdfTotal['total']
407 RDF.units = 'unitless'
408
409 TCF = outputFile.createVariable('tcf-total', N.Float, ('NRVALUES',))
410 TCF[:] = tcfTotal['total']
411 TCF.units = 'unitless'
412
413 asciiVar = sorted(outputFile.variables.keys())
414
415 outputFile.close()
416
417 self.toPlot = {'netcdf' : self.outputPDF, 'xVar' : 'r', 'yVar' : 'pdf-total'}
418
419
420 convertNetCDFToASCII(inputFile = self.outputPDF,\
421 outputFile = os.path.splitext(self.outputPDF)[0] + '.cdl',\
422 variables = asciiVar)
423
424
425
426
428 """Sets up a Coordination Number analysis.
429
430 A Subclass of nMOLDYN.Analysis.Analysis.
431
432 Constructor: CoordinationNumber(|parameters| = None)
433
434 Arguments:
435
436 - |parameters| -- a dictionnary of the input parameters, or 'None' to set up the analysis without parameters.
437 * trajectory -- a trajectory file name or an instance of MMTK.Trajectory.Trajectory class.
438 * timeinfo -- a string of the form 'first:last:step' where 'first' is an integer specifying the first frame
439 number to consider, 'last' is an integer specifying the last frame number to consider and
440 'step' is an integer specifying the step number between two frames.
441 * rvalues -- a string of the form 'rmin:rmax:dr' where 'rmin' is a float specifying the minimum distance to
442 consider, 'rmax' is a float specifying the maximum distance value to consider and 'dr' is a float
443 specifying the distance increment.
444 * group -- a selection string specifying the groups of atoms that will be used to define the points around which
445 the coordination number will be computed. For each group, there is one point defined as the center of
446 gravity of the group.
447 * subset -- a selection string specifying the atoms to consider for the analysis.
448 * deuteration -- a selection string specifying the hydrogen atoms whose atomic parameters will be those of the deuterium.
449 * cn -- the output NetCDF file name. A CDL version of this file will also be generated with the '.cdl' extension
450 instead of the '.nc' extension.
451 * pyroserver -- a string specifying if Pyro will be used and how to run the analysis.
452
453 Running modes:
454
455 - To run the analysis do: a.runAnalysis() where a is the analysis object.
456 - To estimate the analysis do: a.estimateAnalysis() where a is the analysis object.
457 - To save the analysis to 'file' file name do: a.saveAnalysis(file) where a is the analysis object.
458
459 Comments:
460
461 - This code contains a pyrex function for the distance histogram calculation than enhances significantly its
462 performance.
463 """
464
465
466 inputParametersNames = ('trajectory',\
467 'timeinfo',\
468 'rvalues',\
469 'group',\
470 'subset',\
471 'deuteration',\
472 'cn',\
473 'pyroserver',\
474 )
475
476 shortName = 'CN'
477
478
479 canBeEstimated = True
480
482 """The constructor. Insures that the class can not be instanciated directly from here.
483 """
484 raise Error('This class can not be instanciated.')
485
487 """Initializes the analysis (e.g. parses and checks input parameters, set some variables ...).
488 """
489
490
491 self.parseInputParameters()
492
493 if self.universe.basisVectors() is None:
494 raise Error('The universe must be periodic for this kind of calculation.')
495
496 if self.pyroServer != 'monoprocessor':
497 if self.statusBar is not None:
498 LogMessage('warning', 'The job status bar will be desactivated for compatibility with Pyro module.')
499 self.statusBar = None
500
501 self.buildTimeInfo()
502
503 self.subset = self.subsetSelection(self.universe, self.subset)
504 self.nAtoms = self.subset.numberOfAtoms()
505
506 self.deuteration = self.deuterationSelection(self.universe, self.deuteration)
507 self.deuteration = Collection(list(set(self.subset) & set(self.deuteration)))
508
509 self.group = self.groupSelection(self.universe, self.group)
510 self.nGroups = len(self.group)
511
512 self.preLoadTrajectory('frame')
513
514
515 self.population = {}
516
517 for at in self.deuteration.atomList():
518 if self.population.has_key('D'):
519 self.population['D'] += 1.0
520
521 else:
522 self.population['D'] = 1.0
523
524 for at in set(self.subset).difference(set(self.deuteration)):
525 if self.population.has_key(at.symbol):
526 self.population[at.symbol] += 1.0
527
528 else:
529 self.population[at.symbol] = 1.0
530
531
532 self.speciesList = sorted(self.population.keys())
533
534
535 self.nSpecies = len(self.speciesList)
536
537
538 self.indexes = N.zeros((self.nAtoms,), typecode = N.Int)
539
540
541 self.species = N.zeros((self.nAtoms,), typecode = N.Int)
542
543 self.molIdList = []
544
545 self.molecules = N.zeros((self.nAtoms,), typecode = N.Int)
546
547 self.rIndex = {}
548 comp = 0
549 for at in self.subset:
550 self.rIndex[at.index] = comp
551 comp += 1
552
553 for at in self.deuteration.atomList():
554 ind = self.rIndex[at.index]
555 self.species[ind] = self.speciesList.index('D')
556 self.indexes[ind] = at.index
557 molId = id(at.topLevelChemicalObject())
558 if molId not in self.molIdList:
559 self.molIdList.append(molId)
560 self.molecules[ind] = self.molIdList.index(molId)
561
562 for at in set(self.subset).difference(set(self.deuteration)):
563 ind = self.rIndex[at.index]
564 self.species[ind] = self.speciesList.index(at.symbol)
565 self.indexes[ind] = at.index
566 molId = id(at.topLevelChemicalObject())
567 if molId not in self.molIdList:
568 self.molIdList.append(molId)
569 self.molecules[ind] = self.molIdList.index(molId)
570
571 self.groupIndexes = []
572 self.gMolId = []
573 for g in self.group:
574 firstAtom = g.atomList()[0]
575 self.gMolId.append(self.molIdList.index(id(firstAtom.topLevelChemicalObject())))
576 self.groupIndexes.append(N.array([at.index for at in g.atomList()]))
577
578
579 self.hIntra = N.zeros((self.nSpecies,self.nBins,self.nFrames), typecode = N.Float)
580
581 self.hInter = N.zeros((self.nSpecies,self.nBins,self.nFrames), typecode = N.Float)
582
583 if self.pyroServer != 'monoprocessor':
584
585 delattr(self, 'trajectory')
586
587 - def calc(self, frameIndex, trajname):
588 """Calculates the contribution for one frame.
589
590 @param frameIndex: the index of the frame in |self.frameIndexes| array.
591 @type frameIndex: integer.
592
593 @param trajname: the name of the trajectory file name.
594 @type trajname: string
595 """
596
597 frame = self.frameIndexes[frameIndex]
598
599 if self.preLoadedTraj is None:
600 if self.pyroServer == 'monoprocessor':
601 t = self.trajectory
602
603 else:
604
605 t = Trajectory(None, trajname, 'r')
606
607 conf = t.configuration[frame]
608 self.universe.setConfiguration(conf)
609
610 else:
611
612 self.universe.setConfiguration(self.preLoadedTraj[frameIndex])
613
614 directCell = N.ravel(N.array([v for v in self.universe.basisVectors()]))
615 reverseCell = N.ravel(N.transpose(N.array([v for v in self.universe.reciprocalBasisVectors()])))
616
617 hIntraTemp = N.zeros((self.nSpecies,self.nBins), typecode = N.Float)
618 hInterTemp = N.zeros((self.nSpecies,self.nBins), typecode = N.Float)
619 for gInd in range(len(self.group)):
620 coordinationnumber(self.universe.contiguousObjectConfiguration().array,\
621 directCell,\
622 reverseCell,\
623 self.groupIndexes[gInd],\
624 self.gMolId[gInd],\
625 self.indexes,\
626 self.molecules,\
627 self.species,\
628 hIntraTemp,\
629 hInterTemp,\
630 self.rMin,\
631 self.rMax,\
632 self.dr)
633
634 if self.preLoadedTraj is not None:
635 del self.preLoadedTraj[frameIndex]
636
637 return hIntraTemp, hInterTemp
638
640 """
641 """
642
643 frame = self.frameIndexes[frameIndex]
644 index = list(self.frameIndexes).index(frame)
645
646 N.add(self.hIntra[:,:,index], x[0], self.hIntra[:,:,index])
647 N.add(self.hInter[:,:,index], x[1], self.hInter[:,:,index])
648
650 """Finalizes the calculations (e.g. averaging the total term, output files creations ...).
651 """
652
653 coordNumberIntra = {}
654 coordNumberInter = {}
655 coordNumberTotal = {'total' : N.zeros((self.nBins, self.nFrames), typecode = N.Float)}
656 for i in range(self.nSpecies):
657 a = self.speciesList[i]
658
659 coordNumberIntra[a] = self.hIntra[i,:,:]/float(self.nGroups)
660 coordNumberInter[a] = self.hInter[i,:,:]/float(self.nGroups)
661
662 coordNumberTotal[a] = coordNumberIntra[a] + coordNumberInter[a]
663
664 N.add(coordNumberTotal['total'], coordNumberTotal[a], coordNumberTotal['total'])
665
666
667 outputFile = NetCDFFile(self.outputCN, 'w')
668 outputFile.title = self.__class__.__name__
669 outputFile.jobinfo = self.information + '\nOutput file written on: %s\n\n' % asctime()
670
671
672 outputFile.createDimension('NRVALUES', int(self.nBins))
673 outputFile.createDimension('NTIMES', self.nFrames)
674
675
676
677 RVALUES = outputFile.createVariable('r', N.Float, ('NRVALUES',))
678 RVALUES[:] = self.rValues
679 RVALUES.units = 'nm'
680
681 TIMES = outputFile.createVariable('time', N.Float, ('NTIMES',))
682 TIMES[:] = self.times[:]
683 TIMES.units = 'ps'
684
685 for k in coordNumberIntra.keys():
686
687 CNINTRA = outputFile.createVariable('cn-%s-intra' % k, N.Float, ('NRVALUES','NTIMES'))
688 CNINTRA[:] = coordNumberIntra[k]
689 CNINTRA.units = 'unitless'
690
691 CNINTER = outputFile.createVariable('cn-%s-inter' % k, N.Float, ('NRVALUES','NTIMES'))
692 CNINTER[:] = coordNumberInter[k]
693 CNINTER.units = 'unitless'
694
695 CNTOTAL = outputFile.createVariable('cn-%s' % k, N.Float, ('NRVALUES','NTIMES'))
696 CNTOTAL[:] = coordNumberTotal[k]
697 CNTOTAL.units = 'unitless'
698
699 CNCUMULINTRA = outputFile.createVariable('cn-cumul-%s-intra' % k, N.Float, ('NRVALUES','NTIMES'))
700 CNCUMULINTRA[:] = N.cumsum(coordNumberIntra[k])
701 CNCUMULINTRA.units = 'unitless'
702
703 CNCUMULINTER = outputFile.createVariable('cn-cumul-%s-inter' % k, N.Float, ('NRVALUES','NTIMES'))
704 CNCUMULINTER[:] = N.cumsum(coordNumberInter[k])
705 CNCUMULINTER.units = 'unitless'
706
707 CNCUMULTOTAL = outputFile.createVariable('cn-cumul-%s' % k, N.Float, ('NRVALUES','NTIMES'))
708 CNCUMULTOTAL[:] = N.cumsum(coordNumberTotal[k])
709 CNCUMULTOTAL.units = 'unitless'
710
711 CNAVGINTRA = outputFile.createVariable('cn-timeavg-%s-intra' % k, N.Float, ('NRVALUES',))
712 CNAVGINTRA[:] = coordNumberIntra[k].sum(1)/float(self.nFrames)
713 CNAVGINTRA.units = 'unitless'
714
715 CNAVGINTER = outputFile.createVariable('cn-timeavg-%s-inter' % k, N.Float, ('NRVALUES',))
716 CNAVGINTER[:] = coordNumberInter[k].sum(1)/float(self.nFrames)
717 CNAVGINTER.units = 'unitless'
718
719 CNAVGTOTAL = outputFile.createVariable('cn-timeavg-%s' % k, N.Float, ('NRVALUES',))
720 CNAVGTOTAL[:] = coordNumberTotal[k].sum(1)/float(self.nFrames)
721 CNAVGTOTAL.units = 'unitless'
722
723 CNAVGCUMULINTRA = outputFile.createVariable('cn-timeavg-cumul-%s-intra' % k, N.Float, ('NRVALUES',))
724 CNAVGCUMULINTRA[:] = N.cumsum(coordNumberIntra[k].sum(1))/float(self.nFrames)
725 CNAVGCUMULINTRA.units = 'unitless'
726
727 CNAVGCUMULINTER = outputFile.createVariable('cn-timeavg-cumul-%s-inter' % k, N.Float, ('NRVALUES',))
728 CNAVGCUMULINTER[:] = N.cumsum(coordNumberInter[k].sum(1))/float(self.nFrames)
729 CNAVGCUMULINTER.units = 'unitless'
730
731 CNAVGCUMULTOTAL = outputFile.createVariable('cn-timeavg-cumul-%s' % k, N.Float, ('NRVALUES',))
732 CNAVGCUMULTOTAL[:] = N.cumsum(coordNumberTotal[k].sum(1))/float(self.nFrames)
733 CNAVGCUMULTOTAL.units = 'unitless'
734
735 CNAVG = outputFile.createVariable('cn-timeavg', N.Float, ('NRVALUES',))
736 CNAVG[:] = coordNumberTotal['total'].sum(1)/float(self.nFrames)
737 CNAVG.units = 'unitless'
738
739 CNAVGCUMUL = outputFile.createVariable('cn-timeavg-cumul', N.Float, ('NRVALUES',))
740 CNAVGCUMUL[:] = N.cumsum(coordNumberTotal['total'].sum(1))/float(self.nFrames)
741 CNAVGCUMUL.units = 'unitless'
742
743 asciiVar = sorted(outputFile.variables.keys())
744
745 outputFile.close()
746
747 self.toPlot = {'netcdf' : self.outputCN, 'xVar' : 'r', 'yVar' : 'cn-total-tavg-cumul'}
748
749
750 convertNetCDFToASCII(inputFile = self.outputCN,\
751 outputFile = os.path.splitext(self.outputCN)[0] + '.cdl',\
752 variables = asciiVar)
753
754
755
756
758 """Set up a Screw Fit analysis.
759
760 A Subclass of nMOLDYN.Analysis.Analysis.
761
762 Constructor: ScrewFit(|parameters| = None)
763
764 Arguments:
765
766 - |parameters| -- a dictionnary of the input parameters, or 'None' to set up the analysis without parameters.
767 * trajectory -- a trajectory file name or an instance of MMTK.Trajectory.Trajectory class.
768 * timeinfo -- a string of the form 'first:last:step' where 'first' is an integer specifying the first frame
769 number to consider, 'last' is an integer specifying the last frame number to consider and
770 'step' is an integer specifying the step number between two frames.
771 * sfa -- the output NetCDF file name. A CDL version of this file will also be generated with the '.cdl' extension
772 instead of the '.nc' extension.
773 * pyroserver -- a string specifying if Pyro will be used and how to run the analysis.
774
775 Running modes:
776
777 - To run the analysis do: a.runAnalysis() where a is the analysis object.
778 - To estimate the analysis do: a.estimateAnalysis() where a is the analysis object.
779 - To save the analysis to 'file' file name do: a.saveAnalysis(file) where a is the analysis object.
780
781 Comments:
782
783 - This code is based on a first implementation made by Paolo Calligari.
784
785 - For more details: Kneller, G.R., Calligari, P. Acta Crystallographica , D62, 302-311
786 """
787
788
789 inputParametersNames = ('trajectory',\
790 'timeinfo',\
791 'sfa',\
792 'pyroserver',\
793 )
794
795 shortName = 'SFA'
796
797
798 canBeEstimated = True
799
801 """The constructor. Insures that the class can not be instanciated directly from here.
802 """
803 raise Error('This class can not be instanciated.')
804
806 """Initializes the analysis (e.g. parses and checks input parameters, set some variables ...).
807 """
808
809
810 self.parseInputParameters()
811
812 if self.pyroServer != 'monoprocessor':
813 if self.statusBar is not None:
814 LogMessage('warning', 'The job status bar will be desactivated for compatibility with Pyro module.')
815 self.statusBar = None
816
817
818 self.buildTimeInfo()
819
820 self.preLoadTrajectory('frame')
821
822
823
824 self.orientDist = {}
825 self.helixRadius = {}
826 self.straightness = {}
827 foundChain = False
828 for obj in self.universe.objectList():
829 if isinstance(obj, (PeptideChain, Protein)):
830 for chain in obj:
831 foundChain = True
832 self.orientDist[chain] = N.zeros((self.nFrames,len(chain) - 1), typecode = N.Float)
833 self.helixRadius[chain] = N.zeros((self.nFrames,len(chain) - 3), typecode = N.Float)
834 self.straightness[chain] = N.zeros((self.nFrames,len(chain) - 4), typecode = N.Float)
835
836
837 if not foundChain:
838 raise Error('The universe must contains at least one peptide chain to perform a screw fit analysis.')
839
840 if self.pyroServer != 'monoprocessor':
841
842 delattr(self, 'trajectory')
843
844 - def calc(self, frameIndex, trajname):
845 """Calculates the contribution for one frame.
846
847 @param frameIndex: the index of the frame in |self.frameIndexes| array.
848 @type frameIndex: integer.
849
850 @param trajname: the name of the trajectory file name.
851 @type trajname: string
852 """
853
854 frame = self.frameIndexes[frameIndex]
855
856 if self.preLoadedTraj is None:
857 if self.pyroServer == 'monoprocessor':
858 t = self.trajectory
859
860 else:
861
862 t = Trajectory(None, trajname, 'r')
863
864 conf = t.configuration[frame]
865 self.universe.setConfiguration(conf)
866
867 else:
868
869 self.universe.setConfiguration(self.preLoadedTraj[frameIndex])
870
871 od = {}
872 hr = {}
873 s = {}
874 for obj in self.universe.objectList():
875 if isinstance(obj, (PeptideChain, Protein)):
876 for chain in obj:
877 od[chain] = N.zeros((len(chain) - 1), typecode = N.Float)
878 hr[chain] = N.zeros((len(chain) - 3), typecode = N.Float)
879 s[chain] = N.zeros((len(chain) - 4), typecode = N.Float)
880
881 for obj in self.universe.objectList():
882 if isinstance(obj, (PeptideChain, Protein)):
883 for chain in obj:
884 od[chain][:] = self.angularDistance(chain)
885 hr[chain][:], s[chain][:] = self.screwMotionAnalysis(chain)
886
887 if self.preLoadedTraj is not None:
888 del self.preLoadedTraj[frameIndex]
889
890 return od, hr, s
891
893 """
894 """
895
896 for obj in self.universe.objectList():
897 if isinstance(obj, (PeptideChain, Protein)):
898 for chain in obj:
899 self.orientDist[chain][frameIndex,:] = x[0][chain][:]
900 self.helixRadius[chain][frameIndex,:], self.straightness[chain][frameIndex,:] = x[1][chain][:], x[2][chain][:]
901
903 """Finalizes the calculations (e.g. averaging the total term, output files creations ...).
904 """
905
906 outputFile = NetCDFFile(self.outputSFA, 'w', self.information)
907 outputFile.title = self.__class__.__name__
908 outputFile.jobinfo = self.information + '\nOutput file written on: %s\n\n' % asctime()
909
910 outputFile.createDimension('NTIMES', self.nFrames)
911
912 TIMES = outputFile.createVariable('time', N.Float, ('NTIMES',))
913 TIMES[:] = self.times
914 TIMES.units = 'ps'
915
916 for chain in self.straightness.keys():
917
918 sequence = [r.sequence_number for r in chain.residues()]
919 name = chain.fullName().upper()
920
921 straightnessY = name + '_STRAIGHTNESS_Y'
922 helixRadiusY = name + '_HELIXRADIUS_Y'
923 orientDistY = name + '_ORIENTDIST_Y'
924
925 outputFile.createDimension(straightnessY, len(chain)-4)
926 outputFile.createDimension(helixRadiusY, len(chain)-3)
927 outputFile.createDimension(orientDistY, len(chain)-1)
928
929 STRAIGHTNESS_Y = outputFile.createVariable(straightnessY.lower(), N.Int, (straightnessY,))
930 HELIXRADIUS_Y = outputFile.createVariable(helixRadiusY.lower(), N.Int, (helixRadiusY,))
931 ORIENTDIST_Y = outputFile.createVariable(orientDistY.lower(), N.Int, (orientDistY,))
932
933 for r in range(len(chain)-4):
934 res = sequence[r]
935 STRAIGHTNESS_Y[r] = res
936 HELIXRADIUS_Y[r] = res
937 ORIENTDIST_Y[r] = res
938
939 res = sequence[r+1]
940 HELIXRADIUS_Y[r+1] = res
941 ORIENTDIST_Y[r+1] = res
942
943 res = sequence[r+2]
944 ORIENTDIST_Y[r+2] = res
945
946 res = sequence[r+3]
947 ORIENTDIST_Y[r+3] = res
948
949 STRAIGHTNESS = outputFile.createVariable(name + '_straightness', N.Float, ('NTIMES',straightnessY))
950 STRAIGHTNESS[:,:] = self.straightness[chain][:,:]
951
952 HELIXRADIUS = outputFile.createVariable(name + '_helixradius' , N.Float, ('NTIMES',helixRadiusY))
953 HELIXRADIUS.units = 'nm'
954 HELIXRADIUS[:,:] = self.helixRadius[chain][:,:]
955
956 ORIENTDIST = outputFile.createVariable(name + '_orientdist' , N.Float, ('NTIMES',orientDistY))
957 ORIENTDIST[:,:] = self.orientDist[chain][:,:]
958
959 asciiVar = sorted(outputFile.variables.keys())
960
961 outputFile.close()
962
963 self.toPlot = None
964
965
966 convertNetCDFToASCII(inputFile = self.outputSFA,\
967 outputFile = os.path.splitext(self.outputSFA)[0] + '.cdl',\
968 variables = asciiVar)
969
971 """
972 Returns the complete matrix of quaternions compatibles with linear trasformation.|conf1| is
973 the reference configuration. |point_ref| is the reference point about which the fit is calculated
974 """
975
976 if conf1.universe != peptide.universe():
977 raise Error("conformation is for a different universe")
978
979 if conf2 is None:
980 conf1, conf2 = conf2, conf1
981
982 else:
983 if conf2.universe != peptide.universe():
984 raise Error('conformation is for a different universe')
985
986 ref = conf1
987 conf = conf2
988 weights = peptide.universe().masses()
989 weights = weights/peptide.mass()
990 ref_cms = point_ref.position().array
991 pos = N.zeros((3,), typecode = N.Float)
992 pos = point_ref.position(conf).array
993 possq = 0.
994 cross = N.zeros((3, 3), typecode = N.Float)
995
996 for a in peptide.atomList():
997 r = a.position(conf).array - pos
998 r_ref = a.position(ref).array-ref_cms
999 w = weights[a]
1000 possq = possq + w*N.add.reduce(r*r) + w*N.add.reduce(r_ref*r_ref)
1001 cross = cross + w*r[:, N.NewAxis]*r_ref[N.NewAxis, :]
1002
1003 k = N.zeros((4, 4), typecode = N.Float)
1004 k[0, 0] = -cross[0, 0] - cross[1, 1] - cross[2, 2]
1005 k[0, 1] = cross[1, 2] - cross[2, 1]
1006 k[0, 2] = cross[2, 0] - cross[0, 2]
1007 k[0, 3] = cross[0, 1] - cross[1, 0]
1008 k[1, 1] = -cross[0, 0] + cross[1, 1] + cross[2, 2]
1009 k[1, 2] = -cross[0, 1] - cross[1, 0]
1010 k[1, 3] = -cross[0, 2] - cross[2, 0]
1011 k[2, 2] = cross[0, 0] - cross[1, 1] + cross[2, 2]
1012 k[2, 3] = -cross[1, 2] - cross[2, 1]
1013 k[3, 3] = cross[0, 0] + cross[1, 1] - cross[2, 2]
1014
1015 for i in range(1, 4):
1016 for j in range(i):
1017 k[i, j] = k[j, i]
1018
1019 k = 2.*k
1020 for i in range(4):
1021 k[i, i] = k[i, i] + possq - N.add.reduce(pos*pos)
1022
1023 e, v = linalg.eigh(k)
1024 emin = e.argmin()
1025
1026 v = v[emin]
1027 if v[0] < 0: v = -v
1028 if e[emin] <= 0.:
1029 rms = 0.
1030 else:
1031 rms = N.sqrt(e[emin])
1032
1033 if matrix:
1034 emax = e.argmax()
1035 QuatMatrix = v
1036
1037 return Quaternion.Quaternion(QuatMatrix),v, e, e[emin],e[emax], rms
1038
1039 else:
1040 return Quaternion.Quaternion(v), Vector(ref_cms), Vector(pos), rms
1041
1047
1049
1050 universe = InfiniteUniverse()
1051
1052 Crefpos = Vector(0.00,0.00,0.00)
1053 Orefpos = Vector(0.00,0.00,0.123)
1054 Nrefpos = Vector(0.111,0.00,-0.0728)
1055
1056 Catom = Atom('C', position = Crefpos)
1057 Oatom = Atom('O', position = Orefpos)
1058 Natom = Atom('N', position = Nrefpos)
1059
1060 plane = Collection()
1061 plane.addObject(Catom)
1062 plane.addObject(Oatom)
1063 plane.addObject(Natom)
1064
1065 universe.addObject(plane)
1066
1067 refconf = copy.copy(universe.configuration())
1068
1069 chainLength = len(chain) - 1
1070 orientDist = N.zeros((chainLength,),typecode = N.Float)
1071
1072 for resInd in range(chainLength):
1073
1074 Cpos = chain[resInd].peptide.C.position()
1075 Opos = chain[resInd].peptide.O.position()
1076 Npos = chain[resInd+1].peptide.N.position()
1077
1078 Catom.setPosition(Cpos)
1079 Oatom.setPosition(Opos)
1080 Natom.setPosition(Npos)
1081
1082 distV = universe.distanceVector(Cpos,Crefpos)
1083 Catom.translateBy(distV)
1084 Oatom.translateBy(distV)
1085 Natom.translateBy(distV)
1086
1087 Qm, v, e, emin, emax, rms = self.findQuaternionMatrix(plane, Catom, refconf, matrix = True)
1088
1089 rmsD = plane.rmsDifference(refconf)
1090
1091 orientDist[resInd] = rmsD/ N.sqrt(emax)
1092
1093 refconf = copy.copy(universe.configuration())
1094
1095 return orientDist
1096
1098
1099 Crefpos = chain[0].peptide.C.position()
1100 Orefpos = chain[0].peptide.O.position()
1101 Nrefpos = chain[1].peptide.N.position()
1102
1103 Catom = Atom('C', position=Crefpos)
1104 Oatom = Atom('O', position=Orefpos)
1105 Natom = Atom('N', position=Nrefpos)
1106
1107 plane = Collection()
1108 plane.addObject(Catom)
1109 plane.addObject(Oatom)
1110 plane.addObject(Natom)
1111
1112 universe = InfiniteUniverse()
1113 universe.addObject(plane)
1114
1115 refconf = copy.copy(universe.configuration())
1116
1117 chainLength = len(chain)
1118
1119 DatiAsse=[]
1120 for i in range(chainLength - 2):
1121
1122 Cpos = chain[i+1].peptide.C.position()
1123 Opos = chain[i+1].peptide.O.position()
1124 Npos = chain[i+2].peptide.N.position()
1125
1126 Catom.setPosition(Cpos)
1127 Oatom.setPosition(Opos)
1128 Natom.setPosition(Npos)
1129
1130 newconf = copy.copy(universe.configuration())
1131
1132 TL1, ROT, TL2, rms = self.findGenericTransformation(plane, Catom, refconf, newconf)
1133 Tr = TL1*ROT*TL2
1134 DatiAsse.append(Tr.screwMotion())
1135 refconf = copy.copy(universe.configuration())
1136
1137 asseIP = []
1138 for i in range(len(DatiAsse)):
1139 if DatiAsse[i][1].length() != 0.:
1140 asseIP.append(Line(DatiAsse[i][0],DatiAsse[i][1]))
1141 else:
1142 asseIP.append(Line(DatiAsse[i][0],DatiAsse[i+1][1]))
1143
1144 helixRadius = N.zeros((chainLength-3,), typecode = N.Float)
1145
1146 straightness = N.zeros((chainLength-4,), typecode = N.Float)
1147
1148 for resInd in range(chainLength-3):
1149
1150 Ca = chain[resInd].peptide.C.position()
1151 Ca1 = chain[resInd+1].peptide.C.position()
1152 Ca2 = chain[resInd+2].peptide.C.position()
1153
1154 helixRadius[resInd]= asseIP[resInd].distanceFrom(Ca)
1155
1156 if resInd <= len(DatiAsse)-3:
1157 projCa = asseIP[resInd].projectionOf(Ca)
1158 projCa1 = asseIP[resInd+1].projectionOf(Ca1)
1159 projCa2 = asseIP[resInd+2].projectionOf(Ca2)
1160
1161 pscal = (projCa - projCa1)*(projCa1 - projCa2)
1162 mod1 = (projCa - projCa1).length()
1163 mod2 = (projCa1 - projCa2).length()
1164
1165 straightness[resInd] = pscal/(mod1*mod2)
1166
1167 return helixRadius, straightness
1168
1169
1170
1171
1173 """Sets up a Spatial Density analysis.
1174
1175 A Subclass of nMOLDYN.Analysis.Analysis.
1176
1177 Constructor: SpatialDensity(|parameters| = None)
1178
1179 Arguments:
1180
1181 - |parameters| -- a dictionnary of the input parameters, or 'None' to set up the analysis without parameters.
1182 * trajectory -- a trajectory file name or an instance of MMTK.Trajectory.Trajectory class.
1183 * timeinfo -- a string of the form 'first:last:step' where 'first' is an integer specifying the first frame
1184 number to consider, 'last' is an integer specifying the last frame number to consider and
1185 'step' is an integer specifying the step number between two frames.
1186 * rvalues -- a string of the form 'rmin:rmax:dr' where 'rmin' is a float specifying the minimum distance to
1187 consider, 'rmax' is a float specifying the maximum distance value to consider and 'dr' is a float
1188 specifying the distance increment.
1189 * group -- a selection string specifying the groups of atoms that will be used to define the points around which
1190 the coordination number will be computed. For each group, there is one point defined as the center of
1191 gravity of the group.
1192 * atomorder -- a string of the form 'atom1,atom2,atom3' where 'atom1', 'atom2' and 'atom3' are
1193 respectively the MMTK atom names of the atoms in the way they should be ordered.
1194 * target -- a selection string specifying the groups of atoms that will be used to define the points around which
1195 the coordination number will be computed. For each group, there is one point defined as the center of
1196 gravity of the group.
1197 * sd -- the output NetCDF file name. A CDL version of this file will also be generated with the '.cdl' extension
1198 instead of the '.nc' extension.
1199 * pyroserver -- a string specifying if Pyro will be used and how to run the analysis.
1200
1201 Running modes:
1202
1203 - To run the analysis do: a.runAnalysis() where a is the analysis object.
1204 - To estimate the analysis do: a.estimateAnalysis() where a is the analysis object.
1205 - To save the analysis to 'file' file name do: a.saveAnalysis(file) where a is the analysis object.
1206
1207 Comments:
1208
1209 - This code contains a pyrex function for the distance histogram calculation than enhances significantly its
1210 performance.
1211 """
1212
1213
1214 inputParametersNames = ('trajectory',\
1215 'timeinfo',\
1216 'rvalues',\
1217 'thetavalues',\
1218 'phivalues',\
1219 'triplet',\
1220 'atomorder',\
1221 'group',\
1222 'sd',\
1223 'pyroserver',\
1224 )
1225
1226 shortName = 'SD'
1227
1228
1229 canBeEstimated = True
1230
1232 """The constructor. Insures that the class can not be instanciated directly from here.
1233 """
1234 raise Error('This class can not be instanciated.')
1235
1237 """Initializes the analysis (e.g. parses and checks input parameters, set some variables ...).
1238 """
1239
1240
1241 self.parseInputParameters()
1242
1243 if self.universe.basisVectors() is None:
1244 raise Error('The universe must be periodic for this kind of calculation.')
1245
1246 if self.pyroServer != 'monoprocessor':
1247 if self.statusBar is not None:
1248 LogMessage('warning', 'The job status bar will be desactivated for compatibility with Pyro module.')
1249 self.statusBar = None
1250
1251 self.buildTimeInfo()
1252
1253 self.triplet = self.groupSelection(self.universe, self.triplet)
1254 self.triplet = [g for g in self.triplet if g.numberOfAtoms() == 3]
1255
1256 if self.atomOrder is None:
1257 self.triplet = [sorted(g, key = operator.attrgetter('name')) for g in self.triplet]
1258
1259 else:
1260 orderedTriplets = []
1261 for triplet in self.triplet:
1262 tr = []
1263 for atName in self.atomOrder:
1264 found = False
1265 for at in triplet:
1266 if at.name == atName:
1267 found = True
1268 tr.append(at)
1269 if len(tr) == 3:
1270 orderedTriplets.append(tr)
1271 tr = []
1272 else:
1273 if not found:
1274 raise Error('No atom corresponding %s could be found.' % atName)
1275
1276 self.triplet = orderedTriplets
1277
1278 if len(self.triplet) == 0:
1279 raise Error('Your group selection does not contain any triplet.')
1280
1281 self.nTriplets = len(self.triplet)
1282
1283 self.group = self.groupSelection(self.universe, self.group)
1284
1285 self.nGroups = len(self.group)
1286
1287 self.preLoadTrajectory('frame')
1288
1289 self.SD = N.zeros((self.nBins,self.nThetaBins,self.nPhiBins), N.Float)
1290
1291
1292
1293 if self.pyroServer != 'monoprocessor':
1294
1295 delattr(self, 'trajectory')
1296
1297 - def calc(self, frameIndex, trajname):
1298 """Calculates the contribution for one frame.
1299
1300 @param frameIndex: the index of the frame in |self.frameIndexes| array.
1301 @type frameIndex: integer.
1302
1303 @param trajname: the name of the trajectory file name.
1304 @type trajname: string
1305 """
1306
1307 frame = self.frameIndexes[frameIndex]
1308
1309 if self.preLoadedTraj is None:
1310 if self.pyroServer == 'monoprocessor':
1311 t = self.trajectory
1312
1313 else:
1314
1315 t = Trajectory(None, trajname, 'r')
1316
1317 conf = t.configuration[frame]
1318 self.universe.setConfiguration(self.universe.contiguousObjectConfiguration(Collection(self.group), conf))
1319
1320 else:
1321
1322 self.universe.setConfiguration(self.universe.contiguousObjectConfiguration(Collection(self.group), self.preLoadedTraj[frameIndex]))
1323
1324 comList = []
1325 for g in self.group:
1326 comList.append(g.centerOfMass())
1327
1328 sphericalBinList = []
1329
1330 for t in self.triplet:
1331
1332 tcom = Collection(t).centerOfMass()
1333 vi, vj, vk = self.constructBasisFromAtoms(t)
1334
1335 for com in comList:
1336
1337 if com == tcom:
1338 continue
1339
1340 r = (tcom - com).length()
1341
1342 if (r < self.rMin) or (r > self.rMax):
1343 continue
1344
1345 x, y, z = changeBasis(tcom,com,vi,vj,vk)
1346
1347
1348
1349
1350
1351
1352
1353
1354
1355
1356
1357 theta = N.arccos(z/r)
1358
1359 if (theta < self.thetaMin) or (theta > self.thetaMax):
1360 continue
1361
1362
1363 phi = N.arctan2(y,x)
1364
1365 if (phi < self.phiMin) or (phi > self.phiMax):
1366 continue
1367
1368 rBin = int((r - self.rMin)/self.dr)
1369 thetaBin = int((theta - self.thetaMin)/self.dTheta)
1370 phiBin = int((phi - self.phiMin)/self.dPhi)
1371
1372 sphericalBinList.append((rBin, thetaBin, phiBin))
1373
1374 if self.preLoadedTraj is not None:
1375 del self.preLoadedTraj[frameIndex]
1376
1377 return sphericalBinList
1378
1379 - def combine(self, frameIndex, x):
1380 """
1381 """
1382 for rBin, thetaBin, phiBin in x:
1383
1384 self.SD[rBin, thetaBin, phiBin] += 1.0
1385
1387 """Finalizes the calculations (e.g. averaging the total term, output files creations ...).
1388 """
1389
1390
1391
1392
1393
1394
1395 self.SD /= float(self.nFrames*self.nTriplets*self.nGroups)
1396
1397
1398 outputFile = NetCDFFile(self.outputSD, 'w')
1399 outputFile.title = self.__class__.__name__
1400 outputFile.jobinfo = self.information + '\nOutput file written on: %s\n\n' % asctime()
1401
1402
1403 outputFile.createDimension('NTHETAVALUES', self.nThetaBins)
1404 outputFile.createDimension('NPHIVALUES', self.nPhiBins)
1405
1406
1407 THETAS = outputFile.createVariable('theta', N.Float, ('NTHETAVALUES',))
1408 THETAS[:] = self.thetaValues
1409 THETAS.units = 'deg'
1410
1411 PHIS = outputFile.createVariable('phi', N.Float, ('NPHIVALUES',))
1412 PHIS[:] = self.phiValues
1413 PHIS.units = 'deg'
1414
1415 for rIndex in range(len(self.rValues)):
1416
1417 r = self.rValues[rIndex]
1418
1419 SD = outputFile.createVariable('sd_r%snm' % round(r,3), N.Float, ('NTHETAVALUES','NPHIVALUES'))
1420 SD[:,:] = self.SD[rIndex,:,:]
1421 SD.units = 'unitless'
1422
1423 asciiVar = sorted(outputFile.variables.keys())
1424
1425 outputFile.close()
1426
1427 self.toPlot = None
1428
1429
1430 convertNetCDFToASCII(inputFile = self.outputSD,\
1431 outputFile = os.path.splitext(self.outputSD)[0] + '.cdl',\
1432 variables = asciiVar)
1433
1435 """This method construct a set of three oriented orthonormal axes i, j, k from a triplet of atoms
1436 such as (i,j,k) forms a clockwise orthonormal basis.
1437 If a1, a2 and a3 stand respectively for the three atoms of the triplet then:
1438 vector1 = (vector(a1,a2)_normalized + vector(a1,a3)_normalized)_normalized
1439 vector3 = (vector1 ^ vector(a1,a3))_normalized and correclty oriented
1440 vector2 = (vector3 ^ vector1)_normalized
1441
1442 @param triplet: the triplet of atoms.
1443 @type triplet: a list of three MMTK Atoms
1444
1445 @return: the three axis.
1446 @rtype: a list of three Scientific Vector
1447
1448 """
1449
1450 a1, a2, a3 = triplet
1451 v1 = self.universe.distanceVector(a1.position(), a2.position()).normal()
1452 v2 = self.universe.distanceVector(a1.position(), a3.position()).normal()
1453
1454 i = (v1 + v2).normal()
1455 k = i.cross(v2).normal()
1456 if k*a1.position() < 0:
1457 k = -k
1458 j = k.cross(i).normal()
1459
1460 return i, j, k
1461