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

# Source Code for Module nMOLDYN.Analysis.Analysis

       1  """
       2  This modules implements the base class for all the analysis available in nMOLDYN.
       3  """
       4
       5  # The python distribution modules
       6  from distutils.version import LooseVersion
       7  import getpass
       8  import os
       9  import platform
      10  from random import sample, uniform
      11  import re
      12  import sys
      13  from time import asctime
      14  from timeit import default_timer
      15
      16  # The ScientificPython modules
      17  from Scientific import N
      18  from Scientific.Geometry import Tensor, Vector
      19
      20  # The MMTK distribution modules
      21  from MMTK import Atom, AtomCluster, Molecule
      22  from MMTK.Collections import Collection
      23  from MMTK.NucleicAcids import NucleotideChain
      24  from MMTK.ParticleProperties import ParticleScalar
      25  from MMTK.PDB import PDBConfiguration
      26  from MMTK.Proteins import PeptideChain, Protein
      27  from MMTK.Trajectory import Trajectory, TrajectorySet
      28  from MMTK import Units
      29
      30  # The nMOLDYN modules
      31  import nMOLDYN
      32  from nMOLDYN.Core.Chemistry import belongToAnAmine, belongToAHydroxy, belongToAMethyl, belongToAThiol
      33  from nMOLDYN.Core.Error import Error
      34  from nMOLDYN.Core.IOFiles import TemporaryFile
      35  from nMOLDYN.Core.Logger import LogMessage
      36  from nMOLDYN.Core.Mathematics import randomVector
      37  from nMOLDYN.Core.Preferences import PREFERENCES
      38
      39  # Dictionnary describing the residus contents of each chemical family (e.g. hydrophobic ...)
      40  residusChemFamily = {'acidic'      : ('Asp','Glu'),
      41                       'aliphatic'   : ('Ile','Leu','Val'),
      42                       'aromatic'    : ('His','Phe','Trp','Tyr'),
      43                       'basic'       : ('Arg','His','Lys'),
      44                       'charged'     : ('Arg','Asp','Glu','His','Lys'),
      45                       'hydrophobic' : ('Ala','Cys','Cyx','Gly','His','Ile','Leu','Lys','Met','Phe','Thr','Trp','Tyr','Val'),
      46                       'polar'       : ('Arg','Asn','Asp','Cys','Gln','Glu','His','Lys','Ser','Thr','Trp','Tyr'),
      47                       'small'       : ('Ala','Asn','Asp','Cys','Cyx','Gly','Pro','Ser','Thr','Val')}
      48
      49  # Package path
      50  nmoldyn_package_path = os.path.dirname(os.path.split(__file__)[0])
      51
      52 -class Analysis:
      53      """Base class for all analysis defined in nMOLDYN.
      54
      55      The class Analysis is an abstract-base-class that defines attributes and methods common to all the analysis
      56      available in nMOLDYN. To set up an analysis object, use one of its subclass.
      57      """
      58
      59
      60 -    def __init__(self, parameters = None, statusBar = None):
      61          """
      62          The constructor.
      63
      64          @param parameters: a dictionnary that contains parameters of the selected analysis.
      65          @type parameters: dict
      66
      67          @param statusBar: if not None, an instance of nMOLDYN.GUI.Widgets.StatusBar. Will attach a status bar to the
      68              selected analysis.
      69          @type statusBar: instance of nMOLDYN.GUI.Widgets.StatusBar
      70          """
      71
      72          # This variable will store the the the time taken for the analysis in seconds.
      73          self.chrono = None
      74
      75          try:
      76              self.displayProgressRate = int(PREFERENCES.progress_rate)
      77
      78          except ValueError:
      79              self.displayProgressRate = 10
      80
      81          if (self.displayProgressRate <= 1) or (self.displayProgressRate >= 100):
      82              self.displayProgressRate = 10
      83
      84          # In case of an anlysis run from the GUI, this variable will be the statusbar displaying hte job progress.
      85          if statusBar is None:
      86              self.statusBar = None
      87          else:
      88              self.statusBar = statusBar
      89
      90          if isinstance(parameters, dict):
      91              self.parameters = parameters
      92
      93          else:
      94              self.parameters = None
      95
      96 -    def setInputParameters(self, parameters):
      97          """Sets the input parameters dictionnary.
      98          """
      99
     100          if isinstance(parameters, dict):
     101              self.parameters = parameters
     102          else:
     103              raise Error('Bad input for parameters constructor argument. It must be a dictionnary')
     104
     105
     106 -    def parseInputParameters(self):
     107          """Parses the input parameters stored in |parameters| dictionnary.
     108
     109          @return: a dictionnary of the parsed parameters.
     110          @rtype: dict
     111          """
     112
     113          # First check for the occurence of some general parameters.
     114          if not self.parameters.has_key('version'):
     115              self.parameters['version'] = 'unknown'
     116
     117          if self.canBeEstimated:
     118              if not self.parameters.has_key('estimate'):
     119                  self.parameters['estimate'] = 'no'
     120          else:
     121              if self.parameters.has_key('estimate'):
     122                  if self.parameters['estimate'] == 'yes':
     123                      raise Error('%s analysis can not be estimated.' % self.__class__.__name__.split('_')[0])
     124              else:
     125                  self.parameters['estimate'] = 'no'
     126
     127          # First, checks whether some input parameters are not used by the analysis. If so, just raises a warning.
     128          for p in self.parameters.keys():
     129              if p.lower() in ['version','estimate']:
     130                  continue
     131
     132              if p.lower() not in [v.lower() for v in self.inputParametersNames]:
     133                  LogMessage('warning','%s is not a useful parameter for %s analysis.' % (p, self.__class__.__name__),['file','console'])
     134                  del self.parameters[p]
     135
     136          # Then, checks whether some analysis necessary parameters are missing from the input parameters.
     137          # If so, raises an error.
     138          for p in self.inputParametersNames:
     139              if p.lower() not in [k.lower() for k in self.parameters.keys()]:
     140                  raise Error('%s input parameter is missing for %s analysis.' % (p, self.__class__.__name__))
     141
     142          # Loop over the parameters names necessary for the analysis.
     143          for pName in self.parameters:
     144
     145              # Their corresponding input value.
     146              pValue = self.parameters[pName]
     147              pName = pName.lower()
     148
     149              if pName == 'armodelorder':
     150                  self.arModelOrder = int(pValue)
     151
     152              elif pName == 'atomorder':
     153
     154                  if pValue is None:
     155                      self.atomOrder = None
     156
     157                  elif isinstance(pValue, str):
     158
     159                      if pValue.strip().lower() == 'no':
     160                          self.atomOrder = None
     161
     162                      else:
     163                          try:
     164                              self.atomOrder = [p.strip() for p in pValue.split(',')]
     165                              if not self.atomOrder:
     166                                  raise
     167
     168                          except:
     169                              raise Error('Error when parsing % parameter: must be comma-separated strings.' % pName)
     170
     171                  else:
     172                      raise Error('Error when parsing %s parameter: must be a string.' % pName)
     173
     174              elif pName == 'bond':
     175
     176                  if pValue is None:
     177                      self.bond = 'all'
     178                  else:
     179                      self.bond = pValue
     180
     181              elif pName == 'deuteration':
     182
     183                  if pValue is None:
     184                      self.deuteration = 'no'
     185                  else:
     186                      self.deuteration = pValue
     187
     188              elif pName == 'differentiation':
     189
     190                  try:
     191                      self.differentiation = int(pValue)
     192
     193                  except:
     194                      raise Error('Error when parsing %s parameter: must be an integer.' % pName)
     195
     196              elif pName == 'estimate':
     197                  try:
     198                      p = pValue.strip().lower()
     199                      if p == 'yes':
     200                          self.estimate = True
     201                      elif p == 'no':
     202                          self.estimate = False
     203                      else:
     204                          raise
     205
     206                  except:
     207                      raise Error('Error when parsing %s parameter: must be an "yes" or "no".' % pName)
     208
     209              # 'fftWindow' input parameter name.
     210              elif pName == 'fftwindow':
     211                  try:
     212                      self.fftWindow = float(pValue)
     213                      if (self.fftWindow <= 0.0) or (self.fftWindow >= 100.0):
     214                          raise
     215
     216                  except:
     217                      raise Error('Error when parsing %s parameter: must be a float in ]0.0,100.0[.' % pName)
     218
     219              # 'filter' input parameter name.
     220              elif pName == 'filter':
     221                  try:
     222                      self.freqMin, self.freqMax = [float(v) for v in pValue.split(':')]
     223
     224                      if self.freqMin >= self.freqMax:
     225                          raise
     226
     227                  except:
     228                      raise Error('Error when parsing %s parameter: must be "float:float".' % pName)
     229
     230              elif pName == 'group':
     231
     232                  if pValue is None:
     233                      self.group = 'all'
     234                  else:
     235                      self.group = pValue
     236
     237              elif pName == 'normalize':
     238                  try:
     239                      p = pValue.strip().lower()
     240                      if p == 'yes':
     241                          self.normalize = True
     242                      elif p == 'no':
     243                          self.normalize = False
     244                      else:
     245                          raise
     246
     247                  except:
     248                      raise Error('Error when parsing %s parameter: must be an "yes" or "no".' % pName)
     249
     250              elif pName == 'phivalues':
     251
     252                  if not isinstance(pValue, str):
     253                      raise Error('Error when parsing %s parameter: must be a string.' % pName)
     254
     255                  try:
     256
     257                      self.phiMin, self.phiMax, self.dPhi = [float(v) for v in pValue.split(':')]
     258
     259                      if self.phiMin < -180.0:
     260                          raise
     261
     262                      if self.phiMax > 180.0:
     263                          raise
     264
     265                      if self.phiMax <= self.phiMin:
     266                          raise
     267
     268                      if self.dPhi > (self.phiMax - self.phiMin):
     269                          raise
     270
     271                      self.phiMin *= Units.deg
     272                      self.phiMax *= Units.deg
     273                      self.dPhi *= Units.deg
     274
     275                      phi = []
     276                      phi0 = self.phiMin
     277                      while phi0 <= self.phiMax + 1.0e-6:
     278                          phi.append(phi0)
     279                          phi0 += self.dPhi
     280
     281                      if not phi:
     282                          raise
     283
     284                      self.nPhiBins = len(phi) - 1
     285
     286                      # The rMax is readjusted to fit the last r values of the list.
     287                      self.phiMax = phi[-1]
     288
     289                      self.phiValues = self.phiMin + self.dPhi/2.0 + self.dPhi*N.arange(self.nPhiBins)
     290
     291                  except:
     292                      raise Error('Error when parsing %s parameter: wrong format for the string.' % pName)
     293
     294              # 'projection' input parameter name.
     295              # The format must be 'float,float,float'
     296              elif pName == 'projection':
     297
     298                  if pValue is None:
     299                      self.projection = None
     300
     301                  elif isinstance(pValue, str):
     302
     303                      if pValue.strip().lower() == 'no':
     304                          self.projection = None
     305
     306                      else:
     307
     308                          try:
     309                              self.projection = Vector([float(v) for v in pValue.split(',')]).normal()
     310
     311                          except:
     312                              raise Error('Error when parsing % parameter: must be "float,float,float".' % pName)
     313
     314                  else:
     315                      raise Error('Error when parsing %s parameter: must be a string.' % pName)
     316
     317              elif pName == 'pyroserver':
     318
     319                  if pValue is None:
     320                      self.pyroServer = 'monoprocessor'
     321                      self.pyroNodes = {}
     322
     323                  elif isinstance(pValue, str):
     324
     325                      try:
     326                          if pValue.strip().lower() in ['monoprocessor', 'no']:
     327                              self.pyroServer = 'monoprocessor'
     328                              self.pyroNodes = {}
     329
     330                          elif pValue[:14].strip().lower() == 'multiprocessor':
     331                              hostName, nAllocatedProcs = pValue.split('::')[1].split(':')
     332                              self.pyroServer = 'multiprocessor'
     333                              self.pyroNodes = {hostName.strip() : int(nAllocatedProcs.strip())}
     334
     335                          elif pValue[:7].strip().lower() == 'cluster':
     336                              raise Error('Pyroserver in cluster configuration not implemented yet.')
     337                              self.pyroServer = 'cluster'
     338                              self.pyroNodes = {}
     339                              self.pyroNodes.update([[vv[0].strip(),int(vv[1].strip())] for vv in [v.split(':') for v in pValue.split('::')[1].split(',')]])
     340
     341                      except:
     342                          raise Error('Error when parsing %s parameter string in %s case.' % (pName,pValue.strip().lower()))
     343
     344              # 'qShellValues' input parameter name.
     345              # It can be either directly a list (or tuple) of float either a formatted string from which a list of float will be constructed.
     346              # The string format is "qmin1:qmax1:qstep1,qmin2:qmax2:qstep2,qmin3:qmax3:qstep3 ..."
     347              elif pName == 'qshellvalues':
     348
     349                  try:
     350
     351                      if not isinstance(pValue, (list, tuple, str)):
     352                          raise
     353
     354                      if isinstance(pValue,(list,tuple)):
     355                          self.qShellValues = [float(v) for v in pValue]
     356
     357                      elif isinstance(pValue, str):
     358                          self.qShellValues = []
     359                          for v in pValue.split(';'):
     360                              qMin, qMax, dq = [float(vv) for vv in v.split(':')]
     361
     362                              if qMin < 0:
     363                                  raise
     364
     365                              if qMax <= qMin:
     366                                  raise
     367
     368                              if dq <= 0.0:
     369                                  raise
     370
     371                              if dq > (qMax - qMin):
     372                                  raise
     373
     374                              q0 = qMin
     375                              while q0 <= qMax + 1.0e-6:
     376                                  self.qShellValues.append(q0)
     377                                  q0 += dq
     378
     379                      self.qShellValues = sorted(list(set(self.qShellValues)))
     380                      if not self.qShellValues:
     381                          raise
     382
     383                  except:
     384                      raise Error('Error when parsing %s parameter".' % pName)
     385
     386              elif pName == 'qshellwidth':
     387                  try:
     388                      self.qShellWidth = float(pValue)
     389                      if self.qShellWidth <= 0.0:
     390                          raise
     391
     392                  except:
     393                      raise Error('Error when parsing %s parameter: must be a float.' % pName)
     394
     395              # 'qVectorsDirection' input parameter name.
     396              elif pName == 'qvectorsdirection':
     397
     398                  if pValue is None:
     399                      self.qVectorsDirection = None
     400
     401                  elif isinstance(pValue, str):
     402
     403                      if pValue.strip().lower() == 'no':
     404                          self.qVectorsDirection = None
     405
     406                      else:
     407                          try:
     408                              self.qVectorsDirection = [Vector([float(vv.strip()) for vv in v.strip().split(',')]).normal() for v in pValue.split(';')]
     409                          except:
     410                              raise Error('Error when parsing %s parameter: each direction must be of the form "x,y,z".' % pName)
     411
     412                  else:
     413                      raise Error('Error when parsing %s parameter: must be a string.' % pName)
     414
     415
     416              elif pName == 'qvectorsgenerator':
     417
     418                  try:
     419                      p = pValue.strip().lower()
     420                      if p not in ['3d isotropic', '2d isotropic', 'anisotropic']:
     421                              raise
     422
     423                      self.qVectorsGenerator = p
     424
     425                  except:
     426                      raise Error('Error when parsing %s parameter: must be one of "3d isotropic", "2d isotropic" or "anisotropic".' % pName)
     427
     428              elif pName == 'qvectorspershell':
     429
     430                  try:
     431                      self.qVectorsPerShell = int(pValue)
     432                      if self.qVectorsPerShell <= 0:
     433                          raise
     434
     435                  except:
     436                      raise Error('Error when parsing %s parameter: must be an integer.' % pName)
     437
     438              elif pName == 'referencedirection':
     439                  if pValue is not None:
     440                      try:
     441                          self.referenceDirection = Vector([float(v) for v in pValue.split(',')]).normal()
     442
     443                      except:
     444                          raise Error('Error when parsing %s parameter: must be "float,float,float".' % pName)
     445
     446                  else:
     447                      self.referenceDirection = None
     448
     449              elif pName == 'referenceframe':
     450                  self.referenceFrame = int(pValue) - 1
     451
     452              elif pName == 'removetranslation':
     453                  try:
     454                      p = pValue.strip().lower()
     455                      if p == 'yes':
     456                          self.removeTranslation = True
     457                      elif p == 'no':
     458                          self.removeTranslation = False
     459                      else:
     460                          raise
     461
     462                  except:
     463                      raise Error('Error when parsing %s parameter: must be an "yes" or "no".' % pName)
     464
     465              elif pName == 'rvalues':
     466
     467                  if not isinstance(pValue, str):
     468                      raise Error('Error when parsing %s parameter: must be a string.' % pName)
     469
     470                  try:
     471
     472                      self.rMin, self.rMax, self.dr = [float(v) for v in pValue.split(':')]
     473
     474                      if self.rMin < 0.0:
     475                          raise
     476
     477                      if self.rMax <= self.rMin:
     478                          raise
     479
     480                      if self.dr > (self.rMax - self.rMin):
     481                          raise
     482
     483                      r = []
     484                      r0 = self.rMin
     485                      while r0 <= self.rMax + 1.0e-6:
     486                          r.append(r0)
     487                          r0 += self.dr
     488
     489                      if not r:
     490                          raise
     491
     492                      self.nBins = len(r) - 1
     493
     494                      # The rMax is readjusted to fit the last r values of the list.
     495                      self.rMax = r[-1]
     496
     497                      self.rValues = self.rMin + self.dr/2.0 + self.dr*N.arange(self.nBins)
     498
     499                  except:
     500                      raise Error('Error when parsing %s parameter: wrong format for the string.' % pName)
     501
     502              elif pName == 'stepwiserbt':
     503                  try:
     504                      p = pValue.strip().lower()
     505                      if p == 'yes':
     506                          self.stepwiseRBT = True
     507                      elif p == 'no':
     508                          self.stepwiseRBT = False
     509                      else:
     510                          raise
     511
     512                  except:
     513                      raise Error('Error when parsing %s parameter: must be an "yes" or "no".' % pName)
     514
     515              elif pName == 'subset':
     516
     517                  if pValue is None:
     518                      self.subset = 'all'
     519                  else:
     520                      self.subset = pValue
     521
     522              elif pName == 'thetavalues':
     523
     524                  if not isinstance(pValue, str):
     525                      raise Error('Error when parsing %s parameter: must be a string.' % pName)
     526
     527                  try:
     528
     529                      self.thetaMin, self.thetaMax, self.dTheta = [float(v) for v in pValue.split(':')]
     530
     531                      if self.thetaMin < 0.0:
     532                          raise
     533
     534                      if self.thetaMin > 180.0:
     535                          raise
     536
     537                      if self.thetaMax <= self.thetaMin:
     538                          raise
     539
     540                      if self.dTheta > (self.thetaMax - self.thetaMin):
     541                          raise
     542
     543                      self.thetaMin *= Units.deg
     544                      self.thetaMax *= Units.deg
     545                      self.dTheta *= Units.deg
     546
     547                      theta = []
     548                      theta0 = self.thetaMin
     549                      while theta0 <= self.thetaMax + 1.0e-6:
     550                          theta.append(theta0)
     551                          theta0 += self.dTheta
     552
     553                      if not theta:
     554                          raise
     555
     556                      self.nThetaBins = len(theta) - 1
     557
     558                      # The rMax is readjusted to fit the last r values of the list.
     559                      self.thetaMax = theta[-1]
     560
     561                      self.thetaValues = self.thetaMin + self.dTheta/2.0 + self.dTheta*N.arange(self.nThetaBins)
     562
     563                  except:
     564                      raise Error('Error when parsing %s parameter: wrong format for the string.' % pName)
     565
     566              elif pName == 'timeinfo':
     567
     568                  try:
     569                      self.first, self.last, self.skip = [int(v) for v in pValue.split(':')]
     570                      self.first = self.first - 1
     571
     572                  except:
     573                      raise Error('Error when parsing %s parameter: must be "int:int:int".' % pName)
     574
     575              elif pName == 'temperature':
     576                  try:
     577                      self.temperature = float(pValue)
     578                      if self.temperature <= 0.0:
     579                          raise
     580
     581                  except:
     582                      raise Error('Error when parsing %s parameter: must be a float > 0.' % pName)
     583
     584                      # 'trajectory' input parameter name.
     585
     586              elif pName == 'trajectory':
     587                  # If the value is a string, then load the corresponding MMTK object and stores it.
     588                  if isinstance(pValue, str):
     589                      # The extension of the trajectory to load.
     590                      ext = os.path.splitext(pValue)[1]
     591
     592                      # Trajectory set --> loads a TrajectorySet MMTK object.
     593                      if ext == '.ncs':
     594                          try:
     595                              # The trajectory set is opened for reading and its contents sent to |trajSet| list.
     596                              trajSetFile = open(pValue, 'r')
     597                              trajSet = [t.strip() for t in trajSetFile.readlines()]
     598                              trajSetFile.close()
     599
     600                              # Load the whole trajectory set.
     601                              self.trajectory = TrajectorySet(None, trajSet)
     602
     603                          except:
     604                              raise Error('Error when parsing %s parameter: unable to open trajectory set file %s.' % pValue)
     605
     606                      # Trajectory --> load a Trajectory.
     607                      elif ext == '.nc':
     608                          try:
     609                              # Load the whole trajectory set.
     610                              self.trajectory = Trajectory(None, pValue, 'r')
     611
     612                          except IOError:
     613                              raise Error('Error when parsing %s parameter: unable to open trajectory file %s.' % pValue)
     614
     615                      # Otherwise throw an error.
     616                      else:
     617                          raise Error('Error when parsing %s parameter: not a suitable Trajectory or TrajectorySet file name.' % pValue)
     618
     619                  # If the value is already a Trajectory MMTK object, then stores it.
     620                  elif isinstance(pValue, Trajectory):
     621                      self.trajectory = pValue
     622
     623                  else:
     624                      raise Error('Error when parsing %s parameter: neither a trajectory or trajectory set file name neither an instance of Trajectory class.' % pValue)
     625
     626                  self.trajectoryFilename = self.trajectory.filename
     627
     628                  self.universe = self.trajectory.universe
     629
     630                  # The universe contents.
     631                  setUniverseContents(self.universe)
     632
     633              elif pName == 'triplet':
     634
     635                  if pValue is None:
     636                      self.triplet = 'all'
     637                  else:
     638                      self.triplet = pValue
     639
     640              elif pName == 'version':
     641                  try:
     642                      if pName == 'unknown':
     643                          self.version = LooseVersion(vstring = '0.0')
     644                      else:
     645                          self.version = LooseVersion(vstring = self.parameters['version'])
     646
     647                  except:
     648                      raise Error('Error when parsing %s parameter: must be a string of the form "int.int.int".' % pName)
     649
     650              elif pName == 'weights':
     651                  self.weights = pValue
     652
     653              elif pName == 'wignerindexes':
     654
     655                  try:
     656                      self.wignerIndexes = [int(v) for v in pValue.split(',')]
     657                      if len(self.wignerIndexes) != 3:
     658                          raise
     659
     660                  except:
     661                      raise Error('Error when parsing %s parameter: must be "int:int:int".' % (pName))
     662
     663              elif pName in ['ac', 'ados', 'ara', 'avacf', 'cn', 'comt', 'dcsf', 'dcsfar', 'disf', 'disfar',\
     664                              'disfg', 'dos', 'eisf', 'gmft', 'isfg', 'msd', 'op', 'opcm', 'pbft', 'pdf', 'qha',\
     665                              'rbt','rcf', 'rmsd', 'rog', 'sd', 'sfa', 'scsf', 'sscsf', 'vacf']:
     666
     667                  setattr(self, 'output'+ pName.upper(), pValue)
     668
     669 -    def buildTimeInfo(self):
     670          """Builds some attributes related to the frame selection string. They will be used to define at which
     671          times a given analysis should be run.
     672          """
     673
     674          trajLength = len(self.trajectory)
     675
     676          if self.first < 0:
     677              self.first = 0
     678              LogMessage('warning', 'The first step must be an integer >= 1.\n\
     679              It will be set to 1 for the running analysis.',['console'])
     680
     681          if (self.last <= self.first) or (self.last > trajLength):
     682              self.last = len(self.trajectory)
     683              LogMessage('warning', 'The last step must be an integer in [1,%s].\n\
     684              It will be set to %s for the running analysis.' % (trajLength, trajLength), ['console'])
     685
     686          if (self.skip <= 0) or (self.skip >= trajLength):
     687              self.skip = 1
     688              LogMessage('warning', 'The skip step must be an integer in [1,%s].\n\
     689              It will be set to 1 for the running analysis.' % trajLength,['console'])
     690
     691          self.frameIndexes = N.arange(self.first, self.last, self.skip)
     692          if len(self.frameIndexes) == 0:
     693              self.first = 1
     694              self.last = trajLength
     695              self.skip = 1
     696              self.frameIndexes = N.arange(self.first,self.last,self.skip)
     697              LogMessage('warning','The frame selection is empty.\
     698              It will be set to %s:%s:%s for the running analysis.' % (1,trajLength,1),['console'])
     699
     700          if hasattr(self.trajectory,'time'):
     701              t = self.trajectory.time[self.first:self.last:self.skip]
     702
     703          else:
     704              t = range(trajLength)[self.first:self.last:self.skip]
     705
     706          if len(t) == 1:
     707              self.dt = 1.0
     708          else:
     709              if t[0] == t[1]:
     710                  self.dt = 1.0
     711              else:
     712                  self.dt = t[1] - t[0]
     713
     714          self.times = self.dt * N.arange(len(t))
     715
     716          self.nFrames = len(self.times)
     717
     718 -    def preLoadTrajectory(self, structure, differentiation = 1):
     719          """
     720          """
     721
     722          if '_estimate' in self.__class__.__name__:
     723              self.preLoadedTraj = None
     724              return
     725
     726          try:
     727              LogMessage('info', 'Trying to preload the trajectory. Please wait ...', ['console'])
     728              self.preLoadedTraj = {}
     729              if structure == 'atom':
     730                  test = N.zeros((self.nAtoms, self.nFrames, 3), N.Float)
     731                  del test
     732                  for atom in self.subset:
     733                      if differentiation == 0:
     734                          self.preLoadedTraj[atom] = self.trajectory.readParticleTrajectory(atom, first = self.first, last = self.last, skip = self.skip, variable = 'velocities').array
     735                      else:
     736                          self.preLoadedTraj[atom] = self.trajectory.readParticleTrajectory(atom, first = self.first, last = self.last, skip = self.skip).array
     737              elif structure == 'frame':
     738                  test = N.zeros((self.nFrames, self.nAtoms, 3), N.Float)
     739                  del test
     740                  for frameIndex in range(self.nFrames):
     741                      frame = self.frameIndexes[frameIndex]
     742                      self.preLoadedTraj[frameIndex] = self.trajectory.configuration[frame]
     743              LogMessage('info', 'Enough memory to preload the trajectory.', ['console'])
     744
     745          except:
     746              self.preLoadedTraj = None
     747              LogMessage('info', 'Not enough memory to preload the trajectory. The analysis will run with runtime %s NetCDF extraction.' % structure, ['console'])
     748
     749 -    def saveAnalysis(self, filename):
     750          """Saves the settings of an analysis to an output file.
     751
     752          @param filename: the name of the output file. If the extension is '.nmi' the output file will be a
     753              nMOLDYN input script otherwise the output file will be a nMOLDYN autostart script.
     754          @type: filename: string
     755          """
     756
     757          if filename[-4:].lower() == '.nmi':
     758
     759              # Opens the python script file for writing.
     760              pScript = open(filename, 'w')
     761              # Write the input file header.
     762              pScript.write('################################################################\n')
     763              pScript.write('# This is an automatically generated python-nMOLDYN run script #\n')
     764              pScript.write('################################################################\n\n')
     765
     766              # The line that will create an instance of the analysis child class.
     767              pScript.write(('\nanalysis = %s\n\n') % (self.__class__.__name__))
     768
     769              pScript.write('# General parameters\n')
     770              pScript.write(('version = "%s"\n') % nMOLDYN.__version__)
     771              pScript.write(('estimate = "%s"\n\n') % self.parameters['estimate'])
     772
     773              pScript.write('# Analysis-specific parameters\n')
     774              # Loop over the individual widget names.
     775              for pName in sorted(self.inputParametersNames):
     776                  pValue = self.parameters[pName]
     777                  # Special case of strings. Some quotes are added to the value.
     778                  if isinstance(pValue, str):
     779                      pScript.write(('%s = "%s"\n') % (pName, pValue))
     780                  else:
     781                      pScript.write(('%s = %s\n') % (pName, pValue))
     782
     783              pScript.close()
     784
     785          else:
     786
     787              # Opens the python script file for writing.
     788              pScript = open(filename, 'w')
     789              pScript.write('#!%s\n' % sys.executable)
     790
     791              # Write the input file header.
     792              pScript.write('################################################################\n')
     793              pScript.write('# This is an automatically generated python-nMOLDYN run script #\n')
     794              pScript.write('################################################################\n\n')
     795              # Write the 'import' lines.
     796              pScript.write('from MMTK import *\n')
     797              pScript.write(('from %s import %s\n\n') % (self.__module__, self.__class__.__name__))
     798
     799              # Writes the line that will initialize the parameters dictionnary.
     800              pScript.write('parameters = {}\n\n')
     801
     802              pScript.write('# General parameters\n')
     803              pScript.write('parameters[\'version\'] = "%s"\n' % nMOLDYN.__version__)
     804              pScript.write('parameters[\'estimate\'] = "%s"\n\n' % self.parameters['estimate'])
     805
     806              pScript.write('# Analysis-specific parameters\n')
     807              # Loop over the individual widget names.
     808              for pName in sorted(self.inputParametersNames):
     809                  pValue = self.parameters[pName]
     810                  # Special case of strings. Some quotes are added to the value.
     811                  if isinstance(pValue, str):
     812                      pScript.write(('parameters[\'%s\'] = "%s"\n') % (pName, pValue))
     813                  else:
     814                      pScript.write(('parameters[\'%s\'] = %s\n') % (pName, pValue))
     815
     816              # The line that will create an instance of the analysis child class.
     817              pScript.write(('\nanalysis = %s(parameters)') % (self.__class__.__name__))
     818              # The line that will actually launch the analysis.
     819              pScript.write('\nanalysis.runAnalysis()')
     820              pScript.close()
     821
     822              # Make the script executable on posix systems.
     823              if os.name == 'posix':
     824                  os.system('chmod u+x %s' % filename)
     825
     826 -    def runAnalysis(self):
     827          """Runs an analysis.
     828
     829          @return: a dictionnary of the form {'days' : d, 'hours' : h, 'minutes' : m, 'seconds' : s} specifying the
     830              time the analysis took in dayx, hours, minutes and seconds.
     831          @rtype: dict
     832          """
     833
     834          # The job progress in percentage.
     835          self.jobCounter = 0
     836
     837          # Initialize the analysis.
     838          self.initialize()
     839
     840          # Displays the informations about the running analysis.
     841          self.buildJobInfo()
     842
     843          # Checks the nMOLDYN version of the running analysis
     844          if self.version < LooseVersion(vstring = nMOLDYN.__version__):
     845              # The elapsed time for the (estimated) analysis is written in the log file.
     846              LogMessage('warning','The nMOLDYN version number of the script is unknown or deprecated.\n\
     847  Some analysis keywords name/values may have changed. Please check in case of troubles.',['file','console'])
     848
     849          # The analysis is actually processed here.
     850          self.internalRun()
     851
     852          # The ending date of the analysis.
     853          LogMessage('info','Job finished on '+asctime()+'.\n',['console','file'])
     854
     855          # The time taken for the analysis.
     856          timeTakenForAnalysis = self.analysisTime(self.chrono)
     857
     858          # The log message that displays the time taken for the analysis is a little bit different
     859          # depending on whether we are in estimate mode or not.
     860          if self.estimate:
     861              # The elapsed time for the (estimated) analysis is written in the log file.
     862              LogMessage('info','Estimated time for analysis: %(days)s days %(hours)s hours %(minutes)s minutes %(seconds)s seconds.\n' % timeTakenForAnalysis,['file','console'])
     863          else:
     864              # The elapsed time for the (estimated) analysis is written in the log file.
     865              LogMessage('info','Elapsed time for analysis: %(days)s days %(hours)s hours %(minutes)s minutes %(seconds)s seconds.\n' % timeTakenForAnalysis,['file','console'])
     866
     867          return timeTakenForAnalysis
     868
     869 -    def updateJobProgress(self, norm):
     870          """Check the progress of the running analysis and displays periodically on the console and the logfile how far is the analysis.
     871          Called each time a step of an analysis loop is achieved.
     872
     873          @param norm: the maximum number of steps of the analysis.
     874          @type: norm: integer
     875          """
     876
     877          # computes the old percentage
     878          t = int(100*self.jobCounter/norm)
     879
     880          oldProgress = (t/self.displayProgressRate)*self.displayProgressRate
     881
     882          # A new step has been done
     883          self.jobCounter += 1
     884
     885          # computes the new percentage
     886          t = int(100*self.jobCounter/norm)
     887          newProgress = (t/self.displayProgressRate)*self.displayProgressRate
     888
     889          if newProgress != oldProgress:
     890              LogMessage('info','Progress rate = %3d %%' % newProgress,['file','console'])
     891
     892          if self.statusBar is not None:
     893              self.statusBar.setValue(t)
     894
     895 -    def buildJobInfo(self):
     896          """Display on the console and in the log file the main ifnormation about the analysis to run.
     897          """
     898
     899          self.information = ""
     900
     901          # The ascii time when the job was launched.
     902          jobCreationTime = asctime()
     903
     904          self.information += '#'*90 +'\n'
     905          self.information += 'Job information for %s analysis.\n' % self.__class__.__name__
     906          self.information += '#'*90 +'\n\n'
     907
     908          # Some information about the analysis that will be written as a header of the NetCDF output file.
     909          self.information += 'Job launched on: %s\n\n' % jobCreationTime
     910          self.information += 'General informations\n'
     911          self.information += '--------------------\n'
     912          self.information += 'User: %s\n' % getpass.getuser()
     913          self.information += 'OS: %s-%s\n' % (platform.system(), platform.release())
     914          self.information += 'Processor: %s\n' % platform.machine()
     915          self.information += 'nMOLDYN version: %s\n' % self.parameters['version']
     916          self.information += 'Estimate run: %s\n\n' % self.parameters['estimate']
     917
     918          # Some informations for the logfile an output file.
     919          self.information += 'Parameters\n'
     920          self.information += '----------\n'
     921          self.information += '\n'.join(['%s = %s' % (k,self.parameters[k]) for k in self.parameters if k not in ['version','estimate']])
     922          self.information += '\n\n\n'
     923
     924          # Some informations about the various selections.
     925          if hasattr(self, 'subset'):
     926              self.information += 'Subset selection\n'
     927              self.information += '----------------\n'
     928              self.information += 'Number of atoms selected for analysis: %s\n\n' % len(self.subset)
     929
     930          if hasattr(self, 'deuteration'):
     931              self.information += 'Deuteration selection\n'
     932              self.information += '---------------------\n'
     933              self.information += 'Number of atoms selected for deuteration: %s\n\n' % len(self.deuteration)
     934
     935          if hasattr(self, 'group'):
     936              self.information += 'Group selection\n'
     937              self.information += '---------------------\n'
     938              self.information += 'Number of groups selected: %s\n\n' % len(self.group)
     939
     940          if hasattr(self, 'triplet'):
     941              self.information += 'Triplet selection\n'
     942              self.information += '---------------------\n'
     943              self.information += 'Number of triplets selected: %s\n\n' % len(self.triplet)
     944
     945          if hasattr(self, 'bond'):
     946              self.information += 'Bond selection\n'
     947              self.information += '---------------------\n'
     948              self.information += 'Number of bonds selected: %s\n\n' % len(self.bond)
     949
     950          self.information += '\n'
     951          self.information += 'Job status\n'
     952          self.information += '----------\n'
     953
     954          [LogMessage('info', l, ['file', 'console']) for l in self.information.splitlines()]
     955
     956 -    def analysisTime(self, time):
     957          """Converts a time in second in days, hours, minutes and seconds.
     958
     959          @param time: the time (in seconds) to convert.
     960          @type time: integer.
     961
     962          @return: a dictionnary of the form {'days' : d, 'hours' : h, 'minutes' : m, 'seconds' : s}
     963              where d, h, m and s are integers resulting respectively from the conversion of |time| in days, hours, minutes
     964              and seconds.
     965          @rtype: dict
     966          """
     967
     968          if time is None:
     969              return None
     970
     971          else:
     972              converted = {'days' : 0, 'hours' : 0, 'minutes' : 0, 'seconds' : 0}
     973
     974              converted['days'], time = divmod(time, 86400)
     975              converted['hours'], time = divmod(time, 3600)
     976              converted['minutes'], converted['seconds'] = divmod(time, 60)
     977
     978              return converted
     979
     980 -    def weightingScheme(self, universe, atoms, deuter, scheme = 'equal'):
     981          """Returns the weights of |atoms| MMTK collection of |universe| MMTK universe using the
     982          weighting scheme |scheme|.
     983
     984          @param universe: the MMTK universe.
     985          @type universe: instance of MMTK.Universe
     986
     987          @param atoms: the atoms to take into account when defining the weights.
     988          @type atoms: instance of MMTK.Collections.Collection
     989
     990          @param deuter: the hydrogen atoms that will be parametrized as deuterium atoms.
     991          @type deuter: instance of MMTK.Collections.Collection
     992
     993          @param scheme: a string equal to 'equal', 'mass', 'coherent' , 'incoherent' or 'atomicNumber' that specifies the weighting
     994              scheme to use.
     995          @type scheme: string
     996
     997          @return: the weights of the selected atoms.
     998          @rtype: an instance of MMTK.ParticleProperties.ParticledScalar
     999          """
    1000
    1001          try:
    1002              if scheme == 'equal':
    1003                  # MMTK ParticleScalar object containing 1 for each atom of the universe
    1004                  weights = ParticleScalar(universe, N.ones(universe.numberOfAtoms(), N.Float))
    1005
    1006              elif scheme == 'mass':
    1007                  # The atomic masses of the universe.
    1008                  weights = universe.masses()
    1009
    1010                  for atom in deuter:
    1011                      if abs(atom._mass - atom.mass_deut) < 1.e-10:
    1012                          LogMessage('warning','The atom %s is not sensitive to deuteration.' % str(atom),['file'])
    1013                      else:
    1014                          weights[atom] = atom.mass_deut
    1015
    1016              elif scheme == 'atomicNumber':
    1017                  try:
    1018                      # MMTK ParticleScalar object containing 1 for each atom of the universe
    1019                      weights = ParticleScalar(universe)
    1020                      for at in universe.atomList():
    1021                          weights[at] = at.atomic_number
    1022                  except KeyError:
    1023                      raise Error('One or more element with unknown atomic number in your selection.')
    1024
    1025              elif scheme == 'coherent':
    1026                  # getAtomScalarArray is a wrap for getParticleScalar function of Universe.py module
    1027                  # MMTK ParticleScalar object containing the bcoh.
    1028                  weights = universe.getAtomScalarArray('b_coherent')
    1029
    1030                  # The b value of the hydrogen to be deuterated are changed to b value of deuterium.
    1031                  for atom in deuter:
    1032                      if abs(atom.b_coherent - atom.b_coherent_deut) < 1.e-10:
    1033                          LogMessage('warning','The atom %s is not sensitive to deuteration.' % str(atom),['file'])
    1034                      else:
    1035                          weights[atom] = atom.b_coherent_deut
    1036
    1037              elif scheme == 'incoherent':
    1038                  # getAtomScalarArray is a wrap for getParticleScalar function of Universe.py module
    1039                  # MMTK ParticleScalar object containing the bincoh.
    1040                  weights = universe.getAtomScalarArray('b_incoherent')
    1041
    1042                  # The b value of the hydrogen to be deuterated are changed to b value of deuterium.
    1043                  for atom in deuter:
    1044                      if abs(atom.b_incoherent - atom.b_incoherent_deut) < 1.e-10:
    1045                          LogMessage('warning','The atom %s is not sensitive to deuteration.' % str(atom),['file'])
    1046                      else:
    1047                          weights[atom] = atom.b_incoherent_deut
    1048
    1049                  # b^2
    1050                  weights = weights*weights
    1051
    1052              else:
    1053                  raise Error('%s weighting scheme is unknown for %s analysis.' % (scheme, self.__class__.__name__))
    1054
    1055              # If there is an atom selection. A boolean mask is applied on the atoms to select.
    1056              if atoms != universe:
    1057                  weights = weights*atoms.booleanMask()
    1058
    1059              # The normalization factor is different for the coherent weighting scheme.
    1060              if scheme == 'coherent':
    1061                  weights = weights/N.sqrt((weights*weights).sumOverParticles())
    1062              else:
    1063                  weights = weights/weights.sumOverParticles()
    1064
    1065          except:
    1066              raise Error('Problem when setting up %s weighting scheme for %s analysis.' % (scheme, self.__class__.__name__))
    1067          else:
    1068              return weights
    1069
    1070 -    def subsetSelection(self, universe, selection):
    1071          """Returns a MMTK collection of atoms that matches |selection| selection string. Used to apply an analysis
    1072          to a subset of atoms.
    1073
    1074          @param universe: the universe on which the selection will be performed.
    1075          @type universe: instance of MMTK.Universe
    1076
    1077          @param selection: the selection string that will define the atoms to select.
    1078          @type selection: string
    1079
    1080          @return: a MMTK Collection of the atoms that matches |selection| selection string.
    1081          @rtype: instance of MMTK.Collections.Collection
    1082          """
    1083
    1084          if not isinstance(selection, (Collection, list, tuple, str)):
    1085              raise Error('Wrong format for the subset selection. Must be a string or a MMTK Collection' % str(selection))
    1086
    1087          if isinstance(selection, Collection):
    1088              return selection
    1089
    1090          elif isinstance(selection, (list, tuple)):
    1091              try:
    1092                  subsetSelection = Collection(selection)
    1093                  return subsetSelection
    1094
    1095              except:
    1096                  raise Error('The results of a subset selection must be a MMTK Collection.')
    1097
    1098          # If selection dictionnary is empty or equal to None then, by default, returns the whole universe as the subset.
    1099          elif selection.lower() == 'all':
    1100              return Collection(universe.atomList())
    1101
    1102          else:
    1103
    1104              # If the selection string contains "filename name-of-file" then it will be a selection from a file.
    1105              if selection.lower().strip()[0:8] == 'filename':
    1106                  subsetSelection = self.__parseFileSelectionString(universe, selection, 'subset')
    1107
    1108              elif selection.lower().strip()[0:10] == 'expression':
    1109                  subsetSelection = self.__parseExpressionSelectionString(universe, selection, 'subset')
    1110
    1111              # Otherwise the selection will be from an object.
    1112              else:
    1113                  subsetSelection = self.__parseObjectSelectionString(universe, selection)
    1114
    1115              try:
    1116                  subsetSelection = Collection(list(subsetSelection))
    1117
    1118              except:
    1119                  raise Error('The results of a subset selection must be a MMTK Collection.')
    1120
    1121              return subsetSelection
    1122
    1123 -    def deuterationSelection(self, universe, selection):
    1124          """Returns a MMTK collection of atoms that matches |selection| selection string. Used to switch the parameters
    1125          of a subset (or all) of hydrogen atoms to the parameters of deuterium in order to simulate deuterated system.
    1126
    1127          @param universe: the universe on which the selection will be performed.
    1128          @type universe: instance of MMTK.Universe
    1129
    1130          @param selection: the selection string that will define the atoms to select.
    1131          @type selection: string
    1132
    1133          @return: a MMTK Collection of the atoms that matches |selection| selection string.
    1134          @rtype: instance of MMTK.Collections.Collection
    1135          """
    1136
    1137          if not isinstance(selection, (Collection, list, tuple, str)):
    1138              raise Error('Wrong format for the deuteration selection. Must be a string or a MMTK Collection')
    1139
    1140          if isinstance(selection, Collection):
    1141              return selection
    1142
    1143          elif isinstance(selection, (list, tuple)):
    1144              try:
    1145                  deuterationSelection = Collection(selection)
    1146                  return deuterationSelection
    1147
    1148              except:
    1149                  raise Error('The results of a deuteration selection must be a MMTK Collection.')
    1150
    1151          elif selection.lower() == 'no':
    1152              return Collection()
    1153
    1154          elif selection.lower() == 'all':
    1155              return Collection([at for at in universe.atomList() if at.type.name == 'hydrogen'])
    1156
    1157          else:
    1158
    1159              # If the selection string contains "filename name-of-file" then it will be a selection from a file.
    1160              if selection.lower().strip()[0:8] == 'filename':
    1161                  deuterationSelection = self.__parseFileSelectionString(universe, selection, 'deuteration')
    1162
    1163              elif selection.lower().strip()[0:10] == 'expression':
    1164                  deuterationSelection = self.__parseExpressionSelectionString(universe, selection, 'deuteration')
    1165
    1166              # Otherwise the selection will be from an object.
    1167              else:
    1168                  deuterationSelection = self.__parseObjectSelectionString(universe, selection)
    1169
    1170              deuterationSelection = [at for at in deuterationSelection if at.type.name == 'hydrogen']
    1171
    1172              try:
    1173                  deuterationSelection = Collection(deuterationSelection)
    1174
    1175              except:
    1176                  raise Error('The results of a deuteration selection must be a MMTK Collection.')
    1177
    1178              return deuterationSelection
    1179
    1180 -    def groupSelection(self, universe, selection):
    1181          """Returns a list of MMTK collections where each collection defines a group on which will be applied collectively
    1182          an analysis.
    1183
    1184          @param universe: the universe on which the selection will be performed.
    1185          @type universe: instance of MMTK.Universe
    1186
    1187          @param selection: the selection string that will define the contents of each group.
    1188          @type selection: string
    1189
    1190          @return: a list of MMTK Collection where each collection defines a group..
    1191          @rtype: list
    1192          """
    1193
    1194          if not isinstance(selection, (list, tuple, str)):
    1195              raise Error('Wrong format for the group selection. It must be a list or tuple of MMTK Collections or a group selection string.')
    1196
    1197          if isinstance(selection, (list, tuple)):
    1198              for el in selection:
    1199                  if not isinstance(el, Collection):
    1200                      raise Error('Wrong format for the group selection. It must be a list or tuple of MMTK Collections or a group selection string.')
    1201
    1202              return list(selection)
    1203
    1204          elif selection.lower() == 'all':
    1205              groupSelection = []
    1206              for obj in universe.objectList():
    1207                  if isinstance(obj,(PeptideChain, Protein, NucleotideChain)):
    1208                      for r in obj.residues():
    1209                          groupSelection.append(r.atomList())
    1210
    1211                  else:
    1212                      groupSelection.append(obj.atomList())
    1213
    1214              return [Collection(list(g)) for g in groupSelection]
    1215
    1216          else:
    1217
    1218              groupSelection = []
    1219              # If the selection string contains "filename name-of-file" then it will be a selection from a file.
    1220              if selection.lower().strip()[0:8] == 'filename':
    1221                  groupSelection = self.__parseFileSelectionString(universe, selection, 'group')
    1222
    1223              elif selection.lower().strip()[0:10] == 'expression':
    1224                  groupSelection = self.__parseExpressionSelectionString(universe, selection, 'group')
    1225
    1226              # Otherwise the selection will be from an object.
    1227              else:
    1228                  selectionStringLists = [v.strip() for v in re.split('group\d+', selection) if v != '']
    1229
    1230                  for selString in selectionStringLists:
    1231
    1232                      groupingLevel, sString = [v.strip() for v in re.split('groupinglevel\s+(\w+)\s+', selString) if v != '']
    1233
    1234                      grp = self.__parseObjectSelectionString(universe, sString)
    1235
    1236                      self.__buildGroup(grp, groupSelection, groupingLevel)
    1237
    1238              try:
    1239                  groupSelection = [Collection(list(g)) for g in groupSelection]
    1240
    1241              except:
    1242                  raise Error('Wrong format for the group selection. It must be a list or tuple of MMTK Collections or a group selection string.')
    1243
    1244              return groupSelection
    1245
    1246 -    def __buildGroup(self, g, gSelection, gLevel):
    1247          """
    1248          """
    1249
    1250          gDict= {}
    1251
    1252          if gLevel.lower() == 'chain':
    1253              for at in g:
    1254                  chain = at.parent.parent.parent
    1255                  if gDict.has_key(chain):
    1256                      gDict[chain].append(at)
    1257                  else:
    1258                      gDict[chain] = [at]
    1259
    1260          elif gLevel.lower() == 'residue':
    1261              for at in g:
    1262                  residue = at.parent.parent
    1263                  if gDict.has_key(residue):
    1264                      gDict[residue].append(at)
    1265                  else:
    1266                      gDict[residue] = [at]
    1267
    1268          elif gLevel.lower() in ['molecule', 'cluster', 'nucleicacid', 'protein']:
    1269              for at in g:
    1270                  molecule = at.topLevelChemicalObject()
    1271                  if gDict.has_key(molecule):
    1272                      gDict[molecule].append(at)
    1273                  else:
    1274                      gDict[molecule] = [at]
    1275
    1276          elif gLevel.lower() == 'atom':
    1277              for at in g:
    1278                  gDict[at] = at
    1279
    1280          elif gLevel.lower() == 'amine':
    1281              for at in g:
    1282                  amine = belongToAnAmine(at)
    1283                  if amine is not None:
    1284                      if gDict.has_key(amine):
    1285                          gDict[amine].append(at)
    1286                      else:
    1287                          gDict[amine] = [at]
    1288
    1289          elif gLevel.lower() == 'hydroxy':
    1290              for at in g:
    1291                  hydroxy = belongToAHydroxy(at)
    1292                  if hydroxy is not None:
    1293                      if gDict.has_key(hydroxy):
    1294                          gDict[hydroxy].append(at)
    1295                      else:
    1296                          gDict[hydroxy] = [at]
    1297
    1298          elif gLevel.lower() == 'methyl':
    1299              for at in g:
    1300                  methyl = belongToAMethyl(at)
    1301                  if methyl is not None:
    1302                      if gDict.has_key(methyl):
    1303                          gDict[methyl].append(at)
    1304                      else:
    1305                          gDict[methyl] = [at]
    1306
    1307          elif gLevel.lower() == 'thiol':
    1308              for at in g:
    1309                  thiol = belongToAnHydroxy(at)
    1310                  if thiol is not None:
    1311                      if gDict.has_key(thiol):
    1312                          gDict[thiol].append(at)
    1313                      else:
    1314                          gDict[thiol] = [at]
    1315
    1316          elif gLevel.lower() == 'default':
    1317              for at in g:
    1318                  obj = at.topLevelChemicalObject()
    1319
    1320                  if isinstance(obj, Atom):
    1321                      gDict[at] = at
    1322                  elif isinstance(obj, (AtomCluster,Molecule)):
    1323                      if gDict.has_key(obj):
    1324                          gDict[obj].append(at)
    1325                      else:
    1326                          gDict[obj] = [at]
    1327                  elif isinstance(obj, (NucleotideChain,PeptideChain,Protein)):
    1328                      residue = at.parent.parent
    1329                      if gDict.has_key(residue):
    1330                          gDict[residue].append(at)
    1331                      else:
    1332                          gDict[residue] = [at]
    1333
    1334          for v in gDict.values():
    1335              gSelection.append(set(v))
    1336
    1337 -    def __parseExpressionSelectionString(self, universe, selectionString, selectionMode):
    1338          """Parses a nMOLDYN selection string that starts with 'expression' keyword.
    1339
    1340          @param universe: the universe on which the selection will be performed.
    1341          @type universe: instance of MMTK.Universe
    1342
    1343          @param selectionString: the selection string to parse.
    1344          @type selectionString: string
    1345
    1346          @param selectionMode: a string being one of 'subset', 'deuteration' or 'group' that will specify which kind of
    1347              selection should be performed.
    1348          @type selectionMode: string
    1349
    1350          @return: will depend on |selectionMode|.
    1351          @rtype: Python set for |selectionMode| = 'subset' or 'deuteration' and dict for |selectionMode| = 'group'
    1352          @return: a set of the atoms that match |selectionString| selection string.
    1353          @rtype: set
    1354          """
    1355
    1356          try:
    1357              expression = ' '.join(selectionString.split()[1:])
    1358              exec(expression)
    1359          except:
    1360              raise Error('Could not parse the expression of the selection string.')
    1361
    1362          if not isinstance(selection,list):
    1363              raise Error('The expression of the selection string must provide a list.')
    1364
    1365          if selectionMode in ['subset', 'deuteration']:
    1366
    1367              selection = set(selection)
    1368
    1369          elif selectionMode == 'group':
    1370              pass
    1371
    1372          if not selection:
    1373              raise Error('The selection associated with\n%s\nselection string gave an empty selection.\nSomething might be wrong with it.' % selectionString)
    1374
    1375          return selection
    1376
    1377 -    def __parseFileSelectionString(self, universe, selectionString, selectionMode):
    1378          """Parses a nMOLDYN Selection file (nms file) in order to perform a subset, a deuteration or a group
    1379          selection.
    1380
    1381          @param universe: the universe on which the selection will be performed.
    1382          @type universe: instance of MMTK.Universe
    1383
    1384          @param selectionString: the selection string to parse.
    1385          @type selectionString: string
    1386
    1387          @param selectionMode: a string being one of 'subset', 'deuteration' or 'group' that will specify which kind of
    1388              selection should be performed.
    1389          @type selectionMode: string
    1390
    1391          @return: will depend on |selectionMode|.
    1392          @rtype: Python set for |selectionMode| = 'subset' or 'deuteration' and dict for |selectionMode| = 'group'
    1393          """
    1394
    1395          # First try to open the nms file that contains the informations about the atom to select.
    1396          try:
    1397              filename = selectionString.split()[1].strip()
    1398              nmsFile = open(filename, 'r')
    1399              exec nmsFile
    1400              nmsFile.close()
    1401          except:
    1402              raise Error('Unable to open the file %s.' % filename)
    1403
    1404          try:
    1405              pdb = pdb.strip()
    1406              if not os.path.exists(pdb):
    1407                  raise Error('The pdb file %s does not exist.' % pdb)
    1408
    1409          except NameError:
    1410              raise Error('The "pdb" field was not set in the %s nms file.' % filename)
    1411
    1412          # Sets up the output variable depending on the selected selection type.
    1413          if selectionMode in ['subset', 'deuteration']:
    1414              try:
    1415                  selectedAtoms = [int(v) for v in eval(selectionMode)]
    1416              except:
    1417                  raise Error('The %s selection field must be of the form %s = [int,int,int...]' % (selectionMode,selectionMode))
    1418
    1419              selection = set()
    1420
    1421          elif selectionMode == 'group':
    1422              try:
    1423                  selectedAtoms = [[int(vv) for vv in v] for v in eval(selectionMode)]
    1424
    1425              except :
    1426                      raise Error('The group selection field must be of the form group = [[int,int,...],[int,int,...],[int,int,...]...]')
    1427
    1428              selection = [set()]*len(selectedAtoms)
    1429
    1430          else:
    1431              raise Error('Unknown selection type: %s.' % selectionMode)
    1432
    1433          try:
    1434              atomPos = []
    1435              # Loop over all the atom of the universe.
    1436              for atom in universe.atomList():
    1437                  atomPos.append([atom, [round(v,2) for v in self.trajectory.configuration[0][atom]]])
    1438
    1439              # The PDB file is opened for reading.
    1440              pdb = PDBConfiguration(pdb)
    1441
    1442              if selectionMode in ['subset','deuteration']:
    1443                  # Loop over all the residues found in the pdb file. Here residues are defined in a broader term.
    1444                  # If the system is not a protein, residues referes to a single PDB line.
    1445                  for res in pdb.residues:
    1446                      # Loop over all the atoms of |res|.
    1447                      for at in res.atom_list:
    1448                          # The atom should be selected.
    1449                          if at.properties['serial_number'] in selectedAtoms:
    1450                              pos = [round(v,2) for v in at.position]
    1451                              # Loop over all the atom of the universe.
    1452                              for atom in atomPos:
    1453                                  # Find the one that matches the coordinates of atom |at|.
    1454                                  if atom[1] == pos:
    1455                                      # And add it to the selection.
    1456                                      selection.add(atom[0])
    1457                                      break
    1458
    1459              elif selectionMode == 'group':
    1460                  # Loop over all the residues found in the pdb file. Here residues are definer in a broader term.
    1461                  # If the system is not a protein, residues referes to a single PDB line.
    1462                  for res in pdb.residues:
    1463                      # Loop over all the atoms of |res|.
    1464                      for at in res.atom_list:
    1465                          # Loop over all the selected atom pairs tuples.
    1466                          for p in selectedAtoms:
    1467                              # The pdb atom should be selected.
    1468                              if at.properties['serial_number'] in p:
    1469                                  pos = [round(v,2) for v in at.position]
    1470                                  # Loop over all the atom of the universe.
    1471                                  for atom in atomPos:
    1472                                      # Find the one that matches the coordinates of atom |at|.
    1473                                      if atom[1] == pos:
    1474                                          selection[selectedAtoms.index(p)].add(atom[0])
    1475                                          break
    1476
    1477          except:
    1478              raise Error('Could not parse properly for selection the PDB file associated with %s selection file.' % filename)
    1479
    1480          if not selection:
    1481              raise Error('The selection from file %s gave an empty selection. Something might be wrong with it.' % filename)
    1482
    1483          return selection
    1484
    1485 -    def __parseObjectSelectionString(self, universe, selectionString):
    1486          """Parses a selection string.
    1487
    1488          @param universe: the universe on which the selection will be performed.
    1489          @type universe: instance of MMTK.Universe
    1490
    1491          @param selectionString: the selection strin to parse.
    1492          @type selectionString: string
    1493
    1494          @return: a set of the atoms that match |selectionString| selection string.
    1495          @rtype: set
    1496          """
    1497
    1498          selStr = re.sub('\(',' ( ',selectionString)
    1499          selStr = re.sub('\)',' ) ',selStr)
    1500
    1501          l = [v.strip() for v in re.split('\s+|,',selStr) if v.strip()]
    1502
    1503          processedList = []
    1504          selectionValue = []
    1505
    1506          try:
    1507              while True:
    1508
    1509                  if not l:
    1510                      break
    1511
    1512                  l0 = l[0].lower()
    1513                  if l0 == 'objectname':
    1514                      objectName = l[1]
    1515                      del l[0:2]
    1516                      objectClass = universe.nmoldyncontents[objectName]['objectclass'].lower()
    1517                      continue
    1518
    1519                  elif l0 in universe.nmoldyncontents[objectName].keys():
    1520                      selectionKeyword = l0
    1521                      del l[0]
    1522
    1523                      while True:
    1524
    1525                          ll0 = l[0].lower()
    1526                          if ll0 == 'objectname':
    1527                              raise
    1528
    1529                          elif ll0 in universe.nmoldyncontents[objectName].keys():
    1530                              raise
    1531
    1532                          elif ll0 in ['and','or',')']:
    1533                              if not selectionValue:
    1534                                  raise
    1535                              indexes = self.__retrieveAtomIndexes(universe, objectClass, objectName, selectionKeyword, selectionValue)
    1536                              processedList.append(str(set(indexes)))
    1537                              selectionValue = []
    1538                              break
    1539
    1540                          elif ll0 == '(':
    1541                              raise
    1542
    1543                          else:
    1544                              selectionValue.append(ll0)
    1545                              del l[0]
    1546                              if not l:
    1547                                  indexes = self.__retrieveAtomIndexes(universe, objectClass, objectName, selectionKeyword, selectionValue)
    1548                                  processedList.append(str(set(indexes)))
    1549                                  selectionValue = []
    1550                                  break
    1551
    1552                  elif l0 in ['and','or','(']:
    1553                      processedList.append(l0)
    1554                      del l[0]
    1555                      if not l:
    1556                          raise
    1557
    1558                  elif l0 == ')':
    1559                      processedList.append(l0)
    1560                      del l[0]
    1561
    1562                  else:
    1563                      raise
    1564
    1565          except:
    1566              raise Error('Badly formatted selection string.')
    1567
    1568          processedString = ''.join(processedList)
    1569          processedString = processedString.replace('and', '&')
    1570          processedString = processedString.replace('or', '|')
    1571
    1572          # The MMTK objects corresponding to the selection indexes are returned.
    1573          selection = set([universe.indexToAtoms[ind] for ind in eval(processedString)])
    1574
    1575          if not selection:
    1576              raise Error('The selection associated with\n%s\nselection string gave an empty selection.\nSomething might be wrong with it.' % selectionString)
    1577
    1578          return selection
    1579
    1580 -    def __retrieveAtomIndexes(self, universe, objectClass, objectName, selectionKeyword, selectionValue):
    1581          """Retrieves the MMTK index of the atoms matching |objectClass| MMTK chemical object class,
    1582             |objectName| nMOLDYN name, |selectionKeyword| selection keyword and |selectionValue| selection value.
    1583
    1584             @param universe: the universe on which the selection will be performed.
    1585             @type universe: instance of MMTK.Universe
    1586
    1587             @param objectClass: the MMTK chemical object to match.
    1588             @type objectClass: string
    1589
    1590             @param objectName: the nMOLDYN name to match.
    1591             @type objectName: string
    1592
    1593             @param selectionKeyword: the selection keyword to match.
    1594             @type selectionKeyword: string
    1595
    1596             @param selectionValue: the selection value to match.
    1597             @type selectionValue: string
    1598
    1599             @return: a list of MMTK indexes of the selected atoms.
    1600             @rtype: list
    1601          """
    1602
    1603          if objectClass == 'allclass':
    1604              # The atom parser.
    1605              indexes = self.__allClassParser(universe, selectionKeyword, selectionValue)
    1606
    1607          elif objectClass == 'atom':
    1608              # The atom parser.
    1609              indexes = self.__atomParser(universe, objectName.lower(), selectionKeyword, selectionValue)
    1610
    1611          elif objectClass == 'atomcluster':
    1612              # The atom cluster parser.
    1613              indexes = self.__atomClusterParser(universe, objectName.lower(), selectionKeyword, selectionValue)
    1614
    1615          elif objectClass == 'molecule':
    1616              # The molecule parser.
    1617              indexes = self.__moleculeParser(universe, objectName.lower(), selectionKeyword, selectionValue)
    1618
    1619          elif objectClass == 'nucleotidechain':
    1620              # The nucleotide chain parser.
    1621              indexes = self.__nucleotideChainParser(universe, objectName.lower(), selectionKeyword, selectionValue)
    1622
    1623          elif objectClass == 'peptidechain':
    1624              # The peptide chain parser.
    1625              indexes = self.__peptideChainParser(universe, objectName.lower(), selectionKeyword, selectionValue)
    1626
    1627          elif objectClass == 'protein':
    1628              # The protein parser.
    1629              indexes = self.__proteinParser(universe, objectName.lower(), selectionKeyword, selectionValue)
    1630
    1631          else:
    1632              raise Error('The chemical class %s is not defined in MMTK.' % objectClass)
    1633
    1634          return indexes
    1635
    1636 -    def __allClassParser(self, universe, selectionKeyword, selectionValue):
    1637          """Retrieves the MMTK indexes of the atoms whose nMOLDYN name is '*' matching |selectionKeyword|
    1638             selection keyword and |selectionValue| selection value.
    1639
    1640             @param universe: the universe on which the selection will be performed.
    1641             @type universe: instance of MMTK.Universe
    1642
    1643             @param selectionKeyword: the selection keyword to match.
    1644             @type selectionKeyword: string
    1645
    1646             @param selectionValue: the selection value to match.
    1647             @type selectionValue: string
    1648
    1649             @return: a list of MMTK indexes of the selected atoms.
    1650             @rtype: list
    1651          """
    1652
    1653          selectedObjects = universe.objectList()
    1654
    1655          selection = []
    1656
    1657          if selectionKeyword == 'atomelement':
    1658              if '*' in selectionValue:
    1659                  for obj in selectedObjects:
    1660                      selection.extend([at.index for at in obj.atomList()])
    1661              else:
    1662                  selection = []
    1663                  for v in selectionValue:
    1664                      for obj in selectedObjects:
    1665                          selection.extend([at.index for at in obj.atomList() if at.type.name.strip().lower() == v])
    1666
    1667          return selection
    1668
    1669 -    def __atomParser(self, universe, objectName, selectionKeyword, selectionValue):
    1670          """Retrieves the MMTK index of the atoms whose MMTK chemical object class is 'Atom' matching
    1671             |objectName| nMOLDYN name, |selectionKeyword| selection keyword and |selectionValue| selection value.
    1672
    1673             @param universe: the universe on which the selection will be performed.
    1674             @type universe: instance of MMTK.Universe
    1675
    1676             @param objectName: the nMOLDYN name to match.
    1677             @type objectName: string
    1678
    1679             @param selectionKeyword: the selection keyword to match.
    1680             @type selectionKeyword: string
    1681
    1682             @param selectionValue: the selection value to match.
    1683             @type selectionValue: string
    1684
    1685             @return: a list of MMTK indexes of the selected atoms.
    1686             @rtype: list
    1687          """
    1688
    1689          selection = []
    1690
    1691          if selectionKeyword == 'name':
    1692              # If the '*' selection keyword is in the 'name'-key values then the selection is made of all the objects whose
    1693              # nmoldynname is 'selectionKey'.
    1694              if '*' in selectionValue:
    1695                  selection = [obj.index for obj in universe.objectList() if obj.nmoldynname.strip().lower() == objectName]
    1696
    1697              # Otherwise, only keeps the objects whose nmoldynname is 'selectionKey' AND name matches the selection keyword
    1698              # processed.
    1699              else:
    1700                  selection = []
    1701                  for v in selectionValue:
    1702                      selection.extend([obj.index for obj in universe.objectList() if obj.name.strip().lower() == v])
    1703
    1704          return selection
    1705
    1706 -    def __atomClusterParser(self, universe, objectName, selectionKeyword, selectionValue):
    1707          """Retrieves the MMTK index of the atoms whose MMTK chemical object class is 'AtomCluster' matching
    1708             |objectName| nMOLDYN name, |selectionKeyword| selection keyword and |selectionValue| selection value.
    1709
    1710             @param universe: the universe on which the selection will be performed.
    1711             @type universe: instance of MMTK.Universe
    1712
    1713             @param objectName: the nMOLDYN name to match.
    1714             @type objectName: string
    1715
    1716             @param selectionKeyword: the selection keyword to match.
    1717             @type selectionKeyword: string
    1718
    1719             @param selectionValue: the selection value to match.
    1720             @type selectionValue: string
    1721
    1722             @return: a list of MMTK indexes of the selected atoms.
    1723             @rtype: list
    1724          """
    1725
    1726          selection = []
    1727
    1728          # If the selection dictionnary contains the 'name'-key, then just consider for further selection a subset of objects
    1729          # that will match the selection 'name'-key corresponding values.
    1730          if selectionKeyword == 'name':
    1731              # If the '*' selection keyword is in the 'name'-key values then the subset is made of all the objects whose
    1732              # nmoldynname is 'selectionKey'.
    1733              if '*' in selectionValue:
    1734                  selectedObjects = set([obj for obj in universe.objectList() if obj.nmoldynname.strip().lower() == objectName])
    1735
    1736              # Otherwise, only keeps the objects whose nmoldynname is 'selectionKey' AND name matches the selection keyword
    1737              # processed.
    1738              else:
    1739                  selectedObjects = set()
    1740                  for v in selectionValue:
    1741                      selectedObjects.update([obj for obj in universe.objectList() if obj.name.strip().lower() == v])
    1742
    1743              for obj in selectedObjects:
    1744                  selection.extend([at.index for at in obj.atomList()])
    1745
    1746          else:
    1747              selectedObjects = set([obj for obj in universe.objectList() if obj.nmoldynname.strip().lower() == objectName])
    1748
    1749          if selectionKeyword == 'atomname':
    1750              for v in selectionValue:
    1751                  if v == '*':
    1752                      for obj in selectedObjects:
    1753                          selection.extend([at.index for at in obj.atomList()])
    1754                  else:
    1755                      for obj in selectedObjects:
    1756                          selection.extend([at.index for at in obj.atomList() if at.fullName().strip().lower() == v])
    1757
    1758          elif selectionKeyword == 'atomelement':
    1759              for v in selectionValue:
    1760                  if v == '*':
    1761                      for obj in selectedObjects:
    1762                          selection.extend([at.index for at in obj.atomList()])
    1763                  else:
    1764                      for obj in selectedObjects:
    1765                          selection.extend([at.index for at in obj.atomList() if at.type.name.strip().lower() == v])
    1766
    1767          return selection
    1768
    1769 -    def __moleculeParser(self, universe, objectName, selectionKeyword, selectionValue):
    1770          """Retrieves the MMTK index of the atoms whose MMTK chemical object class is 'Molecule' matching
    1771             |objectName| nMOLDYN name, |selectionKeyword| selection keyword and |selectionValue| selection value.
    1772
    1773             @param universe: the universe on which the selection will be performed.
    1774             @type universe: instance of MMTK.Universe
    1775
    1776             @param objectName: the nMOLDYN name to match.
    1777             @type objectName: string
    1778
    1779             @param selectionKeyword: the selection keyword to match.
    1780             @type selectionKeyword: string
    1781
    1782             @param selectionValue: the selection value to match.
    1783             @type selectionValue: string
    1784
    1785             @return: a list of MMTK indexes of the selected atoms.
    1786             @rtype: list
    1787          """
    1788
    1789          selection = []
    1790          # If the selection dictionnary contains the 'name'-key, then just consider for further selection a subset of objects
    1791          # that will match the selection 'name'-key corresponding values.
    1792          if selectionKeyword == 'name':
    1793              # If the '*' selection keyword is in the 'name'-key values then the subset is made of all the objects whose
    1794              # nmoldynname is 'selectionKey'.
    1795              if '*' in selectionValue:
    1796                  selectedObjects = set([obj for obj in universe.objectList() if obj.nmoldynname.strip().lower() == objectName])
    1797              # Otherwise, only keeps the objects whose nmoldynname is 'selectionKey' AND name matches the selection keyword
    1798              # processed.
    1799
    1800              else:
    1801                  selectedObjects = set()
    1802                  for v in selectionValue:
    1803                      selectedObjects.update([obj for obj in universe.objectList() if obj.name.strip().lower() == v])
    1804              for obj in selectedObjects:
    1805                  selection.extend([at.index for at in obj.atomList()])
    1806
    1807          else:
    1808              selectedObjects = set([obj for obj in universe.objectList() if obj.nmoldynname.strip().lower() == objectName])
    1809
    1810          if selectionKeyword == 'atomname':
    1811              for v in selectionValue:
    1812                  if v == '*':
    1813                      for obj in selectedObjects:
    1814                          selection.extend([at.index for at in obj.atomList()])
    1815                  else:
    1816                      for obj in selectedObjects:
    1817                          selection.extend([at.index for at in obj.atomList() if at.fullName().strip().lower() == v])
    1818
    1819          elif selectionKeyword == 'atomelement':
    1820              for v in selectionValue:
    1821                  if v == '*':
    1822                      for obj in selectedObjects:
    1823                          selection.extend([at.index for at in obj.atomList()])
    1824                  else:
    1825                      for obj in selectedObjects:
    1826                          selection.extend([at.index for at in obj.atomList() if at.type.name.strip().lower() == v])
    1827
    1828          elif selectionKeyword == 'chemfragment':
    1829              for v in selectionValue:
    1830
    1831                  if v == 'amine':
    1832                      for obj in selectedObjects:
    1833                          nitrogens = [at for at in obj.atomList() if at.type.name.strip().lower() == 'nitrogen']
    1834                          for n in nitrogens:
    1835                              neighbours = n.bondedTo()
    1836                              hydrogens = [neigh for neigh in neighbours if neigh.type.name.strip().lower() == 'hydrogen']
    1837                              # The amine.
    1838                              if len(hydrogens) == 2:
    1839                                  selection.extend([n.index] + [h.index for h in hydrogens])
    1840                              # The ammonium.
    1841                              elif (len(hydrogens) == 3) and (obj.numberOfAtoms() == 4):
    1842                                  selection.extend([n.index] + [h.index for h in hydrogens])
    1843
    1844                  elif v == 'hydroxy':
    1845                      for obj in selectedObjects:
    1846                          oxygens = [at for at in obj.atomList() if at.type.name.strip().lower() == 'oxygen']
    1847                          for o in oxygens:
    1848                              neighbours = o.bondedTo()
    1849                              hydrogens = [neigh for neigh in neighbours if neigh.type.name.strip().lower() == 'hydrogen']
    1850                              # The hydroxy.
    1851                              if len(hydrogens) == 1:
    1852                                  selection.extend([o.index] + [h.index for h in hydrogens])
    1853                              # The water.
    1854                              elif (len(hydrogens) == 2) and (obj.numberOfAtoms() == 3):
    1855                                  selection.extend([o.index] + [h.index for h in hydrogens])
    1856
    1857                  elif v == 'methyl':
    1858                      for obj in selectedObjects:
    1859                          carbons = [at for at in obj.atomList() if at.type.name.strip().lower() == 'carbon']
    1860                          for c in carbons:
    1861                              neighbours = c.bondedTo()
    1862                              hydrogens = [neigh for neigh in neighbours if neigh.type.name.strip().lower() == 'hydrogen']
    1863                              # The methyl
    1864                              if len(hydrogens) == 3:
    1865                                  selection.extend([c.index] + [h.index for h in hydrogens])
    1866                              # The methane
    1867                              elif (len(hydrogens) == 4) and (obj.numberOfAtoms() == 5):
    1868                                  selection.extend([c.index] + [h.index for h in hydrogens])
    1869
    1870                  elif v == 'thiol':
    1871                      for obj in selectedObjects:
    1872                          sulphurs = [at for at in obj.atomList() if at.type.name.strip().lower() in ['sulphur','sulfur']]
    1873                          for s in sulphurs:
    1874                              neighbours = s.bondedTo()
    1875                              hydrogens = [neigh for neigh in neighbours if neigh.type.name.strip().lower() == 'hydrogen']
    1876                              # The thiol.
    1877                              if len(hydrogens) == 1:
    1878                                  selection.extend([o.index] + [h.index for h in hydrogens])
    1879                              # The SH2. Quite unusual ...
    1880                              elif (len(hydrogens) == 2) and (obj.numberOfAtoms() == 3):
    1881                                  selection.extend([s.index] + [h.index for h in hydrogens])
    1882
    1883          return selection
    1884
    1885 -    def __nucleotideChainParser(self, universe, objectName, selectionKeyword, selectionValue):
    1886          """Retrieves the MMTK index of the atoms whose MMTK chemical object class is 'NucleotideChain' matching
    1887             |objectName| nMOLDYN name, |selectionKeyword| selection keyword and |selectionValue| selection value.
    1888
    1889             @param universe: the universe on which the selection will be performed.
    1890             @type universe: instance of MMTK.Universe
    1891
    1892             @param objectName: the nMOLDYN name to match.
    1893             @type objectName: string
    1894
    1895             @param selectionKeyword: the selection keyword to match.
    1896             @type selectionKeyword: string
    1897
    1898             @param selectionValue: the selection value to match.
    1899             @type selectionValue: string
    1900
    1901             @return: a list of MMTK indexes of the selected atoms.
    1902             @rtype: list
    1903          """
    1904
    1905          selection = []
    1906
    1907          # If the selection dictionnary contains the 'name'-key, then just consider for further selection a subset of objects
    1908          # that will match the selection 'name'-key corresponding values.
    1909          if selectionKeyword == 'name':
    1910              # If the '*' selection keyword is in the 'name'-key values then the subset is made of all the objects whose
    1911              # nmoldynname is 'selectionKey'.
    1912              if '*' in selectionValue:
    1913                  selectedObjects = set([obj for obj in universe.objectList() if obj.nmoldynname.strip().lower() == objectName])
    1914
    1915              # Otherwise, only keeps the objects whose nmoldynname is 'selectionKey' AND name matches the selection keyword
    1916              # processed.
    1917              else:
    1918                  selectedObjects = set([])
    1919                  for v in selectionValue:
    1920                      selectedObjects.update([obj for obj in universe.objectList() if obj.name.strip().lower() == v])
    1921
    1922              for obj in selectedObjects:
    1923                  selection.extend([at.index for at in obj.atomList()])
    1924
    1925          else:
    1926              selectedObjects = set([obj for obj in universe.objectList() if obj.nmoldynname.strip().lower() == objectName])
    1927
    1928          if selectionKeyword == 'atomname':
    1929              for v in selectionValue:
    1930                  if v == '*':
    1931                      for obj in selectedObjects:
    1932                          selection.extend([at.index for at in obj.atomList()])
    1933                  else:
    1934                      for obj in selectedObjects:
    1935                          selection.extend([at.index for at in obj.atomList() if at.fullName().strip().lower() == v])
    1936
    1937          elif selectionKeyword == 'atomtype':
    1938              for v in selectionValue:
    1939                  if v == '*':
    1940                      for obj in selectedObjects:
    1941                          selection.extend([at.index for at in obj.atomList()])
    1942
    1943                  else:
    1944                      for obj in selectedObjects:
    1945                          selection.extend([at.index for at in obj.atomList() if at.name.strip().lower() == v])
    1946
    1947          elif selectionKeyword == 'atomelement':
    1948              for v in selectionValue:
    1949                  if v == '*':
    1950                      for obj in selectedObjects:
    1951                          selection.extend([at.index for at in obj.atomList()])
    1952                  else:
    1953                      for obj in selectedObjects:
    1954                          selection.extend([at.index for at in obj.atomList() if at.type.name.strip().lower() == v])
    1955
    1956          elif selectionKeyword == 'nuclname':
    1957              for v in selectionValue:
    1958                  if v == '*':
    1959                      for obj in selectedObjects:
    1960                          selection.extend([at.index for at in obj.atomList()])
    1961                  else:
    1962                      for obj in selectedObjects:
    1963                          for nucl in obj.residues():
    1964                              if nucl.fullName().strip().lower() == v:
    1965                                  selection.extend([at.index for at in nucl.atomList()])
    1966
    1967          elif selectionKeyword == 'nucltype':
    1968              for v in selectionValue:
    1969                  if v == '*':
    1970                      for obj in selectedObjects:
    1971                          selection.extend([at.index for at in obj.atomList()])
    1972                  else:
    1973                      for obj in selectedObjects:
    1974                          for nucl in obj.residues():
    1975                              if nucl.symbol.strip().lower() == v:
    1976                                  selection.extend([at.index for at in nucl.atomList()])
    1977
    1978          elif selectionKeyword == 'misc':
    1979              for v in selectionValue:
    1980                  if v == 'bases':
    1981                      for obj in selectedObjects:
    1982                          selection.extend([at.index for at in obj.bases().atomList()])
    1983
    1984                  elif v == 'backbone':
    1985                      for obj in selectedObjects:
    1986                          selection.extend([at.index for at in obj.backbone().atomList()])
    1987
    1988          return selection
    1989
    1990 -    def __peptideChainParser(self, universe, objectName, selectionKeyword, selectionValue):
    1991          """Retrieves the MMTK index of the atoms whose MMTK chemical object class is 'PeptideChain' matching
    1992             |objectName| nMOLDYN name, |selectionKeyword| selection keyword and |selectionValue| selection value.
    1993
    1994             @param universe: the universe on which the selection will be performed.
    1995             @type universe: instance of MMTK.Universe
    1996
    1997             @param objectName: the nMOLDYN name to match.
    1998             @type objectName: string
    1999
    2000             @param selectionKeyword: the selection keyword to match.
    2001             @type selectionKeyword: string
    2002
    2003             @param selectionValue: the selection value to match.
    2004             @type selectionValue: string
    2005
    2006             @return: a list of MMTK indexes of the selected atoms.
    2007             @rtype: list
    2008          """
    2009
    2010          selection = []
    2011          # If the selection dictionnary contains the 'name'-key, then just consider for further selection a subset of objects
    2012          # that will match the selection 'name'-key corresponding values.
    2013          if selectionKeyword == 'name':
    2014
    2015              # If the '*' selection keyword is in the 'name'-key values then the subset is made of all the objects whose
    2016              # nmoldynname is 'selectionKey'.
    2017              if '*' in selectionValue:
    2018                  selectedObjects = set([obj for obj in universe.objectList() if obj.nmoldynname.strip().lower() == objectName])
    2019
    2020              # Otherwise, only keeps the objects whose nmoldynname is 'selectionKey' AND name matches the selection keyword
    2021              # processed.
    2022              else:
    2023                  selectedObjects = set([])
    2024                  for v in selectionValue:
    2025                      selectedObjects.update([obj.index for obj in universe.objectList() if obj.name.strip().lower() == v])
    2026
    2027              for obj in selectedObjects:
    2028                  selection.extend([at.index for at in obj.atomList()])
    2029          else:
    2030              selectedObjects = set([obj for obj in universe.objectList() if obj.nmoldynname.strip().lower() == objectName])
    2031
    2032          if selectionKeyword == 'atomname':
    2033              for v in selectionValue:
    2034                  if v == '*':
    2035                      for obj in selectedObjects:
    2036                          selection.extend([at.index for at in obj.atomList()])
    2037                  else:
    2038                      for obj in selectedObjects:
    2039                          selection.extend([at.index for at in obj.atomList() if at.fullName().strip().lower() == v])
    2040
    2041          elif selectionKeyword == 'atomtype':
    2042              for v in selectionValue:
    2043                  if v == '*':
    2044                      for obj in selectedObjects:
    2045                          selection.extend([at.index for at in obj.atomList()])
    2046                  else:
    2047                      for obj in selectedObjects:
    2048                          selection.extend([at.index for at in obj.atomList() if at.name.strip().lower() == v])
    2049
    2050          elif selectionKeyword == 'atomelement':
    2051              for v in selectionValue:
    2052                  if v == '*':
    2053                      for obj in selectedObjects:
    2054                          selection.extend([at.index for at in obj.atomList()])
    2055                  else:
    2056                      for obj in selectedObjects:
    2057                          selection.extend([at.index for at in obj.atomList() if at.type.name.strip().lower() == v])
    2058
    2059          elif selectionKeyword == 'chemfragment':
    2060              for v in selectionValue:
    2061                  if v == 'amine':
    2062                      for obj in selectedObjects:
    2063                          nitrogens = [at for at in obj.atomList() if at.type.name.strip().lower() == 'nitrogen']
    2064                          for n in nitrogens:
    2065                              neighbours = n.bondedTo()
    2066                              hydrogens = [neigh for neigh in neighbours if neigh.type.name.strip().lower() == 'hydrogen']
    2067                              # The amine.
    2068                              if len(hydrogens) == 2:
    2069                                  selection.extend([n.index] + [h.index for h in hydrogens])
    2070                              # The ammonium.
    2071                              elif (len(hydrogens) == 3) and (obj.numberOfAtoms() == 4):
    2072                                  selection.extend([n.index] + [h.index for h in hydrogens])
    2073
    2074                  elif v == 'c_alphas':
    2075                      for obj in selectedObjects:
    2076                          selection.extend([at.index for at in obj.atomList() if at.name.strip().lower() == 'c_alpha'])
    2077
    2078                  elif v == 'hydroxy':
    2079                      for obj in selectedObjects:
    2080                          oxygens = [at for at in obj.atomList() if at.type.name.strip().lower() == 'oxygen']
    2081                          for o in oxygens:
    2082                              neighbours = o.bondedTo()
    2083                              hydrogens = [neigh for neigh in neighbours if neigh.type.name.strip().lower() == 'hydrogen']
    2084                              # The hydroxy.
    2085                              if len(hydrogens) == 1:
    2086                                  selection.extend([o.index] + [h.index for h in hydrogens])
    2087                              # The water.
    2088                              elif (len(hydrogens) == 2) and (obj.numberOfAtoms() == 3):
    2089                                  selection.extend([o.index] + [h.index for h in hydrogens])
    2090
    2091                  elif v == 'methyl':
    2092                      for obj in selectedObjects:
    2093                          carbons = [at for at in obj.atomList() if at.type.name.strip().lower() == 'carbon']
    2094                          for c in carbons:
    2095                              neighbours = c.bondedTo()
    2096                              hydrogens = [neigh for neigh in neighbours if neigh.type.name.strip().lower() == 'hydrogen']
    2097                              # The methyl
    2098                              if len(hydrogens) == 3:
    2099                                  selection.extend([c.index] + [h.index for h in hydrogens])
    2100                              # The methane
    2101                              elif (len(hydrogens) == 4) and (obj.numberOfAtoms() == 5):
    2102                                  selection.extend([c.index] + [h.index for h in hydrogens])
    2103
    2104                  elif v == 'thiol':
    2105                      for obj in selectedObjects:
    2106                          sulphurs = [at for at in obj.atomList() if at.type.name.strip().lower() in ['sulphur', 'sulfur']]
    2107                          for s in sulphurs:
    2108                              neighbours = s.bondedTo()
    2109                              hydrogens = [neigh for neigh in neighbours if neigh.type.name.strip().lower() == 'hydrogen']
    2110                              # The thiol.
    2111                              if len(hydrogens) == 1:
    2112                                  selection.extend([o.index] + [h.index for h in hydrogens])
    2113                              # The SH2. Quite unusual ...
    2114                              elif (len(hydrogens) == 2) and (obj.numberOfAtoms() == 3):
    2115                                  selection.extend([s.index] + [h.index for h in hydrogens])
    2116
    2117          elif selectionKeyword == 'resname':
    2118              for v in selectionValue:
    2119                  if v == '*':
    2120                      for obj in selectedObjects:
    2121                          selection.extend([at.index for at in obj.atomList()])
    2122                  else:
    2123                      for obj in selectedObjects:
    2124                          for res in obj.residues():
    2125                              if res.fullName().strip().lower() == v:
    2126                                  selection.extend([at.index for at in res.atomList()])
    2127
    2128          elif selectionKeyword == 'restype':
    2129              for v in selectionValue:
    2130                  if v == '*':
    2131                      for obj in selectedObjects:
    2132                          selection.extend([at.index for at in obj.atomList()])
    2133                  else:
    2134                      for obj in selectedObjects:
    2135                          for res in obj.residues():
    2136                              if res.symbol.strip().lower() == v:
    2137                                  selection.extend([at.index for at in res.atomList()])
    2138
    2139          elif selectionKeyword == 'resclass':
    2140              for v in selectionValue:
    2141                  for obj in selectedObjects:
    2142                      for res in obj.residues():
    2143                          if res.symbol.strip().lower() in residusChemFamily[v]:
    2144                              selection.extend([at.index for at in res.atomList()])
    2145
    2146          elif selectionKeyword == 'misc':
    2147              for v in selectionValue:
    2148                  if v == 'sidechains':
    2149                      for obj in selectedObjects:
    2150                          selection.extend([at.index for at in obj.sidechains().atomList()])
    2151
    2152                  elif v == 'backbone':
    2153                      for obj in selectedObjects:
    2154                          selection.extend([at.index for at in obj.backbone().atomList()])
    2155
    2156          return selection
    2157
    2158 -    def __proteinParser(self, universe, objectName, selectionKeyword, selectionValue):
    2159          """Retrieves the MMTK index of the atoms whose MMTK chemical object class is 'Protein' matching
    2160             |objectName| nMOLDYN name, |selectionKeyword| selection keyword and |selectionValue| selection value.
    2161
    2162             @param universe: the universe on which the selection will be performed.
    2163             @type universe: instance of MMTK.Universe
    2164
    2165             @param objectName: the nMOLDYN name to match.
    2166             @type objectName: string
    2167
    2168             @param selectionKeyword: the selection keyword to match.
    2169             @type selectionKeyword: string
    2170
    2171             @param selectionValue: the selection value to match.
    2172             @type selectionValue: string
    2173
    2174             @return: a list of MMTK indexes of the selected atoms.
    2175             @rtype: list
    2176          """
    2177
    2178          selection = []
    2179          # If the selection dictionnary contains the 'name'-key, then just consider for further selection a subset of objects
    2180          # that will match the selection 'name'-key corresponding values.
    2181          if selectionKeyword == 'name':
    2182              # If the '*' selection keyword is in the 'name'-key values then the subset is made of all the objects whose
    2183              # nmoldynname is 'selectionKey'.
    2184              if '*' in selectionValue:
    2185                  selectedObjects = set([obj for obj in universe.objectList() if obj.nmoldynname.strip().lower() == objectName])
    2186              # Otherwise, only keeps the objects whose nmoldynname is 'selectionKey' AND name matches the selection keyword
    2187              # processed.
    2188              else:
    2189                  selectedObjects = set([])
    2190                  for v in selectionValue:
    2191                      selectedObjects.update([obj.index for obj in universe.objectList() if obj.name.strip().lower() == v])
    2192
    2193              for obj in selectedObjects:
    2194                  selection.extend([at.index for at in obj.atomList()])
    2195          else:
    2196              selectedObjects = set([obj for obj in universe.objectList() if obj.nmoldynname.strip().lower() == objectName])
    2197
    2198
    2199          if selectionKeyword == 'atomname':
    2200              for v in selectionValue:
    2201                  if v == '*':
    2202                      for obj in selectedObjects:
    2203                          selection.extend([at.index for at in obj.atomList()])
    2204                  else:
    2205                      for obj in selectedObjects:
    2206                          selection.extend([at.index for at in obj.atomList() if at.fullName().strip().lower() == v])
    2207
    2208          elif selectionKeyword == 'atomtype':
    2209              for v in selectionValue:
    2210                  if v == '*':
    2211                      for obj in selectedObjects:
    2212                          selection.extend([at.index for at in obj.atomList()])
    2213                  else:
    2214                      for obj in selectedObjects:
    2215                          selection.extend([at.index for at in obj.atomList() if at.name.strip().lower() == v])
    2216
    2217          elif selectionKeyword == 'atomelement':
    2218              for v in selectionValue:
    2219                  if v == '*':
    2220                      for obj in selectedObjects:
    2221                          selection.extend([at.index for at in obj.atomList()])
    2222                  else:
    2223                      for obj in selectedObjects:
    2224                          selection.extend([at.index for at in obj.atomList() if at.type.name.strip().lower() == v])
    2225
    2226          elif selectionKeyword == 'chemfragment':
    2227              for v in selectionValue:
    2228
    2229                  if v == 'amine':
    2230                      for obj in selectedObjects:
    2231                          nitrogens = [at for at in obj.atomList() if at.type.name.strip().lower() == 'nitrogen']
    2232                          for n in nitrogens:
    2233                              neighbours = n.bondedTo()
    2234                              hydrogens = [neigh for neigh in neighbours if neigh.type.name.strip().lower() == 'hydrogen']
    2235                              # The amine.
    2236                              if len(hydrogens) == 2:
    2237                                  selection.extend([n.index] + [h.index for h in hydrogens])
    2238                              # The ammonium.
    2239                              elif (len(hydrogens) == 3) and (obj.numberOfAtoms() == 4):
    2240                                  selection.extend([n.index] + [h.index for h in hydrogens])
    2241
    2242                  elif v == 'c_alphas':
    2243                      for obj in selectedObjects:
    2244                          selection.extend([at.index for at in obj.atomList() if at.name.strip().lower() == 'c_alpha'])
    2245
    2246                  elif v == 'hydroxy':
    2247                      for obj in selectedObjects:
    2248                          oxygens = [at for at in obj.atomList() if at.type.name.strip().lower() == 'oxygen']
    2249                          for o in oxygens:
    2250                              neighbours = o.bondedTo()
    2251                              hydrogens = [neigh for neigh in neighbours if neigh.type.name.strip().lower() == 'hydrogen']
    2252                              # The hydroxy.
    2253                              if len(hydrogens) == 1:
    2254                                  selection.extend([o.index] + [h.index for h in hydrogens])
    2255                              # The water.
    2256                              elif (len(hydrogens) == 2) and (obj.numberOfAtoms() == 3):
    2257                                  selection.extend([o.index] + [h.index for h in hydrogens])
    2258
    2259                  elif v == 'methyl':
    2260                      for obj in selectedObjects:
    2261                          carbons = [at for at in obj.atomList() if at.type.name.strip().lower() == 'carbon']
    2262                          for c in carbons:
    2263                              neighbours = c.bondedTo()
    2264                              hydrogens = [neigh for neigh in neighbours if neigh.type.name.strip().lower() == 'hydrogen']
    2265                              # The methyl
    2266                              if len(hydrogens) == 3:
    2267                                  selection.extend([c.index] + [h.index for h in hydrogens])
    2268                              # The methane
    2269                              elif (len(hydrogens) == 4) and (obj.numberOfAtoms() == 5):
    2270                                  selection.extend([c.index] + [h.index for h in hydrogens])
    2271
    2272                  elif v == 'thiol':
    2273                      for obj in selectedObjects:
    2274                          sulphurs = [at for at in obj.atomList() if at.type.name.strip().lower() in ['sulphur','sulfur']]
    2275                          for s in sulphurs:
    2276                              neighbours = s.bondedTo()
    2277                              hydrogens = [neigh for neigh in neighbours if neigh.type.name.strip().lower() == 'hydrogen']
    2278                              # The thiol.
    2279                              if len(hydrogens) == 1:
    2280                                  selection.extend([o.index] + [h.index for h in hydrogens])
    2281                              # The SH2. Quite unusual ...
    2282                              elif (len(hydrogens) == 2) and (obj.numberOfAtoms() == 3):
    2283                                  selection.extend([s.index] + [h.index for h in hydrogens])
    2284
    2285          elif selectionKeyword == 'resname':
    2286              for v in selectionValue:
    2287                  if v == '*':
    2288                      for obj in selectedObjects:
    2289                          selection.extend([at.index for at in obj.atomList()])
    2290                  else:
    2291                      for obj in selectedObjects:
    2292                          for res in obj.residues():
    2293                              if res.fullName().strip().lower() == v:
    2294                                  selection.extend([at.index for at in res.atomList()])
    2295
    2296          elif selectionKeyword == 'restype':
    2297              for v in selectionValue:
    2298                  if v == '*':
    2299                      for obj in selectedObjects:
    2300                          selection.extend([at.index for at in obj.atomList()])
    2301                  else:
    2302                      for obj in selectedObjects:
    2303                          for res in obj.residues():
    2304                              if res.symbol.strip().lower() == v:
    2305                                  selection.extend([at.index for at in res.atomList()])
    2306
    2307          elif selectionKeyword == 'resclass':
    2308              for v in selectionValue:
    2309                  for obj in selectedObjects:
    2310                      for res in obj.residues():
    2311                          if res.symbol.strip().lower() in residusChemFamily[v]:
    2312                              selection.extend([at.index for at in res.atomList()])
    2313
    2314          elif selectionKeyword == 'chainname':
    2315              for v in selectionValue:
    2316                  if v == '*':
    2317                      for obj in selectedObjects:
    2318                          selection.extend([at.index for at in obj.atomList()])
    2319                  else:
    2320                      for obj in selectedObjects:
    2321                          for chain in obj:
    2322                              if chain.name.strip().lower() == v:
    2323                                  selection.extend([at.index for at in chain.atomList()])
    2324
    2325          elif selectionKeyword == 'misc':
    2326              for v in selectionValue:
    2327                  if v == 'sidechains':
    2328                      for obj in selectedObjects:
    2329                          selection.extend([at.index for at in obj.sidechains().atomList()])
    2330
    2331                  elif v == 'backbone':
    2332                      for obj in selectedObjects:
    2333                          selection.extend([at.index for at in obj.backbone().atomList()])
    2334
    2335          return selection
    2336
    2337 -def setUniverseContents(universe):
    2338      """Sets the contents of each object found in the universe.
    2339
    2340      @param universe: the MMTK universe to look in.
    2341      @type universe: a instance of MMTK.Universe.
    2342      """
    2343
    2344      if hasattr(universe, 'nmoldyncontents'):
    2345          return
    2346
    2347      isotope = {}
    2348      universe.indexToAtoms = {}
    2349      # Loop over the object found in the universe.
    2350      for obj in universe.objectList():
    2351          if not isinstance(obj,(Atom, AtomCluster, Molecule, NucleotideChain, PeptideChain, Protein)):
    2352              continue
    2353
    2354          for at in obj.atomList():
    2355              try:
    2356                  universe.indexToAtoms[at.index] = at
    2357              except:
    2358                  pass
    2359
    2360          try:
    2361              # Check if the object is referenced in the MMTK database. In that case, it has a name.
    2362              if obj.type.name:
    2363                  obj.nmoldynname = obj.type.name
    2364              else:
    2365                  raise
    2366
    2367          except:
    2368
    2369              if isinstance(obj, Atom):
    2370                  obj.nmoldynname = 'A'
    2371              elif isinstance(obj, AtomCluster):
    2372                  obj.nmoldynname = 'AC'
    2373              elif isinstance(obj, Molecule):
    2374                  obj.nmoldynname = 'M'
    2375              elif isinstance(obj, NucleotideChain):
    2376                  obj.nmoldynname = 'NC'
    2377              elif isinstance(obj, PeptideChain):
    2378                  obj.nmoldynname = 'PC'
    2379              elif isinstance(obj, Protein):
    2380                  obj.nmoldynname = 'P'
    2381
    2382              obj.nmoldynname += str(obj.numberOfAtoms())
    2383
    2384          # Case of the isotopes.
    2385          obj.nmoldynnamelist = sorted([at.name for at in obj.atomList()])
    2386          if not isotope.has_key(obj.nmoldynname):
    2387              isotope[obj.nmoldynname] = [obj.nmoldynnamelist]
    2388
    2389          else:
    2390              if not obj.nmoldynnamelist in isotope[obj.nmoldynname]:
    2391                  isotope[obj.nmoldynname].append(obj.nmoldynnamelist)
    2392
    2393      universe.nmoldyncontents = {}
    2394
    2395      # Loop over the object found in the universe.
    2396      for obj in universe.objectList():
    2397
    2398          if not isinstance(obj,(Atom, AtomCluster, Molecule, NucleotideChain, PeptideChain, Protein)):
    2399              continue
    2400
    2401          if len(isotope[obj.nmoldynname]) > 1:
    2402              LogMessage('info','Isotope found for %s' % obj.nmoldynname,['file','console'])
    2403              iso = str(isotope[obj.nmoldynname].index(obj.nmoldynnamelist) + 1)
    2404              obj.nmoldynname += '_iso' + iso
    2405              obj.writeToFile(os.path.join(PREFERENCES.outputfile_path,obj.nmoldynname + '.pdb'), format = 'pdb')
    2406
    2407          delattr(obj, 'nmoldynnamelist')
    2408
    2409          if not obj.name:
    2410              obj.name = obj.nmoldynname
    2411
    2412          # If the objet type name is already a key of the |universeContents| dictionnary then just add the name of
    2413          # the object to the 'name' set and increment the 'number' counter.
    2414          if universe.nmoldyncontents.has_key(obj.nmoldynname):
    2415              universe.nmoldyncontents[obj.nmoldynname]['name'].add(obj.name)
    2416              universe.nmoldyncontents[obj.nmoldynname]['number'] += 1
    2417
    2418          else:
    2419              universe.nmoldyncontents[obj.nmoldynname] = {'name' : set(['*', obj.name]) , 'number' : 1}
    2420
    2421              # Case of an object of |Atom| class.
    2422              # No additional selection fields:
    2423              if isinstance(obj, Atom):
    2424                  universe.nmoldyncontents[obj.nmoldynname]['objectclass'] = 'Atom'
    2425                  universe.nmoldyncontents[obj.nmoldynname]['groupinglevel'] = [0, ['atom',]]
    2426
    2427              # Case of an object of |AtomCluster| class.
    2428              # Two additional selection fields:
    2429              #     -atomname:    the full name of the atom.
    2430              #     -atomelement: the element of the atom.
    2431              elif isinstance(obj, AtomCluster):
    2432                  universe.nmoldyncontents[obj.nmoldynname]['objectclass'] = 'AtomCluster'
    2433                  universe.nmoldyncontents[obj.nmoldynname]['groupinglevel'] = [1,  ['atom','cluster']]
    2434                  universe.nmoldyncontents[obj.nmoldynname]['atomname'] = set(['*'])
    2435                  universe.nmoldyncontents[obj.nmoldynname]['atomelement'] = set(['*'])
    2436                  for atom in obj.atomList():
    2437                      universe.nmoldyncontents[obj.nmoldynname]['atomname'].add(atom.fullName())
    2438                      universe.nmoldyncontents[obj.nmoldynname]['atomelement'].add(atom.type.name)
    2439
    2440                  universe.nmoldyncontents[obj.nmoldynname]['atomname'] = sorted(universe.nmoldyncontents[obj.nmoldynname]['atomname'])
    2441                  universe.nmoldyncontents[obj.nmoldynname]['atomelement'] = sorted(universe.nmoldyncontents[obj.nmoldynname]['atomelement'])
    2442
    2443              # Case of an object of |Molecule| class.
    2444              # Three additional selection fields:
    2445              #     -atomname:     the full name of the atom.
    2446              #     -atomelement:  the element of the atom.
    2447              #     -chemfragment: the chemical class of the atom.
    2448              elif isinstance(obj, Molecule):
    2449                  universe.nmoldyncontents[obj.nmoldynname]['objectclass'] = 'Molecule'
    2450                  universe.nmoldyncontents[obj.nmoldynname]['groupinglevel'] = [5, ['atom','amine', 'hydroxy', 'methyl', 'thiol', 'molecule']]
    2451                  universe.nmoldyncontents[obj.nmoldynname]['atomname'] = set(['*'])
    2452                  universe.nmoldyncontents[obj.nmoldynname]['atomelement'] = set(['*'])
    2453                  universe.nmoldyncontents[obj.nmoldynname]['chemfragment'] = ['amine', 'hydroxy', 'methyl', 'thiol']
    2454
    2455                  for atom in obj.atomList():
    2456                      universe.nmoldyncontents[obj.nmoldynname]['atomname'].add(atom.fullName())
    2457                      universe.nmoldyncontents[obj.nmoldynname]['atomelement'].add(atom.type.name)
    2458
    2459                  universe.nmoldyncontents[obj.nmoldynname]['atomname'] = sorted(universe.nmoldyncontents[obj.nmoldynname]['atomname'])
    2460                  universe.nmoldyncontents[obj.nmoldynname]['atomelement'] = sorted(universe.nmoldyncontents[obj.nmoldynname]['atomelement'])
    2461
    2462              # Case of an object of |NucleotideChain| class.
    2463              # Six additional selection fields:
    2464              #     -atomname:     the full name of the atom.
    2465              #     -atomtype:     the type of the atom.
    2466              #     -atomelement:  the element of the atom.
    2467              #     -nuclname:     the name of the nucleotide.
    2468              #     -nucltype:     the type of the nucleotide.
    2469              #     -chemfragment: some miscellaneous selection keywords.
    2470              elif isinstance(obj, NucleotideChain):
    2471                  universe.nmoldyncontents[obj.nmoldynname]['objectclass'] = 'NucleotideChain'
    2472                  universe.nmoldyncontents[obj.nmoldynname]['groupinglevel'] = [2, ['atom', 'amine', 'residue', 'nucleicacid']]
    2473                  universe.nmoldyncontents[obj.nmoldynname]['atomname'] = ['*']
    2474                  universe.nmoldyncontents[obj.nmoldynname]['atomtype'] = set(['*'])
    2475                  universe.nmoldyncontents[obj.nmoldynname]['atomelement'] = set(['*'])
    2476                  universe.nmoldyncontents[obj.nmoldynname]['nuclname'] = ['*']
    2477                  universe.nmoldyncontents[obj.nmoldynname]['nucltype'] = set(['*'])
    2478                  universe.nmoldyncontents[obj.nmoldynname]['misc'] = ['backbone', 'bases']
    2479
    2480                  # Loop over the residues of the chain.
    2481                  for nucl in obj.residues():
    2482                      universe.nmoldyncontents[obj.nmoldynname]['nuclname'].append(nucl.fullName())
    2483                      universe.nmoldyncontents[obj.nmoldynname]['nucltype'].add(nucl.symbol)
    2484
    2485                      # Loop over the atoms of the residue.
    2486                      for atom in nucl.atomList():
    2487                          universe.nmoldyncontents[obj.nmoldynname]['atomname'].append(atom.fullName())
    2488                          universe.nmoldyncontents[obj.nmoldynname]['atomtype'].add(atom.name)
    2489                          universe.nmoldyncontents[obj.nmoldynname]['atomelement'].add(atom.type.name)
    2490
    2491                  universe.nmoldyncontents[obj.nmoldynname]['atomtype'] = sorted(universe.nmoldyncontents[obj.nmoldynname]['atomtype'])
    2492                  universe.nmoldyncontents[obj.nmoldynname]['atomelement'] = sorted(universe.nmoldyncontents[obj.nmoldynname]['atomelement'])
    2493
    2494                  universe.nmoldyncontents[obj.nmoldynname]['nucltype'] = sorted(universe.nmoldyncontents[obj.nmoldynname]['nucltype'])
    2495
    2496              # Case of an object of |PeptideChain| class.
    2497              # Eight additional selection fields:
    2498              #     -atomname:     the full name of the atom.
    2499              #     -atomtype:     the type of the atom.
    2500              #     -atomelement:  the element of the atom.
    2501              #     -chemfragment: the chemical class of the atom.
    2502              #     -resname:      the name of the residue.
    2503              #     -restype:      the type of the residue.
    2504              #     -resclass:     the chemical class of the residue.
    2505              #     -chemfragment: some miscellaneous selection keywords.
    2506              elif isinstance(obj, PeptideChain):
    2507                  universe.nmoldyncontents[obj.nmoldynname]['objectclass'] = 'PeptideChain'
    2508                  universe.nmoldyncontents[obj.nmoldynname]['groupinglevel'] = [5, ['atom', 'amine', 'hydroxy', 'methyl', 'thiol', 'residue','chain']]
    2509                  universe.nmoldyncontents[obj.nmoldynname]['atomname'] = ['*']
    2510                  universe.nmoldyncontents[obj.nmoldynname]['atomtype'] = set(['*'])
    2511                  universe.nmoldyncontents[obj.nmoldynname]['atomelement'] = set(['*'])
    2512                  universe.nmoldyncontents[obj.nmoldynname]['chemfragment'] = ['amine', 'c_alphas', 'hydroxy', 'methyl', 'thiol']
    2513                  universe.nmoldyncontents[obj.nmoldynname]['resname'] = ['*']
    2514                  universe.nmoldyncontents[obj.nmoldynname]['restype'] = set(['*'])
    2515                  universe.nmoldyncontents[obj.nmoldynname]['resclass'] = sorted(residusChemFamily.keys())
    2516                  universe.nmoldyncontents[obj.nmoldynname]['misc'] = ['backbone', 'sidechains']
    2517
    2518                  # Loop over the residues of the chain.
    2519                  for residue in obj.residues():
    2520                      universe.nmoldyncontents[obj.nmoldynname]['resname'].append(residue.fullName())
    2521                      universe.nmoldyncontents[obj.nmoldynname]['restype'].add(residue.symbol)
    2522
    2523                      # Loop over the atoms of the residue.
    2524                      for atom in residue.atomList():
    2525                          universe.nmoldyncontents[obj.nmoldynname]['atomname'].append(atom.fullName())
    2526                          universe.nmoldyncontents[obj.nmoldynname]['atomtype'].add(atom.name)
    2527                          universe.nmoldyncontents[obj.nmoldynname]['atomelement'].add(atom.type.name)
    2528
    2529                  universe.nmoldyncontents[obj.nmoldynname]['atomtype'] = sorted(universe.nmoldyncontents[obj.nmoldynname]['atomtype'])
    2530                  universe.nmoldyncontents[obj.nmoldynname]['atomelement'] = sorted(universe.nmoldyncontents[obj.nmoldynname]['atomelement'])
    2531
    2532                  universe.nmoldyncontents[obj.nmoldynname]['restype'] = sorted(universe.nmoldyncontents[obj.nmoldynname]['restype'])
    2533
    2534              # Case of an object of |Protein| class.
    2535              # Nine additional selection fields:
    2536              #     -atomname:     the full name of the atom.
    2537              #     -atomtype:     the type of the atom.
    2538              #     -atomelement:  the element of the atom.
    2539              #     -chemfragment: the chemical class of the atom.
    2540              #     -resname:      the name of the residue.
    2541              #     -restype:      the type of the residue.
    2542              #     -resclass:     the chemical class of the residue.
    2543              #     -chainname:    the name of the chain.
    2544              #     -chemfragment: some miscellaneous selection keywords.
    2545              elif isinstance(obj, Protein):
    2546                  universe.nmoldyncontents[obj.nmoldynname]['objectclass'] = 'Protein'
    2547                  universe.nmoldyncontents[obj.nmoldynname]['groupinglevel'] = [5, ['atom', 'amine', 'hydroxy', 'methyl', 'thiol', 'residue','chain', 'protein']]
    2548                  universe.nmoldyncontents[obj.nmoldynname]['atomname'] = ['*']
    2549                  universe.nmoldyncontents[obj.nmoldynname]['atomtype'] = set(['*'])
    2550                  universe.nmoldyncontents[obj.nmoldynname]['atomelement'] = set(['*'])
    2551                  universe.nmoldyncontents[obj.nmoldynname]['chemfragment'] = ['amine', 'c_alphas', 'hydroxy', 'methyl', 'thiol']
    2552                  universe.nmoldyncontents[obj.nmoldynname]['resname'] = ['*']
    2553                  universe.nmoldyncontents[obj.nmoldynname]['restype'] = set(['*'])
    2554                  universe.nmoldyncontents[obj.nmoldynname]['resclass'] = sorted(residusChemFamily.keys())
    2555                  universe.nmoldyncontents[obj.nmoldynname]['chainname'] = ['*']
    2556                  universe.nmoldyncontents[obj.nmoldynname]['misc'] = ['backbone', 'sidechains']
    2557
    2558                  # Loop over the chains of the protein.
    2559                  for chain in obj:
    2560                      # Append the chain name.
    2561                      universe.nmoldyncontents[obj.nmoldynname]['chainname'].append(chain.fullName())
    2562
    2563                      # Loop over the residues of the chain.
    2564                      for residue in chain.residues():
    2565                          universe.nmoldyncontents[obj.nmoldynname]['resname'].append(residue.fullName())
    2566                          universe.nmoldyncontents[obj.nmoldynname]['restype'].add(residue.symbol)
    2567
    2568                          # Loop over the atoms of the residue.
    2569                          for atom in residue.atomList():
    2570                              universe.nmoldyncontents[obj.nmoldynname]['atomname'].append(atom.fullName())
    2571                              universe.nmoldyncontents[obj.nmoldynname]['atomtype'].add(atom.name)
    2572                              universe.nmoldyncontents[obj.nmoldynname]['atomelement'].add(atom.type.name)
    2573
    2574                  universe.nmoldyncontents[obj.nmoldynname]['atomtype'] = sorted(universe.nmoldyncontents[obj.nmoldynname]['atomtype'])
    2575                  universe.nmoldyncontents[obj.nmoldynname]['atomelement'] = sorted(universe.nmoldyncontents[obj.nmoldynname]['atomelement'])
    2576
    2577                  universe.nmoldyncontents[obj.nmoldynname]['restype'] = sorted(universe.nmoldyncontents[obj.nmoldynname]['restype'])
    2578
    2579      # The 'name' set value is converted to a sorted list.
    2580      for objTypeName in universe.nmoldyncontents.keys():
    2581          universe.nmoldyncontents[objTypeName]['name'] = sorted(universe.nmoldyncontents[objTypeName]['name'])
    2582
    2583      # Special case of * key used to selected all the nmoldyn type found in the universe.
    2584      universe.nmoldyncontents['*'] = {'atomelement' : set(['*']),\
    2585                                      'number' : universe.numberOfAtoms(),\
    2586                                      'objectclass' : 'AllClass',\
    2587                                      'groupinglevel' : [0,['default',]]}
    2588      for atom in universe.atomList():
    2589          universe.nmoldyncontents['*']['atomelement'].add(atom.type.name)
    2590      universe.nmoldyncontents['*']['atomelement'] = sorted(universe.nmoldyncontents['*']['atomelement'])
    2591
    2592 -class QVectors:
    2593      """Generates a set of QVectors within a given shell.
    2594      """
    2595 -    def __init__(self, universe, generator, qRadii, dq, qVectorsPerShell, qVectorsDirection = None):
    2596          """The constructor.
    2597
    2598          @param universe: the MMTK universe used to define the reciprocal space.
    2599          @type universe: a MMTK.Universe subclass object
    2600
    2601          @param generator: a string being one of '3d isotropic', '2d isotropic' or 'anistropic' the way the q-vectors
    2602              should be generated.
    2603          @type generator: string
    2604
    2605          @param qRadii: a list of floats specifying the radii of the shell in which the q vectors have to be
    2606              generated.
    2607          @type qRadii: list
    2608
    2609          @param dq: a float specifying the width of a qhsell defined as [|qRadius| - dq/2,|qRadius|+dq/2].
    2610          @type dq: float
    2611
    2612          @param qVectorsPerShell: an integer specifying the number of q-vectors to generate for each shell.
    2613          @type qVectorsPerShell: integer
    2614
    2615          @param qVectorsDirection: a list of Scientific.Geometry.Vector objects specifying the directions along which
    2616              the q-vectors should be generated. If None, the q-vectors generation will be isotropic.
    2617          @type qVectorsDirection: list
    2618          """
    2619
    2620          self.universe = universe
    2621          self.generator = generator
    2622          self.dq = dq
    2623          self.qVectorsPerShell = qVectorsPerShell
    2624          self.qVectorsDirection = qVectorsDirection
    2625
    2626          if self.generator == '3d isotropic':
    2627              self.qVectorsDirection = None
    2628
    2629          else:
    2630
    2631              if self.qVectorsDirection is None:
    2632                  raise Error('Directions must be given for %s q vectors generation.' % self.generator)
    2633
    2634              if self.generator == '2d isotropic':
    2635                  if len(self.qVectorsDirection) != 2:
    2636                      raise Error('Two vectors must be given to generate q vectors in a plane.')
    2637
    2638                  v1, v2 = self.qVectorsDirection
    2639                  if abs(v1.angle(v2)) < 1.e-10:
    2640                      raise Error('The two vectors are colinear.')
    2641
    2642              self.qVectorsDirection = qVectorsDirection
    2643
    2644          self.qRadii = []
    2645          self.qMin = []
    2646          self.qMax = []
    2647          # Loop over all the sphere radius on which the q vectors have to be generated.
    2648          for qRadius in sorted(qRadii):
    2649              # The case q = 0 produces untrustable results. So skip it.
    2650              if qRadius <= 0.0:
    2651                  LogMessage('warning','No Q vector could be generated for shell with radius Q = %s.' % qRadius,['file','console'])
    2652                  continue
    2653              else:
    2654                  self.qRadii.append(qRadius)
    2655                  self.qMin.append(max(qRadius - 0.5*self.dq,0.0))
    2656                  self.qMax.append(qRadius + 0.5*self.dq)
    2657
    2658          # Get the reciprocal basis of the elementary cell. If universe is non-periodic, basis = None
    2659          self.reciprocalBasis = universe.reciprocalBasisVectors()
    2660
    2661          # Non-periodic universe.
    2662          if self.reciprocalBasis is None:
    2663              self.directMatrix = None
    2664              self.reciprocalCellVolume = None
    2665
    2666          # Periodic universe.
    2667          else:
    2668              self.reciprocalBasis = [2.0*N.pi*v for v in self.reciprocalBasis]
    2669              self.directMatrix = N.array([N.array(v) for v in self.universe.basisVectors()])
    2670              self.directMatrix = Tensor(self.directMatrix)/(2.*N.pi)
    2671              self.reciprocalCellVolume = self.reciprocalBasis[0]*(self.reciprocalBasis[1].cross(self.reciprocalBasis[2]))
    2672
    2673          self.qVectors = {}
    2674
    2675          # Loop over all the sphere radius on which the q vectors have to be generated.
    2676          for index in range(len(self.qRadii)):
    2677              qRadius = self.qRadii[index]
    2678              self.qVectors[qRadius] = []
    2679              qMin = self.qMin[index]
    2680              qMax = self.qMax[index]
    2681
    2682              if self.generator in ['2d isotropic', '3d isotropic']:
    2683                  self.__isotropicGeneration(qRadius, qMin, qMax)
    2684
    2685              elif self.generator == '2d isotropic':
    2686                  self.__planarIsotropicGeneration(qRadius, qMin, qMax)
    2687
    2688              elif self.generator == 'anisotropic':
    2689                  self.__anisotropicGeneration(qRadius, qMin, qMax)
    2690
    2691              if len(self.qVectors[qRadius]) == 0:
    2692                  LogMessage('warning','No Q vector could be generated for shell with radius Q = %s.' % qRadius,['file','console'])
    2693                  del self.qVectors[qRadius]
    2694                  continue
    2695
    2696          if not self.qVectors:
    2697              raise Error('No Q Vectors generated. Please check your settings.')
    2698
    2699          self.statistics = N.zeros((len(self.qVectors), 8), typecode = N.Int)
    2700          comp = 0
    2701          for qRadius in sorted(self.qVectors.keys()):
    2702              qVectors = self.qVectors[qRadius]
    2703              for qVect in qVectors:
    2704                  if qVect[0] > 0:
    2705                      if qVect[1] > 0:
    2706                          # X+Y+Z+
    2707                          if qVect[2] > 0:
    2708                              self.statistics[comp, 0] += 1
    2709                          # X+Y+Z-
    2710                          else:
    2711                              self.statistics[comp, 1] += 1
    2712                      else:
    2713                          # X+Y-Z+
    2714                          if qVect[2] > 0:
    2715                              self.statistics[comp, 2] += 1
    2716                          # X+Y-Z-
    2717                          else:
    2718                              self.statistics[comp, 3] += 1
    2719                  else:
    2720                      if qVect[1] > 0:
    2721                          # X-Y+Z+
    2722                          if qVect[2] > 0:
    2723                              self.statistics[comp, 4] += 1
    2724                          # X-Y+Z-
    2725                          else:
    2726                              self.statistics[comp, 5] += 1
    2727                      else:
    2728                          # X-Y-Z+
    2729                          if qVect[2] > 0:
    2730                              self.statistics[comp, 6] += 1
    2731                          # X-Y-Z-
    2732                          else:
    2733                              self.statistics[comp, 7] += 1
    2734              comp += 1
    2735
    2736 -    def __isotropicGeneration(self, qRadius, qMin, qMax):
    2737          if self.reciprocalCellVolume is None:
    2738              self.__makeRandomQList()
    2739
    2740          else:
    2741              reciprocalShellVolume = 4.0*N.pi*qRadius**2 * (qMax - qMin)
    2742              nCells = reciprocalShellVolume/self.reciprocalCellVolume
    2743              if nCells > 500:
    2744                  LogMessage('info', 'Random Q vectors generation for shell with radius Q = %s.' % qRadius, ['file','console'])
    2745                  self.__makeRandomQList(qRadius, qMin, qMax)
    2746              else:
    2747                  LogMessage('info', 'Explicit Q vectors generation for shell with radius Q = %s.' % qRadius, ['file','console'])
    2748                  self.__makeExplicitQList(qRadius, qMin, qMax)
    2749
    2750 -    def __anisotropicGeneration(self, qRadius, qMin, qMax):
    2751          LogMessage('info', 'Explicit Q vectors generation for shell with radius Q = %s.' % qRadius, ['file','console'])
    2752          self.__makeExplicitQList(qRadius, qMin, qMax)
    2753
    2754 -    def __makeRandomQList(self, qRadius, qMin, qMax):
    2755
    2756          # Generates randomly, 'number' q vectors.
    2757          # This counter is used to avoid infinite loop that can occur during the q vector generation.
    2758          redundancyCounter = 0
    2759          while len(self.qVectors[qRadius]) < self.qVectorsPerShell:
    2760              # Generate a normalized random vector.
    2761              # The actual q vector.
    2762              q = randomVector(self.qVectorsDirection)*uniform(qMin, qMax)
    2763              if self.reciprocalBasis is not None:
    2764                  # Its projection on the reciprocal basis. If the vector is not an exact reciprocal space vector, it is
    2765                  # rounded to the lower corner of its corresponding reciprocal cell.
    2766                  q = [round(v) for v in N.array(self.directMatrix*q)]
    2767                  # The approximated reciprocal space vector.
    2768                  q = q[0]*self.reciprocalBasis[0] + q[1]*self.reciprocalBasis[1] + q[2]*self.reciprocalBasis[2]
    2769
    2770              # Checks whether the length of q vector is in the user-defined length interval [lower, upper]
    2771              # If not, regenerates a new q vector
    2772              if qMin < q.length() <= qMax:
    2773                  # The generated q vector is already in the q vector list so skip it.
    2774                  if q in self.qVectors[qRadius]:
    2775                      redundancyCounter += 1
    2776                      # After 5000 trials to generate a new q vector, gives up and returns the q vector list in its current state.
    2777                      if redundancyCounter == 5000:
    2778                          return
    2779                      continue
    2780                  self.qVectors[qRadius].append(q)
    2781                  # A new q vector could be found, so reinit the redundancy counter.
    2782                  redundancyCounter = 0
    2783              else:
    2784                  redundancyCounter += 1
    2785                  # After 5000 trials to generate a new q vector, gives up and returns the q vector list in its current state.
    2786                  if redundancyCounter == 5000:
    2787                      return
    2788
    2789 -    def __makeExplicitQList(self, qRadius, qMin, qMax):
    2790
    2791          # isotropic case
    2792          if self.generator == '3d isotropic':
    2793              upperLimit = [int(qMax/v.length()) for v in self.reciprocalBasis]
    2794              # Computes all the possible combination of basis vectors and keep the ones that falls in the selected shell.
    2795              for l in range(-upperLimit[0], upperLimit[0]+1):
    2796                  for m in range(-upperLimit[1], upperLimit[1]+1):
    2797                      for n in range(-upperLimit[2], upperLimit[2]+1):
    2798                          if l == m == n == 0:
    2799                              continue
    2800
    2801                          q = l*self.reciprocalBasis[0] + m*self.reciprocalBasis[1] + n*self.reciprocalBasis[2]
    2802
    2803                          if qMin < q.length() <= qMax:
    2804                              if q in self.qVectors[qRadius]:
    2805                                  continue
    2806                              self.qVectors[qRadius].append(q)
    2807
    2808          elif self.generator == '2d isotropic':
    2809              upperLimit = [int(qMax/v.length()) for v in self.qVectorsDirection]
    2810              for l in range(-upperLimit[0], upperLimit[0]+1):
    2811                  for m in range(-upperLimit[1], upperLimit[1]+1):
    2812                      if l == m == 0:
    2813                          continue
    2814
    2815                      # The actual q vector.
    2816                      q = l*self.qVectorsDirection[0] + m*self.qVectorsDirection[1]
    2817
    2818                      # Its projection on the reciprocal basis vector. If the vector is not an exact reciprocal space vector, it is
    2819                      # rounded to the lower corner of its corresponding reciprocal cell.
    2820                      q = [round(v) for v in N.array(self.directMatrix*q)]
    2821                      # The approximated reciprocal space vector.
    2822                      q = q[0]*self.reciprocalBasis[0] + q[1]*self.reciprocalBasis[1] + q[2]*self.reciprocalBasis[2]
    2823
    2824                      if qMin < q.length() <= qMax:
    2825                          if q in self.qVectors[qRadius]:
    2826                              continue
    2827                      self.qVectors[qRadius].append(q)
    2828
    2829          elif self.generator == 'anisotropic':
    2830
    2831              for qVect in self.qVectorsDirection:
    2832                  upperLimit = int(qMax/qVect.length())
    2833                  for l in range(-upperLimit, upperLimit + 1):
    2834                      if l == 0:
    2835                          continue
    2836
    2837                      # Generate a normalized random vector.
    2838                      # The actual q vector.
    2839                      q = l*qVect
    2840
    2841                      # Its projection on the reciprocal basis vector. If the vector is not an exact reciprocal space vector, it is
    2842                      # rounded to the lower corner of its corresponding reciprocal cell.
    2843                      q = [round(v) for v in N.array(self.directMatrix*q)]
    2844                      # The approximated reciprocal space vector.
    2845                      q = q[0]*self.reciprocalBasis[0] + q[1]*self.reciprocalBasis[1] + q[2]*self.reciprocalBasis[2]
    2846
    2847                      if qMin < q.length() <= qMax:
    2848                          if q in self.qVectors[qRadius]:
    2849                              continue
    2850                          self.qVectors[qRadius].append(q)
    2851
    2852          # If there are more q vectors generated than the user-defined number of q vectors per shell
    2853          # then pick up randomly |number| q vectors among the generated q vectors.
    2854          if len(self.qVectors[qRadius]) > self.qVectorsPerShell:
    2855              LogMessage('warning', '%s qVectors generated for shell with radius Q = %s. Will be reduced to %s.' % (len(self.qVectors[qRadius]), qRadius, self.qVectorsPerShell), ['file','console'])
    2856              self.qVectors[qRadius] = sample(self.qVectors[qRadius], self.qVectorsPerShell)
    2857