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

Source Code for Module nMOLDYN.Analysis.Structure

   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  # The python distribution modules 
  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  # The ScientificPython modules 
  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  # The MMTK distribution modules 
  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  # The nMOLDYN modules 
  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  # PAIR DISTRIBUTION FUNCTION ANALYSIS 
  46  ##################################################################################### 
47 -class PairDistributionFunction(Analysis):
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 # The minimal set of input parameters names necessary to perform the analysis. 85 inputParametersNames = ('trajectory',\ 86 'timeinfo',\ 87 'rvalues',\ 88 'subset',\ 89 'deuteration',\ 90 'weights',\ 91 'pdf',\ 92 'pyroserver',\ 93 ) 94 95 shortName = 'PDF' 96 97 # Indicate whether this analysis can be estimated or not. 98 canBeEstimated = True 99
100 - def __init__(self):
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
105 - def initialize(self):
106 """Initializes the analysis (e.g. parses and checks input parameters, set some variables ...). 107 """ 108 109 # The input parameters are parsed. 110 self.parseInputParameters() 111 112 # Some additional checkings. 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 # The number of atoms for each specie found in the universe. 134 self.population = {} 135 136 # The weight corresponding to each specie found in the universe.. 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 # The list of the species name found in the universe. 160 self.speciesList = sorted(self.population.keys()) 161 162 # The number of species found in the universe. 163 self.nSpecies = len(self.speciesList) 164 165 # The indexes of the selected atoms. 166 self.indexes = N.zeros((self.nAtoms,), typecode = N.Int) 167 168 # The id of the specie to which belong each selected atom. 169 self.species = N.zeros((self.nAtoms,), typecode = N.Int) 170 171 self.molIdList = [] 172 # The id of the molecule to which belongs each selected atoms. 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 # The histogram of the intramolecular distances. 200 self.hIntra = N.zeros((self.nSpecies,self.nSpecies,self.nBins), typecode = N.Float) 201 # The histogram of the intermolecular distances. 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 # The attribute trajectory is removed because it can not be pickled by Pyro. 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 # Load the whole trajectory set. 228 t = Trajectory(None, trajname, 'r') 229 230 conf = t.configuration[frame] 231 self.universe.setConfiguration(conf) 232 233 else: 234 # The configuration is updated to the current trajectory step. 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
264 - def combine(self, frameIndex, x):
265 """ 266 """ 267 268 self.averageDensity += self.nAtoms/x[0] 269 270 # The temporary distance histograms are normalized by the volume. This is done for each step because the 271 # volume can variate during the MD (e.g. NPT conditions). This volume is the one that intervene in the density 272 # calculation. 273 N.add(self.hIntra, x[1], self.hIntra) 274 N.add(self.hInter, x[2], self.hInter)
275
276 - def finalize(self):
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 # The intramolecular TCF. 286 tcfIntra = {} 287 # The intermolecular TCF. 288 tcfInter = {} 289 # The total TCF. 290 tcfTotal = {'total' : N.zeros((self.nBins), typecode = N.Float)} 291 292 # The intramolecular RDF. 293 rdfIntra = {} 294 # The intermolecular RDF. 295 rdfInter = {} 296 # The total RDF. 297 rdfTotal = {'total' : N.zeros((self.nBins), typecode = N.Float)} 298 299 # The intramolecular RDF. 300 pdfIntra = {} 301 # The intermolecular RDF. 302 pdfInter = {} 303 # The total RDF. 304 pdfTotal = {'total' : N.zeros((self.nBins), typecode = N.Float)} 305 for i in range(self.nSpecies): 306 a = self.speciesList[i] 307 # The number of atoms of specie a. 308 nA = self.population[a] 309 # The weight corresponding to specie a. 310 wA = self.symToWeight[a] 311 312 for j in range(i, self.nSpecies): 313 b = self.speciesList[j] 314 # The number of atoms of specie b. 315 nB = self.population[b] 316 # The weight corresponding to specie b. 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 # The NetCDF output file is opened for writing. 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 # Some dimensions are created. 355 outputFile.createDimension('NRVALUES', int(self.nBins)) 356 357 # Creation of the NetCDF output variables. 358 # The Q. 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 # Creates an ASCII version of the NetCDF output file. 420 convertNetCDFToASCII(inputFile = self.outputPDF,\ 421 outputFile = os.path.splitext(self.outputPDF)[0] + '.cdl',\ 422 variables = asciiVar)
423 424 ##################################################################################### 425 # COORDINATION NUMBER ANALYSIS 426 #####################################################################################
427 -class CoordinationNumber(Analysis):
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 # The minimal set of input parameters names necessary to perform the analysis. 466 inputParametersNames = ('trajectory',\ 467 'timeinfo',\ 468 'rvalues',\ 469 'group',\ 470 'subset',\ 471 'deuteration',\ 472 'cn',\ 473 'pyroserver',\ 474 ) 475 476 shortName = 'CN' 477 478 # Indicate whether this analysis can be estimated or not. 479 canBeEstimated = True 480
481 - def __init__(self):
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
486 - def initialize(self):
487 """Initializes the analysis (e.g. parses and checks input parameters, set some variables ...). 488 """ 489 490 # The input parameters are parsed. 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 # The number of atoms for each specie found in the universe. 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 # The list of the species name found in the universe. 532 self.speciesList = sorted(self.population.keys()) 533 534 # The number of species found in the universe. 535 self.nSpecies = len(self.speciesList) 536 537 # The indexes of the selected atoms. 538 self.indexes = N.zeros((self.nAtoms,), typecode = N.Int) 539 540 # The id of the specie to which belong each selected atom. 541 self.species = N.zeros((self.nAtoms,), typecode = N.Int) 542 543 self.molIdList = [] 544 # The id of the molecule to which belongs each selected atoms. 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 # The histogram of the intramolecular distances. 579 self.hIntra = N.zeros((self.nSpecies,self.nBins,self.nFrames), typecode = N.Float) 580 # The histogram of the intermolecular distances. 581 self.hInter = N.zeros((self.nSpecies,self.nBins,self.nFrames), typecode = N.Float) 582 583 if self.pyroServer != 'monoprocessor': 584 # The attribute trajectory is removed because it can not be pickled by Pyro. 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 # Load the whole trajectory set. 605 t = Trajectory(None, trajname, 'r') 606 607 conf = t.configuration[frame] 608 self.universe.setConfiguration(conf) 609 610 else: 611 # The configuration is updated to the current trajectory step. 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
639 - def combine(self, frameIndex, x):
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
649 - def finalize(self):
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 # The NetCDF output file is opened for writing. 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 # Some dimensions are created. 672 outputFile.createDimension('NRVALUES', int(self.nBins)) 673 outputFile.createDimension('NTIMES', self.nFrames) 674 675 # Creation of the NetCDF output variables. 676 # The Q. 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 # Creates an ASCII version of the NetCDF output file. 750 convertNetCDFToASCII(inputFile = self.outputCN,\ 751 outputFile = os.path.splitext(self.outputCN)[0] + '.cdl',\ 752 variables = asciiVar)
753 754 ##################################################################################### 755 # SCREW FIT ANALYSIS 756 #####################################################################################
757 -class ScrewFitAnalysis(Analysis):
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 # The minimal set of input parameters names necessary to perform the analysis. 789 inputParametersNames = ('trajectory',\ 790 'timeinfo',\ 791 'sfa',\ 792 'pyroserver',\ 793 ) 794 795 shortName = 'SFA' 796 797 # Indicate whether this analysis can be estimated or not. 798 canBeEstimated = True 799
800 - def __init__(self):
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
805 - def initialize(self):
806 """Initializes the analysis (e.g. parses and checks input parameters, set some variables ...). 807 """ 808 809 # The input parameters are parsed. 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 # Some additional checkings. 818 self.buildTimeInfo() 819 820 self.preLoadTrajectory('frame') 821 822 # The dictionnaries that will store respectively the orientational distance, the helix radius and the 823 # straightness for each chain found in the universe. 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 # Case where no protein or peptide chain was found in the universe. 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 # The attribute trajectory is removed because it can not be pickled by Pyro. 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 # Load the whole trajectory set. 862 t = Trajectory(None, trajname, 'r') 863 864 conf = t.configuration[frame] 865 self.universe.setConfiguration(conf) 866 867 else: 868 # The configuration is updated to the current trajectory step. 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
892 - def combine(self, frameIndex, x):
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
902 - def finalize(self):
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 # Creates an ASCII version of the NetCDF output file. 966 convertNetCDFToASCII(inputFile = self.outputSFA,\ 967 outputFile = os.path.splitext(self.outputSFA)[0] + '.cdl',\ 968 variables = asciiVar)
969
970 - def findQuaternionMatrix(self, peptide, point_ref, conf1, conf2 = None, matrix = True):
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
1042 - def findGenericTransformation(self,peptide, point_ref,conf1, conf2 = None):
1043 1044 q, cm1, cm2, rms = self.findQuaternionMatrix(peptide,point_ref,conf1,conf2,False) 1045 1046 return Transformation.Translation(cm2), q.asRotation(), Transformation.Translation(-cm1), rms
1047
1048 - def angularDistance(self, chain):
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
1097 - def screwMotionAnalysis(self, chain):
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 # SPATIAL DENSITY ANALYSIS 1171 #####################################################################################
1172 -class SpatialDensity(Analysis):
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 # The minimal set of input parameters names necessary to perform the analysis. 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 # Indicate whether this analysis can be estimated or not. 1229 canBeEstimated = True 1230
1231 - def __init__(self):
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
1236 - def initialize(self):
1237 """Initializes the analysis (e.g. parses and checks input parameters, set some variables ...). 1238 """ 1239 1240 # The input parameters are parsed. 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 # self.testLAMP = N.zeros((20,20,20), N.Float) 1292 1293 if self.pyroServer != 'monoprocessor': 1294 # The attribute trajectory is removed because it can not be pickled by Pyro. 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 # Load the whole trajectory set. 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 # The configuration is updated to the current trajectory step. 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 # The group is not considered if it falls outside the desired r interval. 1342 if (r < self.rMin) or (r > self.rMax): 1343 continue 1344 1345 x, y, z = changeBasis(tcom,com,vi,vj,vk) 1346 1347 # xBin = int((x + 0.5)/0.05) 1348 # yBin = int((y + 0.5)/0.05) 1349 # zBin = int((z + 0.5)/0.05) 1350 1351 # try: 1352 # self.testLAMP[xBin, yBin, zBin] += 1.0 1353 # except: 1354 # continue 1355 1356 # The spherical theta is computed 1357 theta = N.arccos(z/r) 1358 # The group is not considered if it falls outside the desired theta interval 1359 if (theta < self.thetaMin) or (theta > self.thetaMax): 1360 continue 1361 1362 # The spherical phi is computed 1363 phi = N.arctan2(y,x) 1364 # The group is not considered if it falls outside the desired phi interval 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 # The Spatial density is incremented by one unit. 1384 self.SD[rBin, thetaBin, phiBin] += 1.0
1385
1386 - def finalize(self):
1387 """Finalizes the calculations (e.g. averaging the total term, output files creations ...). 1388 """ 1389 1390 # for k in range(self.testLAMP.shape[2]): 1391 # for j in range(self.testLAMP.shape[1]): 1392 # for i in range(self.testLAMP.shape[0]): 1393 # print self.testLAMP[i,j,k] 1394 1395 self.SD /= float(self.nFrames*self.nTriplets*self.nGroups) 1396 1397 # The NetCDF output file is opened for writing. 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 # Some dimensions are created. 1403 outputFile.createDimension('NTHETAVALUES', self.nThetaBins) 1404 outputFile.createDimension('NPHIVALUES', self.nPhiBins) 1405 1406 # Creation of the NetCDF output variables. 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 # Creates an ASCII version of the NetCDF output file. 1430 convertNetCDFToASCII(inputFile = self.outputSD,\ 1431 outputFile = os.path.splitext(self.outputSD)[0] + '.cdl',\ 1432 variables = asciiVar)
1433
1434 - def constructBasisFromAtoms(self, triplet):
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