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

# Source Code for Module nMOLDYN.Analysis.NMR

      1  """Collections of classes for the determination of NMR-related properties.
      2
      3  Classes:
      4
      5      * OrderParameter             : sets up an order parameter analysis.
      6      * OrderParameterContactModel : sets up an order parameter analysis using the contact model approach.
      7  """
      8
      9  # The python distribution modules
     10  import copy
     11  import math
     12  from numpy import linalg
     13  import operator
     14  import os
     15  import sys
     16  from time import asctime
     17  from timeit import default_timer
     18
     19  # The ScientificPython modules
     20  from Scientific import N
     21  from Scientific.Geometry import Vector
     22  from Scientific.IO.NetCDF import NetCDFFile
     23
     24  # The MMTK distribution modules
     25  from MMTK.Collections import Collection
     26  from MMTK import Units
     27  from MMTK.NucleicAcids import NucleotideChain
     28  from MMTK.Proteins import PeptideChain, Protein
     29
     30  # The nMOLDYN modules
     31  from nMOLDYN.Analysis.Analysis import Analysis
     32  from nMOLDYN.Core.Error import Error
     33  from nMOLDYN.Core.IOFiles import convertNetCDFToASCII
     34  from nMOLDYN.Core.Logger import LogMessage
     35  from nMOLDYN.Core.Mathematics import correlation
     36  from nMOLDYN_S2 import computeS2
     37
     38  #####################################################################################
     39  # ORDER PARAMETER ANALYSIS
     40  #####################################################################################
     41 -class OrderParameter(Analysis):
     42      """Sets up an order parameter analysis.
     43
     44      A Subclass of nMOLDYN.Analysis.Analysis.
     45
     46      Constructor: OrderParameter(|parameters| = None)
     47
     48      Arguments:
     49
     50          - |parameters| -- a dictionnary of the input parameters, or 'None' to set up the analysis without parameters.
     51              * trajectory         -- a trajectory file name or an instance of MMTK.Trajectory.Trajectory class.
     52              * timeinfo           -- a string of the form 'first:last:step' where 'first' is an integer specifying the first frame
     53                                      number to consider, 'last' is an integer specifying the last frame number to consider and
     54                                      'step' is an integer specifying the step number between two frames.
     55              * group              -- a selection string specifying the groups of atoms that will define the vectors on which the
     56                                      analysis will be computed. Each group must contain two and only two atoms.
     57              * atomorder          -- a string of the form 'atom1,atom2,atom3' where 'atom1', 'atom2' and 'atom3' are
     58                                      respectively the MMTK atom names of the atoms in the way they should be ordered.
     59              * op                 -- the output NetCDF file name. A CDL version of this file will also be generated with the '.cdl' extension
     60                                      instead of the '.nc' extension.
     61              * pyroserver         -- a string specifying if Pyro will be used and how to run the analysis.
     62
     63      Running modes:
     64
     65          - To run the analysis do: a.runAnalysis() where a is the analysis object.
     66          - To estimate the analysis do: a.estimateAnalysis() where a is the analysis object.
     67          - To save the analysis to 'file' file name do: a.saveAnalysis(file) where a is the analysis object.
     68
     69
     70      Comments:
     71
     72          - This code is based on a first implementation made by Vania Calandrini.
     73      """
     74
     75      # The minimal set of input parameters names necessary to perform the analysis.
     76      inputParametersNames = ('trajectory',\
     77                              'timeinfo',\
     78                              'bond',\
     79                              'atomorder',\
     80                              'op',\
     81                              'pyroserver',\
     82                              )
     83
     84      shortName = 'OP'
     85
     86      # Indicate whether this analysis can be estimated or not.
     87      canBeEstimated = True
     88
     89 -    def __init__(self):
     90          """The constructor. Insures that the class can not be instanciated directly from here.
     91          """
     92          raise Error('This class can not be instanciated directly.')
     93
     94 -    def initialize(self):
     95          """Initializes the analysis (e.g. parses and checks input parameters, set some variables ...).
     96          """
     97
     98          # The input parameters are parsed.
     99          self.parseInputParameters()
    100
    101          if self.pyroServer != 'monoprocessor':
    102              if self.statusBar is not None:
    103                  LogMessage('warning', 'The job status bar will be desactivated for compatibility with Pyro module.')
    104                  self.statusBar = None
    105
    106          self.buildTimeInfo()
    107
    108          self.bond = self.groupSelection(self.universe, self.bond)
    109          self.bond = [g for g in self.bond if g.numberOfAtoms() == 2]
    110
    111          if self.atomOrder is None:
    112              self.bond = [sorted(g, key = operator.attrgetter('name')) for g in self.bond]
    113
    114          else:
    115              orderedBonds = []
    116              for g in self.bond:
    117                  gr = []
    118                  for atName in self.atomOrder:
    119                      found = False
    120                      for at in g:
    121                          if at.name == atName:
    122                              found = True
    123                              gr.append(at)
    124                              if len(gr) == 2:
    125                                  orderedBonds.append(gr)
    126                                  gr = []
    127                      else:
    128                          if not found:
    129                              raise Error('No atom corresponding %s could be found.' % atName)
    130
    131              self.bond = orderedBonds
    132
    133          if len(self.bond) == 0:
    134              raise Error('Your group selection does not contain any bonds.')
    135
    136          self.group = self.bond
    137          self.nBonds = self.nGroups = len(self.bond)
    138
    139          self.preLoadTrajectory('atom')
    140
    141          self.referenceDirection = Vector(0.0,0.0,1.0)
    142
    143          # The results are stored in dictionnary because the order in which the bonds are treated does not
    144          # always follow the structure.
    145          self.P2 = {}
    146          self.S2 = {}
    147
    148          if self.pyroServer != 'monoprocessor':
    149              # The attribute trajectory is removed because it can not be pickled by Pyro.
    150              delattr(self, 'trajectory')
    151
    152 -    def calc(self, bondIndex, trajname):
    153          """Calculates the contribution for one group.
    154
    155          @param bondIndex: the index of the group in |self.bond| list.
    156          @type bondIndex: integer.
    157
    158          @param trajname: the name of the trajectory file name.
    159          @type trajname: string
    160          """
    161
    162          # The atoms forming the bond.
    163          at1, at2 = self.bond[bondIndex]
    164
    165          if self.preLoadedTraj is None:
    166              if self.pyroServer == 'monoprocessor':
    167                  t = self.trajectory
    168
    169              else:
    170                  # Load the whole trajectory set.
    171                  t = Trajectory(None, trajname, 'r')
    172
    173              at1Traj = t.readParticleTrajectory(at1, first = self.first, last = self.last, skip = self.skip)
    174              at2Traj = t.readParticleTrajectory(at2, first = self.first, last = self.last, skip = self.skip)
    175
    176          else:
    177              at1Traj = self.preLoadedTraj[at1]
    178              at2Traj = self.preLoadedTraj[at2]
    179
    180          costheta = N.zeros((self.nFrames,), N.Float)
    181          sinphi   = N.zeros((self.nFrames,), N.Float)
    182          cosphi   = N.zeros((self.nFrames,), N.Float)
    183          sintheta = N.zeros((self.nFrames,), N.Float)
    184
    185          for frameIndex in range(self.nFrames):
    186
    187              pVect = Vector(self.universe.distanceVector(at1Traj[frameIndex], at2Traj[frameIndex]))
    188              pVect = pVect.normal()
    189
    190              costheta[frameIndex] = pVect * self.referenceDirection
    191              sintheta[frameIndex] = pVect.cross(self.referenceDirection).length()
    192              cosphi[frameIndex] = (pVect[0]/sintheta[frameIndex])
    193              sinphi[frameIndex] = (pVect[1]/sintheta[frameIndex])
    194
    195          tr2         = 3.0*costheta**2 - 1.
    196          cos2phi     = 2.0*cosphi**2-1.
    197          sin2phi     = 2.0*sinphi*cosphi
    198          cossintheta = costheta*sintheta
    199          sintheta_sq = sintheta**2
    200
    201          # calcul de <P2(cos(theta))> en termes de somme de fonctions de correlation
    202          # d'harmoniques spheriques (theoreme d'addition des harmoniques spheriques).
    203          P2 = (0.25*correlation(tr2) + 3.00*correlation(cosphi*cossintheta) + \
    204                3.00*correlation(sinphi*cossintheta) + 0.75*correlation(cos2phi*sintheta_sq) + \
    205                0.75*correlation(sin2phi*sintheta_sq))
    206
    207          # calcul du parametre d'ordre S^2 (limite pour t->infini de <P2(cos(theta))>).
    208          S2 = (0.75 * (N.sum(cos2phi*sintheta_sq)**2 + N.sum(sin2phi*sintheta_sq)**2) + \
    209                3.00 * (N.sum(cosphi*cossintheta)**2 + N.sum(sinphi*cossintheta)**2) +
    210                0.25 * N.sum(tr2)**2) / self.nFrames**2
    211
    212          return P2, S2
    213
    214 -    def combine(self, bondIndex, x):
    215          """
    216          """
    217          # The atoms forming the bond.
    218          at1, at2 = self.bond[bondIndex]
    219          if isinstance(at1.topLevelChemicalObject(), (Protein, PeptideChain, NucleotideChain)):
    220              bondNumber = at1.parent.parent.sequence_number
    221          else:
    222              bondNumber = bondIndex + 1
    223
    224          # calcul de <P2(cos(theta))> en termes de somme de fonctions de correlation
    225          # d'harmoniques spheriques (theoreme d'addition des harmoniques spheriques).
    226          self.P2[bondNumber] = x[0]
    227
    228          # calcul du parametre d'ordre S^2 (limite pour t->infini de <P2(cos(theta))>).
    229          self.S2[bondNumber] = x[1]
    230
    231 -    def finalize(self):
    232          """Finalizes the calculations (e.g. averaging the total term, output files creations ...).
    233          """
    234
    235          outputFile = NetCDFFile(self.outputOP, 'w')
    236          outputFile.title = self.__class__.__name__
    237          outputFile.jobinfo = self.information + '\nOutput file written on: %s\n\n' % asctime()
    238
    239          outputFile.createDimension('NBONDS', self.nBonds)
    240          outputFile.createDimension('NTIMES', self.nFrames)
    241
    242          TIMES = outputFile.createVariable('time', N.Float, ('NTIMES',))
    243          TIMES[:] = self.times
    244          TIMES.units = 'ps'
    245
    246          BONDS = outputFile.createVariable('bond_id', N.Int, ('NBONDS',))
    247          P2 = outputFile.createVariable('p2', N.Float, ('NBONDS','NTIMES'))
    248          P2AVG = outputFile.createVariable('p2-bondavg', N.Float, ('NTIMES',))
    249          S2   = outputFile.createVariable('s2', N.Float, ('NBONDS',))
    250          p2Avg = N.zeros((self.nFrames),N.Float)
    251          comp = 0
    252          for ind in sorted(self.S2.keys()):
    253              BONDS[comp] = ind
    254              S2[comp] = self.S2[ind]
    255              P2[comp,:] = self.P2[ind]
    256              N.add(p2Avg,self.P2[ind],p2Avg)
    257              comp += 1
    258
    259          P2AVG[:] = p2Avg/float(self.nBonds)
    260
    261          asciiVar = sorted(outputFile.variables.keys())
    262
    263          outputFile.close()
    264
    265          self.toPlot = {'netcdf' : self.outputOP, 'xVar' : 'pair', 'yVar' : 'S2'}
    266
    267          # Creates an ASCII version of the NetCDF output file.
    268          convertNetCDFToASCII(inputFile = self.outputOP,\
    269                               outputFile = os.path.splitext(self.outputOP)[0] + '.cdl',\
    270                               variables = asciiVar)
    271
    272  #####################################################################################
    273  # ORDER PARAMETER ANALYSIS USING CONTACT MODEL
    274  #####################################################################################
    275 -class OrderParameterContactModel(Analysis):
    276      """Sets up an order parameter analysis using the contact model .
    277
    278      A Subclass of nMOLDYN.Analysis.Analysis.
    279
    280      Constructor: OrderParameterContactModel(|parameters| = None)
    281
    282      Arguments:
    283
    284          - |parameters| -- a dictionnary of the input parameters, or 'None' to set up the analysis without parameters.
    285              * trajectory -- a trajectory file name or an instance of MMTK.Trajectory.Trajectory class.
    286              * timeinfo   -- a string of the form 'first:last:step' where 'first' is an integer specifying the first frame
    287                              number to consider, 'last' is an integer specifying the last frame number to consider and
    288                              'step' is an integer specifying the step number between two frames.
    289              * opcm       -- the output NetCDF file name. A CDL version of this file will also be generated with the '.cdl' extension
    290                              instead of the '.nc' extension.
    291              * pyroserver -- a string specifying if Pyro will be used and how to run the analysis.
    292
    293      Running modes:
    294
    295          - To run the analysis do: a.runAnalysis() where a is the analysis object.
    296          - To estimate the analysis do: a.estimateAnalysis() where a is the analysis object.
    297          - To save the analysis to 'file' file name do: a.saveAnalysis(file) where a is the analysis object.
    298
    299      Comments:
    300
    301          - This code is adapted from the s2predict code developped by F. Zhang and R. Bruschweiler and available in:
    302            http://nmr.clarku.edu/software/S2/s2predict.html
    303
    304          - For more details about the method: Zhang, F., Bruschweiler, R. J. AM. Chem. Soc. 2002, 124, 12654-12655.
    305      """
    306
    307      # The minimal set of input parameters names necessary to perform the analysis.
    308      inputParametersNames = ('trajectory',\
    309                              'timeinfo',\
    310                              'opcm',\
    311                              'pyroserver',\
    312                              )
    313
    314      shortName = 'OPCM'
    315
    316      # Indicate whether this analysis can be estimated or not.
    317      canBeEstimated = True
    318
    319 -    def __init__(self):
    320          """The constructor. Insures that the class can not be instanciated directly from here.
    321          """
    322          raise Error('This class can not be instanciated directly.')
    323
    324 -    def initialize(self):
    325          """Initializes the analysis (e.g. parses and checks input parameters, set some variables ...).
    326          """
    327
    328          # The input parameters are parsed.
    329          self.parseInputParameters()
    330
    331          if self.pyroServer != 'monoprocessor':
    332              if self.statusBar is not None:
    333                  LogMessage('warning', 'The job status bar will be desactivated for compatibility with Pyro module.')
    334                  self.statusBar = None
    335
    336          self.buildTimeInfo()
    337
    338          self.preLoadTrajectory('frame')
    339
    340          if self.universe.basisVectors() is None:
    341              raise Error('The universe must be periodic for this kind of calculation.')
    342
    343          self.S2 = {}
    344          self.sequence = {}
    345          self.hnLookup = {}
    346          for obj in self.universe.objectList():
    347              if isinstance(obj, (PeptideChain, Protein)):
    348                  self.sequence[obj] = []
    349                  self.hnLookup[obj] = []
    350                  comp = 1
    351                  for r in obj.residues()[1:]:
    352                      if r.symbol.lower() != 'pro':
    353                          self.sequence[obj].append(r.sequence_number)
    354                          self.hnLookup[obj].append(comp)
    355                      comp += 1
    356                  self.S2[obj] = N.zeros((len(self.sequence[obj]),self.nFrames), N.Float)
    357
    358          # Case where no protein or peptide chain was found in the universe.
    359          if not self.sequence:
    360              raise Error('The universe must contains at least one peptide chain to perform the analysis.')
    361
    362          if self.pyroServer != 'monoprocessor':
    363              # The attribute trajectory is removed because it can not be pickled by Pyro.
    364              delattr(self, 'trajectory')
    365
    366 -    def calc(self, frameIndex, trajname):
    367          """Calculates the contribution for one group.
    368
    369          @param frameIndex: the index of the frame in |self.frameIndexes| array.
    370          @type frameIndex: integer.
    371
    372          @param trajname: the name of the trajectory file name.
    373          @type trajname: string
    374          """
    375
    376          frame = self.frameIndexes[frameIndex]
    377
    378          if self.preLoadedTraj is None:
    379              if self.pyroServer == 'monoprocessor':
    380                  t = self.trajectory
    381
    382              else:
    383                  # Load the whole trajectory set.
    384                  t = Trajectory(None, trajname, 'r')
    385
    386              conf = t.configuration[frame]
    387              self.universe.setConfiguration(conf)
    388
    389          else:
    390              # The configuration is updated to the current trajectory step.
    391              self.universe.setConfiguration(self.preLoadedTraj[frameIndex])
    392
    393          directCell = N.ravel(N.array([v for v in self.universe.basisVectors()]))
    394          reverseCell = N.ravel(N.transpose(N.array([v for v in self.universe.reciprocalBasisVectors()])))
    395
    396          S2 = {}
    397          for obj in self.universe.objectList():
    398              if isinstance(obj, (PeptideChain, Protein)):
    399                  S2[obj] = N.zeros((len(self.sequence[obj])), N.Float)
    400
    401          for obj in self.sequence.keys():
    402              for v in range(len(self.hnLookup[obj])):
    403                  aa_no = self.hnLookup[obj][v]
    404                  aa_no_minus = aa_no - 1
    405
    406                  H = obj.backbone()[aa_no].H
    407                  O = obj.backbone()[aa_no_minus].O
    408
    409                  heavyAtoms = Collection()
    410                  for m in range(len(obj.residues())):
    411                      if (m != aa_no) and (m != aa_no_minus):
    412                          for atom in obj.residues()[m].atomList():
    413                              if atom.symbol != 'H':
    414                                  heavyAtoms.addObject(atom)
    415
    416                  # The indexes of the selected atoms.
    417                  indexes = N.zeros((heavyAtoms.numberOfAtoms(),), typecode = N.Int)
    418                  comp = 0
    419                  for at in heavyAtoms.atomList():
    420                      indexes[comp] = at.index
    421                      comp += 1
    422
    423                  val = computeS2(H.index,\
    424                                  O.index,\
    425                                  self.universe.contiguousObjectConfiguration().array,\
    426                                  directCell,\
    427                                  reverseCell,\
    428                                  indexes,\
    429                                  Units.Ang)
    430
    431                  S2[obj][v] = val
    432
    433          if self.preLoadedTraj is not None:
    434              del self.preLoadedTraj[frameIndex]
    435
    436          return S2
    437
    438 -    def combine(self, frameIndex, x):
    439          """
    440          """
    441          for obj in self.universe.objectList():
    442              if isinstance(obj, (PeptideChain, Protein)):
    443                  self.S2[obj][:,frameIndex] = x[obj]
    444
    445 -    def finalize(self):
    446          """Finalizes the calculations (e.g. averaging the total term, output files creations ...).
    447          """
    448
    449          # The NetCDF output file is opened for writing.
    450          outputFile       = NetCDFFile(self.outputOPCM, 'w')
    451          outputFile.title = self.__class__.__name__
    452          outputFile.jobinfo = self.information + '\nOutput file written on: %s\n\n' % asctime()
    453
    454          outputFile.createDimension('TIME', self.nFrames)
    455          TIMES = outputFile.createVariable('time', N.Float, ('TIME',))
    456          TIMES[:] = self.times[:]
    457          TIMES.units = 'ps'
    458
    459          for obj in self.sequence:
    460              outputFile.createDimension('SEQ%s' % obj.name, len(self.sequence[obj]))
    461
    462              SEQUENCE = outputFile.createVariable('%s_sequence' % obj.name, N.Int, ('SEQ%s' % obj.name,))
    463              SEQUENCE[:] = N.array(self.sequence[obj])
    464              SEQUENCE.units = 'unitless'
    465
    466              S2 = outputFile.createVariable('%s_s2' % obj.name, N.Float, ('SEQ%s' % obj.name,'TIME'))
    467              S2[:] = self.S2[obj][:,:]
    468              S2.units = 'unitless'
    469
    470              S2AVG = outputFile.createVariable('%s_s2_timeavg' % obj.name, N.Float, ('SEQ%s' % obj.name,))
    471              S2AVG[:] = self.S2[obj][:,:].sum(1)/float(self.nFrames)
    472              S2AVG.units = 'unitless'
    473
    474          asciiVar = sorted(outputFile.variables.keys())
    475
    476          outputFile.close()
    477
    478          self.toPlot = None
    479
    480          # Creates an ASCII version of the NetCDF output file.
    481          convertNetCDFToASCII(inputFile = self.outputOPCM,\
    482                               outputFile = os.path.splitext(self.outputOPCM)[0] + '.cdl',\
    483                               variables = asciiVar)
    484