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